解析雨量を Pythonで読み、作図する (試作版)


作図プログラム

気象庁の解析雨量 (2006年以降) のデータを、Python で読み、 Python のパッケージ「Cartopy」をつかって地図上に作図するプログラムを作成した。

(Python での作図をはじめたばかりのわたしが書いてとりあえず動いたという段階のもので、 変なところがあるかもしれない。)

結果の画像の例 (2019年10月25日 日本時間 12:00 の1時間雨量)

rr201910251200jst.png

プログラム例文

[2021-06-02] 1日ごとの値を作図するプログラム例を追加しました。 また、1時間ごとの値を作図するプログラム例を改訂しました (実質の処理はかわらず、プログラム例としての見とおしをよくするための改訂です)。

[2021-06-02] 日あるいは月の切れめが世界時と日本標準時でずれるところの処理で失敗していたところがあり、 修正しましたが、まだ問題があるかもしれません。年の切れめの処理はまだ未確認です。

[2021-06-05] 2020年のデータを処理できるようにプログラムを改訂しました。

気象庁の解析雨量のデータ形式は、2005年のぶんまでは独自の形式、 2006年からはWMO (世界気象機関) できめられた GRIB version 2 (GRIB2 ともいう) である。

解析雨量の値は、本来、1時間降水量 (mm/h) である。 (30分間隔のばあいも、その時刻までの1時間の降水量である。) しかし、有効数字があまり多くないので、1バイト (8ビット) で表現できる 0〜255のレベルでおきかえられ、 さらに、となりどうしの画素の値が同じであるばあいに連長 (run length) 法によって圧縮されている。

ここでのプログラム例文は、GRIB2 形式が採用された 2006 年以後のデータを対象としている。 2005年以前の独自形式には対応していない。 また、2006年1月1日のデータをまだうまくあつかえない。 (世界時で前年に属する観測値をふくむためである。 GRIB2形式のデータは用意されているのだが、わたしのプログラムがまだ対応していない。) 他の年の1月1日についても、うまくいかないことがあるかもしれない。

GRIB2形式のデータでは、日時は世界時 (UTC) であたえられているが、 利用者は日時を日本標準時 (JST) であたえられるようにし、 図には UTC と JST の両方の日時を表示するようにした。 画像ファイル出力のファイル名には JST をつかい、それを明示するようにした。 (Python の datetime パッケージの機能をつかった。)

入力データとして、年をあらわす4けたの数字のディレクトリの下に、 その年のデータをふくむ気象業務支援センターのDVDディスクの内容がそのディレクトリ構成で置かれていることを前提とする。 (ただし、2013年から2019年までのデータは、年をあらわす4けたの数字のディレクトリの下に、 直接データファイルが置かれていることを前提とする。 ディレクトリ構成のちがいには、プログラムで対処する。 [この部分、2021-06-05 改訂])

つかいかた

プログラムやデータを置くディレクトリ名の文字としては、ASCIIの1バイト文字コードで 表現できるアルファベット・数字・算数記号にかぎったほうがよいだろう。 日本語の漢字・かななどがあるとうまくいかないかもしれない。

この例文を各自の目的にやくだてるには、 プログラムを編集して、文字を書きかえる必要がある。 (日時の情報は、人がキーボードからうちこむようにしたり、 別のファイルから読みこむようにすることもできるが、 そのための書きかえかたの説明は、まだ用意していない。)

例文プログラムの「##### Main program #####」というコメントのあとにある 「datapath」という変数に代入する値を、自分がデータを置いたディレクトリ名に変える。 (ディレクトリ区切り文字は、Windows ではふつうバックスラッシュか円マーク ('\') だが、 ここでは Linux と共通にするためにスラッシュ ('/') にしている。このほうがよいと思う。)

対象日時を指定する。

図形の表示さきを指定する。

対象地域を指定する。

数値レベルと図の色わけとの関係は、draw_a_contour_map 関数の中の levels と colors という リスト変数のなかみで指定している。現在あたえている色のくぎりの数値レベルはつぎのとおり。 (あまり深く考えたものではなく、ためしに表示した値がおさまるようにした。)

なお、解析雨量の値が欠損になっているばあいは、関数 load_jmara_grib2 の中であたえる値は (@vcpf さんのプログラムをひきついだので) -10 である。 メインプログラムで、降水量の値を 0.1 mm 単位から mm 単位になおすために 10 でわっているので、 欠損をあらわす値は -1 になる。表示の際はこれを考慮に入れる。 (いまの配色では、値が -1 (欠損) ならば薄い赤、0 (無降水) ならば白、0.4 (微量の降水) ならば薄い黄色としてある。)


時系列書き出しプログラム [2021-06-03 追加, 2021-06-05 改訂]

