PyGMT で海底地形図を描く

以前に「Python で海底地形図を描く」という記事を書きました。この時は matplotlib でデータをゴリゴリいじっていたのですが、pygmt というライブラリを使うことでスマートにプログラムできることが分かったのでまとめておきます。

準備

環境構築

公式サイトに従って、以下でインストールします。ライブラリの依存関係があるようなので、仮想環境を作るのが良いと思います。

conda install --channel conda-forge pygmt

2024.9.27 追記：
初心者の方向けに，Google Colaboratory を用いた環境構築の方法を書きました．環境構築でトラブってしまった方はまずは以下を試してみてください．

【初心者向け】Python + GMT の (たぶん) 一番簡単な環境構築｜mim
https://note.com/mim_2020/n/n646cfd0288d7

グリッドデータの入手

データセットをダウンロードします。以前の記事と同様に、NOAA のサイトから、ETOPO1_Bed_g_gmt4.grd.gz というファイルをダウンロード・解凍します。

チュートリアル

解凍したグリッドデータと同じフォルダに test.ipynb をつくり、公式サイトにある以下のプログラムを試してみましょう。

import pygmt
pygmt.show_versions()

PyGMT information:
 version: v0.12.0
System information:
 python: 3.12.4 …

出力（2024.8.6 時点）

fig = pygmt.Figure()
fig.coast(projection="H10c", region="g", frame=True, land="gray")
fig.show()

3 行で世界地図が書けてしまいました

地形図の作成

以下のセルを実行します。ダウンロード・解凍したグリッドデータがプログラムコードと同じフォルダにある想定です。

grid = 'ETOPO1_Bed_g_gmt4.grd' #グリッドファイルへのパス
region = [135, 155, 5, 25]

fig = pygmt.Figure()
fig.grdimage(grid, 
             cmap = 'haxby', 
             region = region,
             frame = ['WSen', 'xaf', 'yaf']
             )
fig.colorbar()
fig.show()

地形図が描画できました。ここに情報を足していきます。

# マリアナ海溝最深部
ChallengerDeep = {"lat": 11.373333, 
                  "lon": 142.591667}
# 点を打つ
fig.plot(x = ChallengerDeep["lon"], 
         y = ChallengerDeep["lat"], 
         fill = 'black',
         style = 'a0.7c',
         pen = '1p,white'
         )
# 文字を書く
fig.text(text = 'Challenger Deep',
         x = ChallengerDeep['lon'] + 0.5,
         y = ChallengerDeep['lat'],
         font = '12p,black',
         justify='LM' # Left Middle
         )

fig.show()

180°をまたぐ場合

緯度経度をテキストで指定してあげるだけです。

region = ['120E', "120W", '60S', '60N']

fig = pygmt.Figure()
fig.grdimage(grid, 
             cmap = 'haxby', 
             region = region,
             frame = ['WSen', 'xaf', 'yaf'],
             projection = 'M12c',
             )
fig.colorbar()
fig.show()

これまでのプログラムと比べると 10 分の 1 くらいのコード量になっているかと思います。

補足

緯度経度を分-秒表記で指定したい場合は、以下のように入力します。

region = ['140:30', '142:15', '10:45:30', '12:50:15']

fig = pygmt.Figure()
fig.grdimage(grid, 
             cmap = 'haxby', 
             region = region,
             frame = ['WSen', 'xaf', 'yaf']
             )
fig.colorbar()
fig.show()

以前のプログラムと比較すると、 10 分の 1 くらいのコードで地形図が描画できるようになりました。

参考

公式サイト：https://www.pygmt.org/latest/

PyGMT-HOWTO: https://tktmyd.github.io/pygmt-howto-jp/pygmt.html