国立研究開発法人防災科学技術研究所 水・土砂防災研究部門
国立研究開発法人防災科学技術研究所 水・土砂防災研究部門
トップ 一覧 検索 ヘルプ RSS ログイン

国土地理院地図の上に1変数と風ベクトル水平断面図を書く(xarray)

国土地理院の地図(タイル)を取得し、その上に分布を半透明にさせて作図する

ver5.tar.gz(16)
ver6.tar.gz(14)
ver7.tar.gz(21)

ver5はver4の上位互換である。またver4のDHDL/subprog.pyのcheck_variableサブルーチンで
水平断面だけを書く場合、コマンド引数X0_VCSが定義されていないというエラーが出てかけなかった問題を解決した。

 追加機能(ver4からver5)

指定した緯度経度の、指定したzoom-levelでの地図を国土地理院から画像タイルを複数ダウンロードして、
指定した緯度経度分を合成したり、切り出したりするgsitile.pyを作成し、ベースとなるPNGを作図し、その上に
雨等の分布を示すことができるようになった。
mainプログラムはmain_H1UV_GSItile.py です。1枚の水平断面(?変数と風を記載可能)を国土地理院の地図の上に重ねて書きます。

 追加機能(ver5からver6)

ver6はver5の上位互換。共通化できるプログラムをsubprogsに含めました。また入力のNETCDFデータが
これまでは4次元(TIME,LEV,LAT,LON)を想定していたのに対して、2次元以上であればどんな入力でも対応できるようになりました。
例えば、TOPO(LAT, LON)のような高度と時間の次元を持たないものや、dual解析のような3次元(LEV,LAT,LON)のものであったり、
地表面の時間変化SU(TIME,LAT,LON)のような変数も統一して扱えるようにしました。また、get1var.pyはget1varUV.pyで統一することができ、
get1var.pyを廃止しました。基本的な使い方は全く同じです。ver6の方が様々な入力データに対応できるようになったということです。
また、setcolors.pyでは暖色系の色のみの色調もCT=3で選択できるようになりました。

 追加機能(ver6からver7)

アニメーションを作る際の国土地理院サーバへの負荷を減らす工夫をした。スナップショットをたくさんつくり、

convert -loop 0 - delay *.png anime.gif

のコマンドで最終的にGIFアニメを作るときには、main_H1V1UV_GSItile.pyを何回も実施することになる。その際に
バックグラウンドになる国土地理院地図はかならず同じ領域だと思うので、2枚目以降、ダウンロードと作成が不要になる。
-ANM をオプションをつけると、1回目にbackground.pngを背景地図として作成し、2回目以降はこのbackgroud.pngをコピーして
使うようになる。また、-ANM をつけると、plt.show()を省略することで、図を表示せずに図の作成を淡々とやってくれる。
animateで図を作る場合、PNGで作るのではなく、-Oで指定するファイルはJPGがオススメ。pngcrushなどで最適化する等
しないと、matolotlibでつくるPNGファイルはdisplayやeogコマンドで表示しようとすると大変遅くなる。pngcrush -brute 元.png 新.pngで
ファイルを最適化するとようやく、displayコマンドなどでスピーディーに表示ができるようになる。
また、国土地理院地図を示す相手はおそらく、日本人だろう。時刻をUTCで表示すると、難しくなってしまうため、-JPN オプションを
追加した。-JPN をつけると、時間をJSTで表示してくれる。また、-ARオプション(凡例を自動化しない)と-VTオプション(vertically interpolation)の
オプションは以前は、なにか数字を入れないとエラーになっていたが、このバージョンから単に-ARや -VT をつけるだけでよい。
同様に単に"-XXX"のようにオプションだけつけておけばいいオプションは、-AMNと-JPN である。 -AMNは国土地理院地図を表示するときだけに
使われる機能である。-JPNはすべてのmainプログラムで使用が可能となっている。

  使用するライブラリー

pillowを利用しますので、conda install pillow でインストールが必要

 使い方 (ver5 ver6)

pytyon main_H1V1UV_GSItile.py NETCDFファイル名 変数名 
-W 風オプション(0: ベクトル 1:矢羽。鉛直断面と水平断面を両方作図するこのバージョンではベクトルのみ対応)
-WT 風間引きオプション(整数を設定すると、その整数分だけデータをスキップして間引く)
-CT カラーパレット番号を入れることで色パターンを設定します. 0:青から赤へ 1: 気象庁降雨分布 2:白から青へ という3パターンあります
-AR 0 AutoRangeを行わないオプション(-AR を記載しないと自動的に凡例の範囲を検索する)
-O 図をファイルとして保存するオプション。保存ファイル名を指定する。 今回はPNGしか受け付けない。
-T 表示する時刻を指定(大文字のTを指定すると時刻をYYYY-MM-DD-HH:MM:SSのフォーマットで指定)
-t 表示する時刻を指定(小文字のTを指定するとnetcdfの中での初期時刻からのスキップする時刻数を整数で指定)
   時刻を指定しない場合は、インタラクティブに時刻をpytyonが聞いてくれるので、入力する
