線状降水帯予測プロトタイプを作ってみる【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