2025年12月 – Linux & Android Dialy

時間到達圏（isochrone）生成関数から作成

G2（浸水回避済みグラフ）と自宅座標から徒歩Y分の到達圏ポリゴンを作る

到達圏は「到達可能エッジを太らせて面化」する堅牢方式
（alpha-shape不要、Shapely/GeoPandas 標準で完結）

前提：G2（到達可能道路グラフ）、home_wgs、GPKG、METRIC=”EPSG:32654″ が既に定義済み

# === Isochrone 生成ユーティリティ ===
import math
import geopandas as gpd
import networkx as nx
import osmnx as ox
from shapely.ops import unary_union

def add_edge_time_weight(G, speed_kmh: float = 4.8, attr_name: str = "time_min"):
    """
    エッジに 'time_min'（分）を付与（長さ[m] / 速度[km/h]）.
    既にある場合は上書き。
    """
    fac = speed_kmh * 1000 / 60.0  # m / min
    for u, v, k, data in G.edges(keys=True, data=True):
        length_m = float(data.get("length", 0.0))
        data[attr_name] = length_m / fac
    return G

def nearest_node_from_wgs(G, lat: float, lon: float):
    """WGS84の(lat,lon)から最近傍ノードIDを返す（osmnxのAPI差異に両対応）。"""
    try:
        return ox.distance.nearest_nodes(G, X=lon, Y=lat)
    except AttributeError:
        return ox.nearest_nodes(G, X=lon, Y=lat)

def isochrone_polygon(
    G,
    origin_node: int,
    cutoff_min: float,
    proj_crs: str = "EPSG:32654",
    buffer_m: float = 35.0,
    weight_attr: str = "time_min",
):
    """
    到達可能サブグラフのエッジを buffer で面化 → マルチポリゴン化。
    返り値: (isochrone_gdf[EPSG:4326], reachable_edges_gdf[EPSG:4326])

    - cutoff_min: 到達時間しきい値（分）
    - buffer_m: エッジに与えるバッファ[m]（道路幅+余裕）
    """
    # dijkstra で到達ノード
    dist = nx.single_source_dijkstra_path_length(G, origin_node, cutoff=cutoff_min, weight=weight_attr)
    reachable_nodes = set(dist.keys())
    if not reachable_nodes:
        raise RuntimeError("到達ノードが空です。cutoff_min が小さすぎる可能性。")

    # 到達エッジ抽出
    nodes_gdf, edges_gdf = ox.graph_to_gdfs(G, nodes=True, edges=True)
    idx = edges_gdf.reset_index()[["u", "v", "key"]].apply(tuple, axis=1)
    mask = [(u in reachable_nodes and v in reachable_nodes) for (u, v, k) in idx]
    reach_edges = edges_gdf.loc[mask].copy()
    if reach_edges.empty:
        raise RuntimeError("到達エッジが空です。グラフが疎か cutoff_min が小さい可能性。")

    # 面化（project→buffer→union→修正→逆投影）
    reach_proj = reach_edges.to_crs(proj_crs)
    merged = unary_union(reach_proj.buffer(buffer_m).geometry.values)
    # geometry 修正（buffer(0) に相当）
    poly = gpd.GeoSeries([merged], crs=proj_crs).buffer(0).unary_union

    # 出力（4326）
    iso_gdf = gpd.GeoDataFrame({"mode": ["walk"], "minutes": [float(cutoff_min)]},
                               geometry=[poly], crs=proj_crs).to_crs(4326)
    reach_edges_4326 = reach_proj.to_crs(4326)
    return iso_gdf, reach_edges_4326

def subtract_flood(iso_gdf: gpd.GeoDataFrame, flood_union_gdf: gpd.GeoDataFrame, proj_crs: str = "EPSG:32654"):
    """
    isochrone（面）から津波浸水（union）を差し引き → 「到達可能 × 非浸水」エリア。
    返り値: safe_area_now_gdf[EPSG:4326]
    """
    iso_m = iso_gdf.to_crs(proj_crs)
    fl_m  = flood_union_gdf.to_crs(proj_crs)
    safe_geom = iso_m.geometry.iloc[0].difference(fl_m.unary_union)
    safe = gpd.GeoDataFrame({"minutes":[iso_gdf["minutes"].iloc[0]]}, geometry=[safe_geom], crs=proj_crs).to_crs(4326)
    return safe

実行すると

/var/folders/db/tkl5w9dd3kn3h53ctyl4s6180000gn/T/ipykernel_72100/1981613481.py:59: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  poly = gpd.GeoSeries([merged], crs=proj_crs).buffer(0).unary_union
/var/folders/db/tkl5w9dd3kn3h53ctyl4s6180000gn/T/ipykernel_72100/1981613481.py:74: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  safe_geom = iso_m.geometry.iloc[0].difference(fl_m.unary_union)
<folium.features.GeoJson at 0x178279f10>

警告は Shapely/GeoPandas 2.x 系の非推奨によるもの
unary_union は廃止予定 → union_all() に置き換えればOK

変更すると

/var/folders/db/tkl5w9dd3kn3h53ctyl4s6180000gn/T/ipykernel_72100/832098648.py:59: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  poly = gpd.GeoSeries([merged], crs=proj_crs).buffer(0).unary_union
<folium.features.GeoJson at 0x17dfdf350>

    poly = gpd.GeoSeries([merged], crs=proj_crs).buffer(0).unary_union

も修正する

poly = gpd.GeoSeries([merged], crs=proj_crs).buffer(0).union_all()

merged = unary_union(reach_proj.buffer(buffer_m).geometry.values)

も修正

buf = reach_proj.buffer(buffer_m)
merged = gpd.GeoSeries(buf, crs=proj_crs).union_all()

に書き換える

これでエラーが消える

重みtime_minはエッジ長さ/歩行速度。最短経路と同様にDijkstraで「時間≦Y分」のノード集合を取り、両端が到達可能なエッジのみ抽出→buffer→unionで面化しています。
buffer_m は道路幅＋余裕のパラメータ。市街地なら 30–50m で見た目・連続性が良好です。
非浸水差引は isochrone – flood_union。あなたの既存 flood_union（津波浸水 union）をそのまま使えます。
保存は pyogrio + OVERWRITE に統一（毎回クリーン）。

次は

# ② 津波「到達時間・高さ」簡易推定モジュール

**目的**: 「あと何分で来る？」に即答。初期は簡易/安全側、後で精緻化。
**ステップ**

* v0（当日実用版）

* **深海近似の伝播速度** (c ≈ √(g·h)) を **一定深さ代表値**で近似（安全側に遅めに見積もる）
* 震源→沿岸/自宅最寄り海岸線までの距離から到達時間を算出
* 高さは情報源（JMA速報の区分）を採用、未入手なら「不明→避難優先」のメッセージ
* v1（精緻化）

* 沿岸ごとの**水深ラスタ**（GEBCO/J-EGG等）から距離積分で可変速度伝播
* 複数反射・地形収束は扱わず「防災判断に使える保守的な上限/下限」を出す
* 出力

* `arrival_time_min`, `height_class`, `confidence` を `route_latest`/`sim_runs` に書込

使用する用途のコンセプト

✅ あなたの用途の要件
* 精度：震源→沿岸までの到達時間・波高推定の入力に使うので、誤差が大きすぎると危険。
* 範囲：主に日本周辺で十分。全球は必須ではない。
* 更新頻度：将来的に OSM/防災機関の最新データと連携したい。
* GraphDB連携：避難所・道路・避難経路・浸水域をノード/リレーションで管理。
* 利用イメージ：
* 「地震が来たら最寄りの避難所はどこ？」 → 避難所ノード＋道路グラフ
* 「非常食の賞味期限が近いものは？」 → 在庫管理ノード
* 「津波到達までに何分ある？」 → 震源ノード→海岸ノードのコスト距離

✅ データセットごとの特徴
1. 国土数値情報「海岸線データ」
* 日本限定。行政データとして信頼性が高い。
* GML/Shapefile 提供。
* 避難計画（住民向け防災マップ）に最も適する。 ➡ 推奨（基準データにする）
2. OSM Coastlines / Land polygons
* 更新頻度が高い（週次〜月次）。
* 最新の人工改変（埋立地・護岸）も比較的反映。
* 避難経路や実際の生活圏との整合性が良い。 ➡ 補助データ（現況把握用）
3. Natural Earth
* 全球対応、軽量。
* 精度は粗い（1:10mでも数百 m〜km の誤差）。
* 可視化や概念実証には便利。 ➡ プロトタイプ用のみ
4. GSHHG
* 全球対応、高解像度あり。
* ファイルが大きく、取り回しに手間。
* 学術研究向き。 ➡ 専門シミュレーション以外では過剰

✅ 推奨構成（防災マップ＋GraphDB向け）
1. 海岸線：
* 基本 → 国土数値情報（海岸線）
* 最新の現況補正 → OSM Coastlines
2. 標高/地形：
* 陸域 DEM → 国土地理院 DEM / ALOS / SRTM
* 海域 DEM → GEBCO
3. 避難所・道路：
* 避難所 → 自治体オープンデータ（CSV/GeoJSON）
* 道路 → OSMnx で抽出
4. GraphDBに格納：
* ノード：震源・海岸・浸水域・避難所・家庭・在庫
* リレーション：到達時間・避難経路・所要時間・備蓄品管理

✅ まとめ
ご提示の用途（家族避難＋地域防災＋GraphDB管理）なら：
* **国土数値情報（海岸線）**をメイン
* OSM Coastlinesで現況更新
* DEM（陸域＋海域）＋道路＋避難所を組み合わせ
* GraphDBでノード/リレーションを整理して「質問に答えられる防災知識ベース」を構築

これらを考慮して次に作る機能は？

いいね、その前提なら“すぐ効く順”でこれを作るのが最短です。各項目は**完成定義(DoD)**と**最小コード断片**まで落とし込みました。

# 1) 震源→到達時間v0（安全側近似）

**目的**: 「あと何分？」に即答（まずは簡易モデル）。
**DoD**: `route_latest` に `arrival_time_min`, `height_class`, `confidence` を追記して保存。
**やること**

* 震源点(緯度経度)と自宅最寄り海岸線点の距離`D`を計算
* 代表水深`h_rep`で c≈√(g·h_rep)（保守的に**遅め**で見積り）
* `arrival_time_min = D / c / 60`（分）
* 波高は「不明なら ‘unknown’」、速報値が入手できたら `height_class` に格納

def arrival_time_min_km(d_km, h_rep=3000):  # m
    import math
    g=9.81; c=math.sqrt(g*h_rep) # m/s
    return (d_km*1000)/c/60

