気象庁の全国合成レーダーデータを、 気象業務支援センターからDVDディスクでとりよせました。
DVDディスクにはつぎのものがふくまれています。 格子は、緯度・経度それぞれの等間隔ですが、間隔のおよその長さでよばれています。
【2024年3月以後のものは、「エコー強度、1 km メッシュ、5分間隔」だけになっています。】
おなじデータが、京都大学 生存圏研究所 のデータベースのウェブサイトから、 観測の翌日に公開され、2003年6月以後つづけて蓄積されています。 ここからダウンロードしたものを増田のプログラム例で作図するには、 プログラム中のディレクトリ名のところを変更する必要があります。
5分間隔のエコー強度データについては、 ディレクトリ名・ファイル名だけの変更 (ファイル名の Prr10lv の部分を Prr05lv にかえる、など) で作図できることを確認していますが、 このページで紹介するプログラム例には、まだ反映していません。
10分間隔のエコー強度を読んでみました。 解析雨量のデータと、ドキュメントファイルを見くらべてみると (解析雨量のドキュメントがテキストファイル、全国レーダーはPDFであり、 記録形式に関する整数値が、解析雨量では16進数、全国レーダーでは10進数で書かれているなどの ちがいがあって、まだきちんと照合できていないのですが)、 フォーマットは GRIB2 準拠であるだけでなく、 格子点数、連長 (run length) 圧縮も同じであり、 さらに、エコー強度の数量も、1時間降水量 (10分降水量ではない) に換算されていて、 読むファイル名のきめかただけなおせば、 解析雨量を読むプログラムがそのまま動作しました。
ただし、降水強度 (mm/時) に換算するのに、解析雨量のばあいは 10 でわればよかったのですが、 全国合成レーダーのばあいは、100 でわります。 【わたしが 2021年7月に作成したプログラムは、解析雨量のばあいと同じ係数にしていたので、 表示されていた値は、正しい値の 10 倍の値を表示するものになっていました。】
[解析雨量を読み、作図する] のページの例文を 全国合成レーダーのエコー強度を読むように修正し、 10分ごとのエコー強度を、降水量階級別の色わけで表示するプログラムにして、 動かしてみました。
【データの置きかたは、たとえば2020年6月ぶんのDVDディスクの内容を datapath で指定したディレクトリ (hamana では "/raid/zenkoku_radar/") の下の「202006/」に置いた形を想定しています。】
データの日時は、世界時 (UTC) であたえられていますが、 わたしのプログラムは 日本標準時 (JST) の日時でループをまわすことにしました。
【気象業務支援センターから提供されているディスクは、月ごとに1枚になっていますが、 月の切れめは 世界時 (UTC)の日時できめられており、 (解析雨量のディスクとちがって) 日本時間 (JST) は配慮されていません。 5月のデータをもっていないとき、 6月 1日 0時 (JST) のデータを作図しようとすると、エラーになります。 そのばあい、6月1日については、日本時間 9時 (UTC 0時) からはじめるべきです。 この必要性は、どの年月のディスクを持っているかによるので、 変更点はコメントにしてあります。】
プログラムでは、全領域 (日本全国) のデータを配列に読みこみますが、 draw_a_contour_map() には、表示する領域に対応する部分配列をわたすことにしました。 これで計算機の負荷がだいぶちいさくなったと思います。
雨雲を追跡するために、エコー強度の局所的極大の位置を知りたくなります。
NumPy の関数 argmax() で、配列のうちで最大値がどこにあるかを知ることができます。 その情報を 2次元配列での x方向、y方向の位置 としてつかうには unravel_index() を併用します。 つぎのウェブページを参考にしました。
https://note.nkmk.me/python-numpy-argmax-argmin/ note.nkmk.me NumPy配列ndarrayの最大値・最小値のインデックス（位置）を取得 Posted: 2019-04-15 / Modified: 2019-10-29
試行錯誤しましたが、ここに置くプログラムはつぎのような演算によるものです。 (2*jdy+1)×(2*idx+1)個の格子点 (例では idx = jdy = 20 としたので、41×41個の格子点、つまり約 20 km 四方) からなる template の中心の位置を順にずらしながら、 そのうちの最大値をもとめ、最大値が template の中心の格子点にあり、 その値が、あらかじめきめた しきい値 thresh (例では 20 mm/hr とした) 以上のばあいを 抽出しました。
そして、地図上で、極大の位置に○印、その付近にその値を表示するようにしました。 (色は試行錯誤しましたが、 背景の色がなんであっても字が見えるようにするのはむずかしいです。 ここでは magenta にしてみました。)
たとえば、201 [mm/hr] という値が表示されているならば、 10分間に 33.5 (= 201 / 6) mm の雨がふったことになります。 ただし、エコー強度から降水強度への換算は、多数のシーンを重ねてなりたつのであって、 1時刻について換算した値には大きな不確かさがあるとみるべきでしょう。
画像ファイル名は、"rei" につづいて日時の数字列、それから "jst_" region "_local_max.png" のようにつきます。 ただし region のところにはプログラムの region という変数の値としてあたえた 文字列がはいります。これは地図の緯度経度範囲 (いまのところ極大値探索の緯度経度範囲も同じ) を指定する符丁です。
局所的極大を、図に表示するだけでなく、 CSVファイルに書きだすこともします。
まだ、自動追跡はできていません。
手動追跡といえるものの手順を説明します。 過去のプログラム例が [旧版教材ページ]に ありますが、今回のプログラムの出力ファイル形式にはまだあわせていません。