解析雨量のGRIB2形式データを読みこむところは共通で、作図するのではなく、 地点をえらんで、その位置の降水量の値の時系列をテキストファイルに書き出すプログラムを試作している。 こちらは、Cartopy も matplotlib もつかわない。

まず、指定された緯度・経度にいちばん近い ます目 の降水量の1時間ごとの値を書き出すプログラムをつくった。

月単位の期間を指定して、その全部の日の全部の時刻の値を書き出すことにした。 ただし、年にまたがる期間は指定できない (それがほしいばあいは複数回にわけて実行する)。

つかいかた

主プログラムのはじめ (スクリプトの中では途中になるが、「##### Main Program #####」というコメントのはいった直後) の、 つぎの変数の値をあたえている代入文を、ほしい位置の緯度・経度に書きかえる。

すこしさきの「# 対象日付」というコメントのあとのつぎの変数の値をあたえている代入文を、 ほしい日付に書きかえる。

「datapath」という変数に代入する値を、自分がデータを置いたディレクトリ名に変える。

スクリプトを実行する。(Anaconda Powershell Propmpt ならば、 コマンドラインで「python ku_timeser_hourly_v20210602.py」)

出力データの形

つぎのようなテキストファイルができる。

日降水量の時系列書き出しプログラム

指定された緯度・経度にいちばん近い ます目 の降水量の1日ごとの値を書き出すプログラムをつくった。

日本時間の0時から翌日0時までの集計値である。 欠損値は単純に集計に入れない。 (もし観測されていたとしたら、降水量はもっと大きかった可能性がある。) 集計につかわれた有効な観測値の個数を書き出しておく。

作業手順は1時間ごとの値を書き出すプログラムと同様である。

結果のファイル名も同様で、hourly のかわりに daily である。

結果のファイルの列の意味はつぎのとおり。


時間降水量の領域平均値の時系列書き出しプログラム [2021-06-18 追加]

緯度・経度座標の長方形領域で平均した降水量の1時間ごとの値を書き出すプログラムをつくった。 ひとまず バージョン1 を出しておく。今後、改訂する。

つぎのような不満がのこるプログラムである。

つかいかた

次のところを書きなおす。

スクリプトを実行する。(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度の間の値である。

結果のファイルの列の意味はつぎのとおり。

今後の作業についての構想

複数の地点の値がほしいとき、実行能率の面では、毎時の全国のデータを読むごとに 複数の地点の値を書き出すように、プログラムを改造したほうがよい。 ただし、そうするとプログラムがいくらか複雑になり、プログラム保守の面では、うまくないかもしれない。

領域平均を書き出すプログラムには上にのべたように不満があるので、改訂したい。


実行環境

ここにあげたプログラム例は、MS Windows 日本語版で実行することを想定し、Windows 10 上でつくっている。

Windows での実行は、Anaconda をインストールし、その中で Cartopy をいれたうえで、 Anaconda Powershell Prompt から「 python なになに.py 」で起動するか、 Anaconda といっしょにインストールされる Spyder に読みこんで実行するのがよさそうだ。 (わたしは Anaconda Powershell Prompt から実行している。)

[Windows で Pythonをつかう準備 (Anaconda, cartopy のインストール) と 実行]

【Anaconda はディスクを大量に消費するので (標準的インストールをしただけで 3 ギガバイト)、 最小版の Miniconda を入れて、必要に応じてパッケージを追加していったほうがよいかもしれない。 また、(Python言語処理系はフリーだが) Anacondaのライセンスは、個人利用や教育利用は無料だが、 営利企業での利用は有料になることがあるので、とくに学生が卒業後に仕事でつかおうとするときは 注意が必要だ。】

プログラムの編集には、テキストエディタをつかう。 Windows 10の「メモ帳」でもなんとかなるだろうが、 もう少し高度なものをいれたほうがよいかもしれない。 (例文の日本語文字コードは UTF-8 にしてある。古い Windows のメモ帳だと文字化けするかもしれない。) Windowsの「エクスプローラ」で、「.py」でおわるファイル名が Python と関係づけられていると、 エクスプローラでファイルを単純に左クリックした際にそれを Python で実行しようとするが、うまくいかない。 編集したいときは、右クリックで「プログラムから開く」を選択し、そこからテキストエディタを開くようにする。 関係づけを変えることができれば、左クリックでテキストエディタが起動するようにしておくのがよいと思う。


参考にした情報源

解析雨量 (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データから気温の水平分布を描画する

【ウェブ上でキーワード検索をかけたり、それで見つかったものからいもづる式にたどったりして、 ブログ記事などかなり多くの記事をながめたうえで読むものをえらんだのだが、 たまたま qiita.com のページを参照したものがつづいた。 このサイトはたしかに有用な記事があるが、このサイトの記事がみんな信頼できるわけでもない。】


[2021-04-27 までの作業覚え書き] (うちわむけメモ)


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