damyarou

python, GMT などのプログラム

GMT ETOPO1で地形図作成

記事の最後に行く

はじめに

この記事は、Qiita に投稿した以下の記事の再アレンジ版です。

ETOPO1 は,National Geophysical Data Center (Colorado) により作成された,全球の1分刻みの地形モデルです。 ちょっとミャンマー国周辺の地形を見たかったので、ETOPO1 のデータを GMT で可視化してみました。

作業環境は以下の通り。使っているプログラムは GMT (Generic Mapping Tools) と bash のみです。

参考サイト

まずは作例を

ETOPO1 のデータから,標高 500m ピッチで色分けし,描画領域の主要都市をプロットしました。 マスターカラーテーブルは Qiitaの記事では no_green を使いましたが、今回はより地形図に合うと思った haxby を使いました。 陸地側標高は、チベット高原まで含むので、標高5000mまでの表示にしました。海側は陸地と対称に-5000mまでです。 国境線を黒い点線で、主要河川を青い実線で示しています。

f:id:damyarou:20190519083940j:plain

ETOPO1のデータを取ってくる

ETOPO1 のサイトはこれ.ETOPO1 Ice Surfaceのgrid-registered: netCDF 形式を選択.

データのダウンロードは以下より.ETOPO1_Ice_g_gmt4.grd.gz をダウンロードし作業したいディレクトリに展開しておく.

Name Last modified   Size
Parent Directory        -
ETOPO1_Ice_g_gdal.grd.gz    2011-03-22 20:52    377M
ETOPO1_Ice_g_gmt4.grd.gz    2013-01-02 18:00    377M  <= これをダウンロードして展開
readme_etopo1_netcdf.txt    2010-09-23 18:29    1.1K

ETOPO1_Ice_g_gmt4.grd.gz を展開すると ETOPO1_Ice_g_gmt4.grd (933.5MB) が得られますので、作図にはこれを用います。

コマンドの説明

変数

src=ETOPO1_Ice_g_gmt4.grd #元になるグリッドファイル
range=90/105/15/30 # 作図範囲(最小経度 / 最大経度 / 最小緯度 / 最大緯度)
size=15            # 出力画像幅 (cm)
file=myanmer       # 作業ファイル名
fig=fig_$file.eps  # 出力画像ファイル名

grdcut

gmt grdcut $src -R$range -G$file.grd
  • $src:入力グリッドファイル
  • $range:切り取る範囲
  • $file.grd:出力グリッドファイル

makecpt

gmt makecpt -Chaxby -T-5000/5000/500 > $file.cpt
  • haxby:マスターカラーテーブル
  • -T-5000/5000/500:カラーテーブルの範囲と増分(標高-5000mから+5000mまで500mピッチ)
  • $file.cpt:出力CPTファイル名

grdimage

gmt grdimage $file.grd -R$range -JM$size -C$file.cpt -P -K > $fig
  • $file.grd:入力グリッドファイル
  • -R$range:描画範囲
  • -JM$size:プロジェクション(メルカトル)と画像の大きさ(cm)
  • -C$file.cpt:使用するCPTファイル

grdcontour

gmt grdcontour $file.grd -R -J -C500 -W0 -L0/5000 -A500 -O -K >> $fig
  • $file.grd:入力グリッドファイル
  • -C500:等高線ピッチ
  • -W0:等高線の線の太さはゼロ(描かない).指定しないとデフォルトで描かれる
  • -L0/5000:等高線を描く範囲(標高0m以上5000mまで等高線を描く)
  • -A500:等高線に描く数値のピッチ

psscale

gmt psscale -Ba2500g500f500 -C$file.cpt -D6c/-1c/6c/0.3ch -O -K >> $fig
  • -Ba2500g500f500:スケールに描く体裁
  • a:文字を2500ピッチで,g:グリッド線を500ピッチで,f:補助線を500ピッチで描く
  • -C$file.cpt:用いるCPTファイル名
  • -D6c/-1c/6c/0.3ch:スケールの描画位置,長さ,幅,水平スケールの指定
    ここでは,x方向6cm,y方向-1cmを中心に長さ6cm,幅0.3cmのスケールを描画.hの指定で水平スケールとなる.

pscoast

gmt pscoast -R -J -N1/2p,#000000,2_2:0 -I1/1p,#0000ff -I2/1p,#0000ff -O -K >> $fig
  • pscoast:海岸線の描画
  • -N1/2p,#000000,2_2:0:境界線の描画(N1:国境,線の太さ2pt,線の色は黒,線種は点線)
  • -I1/1p,#0000ff:主要河川
  • -I2/1p,#0000ff:付加的主要河川

psxy

描画する都市の座標に、白抜き丸を打ちます。 座標は文字情報のファイルに格納されており、awk により、経度・緯度のみの情報をGMTに渡しています。

