GMT ver.5で、緯度経度がわかっている複数の点を全て地図中にプロットしようと思います。メルカトル図法を使用するならば、例えば
gmt psxy -JM20c -Rwest/east/south/north -B
のように地図の幅を20cmとして描画することができます。ただし地図の高さが何センチメートルになるかは指定できません。緯度と経度は縮尺を介して従属関係にあるので当然と言えば当然ですが、この為にしばしば図-1のように思いもよらない縦長の絵になってしまうことがあります。
画面では余り気になりませんが、A4ランドスケープで描画した場合は図面の上半が切れてしまいます。そこで、図の幅と高さを指定して緯度経度の範囲を再決定するシェルスクリプトを作成してみました。リスト-1に示します。
指定した緯度経度では、所望のアスペクト比よりも縦長になる場合、東西に範囲を広げて自動的に指定したアスペクト比となるようにします(図-2の例)。逆に横長になってしまった場合は、南北に範囲を拡張します。
メルカトル図法以外には対応していませんし、 緯度・経度の比を求めるのに球面近似を使用しているなど不十分な点はありますが、とりあえず見れる図面を描くためには便利なツールだと思っています。
#!/bin/bash sdoc=\ " gmt_rangeset - 図面の大きさから\n"\ " 緯度経度の範囲を\n"\ " を設定する \n"\ " ※メルカトル図法(-JMオプション)\n"\ " に対応します \n"\ " \n"\ " [ require] \n"\ " width 地図の幅(単位:cm) \n"\ " height 地図の高さ(単位:cm) \n"\ " range 経度・緯度の最大最小値 \n"\ " -Rwest/east/south/northで指定\n"\ " \n"\ " [ 出力 ] \n"\ " -JM{width}c -Rw/e/s/n \n"\ " 図の高さがheight(cm)より小さい\n"\ " ことが予想される場合は \n"\ " n=north+α s=south+α して \n"\ " 図の高さをheight(cm)になるよう\n"\ " 調節する \n"\ " 逆に大きくなる場合は \n"\ " e=east+α w=west-α とする \n"\ " \n" export sdoc width=`getparstr $# "$*" "width"` height=`getparstr $# "$*" "height"` range=`getparstr $# "$*" "range"` if [ -z $width ]; then usage;exit; fi if [ -z $height ]; then usage;exit; fi if [ -z $range ]; then usage;exit; fi if [ `echo $range | cut -c1-2` != "-R" ]; then usage;exit; fi xmin=`echo $range | sed -e "s/-R//g" | awk -F"/" '{print $1}'` xmax=`echo $range | sed -e "s/-R//g" | awk -F"/" '{print $2}'` ymin=`echo $range | sed -e "s/-R//g" | awk -F"/" '{print $3}'` ymax=`echo $range | sed -e "s/-R//g" | awk -F"/" '{print $4}'` ymean=`echo $range | sed -e "s/-R//g" | awk -F"/" '{print 0.5*($3+$4)}'` prjx=$width prjy=$height temp=temp.txt scilab -nw << _SCI_ xmin=${xmin};xmax=${xmax};ymin=${ymin};ymax=${ymax}; height=1/cos($ymean*%pi/180)*(ymax-ymin)/(xmax-xmin) * $prjx; ratio=height/$prjy; fd=mopen('$temp','w'); if ratio<1, ymax=$ymean+($ymax-$ymin)/ratio/2; ymin=$ymean-($ymax-$ymin)/ratio/2; mfprintf(fd,'-JM%sc ',"${prjx}"); mfprintf(fd,'-R%f/%f/%f/%f\n',xmin,xmax,ymin,ymax); else xmin=0.5*$xmin*(1+ratio)+0.5*$xmax*(1-ratio); xmax=0.5*$xmin*(1-ratio)+0.5*$xmax*(1+ratio); mfprintf(fd,'-JM%sc ',"${prjx}"); mfprintf(fd,'-R%f/%f/%f/%f\n',xmin,xmax,ymin,ymax); end mclose(fd); _SCI_ cat $temp; rm $temp
#!/bin/bash cat <<+ > GPSsite.txt 利尻 960501 45.1377 141.1671
滝上 960504 44.1919 143.0768
愛別 960506 43.9084 142.5777
旭川 960508 43.7385 142.4096
北竜 960509 43.7408 141.8738
上川 960510 43.7714 142.9030
富良野 960514 43.3363 142.3953
三笠 960516 43.2484 141.8909
小樽2 960517 43.2095 140.8606
新得 960518 43.1655 142.8097
江別 960520 43.0776 141.5401
帯広 960521 42.9389 143.1706
恵庭 960522 42.8844 141.5774
千歳 960523 42.7732 141.4071
虻田 960525 42.5514 140.7676
伊達 960526 42.4739 140.8764
奥尻2 960527 42.0605 139.4462
砂原 960528 42.1231 140.6669
七飯 960529 41.9767 140.7154
上ノ国 960530 41.8030 140.0712
+ proj=-JM10c range=`\ cat GPSsite.txt |\ awk '{print $4,$3}' |\ gmtinfo -I0.5/0.5 ` range=`gmt_rangeset width=20 height=15 range=$range | awk '{print $1,$2}'` proj="" echo $proj $range gmt psbasemap $proj $range -Baf -K > map.ps gmt pscoast -J -R -B -Ggray -Df -O -K >> map.ps cat GPSsite.txt |\ awk '{print $4,$3}' |\ gmt psxy -J -R -B -Ss0.5 -Gred -O -K >> map.ps cat GPSsite.txt |\ awk '{print $4,$3,$1}' |\ gmt pstext -J -R -O -F+jMR+a0.0+f12p,GothicBBB-Medium-UniJIS-UTF8-H,0/0/0 -Dj0.25c >> map.ps