気象庁の解析雨量を Pythonで読み、作図する
序論
気象庁の解析雨量のデータを、Python で読み、
Python のパッケージ「Cartopy」をつかって地図上に作図するプログラムを作成した。
解析雨量のデータは1995年のものからあるが、
ここでは、2006年以後のデータの処理にかぎって紹介する。
解析雨量とはなにかについては、気象庁のウェブページを見てほしい。
大まかにいえば、多数のレーダー観測で得られたエコー強度を、
雨量計による降水強度 (短時間の降水量) の観測値を参照して降水強度になおし、
緯度経度それぞれ等間隔の格子データにしたものである。
(1999年以後のデータでは) 格子間隔は「1 km」とされるが、
正確には、緯度 (1/120) 度、経度 (1/80) 度である。
数量は、1時間降水量であり、1時間の時間区間のおわりの時刻の値としてしめされている。
(最近の年は30分間隔の値があるが、30分ずつ重複した1時間降水量である。)
その開発にかかわった人の著書のうちの解析雨量のところ (数ページ) も参考になる。
- 牧原 康隆, 2020: 『気象防災の知識と実践』。朝倉書店。
[読書メモ]
気象庁のウェブサイトからオンライン公開されているのは「防災情報」のうちの
「今後の雨」にふくまれた最近12時間の実況だけのようである。
気象業務支援センターからオフライン (DVDディスク) で有料配布されており、
立正大学で購入したものを、研究室の Linux マシンのディスク上に置いてある。
気象庁の解析雨量のデータ形式は、2005年のぶんまでは独自の形式、
2006年からはWMO (世界気象機関) できめられた GRIB 第2版 (GRIB2 ともいわれる) である。
解析雨量の値は、1時間降水量 (mm/h) であるが、
有効数字があまり多くないので、1バイト (8ビット) で表現できる 0〜255のレベルでおきかえられ、
さらに、となりどうしの画素の値が同じであるばあいに連長 (run length) 法によって圧縮された
形で収録されている。
(わたしは GRIB2 のデータ形式の概略を理解したものの、解読プログラムを書けるほどよく理解していないが、)
次のウェブ記事で、解析雨量にかぎってGRIB2を解読して地図上に作図する Python プログラムが紹介されており、
そのGRIB2解読部分をそのままつかって、プログラム例を作成した。
- [外] (Qiita) 気象庁解析雨量(GRIB2形式)をPythonで扱う (@vpcf, 2020年07月19日) https://qiita.com/vpcf/items/b680f504cfe8b6a64222
ただし、作図については、@vcpf さんのプログラム例では Basemap というパッケージをつかっていたが、
ここでは Cartopy をつかうように変更した。
研究室 Linux マシンでのデータの配置
hamana と niigata には、
2025-12-06 現在、1995年から2025年6月末までの日付のデータがある。
/raid/kaiseki_uryo/ の下に、
つぎのようなディレクトリ編成でデータを配置している。
(ただし、niigata での /raid/kaiseki_uryo は、/diskb/kaiseki_uryo に symbolic link でつけた別名である。)
- /raid/kaiseki_uryo/
- 1995-1998/
- 1999/
- ...
- 2005/
- 2006/
- ...
- 2024/
この 2024/ の下はつぎのようになっている。
2024年のデータの大部分は 2024/DATA/2024/ の下の、月、日の階層のサブディレクトリの中にある。
2024/DATA/2023/ の下には、日本標準時 (JST) では 2024年だが 世界時 (UTC) ではまだ2023年である
9時間のぶんのデータがある。
- 2024/
- DATA/
- 2023/
- 12/
- 31/
- Z__C_RJTD_20231231150000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin
- ...
- Z__C_RJTD_20231231233000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin
- 2024/
- 01/
- 01/
- Z__C_RJTD_20240101000000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin
- ...
- Z__C_RJTD_20240101233000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin
- ...
- 31/
- ...
- 12/
- DOC/ ... 説明文書
- RAV_GRIB/ ... データを読むプログラム (MS Windows でうごくらしいが、未確認)
- SAMPLE_C/ ... データを読むソースプログラム例 (C言語で書かれている)
2006年以後の他の年のディレクトリ構成も、原則としてこれと同じである。
ただし、2013年から2019年までは、たとえば、上記のとおりならば「2013/DATA/2013/」に
くるべきものを「2013/」のところに置いてしまっている。
これは、立正大学でのデータ保管の際にたまたまそうなってしまったのであって、
気象庁あるいは気象業務支援センターに由来するものではない。
増田のプログラム例はこのディレクトリ構成に対応するように書いているが、
気象業務支援センターのディスクにもとづいてあらたにデータを構築するならば、
2013-2019年を特別あつかいしないように書きかえるべきである。
降水量の分布を地図上に作図するプログラム
2025-05-04 現在、niigata と hamana の /home/masudako/ex/kaiseki_uryo にある。
動作確認は niigata でおこなった。
プログラムのファイル名には「_」 (下線、アンダースコア) をふくんでおり、空白文字はふくんでいない。
2006年以後の解析雨量データの日時は世界時 (UTC) であるが、
例文プログラムでの日時指定や出力画像ファイル名などの日時は、
気象庁が日本語ウェブサイトで公開している地上観測の降水量にあわせて、
日本標準時 (JST) をつかうことにした。
時間降水量の図の表題文字列には、JST と UTC の両方をしめす。
集計値のばあいは JST だけにしてしまったばあいもある。
- 1時間ごとの降水量 (時間降水量) をぬりわけで表示する。
- 3時間ごとに集計した降水量 (3時間降水量) をぬりわけで表示する。
- 1日ごとに集計した降水量 (日降水量)をぬりわけで表示する。
- 複数日の期間で集計した降水量をぬりわけで表示する。単位時間あたりではない。
- 指定した時間区間について集計した降水量をぬりわけで表示する。単位時間あたりではない。
- 複数年について同じ日付の複数日の期間で集計した降水量を、年数でわってぬりわけで表示する。
日数ではわらない。
つぎの注は、ここにあげた例文に共通である。
- プログラム (スクリプト) ファイル中では、日本語文字はコメントだけにふくまれ、その文字コードは UTF-8 になっている。
ブラウザでは文字化けするかもしれないが、
ブラウザに「テキストエンコーディングを修復」の機能があればそれをつかって読める。
なければ、(右クリックから) ダウンロードしてパソコンで読むか、
niigata や hamana 上の emacs で読んでほしい。
- 緯度経度範囲は、lonmin, lonmax, latmin, latmax にあたえる値で指定する。
いくつかの緯度経度範囲に略称をつけ、画像ファイル名にもそれをふくめることにした。
仮の代表例として、北緯 34.5～38.5度、東経136.5～141.5度の領域
(関東地方のほか、中部地方の大部分と東北地方の南部をふくむ) を、
(うまい表現がなくてくるしまぎれだが) 「gkanto」 ("Greater Kanto", 「関東広域」) とした。
- 図は、savefig が True ならば、imagefilepath で指定したところ (例のまま変更しなければ、"../image/") に png 画像ファイルで保存される。
ファイル名には 日本標準時 (JST) で日時がはいるようにしてある。
savefig が False ならば X の画面に表示する。
- 海岸線のほか、都道府県の境をいれている。
その方法は [日本の都道府県境をかく / 中国の「省」の境をかく] に書いた。niigata・hamana 以外の計算機で実行するときは、あらかじめ海岸線データを用意する必要がある。
なお、日本の国土数値情報をつかって、市町村までの行政界を入れることもできる。
そのための変更箇所はコメント文にしてある。
- 数値レベルと図の色わけとの関係は、draw_a_contour_map 関数の中の
levels と colors というリスト変数のなかみで指定している。
colors にあたえる文字列は Matplotlib の List of named colors にあるものである。
(色をその範囲にかぎらずに16進数で指定したほうがよいかもしれない。)
levels の値はあまり深く考えたものではなく、ためしに表示した値がおさまるようにしただけなので、
適宜変更してつかってほしい。
- 数値レベルと図の色わけとの関係について、levelsの下限より下、上限より上の値について
色を指定するしくみ「cmap.set_under」「cmap.set_over」が、
現在 niigata や hamana にはいっているバージョンの matplotlib では、なぜか動作しないので、
levels の下限と上限にデータ中にありそうもない極端な値を入れることによって問題を回避している。
例文プログラムはたびたび改訂することになるだろう。
とくに、複数のプログラムで同じことをやっているところは、
ひとつのファイルに書いて各プログラムからそれを共有するようにしたいのだが、
それを実現する方針はまだたっていない。
降水量の時系列の数値をとりだすプログラム
- 時別値を、緯線経線でかこまれた長方形領域で平均した値を、CSVファイルに書き出す。
(図形出力はしない。下にしめす ku_graph_timeser のプログラムをすこし修正すれば
このCSVを読んでグラフにすることができる。)
- [ku_areamean_timeser_hourly_v20251107.py]
- 日時のうち年は固定。開始と終了の月・日・時を日本時間で指定する。
「時」は「24時」も指定可能にしてあるが、出力の日時表示は翌日の0時となる。
- 領域は、南端・北端の緯度、西端・東端の経度を指定する。
「度、分、秒」で指定可能だが、「度」の小数部つきで指定して「分、秒」を0にしてもよい。
- 出力ファイル名には、開始・終了の日付と、領域の端の緯度・経度がはいるようにした。
- 出力ファイルは CSV。便宜上、列を表計算ソフトウェアの方式で示す。
- A列が日付 ("2024-12-31 23:00" のような形、内容は日本時間)。
- B列が降水量時別値 (単位はmm/h)、ただしデータ欠損は -1 とする。
- C列は長方形領域にふくまれた格子点のうち有効な値があったものの割合 (全部有効ならば 1.0)。
- 1つの格子点の時別値を、CSVファイルに書き出す。
- [ku_timeser_hourly_v20251124.py]
- 日時のうち年・月は固定。開始と終了の日・時を日本時間で指定する。
- 緯度・経度は、
「度、分、秒」で指定可能だが、「度」の小数部つきで指定して「分、秒」を0にしてもよい。
指定した緯度・経度にいちばん近い格子点の値を出力する。
出力ファイル名には、(指定した値ではなく) 格子点の緯度・経度を小数点下3けたで近似した値がはいる。
- 出力ファイルはCSVで、A列が日付、B列が降水量時別値である。データ欠損は「-1.0」と表示される。
1つの時系列をグラフに作図する。
しきい値をこえた日時と最大値を書き出す
緯線経線の長方形領域の降水量の格子点値を書き出す
- [ku_kiridasi_v20251205.py]
- 1時間ごとの解析雨量のうち、指定した緯線・経線でかこまれた長方形の内側の格子点値を、
単純バイナリファイル (64ビット浮動小数点数、増田の符丁「r8l」[最後の字は小文字のエル]) で書き出す。
- 開始日は出力ファイル名にふくまれるが、その他の日時情報は書き出されない。
時刻の範囲は日本時間1時から始めて24時までに固定している。
- 対象となった格子点の緯度と経度はそれぞれ lat_.r8l, lon_.r8l という名の
単純バイナリファイル名に書き出される。
(プログラムの region に地域名をあらわす文字列を指定したばあいは、
その文字列が出力ファイル名の「_」(アンダースコア) と「.」(ピリオド) のあいだにはいる。)
- [ku_map_kiridasareta_v20251207.py]
- 上の kiridasi のプログラムで書きだされたファイルを読み、ku_map_hourly のプログラムと同様に作図する。
土壌雨量指数を計算する [2025-12-07, 改訂 2025-12-25]
「土壌雨量指数」は、気象庁が大雨警報のために計算しているものである。
計算の方式は水文の流出解析のためのタンクモデルと同様で、
土壌雨量指数はタンクにたまっている水の量にあたる。
しかし、目的が水文解析ではなく全国くまなく警報を出すことなので、
モデルパラメータは全国一律に決められている。
(警報のしきい値は地域によって変えられている。)
説明は気象庁ウェブサイトのつぎのページにある。
- [calc_swigrid_ku_v20251223.py]
- 指定した緯線・経線でかこまれた長方形の内側の格子点のそれぞれについて
土壌雨量指数を計算し、1時間きざみで、単純バイナリファイルに書き出す。
緯度経度や日時の指定のしかたは、上の ku_kiridasi と同様である。
- 気象庁の土壌雨量指数の入力データは解析雨量の速報値あるいは予報値の10分値だが、
ここでは、解析雨量(確定版)の1時間値を、10分ごとの6ステップくりかえしあたえることにした。
- 現在のプログラム例では、タンクの水の量の初期値は 0 になっている。
したがって、まとまった量の雨がふるまえの土壌雨量指数は過小評価になっている。
- [map_swigrid_hourly_v20251223.py]
- 上のプログラムで書きだされたファイルの土壌雨量指数を読み、地図上に作図する。
- (数量と色との対応は、解析雨量の1時間値用のままで、まだ土壌雨量指数用に検討されていない。)
- シミュレーションの初期の期間の数値を作図からは省略したいことがあるので、
calc_swigrid の計算の最初の日時を「from」のついた変数で、
作図する最初と最後の日時を「start, end」のついた変数で指定するようにした。
2021-2023年度の教材ページ。機能のほとんどはこのページで紹介する例文にひきついだが、
事例の図や補足的説明はひきついでいないことがあるので、旧ページも残しておく。
2022-2023年度の教材ページ。わたしの時間がとれれば2025年度のスタイルに修正したい。
解析雨量の2005年以前のデータ形式は日本の気象庁の独自のものであり、日時は JST で収録されている。
それを処理するプログラムは2006年度以後のデータ用とは別につくる。
(2021-04-27 .. 2023-09-12) 2025-05-04 .. 2025-05-13, 2025-12-07, 2025-12-25
増田 耕一 (MASUDA Kooiti)