防災GISデータを全国1kmメッシュへ集計する ― 防災Web地図 Step 2

Step 1では、e-Statの3次メッシュ境界を全国分まとめ、共通の1kmメッシュをGeoParquetにしました。

今回は、そのメッシュへ洪水浸水想定区域と土砂災害警戒区域を重ねます。

最初は、公式データを読み込んで1kmメッシュ別の危険度だけ作ればよいと考えていました。ただ、集計結果だけ残すと、後で「どの区域がこの値になったのか」を追えません。土砂災害警戒区域のように1kmメッシュよりかなり細かなデータでは、元の形も必要です。

そこでStep 2では、次の2種類を作ります。

元の区域を残したGeoParquet
  詳細表示、原本との照合、再集計に使う

1kmメッシュへ集計したGeoParquet
  全国表示、PMTiles、DuckDBの集計に使う

Step 1と同様に、入力を確認してから加工し、GeoParquetを書き出した後に必ず読み戻し検査を行う流れです。Step 2では、その間に「空間結合」「交差geometry」「重複区域のunion」「面積率の計算」が加わります。

公式防災GISデータ
    ↓
安全に展開
    ↓
ファイル・レイヤー・属性を確認
    ↓
JSON設定に従って正規化
    ├── 元形状GeoParquet
    │
    └── Step 1の1kmメッシュと空間結合
    │      ↓
    │   交差geometry作成
    │      ↓
    │   重複区域をunion
    │      ↓
    │   1kmメッシュ単位で集計
    │      ↓
    │   1km集計GeoParquet
    │
    └────────────┐
                 ↓
             読み戻し検査
                 ↓
       manifest・catalog・検査報告

この記事は、次の5段階で防災Web地図を作るうちのStep 2です。

防災Web地図を作る

Step 1 全国1kmメッシュをGeoParquetへ変換
Step 2 防災区域を正規化し、1kmメッシュへ集計
Step 3 TippecanoeでPMTilesを作成
Step 4 NginxとOpenLayersで静的表示
Step 5 DuckDBとFastAPIで詳細・統計APIを追加

今回の処理範囲

Step 2では、地図表示用のタイルはまだ作りません。GeoParquetを検査済みの中間成果物として固め、Step 3へ渡すところまでです。

使用するデータ

実装例には、国土数値情報の次の2データを入れました。

洪水浸水想定区域（河川単位）
土砂災害警戒区域

洪水浸水想定区域

2026年8月に確認した時点では、河川単位の洪水浸水想定区域は2025年度版が最新です。2026年5月に追加データが公開されています。

データは、次の5カテゴリに分かれています。

計画規模
想定最大規模
浸水継続時間
家屋倒壊等氾濫想定区域（氾濫流）
家屋倒壊等氾濫想定区域（河岸浸食）

今回は、想定最大規模の浸水深ランクを使います。

Shapefileの主な属性は次です。

属性	内容
`A31a_201`	河川番号
`A31a_202`	河川名
`A31a_203`	河川管理者番号
`A31a_204`	河川管理者
`A31a_205`	浸水深ランク

計画規模を追加するときは、同じテーブルへ混ぜず、別のscenario_idで扱います。

dataset_id   flood_a31a_2025
scenario_id  maximum

dataset_id   flood_a31a_2025
scenario_id  planned

URLやファイル名だけで区別せず、成果物の中にもシナリオを残します。

洪水浸水想定区域の説明には、前提となる降雨を超える雨、支川、高潮、内水などによって、区域外でも浸水する可能性があると記載されています。集計結果がない場所を「安全」と表示しないことも、Step 2のデータ契約に含めました。

土砂災害警戒区域

土砂災害警戒区域も、2025年度版が最新です。データ基準日は2025年8月1日で、2026年5月に更新されています。

主な属性は次です。

属性	内容
`A33_001`	現象の種類
`A33_002`	区域区分
`A33_003`	都道府県コード
`A33_004`	区域番号
`A33_005`	区域名
`A33_006`	所在地
`A33_007`	告示日
`A33_008`	特別警戒未指定フラグ

現象の種類は、急傾斜地の崩壊、土石流、地すべりです。

区域区分には、指定済みと指定前があります。

1  土砂災害警戒区域（指定済）
2  土砂災害特別警戒区域（指定済）
3  土砂災害警戒区域（指定前）
4  土砂災害特別警戒区域（指定前）

