地形 (標高) データを Pythonで読み、作図する (試作版)

[解析雨量を Pythonで読み、作図する (試作版)]で 図化した雨量の分布を、地形とくらべてみるために、 地面の標高のデータを、雨量の図と同じ図法・同じ大きさの地図に表示してみようと思った。

対象地域が日本だから、日本の国土地理院の標高データがよいと思ったのだが、 つぎのような状況で、つかいにくかった。

世界の陸の標高と海底の水深をまとめたデータとしては、 アメリカ合衆国の NOAA の 旧 National Geophysical Data Center (現在は National Centers for Environmental Information の一部) がつくったものが知られている。 いま代表的なデータセットは ETOPO1 とよばれるもので、 「ETOPO1 Global Relief Model」 https://www.ngdc.noaa.gov/mgg/global/ のページにある。 氷床のところで、氷の表面を採用したものと、氷の底の岩石の表面を採用したものとの2種類の データセットをつくっている。 測地系は「WGS84」で、日本が2002年に採用した「世界測地系」と厳密に同じではないがかなり近い。 データのダウンロードは https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ からできる。 格子点の位置に値がある grid-registered と、格子線でかこまれた升目に値がある cell-registered の データがつくられているが、grid-registered のほうがもとである。 データ形式は、netCDF、geoTIFF などもえらべるが、ここでは単純バイナリにした。 grid-registered の単純バイナリのデータは、4バイト浮動小数点数、下位バイト先行 (little endian)で、 格子点数は 東西 21601、南北 10801 である。これをダウンロードした。

このうち、関東甲信越をふくむ 北緯 34.5~38.5度、東経 136.5~141.5度 の範囲を切り出した。 (切り出しは GrADS でやった [注]。まず手作業で GrADS の ctl ファイルをつくり、 GrADS の中で緯度経度範囲を指定して「set gxout fwrite」でファイルに書きだした。) もとの ETOPO1 のデータの格子点のならびは北からだったが、切り出したものは (GrADSの default の) 南からになっている。 格子点数は 東西 301、南北 241 になった。

【[注] バイナリデータとして読んでその一部分をバイナリデータとして書きだすだけなので、 Python でもできるのだが、わたしが切り出すべき部分を正しくかぞえられるか自信がなく、 GrADS で 中間結果を図示しながらやるほうがまちがえにくいと判断した。】

単純バイナリデータを Python でよむ方法は、つぎのウェブページを参考にした。 (big endian のデータむけに「>」となっているところを、little endian なので「<」に変えた。)

https://qiita.com/mhangyo/items/76db7c6a6ebba6cf4330
@mhangyoが2017年08月07日に作成
Pythonでbig endianのbinaryを読み込んでndarrayに変換する

Cartopy を利用した作図は、解析雨量のばあいと基本的に同じで、手本はつぎのウェブページである。

https://qiita.com/wm-ytakano/items/1e6a5d250fe0d079b08e
@wm-ytakanoが2021年02月14日に更新
PythonでGSMのGPVデータから気温の水平分布を描画する

プログラム例

プログラム例で読んでいるデータ (ETOPO1 から部分をきりだしたもの)

結果の画像

etopo1_around_kanto.png

(海岸線は Cartopy によって Natural Earth というサイトから取得されたデータがつかわれているが、 ここでは、解析雨量の図につかった「50 m」ではあらすぎると判断して「10 m」を採用した。)

標高 0 m をさかいに色をかえたが、ETOPO1 には陸と海の区別がないので、 海面より低い陸地が海と同様な色になり、湖が部分的に海と同様な色になって、 平野部の地図としてはちょっと気もちのよくないところがある。 山地の分布をみるにはさしつかえないだろう。

[2021-12-09 補足] 実行環境によっては、海が青でなく陸の低地と同じ緑でぬられてしまうことがある。 そこでうごいている Cartopy のバージョンか、matplotlib のバージョンでは、 mappable.cmap.set_under(), mappable.cmap.set_over() が無効になっているらしい。 そのばあいの対策として、0 m 未満、1600 m 以上の区間も levels にふくめ、 それに対応する色も colors にふくめてぬりわけるように修正したバージョンをつくった。


2021-05-05, 2021-05-11, 2021-12-09
増田 耕一 (MASUDA Kooiti)
[計算機を利用した気象分野の実習の教材] [教材目録ページ] [日本語入り口ページ] [ホームページ]
[卒論指導のページ] [3年生セミナーのページ]