# 保存時に route_latest の属性へ追記
route_gdf["arrival_time_min"] = [arrival_time_min_km(D_km)]
route_gdf["height_class"] = ["unknown"]
route_gdf["confidence"] = ["low"]  # v0は保守的
```

# 2) 「安全エリア」ポリゴン生成（反転×等時間到達圏）

**目的**: 地図で「今向かうべき範囲」を1枚で示す。
**DoD**: `safe_area_now` と `isochrone_walk_min` を GPKG に**上書き**保存し、Foliumに重ねる。
**やること**

* 浸水想定を反転（`safe = study_area – tsunami_inundation`）
* 歩行速度4.8km/hで `Y = min(arrival_time*0.7, 20)` 分の等時間到達圏（graph isochrone）
* `safe_area_now = safe ∩ isochrone`

（等時間到達圏は既存グラフ`G2`から到達ノード抽出→凸包/マルチポリゴン化でOK）

—

# 3) 自動実行パイプライン（CLI化）

**目的**: トリガーが来たら無人で回る。
**DoD**: `python evac_pipeline.py –home 34.728,137.978 –quake 33.0,138.5 –radius 12000` で
`roads_safe`,`shelters_safe`,`route_latest`,`safe_area_now` が更新され、LINE/音声に通知。

**雛形**

# evac_pipeline.py
import argparse
def main():
    ap=argparse.ArgumentParser()
    ap.add_argument("--home", required=True)  # "lat,lon"
    ap.add_argument("--quake", required=True) # "lat,lon"
    ap.add_argument("--radius", type=int, default=12000)
    ap.add_argument("--notify", default="line,voice")
    args=ap.parse_args()
    # 1) データ更新→ 2) ルート再計算→ 3) 到達時間→ 4) 安全エリア→ 5) 保存→ 6) 通知
    # （既存コードを関数化して順番に呼ぶ）
if __name__=="__main__":
    main()

—

# 4) GraphDB(Neo4j) スキーマ初期化

**目的**: 「質問に答えられる」グラフへ。
**DoD**: `Schema v1` が作成され、`Route(latest=true)` が常に1本だけ存在。

**Cypher（最小）**

CREATE CONSTRAINT shelter_id IF NOT EXISTS FOR (s:Shelter) REQUIRE s.sid IS UNIQUE;
CREATE CONSTRAINT route_id   IF NOT EXISTS FOR (r:Route)   REQUIRE r.rid IS UNIQUE;

UNWIND $shelters AS s
MERGE (sh:Shelter {sid:s.sid})
  SET sh.name=s.name, sh.lat=s.lat, sh.lon=s.lon, sh.safe=s.safe;

MERGE (h:Home {hid:'home-1'}) SET h.lat=$home.lat, h.lon=$home.lon;

WITH $route AS r
MERGE (rt:Route {rid:r.rid})
  SET rt.meters=r.meters, rt.eta_walk=r.eta_walk, rt.eta_walk_safe=r.eta_walk_safe,
      rt.arrival_time_min=r.arrival_time_min, rt.height_class=r.height_class,
      rt.created_at=r.created_at, rt.latest=true;

MATCH (h:Home {hid:'home-1'}), (best:Shelter {sid:$best_sid})
MERGE (h)-[:SAFE_PATH_TO {mode:'walk'}]->(best);

MATCH (x:Route {latest:true}) WHERE x.rid <> $route.rid
SET x.latest=false;

—

# 5) 緊急ダッシュボード（Folium MVP）

**目的**: 家族が1枚で把握。
**DoD**: `outputs/evac.html` に以下が出る。

* `route_latest`（距離/徒歩ETA/安全係数/推定到達までX分）
* `safe_area_now`（半透明塗り）
* 「Googleマップで開く」直リンク

**追加断片**

folium.GeoJson(safe_area_now.__geo_interface__, name="安全エリア", style_function=lambda x: {
    "fillOpacity":0.25, "weight":0
}).add_to(m)

gmap_url=f"https://www.google.com/maps/dir/?api=1&destination={dest_pt.y},{dest_pt.x}&travelmode=walking"
folium.Marker([dest_pt.y, dest_pt.x],
    popup=f"<a href='{gmap_url}' target='_blank'>Googleマップで開く</a>",
    icon=folium.Icon(color='blue')).add_to(m)

—

# 6) 期限ウォッチャ（非常食）

**目的**: 「切れてた」を無くす。
**DoD**: SQLite(またはNeo4j)から**14日以内期限**を抽出し、朝9時にLINE/音声通知。

**SQL 例**

SELECT name, qty, expiry_date
FROM supplies
WHERE DATE(expiry_date) <= DATE('now','+14 day')
ORDER BY expiry_date ASC;

—

# 7) イベントトリガー（擬似EEW→本番）

**目的**: 自動化の完成。
**DoD**: ダミーの震源JSONを受けて③のCLIを起動→将来EEW/津波情報に差し替え。

**擬似Webhook**

curl -X POST http://localhost:8080/quake \
  -H "Content-Type: application/json" \
  -d '{"lat":33.0,"lon":138.5,"mag":7.6,"depth_km":20}'

—

## 実装順（2スプリント想定）

1. **③ CLI化** → 既存処理の関数化・引数化
2. **1) 到達時間v0** → `route_latest` 属性追加
3. **2) 安全エリア** → 反転×等時間到達圏を保存＆Folium表示
4. **5) ダッシュボード** → GMapリンク・要点表示
5. **4) GraphDB投入** → Route/Shelter/関係だけ先行
6. **6) 期限ウォッチャ** → 朝9時通知
7. **7) 擬似EEWトリガー** → 自動化通し

必要なら、上記の**等時間到達圏（isochrone）生成関数**と**Neo4j投入用の具体的なPythonコード**をすぐ出します。どこから書くか一言くれたら、その部分を一気に仕上げます。

ということでまずは
上記の等時間到達圏（isochrone）生成関数から作成

次に作る機能について

これで

使用する用途のコンセプトまた震源地から津波が到達するまでの時間と津波の高さを計算し GIS情報を使い避難するエリアを算出しルート案内できるようにしたい家族の避難経路・地域防災マップをGraphDBで整理避難所や重要連絡先も紐付け使い方例: 「地震が来たら最寄りの避難所はどこ？」「非常食の賞味期限が近いものは？」などこの用途で使う
また
防災・非常時情報管理
やること:
非常持ち出し品・家族の避難経路・地域防災マップをGraphDBで整理
避難所や重要連絡先も紐付け
使い方例:
「地震が来たら最寄りの避難所はどこ？」
「非常食の賞味期限が近いものは？」
メリット:
緊急時に必要な情報を瞬時に取得
をやりたい

またYoutubeライブカメラなどで現地の情報をgeminiなどでマルチモーダルを使うことで解説できるはず。
また震源地から津波が到達するまでの時間と津波の高さを計算し
QGIS情報を使い避難するエリアを算出しルート案内できるようにしたい

地震速報などが出たときや揺れを検知したタイミングで自動動作させたい

この場合、次に作る機能は？

最高です。そのコンセプトなら、まずは「自動で動く骨格」と「意思決定に直結する出力」を先に固めましょう。次に作る機能を、優先度つきの“最短実用ルート”で並べます。すぐ作れる粒度で書きました。

# ① 地震トリガー → 自動実行（最優先）

**目的**: 揺れ/速報を検知したら、対話なしで「最新の安全経路」と「要点ダイジェスト」を家族に配信。
**やること**

* トリガー入力（どれか1つでOK、複数あれば冗長化）

* EEW（緊急地震速報）/ 津波情報のフィード監視（RSS/JSON/PUSH）
* ローカル加速度センサ（Raspberry Pi + MEMS）でしきい値超え検知
* **自動パイプライン起動**：

1. 震源・規模を取得 → 2) 津波到達目安・高さの推定 → 3) 浸水想定×道路を再マスク → 4) `route_latest` を再計算 → 5) 配信
* **LINE/音声配信**（あなたの既存基盤を活用）

* VoiceVox で読み上げ（「徒歩X分、安全側Y分。最寄りは◯◯小学校」）
* LINE Bot へ地図URL（`outputs/evac.html`）＋テキスト要約

> 実装ヒント：`quake_watcher.py`（常駐）→ イベント時に `evac_pipeline.py –mode=auto –home ` を起動

—

# ② 津波「到達時間・高さ」簡易推定モジュール

**目的**: 「あと何分で来る？」に即答。初期は簡易/安全側、後で精緻化。
**ステップ**

* v0（当日実用版）

* `arrival_time_min`, `height_class`, `confidence` を `route_latest`/`sim_runs` に書込

—

# ③ 「安全エリア」ポリゴン生成（QGIS/GeoPandas）

**目的**: 地図上で「今向かうべきエリア」を一目で示す。
**やること**

* 浸水想定（`tsunami_inundation`）の**反転**＋標高/堤防レイヤ（任意）で**安全領域**を生成
* **到達時間で内側制限**（例：徒歩Y分の等時間到達圏 isochrone）× 安全領域 = 「到達可能×安全」領域
* GPKGレイヤ例：`safe_area_now`（上書き）、`isochrone_walk_Ymin`（上書き）

—

# ④ ルート再計算の自動化（既存コードの自動再実行）

**目的**: クリックなしで最新ルートが常に用意される。
**やること**

* 既に作成済みの `roads_safe` / `shelters_safe` / `route_latest` 更新処理を

* **パラメトリックCLI化**：`evac_pipeline.py –home 34.728,137.978 –radius 12000`
* **一時ファイル→置換**でアトミックに上書き（停電等でも壊れない）
* **キャッシュキー**（あなたの `rid` 実装）を、**震源メタ**・**津波メタ**も内包して無効化制御

—

# ⑤ 家族向け「緊急ダッシュボード」（最短 MVP）

**目的**: 家族が迷わず見る1枚。
**やること**

* `outputs/evac.html` に以下を追記（Folium簡易UIでOK）

* 現在の `route_latest`（距離/ETA/安全係数）
* 「推定到達まで X 分 / 高さ区分：◯◯（安全側）」
* 「最寄り安全避難所：◯◯」「次点：◯◯」
* **赤ボタン：『ナビ開始（Googleマップ）』**（目的地座標埋め込みURL）
* LINEにも同じ情報をテンプレで送信（短文＋URL）

—

# ⑥ GraphDB スキーマ拡張（避難・備蓄）

**目的**: 「何を、どこに、誰が、いつまで」をグラフで統合。
**ノード**

* `Person`（家族）/ `Shelter` / `Place(POI)` / `Contact`（電話・LINE）
* `HazardZone` / `Route`（`latest`/`run_id`）/ `Supply`（非常食・水・電池）
**リレーション**
* `LIVES_AT`, `NEAR`, `SAFE_PATH_TO`, `INSIDE`, `HAS_CONTACT`, `HAS_SUPPLY`, `BEST_OPTION`
**最低限のプロパティ**
* `Route`: meters, eta_min_walk, eta_min_walk_safe, created_at, arrival_time_min, height_class
* `Supply`: name, qty, expiry_date, storage_place
**使い方例**
* 「最寄りの安全避難所」→ `MATCH (home)-[:NEAR]->(s: Shelter) WHERE s.safe=true RETURN s LIMIT 1`
* 「賞味期限が近い」→ `MATCH (sup:Supply) WHERE sup.expiry_date < date()+duration('P14D') RETURN sup` --- # ⑦ 非常食・備蓄の賞味期限ウォッチャ **目的**: 「切らしてた/期限切れ」をなくす。 **やること** * `supplies.db`（SQLite） or GraphDB に `Supply` 登録 * 毎朝 9時に期限14日以内を抽出 → LINE/音声通知 * 地震トリガー時は**即時**に「持ち出し優先リスト」を配信（例：水/ヘッドライト/救急） --- # ⑧ ライブカメラ要約（マルチモーダル） **目的**: 「今、現地はどう見える？」を数十秒で要約。 **やること** * YouTube Live の**固定リスト**（地域局・河川・港湾）を管理 * イベント時：最新サムネイル or 数秒キャプチャを**LLMに要約**（水位・人流・通行可否など） * 要約に「**不確実性/注意喚起**」を必ず付け、安全側誘導 --- # ⑨ 訓練モード & 失敗時フォールバック **目的**: いざ本番でも動く。 **やること** * `--drill` フラグで**擬似震源**を流し、全自動ルート→配信まで通す * 外部データ失敗時は**ローカルのみ**で保守的判定（津波到達は短め、浸水は広め） --- # ⑩ すぐに書き足すスクリプト断片（抜粋） **自動起動フック（例・擬似）** [shell] # systemd 例 [Unit] Description=Quake watcher After=network-online.target [Service] ExecStart=/usr/bin/python3 /opt/evac/quake_watcher.py Restart=always [Install] WantedBy=multi-user.target ``` **パイプライン入口** ```bash python evac_pipeline.py \ --home 34.728,137.978 \ --radius 12000 \ --mode auto \ --notify line,voice ``` [/shell] ## 今日からの実装順（最短2–3スプリント想定） 1. **①トリガー→自動実行**（最低限：ローカル加速度しきい値 or ダミーEEW） 2. **④ルート再計算のCLI化**（既存コードを関数化→引数で動く） 3. **⑤ダッシュボードMVP**（Foliumに要点3行＋ナビボタン） 4. **②到達時間v0**（一定速度・安全側近似）＋`route_latest`に追記 5. **⑦賞味期限ウォッチャ**（朝9時通知） 6. **⑥GraphDBスキーマ**（`Route`と`Supply`だけ先に） 7. **⑧ライブ要約**（固定ライブ2–3本から） --- 必要なら、この骨格に合わせて**`evac_pipeline.py` の雛形**（CLI引数・モジュール分割）と **Neo4jの初期Cypher** もすぐ書きます。どこから着手するかだけ言ってくれたら、実装用のコードをそのまま渡します。がGPTの答え

コンセプトに基づく改良

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString
import osmnx as ox
import networkx as nx
import folium
import os

GPKG = "outputs/shizuoka_hazard.gpkg"

# 1) 既存レイヤ読み込み（両方とも 4326）
coast   = gpd.read_file(GPKG, layer="coastline")
tsunami = gpd.read_file(GPKG, layer="tsunami_inundation")

# ----- パラメータ -----
# 自宅座標（lat, lon）←必ず置き換え
home_wgs = (34.728, 137.978)
# OSM道路の取得半径[m]（到達不能なら大きく：15000, 20000…）
GRAPH_RADIUS_M = 12000
# 避難所CSV（cp932想定）
SHELTER_CSV = "data/20250801_指定緊急避難場所一覧.CSV"
# ----------------------

# 2) 避難所CSV→点化（列名は自動検出）
df = pd.read_csv(SHELTER_CSV, encoding="cp932")
lat_candidates = 
lon_candidates = 
if not lat_candidates or not lon_candidates:
    raise ValueError("CSVに緯度/経度列が見つかりません。列名を教えてください。")
lat_col, lon_col = lat_candidates[0], lon_candidates[0]
df[lat_col] = pd.to_numeric(df[lat_col], errors="coerce")
df[lon_col] = pd.to_numeric(df[lon_col], errors="coerce")
df = df.dropna(subset=[lat_col, lon_col])
shelters = gpd.GeoDataFrame(df.copy(),
    geometry=gpd.points_from_xy(df[lon_col], df[lat_col], crs="EPSG:4326")
)

# 3) 安全な避難所 = 浸水域と非交差
#   （unionはメモリ負荷になるので、空間結合で判定）
safe_join = gpd.sjoin(shelters, tsunami[["geometry"]], how="left", predicate="intersects")
shelters["safe"] = safe_join["index_right"].isna()
safe_shelters = shelters[shelters["safe"]].copy()

# --- 4) OSM道路取得（徒歩網） ---
G = ox.graph_from_point(home_wgs, dist=GRAPH_RADIUS_M, network_type="walk", simplify=True, retain_all=False)
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True)

# u,v,key を“列”にする（重要）
edges = edges.reset_index()

# CRSを4326に
if edges.crs is None:
    edges = edges.set_crs(4326)
elif edges.crs.to_epsg() != 4326:
    edges = edges.to_crs(4326)
# --- 5) 浸水域と交差する道路を除去（sjoin→集計で一意化） ---
left = edges[["u","v","key","geometry"]].copy()   # edges は reset_index() 済み想定
right = tsunami[["geometry"]]

tag = gpd.sjoin(left, right, how="left", predicate="intersects")

blocked = (
    tag.assign(hit=tag["index_right"].notna())
       .groupby(["u","v","key"], as_index=False)["hit"].max()
       .rename(columns={"hit":"blocked"})
)

edges2 = edges.merge(blocked, on=["u","v","key"], how="left")
edges2["blocked"] = edges2["blocked"].fillna(False)

edges_keep = edges2.loc[~edges2["blocked"]].copy()

# ★ ここが重要：graph_from_gdfs に渡す前に MultiIndex に戻す
#   （key が float になっていたら int にしておくと安心）
if edges_keep["key"].isna().any():
    edges_keep["key"] = edges_keep["key"].fillna(0)

edges_keep["key"] = edges_keep["key"].astype(int)
edges_keep = edges_keep.set_index(["u","v","key"]).sort_index()

# --- 6) サブグラフ再構築（通行可のみ） ---
used_nodes = pd.unique(edges_keep.reset_index()[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()

# nodes_keep の geometry 列が欠けていたら補完
if "geometry" not in nodes_keep or nodes_keep["geometry"].isna().any():
    nodes_keep = nodes_keep.set_geometry(
        gpd.points_from_xy(nodes_keep["x"], nodes_keep["y"]), crs=edges_keep.crs
    )

# OSMnx のバージョン差に対応
try:
    G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))
except AttributeError:
    from osmnx import utils_graph as ug
    G2 = ug.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))


# 7) 出発・目的地ノード（最寄り“安全”避難所）
home_pt = Point(home_wgs[1], home_wgs[0])
# 最寄り安全避難所（直線距離）…なければ全体から選ぶ
TARGET = safe_shelters if len(safe_shelters) else shelters
METRIC = "EPSG:32654"# 静岡周辺 UTM
hp = gpd.GeoSeries([home_pt], crs=4326).to_crs(METRIC).iloc[0]
tp = TARGET.to_crs(METRIC).copy()
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = TARGET.loc[tp["__dist_m__"].idxmin(), "geometry"]

# ▶ 互換呼び分け
try:
    orig = ox.distance.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)
except AttributeError:
    # 新API名（環境によってはこちらが有効）
    orig = ox.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)


# 8) 最短経路（通行可能グラフ）
try:
    route = nx.shortest_path(G2, orig, dest, weight="length")
except nx.NetworkXNoPath:
    raise RuntimeError("到達経路なし。GRAPH_RADIUS_Mを増やすか、別の避難所で再試行を。")

# 9) GPKGへ追記
#    既にある場合は上書きしたいなら driverオプションを切り替える（ここは追記）
shelters.to_file(GPKG, layer="shelters_all", driver="GPKG")
safe_shelters.to_file(GPKG, layer="shelters_safe", driver="GPKG")
edges.to_file(GPKG, layer="roads_all", driver="GPKG")
edges_keep.to_file(GPKG, layer="roads_safe", driver="GPKG")

# 10) Folium地図
m = folium.Map(location=[home_wgs[0], home_wgs[1]], zoom_start=12, control_scale=True)
folium.GeoJson(tsunami.__geo_interface__, name="津波浸水想定(A40)").add_to(m)
folium.GeoJson(edges_keep.__geo_interface__, name="道路(浸水回避)").add_to(m)

# 経路を線で描画
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(route_latlon, weight=6, opacity=0.9, tooltip="避難経路").add_to(m)

# 出発・目的地
folium.Marker([home_wgs[0], home_wgs[1]], popup="現在地", icon=folium.Icon(color="green")).add_to(m)
folium.Marker([dest_pt.y, dest_pt.x], popup="避難所", icon=folium.Icon(color="blue")).add_to(m)

# 避難所レイヤ（安全のみ強調）
folium.GeoJson(safe_shelters.__geo_interface__, name="安全な避難所",
               marker=folium.CircleMarker(radius=4)).add_to(m)
folium.LayerControl().add_to(m)

# 保存（任意）
m.save("outputs/evac.html")
print("✓ 完了: GPKG追記 & outputs/evac.html を作成")

の結果

---------------------------------------------------------------------------
CPLE_AppDefinedError                      Traceback (most recent call last)
File pyogrio/_io.pyx:2344, in pyogrio._io.create_ogr_dataset_layer()

File pyogrio/_err.pyx:218, in pyogrio._err.check_pointer()

CPLE_AppDefinedError: Layer roads_all already exists, CreateLayer failed. Use the layer creation option OVERWRITE=YES to replace it.

During handling of the above exception, another exception occurred:

DataLayerError                            Traceback (most recent call last)
Cell In[2], line 129
    127 shelters.to_file(GPKG, layer="shelters_all", driver="GPKG")
    128 safe_shelters.to_file(GPKG, layer="shelters_safe", driver="GPKG")
--> 129 edges.to_file(GPKG, layer="roads_all", driver="GPKG")
    130 edges_keep.to_file(GPKG, layer="roads_safe", driver="GPKG")
    132 # 10) Folium地図

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/geodataframe.py:1637, in GeoDataFrame.to_file(self, filename, driver, schema, index, **kwargs)


   1543 """Write the ``GeoDataFrame`` to a file.
   1544 
   1545 By default, an ESRI shapefile is written, but any OGR data source
   (...)   1633 
   1634 """
   1635 from geopandas.io.file import _to_file
-> 1637 _to_file(self, filename, driver, schema, index, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/io/file.py:731, in _to_file(df, filename, driver, schema, index, mode, crs, engine, metadata, **kwargs)
    728     raise ValueError(f"'mode' should be one of 'w' or 'a', got '{mode}' instead")
    730 if engine == "pyogrio":
--> 731     _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs)
    732 elif engine == "fiona":
    733     _to_file_fiona(df, filename, driver, schema, crs, mode, metadata, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/io/file.py:793, in _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs)
    790 if not df.columns.is_unique:
    791     raise ValueError("GeoDataFrame cannot contain duplicated column names.")