今回は、指定済みの1と2だけを集計します。

また、このデータには面と線が含まれます。面積率を求めるStep 2では面レイヤーだけを使い、線レイヤーは検査記録へ残したうえで除外します。

都道府県によって利用条件が異なる場合があります。公開ページに記載された各都道府県の条件とデータ時点は、ダウンロード時に必ず確認します。

1kmメッシュだけを残さない

防災区域を1kmメッシュへ集計すると、広域表示はかなり扱いやすくなります。

一方、1kmメッシュだけでは分からなくなる情報もあります。

区域がメッシュのどこにあるか
住宅地側か山側か
狭い特別警戒区域が含まれているか
複数河川の区域がどこで重なっているか

そのため、出力を次のように分けました。

workspace/artifacts/step2/
├─ native/
│  ├─ flood_a31a_2025__maximum.parquet
│  └─ landslide_a33_2025__designated.parquet
│
├─ summary_1km/
│  ├─ flood_a31a_2025__maximum.parquet
│  └─ landslide_a33_2025__designated.parquet
│
├─ dataset_catalog_fragment.json
├─ manifest.json
└─ validation_report.json

nativeには元のポリゴンを残します。

summary_1kmには、最大クラス、影響面積、面積率などを保存します。

実行環境

Docker Composeで実行します。

参考のためにZIPファイルを置いときますので、ダウンロードして使ってください。

Python 3.13
GeoPandas
Pyogrio
PyArrow
Shapely 2
PyProj

処理はbuilderコンテナだけで行います。Step 2の段階では、FastAPIやNginxは起動しません。

services:
  builder:
    build:
      context: .
      dockerfile: builder/Dockerfile
    working_dir: /app
    environment:
      PYTHONPATH: /app/src
      TZ: Asia/Tokyo
    volumes:
      - ./:/app

プロジェクト全体を/appへマウントするため、設定やPythonソースを修正しても、イメージを毎回作り直す必要はありません。依存ライブラリを変えた場合だけ再ビルドします。

フォルダー構成

japan_disaster_pmtiles_step2_r001/
├─ README.md
├─ VALIDATION.md
├─ japan_disaster_pmtiles_step2_qiita_r001.md
├─ docker-compose.yml
├─ Makefile
├─ .env.example
├─ pyproject.toml
│
├─ builder/
│  ├─ Dockerfile
│  └─ requirements.txt
│
├─ config/
│  ├─ hazard_layers.demo.json
│  └─ hazard_layers.example.json
│
├─ contracts/
│  ├─ artifact_manifest.schema.json
│  ├─ dataset_catalog_fragment.schema.json
│  ├─ hazard_layers.schema.json
│  └─ validation_report.schema.json
│
├─ src/hazard_step2/
│  ├─ aggregation.py
│  ├─ archives.py
│  ├─ catalog.py
│  ├─ cli.py
│  ├─ config.py
│  ├─ demo.py
│  ├─ errors.py
│  ├─ hashing.py
│  ├─ inspect_output.py
│  ├─ json_io.py
│  ├─ logging_utils.py
│  ├─ mesh.py
│  ├─ metadata.py
│  ├─ normalization.py
│  ├─ pipeline.py
│  ├─ sources.py
│  ├─ validation.py
│  └─ writer.py
│
├─ tests/
├─ docs/
├─ scripts/
└─ workspace/
   ├─ artifacts/step1/
   ├─ artifacts/step2/
   ├─ downloads/step2/
   ├─ raw/hazards/
   ├─ demo/
   └─ logs/

Step 1の成果物は次へ配置します。

workspace/artifacts/step1/
├─ mesh_master_1km.parquet
└─ manifest.json

Step 2はmanifest.jsonのstepとstatusを確認します。

step    step1
status  passed

検査に合格していないメッシュを、そのまま空間集計へ使わないためです。

属性名はJSONで対応付ける

防災データは、種類によって属性名がかなり違います。

洪水では河川番号や浸水深ランク、土砂災害では現象種別や区域区分を使います。これをPythonへ直接書くと、データを1つ増やすたびに処理本体を変更することになります。

そこで、属性対応をJSONへ分けました。

洪水の例です。

