見出し画像

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

前回の続きです.前回は,ETOPO1 のデータセットから指定した緯度経度の範囲のデータだけを取り出して,マリアナ海溝周辺の海底地形図を描きました.これを応用して,今度は太平洋の地形図を描いてみようと思います.

ここまで書くと,「描画範囲を変えるだけじゃん」と思われてしまいそうですが,そのままではエラーになります.描画に用いている 標高(&水深)のデータは経度 180° で切れているので,180° 線をまたぐ処理が必要になります.

考え方

オリジナルの地形データの右側に,同じデータセットをもう一つくっつけます.これが作れれば,連続した描画範囲を指定できるので,前回と同様の方法で描画できるはずです.

画像1

コード

最初の方は前回と同じです.まずはライブラリのインポートから.

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()

出力はこんな感じになります.

画像2

ライブラリのバージョン

以下のバージョンで動作確認を行っています.バージョンの確認方法については,こちらの記事も合わせてどうぞ.

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

参考

地球全体を、標高1m毎に別の色で塗り絵して地図を作ってみる。
Python で海底地形図を描く (1)

この記事が気に入ったらサポートをしてみませんか?