--> 793 pyogrio.write_dataframe(df, filename, driver=driver, metadata=metadata, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/geopandas.py:710, in write_dataframe(df, path, layer, driver, encoding, geometry_type, promote_to_multi, nan_as_null, append, use_arrow, dataset_metadata, layer_metadata, metadata, dataset_options, layer_options, **kwargs)
    707         field_data.append(values)
    708         field_mask.append(None)
--> 710 write(
    711     path,
    712     layer=layer,
    713     driver=driver,
    714     geometry=geometry,
    715     field_data=field_data,
    716     field_mask=field_mask,
    717     fields=fields,
    718     crs=crs,
    719     geometry_type=geometry_type,
    720     encoding=encoding,
    721     promote_to_multi=promote_to_multi,
    722     nan_as_null=nan_as_null,
    723     append=append,
    724     dataset_metadata=dataset_metadata,
    725     layer_metadata=layer_metadata,
    726     metadata=metadata,
    727     dataset_options=dataset_options,
    728     layer_options=layer_options,
    729     gdal_tz_offsets=gdal_tz_offsets,
    730     **kwargs,
    731 )

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723, in write(path, geometry, field_data, fields, field_mask, layer, driver, geometry_type, crs, encoding, promote_to_multi, nan_as_null, append, dataset_metadata, layer_metadata, metadata, dataset_options, layer_options, gdal_tz_offsets, **kwargs)
    718 # preprocess kwargs and split in dataset and layer creation options
    719 dataset_kwargs, layer_kwargs = _preprocess_options_kwargs(
    720     driver, dataset_options, layer_options, kwargs
    721 )
--> 723 ogr_write(
    724     path,
    725     layer=layer,
    726     driver=driver,
    727     geometry=geometry,
    728     geometry_type=geometry_type,
    729     field_data=field_data,
    730     field_mask=field_mask,
    731     fields=fields,
    732     crs=crs,
    733     encoding=encoding,
    734     promote_to_multi=promote_to_multi,
    735     nan_as_null=nan_as_null,
    736     append=append,
    737     dataset_metadata=dataset_metadata,
    738     layer_metadata=layer_metadata,
    739     dataset_kwargs=dataset_kwargs,
    740     layer_kwargs=layer_kwargs,
    741     gdal_tz_offsets=gdal_tz_offsets,
    742 )

File pyogrio/_io.pyx:2463, in pyogrio._io.ogr_write()

File pyogrio/_io.pyx:2360, in pyogrio._io.create_ogr_dataset_layer()

DataLayerError: Layer roads_all already exists, CreateLayer failed. Use the layer creation option OVERWRITE=YES to replace it.

となった

原因は
エラーの原因は、そのレイヤー名（roads_all）が GPKG 内にすでに存在しているのに、上書きオプションを付けずに to_file() しているためです。対処は「上書き」するか「追記」するかの2択です。再実行に強いのは上書きです。

いくつか方法を提示されたが、コンセプトに基づき改造しないと後で全部修正になる

使用する用途のコンセプトまた震源地から津波が到達するまでの時間と津波の高さを計算し GIS情報を使い避難するエリアを算出しルート案内できるようにしたい家族の避難経路・地域防災マップをGraphDBで整理避難所や重要連絡先も紐付け使い方例: 「地震が来たら最寄りの避難所はどこ？」「非常食の賞味期限が近いものは？」などこの用途で使うなら今回の場合は（上書き）既存レイヤーを削除してから書く同じレイヤーにレコードを追加どれを推奨するか？

結論から：
* 今回の用途（津波到達時間・高さを計算 → 安全エリア算出 → その時点で最良の避難経路を提示）では、 “計算で毎回作り直す派生レイヤー”は原則《上書き》を推奨です。
* 履歴を残したいテーブル（実行ログ・過去の経路・賞味期限チェック結果など）は《追記》を併用します。
* 《既存レイヤー削除→書き込み》は、スキーマを頻繁に変えるときやエンジン差の不具合回避用のセーフティとして使う位置づけ。

なぜ《上書き》が最適か
あなたのワークフローは「同じ入力（ベースの津波想定・OSM・避難所）→同じ処理→決定論的な派生成果物（roads_safe, shelters_safe, route_latest など）」という再現可能・冪等パイプラインです。 このタイプは毎回クリーンに作り直せることが最重要なので、レイヤー名固定＋上書きが最もシンプルで壊れにくいです。
* 冪等：再実行で結果が混ざらない
* 運用：下流（Folium表示／GraphDB同期）が常に同じレイヤー名を参照
* 事故防止：前回分の残骸や列増減の不一致で壊れない

まとめ（今回どれを勧める？）
* メインは《上書き》一択（派生レイヤーは常に最新をクリーンに再生成）
* 履歴が必要なものだけ《追記》を別レイヤーで併用
* スキーマ変更直後や不具合回避には《削除→書き込み》をスポット的に
この方針にすると、 「地震が来たら最寄りの避難所は？」→ shelters_safe＋route_latest を即参照、 「非常食の賞味期限が近いものは？」→ food_expiry_alerts（追記ログ）や在庫DB を参照、 という形で運用がシンプル＆堅牢になります。

# 9) GPKGへ追記
# #    既にある場合は上書きしたいなら driverオプションを切り替える（ここは追記）
# shelters.to_file(GPKG, layer="shelters_all", driver="GPKG")
# safe_shelters.to_file(GPKG, layer="shelters_safe", driver="GPKG")
# edges.to_file(GPKG, layer="roads_all", driver="GPKG")
# edges_keep.to_file(GPKG, layer="roads_safe", driver="GPKG")

# 9) GPKGへ書き出し（毎回クリーンに上書き）
common_kwargs = dict(
    driver="GPKG",
    engine="pyogrio",
    layer_options={"OVERWRITE": "YES"},  # 既存レイヤーを置き換える
)

shelters.to_file(GPKG, layer="shelters_all",   **common_kwargs)
safe_shelters.to_file(GPKG, layer="shelters_safe", **common_kwargs)
edges.to_file(GPKG,    layer="roads_all",     **common_kwargs)
edges_keep.to_file(GPKG, layer="roads_safe",  **common_kwargs)

というようにコードを変更

これでエラーは解決

しかし

# # === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
# from shapely.geometry import LineString
# import geopandas as gpd

# # 1) ルートのLineString（まだ属性は付けない）
# route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
# route_geom = LineString(route_lonlat)
# route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# # 2) 堅牢チェック（交差ゼロ＆境界近接なし）
# #   union_all が使えない環境向けにフォールバックも用意
# try:
#     flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
# except Exception:
#     flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

# METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
# route_m = route_gdf.to_crs(METRIC)
# flood_m = flood_union.to_crs(METRIC)

# # 交差長[m]
# inter = gpd.overlay(route_m, flood_m, how="intersection")
# inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# # 近接（境界 ±2m 以内をNGにする例）
# # near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
# buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
# near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty


# if inside_len_m > 0 or near:
#     raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
#                        "別避難所 or 半径拡大で再探索してください。")

# # 3) 距離[m]→ETA（バージョン互換）
# try:
#     route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
#     meters = float(route_edges["length"].sum())
# except Exception:
#     meters = 0.0
#     for u, v in zip(route[:-1], route[1:]):
#         ed = G2.get_edge_data(u, v)
#         meters += min(d.get("length", 0.0) for d in ed.values())

# walk_kmh = 4.8
# drive_kmh = 20.0
# eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
# eta_drive_min = meters / (drive_kmh * 1000 / 60)
# eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# # 4) GPKG に“合格ルート”だけ属性付きで保存
# out_gpkg = "outputs/shizuoka_hazard.gpkg"
# route_gdf = gpd.GeoDataFrame(
#     {
#         "meters":[meters],
#         "eta_min_walk":[eta_walk_min],
#         "eta_min_walk_safe":[eta_walk_safe],
#         "eta_min_drive":[eta_drive_min],
#     },
#     geometry=[route_geom], crs="EPSG:4326",
# )
# route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
# print("✅ 交差なし・保存完了")

# # 5) （任意）Foliumに表示するならここで
# route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
# folium.PolyLine(
#     route_latlon, weight=6, opacity=0.9,
#     tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
# ).add_to(m)

# === ルート検証 → ETA → 保存（キャッシュ対応） ===
from shapely.geometry import LineString, Point
import geopandas as gpd
import hashlib, json
from datetime import datetime
from fiona import listlayers
import osmnx as ox
import networkx as nx

GPKG = "outputs/shizuoka_hazard.gpkg"
METRIC = "EPSG:32654"  # 静岡周辺

# 0) 津波 union はキャッシュを優先（なければ作る）
if "tsunami_union" in listlayers(GPKG):
    flood_union = gpd.read_file(GPKG, layer="tsunami_union")
else:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs=tsunami.crs)
    flood_union.to_file(GPKG, layer="tsunami_union", driver="GPKG")

# 0') ハザードの簡易メタ（A40が変わったらルートキャッシュ無効化に使う）
ti = gpd.read_file(GPKG, layer="tsunami_inundation")
hazard_meta = f'{len(ti)}|{",".join(map(str, ti.total_bounds))}'

# 1) キャッシュキー生成（端点・モード・半径・ハザードメタ）
def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000, hazard_meta=None):
    payload = {"home": home_wgs, "dest": dest_wgs, "mode": mode, "r": radius_m, "hz": hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk",
                radius_m=GRAPH_RADIUS_M, hazard_meta=hazard_meta)

# 2) ルートキャッシュにヒットすれば使う
route_gdf = None
use_cached_route = False
if "route_cache" in listlayers(GPKG):
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"] == key]
    if len(hit):
        route_gdf = hit.set_crs(4326)
        meters = float(route_gdf["meters"].iloc[0])
        eta_walk_min  = float(route_gdf["eta_min_walk"].iloc[0])
        eta_walk_safe = float(route_gdf["eta_min_walk_safe"].iloc[0])
        eta_drive_min = float(route_gdf["eta_min_drive"].iloc[0])
        coords = list(route_gdf.geometry.iloc[0].coords)
        route_latlon = [(lat, lon) for lon, lat in coords]  # folium用 (lat,lon)
        use_cached_route = True
        print("↩︎ ルート: キャッシュヒット")

if not use_cached_route:
    # 3) これまで通りの経路（'route'）から LineString を作成
    route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
    route_geom = LineString(route_lonlat)
    route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

    # 4) 堅牢チェック（交差ゼロ＆境界近接なし）
    route_m = route_gdf.to_crs(METRIC)
    flood_m = flood_union.to_crs(METRIC)

    inter = gpd.overlay(route_m, flood_m, how="intersection")
    inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

    buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
    near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

    if inside_len_m > 0 or near:
        raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                           "別避難所 or 半径拡大で再探索してください。")

    # 5) 距離[m]→ETA
    try:
        route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
        meters = float(route_edges["length"].sum())
    except Exception:
        meters = 0.0
        for u, v in zip(route[:-1], route[1:]):
            ed = G2.get_edge_data(u, v)
            meters += min(d.get("length", 0.0) for d in ed.values())

    walk_kmh = 4.8
    drive_kmh = 20.0
    eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
    eta_drive_min = meters / (drive_kmh * 1000 / 60)
    eta_walk_safe = eta_walk_min * 1.3

    # 6) route_latest へ保存（検証済みのみ）
    out_gpkg = GPKG
    route_gdf = gpd.GeoDataFrame(
        {
            "meters":[meters],
            "eta_min_walk":[eta_walk_min],
            "eta_min_walk_safe":[eta_walk_safe],
            "eta_min_drive":[eta_drive_min],
            "created_at":[datetime.utcnow().isoformat(timespec="seconds")+"Z"],
        },
        geometry=[route_geom], crs="EPSG:4326",
    )
    route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
    print("✅ 交差なし・保存完了")

    # 7) ルートキャッシュに追記
    route_cache_rec = route_gdf.copy()
    route_cache_rec["rid"] = key
    try:
        route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG", if_exists="append")
    except TypeError:
        # 古いGeoPandasのためのフォールバック
        if "route_cache" in listlayers(out_gpkg):
            old = gpd.read_file(out_gpkg, layer="route_cache")
            gpd.GeoDataFrame(pd.concat([old, route_cache_rec], ignore_index=True), crs=4326)\
               .to_file(out_gpkg, layer="route_cache", driver="GPKG")
        else:
            route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG")

    # 8) folium用座標
    route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]

# 9) Folium表示（キャッシュ・新規どちらでも）
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

で

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 175
    173 route_cache_rec["rid"] = key
    174 try:
--> 175     route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG", if_exists="append")
    176 except TypeError:
    177     # 古いGeoPandasのためのフォールバック
    178     if "route_cache" in listlayers(out_gpkg):

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/geodataframe.py:1637, in GeoDataFrame.to_file(self, filename, driver, schema, index, **kwargs)
   1543 """Write the ``GeoDataFrame`` to a file.
   1544 
   1545 By default, an ESRI shapefile is written, but any OGR data source
   (...)   1633 
   1634 """
   1635 from geopandas.io.file import _to_file
-> 1637 _to_file(self, filename, driver, schema, index, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/io/file.py:731, in _to_file(df, filename, driver, schema, index, mode, crs, engine, metadata, **kwargs)
    728     raise ValueError(f"'mode' should be one of 'w' or 'a', got '{mode}' instead")
    730 if engine == "pyogrio":
--> 731     _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs)
    732 elif engine == "fiona":
    733     _to_file_fiona(df, filename, driver, schema, crs, mode, metadata, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/io/file.py:793, in _to_file_pyogrio(df, filename, driver, schema, crs, mode, metadata, **kwargs)
    790 if not df.columns.is_unique:
    791     raise ValueError("GeoDataFrame cannot contain duplicated column names.")
--> 793 pyogrio.write_dataframe(df, filename, driver=driver, metadata=metadata, **kwargs)

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/geopandas.py:710, in write_dataframe(df, path, layer, driver, encoding, geometry_type, promote_to_multi, nan_as_null, append, use_arrow, dataset_metadata, layer_metadata, metadata, dataset_options, layer_options, **kwargs)
    707         field_data.append(values)
    708         field_mask.append(None)
--> 710 write(
    711     path,
    712     layer=layer,
    713     driver=driver,
    714     geometry=geometry,
    715     field_data=field_data,
    716     field_mask=field_mask,
    717     fields=fields,
    718     crs=crs,
    719     geometry_type=geometry_type,
    720     encoding=encoding,
    721     promote_to_multi=promote_to_multi,
    722     nan_as_null=nan_as_null,
    723     append=append,
    724     dataset_metadata=dataset_metadata,
    725     layer_metadata=layer_metadata,
    726     metadata=metadata,
    727     dataset_options=dataset_options,
    728     layer_options=layer_options,
    729     gdal_tz_offsets=gdal_tz_offsets,
    730     **kwargs,
    731 )

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:719, in write(path, geometry, field_data, fields, field_mask, layer, driver, geometry_type, crs, encoding, promote_to_multi, nan_as_null, append, dataset_metadata, layer_metadata, metadata, dataset_options, layer_options, gdal_tz_offsets, **kwargs)
    711     warnings.warn(
    712         "'crs' was not provided.  The output dataset will not have "
    713         "projection information defined and may not be usable in other "
    714         "systems.",
    715         stacklevel=2,
    716     )
    718 # preprocess kwargs and split in dataset and layer creation options
--> 719 dataset_kwargs, layer_kwargs = _preprocess_options_kwargs(
    720     driver, dataset_options, layer_options, kwargs
    721 )
    723 ogr_write(
    724     path,
    725     layer=layer,
   (...)    741     gdal_tz_offsets=gdal_tz_offsets,
    742 )

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:514, in _preprocess_options_kwargs(driver, dataset_options, layer_options, kwargs)
    512             layer_kwargs[k] = v
    513         else:
--> 514             raise ValueError(f"unrecognized option '{k}' for driver '{driver}'")
    516 return dataset_kwargs, layer_kwargs

ValueError: unrecognized option 'IF_EXISTS' for driver 'GPKG'

原因は
if_exists=”append” は pyogrio では非対応なのでエラーになっています。
to_file() は engine=”pyogrio” + append=True を使ってください。ついでに route_latest など「毎回作り直す」レイヤーは OVERWRITE で上書きするのが安全です

要点
* pyogrio では if_exists=”append” は使えません → append=True を使う。
* 毎回再計算するレイヤーは layer_options={“OVERWRITE”:”YES”} で確実に上書き。
* 追記でスキーマ差が出たら、読み出し→列合わせ→連結→上書きで堅牢化。

route_latest は毎回クリーンに上書き、
route_cache はappendで追記（pyogrio準拠）、
tsunami_union も必要に応じて上書き耐性で作成、 の3点が組み込んでコードを書き換える

# # === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
# from shapely.geometry import LineString
# import geopandas as gpd

# # 1) ルートのLineString（まだ属性は付けない）
# route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
# route_geom = LineString(route_lonlat)
# route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# # 2) 堅牢チェック（交差ゼロ＆境界近接なし）
# #   union_all が使えない環境向けにフォールバックも用意
# try:
#     flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
# except Exception:
#     flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

# METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
# route_m = route_gdf.to_crs(METRIC)
# flood_m = flood_union.to_crs(METRIC)

# # 交差長[m]
# inter = gpd.overlay(route_m, flood_m, how="intersection")
# inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# # 近接（境界 ±2m 以内をNGにする例）
# # near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
# buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
# near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty


# if inside_len_m > 0 or near:
#     raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
#                        "別避難所 or 半径拡大で再探索してください。")

# # 3) 距離[m]→ETA（バージョン互換）
# try:
#     route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
#     meters = float(route_edges["length"].sum())
# except Exception:
#     meters = 0.0
#     for u, v in zip(route[:-1], route[1:]):
#         ed = G2.get_edge_data(u, v)
#         meters += min(d.get("length", 0.0) for d in ed.values())

# walk_kmh = 4.8
# drive_kmh = 20.0
# eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
# eta_drive_min = meters / (drive_kmh * 1000 / 60)
# eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# # 4) GPKG に“合格ルート”だけ属性付きで保存
# out_gpkg = "outputs/shizuoka_hazard.gpkg"
# route_gdf = gpd.GeoDataFrame(
#     {
#         "meters":[meters],
#         "eta_min_walk":[eta_walk_min],
#         "eta_min_walk_safe":[eta_walk_safe],
#         "eta_min_drive":[eta_drive_min],
#     },
#     geometry=[route_geom], crs="EPSG:4326",
# )
# route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
# print("✅ 交差なし・保存完了")

# # 5) （任意）Foliumに表示するならここで
# route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
# folium.PolyLine(
#     route_latlon, weight=6, opacity=0.9,
#     tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
# ).add_to(m)

# === ルート検証 → ETA → 保存（キャッシュ対応） ===
from shapely.geometry import LineString, Point
import geopandas as gpd
import hashlib, json
from datetime import datetime
from fiona import listlayers
import osmnx as ox
import networkx as nx

GPKG = "outputs/shizuoka_hazard.gpkg"
METRIC = "EPSG:32654"  # 静岡周辺

# 0) 津波 union はキャッシュを優先（なければ作る）
if "tsunami_union" in listlayers(GPKG):
    flood_union = gpd.read_file(GPKG, layer="tsunami_union")
else:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs=tsunami.crs)
    flood_union.to_file(GPKG, layer="tsunami_union", driver="GPKG")