{
  "dataset_id": "flood_a31a_2025",
  "dataset_name": "洪水浸水想定区域（想定最大規模）",
  "hazard_type": "flood",
  "scenario_id": "maximum",
  "source_globs": [
    "**/20_想定最大規模/**/*.shp"
  ],
  "field_mapping": {
    "source_record_id": null,
    "source_code": "A31a_201",
    "source_name": "A31a_202",
    "admin_code": "A31a_203",
    "admin_name": "A31a_204",
    "hazard_class": "A31a_205",
    "hazard_value": null
  },
  "class_mapping": {
    "1": 1,
    "2": 2,
    "3": 3,
    "4": 4,
    "5": 5,
    "6": 6
  }
}

土砂災害の例では、指定前区域を属性フィルターで外します。

{
  "dataset_id": "landslide_a33_2025",
  "hazard_type": "landslide",
  "scenario_id": "designated",
  "source_globs": [
    "**/A33*.shp"
  ],
  "field_mapping": {
    "source_record_id": "A33_004",
    "source_code": "A33_001",
    "source_name": "A33_005",
    "admin_code": "A33_003",
    "hazard_class": "A33_002"
  },
  "attribute_filters": {
    "A33_002": ["1", "2"]
  },
  "class_mapping": {
    "1": 1,
    "2": 2
  }
}

配布データのフォルダー名や属性名が変わった場合は、PythonではなくこのJSONを直します。

ただし、設定例をそのまま信用せず、実際に取得したファイルの属性を先に確認します。

make inspect-source

結果は次へ保存します。

workspace/logs/step2_source_inspection.json

ここで確認するのは次です。

globに一致したファイル
GDALレイヤー名
CRS
geometry型
属性列
地物数
入力ファイル一式のSHA-256

全国処理を始めてから属性不足で止まるより、先に入力構成だけを見る方が修正しやすくなります。

ZIPを安全に展開する

ダウンロードしたZIPは次へ置きます。

workspace/downloads/step2/

展開は次です。

make prepare

archives.pyでは、ZIP内のパスを検査してから1件ずつ展開します。

def _validate_member_name(member_name: str) -> PurePosixPath:
    """ZIP内のパスが展開先の外を指していないか確認する。"""

    normalized = member_name.replace("\\", "/")
    member_path = PurePosixPath(normalized)

    if member_path.is_absolute():
        raise SourceDataError(
            f"ZIP内に絶対パスがあります: {member_name!r}"
        )

    if any(part in {"", ".", ".."} for part in member_path.parts):
        raise SourceDataError(
            f"ZIP内に不正な相対パスがあります: {member_name!r}"
        )

    return member_path

途中で失敗した場合は、半端に展開したフォルダーを削除します。

展開途中のデータを完全な原本と誤認しないためです。

ファイルとレイヤーを事前検査する

Shapefileは、少なくとも次の構成を確認します。

.shp
.shx
.dbf

.prjがない場合はCRSを取得できないため、通常は処理を止めます。

どうしても原本にCRS情報がない場合だけ、設定のsource_crs_overrideを使います。ただし、これは推測で指定する項目ではありません。原典資料で座標系を確認できた場合に限ります。

GeoPackageやGMLのように複数レイヤーを持つ形式では、source_layerを指定します。

暗黙に最初のレイヤーを選ぶと、面ではなく線を読んでしまうことがあるためです。

共通スキーマへ正規化する

元の区域は、次の列へそろえます。

feature_id
hazard_id
dataset_id
hazard_type
scenario_id
hazard_class
hazard_class_label
hazard_value
source_record_id
source_code
source_name
admin_code
admin_name
source_file
source_layer
extra_attributes_json
geometry
bbox

hazard_idは、データセットID、シナリオ、元ファイル、レイヤー、元レコードIDからSHA-256を作り、先頭24文字を使います。

def _build_hazard_id(
    dataset: DatasetConfig,
    source: HazardSource,
    source_record_id: str,
    source_position: int,
) -> str:
    """原本レコードから再現可能な短い識別子を作る。"""

    material = "|".join(
        (
            dataset.dataset_id,
            dataset.scenario_id,
            source.relative_path,
            source.layer,
            source_record_id,
            str(source_position),
        )
    )
    return hashlib.sha256(
        material.encode("utf-8")
    ).hexdigest()[:24]

