皆さんこんにちは。
前回はnetCDFの構造と参照の仕方について簡単に説明を行いました。
今日はデータの詳細と図化をしていきたいと思います。

今見ているブログページを上スクロールしてもらうと
「netCDF内部構造参照結果」というスクショ画像があります。
（※第80回本文）

これを見ると記載されていることですが、次元サイズはｘが21601で、これは-180度から180度までの範囲の東経を経度としています。
yは次元サイズ10801で-90度から90度までの北緯を緯度としています。
これらの緯度(ｙ)経度(ｘ)を格子状で区切った時、1つのセルに当てはめられる整数型の値がｚで、つまり標高(ｍ)ですね(^^)

例として世界マップから日本の位置を抽出して標高マップを表示させてみます。
スクリプトコードは既に個人ブログ等で参考となる資料があったので、今回はmemomemokunさんのコードを基本としてアレンジを加えました。

そう言えば、プログラムコードに関する著作権に関しては丸パクリはダメだけど、元のコードを参考に新たなプログラムを作ったり組み替える分には認められるようです。
ただ実際法律は曖昧でその区分も個人か企業かなど時と場合で変わるでしょうから、トラブルを避けるため最初に規約として提示したほうがいいみたいですね。

ちなみに本ブログのコードは個人・企業問わず誰でも研究開発や趣味、勉強に使用することができます。
（写経をしないとコーディングは上手くなりません！）
商用利用については、公開しているプログラムコードと全く同じものを販売目的で使用するのはご遠慮くだい。
その他は特に規制ありません。

では全コードを表示します。

1～19行目まではnetCDFからデータセットの読み込みを行っています。
そして20行目からは日本の座標位置を取得しています。

#japan:東経127-150度,北緯25-47度
lower_left_x=np.argmin(np.abs(x-127))
lower_left_y=np.argmin(np.abs(y-25))
height=np.argminnp.abs(y-47))
width=np.argmin(np.abs(x-150))

X=X[lower_left_y:height,lower_left_x:width]
Y=Y[lower_left_y:height,lower_left_x:width]
Z=Z[lower_left_y:height,lower_left_x:width]

コード内、小文字のｘ、ｙはそれぞれ-180~180と-90~90の範囲を１次元のfloat型で表しています。
つまりこれは東経と北緯が何度なのか、地図上の位置を示します。
大文字のX、Yは、matplotlib によって実際に標高Zをマッピングさせるため、np.meshgrid (x,y)で与えられた2つの配列から2次元の方眼線を作成します。
今回日本の東経は約127度～150度。北緯は25度～47度で切り取ることにしました。
切り取りの方法はnp.argmin()を使って、四角形の4カ所のポイントの東経北緯に等しいインデックスを抽出します。
後はマッピング用の変数X、Y、Zから日本の座標位置を切り出します。

各標高に対する配色については、元のコードと地理院地図を参考にして少し書き直しました。

# 配色設定
cmap = colors.ListedColormap(['blue', 'royalblue', 'cornflowerblue', 'skyblue', 'powderblue', 'g', 'khaki', 'wheat', 'tan', 'peru', 'chocolate', 'salmon', 'indianred', 'brown', 'snow'])
bounds = [-8000, -4000, -2000, -1000, -500, 0, 100, 300, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]

その結果がこちら・・

上手く切り取れました！

参考までにETOPO1の他にもNOAAにはGSHHGという世界の海岸線や湖沼を含む高解像度地形のベクトルデータも公開されており、こちらも自由にダウンロードすることができます。
ただし、ダウンロードしたgshhg-shp-2.3.7データに関してはモトハシも普段目にしない圧縮形式やファイル名だらけで、どうやらCADで扱うタイプのファイル？なのか詳しいことはよく分かっていません。
(ダムとか湖の解析に役立ちそうなのに・・・)

まだまだ勉強することがたくさんあるなあ～と痛感しましたが、ゆるゆる頑張っていこうと思います。
ではでは！