awk '{print $1,$2}' _city.txt | gmt psxy -R -J -SC0.3 -W2,#000000 -G#ffffff -O -K >> $fig
  • psxy:線を引く
  • -SC0.3 -W2,#000000 -G#ffffff:0.3cmの円(SC0.3)を太さ2pt黒実線で描き(-W2,#000000)、中を白で塗る(-G#ffffff)。

pstext

都市名を描画します。

awk '{print $1+0.2,$2-0.2,$3,$4,$5,$6}' _city.txt | gmt pstext -R -J -F+f+a+j -O -K >> $fig
  • pstest:文字列の描画
  • -F+f+a+j:座標データのあとにフォント情報,描画角度,指定座標に対する位置を指定
  • ここでは、awk を用いて読み込んだ座標に対し、経度は +0.2度、緯度は -0.2度を付加して文字のプロット座標を指定しています。
  x(lon)  y(lat)     f              a  j  text
---------------------------------------------------------
 90.3669 23.7020 12p,Helvetica-Bold 0 ML Dhaka
 91.8000 22.3667 12p,Helvetica-Bold 0 TC Chattogram
 96.1500 16.8000 12p,Helvetica-Bold 0 TC Yangon
 96.1000 19.7500 12p,Helvetica-Bold 0 TC Naypyidaw
 96.0667 21.9667 12p,Helvetica-Bold 0 TC Mandalay
102.6000 17.9667 12p,Helvetica-Bold 0 TC Vientiane
 91.1194 29.6581 12p,Helvetica-Bold 0 ML Lhasa 
102.6830 25.0667 12p,Helvetica-Bold 0 ML Kunming

psbasemap

gmt psbasemap -R -J -Bxa5f1 -Bya5f1 -O >> $fig
  • psbasemap:枠線の描画
  • -Bxa5f1:x軸方向の指定(数値を5単位で,補助線を1単位で描画)
  • -Bya5f1:y軸方向の指定(数値を5単位で,補助線を1単位で描画)

psconvert

epsを余白削除してjpgに変換します。

gmt psconvert $fig -Tj -A
  • $fn に eps のファイル名を入れます.
  • -Tj で jpg に変換.デフォルトでは出力ファイル名は拡張子のみが変わるようです.
  • -A で余白削除.
  • 解像度の指定もできるようですが,デフォルトでは 300dpi のようです.
  • 変換できる画像は,BMP, EPS, JPEG, PDF, PNG, PPM, SVG, TIFF ということです.

スクリプト全文

スクリプトでは最初に _city.txt という名前で、都市名プロット用のデータファイルを作成しています。

# City plot data
cat << EOT > _city.txt
 90.3669 23.7020 12p,Helvetica-Bold 0 ML Dhaka
 91.8000 22.3667 12p,Helvetica-Bold 0 TC Chattogram
 96.1500 16.8000 12p,Helvetica-Bold 0 TC Yangon
 96.1000 19.7500 12p,Helvetica-Bold 0 TC Naypyidaw
 96.0667 21.9667 12p,Helvetica-Bold 0 TC Mandalay
102.6000 17.9667 12p,Helvetica-Bold 0 TC Vientiane
 91.1194 29.6581 12p,Helvetica-Bold 0 ML Lhasa 
102.6830 25.0667 12p,Helvetica-Bold 0 ML Kunming
EOT


# ====================
# GMT command
# ====================
src=ETOPO1_Ice_g_gmt4.grd
range=90/105/15/30 # W/E/S/N (range of drawing)
size=15            # width of drawing (cm)
file=myanmer       # part of drawing name
fig=fig_$file.eps  # output image file name
gmt grdcut $src -R$range -G$file.grd
gmt makecpt -Chaxby -T-5000/5000/500 > $file.cpt
gmt grdimage $file.grd -R$range -JM$size -C$file.cpt -P -K > $fig
gmt grdcontour $file.grd -R -J -C500 -W0 -L0/5000 -A500 -O -K >> $fig
gmt psscale -Ba2500g500f500 -C$file.cpt -D6c/-1c/6c/0.3ch -O -K >> $fig
gmt pscoast -R -J -N1/2p,#000000,2_2:0 -I1/1p,#0000ff -I2/1p,#0000ff -O -K >> $fig
awk '{print $1,$2}' _city.txt | gmt psxy -R -J -SC0.3 -W2,#000000 -G#ffffff -O -K >> $fig
awk '{print $1+0.2,$2-0.2,$3,$4,$5,$6}' _city.txt | gmt pstext -R -J -F+f+a+j -O -K >> $fig
gmt psbasemap -R -J -Bxa5f1 -Bya5f1 -O >> $fig

gmt psconvert $fig -Tj -A

Thank you.

記事の先頭に行く