気象庁の解析雨量を 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-05-04 現在、1995年から2024年末までの日付のデータがある。 /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時間ごとの降水量 (時間降水量) をぬりわけで表示する。
- [ku_map_hourly_v20250509.py] [2025-05-09 バグ修正]
- 日時の指定は、年、月はひとつだけにかぎる。その中で「何日何時」から「何日何時」までかを指定する。
3時間ごとに集計した降水量 (3時間降水量) をぬりわけで表示する。
- [ku_map_3hourly_v20250503.py]
- 日時の指定は、年、月はひとつだけにかぎる。その中で何日から何日までかを指定する。指定された各日の 3, 6, 9, 12, 15, 18, 21, 24 時までの3時間集計値を出力する。
1日ごとに集計した降水量 (日降水量)をぬりわけで表示する。
- [ku_map_daily_v20250503.py]
- 日時の指定は、年、月はひとつだけにかぎる。その中で何日から何日までかを指定する。指定された各日の24時までの24時間集計値をそれぞれ出力する。
複数日の期間で集計した降水量をぬりわけで表示する。単位時間あたりではない。
- [ku_map_multiday_v20250509.py]
- 日時の指定は、何年何月何日から何年何月何日までかを指定する [2025-05-09 改訂]。指定された (24時間の倍数の) 期間での集計値を出力する。
指定した時間区間について集計した降水量をぬりわけで表示する。単位時間あたりではない。
- [ku_map_interval_v20250504.py]
- 日時の指定は、「何年何月何日何時から、何年何月何日何時まで」という形で指定する。集計値には、開始時刻までの1時間の値はふくまれない。
複数年について同じ日付の複数日の期間で集計した降水量を、年数でわってぬりわけで表示する。日数ではわらない。
- [ku_map_multiyear_v20250509.py]
- 日時の指定は、何年から何年まで、何月何日から何月何日までかを指定する。

つぎの注は、ここにあげた例文に共通である。

プログラム (スクリプト) ファイル中では、日本語文字はコメントだけにふくまれ、その文字コードは 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ファイルに書き出す。 (図形出力はしない。このCSVを読んでグラフにするプログラムを追って整備する予定。)
- [ku_areamean_hourly_csv_v20250512.py]
- 日時のうち年・月は固定。開始と終了の日を指定する。日本時間でその日の1時から24時までを書き出す。
- 領域は、南端・北端の緯度、西端・東端の経度を指定する。「度、分、秒」で指定可能だが、「度」の小数部つきで指定して「分、秒」を0にしてもよい。
- 出力ファイル名には、開始・終了の日付と、領域の端の緯度・経度がはいるようにした。
- 出力ファイルは CSV。便宜上、列を表計算ソフトウェアの方式で示す。
  - A列が日付 ("2024-12-31 23:00" のような形、内容は日本時間)。
  - B列が降水量時別値 (単位はmm/h)、ただしデータ欠損は -1 とする。
  - C列は長方形領域にふくまれた格子点のうち有効な値があったものの割合 (全部有効ならば 1.0)。

今後の予定

[旧年度のページ]にあってまだ改訂していないものを、ここに追加する予定。
2005年以前の解析雨量については、旧年度のページ [1995-2005年の解析雨量を読む] があるが、要望があれば今年度のスタイルに改訂したい。 (2005年以前のデータ形式は日本の気象庁の独自のものであり、日時は JST で収録されている。)

(2021-04-27 .. 2023-09-12) 2025-05-04 .. 2025-05-13
増田耕一 (MASUDA Kooiti)

[計算機を利用した気象分野の実習の教材] [立正大学用の教材] [教材目録ページ] [日本語入り口ページ] [ホームページ]へ