[2021-07-09] [解析雨量を Pythonで読み、作図する (試作版)]の続編のようなものです。
気象庁の全国合成レーダーデータを、 気象業務支援センターからDVDディスクでとりよせました。
DVDディスクにはつぎのものがふくまれています。 【[2024-08-22 注] 2024年3月以後のものは、「エコー強度、1 km メッシュ、5分間隔」だけになっています。】 格子は、緯度・経度それぞれの等間隔ですが、間隔のおよその長さでよばれています。
[2024-08-22 補足] おなじデータが、 京都大学 生存圏研究所 のデータベースのウェブサイトから、 観測の翌日に公開され、2003年6月以後つづけて蓄積されています。 ここからダウンロードしたものを増田のプログラム例で作図するには、 プログラム中のディレクトリ名のところを変更する必要があります。 (そのための説明は、今後整備しようと思っています。)
[2021-07-09, 2021-07-29] エコー強度を地図上に作図するプログラムをつくりました。
[2021-07-30, 2021-12-08] エコー強度の局所的極大を抽出して、手作業での追跡をへて、極大の移動経路を作図する、 一連のプログラムを試作しました。(追跡の自動化はできていません。)
[2022-04-28] 2021年7月からこれまでここにしめしていたプログラム例の、エコー強度を降水強度に換算した値について、 わたしが換算定数をまちがえており、正しい値はこれまでプログラムや図の例で示した値の 10分の1です。 2022-04-28 にプログラム例や図をさしかえました。
[2022-04-28] エコー頂高度を読み作図するプログラム例を追加しました。
[2022-11-07] 作図プログラム例の都道府県境の情報の読みかたを [日本の都道府県境をかく / 中国の「省」の境をかく] のページでのべたように 変更しました。 【[2024-08-22 注] 都道府県境のデータが置かれた場所は変数 shapefilepath で指定してあり、 実行する計算機によって変更する必要があるかもしれません。】
[2024-08-22] プログラム例をさしかえました。こまかい改訂をしています。
[2024-08-22 注] 5分間隔のエコー強度データについては、 ディレクトリ名・ファイル名だけの変更 (ファイル名の Prr10lv の部分を Prr05lv にかえる、など) で作図できることを確認していますが、 このページで紹介するプログラム例には、まだ反映していません。
[2024-08-22 注] 現在 niigata や hamana にはいっている matplotlib では、 contourf と colorbar で、cmap.set_under, cmap.set_over が無効になっており、 levels の範囲をひろげるようにプログラムを変更して対処しています。 このページで紹介するプログラム例は、cmap.set_{under,over} が有効なばあいのものです。
ひとまず、10分間隔のエコー強度を読んでみました。 解析雨量のデータと、ドキュメントファイルを見くらべてみると (解析雨量のドキュメントがテキストファイル、全国レーダーはPDFであり、 記録形式に関する整数値が、解析雨量では16進数、全国レーダーでは10進数で書かれているなどの ちがいがあって、まだきちんと照合できていないのですが)、 フォーマットは GRIB2 準拠であるだけでなく、 格子点数、連長 (run length) 圧縮も同じであり、 さらに、エコー強度の数量も、1時間降水量 (10分降水量ではない) に換算されていて、 読むファイル名のきめかただけなおせば、 解析雨量を読むプログラムがそのまま動作しました。 ただし、降水強度 (mm/時) に換算するのに、解析雨量のばあいは 10 でわればよかったのですが、 全国合成レーダーのばあいは、100 でわります。 【わたしが 2021年7月に作成したプログラムは、解析雨量のばあいと同じ係数にしていたので、 表示されていた値は、正しい値の 10 倍の値を表示するものになっていました。】 [2022-04-28 改訂]
[解析雨量を読み、作図する] のページの例文 ku_map_hourly_v20210701.py を、全国合成レーダーのエコー強度を読むように修正し、 そこで便宜上「#狭関東」とよぶことにした領域について、 10分ごとのエコー強度を、降水量階級別の色わけで表示するプログラムにして、 動かしてみました。
【[この項目 2024-08-22 改訂] データの置きかたは、たとえば2020年6月ぶんのDVDディスクの内容を datapath で指定したディレクトリ (hamana では "/raid/zenkoku_radar/") の下の「202006/」に置いた形を想定しています。】
データの日時は、世界時 (UTC) であたえられていますが、 わたしのプログラムは 日本標準時 (JST) の日時でループをまわすことにしました。
[2021-07-29] プログラムを改訂しました。改訂した点はつぎのとおりです。
[2022-11-07] プログラムの県境の作図のところを改訂しました。
[2024-08-22] メインプログラムから draw_contour_map にひきわたす配列の上限を修正しました。
[2021-07-29] 雨雲を追跡するために、エコー強度の局所的極大の位置を知りたくなります。
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
まず、図示する領域全体での最大値の位置を知り、そこに〇印と、 エコー強度 (降水強度換算) の値を表示するプログラムをつくってみました。
(最大値を表示する色は試行錯誤していますが、 背景の色がなんであっても字が見えるようにするのはむずかしいです。 ここでは magenta にしてみました。)
この実行例では 201 [mm/hr] という値が表示されており、 10分間に 33.5 (= 201 / 6) mm の雨がふったことになります。 ただし、エコー強度から降水強度への換算は、多数のシーンを重ねてなりたつのであって、 1時刻について換算した値には大きな不確かさがあるとみるべきでしょう。[2022-04-28 部分修正]
[2021-07-30] 局所的極大をもとめてみました。
試行錯誤しましたが、ここに置くプログラムはつぎのような演算によるものです。 21×21個の格子点からなる template の中心の位置を順にずらしながら、 そのうちの最大値をもとめ、最大値が template の中心の格子点にあり、 その値が、あらかじめきめたしきい値 (ここでは 20 mm/hr とした) 以上のばあいを 抽出しました。
[2021-12-08] 局所的極大を、図に表示するだけでなく、 テキストファイルに書きだすこともするようにプログラムを改造しました。 なお、template の大きさを 41×41個の格子点にし、 処理対象の時間間隔を 10 分間隔にする (材料のデータにある時刻をとばさずに処理する) などの変更をしています。)
この時刻には熊谷をふくむ埼玉県北部にエコーが出ていますが、 極大はそのうち西部の長瀞あたりで、値は 150 mm/hr となっています。 (10分降水量 25 mm となりますが、降水量の値としては大きな不たしかさをもつとみるべきでしょう。) [2022-04-28 部分修正]
[2021-12-08] 書きだした極大値のデータを、緯度・経度でしぼりこみます。 趣旨は、[表データの条件にあてはまる行を抽出 -- 緯度経度が指定範囲内の地点を抽出] の ページに書いたものと同じですが、 zgrei_map_local_max.py が出力したファイル rei_local_max.txt のうちから抽出するようにした プログラムがつぎのものです [2022-04-28 文章修正]。
[2021-12-08] それから、エディタをつかった手作業で、 しぼりこまれた極大値のデータから行をぬきだし、 それぞれの雨雲の経路をあらわすテキストファイルを、 追いかける雨雲の数だけ別々につくります。 そして、つくったテキストファイル名のファイル名を1行に1つずつ書きならべた filelist.txt をつくります。
[2021-12-08] そして、つぎのプログラムを実行すると、 雨雲の経路が地図上に折れ線として表示されます。
[2022-04-28] エコー頂高度のデータを読み、地図上に作図するプログラム例をつくりました。
エコー頂高度の空間分解能は「2.5kmメッシュ」とされており、 くわしくは、経度間隔 1.875 分、緯度間隔 1.5分 (度分秒の分) です。 また、エコー頂高度のデータ値は、エコーがないばあいが 0、 2 km 未満が 1、2-4 km が 3, ... 12-14 km が 13, 14 km 以上が 15 となっています。 そのまま km 単位の値とみなして表示することにします。
[2023-08-29] メインプログラムから draw_contour_map にひきわたす配列の上限を修正しました。
