Python で海底地形図を描く (2)
前回の続きです．前回は，ETOPO1 のデータセットから指定した緯度経度の範囲のデータだけを取り出して，マリアナ海溝周辺の海底地形図を描きました．これを応用して，今度は太平洋の地形図を描いてみようと思います．
ここまで書くと，「描画範囲を変えるだけじゃん」と思われてしまいそうですが，そのままではエラーになります．描画に用いている 標高（&水深）のデータは経度 180° で切れているので，180° 線をまたぐ処理が必要になります．
考え方
オリジナルの地形データの右側に，同じデータセットをもう一つくっつけます．これが作れれば，連続した描画範囲を指定できるので，前回と同様の方法で描画できるはずです．
コード
最初の方は前回と同じです．まずはライブラリのインポートから．
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'] = 14
mpl.rcParams['font.family'] = 'Arial'
mpl.rcParams['pdf.fonttype'] = 42 # pdf の出力をベクター形式にする次に，ETOPO1 のデータセットを読み込みます．経度の -180° と 180° で重複があるので，今回は 180° のデータを除外しておきます．データ点の間隔は，前回よりも粗めに設定しています．
# データセットの読み込み
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)
# 重複している 180 度のデータを除外
X = X[:, :-1]
Y = Y[:, :-1]
Z = Z[:, :-1]
# データを間引く
step = 12
X_rough = X[::step, ::step]
Y_rough = Y[::step, ::step]
Z_rough = Z[::step, ::step]次に，元データの経度に 360° (地球 1 周分) 加えた新たなデータセットを作ります．緯度や標高はそのままです．そして，このデータセットを元データの右側に結合します．
# 右側のデータ
X_right = X_rough + 360
Y_right = Y_rough.copy()
Z_right = Z_rough.copy()
# 結合
X_stack = np.hstack([X_rough, X_right])
Y_stack = np.hstack([Y_rough, Y_right])
Z_stack = np.hstack([Z_rough, Z_right])描画範囲を指定します．180° をまたぐ場合は，right_edge に 360 を加えることで (例えば西経 60° は東経 300° とみなせます)，連続的な描画範囲が指定できるようになります．
# 描画範囲を決定
left_edge, right_edge = 120, -60
bottom_edge, top_edge = -60, 60
# 180°をまたぐ場合の処理
if left_edge > right_edge:
right_edge += 360描画範囲のデータ抽出は，前回と同じ処理です．
# 描画範囲に含まれるかを判定
X_isin_drawRange = (left_edge <= X_stack[0, :])\
* (X_stack[0, :] <= right_edge)
Y_isin_drawRange = (bottom_edge <= Y_stack[:, 0])\
* (Y_stack[:, 0] <= top_edge)
# 描画範囲のデータだけを抽出する
X_drawRange = X_stack[Y_isin_drawRange, :][:, X_isin_drawRange]
Y_drawRange = Y_stack[Y_isin_drawRange, :][:, X_isin_drawRange]
Z_drawRange = Z_stack[Y_isin_drawRange, :][:, X_isin_drawRange]水深ごとの色を設定します．前回と同様，こちらのコードを参考にさせて頂きました．
# カラーバーの設定
color = np.array([[-7000, 'darkblue'],
[-6000, 'blue'],
[-5000, 'deepskyblue'],
[-4000, 'darkseagreen'],
[-3000, 'yellowgreen'],
[-2000, 'orange'],
[-1000, 'bisque'],
[0, 'whitesmoke']])
####################################################
# 水深・標高用カラーマップ作成
# 標高値を0..1の範囲に正規化
height = color[:, 0].astype(np.float32)
heightnorm = preprocessing.minmax_scale(height)
# 正規化した標高値と各色の組みのリストcolornormを生成
colornorm = []
for no, norm in enumerate(heightnorm):
colornorm.append([norm, color[no, 1]])
# colornormに従いグラデーションカラーマップを作成
cmap = colors.LinearSegmentedColormap.from_list('a', colornorm, N=height.max()-height.min()+1)
# カラーマップを割り当てる1m間隔の標高値
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)描画も基本的には前回と同じです．ただし，経度の軸ラベル (xticklabels) は，180° を超えた値のみ 360 を引くという処理をしています．
fig = plt.figure(figsize = (27, 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)
# X 軸を振り直す
x_tick0 = ax.get_xticks()
ax.set_xticklabels(np.where(x_tick0 > 180,
x_tick0 - 360, x_tick0).astype(int))
# 標高バーをつける
cb = fig.colorbar(pcolor, shrink=0.6, aspect=20, extend='both',
ticks = color[:, 0].astype(np.float32),
spacing='proportional')
cb.set_label('Depth')
# 書き出した地形図を保存
plt.savefig("PacificMap.png")
plt.savefig("PacificMap.pdf")
plt.show()出力はこんな感じになります．
ライブラリのバージョン
以下のバージョンで動作確認を行っています．バージョンの確認方法については，こちらの記事も合わせてどうぞ．
numpy : 1.16.2 matplotlib: 3.0.3
sklearn : 0.20.3 netCDF4 : 1.4.2
追記
PyGMT で海底地形図を描くを書きました．よりシンプルなコードで実装したい方はこちらもご覧ください (2024.8.7)．