-z 水平断面を表示する場合の層番号を指定。
-Z 水平断面を表示する、高度をメートル単位で指定
-ZL zoom level 10くらいがちょうどいい。13になると地図が細かすぎる。
-TL 分布の透過度。0が完全に透過(地図だけが見える)1が完全に不透明(地図がまったく見えない)。0.3くらいがいい。

 使い方 (ver7)

地理院地図で表示する場合、

pytyon main_H1V1UV_GSItile.py NETCDFファイル名 変数名 のあとに 
-W 風オプション(0: ベクトル 1:矢羽。鉛直断面と水平断面を両方作図するこのバージョンではベクトルのみ対応)
-WT 風間引きオプション(整数を設定すると、その整数分だけデータをスキップして間引く)
-CT カラーパレット番号を入れることで色パターンを設定します. 0:青から赤へ 1: 気象庁降雨分布 2:白から青へ 3:白から赤へ という4パターンあります
-AR AutoRangeを行わないオプション。つまり手動で凡例の値を決める必要がある。DHDLのsubprogs.pyの中で設定しています(-AR を記載しないと自動的に凡例の範囲を検索する)
-O 図をファイルとして保存するオプション。保存ファイル名を指定する。 今回はPNGとJPGを想定。
-T 表示する時刻を指定(大文字のTを指定すると時刻をYYYY-MM-DD-HH:MM:SSのフォーマットで指定)
-t 表示する時刻を指定(小文字のTを指定するとnetcdfの中での初期時刻からのスキップする時刻数を整数で指定)
   時刻を指定しない場合は、インタラクティブに時刻をpytyonが聞いてくれるので、入力する
-z 水平断面を表示する場合の層番号を指定。
-Z 水平断面を表示する、高度をメートル単位で指定
-ZL zoom level 10くらいがちょうどいい。13になると地図が細かすぎる。
-TL 分布の透過度。0が完全に透過(地図だけが見える)1が完全に不透明(地図がまったく見えない)。0.3くらいがいい。
-AMN アニメション支援オプション。図を表示せず、作成のみ、背景の国土地理院地図は最初の一回だけ作成
-JPN 時刻表記をJSTとする。

水平断面図では

pytyon main_H1V1UV.py NETCDFファイル名 変数名 のあとに 
-W 風オプション(0: ベクトル このバージョンではベクトルのみ対応)
-WT 風間引きオプション(整数を設定すると、その整数分だけデータをスキップして間引く)
-CT カラーパレット番号を入れることで色パターンを設定します. 0:青から赤へ 1: 気象庁降雨分布 2:白から青へ 3:白から赤へ という4パターンあります
-AR AutoRangeを行わないオプション。つまり手動で凡例の値を決める必要がある。DHDLのsubprogs.pyの中で設定しています(-AR を記載しないと自動的に凡例の範囲を検索する)
-O 図をファイルとして保存するオプション。保存ファイル名を指定する。 今回はPNGとJPGを想定。
-T 表示する時刻を指定(大文字のTを指定すると時刻をYYYY-MM-DD-HH:MM:SSのフォーマットで指定)
-t 表示する時刻を指定(小文字のTを指定するとnetcdfの中での初期時刻からのスキップする時刻数を整数で指定)
   時刻を指定しない場合は、インタラクティブに時刻をpytyonが聞いてくれるので、入力する
-z 水平断面を表示する場合の層番号を指定。
-Z 水平断面を表示する、高度をメートル単位で指定
-AMN アニメション支援オプション。図を表示せず、作成のみ、背景の国土地理院地図は最初の一回だけ作成
-JPN 時刻表記をJSTとする。

鉛直断面図では

pytyon main_V1.py NETCDFファイル名 変数名 のあとに 
-W 風オプション(0: ベクトル 1:矢羽。鉛直断面と水平断面を両方作図するこのバージョンではベクトルのみ対応)
-WT 風間引きオプション(整数を設定すると、その整数分だけデータをスキップして間引く)
-CT カラーパレット番号を入れることで色パターンを設定します. 0:青から赤へ 1: 気象庁降雨分布 2:白から青へ 3:白から赤へ という4パターンあります
-AR AutoRangeを行わないオプション。つまり手動で凡例の値を決める必要がある。DHDLのsubprogs.pyの中で設定しています(-AR を記載しないと自動的に凡例の範囲を検索する)
-O 図をファイルとして保存するオプション。保存ファイル名を指定する。 今回はPNGとJPGを想定。
-T 表示する時刻を指定(大文字のTを指定すると時刻をYYYY-MM-DD-HH:MM:SSのフォーマットで指定)
-t 表示する時刻を指定(小文字のTを指定するとnetcdfの中での初期時刻からのスキップする時刻数を整数で指定)
   時刻を指定しない場合は、インタラクティブに時刻をpytyonが聞いてくれるので、入力する