# 0') ハザードの簡易メタ（A40が変わったらルートキャッシュ無効化に使う）
ti = gpd.read_file(GPKG, layer="tsunami_inundation")
hazard_meta = f'{len(ti)}|{",".join(map(str, ti.total_bounds))}'

# 1) キャッシュキー生成（端点・モード・半径・ハザードメタ）
def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000, hazard_meta=None):
    payload = {"home": home_wgs, "dest": dest_wgs, "mode": mode, "r": radius_m, "hz": hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk",
                radius_m=GRAPH_RADIUS_M, hazard_meta=hazard_meta)

# 2) ルートキャッシュにヒットすれば使う
route_gdf = None
use_cached_route = False
if "route_cache" in listlayers(GPKG):
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"] == key]
    if len(hit):
        route_gdf = hit.set_crs(4326)
        meters = float(route_gdf["meters"].iloc[0])
        eta_walk_min  = float(route_gdf["eta_min_walk"].iloc[0])
        eta_walk_safe = float(route_gdf["eta_min_walk_safe"].iloc[0])
        eta_drive_min = float(route_gdf["eta_min_drive"].iloc[0])
        coords = list(route_gdf.geometry.iloc[0].coords)
        route_latlon = [(lat, lon) for lon, lat in coords]  # folium用 (lat,lon)
        use_cached_route = True
        print("↩︎ ルート: キャッシュヒット")

if not use_cached_route:
    # 3) これまで通りの経路（'route'）から LineString を作成
    route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
    route_geom = LineString(route_lonlat)
    route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

    # 4) 堅牢チェック（交差ゼロ＆境界近接なし）
    route_m = route_gdf.to_crs(METRIC)
    flood_m = flood_union.to_crs(METRIC)

    inter = gpd.overlay(route_m, flood_m, how="intersection")
    inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

    buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
    near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

    if inside_len_m > 0 or near:
        raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                           "別避難所 or 半径拡大で再探索してください。")

    # 5) 距離[m]→ETA
    try:
        route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
        meters = float(route_edges["length"].sum())
    except Exception:
        meters = 0.0
        for u, v in zip(route[:-1], route[1:]):
            ed = G2.get_edge_data(u, v)
            meters += min(d.get("length", 0.0) for d in ed.values())

    walk_kmh = 4.8
    drive_kmh = 20.0
    eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
    eta_drive_min = meters / (drive_kmh * 1000 / 60)
    eta_walk_safe = eta_walk_min * 1.3

    # 6) route_latest へ保存（検証済みのみ）
    out_gpkg = GPKG
    route_gdf = gpd.GeoDataFrame(
        {
            "meters":[meters],
            "eta_min_walk":[eta_walk_min],
            "eta_min_walk_safe":[eta_walk_safe],
            "eta_min_drive":[eta_drive_min],
            "created_at":[datetime.utcnow().isoformat(timespec="seconds")+"Z"],
        },
        geometry=[route_geom], crs="EPSG:4326",
    )
    route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
    print("✅ 交差なし・保存完了")

    # 7) ルートキャッシュに追記
    route_cache_rec = route_gdf.copy()
    route_cache_rec["rid"] = key
    try:
        route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG", if_exists="append")
    except TypeError:
        # 古いGeoPandasのためのフォールバック
        if "route_cache" in listlayers(out_gpkg):
            old = gpd.read_file(out_gpkg, layer="route_cache")
            gpd.GeoDataFrame(pd.concat([old, route_cache_rec], ignore_index=True), crs=4326)\
               .to_file(out_gpkg, layer="route_cache", driver="GPKG")
        else:
            route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG")

    # 8) folium用座標
    route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]

# 9) Folium表示（キャッシュ・新規どちらでも）
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

を
# === ルート検証 → ETA → 保存（キャッシュ対応, 上書き＆追記の方針を反映） ===
from shapely.geometry import LineString, Point
import geopandas as gpd
import pandas as pd
import hashlib, json
from datetime import datetime
from fiona import listlayers
import osmnx as ox
import networkx as nx

GPKG = "outputs/shizuoka_hazard.gpkg"
METRIC = "EPSG:32654"  # 静岡周辺

# 0) 津波 union はキャッシュを優先（なければ作る）→ 上書き耐性ありで作成可
if "tsunami_union" in listlayers(GPKG):
    flood_union = gpd.read_file(GPKG, layer="tsunami_union")
else:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs=tsunami.crs)
    flood_union.to_file(
        GPKG, layer="tsunami_union",
        driver="GPKG", engine="pyogrio",
        layer_options={"OVERWRITE": "YES"}  # (3) 上書き耐性
    )

# 0') ハザードの簡易メタ（A40が変わったらルートキャッシュ無効化に使う）
ti = gpd.read_file(GPKG, layer="tsunami_inundation")
hazard_meta = f'{len(ti)}|{",".join(map(str, ti.total_bounds))}'

# 1) キャッシュキー生成（端点・モード・半径・ハザードメタ）
def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000, hazard_meta=None):
    payload = {"home": home_wgs, "dest": dest_wgs, "mode": mode, "r": radius_m, "hz": hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk",
                radius_m=GRAPH_RADIUS_M, hazard_meta=hazard_meta)

# 2) ルートキャッシュにヒットすれば使う
route_gdf = None
use_cached_route = False
if "route_cache" in listlayers(GPKG):
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"] == key]
    if len(hit):
        route_gdf = hit.set_crs(4326)
        meters = float(route_gdf["meters"].iloc[0])
        eta_walk_min  = float(route_gdf["eta_min_walk"].iloc[0])
        eta_walk_safe = float(route_gdf["eta_min_walk_safe"].iloc[0])
        eta_drive_min = float(route_gdf["eta_min_drive"].iloc[0])
        coords = list(route_gdf.geometry.iloc[0].coords)
        route_latlon = [(lat, lon) for lon, lat in coords]  # folium用 (lat,lon)
        use_cached_route = True
        print("↩︎ ルート: キャッシュヒット")

if not use_cached_route:
    # 3) これまで通りの経路（'route'）から LineString を作成
    route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
    route_geom = LineString(route_lonlat)
    route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

    # 4) 堅牢チェック（交差ゼロ＆境界近接なし）
    route_m = route_gdf.to_crs(METRIC)
    flood_m = flood_union.to_crs(METRIC)

    inter = gpd.overlay(route_m, flood_m, how="intersection")
    inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

    buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
    near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

    if inside_len_m > 0 or near:
        raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                           "別避難所 or 半径拡大で再探索してください。")

    # 5) 距離[m]→ETA
    try:
        route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
        meters = float(route_edges["length"].sum())
    except Exception:
        meters = 0.0
        for u, v in zip(route[:-1], route[1:]):
            ed = G2.get_edge_data(u, v)
            meters += min(d.get("length", 0.0) for d in ed.values())

    walk_kmh = 4.8
    drive_kmh = 20.0
    eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
    eta_drive_min = meters / (drive_kmh * 1000 / 60)
    eta_walk_safe = eta_walk_min * 1.3

    # 6) route_latest へ保存（検証済みのみ）→ 上書き
    out_gpkg = GPKG
    route_gdf = gpd.GeoDataFrame(
        {
            "meters":[meters],
            "eta_min_walk":[eta_walk_min],
            "eta_min_walk_safe":[eta_walk_safe],
            "eta_min_drive":[eta_drive_min],
            "created_at":[datetime.utcnow().isoformat(timespec="seconds")+"Z"],
        },
        geometry=[route_geom], crs="EPSG:4326",
    )
    route_gdf.to_file(
        out_gpkg, layer="route_latest",
        driver="GPKG", engine="pyogrio",
        layer_options={"OVERWRITE": "YES"}  # (1) 上書き
    )
    print("✅ 交差なし・保存完了")

    # 7) ルートキャッシュに追記（pyogrio 流儀）
    route_cache_rec = route_gdf.copy()
    route_cache_rec["rid"] = key
    try:
        route_cache_rec.to_file(
            out_gpkg, layer="route_cache",
            driver="GPKG", engine="pyogrio",
            append=True                     # (2) 追記
        )
    except Exception:
        # フォールバック：スキーマ差や古い環境でも確実に追記する
        if "route_cache" in listlayers(out_gpkg):
            old = gpd.read_file(out_gpkg, layer="route_cache")
            cols = sorted(set(old.columns) | set(route_cache_rec.columns))
            old = old.reindex(columns=cols)
            new = route_cache_rec.reindex(columns=cols)
            merged = gpd.GeoDataFrame(pd.concat([old, new], ignore_index=True), crs=4326)
            merged.to_file(
                out_gpkg, layer="route_cache",
                driver="GPKG", engine="pyogrio",
                layer_options={"OVERWRITE": "YES"}  # マージして上書き
            )
        else:
            route_cache_rec.to_file(
                out_gpkg, layer="route_cache",
                driver="GPKG", engine="pyogrio"
            )

    # 8) folium用座標
    route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]

# 9) Folium表示（キャッシュ・新規どちらでも）
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

へ変更

これで動くようになった

!ogrinfo -so outputs/shizuoka_hazard.gpkg route_latest

で動作も問題なし

import geopandas as gpd
gpd.read_file("outputs/shizuoka_hazard.gpkg", layer="route_latest")[["meters","eta_min_walk","eta_min_walk_safe","eta_min_drive"]]

も問題なし

再度 notebook での確認

import sys, os
sys.path.append(os.path.abspath("."))  # プロジェクト直下をパスに追加
# 以後 from evac.config import Cfg が通る

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString
import osmnx as ox
import networkx as nx
import folium
import os

GPKG = "outputs/shizuoka_hazard.gpkg"

# 1) 既存レイヤ読み込み（両方とも 4326）
coast   = gpd.read_file(GPKG, layer="coastline")
tsunami = gpd.read_file(GPKG, layer="tsunami_inundation")

# ----- パラメータ -----
# 自宅座標（lat, lon）←必ず置き換え
home_wgs = (34.728, 137.978)
# OSM道路の取得半径[m]（到達不能なら大きく：15000, 20000…）
GRAPH_RADIUS_M = 12000
# 避難所CSV（cp932想定）
SHELTER_CSV = "data/20250801_指定緊急避難場所一覧.CSV"
# ----------------------

# 2) 避難所CSV→点化（列名は自動検出）
df = pd.read_csv(SHELTER_CSV, encoding="cp932")
lat_candidates = 
lon_candidates = 
if not lat_candidates or not lon_candidates:
    raise ValueError("CSVに緯度/経度列が見つかりません。列名を教えてください。")
lat_col, lon_col = lat_candidates[0], lon_candidates[0]
df[lat_col] = pd.to_numeric(df[lat_col], errors="coerce")
df[lon_col] = pd.to_numeric(df[lon_col], errors="coerce")
df = df.dropna(subset=[lat_col, lon_col])
shelters = gpd.GeoDataFrame(df.copy(),
    geometry=gpd.points_from_xy(df[lon_col], df[lat_col], crs="EPSG:4326")
)

# 3) 安全な避難所 = 浸水域と非交差
#   （unionはメモリ負荷になるので、空間結合で判定）
safe_join = gpd.sjoin(shelters, tsunami[["geometry"]], how="left", predicate="intersects")
shelters["safe"] = safe_join["index_right"].isna()
safe_shelters = shelters[shelters["safe"]].copy()

# --- 4) OSM道路取得（徒歩網） ---
G = ox.graph_from_point(home_wgs, dist=GRAPH_RADIUS_M, network_type="walk", simplify=True, retain_all=False)
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True)

# u,v,key を“列”にする（重要）
edges = edges.reset_index()

# CRSを4326に
if edges.crs is None:
    edges = edges.set_crs(4326)
elif edges.crs.to_epsg() != 4326:
    edges = edges.to_crs(4326)
# --- 5) 浸水域と交差する道路を除去（sjoin→集計で一意化） ---
left = edges[["u","v","key","geometry"]].copy()   # edges は reset_index() 済み想定
right = tsunami[["geometry"]]

tag = gpd.sjoin(left, right, how="left", predicate="intersects")

blocked = (
    tag.assign(hit=tag["index_right"].notna())
       .groupby(["u","v","key"], as_index=False)["hit"].max()
       .rename(columns={"hit":"blocked"})
)

edges2 = edges.merge(blocked, on=["u","v","key"], how="left")
edges2["blocked"] = edges2["blocked"].fillna(False)

edges_keep = edges2.loc[~edges2["blocked"]].copy()

# ★ ここが重要：graph_from_gdfs に渡す前に MultiIndex に戻す
#   （key が float になっていたら int にしておくと安心）
if edges_keep["key"].isna().any():
    edges_keep["key"] = edges_keep["key"].fillna(0)

edges_keep["key"] = edges_keep["key"].astype(int)
edges_keep = edges_keep.set_index(["u","v","key"]).sort_index()

# --- 6) サブグラフ再構築（通行可のみ） ---
used_nodes = pd.unique(edges_keep.reset_index()[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()

# nodes_keep の geometry 列が欠けていたら補完
if "geometry" not in nodes_keep or nodes_keep["geometry"].isna().any():
    nodes_keep = nodes_keep.set_geometry(
        gpd.points_from_xy(nodes_keep["x"], nodes_keep["y"]), crs=edges_keep.crs
    )

# OSMnx のバージョン差に対応
try:
    G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))
except AttributeError:
    from osmnx import utils_graph as ug
    G2 = ug.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))


# 7) 出発・目的地ノード（最寄り“安全”避難所）
home_pt = Point(home_wgs[1], home_wgs[0])
# 最寄り安全避難所（直線距離）…なければ全体から選ぶ
TARGET = safe_shelters if len(safe_shelters) else shelters
METRIC = "EPSG:32654"# 静岡周辺 UTM
hp = gpd.GeoSeries([home_pt], crs=4326).to_crs(METRIC).iloc[0]
tp = TARGET.to_crs(METRIC).copy()
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = TARGET.loc[tp["__dist_m__"].idxmin(), "geometry"]

# ▶ 互換呼び分け
try:
    orig = ox.distance.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)
except AttributeError:
    # 新API名（環境によってはこちらが有効）
    orig = ox.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)


# 8) 最短経路（通行可能グラフ）
try:
    route = nx.shortest_path(G2, orig, dest, weight="length")
except nx.NetworkXNoPath:
    raise RuntimeError("到達経路なし。GRAPH_RADIUS_Mを増やすか、別の避難所で再試行を。")

# 9) GPKGへ追記
#    既にある場合は上書きしたいなら driverオプションを切り替える（ここは追記）
shelters.to_file(GPKG, layer="shelters_all", driver="GPKG")
safe_shelters.to_file(GPKG, layer="shelters_safe", driver="GPKG")
edges.to_file(GPKG, layer="roads_all", driver="GPKG")
edges_keep.to_file(GPKG, layer="roads_safe", driver="GPKG")

# 10) Folium地図
m = folium.Map(location=[home_wgs[0], home_wgs[1]], zoom_start=12, control_scale=True)
folium.GeoJson(tsunami.__geo_interface__, name="津波浸水想定(A40)").add_to(m)
folium.GeoJson(edges_keep.__geo_interface__, name="道路(浸水回避)").add_to(m)

# 経路を線で描画
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(route_latlon, weight=6, opacity=0.9, tooltip="避難経路").add_to(m)

# 出発・目的地
folium.Marker([home_wgs[0], home_wgs[1]], popup="現在地", icon=folium.Icon(color="green")).add_to(m)
folium.Marker([dest_pt.y, dest_pt.x], popup="避難所", icon=folium.Icon(color="blue")).add_to(m)

# 避難所レイヤ（安全のみ強調）
folium.GeoJson(safe_shelters.__geo_interface__, name="安全な避難所",
               marker=folium.CircleMarker(radius=4)).add_to(m)
folium.LayerControl().add_to(m)

# 保存（任意）
m.save("outputs/evac.html")
print("✓ 完了: GPKG追記 & outputs/evac.html を作成")

print("edges_keep index levels:", getattr(edges_keep.index, "nlevels", 1))  # → 3 が正解
print(edges_keep.index.names)  # → ['u','v','key'] が正解

# # ルートの延長[m]は G2 の edge 属性 'length' を合計（OSMnxが付けてくれてます）
# meters = sum(ox.utils_graph.get_route_edge_attributes(G2, route, 'length'))
# walk_kmh = 4.8   # 歩行速度（必要に応じて変更）
# drive_kmh = 20.0 # 災害時の車想定の控えめ値（使うなら）

# eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
# eta_drive_min = meters / (drive_kmh * 1000 / 60)
# print(f"距離: {meters:.0f} m, 徒歩ETA: {eta_walk_min:.1f} 分, 車ETA(仮): {eta_drive_min:.1f} 分")

# ▼ 距離[m]の取り方をバージョン互換で
try:
    # 新しめのOSMnx（推奨）: ルートをGDF化して length を合計
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    # フォールバック: NetworkXからエッジ属性 length を合計
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)  # MultiDiGraph: {key: attr_dict}
        if isinstance(ed, dict) and len(ed):
            meters += min(d.get("length", 0.0) for d in ed.values())