同じ入力と設定から作り直したとき、同じ地物に同じIDを付けるためです。

コード値の表記をそろえる

DBFやGeoJSONからコードを読むと、同じ値でも次のように見える場合があります。

2
2.0
"2"

危険度対応表と照合する前に、表記だけをそろえます。

def normalize_code_key(value: Any) -> str:
    """コード対応表で使う比較用文字列を作る。"""

    if value is None or pd.isna(value):
        return ""

    text = str(value).strip()

    if text.endswith(".0"):
        prefix = text[:-2]
        if prefix and prefix.lstrip("+-").isdigit():
            return prefix

    return text

先頭ゼロには意味がある可能性があるため、整数へ変換して削りません。

geometryを面として整える

土砂災害データには線レイヤーもあります。また、元ポリゴンに自己交差などの不正geometryが含まれる場合もあります。

Step 2では次の順で処理します。

欠損・空geometry
  除外して件数を記録

線・点
  面積集計の対象外として除外

不正Polygon
  repair_invalid=trueならmake_valid

GeometryCollection
  PolygonとMultiPolygonだけを抽出

修復後も不正
  エラーとして停止

中心となる処理です。

def normalize_polygon_geometry(
    geometry: BaseGeometry | None,
    *,
    repair_invalid: bool,
) -> tuple[BaseGeometry | None, bool, str | None]:
    """防災区域geometryを有効なPolygon系へ整える。"""

    if geometry is None or geometry.is_empty:
        return None, False, "empty"

    repaired = False
    candidate = geometry

    if not candidate.is_valid:
        if not repair_invalid:
            raise SourceDataError(
                "不正geometryがあります。"
            )

        candidate = make_valid(candidate)
        repaired = True

    polygonal = _polygonal_part(candidate)

    if polygonal is None or polygonal.is_empty:
        return None, repaired, "non_polygon"

    if not polygonal.is_valid:
        raise SourceDataError(
            "make_valid後も不正なPolygonが残りました。"
        )

    return polygonal, repaired, None

修復件数、空geometryの除外件数、線・点の除外件数はmanifestへ残します。

自動修復したことを隠さないためです。

第1次メッシュ単位で空間結合する

全国分の防災区域と1kmメッシュを一度にoverlayすると、交差後の一時データが大きくなります。

Step 1のメッシュにはmesh1_codeがあるため、第1次メッシュ単位に分けて処理します。

第1次メッシュを1つ取り出す
  ↓
その範囲と交差する防災区域候補を空間索引で取得
  ↓
GeoPandas sjoin
  ↓
実際のintersectionを配列演算
  ↓
面積を持つ交差だけ残す

候補検索と実際のintersectionを分けているのは、境界ボックスだけで候補を絞り、必要な組だけgeometry計算するためです。

joined = gpd.sjoin(
    mesh_chunk,
    hazard_chunk,
    how="inner",
    predicate="intersects",
)

intersections = shapely.intersection(
    joined.geometry.array,
    hazard_geometries,
)

境界に触れただけの地物を除外する

intersectsは、ポリゴン同士が線や点で触れただけでも真になります。

面積集計では、交差面積が設定値以下の組を除外します。

{
  "min_intersection_area_m2": 0.01
}

positive_mask = (
    (~shapely.is_empty(intersections))
    & np.isfinite(intersection_areas)
    & (
        intersection_areas
        > dataset.min_intersection_area_m2
    )
)

0.01平方メートルは危険区域の意味上の閾値ではなく、境界接触や計算誤差を面積として数えないための初期値です。

重複する区域を二重計上しない

河川別の浸水区域をまとめると、同じ1kmメッシュ内でポリゴンが重なることがあります。

各交差面積をそのまま合計すると、次のような結果になります。

メッシュ面積          1.0 km²
河川Aとの交差面積     0.7 km²
河川Bとの交差面積     0.6 km²
単純合計              1.3 km²

この1.3km²は、重なった場所を2回数えています。

今回は、メッシュ内の交差geometryをunionしてから面積を求めます。

for mesh_row, group in pairs.groupby(
    "mesh_row",
    sort=False,
):
    union_geometry = shapely.union_all(
        group["intersection_geometry"].to_numpy()
    )

    affected_area = float(
        shapely.area(union_geometry)
    )

これにより、影響面積率を0～1の範囲に保ちやすくなります。