-X0V 鉛直断面の始点の経度
-Y0V 鉛直断面の始点の緯度
-L0V 鉛直断面の始点からの距離(km)
-A0V 始点から真東を0度とし、反時計回りに回転した方向に鉛直断面を作成する。90度とすると北向き、180度とすると西向き、-90度で南向き
-TP  地形マスクオプション。inputファイルにTOPOという変数で標高が入力されていれば、数字の0を入れればよい。外部ファイルとして-TP ファイル名として同じ座標系の標高のnetcdfファイルを入力できる
-VT  鉛直を線形内挿して、一定間隔の高度とした高度座標系で鉛直断面を表示。-TPと一緒に使うと、地形をマスクしてくれる。気圧座標系の場合、高度座標系に変換する
-JPN 時刻表記をJSTとする。

水平・鉛直断面図では

pytyon main_H1V1UV.py NETCDFファイル名 変数名 のあとに 
-W 風オプション(0: ベクトル 1:矢羽。鉛直断面と水平断面を両方作図するこのバージョンではベクトルのみ対応)
-WT 風間引きオプション(整数を設定すると、その整数分だけデータをスキップして間引く)
-CT カラーパレット番号を入れることで色パターンを設定します. 0:青から赤へ 1: 気象庁降雨分布 2:白から青へ 3:白から赤へ という4パターンあります
-AR AutoRangeを行わないオプション。つまり手動で凡例の値を決める必要がある。DHDLのsubprogs.pyの中で設定しています(-AR を記載しないと自動的に凡例の範囲を検索する)
-O 図をファイルとして保存するオプション。保存ファイル名を指定する。 今回はPNGとJPGを想定。
-T 表示する時刻を指定(大文字のTを指定すると時刻をYYYY-MM-DD-HH:MM:SSのフォーマットで指定)
-t 表示する時刻を指定(小文字のTを指定するとnetcdfの中での初期時刻からのスキップする時刻数を整数で指定)
   時刻を指定しない場合は、インタラクティブに時刻をpytyonが聞いてくれるので、入力する
-X0V 鉛直断面の始点の経度
-Y0V 鉛直断面の始点の緯度
-L0V 鉛直断面の始点からの距離(km)
-A0V 始点から真東を0度とし、反時計回りに回転した方向に鉛直断面を作成する。90度とすると北向き、180度とすると西向き、-90度で南向き
-TP  地形マスクオプション。inputファイルにTOPOという変数で標高が入力されていれば、数字の0を入れればよい。外部ファイルとして-TP ファイル名として同じ座標系の標高のnetcdfファイルを入力できる
-VT  鉛直を線形内挿して、一定間隔の高度とした高度座標系で鉛直断面を表示。-TPと一緒に使うと、地形をマスクしてくれる。気圧座標系の場合、高度座標系に変換する
-z 水平断面を表示する場合の層番号を指定。
-Z 水平断面を表示する、高度をメートル単位で指定
-JPN 時刻表記をJSTとする。

 サンプル

富士山の風下に地形性降雨が予測された分布を国土地理院地図の上に重ねた

python main_H1UV_GSItile.py ../data/JR.dmp201203120000.nc REF3 -t 15 -z 10 -O fuji.png -CT 1 -ZL 10 -W 1 -WT 1 -TL 0.3

fuji.png

 図のサイズの与え方について

DDRW/GSI_mk1figの先頭に

FigsizeFix = 0   # 0-> Fitting for Original PNG size # 1-> Fitting for user's demand size

という変数がある。0にすると、国土地理院のPNGサイズに合わせて作図する。こっちの方がいい。
もしzoomレベルをあげていると地図が細かいので、PNGサイズが大きくないと地図情報がつぶれて読めない。
もし小さい図を作りたいのであれば、細かい地図を書かずにzoom levelを10以下にすべき。その場合、
user が指定する図のサイズに合わせて書くというのも(FigsizeFix=1)ありだと思う。
この作図方法はあまり広い範囲を示すのには適していません。
hoge4.png

 謝辞

国土地理院のタイル画像を使用しています。国土地理院の地図を利用する人は以下のサイトの利用規約に従ってください。
https://maps.gsi.go.jp/development/ichiran.html
「地理院タイル (標準地図タイル)を加工して作成」などと一言謝辞を行う必要があります。

タイルの種類は以下に記載されています。サンプルではstd(標準地図)を使っています
DHDL/GSI_mk1fig.pyの19行目のtile_typeをstdからpale(淡色地図)やblank(白地図)に変更することも可能です。
https://maps.gsi.go.jp/development/ichiran.html#blank

この国土地理院タイルをダウンロードするプログラムは防災科研の前坂さんにperlプログラムを提供していただきました。
perlプログラムを元にchatgptを活用して、pythonに変換しました。両者に感謝します。