GMT6.5のサンプルスクリプトを試す3 basemap~極座標~

GMT（Generic mapping tools）6.5の公式HPに記載されている例・サンプルスクリプトが示されているが、図が示されていないので描いてみる。

GMT6のサンプルスクリプト試す2 basemap~べき乗の続き

引用元
GMT6.5.0 Documentation basemap

極座標表示

gmt basemap -R0/360/0/1000 -JP6i -Bafg -pdf polar

-R:極を中心にプロットさせる場合は0/360にする。
-JP pまたはPで極座標にできる。P6iでは半径3インチ、直径6インチになる。

せっかくだから何かプロットしてみようと思い立って、QZSSの衛星の軌跡を描いてみることにした。
いつもお世話になっている国土地理院の電子基準点データから、適当な点の適当な日の適当なデータをダウンロードしてみる。そして国土地理院が配布しているGNSSデータ前処理ツールringoを使って、人工衛星の方位角と仰角を出力してみた。

ringo rnxcsv 06300010.24o 06300010.24q > 06300010.csv
grep J02 06300010.csv　| awk ‘{ print $4, $5; }‘  > sample4.txt

衛星の軌跡を描くならば、北が0度、中心（天頂）が90度で円周（地表）が0度にするほうがよい。
この場合は、-Jに+aと+fをつければよいらしい。そして範囲を0/90にする。
電子基準点の観測データは３０秒間隔なので、３０秒ごとの衛星の位置をプロットすることにするので、plotを使って小さい赤丸を連続的に描画する。

gmt begin qzss png
gmt basemap -R0/360/0/90 -JP6i+a+f -Bafg
gmt plot sample4.txt -Sc0.05 -Gred
gmt end

QZSSの8の字が描かれた。