# ETA（必要なら速度を調整）
walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
print(f"距離: {meters:.0f} m, 徒歩ETA: {eta_walk_min:.1f} 分, 車ETA(仮): {eta_drive_min:.1f} 分")

# Folium のツールチップにも表示
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分"
).add_to(m)


from shapely.geometry import LineString
import geopandas as gpd

route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]
route_gdf = gpd.GeoDataFrame(
    {"meters":[meters], "eta_min_walk":[eta_walk_min], "eta_min_drive":[eta_drive_min]},
    geometry=[LineString(route_lonlat)], crs="EPSG:4326"
)
route_gdf.to_file("outputs/shizuoka_hazard.gpkg", layer="route_latest", driver="GPKG")

import geopandas as gpd
# ok = not route_gdf.intersects(tsunami.unary_union).iloc[0]
# print("交差なし✅" if ok else "⚠️ 交差あり（パラメータ見直し）")

# 旧: tsunami.unary_union  ← 警告が出る
# 新:
flood_union = tsunami.geometry.union_all()
ok = not route_gdf.intersects(flood_union).iloc[0]
print("交差なし✅" if ok else "⚠️ 交差あり（パラメータ見直し）")

# === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
from shapely.geometry import LineString
import geopandas as gpd

# 1) ルートのLineString（まだ属性は付けない）
route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
route_geom = LineString(route_lonlat)
route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# 2) 堅牢チェック（交差ゼロ＆境界近接なし）
#   union_all が使えない環境向けにフォールバックも用意
try:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
except Exception:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
route_m = route_gdf.to_crs(METRIC)
flood_m = flood_union.to_crs(METRIC)

# 交差長[m]
inter = gpd.overlay(route_m, flood_m, how="intersection")
inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# 近接（境界 ±2m 以内をNGにする例）
# near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty


if inside_len_m > 0 or near:
    raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                       "別避難所 or 半径拡大で再探索してください。")

# 3) 距離[m]→ETA（バージョン互換）
try:
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)
        meters += min(d.get("length", 0.0) for d in ed.values())

walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# 4) GPKG に“合格ルート”だけ属性付きで保存
out_gpkg = "outputs/shizuoka_hazard.gpkg"
route_gdf = gpd.GeoDataFrame(
    {
        "meters":[meters],
        "eta_min_walk":[eta_walk_min],
        "eta_min_walk_safe":[eta_walk_safe],
        "eta_min_drive":[eta_drive_min],
    },
    geometry=[route_geom], crs="EPSG:4326",
)
route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
print("✅ 交差なし・保存完了")

# 5) （任意）Foliumに表示するならここで
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

!ogrinfo -so outputs/shizuoka_hazard.gpkg route_latest

import geopandas as gpd
gpd.read_file("outputs/shizuoka_hazard.gpkg", layer="route_latest")[["meters","eta_min_walk","eta_min_walk_safe","eta_min_drive"]]

までは動く

どこからバグったかを調べたら
まとめ（差し込む位置の目安）
1. 目的地 dest_pt 決定直後に「キャッシュ照会ブロック」を挿入
2. use_cached_route が False なら従来の「計算→検証→保存」ブロックへ進む
3. 保存直後に route_cache へ追記
4. True なら保存はスキップし、地図表示へ直行（または再検証のみ実施）
これで「既知の組み合わせは即復帰／未知は新規計算」の運用ができます。

これはモジュールにして読み込んだ方が良い？
というあたりでバグ発生

つまりここから無理にモジュールしたのが原因で先にプロトタイプを作ればいい

ルート再利用性を高めたいので、ルート結果をキー付きでGPKGにキャッシュする機能を追加

# === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
from shapely.geometry import LineString
import geopandas as gpd

# 1) ルートのLineString（まだ属性は付けない）
route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
route_geom = LineString(route_lonlat)
route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# 2) 堅牢チェック（交差ゼロ＆境界近接なし）
#   union_all が使えない環境向けにフォールバックも用意
try:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
except Exception:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
route_m = route_gdf.to_crs(METRIC)
flood_m = flood_union.to_crs(METRIC)

# 交差長[m]
inter = gpd.overlay(route_m, flood_m, how="intersection")
inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# 近接（境界 ±2m 以内をNGにする例）
# near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty


if inside_len_m > 0 or near:
    raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                       "別避難所 or 半径拡大で再探索してください。")

# 3) 距離[m]→ETA（バージョン互換）
try:
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)
        meters += min(d.get("length", 0.0) for d in ed.values())

walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# 4) GPKG に“合格ルート”だけ属性付きで保存
out_gpkg = "outputs/shizuoka_hazard.gpkg"
route_gdf = gpd.GeoDataFrame(
    {
        "meters":[meters],
        "eta_min_walk":[eta_walk_min],
        "eta_min_walk_safe":[eta_walk_safe],
        "eta_min_drive":[eta_drive_min],
    },
    geometry=[route_geom], crs="EPSG:4326",
)
route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
print("✅ 交差なし・保存完了")

# 5) （任意）Foliumに表示するならここで
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

のコードを
# === ルート検証 → ETA → 保存（キャッシュ対応） ===
from shapely.geometry import LineString, Point
import geopandas as gpd
import hashlib, json
from datetime import datetime
from fiona import listlayers
import osmnx as ox
import networkx as nx

GPKG = "outputs/shizuoka_hazard.gpkg"
METRIC = "EPSG:32654"  # 静岡周辺

# 0) 津波 union はキャッシュを優先（なければ作る）
if "tsunami_union" in listlayers(GPKG):
    flood_union = gpd.read_file(GPKG, layer="tsunami_union")
else:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs=tsunami.crs)
    flood_union.to_file(GPKG, layer="tsunami_union", driver="GPKG")

# 0') ハザードの簡易メタ（A40が変わったらルートキャッシュ無効化に使う）
ti = gpd.read_file(GPKG, layer="tsunami_inundation")
hazard_meta = f'{len(ti)}|{",".join(map(str, ti.total_bounds))}'

# 1) キャッシュキー生成（端点・モード・半径・ハザードメタ）
def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000, hazard_meta=None):
    payload = {"home": home_wgs, "dest": dest_wgs, "mode": mode, "r": radius_m, "hz": hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk",
                radius_m=GRAPH_RADIUS_M, hazard_meta=hazard_meta)

# 2) ルートキャッシュにヒットすれば使う
route_gdf = None
use_cached_route = False
if "route_cache" in listlayers(GPKG):
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"] == key]
    if len(hit):
        route_gdf = hit.set_crs(4326)
        meters = float(route_gdf["meters"].iloc[0])
        eta_walk_min  = float(route_gdf["eta_min_walk"].iloc[0])
        eta_walk_safe = float(route_gdf["eta_min_walk_safe"].iloc[0])
        eta_drive_min = float(route_gdf["eta_min_drive"].iloc[0])
        coords = list(route_gdf.geometry.iloc[0].coords)
        route_latlon = [(lat, lon) for lon, lat in coords]  # folium用 (lat,lon)
        use_cached_route = True
        print("↩︎ ルート: キャッシュヒット")

if not use_cached_route:
    # 3) これまで通りの経路（'route'）から LineString を作成
    route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
    route_geom = LineString(route_lonlat)
    route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

    # 4) 堅牢チェック（交差ゼロ＆境界近接なし）
    route_m = route_gdf.to_crs(METRIC)
    flood_m = flood_union.to_crs(METRIC)

    inter = gpd.overlay(route_m, flood_m, how="intersection")
    inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

    buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
    near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

    if inside_len_m > 0 or near:
        raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                           "別避難所 or 半径拡大で再探索してください。")

    # 5) 距離[m]→ETA
    try:
        route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
        meters = float(route_edges["length"].sum())
    except Exception:
        meters = 0.0
        for u, v in zip(route[:-1], route[1:]):
            ed = G2.get_edge_data(u, v)
            meters += min(d.get("length", 0.0) for d in ed.values())

    walk_kmh = 4.8
    drive_kmh = 20.0
    eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
    eta_drive_min = meters / (drive_kmh * 1000 / 60)
    eta_walk_safe = eta_walk_min * 1.3

    # 6) route_latest へ保存（検証済みのみ）
    out_gpkg = GPKG
    route_gdf = gpd.GeoDataFrame(
        {
            "meters":[meters],
            "eta_min_walk":[eta_walk_min],
            "eta_min_walk_safe":[eta_walk_safe],
            "eta_min_drive":[eta_drive_min],
            "created_at":[datetime.utcnow().isoformat(timespec="seconds")+"Z"],
        },
        geometry=[route_geom], crs="EPSG:4326",
    )
    route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
    print("✅ 交差なし・保存完了")

    # 7) ルートキャッシュに追記
    route_cache_rec = route_gdf.copy()
    route_cache_rec["rid"] = key
    try:
        route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG", if_exists="append")
    except TypeError:
        # 古いGeoPandasのためのフォールバック
        if "route_cache" in listlayers(out_gpkg):
            old = gpd.read_file(out_gpkg, layer="route_cache")
            gpd.GeoDataFrame(pd.concat([old, route_cache_rec], ignore_index=True), crs=4326)\
               .to_file(out_gpkg, layer="route_cache", driver="GPKG")
        else:
            route_cache_rec.to_file(out_gpkg, layer="route_cache", driver="GPKG")

    # 8) folium用座標
    route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]

# 9) Folium表示（キャッシュ・新規どちらでも）
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

に書き換える

津波避難関連のモジュール化(失敗）

津波避難関連のモジュール化

touch config.py

これで定数ひとまとめ

from dataclasses import dataclass

@dataclass(frozen=True)
class Cfg:
    GPKG = "outputs/shizuoka_hazard.gpkg"
    EPSG_WGS84 = 4326
    EPSG_METRIC = 32654           # 静岡周辺
    GRAPH_RADIUS_M = 12000
    WALK_KMH = 4.8
    DRIVE_KMH = 20.0
    SAFETY_FACTOR = 1.3

touch data_io.py

内容は

import geopandas as gpd
from fiona import listlayers
from .config import Cfg

def load_layers():
    coast   = gpd.read_file(Cfg.GPKG, layer="coastline")
    tsunami = gpd.read_file(Cfg.GPKG, layer="tsunami_inundation")
    return coast, tsunami

def get_tsunami_union(tsunami):
    if "tsunami_union" in listlayers(Cfg.GPKG):
        return gpd.read_file(Cfg.GPKG, layer="tsunami_union")
    fu = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs=tsunami.crs)
    fu.to_file(Cfg.GPKG, layer="tsunami_union", driver="GPKG")
    return fu

def write_route_latest(route_gdf):
    route_gdf.to_file(Cfg.GPKG, layer="route_latest", driver="GPKG")

touch hazard.py

これは（検証ロジック）

import geopandas as gpd
def validate_route(route_gdf, flood_union, epsg_metric, near_m=2.0):
    route_m = route_gdf.to_crs(epsg_metric)
    flood_m = flood_union.to_crs(epsg_metric)

    inter = gpd.overlay(route_m, flood_m, how="intersection")
    inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

    buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(near_m), crs=epsg_metric)
    near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty
    return inside_len_m, near

touch network.py

これはOSM取得＆遮断

import osmnx as ox, geopandas as gpd, pandas as pd

def build_safe_graph(center_latlon, radius_m, tsunami):
    G = ox.graph_from_point(center_latlon, dist=radius_m, network_type="walk", simplify=True)
    nodes, edges = ox.graph_to_gdfs(G)
    edges = edges.reset_index()
    if edges.crs is None: edges = edges.set_crs(4326)
    elif edges.crs.to_epsg()!=4326: edges = edges.to_crs(4326)

    tag = gpd.sjoin(edges[["u","v","key","geometry"]], tsunami[["geometry"]],
                    how="left", predicate="intersects")
    blocked = (tag.assign(hit=tag["index_right"].notna())
                  .groupby(["u","v","key"], as_index=False)["hit"].max()
                  .rename(columns={"hit":"blocked"}))
    edges2 = edges.merge(blocked, on=["u","v","key"], how="left")
    edges2["blocked"] = edges2["blocked"].fillna(False)
    edges_keep = edges2.loc[~edges2["blocked"]].copy()
    edges_keep["key"] = edges_keep["key"].fillna(0).astype(int)
    edges_keep = edges_keep.set_index(["u","v","key"]).sort_index()
    used_nodes = pd.unique(edges_keep.reset_index()[["u","v"]].values.ravel())
    nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()
    try:
        G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)
    except AttributeError:
        from osmnx import utils_graph as ug
        G2 = ug.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)
    return G2, nodes_keep, edges_keep

touch routing.py

これは経路・ETA

import osmnx as ox, networkx as nx, geopandas as gpd
from shapely.geometry import LineString
from .config import Cfg

def nearest_node(G, x, y):
    try: return ox.distance.nearest_nodes(G, X=x, Y=y)
    except AttributeError: return ox.nearest_nodes(G, X=x, Y=y)

def route_and_eta(G2, route_nodes):
    try:
        route_edges = ox.routing.route_to_gdf(G2, route_nodes, True)
        meters = float(route_edges["length"].sum())
    except Exception:
        meters = 0.0
        for u,v in zip(route_nodes[:-1], route_nodes[1:]):
            ed = G2.get_edge_data(u,v); meters += min(d.get("length",0.0) for d in ed.values())
    walk = meters/(Cfg.WALK_KMH*1000/60)
    drive= meters/(Cfg.DRIVE_KMH*1000/60)
    return meters, walk, drive

def make_route_gdf(G2, route_nodes, meters, walk, drive):
    latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route_nodes]
    lonlat = [(x,y) for y,x in latlon]
    gdf = gpd.GeoDataFrame(
        {"meters":[meters],"eta_min_walk":[walk],
         "eta_min_walk_safe":[walk*Cfg.SAFETY_FACTOR],
         "eta_min_drive":[drive]},
        geometry=[LineString(lonlat)], crs=4326
    )
    return gdf, latlon

touch cache.py

これはキー生成＆保存

import json, hashlib, geopandas as gpd
from fiona import listlayers
from .config import Cfg

def route_key(home_wgs, dest_wgs, mode="walk", radius_m=Cfg.GRAPH_RADIUS_M, hazard_meta=None):
    payload = {"home":home_wgs, "dest":dest_wgs, "mode":mode, "r":radius_m, "hz":hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

def load_cached_route(rid):
    if "route_cache" not in listlayers(Cfg.GPKG): return None
    gdf = gpd.read_file(Cfg.GPKG, layer="route_cache")
    hit = gdf.loc[gdf["rid"]==rid]
    return hit if len(hit) else None

def save_route_cache(route_gdf, rid):
    route_gdf = route_gdf.copy(); route_gdf["rid"] = rid
    try:
        route_gdf.to_file(Cfg.GPKG, layer="route_cache", driver="GPKG", if_exists="append")
    except TypeError:
        # 旧版対応
        import pandas as pd
        try:
            old = gpd.read_file(Cfg.GPKG, layer="route_cache")
            gpd.GeoDataFrame(pd.concat([old, route_gdf], ignore_index=True), crs=4326)\
               .to_file(Cfg.GPKG, layer="route_cache", driver="GPKG")
        except Exception:
            route_gdf.to_file(Cfg.GPKG, layer="route_cache", driver="GPKG")

 touch viz.py

これは表示するモジュール

import folium, geopandas as gpd

def add_layers_to_map(m, tsunami, edges_keep, route_latlon, tooltip):
    folium.GeoJson(tsunami.__geo_interface__, name="津波浸水想定").add_to(m)
    folium.GeoJson(edges_keep.__geo_interface__, name="道路(浸水回避)").add_to(m)
    folium.PolyLine(route_latlon, weight=6, opacity=0.9, tooltip=tooltip).add_to(m)
    folium.LayerControl().add_to(m)

あとは
Notebook で

from evac.config import Cfg
from evac.data_io import load_layers, get_tsunami_union, write_route_latest
from evac.network import build_safe_graph
from evac.routing import nearest_node, route_and_eta, make_route_gdf
from evac.hazard import validate_route
from evac.cache import route_key, load_cached_route, save_route_cache

coast, tsunami = load_layers()
flood_union = get_tsunami_union(tsunami)

home_wgs = (34.728, 137.978)
G2, nodes_keep, edges_keep = build_safe_graph(home_wgs, Cfg.GRAPH_RADIUS_M, tsunami)

# 目的地は別途選定済みとする: dest_pt
orig = nearest_node(G2, x=home_wgs[1], y=home_wgs[0])
dest = nearest_node(G2, x=dest_pt.x, y=dest_pt.y)

rid = route_key(home_wgs, (dest_pt.y, dest_pt.x))
cached = load_cached_route(rid)
if cached is not None:
    route_gdf = cached.set_crs(4326)
else:
    route_nodes = nx.shortest_path(G2, orig, dest, weight="length")
    meters, walk, drive = route_and_eta(G2, route_nodes)
    route_gdf, route_latlon = make_route_gdf(G2, route_nodes, meters, walk, drive)
    inside_len_m, near = validate_route(route_gdf, flood_union, Cfg.EPSG_METRIC)
    if inside_len_m>0 or near: raise RuntimeError("ルート不適合")
    write_route_latest(route_gdf)
    save_route_cache(route_gdf, rid)
/
を実行したが

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[17], line 1
----> 1 from evac.config import Cfg
      2 from evac.data_io import load_layers, get_tsunami_union, write_route_latest
      3 from evac.network import build_safe_graph

ModuleNotFoundError: No module named 'evac'

これは
evac パッケージがPythonの検索パスに無いだけ

mkdir -p evac
touch evac/__init__.py

でファイルを作成

とりあえず中身は空でも良いらしい

あとは notebookで

import sys, os
sys.path.append(os.path.abspath("."))  # プロジェクト直下をパスに追加
# 以後 from evac.config import Cfg が通る

を記述

しかしうまくいかないため
Main.pyにして実行してみる

mv *.py evac 
mv evac/main.py .
python main.py

これで実行

今度は

Traceback (most recent call last):
  File "/Users/snowpool/aw10s/FamilySafetyMap/main.py", line 16, in <module>
    dest = nearest_node(G2, x=dest_pt.x, y=dest_pt.y)
                              ^^^^^^^
NameError: name 'dest_pt' is not defined

となる

dest_pt（目的地の点ジオメトリ）を作っていないのが原因

修正案は２つ
A) 手早い版（直線距離で最寄りの安全避難所を選ぶ）
B) 推奨版（上位N候補から“ネットワーク距離最短”を選ぶ）

