防災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へ集計したレイヤーも追加できる構成にしていきます。
参考
- 国土数値情報 洪水浸水想定区域データ（河川単位）2025年度版
- 国土数値情報 土砂災害警戒区域データ 2025年度版
- 国土数値情報（土砂災害警戒区域）製品仕様書
- GeoPandas spatial join
- Shapely unary_union / union_all
- GeoParquet 1.1 specification