解析雨量を Pythonで読み、作図する (試作版)
気象庁の解析雨量 (2006年以降) のデータを、Python で読み、
Python のパッケージ「Cartopy」をつかって地図上に作図するプログラムを作成した。
結果の画像の例 (2022年9月18日 日本時間 21:00 の1時間雨量)
序論
解析雨量とはなにかについては、その開発にかかわった人の著書のうちの
解析雨量のところ (数ページ) を見てほしい。
- 牧原 康隆, 2020: 『気象防災の知識と実践』。朝倉書店。
[読書メモ]
気象庁の解析雨量のデータ形式は、2005年のぶんまでは独自の形式、
2006年からはWMO (世界気象機関) できめられた GRIB 第2版 (GRIB2 ともいわれる) である。
解析雨量の値は、本来、1時間降水量 (mm/h) である。
(30分間隔のばあいも、その時刻までの1時間の降水量である。)
しかし、有効数字があまり多くないので、1バイト (8ビット) で表現できる 0〜255のレベルでおきかえられ、
さらに、となりどうしの画素の値が同じであるばあいに連長 (run length) 法によって圧縮されている。
わたしはこれまで解析雨量のデータを GrADS というソフトウェアであつかってきたが、
GrADSで読むためには、連長圧縮をといた中間ファイルをつくらなければならない。
ウェブ検索したら、GRIB2 形式の解析雨量データを Python で中間ファイルをつくらずに読んで作図するプログラム例が
みつかったので、それにならってやってみることにした。
[補足] このウェブページをつくりはじめるまえの作業の覚え書き
データを地図上に作図するプログラム
ここでのプログラム例文は、GRIB2 形式が採用された 2006 年以後のデータを対象としている。
2005年以前の気象庁独自形式には対応していない。
GRIB2形式のデータでは、日時は世界時 (UTC) であたえられているが、
利用者は日時を日本標準時 (JST) であたえられるようにし、
図には UTC と JST の両方の日時を表示するようにした。
画像ファイル出力のファイル名には JST をつかい、それを明示するようにした。
(Python の datetime パッケージの機能をつかった。)
【[2023-09-12 補足] JST をつかわずに UTC のままあつかいたいばあいもある。
その例として、UTC の 1日 (0時から24時まで) の降水量を積算して表示するプログラム例を追加した。】
数値レベルと図の色わけとの関係は、draw_a_contour_map 関数の中の levels と colors という
リスト変数のなかみで指定している。現在あたえている色のくぎりの数値レベルはつぎのとおり。
(あまり深く考えたものではなく、ためしに表示した値がおさまるようにした。)
- hourly, halfhourly では、0, 0.5, 1.5, 2.5, 5, 10, 20, 40, 80 [mm/hr]。
(データから出てくる値は、0 [無降水], 0.4 [微量の降水], 1 以上の整数値である。)
- 3hourly では、0, 0.5, 5, 10, 20, 40, 80, 160, 320 [mm/3hr]。
- daily では、 0, 0.5, 10, 25, 50, 100, 200, 400, 800 [mm/day]。
解析雨量の値が欠損になっているばあいは、関数 load_jmara_grib2 の中であたえる値は
(@vcpf さんのプログラムをひきついだので) -10 である。
メインプログラムで、降水量の値を 0.1 mm 単位から mm 単位になおすために 10 でわっているので、
欠損をあらわす値は -1 になる。表示の際はこれを考慮に入れる。
(いまの配色では、値が -1 (欠損) ならば薄い赤、0 (無降水) ならば白、0.4 (微量の降水) ならば薄い黄色としてある。)
【[2023-06-30 補足] matplotlib のバージョンによって、contourf の set_under, set_over が無効になってしまう
ばあいがある。そのときは、levels の内容を上下にひろげ、colors に set_under, set_over で指定する色もふくめておく。
ku_map_3hourly.py では、そのばあいの変更点を「#*」ではじまるコメント行でふくめてある。】
図に日本語文字を表示するところで、MS Windows の日本語フォントを指定している。
他の OS では変更が必要だろう。
つかいかた
この例文を各自の目的にやくだてるには、
プログラムを編集して、文字を書きかえる必要がある。
(日時の情報は、人がキーボードからうちこむようにしたり、
別のファイルから読みこむようにすることもできるが、
そのための書きかえかたの説明は、まだ用意していない。)
計算機上でデータのある場所を指定する。
- 入力データとして、年をあらわす4けたの数字のディレクトリの下に、
その年のデータをふくむ気象業務支援センターのDVDディスクの内容がそのディレクトリ構成で置かれていることを前提とする。
(ただし、わたしの所属グループのローカルなデータ保存状況を反映して、
2013年から2019年までのデータは、年をあらわす4けたの数字のディレクトリの下に、
直接データファイルが置かれていることを前提としている。
このディレクトリ構成のちがいには、プログラムで対処している。)
- 例文プログラムの「##### Main program #####」というコメントのあとにある
「datapath」という変数に代入する値を、自分がデータを置いたディレクトリ名に変える。
- [2023-06-30 補足] hamana では 「datapath = '/raid/kaiseki_uryo/'」とする。
- [注] ディレクトリ区切り文字は、Windows ではふつうバックスラッシュか円マーク ('\') だが、
ここでは Linux と共通にするためにスラッシュ ('/') にしている。このほうがよいと思う。
- [注] ディレクトリ名の文字は、ASCIIの1バイト文字コードで
表現できるアルファベット・数字・算数記号にかぎったほうがよいだろう。
日本語の漢字・かななどがあるとうまくいかないかもしれない。
対象日時を指定する。
- hourly のプログラムのばあい
- 「#1時刻作図」「#24時刻作図」「#48時刻作図」と注記した行の
いずれかを有効にし、その他のものをコメントにする。
(複数行つづくときは、「ここから #48時刻作図」「ここまで #48時刻作図」
のようなコメントではさんでしめしてあることもある。)
- 「#24時刻作図」: 日本時間 (JST) の1時から24時までの毎時の図を作図する。
(気象庁の地上観測の雨量データの期間のとりかたにあわせた。)
ただし、24時の図のファイル名や画面表示は、翌日の 0時となる。
- 「#48時刻作図」: 日本時間 (JST) の0時30分から24時までの30分ごとの図を作図する。
ただし、24時の図のファイル名や画面表示は、翌日の 0時となる。
- 「#1時刻作図]: 1つの時刻についての図を作図する。
- iyearjst, imonjst, idayjst の値を、自分が表示させたい日付にする。
- 「#1時刻作図] のばあいは、「for ihourjst in [14]」と
「for iminjst in [0]」の角かっこの中の数値を、自分が表示させたい時刻にする。
- daily, 3hourly のプログラムのばあい
- iyearjst, imonjst, idayjststart, idayjstend の値を、自分が表示させたい日付にする。
(idayjststart と idayjstend の値を同じにすれば, daily では 1枚、3hourly では 8枚の図がつくられる。
作図したい期間が月にまたがるばあいは月の切れめでわけて
複数回プログラムを実行することを想定した。)
図形の表示さきを指定する。
- 「draw_a_coutour_map」関数の最後ちかくの
「#画面表示」とマークしてある「plt.show()」を生かすと、画面に図が出る。
- 「画像保存 & 終了」のコメントのあとの「#画像ファイル」とマークしてある4行を生かすと、png 画像ファイルができる。
対象地域を指定する。
- 「#全国」「#本土」「#関東」「#狭関東」と注記した行の
いずれかを有効にし、その他のものをコメントにする。
(複数行つづくときは、「ここから #全国 以外」「ここまで #全国 以外」
のようなコメントではさんでしめしてあることもある。)
- 「#全国」: 解析雨量の全領域 (北緯 20～48度、東経 118～150度) を
作図する。地図の範囲は 北緯 20～50度、東経 115～150度 にしてある。
- 「#本土」: 沖縄、奄美、小笠原諸島をふくまないがそのほかの日本の陸地のほとんどをふくむ領域
(北緯 30～46度、東経 128～146度) を作図する。
- 「#関東」: 関東甲信越をふくむ つぎの緯度経度範囲
(北緯 34.5～38.5度、東経 136.5～141.5度) を作図する。
- 「#狭関東」: 関東をふくむ つぎの緯度経度範囲
(北緯 34.5～37.5度、東経 138.5～141.5度) を作図する。
- [注] あと知恵ではこれを「関東」とするべきだったのだが、さきに「#関東」を設定してしまったので、
くるしまぎれにこのような名まえにした。
- その他の地域を抽出することは、「#本土」「#関東」「#狭関東」と注記した行のいずれかをもとに数値を書きかえればできると思う。
改訂履歴
- [2021-04-27] このウェブページを作成し、1時間ごとの値を作図するプログラム例を置いた。
- [2021-06-02] 1日ごとの値を作図するプログラム例を追加した。
また、1時間ごとの値を作図するプログラム例を改訂した
(実質の処理はかわらず、デバッグと、プログラムの見とおしをよくするための改訂)。
- [2021-06-02] 日あるいは月の切れめが世界時と日本標準時でずれるところの処理で失敗していたところを修正した。
- [2021-06-05] 2020年のデータを処理できるようにした。
- [2021-07-01] 「#関東」というオプションで指定される緯度経度範囲は関東地方だけでなく中部地方の大部分をふくむものだったので、
もっと関東にしぼりこんだ領域を指定するオプション「#狭関東」を追加した。
- [2021-07-29]
- 1時間ごとの値をかく機能と30分ごとの値をかく機能を hourly のプログラムに統合した。
- 都道府県境界をかくようにした。 (ただし #全国 では かかない。)
- draw_a_countour_map に部分配列をひきわたすようにし、計算機の負荷をへらした。
- daily のプログラムでは、1時間値から1日値への集計も、表示対象となる格子点に
ついてだけおこなうようにし、計算機の負荷をへらした。
ただし、配列の大きさは全国のデータがはいるようにとっており、その部分をつかっている。
- カラーバーの数値の単位の文字列「mm/hr」「mm/day」の位置を適切にした。
- 日本標準時で2006年初、世界時で2005年末となる期間をあつかえるようにした。
(2006年のディスクの2005というディレクトリにデータがあるので、そこを読みにいく。)
- dailyのプログラムの1時間値から1日値への集計での欠損値のあつかいを改訂した。
- [2022-11-07] 都道府県境を作図するところを
[日本の都道府県境をかく / 中国の「省」の境をかく] にのべたように修正した。
- [2023-02-18] 地図の範囲の選択に「#本土」を追加した。結果の画像の例を変えた。
- [2023-06-30] 3hourly を追加した。
1時間ずつずれた3時間ごとの積算値を作図する [2023-06-30 追加]
- 日時の指定を、Pandas の date_range で日付のリストを生成して、それをつかって for でループをまわす方法によることにした。
例文ではつぎのようになっている。
- datetime_list = pd.date_range(start='2022-09-23 09:00:00', periods=26, freq='H')
- ユーザーは、start に与える値を実際に計算させたい期間の最初の時刻にすればよい。
period に与えた個数の1時間値を読み込み、それよりも2個すくない3時間積算値を作図する。
- プログラム内部では、読みこんだうちでいちばん新しい雨量の値を 2次元配列 rain に、
2ステップまえの値を rainbakbak、1ステップ前の値を rainbak に代入してもつようにしている。
(ひとつの3次元配列にまとめたほうがきれいなプログラムになると思うが、まだ成功していない。)
時系列書き出しプログラム
[2021-06-03 追加, 2021-06-05 改訂, 2021-07-30 改訂]
解析雨量のGRIB2形式データを読みこむところは共通で、作図するのではなく、
地点をえらんで、その位置の降水量の値の時系列をテキストファイルに書き出すプログラムを試作している。
こちらは、Cartopy も matplotlib もつかわない。
まず、指定された緯度・経度にいちばん近い ます目 の降水量の1時間ごとの値を書き出すプログラムをつくった。
月単位の期間を指定して、その全部の日の全部の時刻の値を書き出すことにした。
ただし、年にまたがる期間は指定できない (それがほしいばあいは複数回にわけて実行する)。
つかいかた
主プログラムのはじめ (スクリプトの中では途中になるが、「##### Main Program #####」というコメントのはいった直後) の、
つぎの変数の値をあたえている代入文を、ほしい位置の緯度・経度に書きかえる。
- latdeg, latmin, latsec ... 緯度 (北緯)。度、分、秒であたえる形になっているが、整数にかぎることはない。
「度に小数部のある値をあたえて分、秒はゼロにする」
「度に整数、分に小数部のある値をあたえて秒はゼロにする」などのあたえかたもできる。
- londeg, lonmin, lonsec ... 経度 (東経)。緯度と同様。
すこしさきの「# 対象日付」というコメントのあとのつぎの変数の値をあたえている代入文を、
ほしい日付に書きかえる。
- iyearjst ... 年 (西暦)
- imonjststart, imonjstend ... はじめとおわりの月。1から12の整数値であたえる。
はじめとおわりを同じ値にすると 1か月ぶんの出力になる。
「datapath」という変数に代入する値を、自分がデータを置いたディレクトリ名に変える。
スクリプトを実行する。(Anaconda Powershell Propmpt ならば、
コマンドラインで「python ku_timeser_hourly_v20210602.py」)
出力データの形
つぎのようなテキストファイルができる。
- ファイル名は、たとえば「rrhourly201906_09_lat36_158_lon139_375.txt」のようにつけられる。
- 「rr」は雨量によくつけられる略号 (本来は rainfall rate 「降水強度」だったのかもしれない)
- 「hourly」は1時間値であること
- 「2019」が指定した年、「06」と「09」が指定した はじめとおわりの月
- 「lat36_158」「lon139_375」は、ぬきだされた画素の中央の位置の緯度・経度。
度の小数点3けたまでで表現し、小数点を下線("_"、アンダースコア)でおきかえた。
指定した緯度・経度に近いがそれと同じではないことに注意。
- ファイルの内容は、空白くぎり (ただし連続した空白は1つの空白と同じとみなす) でもあり、
フィールド値固定でもある、テキストファイルにした。
表計算ソフトウェアにそのまま読みこめると思う。
(ファイルのあたまにデータが何であるかをしめすコメント行を入れたかったのだが、
そうすると表計算ソフトウェアでうまく読みこめないので、それをあきらめて、
最小限の情報をファイル名に入れた。)
- 左から5列めまで (Python流にいえば第0列から第4列、表計算流にいえばA列からE列)は、
年、月、日、時、分 である。
気象庁の地上気象観測データにあわせて、日本時間で、「0時0分」でなく「24時0分」をふくむ形にした。
- (世界時 (UTC) の日時を書き出すこともできるが、いくつもの変更が必要で、
オプションとして用意できていない。)
- (「分」の値は、このプログラムでは「00」にきまっているので、省略してもよいのだが、
30分ごとの値を書きだすように改造するばあいにそなえて、列をつくっておいた。)
- 6列め (Python流にいえば第5列、表計算流にいえばF列) に、雨量の値 (単位は mm/時) が
小数点下1けたの形で書かれる。
データ欠損のところは「-1.0」と表現されるはずであるが、まだ確認していない。
日降水量の時系列書き出しプログラム
[2021-07-30] 時別値から日別値を集計するところでの欠損値のあつかいを改訂した。
指定された緯度・経度にいちばん近い ます目 の降水量の1日ごとの値を書き出すプログラムをつくった。
日本時間の0時から翌日0時までの集計値である。
欠損値は単純に集計に入れない。
(もし観測されていたとしたら、降水量はもっと大きかった可能性がある。)
集計につかわれた有効な観測値の個数を書き出しておく。
作業手順は1時間ごとの値を書き出すプログラムと同様である。
結果のファイル名も同様で、hourly のかわりに daily である。
結果のファイルの列の意味はつぎのとおり。
- 左から3列めまで (Python流にいえば第0列から第2列、表計算流にいえばA列からC列)は、
年、月、日である。
- 左から4列め (Python流にいえば第3列、表計算流にいえばD列) は、日降水量 (mm/日) である。
- 左から5列め (Python流にいえば第4列、表計算流にいえばE列) は、
集計につかわれた有効な観測値の個数である。観測がとまらず正常におこなわれていれば 24 となる。
時間降水量の領域平均値の時系列書き出しプログラム
[2021-06-18 追加, 2021-07-30 改訂]
緯度・経度座標の長方形領域で平均した降水量の1時間ごとの値を書き出すプログラムをつくった。
ひとまず バージョン1 を出しておく。今後、改訂する。
つぎのような不満がのこるプログラムである。
- 観測値が部分的に欠損のばあい。
【これは、満足な解決方法はなく、妥協が必要である。】
このバージョンでは、
長方形領域の画素数の半分以上に有効な値があれば、有効な画素についての平均を書き出し、
半分未満ならば、結果を欠損とすることにした。
ただし、ここで「半分」としたところは、thresh という変数にしてあり、0.5 という値を与えてある。
なお、全画素のうち有効な画素のわりあいも書き出すことにした。すべて有効ならば 1.0 となる。
- 画素が代表する地球上の面積は緯度によってちがうが、
このバージョンでは、画素の面積は一定であるかのように、画素の個数で平均をとっている。
【近似を一段高めるには、地球の形を球で、画素の形を台形で近似して、画素の面積は 緯度の cos に比例するとするのがよさそうだ。】
- 領域平均を計算するところは、二重ループの中に if の判断がはいった形になっており、
実行に時間がかかる。
【Python の numpy の配列演算をうまくつかえばもっと効率よくできるだろう。】
つかいかた
次のところを書きなおす。
- datapath という変数に与える文字列を、
作業するパソコンで解析雨量データを置いたディレクトリ名にする。
- 対象日時。日本標準時で指定する。
iyearjst に年、imonjststart と imonjstend に はじめとおわりの月を指定すると、
その期間の全部の時刻について値を書き出すようにしてある。
- thresh (上に説明した) を変えたければ変える。
- 領域の北端の緯度を latdeg0, latmin0, latsec0 で指定する。
- 領域の南端の緯度を latdeg1, latmin1, latsec1 で指定する。
- 領域の西端の経度を londeg0, lonmin0, lonsec0 で指定する。
- 領域の東端の経度を londeg1, lonmin1, lonsec1 で指定する。
スクリプトを実行する。(Anaconda Powershell Propmpt ならば、
コマンドラインで「python ku_areamean_hourly_v1.py」 [2021-06-19 訂正])
出力データの形
結果のファイル名は、たとえば
「rrhourly202006_09_lat36_250_35_750_lon138_750_139_750.txt」
のようにつけられる。
これは、2020年 6月から9月まで、
緯度 36.25度と35.75度の間、経度 138.75度と139.75度の間の値である。
結果のファイルの列の意味はつぎのとおり。
- 左から3列めまで (Python流にいえば第0列から第2列、表計算流にいえばA列からC列)は、
年、月、日である。
- 左から4列め (Python流にいえば第3列、表計算流にいえばD列) は、領域平均の時間降水量 (mm/時) である。
- 左から5列め (Python流にいえば第4列、表計算流にいえばE列) は、
集計につかわれた有効な画素数の領域の全画素数のうちでのわりあいである。
全画素が有効ならば 1.0 となる。
今後の作業についての構想
複数の地点の値がほしいとき、実行能率の面では、毎時の全国のデータを読むごとに
複数の地点の値を書き出すように、プログラムを改造したほうがよい。
ただし、そうするとプログラムがいくらか複雑になり、プログラム保守の面では、うまくないかもしれない。
領域平均を書き出すプログラムには上にのべたように不満があるので、改訂したい。
実行環境
ここにあげたプログラム例は、MS Windows 日本語版で実行することを想定し、Windows 10 上でつくっている。
Windows での実行は、Anaconda をインストールし、その中で Cartopy をいれたうえで、
Anaconda Powershell Prompt から「 python なになに.py 」で起動するか、
Anaconda といっしょにインストールされる Spyder に読みこんで実行するのがよさそうだ。
(わたしは Anaconda Powershell Prompt から実行している。)
[Windows で Pythonをつかう準備 (Anaconda, cartopy のインストール) と 実行]
標高の等高線を重ねる (2023-04-28 追加)
降水量と地形の関係をみるため、降水量と標高を同じ図に表示したいことがある。
ここでは、降水量の段彩 (contourf() による) と、標高の等高線 (contour() による) をかさねてみることにした。
例としてとった対象領域は、ここで「#関東」とした緯度経度範囲 (関東地方と中部地方の大部分をふくむ) である。
標高データは [地形 (標高) データを Pythonで読み、作図する (試作版)] であつかったのと同じ、
アメリカ合衆国の NOAA NGDC (現 NCEI) で編集された
世界の標高・水深データセット ETOPO1 からこの緯度経度範囲を切り出したものである。
【国土数値情報の「1 km メッシュ」もためしたのだが、この範囲の図にしめすにはそのままではこまかすぎ、
平滑化または間引きが必要なことがわかったので、ひとまず見送っている。】
プログラム例では、等高線を、標高 500 m 以上について、500 m 間隔にひいた。
海岸線と県境を黒の線でいれているので、等高線は赤にしてみた。
このプログラムでは、メインプログラムで標高データを配列 topo によみこむとともにその格子点の経度と緯度の配列 lontopo, lattopoを
用意し、この 3 つの配列を draw_a_coutour_map() にひきわたして、draw_a_contour_map() のなかの降水量の段彩のあとで
標高の等値線をひいている。
【なお、図の大きさを figure(figsize=(12,8)) (横 12インチ、縦 8インチの紙に相当) とし、それにあわせて あちこちの fontsize を調整した。
ただし、このウェブページにおさめる際には width を指定して縮小表示している。】
参考にした情報源
【ウェブ上でキーワード検索をかけたり、それで見つかったものからいもづる式にたどったりして、
ブログ記事などかなり多くの記事をながめたうえで読むものをえらんだのだが、
qiita.com のページを参照したものがわりあい多い。
このサイトはたしかに有用な記事があるが、このサイトの記事がみんな信頼できるわけでもない。】
解析雨量 (GRIB2形式) のデータの解読は、つぎのウェブページのプログラム例によった。
https://qiita.com/vpcf/items/b680f504cfe8b6a64222#%E8%A7%A3%E6%9E%90%E9%9B%A8%E9%87%8F%E3%82%92python%E3%81%A7%E8%AA%AD%E3%81%BF%E3%81%9F%E3%81%84
@vcpfが2020年07月19日に更新
気象庁解析雨量(GRIB2形式)をPythonで扱う
これは作図に、matplotlib にくわえて Basemap というパッケージをつかっている。
しかし、Basemap は開発終了し、その開発者は Cartopy への移行を勧めている。
Cartopy を利用した作図は、つぎのウェブページのプログラム例を参考にして、前のプログラム例を書きかえた。
https://qiita.com/wm-ytakano/items/1e6a5d250fe0d079b08e
@wm-ytakanoが2021年02月14日に更新
PythonでGSMのGPVデータから気温の水平分布を描画する
都道府県境界を入れるところは、つぎのウェブページのプログラム例
(中国の「省」「自治区」の境界を入れるものだが)を参考にして、
Natural Earth というウェブサイトにある行政区境界の地理データをつかっている。
(ただし、2022-11-07 の改訂で、Cartopy の機能で自動的に読みこむのをやめ、
あらかじめダウンロードして編集した境界データをつかうことにした。)
https://geocat-examples.readthedocs.io/en/latest/gallery/Overlays/NCL_overlay_11b.html
GeoCAT-examples
NCL_overlay_11b.py
(C) Copyright 2021, University Corporation for Atmospheric Research. Revision 4e24dbc1.
2021-04-27 .. 2021-06-05, 2021-06-18, 2021-06-19, 2021-07-01, 2021-07-30, 2022-11-07, 2023-02-18, 2023-04-28, 2023-06-30, 2023-09-12
増田 耕一 (MASUDA Kooiti)