Bを実行する

import geopandas as gpd
import networkx as nx
from shapely.geometry import Point

を追記

もともとの

# 目的地は別途選定済みとする: dest_pt
orig = nearest_node(...)
dest = nearest_node(... dest_pt ...)
rid = route_key(...)

は削除

代わりに pick_best_shelter_by_network で route_nodes と dest_pt を得ます。
pick_best_shelter_by_network をまだモジュール化していないなら、main.py に一時定義してから使ってOKです（前に渡した実装をそのまま貼り付け）。 もしくは evac/routing_extras.py を作って関数を置き、from evac.routing_extras import … で読み込む形にしてください。
キャッシュヒット時は route_gdf をそのまま使えるので、再計算・検証はスキップ（必要なら最新 tsunami_union で再検証だけ行うのも可）。

from evac.config import Cfg
from evac.data_io import load_layers, get_tsunami_union, write_route_latest
from evac.network import build_safe_graph
from evac.routing import nearest_node, route_and_eta, make_route_gdf
from evac.hazard import validate_route
from evac.cache import route_key, load_cached_route, save_route_cache

coast, tsunami = load_layers()
flood_union = get_tsunami_union(tsunami)

home_wgs = (34.728, 137.978)
G2, nodes_keep, edges_keep = build_safe_graph(home_wgs, Cfg.GRAPH_RADIUS_M, tsunami)

# 目的地は別途選定済みとする: dest_pt
orig = nearest_node(G2, x=home_wgs[1], y=home_wgs[0])
dest = nearest_node(G2, x=dest_pt.x, y=dest_pt.y)

rid = route_key(home_wgs, (dest_pt.y, dest_pt.x))
cached = load_cached_route(rid)
if cached is not None:
    route_gdf = cached.set_crs(4326)
else:
    route_nodes = nx.shortest_path(G2, orig, dest, weight="length")
    meters, walk, drive = route_and_eta(G2, route_nodes)
    route_gdf, route_latlon = make_route_gdf(G2, route_nodes, meters, walk, drive)
    inside_len_m, near = validate_route(route_gdf, flood_union, Cfg.EPSG_METRIC)
    if inside_len_m>0 or near: raise RuntimeError("ルート不適合")
    write_route_latest(route_gdf)
    save_route_cache(route_gdf, rid)

のコードを

# --- ここまで同じ ---
from evac.config import Cfg
from evac.data_io import load_layers, get_tsunami_union, write_route_latest
from evac.network import build_safe_graph
from evac.routing import nearest_node, route_and_eta, make_route_gdf
from evac.hazard import validate_route
from evac.cache import route_key, load_cached_route, save_route_cache

import geopandas as gpd
import networkx as nx
from shapely.geometry import Point

# pick_best_shelter_by_network をどこかに実装済みなら import
# 例: evac.routing_extras に置いた場合
# from evac.routing_extras import pick_best_shelter_by_network
# 未実装なら、ひとまず前に貼った関数を main.py に一時定義してもOK

coast, tsunami = load_layers()
flood_union = get_tsunami_union(tsunami)

home_wgs = (34.728, 137.978)
G2, nodes_keep, edges_keep = build_safe_graph(home_wgs, Cfg.GRAPH_RADIUS_M, tsunami)

# ====== ここを「目的地は別途選定済み」から 置き換え ======
# 避難所レイヤ（安全優先 → なければ全体）
try:
    TARGET = gpd.read_file(Cfg.GPKG, layer="shelters_safe")
    if TARGET.empty:
        raise IOError
except Exception:
    TARGET = gpd.read_file(Cfg.GPKG, layer="shelters_all")

if TARGET.empty:
    raise RuntimeError("避難所レイヤが空です（shelters_safe / shelters_all）。作成を確認してください。")

# 出発ノード
home_pt = Point(home_wgs[1], home_wgs[0])  # (lon, lat)
orig = nearest_node(G2, x=home_pt.x, y=home_pt.y)

# 上位N候補からネットワーク距離最短の避難所を選定
# ※ pick_best_shelter_by_network(G2, orig, TARGET, N=5) が定義済み前提
route_nodes, meters_tmp, dest_pt = pick_best_shelter_by_network(G2, orig, TARGET, N=5)

# 以降のキャッシュキー生成に dest_pt を使用
rid = route_key(home_wgs, (dest_pt.y, dest_pt.x))
cached = load_cached_route(rid)
# ====== 置き換えここまで ======

if cached is not None:
    route_gdf = cached.set_crs(4326)
else:
    # 既に route_nodes があるので最短経路は再計算不要
    meters, walk, drive = route_and_eta(G2, route_nodes)
    route_gdf, route_latlon = make_route_gdf(G2, route_nodes, meters, walk, drive)
    inside_len_m, near = validate_route(route_gdf, flood_union, Cfg.EPSG_METRIC)
    if inside_len_m > 0 or near:
        raise RuntimeError("ルート不適合")
    write_route_latest(route_gdf)
    save_route_cache(route_gdf, rid)

に変更

しかし

Traceback (most recent call last):
  File "/Users/snowpool/aw10s/FamilySafetyMap/main.py", line 16, in <module>
    dest = nearest_node(G2, x=dest_pt.x, y=dest_pt.y)
                              ^^^^^^^
NameError: name 'dest_pt' is not defined
snowpool@MacBookAir FamilySafetyMap % python main.py 
Traceback (most recent call last):
  File "/Users/snowpool/aw10s/FamilySafetyMap/main.py", line 42, in <module>
    route_nodes, meters_tmp, dest_pt = pick_best_shelter_by_network(G2, orig, TARGET, N=5)
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
NameError: name 'pick_best_shelter_by_network' is not defined

となる

エラー原因は2つ：
1. dest_pt を作る前に nearest_node(… dest_pt …) を呼んでいる（古い行が残っている）
2. pick_best_shelter_by_network をまだ定義/インポートしていない

このため

evac/routing_extras.py

を作成

# evac/routing_extras.py
import numpy as np
import networkx as nx
import osmnx as ox
import geopandas as gpd

def pick_best_shelter_by_network(G2, orig_node, shelters_gdf, N=5, epsg_metric=32654):
    """直線距離の上位N候補から、ネットワーク距離が最短の避難所を選ぶ。"""
    # 出発点の点（lon,lat）
    hp = gpd.GeoSeries([ (G2.nodes[orig_node]["x"], G2.nodes[orig_node]["y"]) ],
                       crs=4326).to_crs(epsg_metric).iloc[0]
    cand = shelters_gdf.to_crs(epsg_metric).copy()
    cand["__d__"] = cand.geometry.distance(hp)
    cand = shelters_gdf.loc[cand.sort_values("__d__").head(N).index].copy()

    best_route = None
    best_len = np.inf
    best_dest_pt = None

    for _, row in cand.iterrows():
        try:
            try:
                dnode = ox.distance.nearest_nodes(G2, X=row.geometry.x, Y=row.geometry.y)
            except AttributeError:
                dnode = ox.nearest_nodes(G2, X=row.geometry.x, Y=row.geometry.y)
            r = nx.shortest_path(G2, orig_node, dnode, weight="length")
            # 経路長[m]
            try:
                meters = float(ox.routing.route_to_gdf(G2, r, True)["length"].sum())
            except Exception:
                meters = 0.0
                for u, v in zip(r[:-1], r[1:]):
                    ed = G2.get_edge_data(u, v)
                    meters += min(d.get("length", 0.0) for d in ed.values())
            if meters < best_len:
                best_len = meters
                best_route = r
                best_dest_pt = row.geometry
        except nx.NetworkXNoPath:
            continue

    if best_route is None:
        raise RuntimeError("到達可能な候補が見つかりません。半径またはNを増やしてください。")
    return best_route, best_len, best_dest_pt

次に
main.py の該当箇所を置き換え

# ❌ 削除
dest = nearest_node(G2, x=dest_pt.x, y=dest_pt.y)
rid = route_key(home_wgs, (dest_pt.y, dest_pt.x))

コードを変更

# main.py
from evac.config import Cfg
from evac.data_io import load_layers, get_tsunami_union, write_route_latest
from evac.network import build_safe_graph
from evac.routing import nearest_node, route_and_eta, make_route_gdf
from evac.hazard import validate_route
from evac.cache import route_key, load_cached_route, save_route_cache
from evac.routing_extras import pick_best_shelter_by_network   # ← 追加

import geopandas as gpd
import networkx as nx
from shapely.geometry import Point

coast, tsunami = load_layers()
flood_union = get_tsunami_union(tsunami)

home_wgs = (34.728, 137.978)
G2, nodes_keep, edges_keep = build_safe_graph(home_wgs, Cfg.GRAPH_RADIUS_M, tsunami)

# 避難所レイヤ（安全があれば優先）
try:
    TARGET = gpd.read_file(Cfg.GPKG, layer="shelters_safe")
    if TARGET.empty:
        raise IOError
except Exception:
    TARGET = gpd.read_file(Cfg.GPKG, layer="shelters_all")
if TARGET.empty:
    raise RuntimeError("避難所レイヤが空です（shelters_safe / shelters_all）。")

# 出発ノード
home_pt = Point(home_wgs[1], home_wgs[0])
orig = nearest_node(G2, x=home_pt.x, y=home_pt.y)

# ✅ 上位Nからネットワーク距離最短を選定（ここで dest_pt が得られる）
route_nodes, meters_tmp, dest_pt = pick_best_shelter_by_network(G2, orig, TARGET, N=5)

# 以降のキャッシュキー生成は dest_pt を使う
rid = route_key(home_wgs, (dest_pt.y, dest_pt.x))
cached = load_cached_route(rid)

if cached is not None:
    print("↩︎ ルート: キャッシュヒット")
    route_gdf = cached.set_crs(4326)
else:
    meters, walk, drive = route_and_eta(G2, route_nodes)
    route_gdf, route_latlon = make_route_gdf(G2, route_nodes, meters, walk, drive)
    inside_len_m, near = validate_route(route_gdf, flood_union, Cfg.EPSG_METRIC)
    if inside_len_m > 0 or near:
        raise RuntimeError("ルート不適合: 浸水域と交差/近接あり")
    write_route_latest(route_gdf)
    save_route_cache(route_gdf, rid)

tsunami_union を保存 → 交差チェックをそれに切替

tsunami_union を保存 → 交差チェックをそれに切替
上位N→ネットワーク最短の選定を導入（経路品質UP）
アイソクローンのレイヤをGPKG/地図に追加（“どこまで行けるか”の可視化）
警告フラグをUIに表示
GraphDB CSV出力で家族アプリへの橋渡し

を行う

1) tsunami_union をGPKGにキャッシュ（毎回のunion計算を省略）
これで気になったのが
このtsunami_union をGPKGにキャッシュしたものは現在地、もしくはすでに計測したもののみ有効？
もし新規の場所で避難することがあれば対応したいので
すでにあるルートならキャッシュ、
そうでないのなら新規ルートを使うようにしたい

結論：tsunami_union は“現在地に依存しない”全域キャッシュです。
静岡県の「津波浸水想定（A40）」レイヤ全体を union したものなので、同じデータ（A40）がソースならどの場所でも再利用可能。新規の場所でもそのまま使えます。再生成が必要なのは、A40を更新・差し替えたときだけです。

A) tsunami_union（全域）のキャッシュ方針
* 有効範囲：A40レイヤ（静岡版）全域。現在地・ルートに非依存。
* 再生成条件：A40が変わった可能性があるとき（例：ファイル更新、feature数やBBOXが変化）。
* 高速化：ルート検証時は tsunami_union をその場のAOIでクリップして使うと軽いです。

例：キャッシュの検証＆必要なら再生成

import geopandas as gpd
from fiona import listlayers

GPKG = "outputs/shizuoka_hazard.gpkg"

tsunami = gpd.read_file(GPKG, layer="tsunami_inundation")
feat_count = len(tsunami)
bbox = tuple(tsunami.total_bounds)  # (xmin,ymin,xmax,ymax)

# 既存 union があるか＆メタデータが一致するか
need_rebuild = True
if "tsunami_union" in listlayers(GPKG):
    union = gpd.read_file(GPKG, layer="tsunami_union")
    # 簡易同一性チェック（必要なら厳格化）
    if "meta_feat" in union.columns and "meta_bbox" in union.columns:
        need_rebuild = not (int(union.loc[0,"meta_feat"]) == feat_count and
                            tuple(map(float, union.loc[0,"meta_bbox"].split(","))) == bbox)

if need_rebuild:
    flood_union = gpd.GeoDataFrame(
        {"meta_feat":[feat_count], "meta_bbox":[",".join(map(str,bbox))]},
        geometry=[tsunami.geometry.union_all()], crs=tsunami.crs
    )
    flood_union.to_file(GPKG, layer="tsunami_union", driver="GPKG")
else:
    flood_union = gpd.read_file(GPKG, layer="tsunami_union")

例：ルート検証時は AOI で軽量化

# route_gdf: ルート LineString (EPSG:4326)
AOI_m = 15000  # 15km
aoi = route_gdf.to_crs("EPSG:32654").buffer(AOI_m).to_crs(4326)
flood_union_aoi = gpd.overlay(flood_union, gpd.GeoDataFrame(geometry=aoi, crs=4326), how="intersection")
# 以後は flood_union_aoi で交差・近接判定（軽い）

B) 「既に計測したルート」はルート側をキャッシュ（任意）
tsunami_union は“全域・一回”でOKですが、ルートの結果は現在地や避難所・パラメータに依存するので、キー付きでGPKGに保存して再利用が有効です。
例：ルート結果のキャッシュ（キー = 端点・モード・半径）

import hashlib, json, geopandas as gpd

def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000):
    payload = {"home":home_wgs, "dest":dest_wgs, "mode":mode, "r":radius_m}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk", radius_m=GRAPH_RADIUS_M)

# 取り出し（あれば使う）
try:
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"]==key]
    if len(hit):
        route_gdf = hit.set_crs(4326)  # そのまま使える
    else:
        raise IOError
except Exception:
    # 計算 → route_gdf を作ったら保存
    route_gdf["rid"] = key
    route_gdf.to_file(GPKG, layer="route_cache", driver="GPKG", mode="a")  # 追記

A40をアップデートしたら、route_cacheは無効化（レイヤを削除、またはキーにA40のメタ情報を混ぜる）するのが安全です。

ということで

ルート再利用性を高めたい場合は、ルート結果をキー付きでGPKGにキャッシュする機能を追加したい

import hashlib, json, geopandas as gpd

def route_key(home_wgs, dest_wgs, mode="walk", radius_m=12000):
    payload = {"home":home_wgs, "dest":dest_wgs, "mode":mode, "r":radius_m}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk", radius_m=GRAPH_RADIUS_M)

# 取り出し（あれば使う）
try:
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"]==key]
    if len(hit):
        route_gdf = hit.set_crs(4326)  # そのまま使える
    else:
        raise IOError
except Exception:
    # 計算 → route_gdf を作ったら保存
    route_gdf["rid"] = key
    route_gdf.to_file(GPKG, layer="route_cache", driver="GPKG", mode="a")  # 追記

これはどこに追記する？
入れる場所は**「目的地（dest_pt）が決まった直後、経路計算（nx.shortest_path）に入る前」です。
ヒットしたらそのまま使い、ミスなら今までどおり計算→検証→保存しつつキャッシュへ追記**します

# ==== ルートキャッシュ（GPKG: route_cache） ====
import hashlib, json, geopandas as gpd
from fiona import listlayers

GPKG = "outputs/shizuoka_hazard.gpkg"

def route_key(home_wgs, dest_wgs, mode="walk", radius_m=GRAPH_RADIUS_M, hazard_meta=None):
    payload = {"home": home_wgs, "dest": dest_wgs, "mode": mode, "r": radius_m, "hz": hazard_meta}
    return hashlib.md5(json.dumps(payload, sort_keys=True).encode()).hexdigest()

