GMTで特定の海域を塗り潰す

描き方を聞かれて少し迷ったのでメモ。
地球科学分野ではあまり使わないテクニックだが、海底地形を必要としない海域図では利用価値がある。

地図上で任意の領域を塗り潰すためには psxy の -G オプションを使えば良い。しかし沿岸部の場合、海岸線での塗り分けが必要になる。このような場合、psxy で領域を描画した後に pscoast で上書きすると良い。

以下に、大阪湾で例を示す。

例では 8 点のみで湾を囲む。海岸線の数値データが手元になく、領域の正確な範囲を作成するのが困難な場合でも簡易的に用いることができる方法となる。必要な緯度経度は Google Maps でポイントをクリックすると得られるもので十分。

湾内を塗りたい場合、以下のようなスクリプトで実現できる。
（※大阪湾の正確な定義に基づく描画ではない）

paint_area.csh

#!/bin/csh -f

#gmtdefaults の設定
gmtset PAPER_MEDIA=A4
gmtset PAGE_ORIENTATION=portrait
gmtset LABEL_FONT_SIZE=18p
gmtset PLOT_DEGREE_FORMAT=ddd:mm:ssF

#出力ファイル
set psfile=osakabay.ps

#今回は、県境データが以下の名前で存在すると仮定
set prefecture=b_pref.txt

#描画オプションの指定
set Region=134.25/135.75/33.75/35.25
set Projection=M15
set Annotation=f.5a.5g.5WSne

#Basemap の描画（グリッドなし）　※グリッドを陸上にも描きたい場合は pscoast の後に記述するべし
psbasemap -R$Region -J$Projection -B$Annotation -K -V > $psfile

#塗りたい領域。簡易的に大阪湾より広いエリアを指定
psxy -R -J -G100/100/100 -m -L -: -K -V -O -f << END >> $psfile
34:36.555, 135: 0.183
34:38.280, 135:2.066
34:45.926, 135:24.653
34:39.187, 135:32.110
34:16.633, 135:19.462
34:17.840, 135:4.181
34:16.122, 134:56.474
34:16.973, 134:45.679
END

#海岸線の描画
pscoast -R -J -Df -W2 -G220 -P -K -V -O >> $psfile

#県境の描画
psxy $prefecture -R -J -W1/0/0/0 -m -O -V >> $psfile

上記のスクリプトを実行すると、次のような結果を得る。