一方、source_feature_countには交差した元地物数を残します。面積は重複を除きますが、複数河川や複数区域が関係していることは確認できます。

面積計算にEPSG:3857を使わない

Web地図はEPSG:3857で表示しますが、面積計算には使いません。

初期設定では、全国を同じ投影条件で概略集計するため、EPSG:6933を使います。

{
  "target_crs": "EPSG:4326",
  "area_crs": "EPSG:6933"
}

target_crsは保存用、area_crsは面積計算用です。

保存
  EPSG:4326

面積計算
  EPSG:6933

Step 3のタイル化
  Tippecanoe側でWeb Mercatorタイルへ変換

EPSG:6933は全国を同じ条件で比較するための選択です。測量成果や法的境界を確定する用途ではありません。

都道府県、流域、市町村など対象範囲が決まっていて、より厳密な面積が必要な場合は、その地域に合うJGD2011平面直角座標系へ変更して比較します。

1km集計に保存する値

feature_id
mesh_code
dataset_id
hazard_type
scenario_id
mesh1_code
has_hazard
max_class
max_class_label
max_value
affected_area_m2
mesh_area_m2
affected_area_ratio
source_feature_count
intersection_count
geometry
bbox

max_classは、メッシュと交差した区域の最大危険度です。

affected_area_ratioは次で求めます。

affected_area_ratio
  = union後の影響面積
  ÷ 1kmメッシュ面積

デフォルトでは、正の面積で交差したメッシュだけを保存します。

{
  "include_unaffected_meshes": false
}

全メッシュを出力する場合はtrueへ変更できます。ただし、未交差を安全とは表記せず、集計対象区域との交差なしとします。

GeoParquetとして保存する

元形状と集計結果は、どちらも同じ条件で保存します。

GeoParquet 1.1.0
geometry encoding  WKB
圧縮               ZSTD
bbox covering      あり

書き込み途中のファイルを成果物として残さないよう、一時ファイルへ保存してから置き換えます。

def write_geoparquet_atomic(
    frame: gpd.GeoDataFrame,
    output_path: Path,
    *,
    row_group_size: int,
    overwrite: bool,
) -> None:
    """GeoParquet 1.1を途中状態を残さずに保存する。"""

    descriptor, temporary_name = tempfile.mkstemp(
        prefix=f".{output_path.stem}.",
        suffix=".tmp.parquet",
        dir=output_path.parent,
    )
    os.close(descriptor)
    temporary_path = Path(temporary_name)

    try:
        frame.to_parquet(
            temporary_path,
            index=False,
            compression="zstd",
            compression_level=6,
            geometry_encoding="WKB",
            write_covering_bbox=True,
            schema_version="1.1.0",
            row_group_size=row_group_size,
        )
        os.replace(temporary_path, output_path)
    finally:
        temporary_path.unlink(missing_ok=True)

既存成果物がある状態では、--overwriteを付けない限り停止します。

まずデモを実行する

実データを準備する前に、架空データで処理を確認します。

cp .env.example .env
make demo

デモでは次を作ります。

1kmメッシュ                30セル
架空の洪水区域             3ポリゴン
架空の土砂災害警戒区域     3ポリゴン

洪水区域には、次の確認ケースを入れています。

2区域が重なるケース
自己交差Polygonをmake_validするケース
複数の危険度クラス

土砂災害区域には、区域コード3の指定前区域を1件入れています。設定の属性フィルターで除外されることを確認します。

デモ成果物は次です。

workspace/artifacts/step2_demo/
├─ native/
├─ summary_1km/
├─ dataset_catalog_fragment.json
├─ manifest.json
└─ validation_report.json

実在する河川や区域を表したものではありません。地図表示や防災判断には使いません。

実データを配置する

Step 1成果物を配置します。

workspace/artifacts/step1/
├─ mesh_master_1km.parquet
└─ manifest.json

国土数値情報から取得したZIPは、展開せずに次へ置きます。

workspace/downloads/step2/
├─ A31a-25_..._SHP.zip
├─ A33-25_..._SHP.zip
└─ ...

実行順です。

make build
make prepare
make inspect-source
make step2
make validate
make inspect
make test

1. builderイメージを作る

make build

2. ZIPを展開する

make prepare

既存の展開先を作り直す場合です。

