Python で海底地形図を描く (1)

最近，上のような海底地形図が Python で描けることを知ったのでまとめておきます．地球全体の標高マップについての記事を参考に，特定の緯度経度だけを切り抜いて描画します．

データセットのダウンロード

NOAA のサイトから，ETOPO1_Bed_g_gmt4.grd.gz というファイルをダウンロードし，解凍します．gz ファイルの解凍はこちらが参考になるかと思います (windows 限定です)．
解凍した ETOPO1_Bed_g_gmt4.grd というファイルを，プログラムコードと同じフォルダに入れておきます．

ライブラリのインポート

今回使うライブラリは以下です．netCDF4 だけはインストールが必要になる方が多いかもしれません．

from netCDF4 import Dataset
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as colors
# from mpl_toolkits.basemap import Basemap
from sklearn import preprocessing

# フォントの指定
mpl.rcParams['font.size'] = 12
mpl.rcParams['font.family'] = 'Arial'
mpl.rcParams['pdf.fonttype'] = 42

最終行の mpl.rcParams['pdf.fonttype'] = 42 は，出力する PDF ファイルをベクター形式にします．これを入れると，イラストレーターで編集ができるようになります．

データの読み込み

標高のグリッドデータを読み込みます．

etopo1 = Dataset('ETOPO1_Bed_g_gmt4.grd','r')
x=etopo1.variables['x'][:]
y=etopo1.variables['y'][:]
Z=etopo1.variables['z'][:]
X, Y = np.meshgrid(x, y)

このデータは，緯度方向 10801 点×経度方向 21601 点という非常に重いデータなので，数点おきにデータをとって使用します．今回は 3 点につき 1 点をとり，0.05° おきのグリッドデータを描画に使います．
解像度をどこまで落とすかは調整が必要になると思うので，データを読み込むセルとは分けておくことをおすすめします．

step = 3
X_rough = X[::step, ::step]
Y_rough = Y[::step, ::step]
Z_rough = Z[::step, ::step]

np.set_printoptions(precision = 3, floatmode='fixed')
print("**** X_rough ****")
print(X_rough[:5, :5])
print("**** Y_rough ****")
print(Y_rough[:5, :5])
print("**** Z_rough ****")
print(Z_rough[:5, :5].astype(float))

グリッドデータの意味

X, Y, Z は同じ shape の NumPy 配列になっていて，ある番地を指定するとその点の緯度・経度・標高が出力されます．緯度経度はそれぞれ北緯，東経が正で，標高の負の値は水深を表します．

XY_index = (2, 3)
print("longitude: {:.3f} (E), latitude: {:.3f} (N), height: {} (m)"\
     .format(X_rough[XY_index], Y_rough[XY_index], Z_rough[XY_index]))

longitude: -179.850 (E), latitude: -89.900 (N), height: -51 (m)

X_rough, Y_rough, Z_rough の 3 行 4 列の値が抽出されています．

描画範囲のデータ抽出

今回は，地球で最も深いマリアナ海溝周辺の地形図を描いてみます．描画範囲は，東経 135~155°・北緯 5~25° くらいでしょうか．

# 描画範囲を決定
left_edge, right_edge = 135, 155
bottom_edge, top_edge = 5, 25

# 描画範囲に含まれるかを判定
X_isin_drawRange = (left_edge <= X_rough[0, :])\
                  * (X_rough[0, :] <= right_edge)
Y_isin_drawRange = (bottom_edge <= Y_rough[:, 0])\
                  * (Y_rough[:, 0] <= top_edge)

# 描画範囲のデータだけを抽出する
X_drawRange = X_rough[Y_isin_drawRange, :][:, X_isin_drawRange]
Y_drawRange = Y_rough[Y_isin_drawRange, :][:, X_isin_drawRange]
Z_drawRange = Z_rough[Y_isin_drawRange, :][:, X_isin_drawRange]

データの抽出の部分が少し汚いコードになっています．[Y_isin_drawRange, X_isin_drawRange] で指定しようとしたらエラーが出てしまいました…

カラーバーの作成

こちらのサイトを参考に，塗り分けだけ変更しました．

color = np.array([[-10000, 'darkblue'],
                 [-9000,  'blue'],
                 [-6000,  'deepskyblue'],
                 [-4000,  'darkseagreen'],
                 [-3000,  'yellowgreen'],
                 [-2000,  'orange'],
                 [-1000,  'bisque'],
                 [0,      'whitesmoke']])

####################################################

height = color[:, 0].astype(np.float32)
heightnorm = preprocessing.minmax_scale(height)
colornorm = []
for no, norm in enumerate(heightnorm):
   colornorm.append([norm, color[no, 1]])
cmap = colors.LinearSegmentedColormap.from_list('a', colornorm, N=height.max()-height.min()+1)
heightbounds = range(int(height.min()), int(height.max()+1))
cmap.set_over("black") # 陸地は黒で塗る
cmap.set_under(color[0, 1])
norm = colors.BoundaryNorm(heightbounds, cmap.N)

地形図の描画

matplotlib で標高データを表示します．せっかくなので，マリアナ海溝の最深部であるチャレンジャー海淵に星を打ってみます．
PDF ファイルの出力には数十秒かかります．満足のいく図が得られるまでは，図の保存はコメントアウトしておくことをおすすめします．

# カラーマップを作成する
fig = plt.figure(figsize=(18.3, 15), 
                facecolor = 'w', edgecolor = 'k')
plt.subplots_adjust(wspace=0.1, hspace=0.1)
ax = fig.add_subplot(1, 1, 1)
pcolor = ax.pcolor(X_drawRange, Y_drawRange, Z_drawRange, 
                  norm = norm, linewidth = 0, cmap = cmap)

# 標高バーをつける
cb = fig.colorbar(pcolor, shrink=0.6, aspect=20, extend='both', 
                 ticks = color[:, 0].astype(np.float32), spacing='proportional')
cb.set_label('Depth')

# マリアナ海溝最深部
ChallengerDeep = {"lat": 11.373333, 
                 "lon": 142.591667}

# サイトを星で表示
ax.scatter(ChallengerDeep["lon"], ChallengerDeep["lat"], 
          s = 500, marker = "*", color = "r", 
          linewidths = 2, edgecolors = "w")

# 描画した地形図を保存
plt.savefig("MarianaTrench.png")
plt.savefig("MarianaTrench.pdf")

plt.show()

こんな感じの海底地形図が出力されます．星のサイズはもう少し大きくすればよかったか…

基本的な描画はできるようになりました．次回以降，もう少し遊んでみようと思います．

追記
Python で海底地形図を描く (2) を書きました (2020/12/26)．
Basemap のインポートは不要だったのでコメントアウトしました (2024.3.4)．

ライブラリのバージョン

今回使用したライブラリのバージョンは以下の通りです．バージョンの確認方法については以前の記事に書いたので，興味があればあわせてどうぞ．

netCDF4 : 1.4.2,    numpy : 1.16.2,
sklearn : 0.20.3,    matplotlib: 3.0.3

参考​

NOAA: ETOPO1 Global Relief Model
netCDF4 公式ページ
ETOPO1グローバル地形データセットで地球の全体地図をPythonで描いてみる
地球全体を、標高1m毎に別の色で塗り絵して地図を作ってみる。
watermark によるライブラリのバージョン表示 (Python)