# A40が変わったら無効化できるように簡易メタを混ぜる（件数＋BBOX）
hazard_meta = None
try:
    if "tsunami_union" in listlayers(GPKG):
        # union保存時にメタ列を持たせていればそれを使う／無ければ tsunami_inundation から作る
        tu = gpd.read_file(GPKG, layer="tsunami_union")
        if {"meta_feat","meta_bbox"} <= set(tu.columns):
            hazard_meta = f'{int(tu.loc[0,"meta_feat"])}|{tu.loc[0,"meta_bbox"]}'
    if hazard_meta is None:
        ti = gpd.read_file(GPKG, layer="tsunami_inundation")
        hazard_meta = f'{len(ti)}|{",".join(map(str, ti.total_bounds))}'
except Exception:
    pass

key = route_key(home_wgs, (dest_pt.y, dest_pt.x), mode="walk", radius_m=GRAPH_RADIUS_M, hazard_meta=hazard_meta)

route_gdf_cached = None
if "route_cache" in listlayers(GPKG):
    cache = gpd.read_file(GPKG, layer="route_cache")
    hit = cache.loc[cache["rid"] == key]
    if len(hit):
        route_gdf_cached = hit.set_crs(4326)
        print("↩︎ ルート: キャッシュヒット")

# ヒット時の扱い：後続の検証/表示にそのまま使えるよう最低限の変数を整備
if route_gdf_cached is not None:
    route_gdf = route_gdf_cached[["meters","eta_min_walk","eta_min_walk_safe","eta_min_drive","geometry"]].copy()
    # route_latlon / meters 等は既存列から復元（既に保存してある前提）
    meters = float(route_gdf["meters"].iloc[0])
    eta_walk_min  = float(route_gdf["eta_min_walk"].iloc[0])
    eta_walk_safe = float(route_gdf["eta_min_walk_safe"].iloc[0])
    eta_drive_min = float(route_gdf["eta_min_drive"].iloc[0])
    # folium描画用
    coords = list(route_gdf.geometry.iloc[0].coords)
    route_latlon = [(lat, lon) for lon, lat in coords]
    # 以降の「堅牢チェック→保存」はスキップしてOK（または最新A40で再チェックしたいならそのまま通しても良い）
    use_cached_route = True
else:
    use_cached_route = False

その後のフローはこうしてください：
* use_cached_route == True のとき：
* すでに route_gdf / meters / eta_* / route_latlon が揃っているので、堅牢チェックと保存をスキップして Folium 表示だけ行う（必要なら最新 tsunami_union で再チェックだけ通してもOK）。
* use_cached_route == False のとき：
* これまでどおり nx.shortest_path → 堅牢チェック → ETA計算 → route_latest 保存。
* さらに route_cache に追記します（下の1行を“保存直後”に追加）。

# route_latest を保存した直後にキャッシュへ追記
route_gdf_cache = route_gdf.copy()
route_gdf_cache["rid"] = key
# append（GeoPandas 0.14+）：if_exists="append" が使えます
try:
    route_gdf_cache.to_file(GPKG, layer="route_cache", driver="GPKG", if_exists="append")
except TypeError:
    # 古い版は append がないので自前で結合して上書き
    if "route_cache" in listlayers(GPKG):
        old = gpd.read_file(GPKG, layer="route_cache")
        gpd.GeoDataFrame(pd.concat([old, route_gdf_cache], ignore_index=True), crs=4326)\
           .to_file(GPKG, layer="route_cache", driver="GPKG")
    else:
        route_gdf_cache.to_file(GPKG, layer="route_cache", driver="GPKG")

まとめ（差し込む位置の目安）
1. 目的地 dest_pt 決定直後に「キャッシュ照会ブロック」を挿入
2. use_cached_route が False なら従来の「計算→検証→保存」ブロックへ進む
3. 保存直後に route_cache へ追記
4. True なら保存はスキップし、地図表示へ直行（または再検証のみ実施）
これで「既知の組み合わせは即復帰／未知は新規計算」の運用ができます。

これはモジュールにして読み込んだ方が良い？

結論：はい。モジュール化して読み込むのが吉です。
理由は「再利用・テスト・差し替え（例：県を静岡→他県）」が一気に楽になるから。

貼り付け用の完成ブロック

貼り付け用の完成ブロック（そのまま置換OK）

# === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
from shapely.geometry import LineString
import geopandas as gpd

# 1) ルートのLineString（まだ属性は付けない）
route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
route_geom = LineString(route_lonlat)
route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# 2) 堅牢チェック（交差ゼロ＆境界近接なし）
#   union_all が使えない環境向けにフォールバックも用意
try:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
except Exception:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
route_m = route_gdf.to_crs(METRIC)
flood_m = flood_union.to_crs(METRIC)

# 交差長[m]
inter = gpd.overlay(route_m, flood_m, how="intersection")
inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# 近接（境界 ±2m 以内をNGにする例）
near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty

if inside_len_m > 0 or near:
    raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                       "別避難所 or 半径拡大で再探索してください。")

# 3) 距離[m]→ETA（バージョン互換）
try:
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)
        meters += min(d.get("length", 0.0) for d in ed.values())

walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# 4) GPKG に“合格ルート”だけ属性付きで保存
out_gpkg = "outputs/shizuoka_hazard.gpkg"
route_gdf = gpd.GeoDataFrame(
    {
        "meters":[meters],
        "eta_min_walk":[eta_walk_min],
        "eta_min_walk_safe":[eta_walk_safe],
        "eta_min_drive":[eta_drive_min],
    },
    geometry=[route_geom], crs="EPSG:4326",
)
route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
print("✅ 交差なし・保存完了")

# 5) （任意）Foliumに表示するならここで
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

を実行すると

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[12], line 26
     23 inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0
     25 # 近接（境界 ±2m 以内をNGにする例）
---> 26 near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
     28 if inside_len_m > 0 or near:
     29     raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
     30                        "別避難所 or 半径拡大で再探索してください。")

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/tools/overlay.py:273, in overlay(df1, df2, how, keep_geom_type, make_valid)
    270     raise ValueError(f"`how` was '{how}' but is expected to be in {allowed_hows}")
    272 if isinstance(df1, GeoSeries) or isinstance(df2, GeoSeries):
--> 273     raise NotImplementedError(
    274         "overlay currently only implemented for GeoDataFrames"
    275     )
    277 if not _check_crs(df1, df2):
    278     _crs_mismatch_warn(df1, df2, stacklevel=3)

NotImplementedError: overlay currently only implemented for GeoDataFrames

となる

原因を調べたら
route_m.buffer(2.0) は GeoSeries なので、
gpd.overlay()（＝GeoDataFrame同士専用）に渡すとエラーになる

バッファをGeoDataFrame化してから overlay するか、sjoinで近接判定に変える

今後も考慮してGPTに質問する

上位N候補の避難所からネットワーク距離最短を自動選定
10/20分アイソクローンの可視化
現在地が浸水想定内かの警告フラグ
GraphDB用にノード/リレーション出力（Home, Shelter, RoadNode/Edge, TsunamiPolygon, ROUTE_TO など）
も今後行うなら

今回変更するのは
overlay するか、sjoin
どちらが適切？

短く結論：
* 基本は sjoin 推奨（速い・メモリ軽い・属性付与しやすい）。
* 幾何そのものが必要なときだけ overlay（交差部分の長さ/面積を測る、クリップして可視化したい等）。
1. 道路を浸水域で“遮断”（真偽判定・フラグ付け）
* sjoin で十分（predicate=”intersects”）。
* 例：edges[“blocked”] = gpd.sjoin(edges, tsunami[[“geometry”]], predicate=”intersects”, how=”left”)[“index_right”].notna()
2. **近接NG（例：2m以内）**のチェック
* ルートを buffer(2m) → sjoin（軽い）。
* 例：near = not gpd.sjoin(route_m.buffer(2), flood_m, predicate=”intersects”).empty
3. 交差“長さ/面積”が要る（ルートが浸水域内に何 m 入ってる？）
* overlay(…, how=”intersection”) で交差ジオメトリを作り、length/area を合計。
* 今回の「inside_len_m」は overlay が最適。
4. AOIクリップや可視化用の切り出し
* 結果ジオメトリが要るので overlay。
5. 上位N避難所のネットワーク最短選定
* 候補抽出は直線距離（to_crs(UTM)）→ N件
* その後は NetworkX で経路長を比較。
* 候補フィルタ（浸水域内除外）は sjoin でOK。
6. 10/20分アイソクローン
* 等時間ポリゴン生成自体はネットワーク処理（eg. 到達可能ノードの凸包/αシェイプ等）。
* 仕上げの「浸水域で塗り分け」や表示用クリップは overlay。
7. GraphDB 連携（ノード/エッジに属性フラグ）
* IN_FLOOD, NEAR_FLOOD, BLOCKED 等のブール属性付与は sjoin が最適。
* 「交差長」など数値を持たせたい関係は overlay で測って格納。

いまのコードへの具体反映
* 継続：inside_len_m の算出 → overlay（長さが欲しいから）
* 修正済：近接チェック → sjoin（軽い）    

buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

この方針なら、今後の「上位N・アイソクローン・GraphDB」までスムーズに拡張できます。

ということなので
Sjoin　に変更する

# === 検証→ETA→保存（検証に通ったルートだけ保存・表示） ===
from shapely.geometry import LineString
import geopandas as gpd

# 1) ルートのLineString（まだ属性は付けない）
route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]  # (lon,lat)
route_geom = LineString(route_lonlat)
route_gdf  = gpd.GeoDataFrame(geometry=[route_geom], crs="EPSG:4326")

# 2) 堅牢チェック（交差ゼロ＆境界近接なし）
#   union_all が使えない環境向けにフォールバックも用意
try:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.geometry.union_all()], crs="EPSG:4326")
except Exception:
    flood_union = gpd.GeoDataFrame(geometry=[tsunami.unary_union], crs="EPSG:4326")

METRIC = "EPSG:32654"  # 静岡周辺（必要なら地域に合わせて変更）
route_m = route_gdf.to_crs(METRIC)
flood_m = flood_union.to_crs(METRIC)

# 交差長[m]
inter = gpd.overlay(route_m, flood_m, how="intersection")
inside_len_m = float(inter.geometry.length.sum()) if len(inter) else 0.0

# 近接（境界 ±2m 以内をNGにする例）
# near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty


if inside_len_m > 0 or near:
    raise RuntimeError(f"ルート不適合: 交差長 {inside_len_m:.2f} m / 近接={near}. "
                       "別避難所 or 半径拡大で再探索してください。")

# 3) 距離[m]→ETA（バージョン互換）
try:
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)
        meters += min(d.get("length", 0.0) for d in ed.values())

walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
eta_walk_safe = eta_walk_min * 1.3  # 安全係数（例：+30%）

# 4) GPKG に“合格ルート”だけ属性付きで保存
out_gpkg = "outputs/shizuoka_hazard.gpkg"
route_gdf = gpd.GeoDataFrame(
    {
        "meters":[meters],
        "eta_min_walk":[eta_walk_min],
        "eta_min_walk_safe":[eta_walk_safe],
        "eta_min_drive":[eta_drive_min],
    },
    geometry=[route_geom], crs="EPSG:4326",
)
route_gdf.to_file(out_gpkg, layer="route_latest", driver="GPKG")
print("✅ 交差なし・保存完了")

# 5) （任意）Foliumに表示するならここで
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分 (安全側 {eta_walk_safe:.1f} 分)"
).add_to(m)

というように

# 近接（境界 ±2m 以内をNGにする例）
# near = not gpd.overlay(route_m.buffer(2.0), flood_m, how="intersection").empty
を
buf_gdf = gpd.GeoDataFrame(geometry=route_m.buffer(2.0), crs=METRIC)
near = not gpd.sjoin(buf_gdf, flood_m, how="inner", predicate="intersects").empty

にすることで

✅ 交差なし・保存完了
<folium.vector_layers.PolyLine at 0x31c488b10>

となった

ここまでで「安全確認→距離/ETA→保存→表示」の最小ループは完成

を行っていく

ルート距離・ETA計算 → GPKG保存 → 追加の堅牢化

1) ルート長とETA（徒歩/車）を出す

# ルートの延長[m]は G2 の edge 属性 'length' を合計（OSMnxが付けてくれてます）
meters = sum(ox.utils_graph.get_route_edge_attributes(G2, route, 'length'))
walk_kmh = 4.8   # 歩行速度（必要に応じて変更）
drive_kmh = 20.0 # 災害時の車想定の控えめ値（使うなら）

eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
print(f"距離: {meters:.0f} m, 徒歩ETA: {eta_walk_min:.1f} 分, 車ETA(仮): {eta_drive_min:.1f} 分")

を実行すると

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 2
      1 # ルートの延長[m]は G2 の edge 属性 'length' を合計（OSMnxが付けてくれてます）
----> 2 meters = sum(ox.utils_graph.get_route_edge_attributes(G2, route, 'length'))
      3 walk_kmh = 4.8   # 歩行速度（必要に応じて変更）
      4 drive_kmh = 20.0 # 災害時の車想定の控えめ値（使うなら）

AttributeError: module 'osmnx' has no attribute 'utils_graph'

これは OSMnx のAPI変更が原因です。utils_graph.get_route_edge_attributes は環境によって存在しない

# ▼ 距離[m]の取り方をバージョン互換で
try:
    # 新しめのOSMnx（推奨）: ルートをGDF化して length を合計
    route_edges = ox.routing.route_to_gdf(G2, route, fill_edge_geometry=True)
    meters = float(route_edges["length"].sum())
except Exception:
    # フォールバック: NetworkXからエッジ属性 length を合計
    meters = 0.0
    for u, v in zip(route[:-1], route[1:]):
        ed = G2.get_edge_data(u, v)  # MultiDiGraph: {key: attr_dict}
        if isinstance(ed, dict) and len(ed):
            meters += min(d.get("length", 0.0) for d in ed.values())

# ETA（必要なら速度を調整）
walk_kmh = 4.8
drive_kmh = 20.0
eta_walk_min  = meters / (walk_kmh  * 1000 / 60)
eta_drive_min = meters / (drive_kmh * 1000 / 60)
print(f"距離: {meters:.0f} m, 徒歩ETA: {eta_walk_min:.1f} 分, 車ETA(仮): {eta_drive_min:.1f} 分")

# Folium のツールチップにも表示
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(
    route_latlon, weight=6, opacity=0.9,
    tooltip=f"距離 {meters:.0f} m / 徒歩 {eta_walk_min:.1f} 分"
).add_to(m)

に変更

これで

距離: 3017 m, 徒歩ETA: 37.7 分, 車ETA(仮): 9.1 分
<folium.vector_layers.PolyLine at 0x33d5e6650>

となった

約3.0 km / 徒歩 ~38分 / 車 ~9分（仮） は現実的な値

とりあえずやることは
1) ルート長とETA（徒歩/車）を出す
2) ルートを GPKG に属性付きで保存
3) “ネットワーク的に最も近い”安全避難所を選ぶ（推奨）
いまは直線距離で1件選んでいます。候補を直線距離Top Nに絞ってから実際の経路長で最短を選ぶと堅牢
4) 追加の安全策（おすすめ）
* 自宅が浸水域内かを先に判定：内なら「即時・垂直避難」をUIに強調し、経路探索は参考扱いに。
5) 動作チェックのミニ・リスト
* evac.html を開いて、赤いルートが青の浸水域と交差していないか。
* shelters_safe に目的地が含まれているか（含まれなければCSV側の位置精度を見直し）。
* QGISで roads_safe と tsunami_inundation を重ね、視覚的にもOKか。

まずは
1) ルートをGPKGに属性付きで保存

from shapely.geometry import LineString
import geopandas as gpd

route_lonlat = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]
route_gdf = gpd.GeoDataFrame(
    {"meters":[meters], "eta_min_walk":[eta_walk_min], "eta_min_drive":[eta_drive_min]},
    geometry=[LineString(route_lonlat)], crs="EPSG:4326"
)
route_gdf.to_file("outputs/shizuoka_hazard.gpkg", layer="route_latest", driver="GPKG")

2) ルートが浸水域を跨いでないか最終チェック

import geopandas as gpd
ok = not route_gdf.intersects(tsunami.unary_union).iloc[0]
print("交差なし✅" if ok else "⚠️ 交差あり（パラメータ見直し）")

実行結果は

/var/folders/db/tkl5w9dd3kn3h53ctyl4s6180000gn/T/ipykernel_49497/3699053437.py:2: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  ok = not route_gdf.intersects(tsunami.unary_union).iloc[0]

交差なし✅

となった

DeprecationWarning は unary_union が将来的に非推奨という通知だけなので、動作上は問題ありません。念のため推奨APIに差し替える

import geopandas as gpd
# ok = not route_gdf.intersects(tsunami.unary_union).iloc[0]
# print("交差なし✅" if ok else "⚠️ 交差あり（パラメータ見直し）")

# 旧: tsunami.unary_union  ← 警告が出る
# 新:
flood_union = tsunami.geometry.union_all()
ok = not route_gdf.intersects(flood_union).iloc[0]
print("交差なし✅" if ok else "⚠️ 交差あり（パラメータ見直し）")

これで警告は消える

これを修正したらGPTから
もう一段だけ堅牢に（任意）の項目が提案された
この部分は
3) 安全バッファ付きETA（渋滞・混雑を考慮）
4) Foliumツールチップに表示
の前に行う処理？
と聞いてみる