./scripts/prepare_archives.sh --overwrite

3. 入力構成を確認する

make inspect-source

ここで属性名や対象レイヤーが合わない場合は、config/hazard_layers.example.jsonを修正します。

4. Step 2を実行する

make step2

既存成果物を置き換える場合だけ次を使います。

./scripts/run_step2.sh --overwrite

5. 読み戻し検査を再実行する

make validate

6. manifestの要点を見る

make inspect

成果物を読み戻して確認する

validation_report.jsonでは、次を検査します。

native GeoParquet

必須列
GeoParquet 1.1
WKB
bbox covering
CRS
feature_idとhazard_idの一意性
hazard_classが1以上
geometryの欠損、空、不正
PolygonまたはMultiPolygon

summary_1km GeoParquet

必須列
mesh_codeが8桁
mesh_codeとfeature_idの一意性
Step 1に存在するmesh_codeだけであること
影響面積率が0～1
影響面積がメッシュ面積以下
has_hazardとmax_classの対応

検査に失敗した場合、manifestのstatusはfailedになります。Step 3ではpassedだけを入力として扱います。

manifestへ残す情報

Step 1メッシュのSHA-256
Step 1 manifestのSHA-256
設定JSONのSHA-256
入力GISファイル一式のSHA-256
対象・除外レイヤー
正規化件数
geometry修復件数
属性フィルター除外件数
影響メッシュ数
候補交差数
正の面積を持つ交差数
成果物のSHA-256
実行環境とライブラリ版
処理時間

元のZIPを更新したときに、どの入力から作った結果か追えるようにしています。

データセットカタログ

dataset_catalog_fragment.jsonは、Step 3とStep 5へ渡すための一覧です。

{
  "dataset_id": "flood_a31a_2025",
  "dataset_name": "洪水浸水想定区域（想定最大規模）",
  "hazard_type": "flood",
  "scenario_id": "maximum",
  "source_version": "2025年度版",
  "reference_date": "2025年度",
  "target_crs": "EPSG:4326",
  "area_crs": "EPSG:6933",
  "artifacts": {
    "native": {
      "path": "...",
      "sha256": "...",
      "row_count": 0
    },
    "summary_1km": {
      "path": "...",
      "sha256": "...",
      "row_count": 0
    }
  }
}

画面の凡例やポップアップで、出典、版、基準日、利用条件を表示するための元情報にも使います。

テスト

make test

主なテストは次です。

JSON設定の矛盾検出
ZIP Slipの拒否
コード値2.0と2の表記統一
自己交差Polygonの修復
線geometryの除外
重複区域の面積二重計上防止
GeoJSONの探索と共通スキーマ変換
デモ全体のGeoParquet統合処理

重複区域のテストでは、1km²のメッシュを2つのポリゴンで重ねます。

assert row["affected_area_m2"] == pytest.approx(
    1_000_000.0
)
assert row["affected_area_ratio"] == pytest.approx(1.0)
assert row["source_feature_count"] == 2

単純合計なら1km²を超える配置ですが、union後の面積は1km²になります。

防災データとしての注意

1kmメッシュ集計は、全国を同じ単位で見比べるための派生データです。

元の指定区域を置き換えるものではありません。

全国・都道府県の傾向を見る
  1kmメッシュ集計

町丁目・街区で区域を確認する
  native GeoParquet

法的な境界や最新指定を確認する
  所管機関の公示図書・最新資料

特に土砂災害警戒区域は、国土数値情報の公開ページでも、概略位置を示す参考図であり、法的図書ではないことが注意事項として示されています。

洪水も、想定条件を超える現象や内水などによって、区域外で浸水する可能性があります。

Web画面では「該当なし」を「安全」と読み替えない表現にします。

Step 2でできたもの

今回の処理で、次の役割が分かれました。

Step 1
  全国共通の1kmメッシュ

Step 2 native
  元の防災区域

Step 2 summary_1km
  全国表示用の1km集計

catalog / manifest / validation
  出典、条件、検査、再現情報

次のStep 3では、summary_1kmと必要に応じてnativeをFlatGeobufへ変換し、TippecanoeでPMTilesを作ります。

低ズームで1kmメッシュを無理に全部描画せず、10kmや20kmへ集計したレイヤーも追加できる構成にしていきます。

参考