全国合成レーダーデータを Pythonで読み、作図する (試作版)


[2022-04-28] 2021年7月からこれまでここにしめしていたプログラム例の、エコー強度を降水強度に換算した値について、 わたしが換算定数をまちがえており、正しい値はこれまでプログラムや図の例で示した値の 10分の1です。 2022-04-28 にプログラム例や図をさしかえました。

[2022-04-28] エコー頂高度を読み作図するプログラム例を追加しました。

[2022-11-07] 作図プログラム例の都道府県境の情報の読みかたを [日本の都道府県境をかく / 中国の「省」の境をかく] のページでのべたように 変更しました。


[2021-07-09] [解析雨量を Pythonで読み、作図する (試作版)]の続編のようなものです。

気象庁の全国合成レーダーデータを、 気象業務支援センターからDVDディスクでとりよせました。 対象期間は 2019年と2020年のそれぞれ6月から9月です。

DVDディスクにはつぎのものがふくまれています。 格子は、緯度・経度それぞれの等間隔ですが、間隔のおよその長さでよばれています。


エコー強度を読み、作図する

ひとまず、10分間隔のエコー強度を読んでみました。 解析雨量のデータと、ドキュメントファイルを見くらべてみると (解析雨量のドキュメントがテキストファイル、全国レーダーはPDFであり、 記録形式に関する整数値が、解析雨量では16進数、全国レーダーでは10進数で書かれているなどの ちがいがあって、まだきちんと照合できていないのですが)、 フォーマットは GRIB2 準拠であるだけでなく、 格子点数、連長 (run length) 圧縮も同じであり、 さらに、エコー強度の数量も、1時間降水量 (10分降水量ではない) に換算されていて、 読むファイル名のきめかただけなおせば、 解析雨量を読むプログラムがそのまま動作しました。 ただし、降水強度 (mm/時) に換算するのに、解析雨量のばあいは 10 でわればよかったのですが、 全国合成レーダーのばあいは、100 でわります。 【わたしが 2021年7月に作成したプログラムは、解析雨量のばあいと同じ係数にしていたので、 表示されていた値は、正しい値の 10 倍の値を表示するものになっていました。】 [2022-04-28 改訂]

「解析雨量を読み、作図する」ページの例文 ku_map_hourly_v20210701.py を、全国合成レーダーのエコー強度を読むように修正し、 そこで便宜上「#狭関東」とよぶことにした領域について、 10分ごとのエコー強度を、降水量階級別の色わけで表示するプログラムにして、 動かしてみました。

(データの置きかたは、たとえば2020年6月ぶんのDVDディスクのうち 「composition/Reflection_Intensity/」の下の全体を、 作業用ディレクトリの横ならびの「202006」ディレクトリの下に置いた形を想定しています。)

[2021-07-29] プログラムを改訂しました。改訂した点はつぎのとおりです。

[2022-11-07] プログラムの県境の作図のところを改訂しました。

rei202006121500jst.png


極大値の位置を知る

[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 にしてみました。)

rei202006121500jst_regional_max.png

この実行例では 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 分間隔にする (材料のデータにある時刻をとばさずに処理する) などの変更をしています。)

rei202006121500jst_local_max.png

この時刻には熊谷をふくむ埼玉県北部にエコーが出ていますが、 極大はそのうち西部の長瀞あたりで、値は 150 mm/hr となっています。 (10分降水量 25 mm となりますが、降水量の値としては大きな不たしかさをもつとみるべきでしょう。) [2022-04-28 部分修正]

[2021-12-08] 書きだした極大値のデータを、緯度・経度でしぼりこみます。 趣旨は、[表データの条件にあてはまる行を抽出 -- 緯度経度が指定範囲内の地点を抽出] の ページに書いたものと同じですが、 zgrei_map_local_max_v20220428.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 単位の値とみなして表示することにします。

reth202006121500jst.png


2021-07-09, 2021-07-13, 2021-07-16, 2021-07-30, 2021-12-08, 2022-04-28, 2022-06-03, 2022-11-07
増田 耕一 (MASUDA Kooiti)
[計算機を利用した気象分野の実習の教材] [教材目録ページ] [日本語入り口ページ] [ホームページ]
[セミナーのページ]