結論：はい、“堅牢チェック”は 3)安全バッファ付きETA・4)ツールチップ表示の前にやるのがベストです。
理由：交差や境界ニアミスが見つかったらその経路は捨てて再探索するので、ETAや表示は検証済みルートに対してだけ出すのが安全

おすすめの処理順
1. ルート算出（orig→dest）
2. 堅牢チェック
* union_all() で津波ポリゴン結合
* 交差の有無 → あれば再探索（別避難所 or 半径拡大）
* 交差長[m]（=0が理想）
* 近接判定（例：ルート1–5mバッファが浸水域に当たるか）
3. （ここで初めて）距離[m]→ETA計算（徒歩/車＋安全係数）
4. Folium描画（ツールチップに距離/ETAを表示）
5. GPKGへ route_latest などを保存

チェック（交差・近接の堅牢判定）→ OKなら保存の順がベスト
不適切なルートをGPKGに残さないため。保存は“合格”ルートだけにしましょう。
とのこと

道路→浸水域回避→最寄り安全避難所への経路の作成

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString
import osmnx as ox
import networkx as nx
import folium
import os

GPKG = "outputs/shizuoka_hazard.gpkg"

# 1) 既存レイヤ読み込み（両方とも 4326）
coast   = gpd.read_file(GPKG, layer="coastline")
tsunami = gpd.read_file(GPKG, layer="tsunami_inundation")

# ----- パラメータ -----
# 自宅座標（lat, lon）←必ず置き換え
home_wgs = (34.728, 137.978)
# OSM道路の取得半径[m]（到達不能なら大きく：15000, 20000…）
GRAPH_RADIUS_M = 12000
# 避難所CSV（cp932想定）
SHELTER_CSV = "data/20250801_指定緊急避難場所一覧.CSV"
# ----------------------

# 2) 避難所CSV→点化（列名は自動検出）
df = pd.read_csv(SHELTER_CSV, encoding="cp932")
lat_candidates = 
lon_candidates = 
if not lat_candidates or not lon_candidates:
    raise ValueError("CSVに緯度/経度列が見つかりません。列名を教えてください。")
lat_col, lon_col = lat_candidates[0], lon_candidates[0]
df[lat_col] = pd.to_numeric(df[lat_col], errors="coerce")
df[lon_col] = pd.to_numeric(df[lon_col], errors="coerce")
df = df.dropna(subset=[lat_col, lon_col])
shelters = gpd.GeoDataFrame(df.copy(),
    geometry=gpd.points_from_xy(df[lon_col], df[lat_col], crs="EPSG:4326")
)

# 3) 安全な避難所 = 浸水域と非交差
#   （unionはメモリ負荷になるので、空間結合で判定）
safe_join = gpd.sjoin(shelters, tsunami[["geometry"]], how="left", predicate="intersects")
shelters["safe"] = safe_join["index_right"].isna()
safe_shelters = shelters[shelters["safe"]].copy()

# 4) OSM道路取得（徒歩網）
G = ox.graph_from_point(home_wgs, dist=GRAPH_RADIUS_M, network_type="walk", simplify=True, retain_all=False)
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True)
if edges.crs is None: edges.set_crs(4326, inplace=True)
elif edges.crs.to_epsg() != 4326: edges = edges.to_crs(4326)

# 5) 浸水域と交差する道路を除去（空間結合で高速化）
tag = gpd.sjoin(edges[["u","v","key","geometry"]], tsunami[["geometry"]],
                how="left", predicate="intersects")
tag["blocked"] = tag["index_right"].notna()
edges_keep = edges.loc[~tag["blocked"].values].copy()

# 6) サブグラフ再構築（通行可のみ）
used_nodes = pd.unique(edges_keep[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()
G2 = ox.utils_graph.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)

# 7) 出発・目的地ノード（最寄り“安全”避難所）
home_pt = Point(home_wgs[1], home_wgs[0])
# 最寄り安全避難所（直線距離）…なければ全体から選ぶ
TARGET = safe_shelters if len(safe_shelters) else shelters
METRIC = "EPSG:32654"  # 静岡周辺 UTM
hp = gpd.GeoSeries([home_pt], crs=4326).to_crs(METRIC).iloc[0]
tp = TARGET.to_crs(METRIC).copy()
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = TARGET.loc[tp["__dist_m__"].idxmin(), "geometry"]

orig = ox.distance.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)

# 8) 最短経路（通行可能グラフ）
try:
    route = nx.shortest_path(G2, orig, dest, weight="length")
except nx.NetworkXNoPath:
    raise RuntimeError("到達経路なし。GRAPH_RADIUS_Mを増やすか、別の避難所で再試行を。")

# 9) GPKGへ追記
#    既にある場合は上書きしたいなら driverオプションを切り替える（ここは追記）
shelters.to_file(GPKG, layer="shelters_all", driver="GPKG")
safe_shelters.to_file(GPKG, layer="shelters_safe", driver="GPKG")
edges.to_file(GPKG, layer="roads_all", driver="GPKG")
edges_keep.to_file(GPKG, layer="roads_safe", driver="GPKG")

# 10) Folium地図
m = folium.Map(location=[home_wgs[0], home_wgs[1]], zoom_start=12, control_scale=True)
folium.GeoJson(tsunami.__geo_interface__, name="津波浸水想定(A40)").add_to(m)
folium.GeoJson(edges_keep.__geo_interface__, name="道路(浸水回避)").add_to(m)

# 経路を線で描画
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(route_latlon, weight=6, opacity=0.9, tooltip="避難経路").add_to(m)

# 出発・目的地
folium.Marker([home_wgs[0], home_wgs[1]], popup="現在地", icon=folium.Icon(color="green")).add_to(m)
folium.Marker([dest_pt.y, dest_pt.x], popup="避難所", icon=folium.Icon(color="blue")).add_to(m)

# 避難所レイヤ（安全のみ強調）
folium.GeoJson(safe_shelters.__geo_interface__, name="安全な避難所",
               marker=folium.CircleMarker(radius=4)).add_to(m)
folium.LayerControl().add_to(m)

# 保存（任意）
m.save("outputs/evac.html")
print("✓ 完了: GPKG追記 & outputs/evac.html を作成")

実行すると

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[1], line 51
     48 elif edges.crs.to_epsg() != 4326: edges = edges.to_crs(4326)
     50 # 5) 浸水域と交差する道路を除去（空間結合で高速化）
---> 51 tag = gpd.sjoin(edges[["u","v","key","geometry"]], tsunami[["geometry"]],
     52                 how="left", predicate="intersects")
     53 tag["blocked"] = tag["index_right"].notna()
     54 edges_keep = edges.loc[~tag["blocked"].values].copy()

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/geodataframe.py:1896, in GeoDataFrame.__getitem__(self, key)
   1890 def __getitem__(self, key):
   1891     """
   1892     If the result is a column containing only 'geometry', return a
   1893     GeoSeries. If it's a DataFrame with any columns of GeometryDtype,
   1894     return a GeoDataFrame.
   1895     """
-> 1896     result = super().__getitem__(key)
   1897     # Custom logic to avoid waiting for pandas GH51895
   1898     # result is not geometry dtype for multi-indexes
   1899     if (
   1900         pd.api.types.is_scalar(key)
   1901         and key == ""
   (...)   1904         and not is_geometry_type(result)
   1905     ):

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pandas/core/frame.py:4108, in DataFrame.__getitem__(self, key)
   4106     if is_iterator(key):
   4107         key = list(key)
-> 4108     indexer = self.columns._get_indexer_strict(key, "columns")[1]
   4110 # take() does not accept boolean indexers
   4111 if getattr(indexer, "dtype", None) == bool:

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pandas/core/indexes/base.py:6200, in Index._get_indexer_strict(self, key, axis_name)
   6197 else:
   6198     keyarr, indexer, new_indexer = self._reindex_non_unique(keyarr)
-> 6200 self._raise_if_missing(keyarr, indexer, axis_name)
   6202 keyarr = self.take(indexer)
   6203 if isinstance(key, Index):
   6204     # GH 42790 - Preserve name from an Index

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pandas/core/indexes/base.py:6252, in Index._raise_if_missing(self, key, indexer, axis_name)
   6249     raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   6251 not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())
-> 6252 raise KeyError(f"{not_found} not in index")

KeyError: "['u', 'v', 'key'] not in index"

原因はコレ：graph_to_gdfs の edges は MultiIndex（u,v,key が“列ではなくインデックス”）なので、edges[[“u”,”v”,”key”,”geometry”]] が KeyError になります。
→ reset_index() で列に起こす＋sjoin の重複行を集計してからフィルタしてください

ということで

# --- 4) OSM道路取得（徒歩網） ---
G = ox.graph_from_point(home_wgs, dist=GRAPH_RADIUS_M, network_type="walk", simplify=True, retain_all=False)
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True)

# u,v,key を“列”にする（重要）
edges = edges.reset_index()

# CRSを4326に
if edges.crs is None:
    edges = edges.set_crs(4326)
elif edges.crs.to_epsg() != 4326:
    edges = edges.to_crs(4326)

# --- 5) 浸水域と交差する道路を除去（sjoin→集計で一意化） ---
# 属性を絞って結合
left = edges[["u","v","key","geometry"]].copy()
right = tsunami[["geometry"]]

tag = gpd.sjoin(left, right, how="left", predicate="intersects")

# エッジごとに「1回でも当たったら blocked=True」
blocked = (
    tag.assign(hit=tag["index_right"].notna())
       .groupby(["u","v","key"], as_index=False)["hit"].max()
       .rename(columns={"hit":"blocked"})
)

# 元の edges に blocked を付与（ない行は False）
edges2 = edges.merge(blocked, on=["u","v","key"], how="left")
edges2["blocked"] = edges2["blocked"].fillna(False)

edges_keep = edges2.loc[~edges2["blocked"]].copy()

# --- 6) サブグラフ再構築 ---
used_nodes = pd.unique(edges_keep[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()
G2 = ox.utils_graph.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)

へコード変更

しかしまだ

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 80
     78 used_nodes = pd.unique(edges_keep[["u","v"]].values.ravel())
     79 nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()
---> 80 G2 = ox.utils_graph.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)
     83 # 7) 出発・目的地ノード（最寄り“安全”避難所）
     84 home_pt = Point(home_wgs[1], home_wgs[0])

AttributeError: module 'osmnx' has no attribute 'utils_graph'

となる

これは
OSMnx のバージョン差です。
utils_graph.graph_from_gdfs は古い呼び方で、いまは ox.graph_from_gdfs（トップレベル）を使います
とのことなので

# --- 6) サブグラフ再構築 ---
used_nodes = pd.unique(edges_keep[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()
G2 = ox.utils_graph.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)

を

# --- 6) サブグラフ再構築（通行可のみ） ---
used_nodes = pd.unique(edges_keep[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()

# ▶ ここを置き換え
try:
    # 新しめの OSMnx（1.1+ 〜 2.x）
    G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))
except AttributeError:
    # 旧版 OSMnx 向けフォールバック
    from osmnx import utils_graph as ug
    G2 = ug.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))

に変更

最寄りノード取得もバージョン差があるので互換にしておく

home_pt = Point(home_wgs[1], home_wgs[0])

TARGET = safe_shelters if len(safe_shelters) else shelters
METRIC = "EPSG:32654"  
hp = gpd.GeoSeries([home_pt], crs=4326).to_crs(METRIC).iloc[0]
tp = TARGET.to_crs(METRIC).copy()
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = TARGET.loc[tp["__dist_m__"].idxmin(), "geometry"]

orig = ox.distance.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)

を

# 7) 出発・目的地ノード
home_pt = Point(home_wgs[1], home_wgs[0])
TARGET = safe_shelters if len(safe_shelters) else shelters
METRIC = "EPSG:32654"
hp = gpd.GeoSeries([home_pt], crs=4326).to_crs(METRIC).iloc[0]
tp = TARGET.to_crs(METRIC).copy()
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = TARGET.loc[tp["__dist_m__"].idxmin(), "geometry"]

# ▶ 互換呼び分け
try:
    orig = ox.distance.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)
except AttributeError:
    # 新API名（環境によってはこちらが有効）
    orig = ox.nearest_nodes(G2, X=home_pt.x, Y=home_pt.y)
    dest = ox.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)

に変更

OSMnx のバージョンは import osmnx as ox; print(ox.__version__) で確認できます。
もし整えたいなら

 pip install -U osmnx

で最新版に寄せるのもアリ
とのこと

--------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[3], line 87
     84 # ▶ ここを置き換え
     85 try:
     86     # 新しめの OSMnx（1.1+ 〜 2.x）
---> 87     G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))
     88 except AttributeError:
     89     # 旧版 OSMnx 向けフォールバック
     90     from osmnx import utils_graph as ug

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/osmnx/convert.py:302, in graph_from_gdfs(gdf_nodes, gdf_edges, graph_attrs)
    265 def graph_from_gdfs(
    266     gdf_nodes: gpd.GeoDataFrame,
    267     gdf_edges: gpd.GeoDataFrame,
    268     *,
    269     graph_attrs: dict[str, Any] | None = None,
    270 ) -> nx.MultiDiGraph:
    271     """
    272     Convert node and edge GeoDataFrames to a MultiDiGraph.
    273 
   (...)    300         The converted MultiDiGraph.
    301     """
--> 302     _validate_node_edge_gdfs(gdf_nodes, gdf_edges)
    304     # drop geometry column from gdf_nodes (since we use x and y for geometry
    305     # information), but warn the user if the geometry values differ from the
    306     # coordinates in the x and y columns. this results in a df instead of gdf.
    307     if gdf_nodes.active_geometry_name is None:  # pragma: no cover

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/osmnx/convert.py:241, in _validate_node_edge_gdfs(gdf_nodes, gdf_edges)
    239 edges_index_levels = 3
    240 check1 = gdf_edges.index.nlevels == edges_index_levels
--> 241 uv = set(gdf_edges.index.get_level_values(0)) | set(gdf_edges.index.get_level_values(1))
    242 check2 = uv.issubset(set(gdf_nodes.index))
    243 if not (check1 and check2):  # pragma: no cover

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pandas/core/indexes/base.py:2102, in Index._get_level_values(self, level)
   2066 def _get_level_values(self, level) -> Index:
   2067     """
   2068     Return an Index of values for requested level.
   2069 
   (...)   2100     Index(['a', 'b', 'c'], dtype='object')
   2101     """
-> 2102     self._validate_index_level(level)
   2103     return self

File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pandas/core/indexes/base.py:2008, in Index._validate_index_level(self, level)
   2003         raise IndexError(
   2004             "Too many levels: Index has only 1 level, "
   2005             f"{level} is not a valid level number"
   2006         )
   2007     if level > 0:
-> 2008         raise IndexError(
   2009             f"Too many levels: Index has only 1 level, not {level + 1}"
   2010         )
   2011 elif level != self.name:
   2012     raise KeyError(
   2013         f"Requested level ({level}) does not match index name ({self.name})"
   2014     )

IndexError: Too many levels: Index has only 1 level, not 2


原因は、ox.graph_from_gdfs が受け取る edges のインデックス形式です。
この関数は edges のインデックスが MultiIndex（u, v, key） であることを前提にしています。先で reset_index() したので RangeIndex になっており、検証でコケています

# --- 5) 浸水域と交差する道路を除去（sjoin→集計で一意化） ---
left = edges[["u","v","key","geometry"]].copy()   # edges は reset_index() 済み想定
right = tsunami[["geometry"]]

tag = gpd.sjoin(left, right, how="left", predicate="intersects")

blocked = (
    tag.assign(hit=tag["index_right"].notna())
       .groupby(["u","v","key"], as_index=False)["hit"].max()
       .rename(columns={"hit":"blocked"})
)

edges2 = edges.merge(blocked, on=["u","v","key"], how="left")
edges2["blocked"] = edges2["blocked"].fillna(False)

edges_keep = edges2.loc[~edges2["blocked"]].copy()

# ★ ここが重要：graph_from_gdfs に渡す前に MultiIndex に戻す
#   （key が float になっていたら int にしておくと安心）
if edges_keep["key"].isna().any():
    edges_keep["key"] = edges_keep["key"].fillna(0)

edges_keep["key"] = edges_keep["key"].astype(int)
edges_keep = edges_keep.set_index(["u","v","key"]).sort_index()

# --- 6) サブグラフ再構築（通行可のみ） ---
used_nodes = pd.unique(edges_keep.reset_index()[["u","v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_nodes)].copy()

# nodes_keep の geometry 列が欠けていたら補完
if "geometry" not in nodes_keep or nodes_keep["geometry"].isna().any():
    nodes_keep = nodes_keep.set_geometry(
        gpd.points_from_xy(nodes_keep["x"], nodes_keep["y"]), crs=edges_keep.crs
    )

# OSMnx のバージョン差に対応
try:
    G2 = ox.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))
except AttributeError:
    from osmnx import utils_graph as ug
    G2 = ug.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=getattr(G, "graph", None))

これで

✓ 完了: GPKG追記 & outputs/evac.html を作成

となって

print("edges_keep index levels:", getattr(edges_keep.index, "nlevels", 1))  # → 3 が正解
print(edges_keep.index.names)  # → ['u','v','key'] が正解

の結果は

edges_keep index levels: 3
['u', 'v', 'key']

となる

次に
“差分だけ”足せば、ルート距離・ETA計算 → GPKG保存 → 追加の堅牢化まで整う