# -*- coding: utf-8 -*- """L45 ため池浸水想定区域 重い空間処理を事前にキャッシュするスクリプト。 本記事 .py 実行時はこのキャッシュ CSV/GPKG を読み込むだけで完走する (要件 S)。 入力: data/extras/L45_pond_inundation/shp/ため池浸水想定区域.shp (143MB, 7223 polys) data/extras/tameike_basic.csv (1.1MB, 6754 rows) data/extras/L15_admin_zones/admin_922_広島県.zip (都市計画区域) 生成ファイル: data/extras/L45_pond_inundation/_cache/diss_per_pond.gpkg ため池ごとに dissolve した浸水想定 (6730 features) data/extras/L45_pond_inundation/_cache/pond_inund_merge.csv ため池属性 + 浸水想定面積 結合表 (6754 rows) data/extras/L45_pond_inundation/_cache/inund_union.gpkg 浸水想定区域全体 unary_union (面積二重カウント排除用) data/extras/L45_pond_inundation/_cache/per_city_stats.csv 市町別 集計 data/extras/L45_pond_inundation/_cache/admin_diss.gpkg 都市計画区域行政区域 (L15 由来) """ from __future__ import annotations import sys, time, zipfile from pathlib import Path ROOT = Path(__file__).resolve().parents[1] DEST = ROOT / "data" / "extras" / "L45_pond_inundation" CACHE = DEST / "_cache" CACHE.mkdir(parents=True, exist_ok=True) import geopandas as gpd import pandas as pd import shapely TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標系第 III 系 (m 単位) SHP_PATH = DEST / "shp" / "ため池浸水想定区域.shp" CSV_PATH = ROOT / "data" / "extras" / "tameike_basic.csv" print("=== L45 cache builder ===", flush=True) # ----------------------------------------------------------------------------- # STEP 1: ため池ごと dissolve (FIELD001 列で集約) # ----------------------------------------------------------------------------- diss_gpkg = CACHE / "diss_per_pond.gpkg" if diss_gpkg.exists(): diss = gpd.read_file(diss_gpkg) print(f" [diss] cache hit: {len(diss)} ため池ポリゴン", flush=True) else: t = time.time() raw = gpd.read_file(SHP_PATH, encoding="cp932") print(f" [diss] read: {len(raw)} polys, {time.time()-t:.1f}s", flush=True) t = time.time() raw = raw.to_crs(TARGET_CRS) # 不正幾何の修復 geoms = raw.geometry.values inv = ~shapely.is_valid(geoms) if inv.sum() > 0: new = [g.buffer(0) if isi else g for g, isi in zip(geoms, inv)] raw = raw.copy() raw["geometry"] = new print(f" [diss] fixed {inv.sum()} invalid polys", flush=True) diss = raw.dissolve(by="FIELD001", as_index=False) diss["area_km2"] = diss.geometry.area / 1e6 # ポリゴン重心 (代表点) で各ため池の地理座標 diss["centroid_x"] = diss.geometry.centroid.x diss["centroid_y"] = diss.geometry.centroid.y print(f" [diss] dissolve: {len(diss)} unique pond polys, {time.time()-t:.1f}s", flush=True) diss.to_file(diss_gpkg, driver="GPKG") print(f" [diss] saved {diss_gpkg.name}", flush=True) # ----------------------------------------------------------------------------- # STEP 2: ため池基本情報 (CSV) と浸水想定面積を結合 # ----------------------------------------------------------------------------- merge_csv = CACHE / "pond_inund_merge.csv" if merge_csv.exists(): merge = pd.read_csv(merge_csv, dtype={"ため池番号": str, "FIELD001": str}) print(f" [merge] cache hit: {len(merge)} rows", flush=True) else: t = time.time() df = pd.read_csv(CSV_PATH, dtype=str) df.columns = [c.strip() for c in df.columns] df["ため池番号"] = df["ため池番号"].str.strip() df["堤高_m"] = pd.to_numeric(df["堤高(m)"] if "堤高(m)" in df.columns else df["堤高（m）"], errors="coerce") df["貯水量_千m3"] = pd.to_numeric(df["貯水量(千m3)"] if "貯水量(千m3)" in df.columns else df["貯水量（千m3）"], errors="coerce") df["緯度"] = pd.to_numeric(df["緯度"], errors="coerce") df["経度"] = pd.to_numeric(df["経度"], errors="coerce") # FIELD001 = ため池番号 で結合 merge = df.merge(diss[["FIELD001", "area_km2", "centroid_x", "centroid_y"]], left_on="ため池番号", right_on="FIELD001", how="left") merge["has_inund"] = merge["area_km2"].notna() merge.to_csv(merge_csv, index=False, encoding="utf-8-sig") print(f" [merge] {len(merge)} rows ({merge['has_inund'].sum()} with SHP, " f"{(~merge['has_inund']).sum()} without SHP), {time.time()-t:.1f}s", flush=True) # ----------------------------------------------------------------------------- # STEP 3: 浸水想定区域 全体 unary_union (面積の二重カウント排除) # ----------------------------------------------------------------------------- union_gpkg = CACHE / "inund_union.gpkg" if union_gpkg.exists(): union_gdf = gpd.read_file(union_gpkg) print(f" [union] cache hit", flush=True) else: t = time.time() # 全 7223 polygon の unary_union if 'raw' not in dir(): raw = gpd.read_file(SHP_PATH, encoding="cp932").to_crs(TARGET_CRS) # fix invalid geoms = raw.geometry.values inv = ~shapely.is_valid(geoms) if inv.sum() > 0: new = [g.buffer(0) if isi else g for g, isi in zip(geoms, inv)] raw = raw.copy() raw["geometry"] = new union_geom = shapely.unary_union(raw.geometry.values) union_gdf = gpd.GeoDataFrame({"name": ["pond_inund_all"]}, geometry=[union_geom], crs=TARGET_CRS) union_gdf["area_km2"] = union_gdf.geometry.area / 1e6 union_gdf.to_file(union_gpkg, driver="GPKG") print(f" [union] saved, total area = {union_gdf['area_km2'].iloc[0]:.2f} km², " f"{time.time()-t:.1f}s", flush=True) # ----------------------------------------------------------------------------- # STEP 4: 都市計画区域行政区域 (L15) と CITY_CD で結合 # ----------------------------------------------------------------------------- admin_diss_gpkg = CACHE / "admin_diss.gpkg" ADMIN_PREF = ROOT / "data" / "extras" / "L15_admin_zones" / "admin_922_広島県.zip" if admin_diss_gpkg.exists(): admin_diss = gpd.read_file(admin_diss_gpkg) print(f" [admin] cache hit: {len(admin_diss)} cities", flush=True) else: with zipfile.ZipFile(ADMIN_PREF) as zf: name = [n for n in zf.namelist() if n.endswith(".geojson")][0] admin = gpd.read_file(f"zip://{ADMIN_PREF}!{name}").to_crs(TARGET_CRS) admin_diss = admin.dissolve(by="CITY_CD", as_index=False) admin_diss["admin_km2"] = admin_diss.geometry.area / 1e6 admin_diss.to_file(admin_diss_gpkg, driver="GPKG") print(f" [admin] saved: {len(admin_diss)} cities", flush=True) # ----------------------------------------------------------------------------- # STEP 5: 市町別統計 (所管市町列ベースで集計) # ----------------------------------------------------------------------------- city_csv = CACHE / "per_city_stats.csv" if city_csv.exists(): city_stats = pd.read_csv(city_csv, encoding="utf-8-sig") print(f" [city] cache hit: {len(city_stats)} 市町", flush=True) else: t = time.time() g = merge.groupby("所管市町").agg( ため池数=("ため池番号", "count"), 浸水想定あり=("has_inund", "sum"), 平均堤高_m=("堤高_m", "mean"), 最大堤高_m=("堤高_m", "max"), 平均貯水量_千m3=("貯水量_千m3", "mean"), 合計貯水量_千m3=("貯水量_千m3", "sum"), 合計浸水想定_km2=("area_km2", "sum"), 平均浸水想定_km2=("area_km2", "mean"), 中央緯度=("緯度", "median"), 中央経度=("経度", "median"), ).reset_index() city_stats = g.sort_values("ため池数", ascending=False) city_stats.to_csv(city_csv, index=False, encoding="utf-8-sig") print(f" [city] {len(city_stats)} 市町, {time.time()-t:.1f}s", flush=True) # ----------------------------------------------------------------------------- # STEP 6: 県全体集計 (重複排除済みの真の浸水想定面積) # ----------------------------------------------------------------------------- print() print("=== サマリ ===") total_pond = len(merge) n_with_shp = merge["has_inund"].sum() sum_indiv = merge["area_km2"].sum() union_area = union_gdf["area_km2"].iloc[0] overlap_ratio = (sum_indiv - union_area) / sum_indiv * 100 if sum_indiv > 0 else 0 print(f" ため池総数: {total_pond}") print(f" 浸水想定 SHP あり: {n_with_shp} ({n_with_shp/total_pond*100:.1f}%)") print(f" 浸水想定 SHP 無し: {total_pond - n_with_shp}") print(f" 個別浸水想定面積の合計: {sum_indiv:.2f} km² (重複あり)") print(f" unary_union 後の面積: {union_area:.2f} km² (真の浸水想定面積)") print(f" 重複面積 (overlap): {sum_indiv - union_area:.2f} km² ({overlap_ratio:.1f}%)") print(f" 所管市町数: {merge['所管市町'].nunique()}") print("=== done ===")