図の大きさから-Rオプションを決定する

GMT ver.5で、緯度経度がわかっている複数の点を全て地図中にプロットしようと思います。メルカトル図法を使用するならば、例えば

 gmt psxy -JM20c -Rwest/east/south/north -B

のように地図の幅を20cmとして描画することができます。ただし地図の高さが何センチメートルになるかは指定できません。緯度と経度は縮尺を介して従属関係にあるので当然と言えば当然ですが、この為にしばしば図-1のように思いもよらない縦長の絵になってしまうことがあります。

f:id:S_E_Hyphen:20180122093924j:plain

図-1 横幅を大きく指定しすぎた例

 

 画面では余り気になりませんが、A4ランドスケープで描画した場合は図面の上半が切れてしまいます。そこで、図の幅と高さを指定して緯度経度の範囲を再決定するシェルスクリプトを作成してみました。リスト-1に示します。

指定した緯度経度では、所望のアスペクト比よりも縦長になる場合、東西に範囲を広げて自動的に指定したアスペクト比となるようにします(図-2の例)。逆に横長になってしまった場合は、南北に範囲を拡張します。

 

f:id:S_E_Hyphen:20180122095009j:plain

図-2 gmt_rangesetにより東西を自動的に拡張した例

 

メルカトル図法以外には対応していませんし、 緯度・経度の比を求めるのに球面近似を使用しているなど不十分な点はありますが、とりあえず見れる図面を描くためには便利なツールだと思っています。

 

リスト-1 gmt_rangeset のリスト
#!/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

 

 

 

リスト-2 gmt_rangeset の使用例
#!/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