線状降水帯予測プロトタイプを作ってみる【r001】Step 3 ー 3時間積算雨量から強雨域候補を抽出する
はじめに
Step 2では、Step 1でそろえた1時間降水量から3時間移動積算雨量を作りました。
実データでは、DIASのMSMについて、
input (39, 505, 481)
output (37, 505, 481)
max 3h 219.719 mm
verify PASS
まで確認できています。
Open-Meteo側も、
input (12, 9, 9)
output (10, 9, 9)
max 3h 11.500 mm
verify PASS
となり、取得元が違っても同じStep 2処理を通せるところまで確認できました。
ここからStep 3へ進みます。
今回の目的は、3時間積算雨量Gridの中から、まとまった強雨域を見つけてPolygonとして保存することです。
処理の流れは次のようにします。
Step 3で初めて「雨域の形」を扱います。
ただし、ここで抽出するものをそのまま「線状降水帯」とは呼びません。
気象庁の線状降水帯に関する情報では、降水量、面積、形状だけでなく、キキクルなどの災害危険度も組み合わせて判定しています。
このプロトタイプでは、そのうち雨量と形状に関する数値条件を参考にして、強雨域候補として扱います。
Step 3のソースファイル(https://godo-tys.jp/downloads/linear_rainband_r001_step3.zip)
参考:
- 線状降水帯予測プロトタイプ r001 全体構成
https://qiita.com/rino_yume/items/171a7009633f8a11681c - Step 1 完了編
https://qiita.com/rino_yume/items/055986af6f1569b0f90c - 気象庁 線状降水帯に関する情報
https://www.jma.go.jp/jma/kishou/know/bosai/kishojoho_senjoukousuitai.html
1. Step 3で扱う条件
2026年時点の気象庁の線状降水帯直前予測では、降水域について次のような数値条件が使われています。
前3時間降水量100 mm以上の領域
面積500 km²以上
長軸 / 短軸比2.5以上
領域内最大3時間降水量145 mm以上
このほか、キキクルによる危険度条件があります。
一方、発生情報では最大3時間降水量の条件は150 mm以上です。
今回作るのは予測側のプロトタイプなので、設定の初期値は直前予測の数値条件を参考に、次のようにしました。
candidate:
input_variable: precipitation_3h
threshold_mm: 100.0
connectivity: 8
min_area_km2: 500.0
min_axis_ratio: 2.5
min_peak_mm: 145.0
ただし、これを気象庁の判定そのものとしては使いません。
本実装では、まず100 mm/3h以上のすべての連結領域を保存します。
そのうえで、
area_km2 >= 500
axis_ratio >= 2.5
max_rainfall_mm >= 145
を満たした場合だけ、
is_candidate = true
とします。
これなら、「100 mm以上の雨域は存在したが、面積が小さかった」「線状性が足りなかった」といった途中の情報も残せます。
2. Step 2の成果物だけを入力にする
Step 3では、Open-MeteoやDIASからデータを取り直しません。
入力はStep 2のmanifestです。
DIASなら、
data/manifests/step2_accumulation/3h/
└─ dias_msm_archive/
└─ run_20260612T2100Z.json
Open-Meteoなら、
data/manifests/step2_accumulation/3h/
└─ openmeteo_jma_msm/
└─ latest_20260920T005852Z.json
です。
Step 3はStep 2 manifestを最初にverifyし、壊れたZarrや変更された入力をそのまま使わないようにしています。
Step 2 manifest
↓
SHA-256確認
↓
Step 2 Zarr再検証
↓
候補抽出
この考え方はStep 1、Step 2から変えていません。
3. まず100 mm以上のセルを取り出す
入力は、
precipitation_3h
です。
各時刻について、
mask = rainfall >= threshold_mm
とします。
既定値では、
threshold_mm = 100
です。
ここでは単にセルを抜き出すだけで、まだ雨域とは考えません。
4. 隣接したセルを連結成分としてまとめる
次に、隣り合った強雨セルを1つのまとまりにします。
今回はscipy.ndimage.label()を使います。
labels, component_count = ndimage.label(
mask,
structure=structure,
)
既定では8近傍です。
┌─┬─┬─┐
│●│●│●│
├─┼─┼─┤
│●│■│●│
├─┼─┼─┤
│●│●│●│
└─┴─┴─┘
中央のセルに対して、上下左右だけでなく斜めも接続とみなします。
降水域は格子に対して斜め方向へ伸びることもあるので、r001では8近傍を既定にしました。
4近傍へ変更したい場合は設定で切り替えられます。
5. 連結成分をPolygonへ変換する
ラベル付けした格子は、そのままではGISで扱いにくいため、Polygonへ変換します。
ここではrasterio.features.shapes()を使います。
格子ラベル
↓
rasterio.features.shapes
↓
Polygon / MultiPolygon
↓
EPSG:4326
Step 1から緯度軸は南→北で統一しています。
一方、rasterioは画像の上側を北として扱うため、Polygon化の直前だけ配列を上下反転しています。
また、格子点をセル中心として扱い、半セル外側をPolygonの境界にしています。
6. 面積をdegree²で計算しない
緯度経度のPolygonに対して、
polygon.area
をそのまま使うと、単位はdegree²になります。
これでは500 km²という条件と比較できません。
Step 3ではpyproj.Geodを使い、WGS84楕円体上の測地面積を求めます。
geod = Geod(ellps="WGS84")
area_m2, _ = geod.geometry_area_perimeter(geometry)
area_km2 = abs(area_m2) / 1_000_000
これで、
area_km2
として保存できます。
7. 長軸・短軸は局所座標で求める
線状性を見るには、
長軸 / 短軸
が必要です。
経度1度の距離は緯度によって変わるため、緯度経度のまま長さを測るのは避けます。
各雨域の重心を中心に、局所的なAEQD座標へ投影します。
WGS84 Polygon
↓
候補重心を中心とするAEQD
↓
Minimum Rotated Rectangle
↓
長軸
短軸
長短軸比
向き
最小回転矩形の4辺から、
major_axis_km
minor_axis_km
axis_ratio
orientation_deg
を計算します。
orientation_degはStep 4の雨域追跡でも使えるように残します。
8. 最大雨量だけでなく途中の特徴量も残す
GeoParquetには、次の属性を保存します。
region_id
run_id
provider_id
valid_time
component_label
cell_count
threshold_mm
max_rainfall_mm
mean_rainfall_mm
area_km2
major_axis_km
minor_axis_km
axis_ratio
orientation_deg
centroid_lon
centroid_lat
meets_min_area
meets_min_axis_ratio
meets_min_peak
is_candidate
geometry
例えば、候補にならなかった場合でも、
meets_min_area = true
meets_min_axis_ratio = false
meets_min_peak = true
のように残ります。
判定結果だけ保存するより、後で閾値を見直しやすくなります。
9. GeoParquetへ保存する
Step 3の成果物はGeoParquetで保存します。
data/processed/candidates/3h/
├─ dias_msm_archive/
│ └─ run_20260612T2100Z.parquet
│
└─ openmeteo_jma_msm/
└─ latest_20260920T005852Z.parquet
GeoParquetにしておけば、後続のStep 4だけでなく、
DuckDB
Python
QGIS
WebGIS
からも扱いやすくなります。
保存しただけで完了にはせず、もう一度読み直して、
- 行数
- CRS
is_candidate
が変わっていないことまで確認します。
10. 確認用PNGも作る
Step 3でも、計算結果を目で確認できるPNGを出します。
data/preview/candidates/3h/
└─ <provider>/
└─ <run_id>_candidates.png
全時刻の中から、最大3時間雨量が最も大きい時刻を代表時刻として選びます。
その時刻の、
3時間積算雨量Grid
強雨域Polygon
候補Polygon
を重ねて表示します。
ここで作る画像はStep 6のWebGISではなく、処理結果を確認するための簡易図です。
11. Step 3もnetworkなしで実行する
Step 3は外部APIへ接続しません。
Step 2のZarrだけを読みます。
そのためComposeは、
step3:
network_mode: none
としています。
Step 2でDockerのbridge address poolが枯渇したことがあったため、外部通信が不要な解析Stepではnetworkを作らない方針にしました。
Step 1 データ取得 networkあり
Step 2 3時間積算 networkなし
Step 3 強雨域抽出 networkなし
と分けています。
12. Step 2のデータを引き継ぐ
Step 2正式完了版のdata/をStep 3へコピーします。
cp -a \
../linear_rainband_r001_step2_complete/data/. \
data/
手元のStep 2ディレクトリがlinear_rainband_r001_step2_fix003の場合は、そちらからコピーしても構いません。
13. Dockerをビルドする
cp .env.example .env
sed -i "s/^LOCAL_UID=.*/LOCAL_UID=$(id -u)/" .env
sed -i "s/^LOCAL_GID=.*/LOCAL_GID=$(id -g)/" .env
docker compose build --no-cache
Step 3では、Step 1とStep 2に加えて、
SciPy
Rasterio
Shapely
PyProj
GeoPandas
PyArrow
を使います。
14. 配布物を検証する
docker compose run --rm step3 \
./scripts/validate_package.sh
検査内容はこれまでと同じです。
dependency consistency
Shell / Python syntax
Ruff format
Ruff check
mypy
pytest
CLI / source contract
package source files
Step 3では、
- 細長い強雨域が候補になる
- コンパクトな強雨域は候補にならない
- 斜め接続を8近傍で同じ雨域にできる
- 100 mm未満なら0件になる
- 不規則な緯度経度格子を拒否する
- NaNを拒否する
- PNGを生成できる
といったテストを追加しています。
15. DIAS MSMで実行する
Step 2で作ったDIASのmanifestを指定します。
docker compose run --rm step3 \
rainband-step3 extract \
--manifest \
/workspace/data/manifests/step2_accumulation/3h/dias_msm_archive/run_20260612T2100Z.json \
--config /workspace/config/step3.yaml
Step 2で、このRunの最大3時間雨量は、
219.719 mm
まで確認できています。
ただし、最大雨量が145 mmを超えているだけでは候補になるとは限りません。
面積500 km²以上
長短軸比2.5以上
も必要です。
したがって、候補数は実際にStep 3を動かして確認します。
16. Step 3の成果物をverifyする
候補抽出後は、
docker compose run --rm step3 \
rainband-step3 verify \
--manifest \
/workspace/data/manifests/step3_candidates/3h/dias_msm_archive/run_20260612T2100Z.json
を実行します。
ここでは、
Step 2 source manifest
GeoParquet SHA-256
PNG SHA-256
GeoParquet再読込
CRS
行数
候補件数
Geometry validity
を確認します。
17. Open-Meteoは「0件」の確認にも使える
Step 2で確認したOpen-Meteo Runは、
最大3時間雨量 = 11.500 mm
でした。
今回の抽出閾値は100 mmなので、このRunでは強雨域は0件になるはずです。
docker compose run --rm step3 \
rainband-step3 extract \
--manifest \
/workspace/data/manifests/step2_accumulation/3h/openmeteo_jma_msm/latest_20260920T005852Z.json \
--config /workspace/config/step3.yaml
0件は失敗ではありません。
region_count = 0
candidate_count = 0
として、空のGeoParquet、PNG、manifestまで保存できることを確認します。
この経路は、「候補がないときも正常に終わるか」という負のテストになります。
18. Step 3の完了条件
今回も、実機で最後まで通ってからStep 3をCOMPLETEにします。
[ ] Docker build PASS
[ ] validate_package.sh 8/8 PASS
DIAS
[ ] Step 2 manifest verify PASS
[ ] 強雨域抽出 PASS
[ ] GeoParquet保存・再読込 PASS
[ ] PNG生成 PASS
[ ] Step 3 manifest生成 PASS
[ ] Step 3 verify PASS
Open-Meteo
[ ] 0件経路 PASS
[ ] 空GeoParquet保存・再読込 PASS
[ ] PNG生成 PASS
[ ] manifest verify PASS
これらが通ったら、r001 Step 3を完了とします。
19. 次は雨域を時間方向に追跡する
Step 3が完成すると、各時刻にPolygonができます。
Step 4では、これを前後時刻でつなぎます。
比較する候補は、
Polygonの重なり
重心間距離
面積
最大雨量
長軸方向
です。
Step 3で特徴量を多めに残したのは、このStep 4へつなげるためでもあります。
まとめ
Step 3では、3時間積算雨量を単なるRasterのまま扱うのではなく、
雨量Grid
↓
強雨セル
↓
連結した雨域
↓
Polygon
↓
特徴量
へ変換します。
ここで重要なのは、最終判定だけを保存しないことです。
面積が不足したのか、線状性が不足したのか、最大雨量が不足したのかを後から確認できる形にしています。
また、このStepで付けるis_candidateは、気象庁の数値条件の一部を参考にした独自候補です。
キキクル等の条件を扱っていないため、公式の線状降水帯判定とは分けて扱います。
次は、時刻ごとに抽出した強雨域を追跡します。
参考資料
-
気象庁 線状降水帯に関する情報
https://www.jma.go.jp/jma/kishou/know/bosai/kishojoho_senjoukousuitai.html -
気象庁 線状降水帯の事例
https://www.data.jma.go.jp/senjo_list/list_senjoukousuitai.html -
SciPy
ndimage.label
https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.label.html -
Rasterio Features
https://rasterio.readthedocs.io/en/stable/topics/features.html -
Shapely
https://shapely.readthedocs.io/ -
pyproj Geod
https://pyproj4.github.io/pyproj/stable/api/geod.html -
GeoPandas GeoParquet
https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.to_parquet.html