次に作る機能について

次に作る機能について

これで

使用する用途のコンセプト また震源地から津波が到達するまでの時間と津波の高さを計算し 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(当日実用版)

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

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

* `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 での確認

再度 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(目的地の点ジオメトリ)を作っていないのが原因

修正案は2つ
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のメタ情報を混ぜる)するのが安全です。

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→保存→表示」の最小ループは完成

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

を行っていく

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

ルート距離・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保存 → 追加の堅牢化まで整う

国土数値情報「海岸線」データ(全国版 Shapefile)

国土数値情報「海岸線」データ(全国版 Shapefile)のURL
を調べたが全国版はない

都道府県別に Shapefile 形式や GML 形式で取得可能な国土数値情報「海岸線」データのダウンロードページがある

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C23.html?utm_source=chatgpt.com
が海岸線データのリンクだけど
最新のデータは2006年度(平成18年度)版

他を調べてみると

の情報があった
https://nlftp.mlit.go.jp/ksj/index.html

これらを元にGPTで調べてみる

📊 スクショにあるデータの役割
* 津波浸水想定区域データ
→ 津波が来た場合の浸水範囲(自治体がハザードマップで使う)。
* 洪水浸水想定区域データ
→ 河川ごと、メッシュごとの洪水想定範囲。
* 高潮・雨水出水・多段階浸水データ
→ 津波以外の沿岸・豪雨災害関連のリスク領域。
* 避難所・避難道路系(行政区、緊急輸送道路など)
→ 避難経路や防災拠点の把握に必須。

🌊 あなたの用途別の適材適所
1. 「海岸線」データ(国土数値情報 2006年版)
* 海と陸の境界を決める基盤。
* 最新ではないが行政基準値。
* → 基準線として使う。
2. スクショの「津波浸水想定区域データ」
* シミュレーション済みの「どこまで水が来るか」。
* → 津波到達後の浸水エリア算出に必須。
3. スクショの「洪水・高潮データ」
* 津波以外のリスクも扱える。
* → 家族の防災マップを広くカバー。
4. 避難所・緊急輸送道路データ
* 実際の「どこに逃げるか」「どう行くか」。
* → GraphDBの避難所ノード・道路ノードに登録。

✅ 推奨データ利用の組み合わせ
* ベース地図:
* 国土数値情報「海岸線」 + OSM Coastlines(更新補完)
* リスク領域:
* 津波浸水想定区域データ(スクショにある)
* 洪水・高潮想定区域データ
* 地形:
* DEM(国土地理院/GEBCO)
* 避難系:
* 避難所・緊急輸送道路・行政区データ
* GraphDB:
* ノード:震源・海岸線・浸水想定区域・避難所・道路・家族
* リレーション:到達時間・避難経路・所要時間・備蓄品

🔑 結論:
スクショにある「津波浸水想定区域データ」や「洪水浸水想定区域データ」は 必ず併用すべきデータ です。
単なる「海岸線」だけでは津波到達範囲はわからないので、
* 海岸線 = 境界線
* 浸水想定区域 = 被害範囲
として 両方組み合わせるのがベスト

まずは居住区域のみで作成して動くようなら適用範囲を拡大していく

pip install geopandas shapely rtree pyproj osmnx networkx fiona

で必要なパッケージのインストール

📦 パッケージ解説
1. GeoPandas
* 役割: pandas の DataFrame に地理情報(ポリゴンやポイントなどのジオメトリ)を扱えるようにした拡張版。
* 特徴:
* shapefile や GeoJSON など GIS ファイルを直接読み書きできる。
* 空間結合(spatial join)、バッファリング(一定距離の範囲作成)、重なり判定などが簡単にできる。
* 用途例: 避難所の位置と自宅の位置から最寄り避難所を検索。

2. Shapely
* 役割: ジオメトリ(点、線、面)の作成・操作・解析を行うライブラリ。
* 特徴:
* ポリゴンの交差判定やバッファリング、距離計算が可能。
* GeoPandas 内部でも Shapely が使われている。
* 用途例: 津波の浸水範囲ポリゴンと建物の位置の交差判定。

3. Rtree
* 役割: 空間インデックスを構築して、高速に空間検索を行うためのライブラリ。
* 特徴:
* 例えば「ある点から半径 1km 以内の避難所を検索」といった処理が高速化できる。
* GeoPandas で .sindex を使うと裏で Rtree が動作。
* 用途例: 半径○km以内の避難場所を効率的に抽出。

4. PyProj
* 役割: 座標系(緯度経度、UTMなど)の変換を行うライブラリ。
* 特徴:
* EPSG コードを指定して座標変換ができる。
* GeoPandas で .to_crs() を使うと内部で PyProj が使われる。
* 用途例: 世界測地系(WGS84, EPSG:4326)から日本の平面直角座標系(EPSG:6677など)への変換。

5. OSMnx
* 役割: OpenStreetMap (OSM) データを Python で取得・解析するライブラリ。
* 特徴:
* 道路ネットワークや建物形状を簡単にダウンロード可能。
* networkx 形式で道路ネットワークを解析できる。
* 用途例: 自宅から避難所までの徒歩ルート検索。

6. NetworkX
* 役割: グラフ(ネットワーク)解析ライブラリ。
* 特徴:
* 最短経路探索、中心性、連結成分などのネットワーク理論を使った解析が可能。
* OSMnx が取得した道路データを解析するのに必須。
* 用途例: 最短避難経路の算出(Dijkstra アルゴリズムなど)。

7. Fiona
* 役割: 各種 GIS ファイルの入出力を行うためのライブラリ。
* 特徴:
* shapefile や GeoJSON を読み込むときに GeoPandas が内部で使用。
* GDAL(GIS の基盤ライブラリ)の Python バインディング的な存在。
* 用途例: 国土地理院の災害情報シェープファイルの読み込み。

🔗 関係性まとめ
* GeoPandas = GIS用 DataFrame(中心的存在)
* 内部で Shapely(ジオメトリ処理)、PyProj(座標変換)、Fiona(ファイルIO)、Rtree(高速検索)を利用。
* OSMnx = OpenStreetMap から道路や建物を取得 → NetworkX で解析。

✅ あなたの防災マップ用途との対応
* 津波浸水範囲と建物位置の判定 → Shapely / GeoPandas
* 避難所の検索(近い順) → Rtree / GeoPandas
* ルート案内 → OSMnx + NetworkX
* 座標変換(地図投影法の切り替え) → PyProj
* データの入出力(shapefile, GeoJSONなど) → Fiona

👉 この7つをまとめてインストールしておくと、ほぼ PythonでGIS + 避難経路システムを自作できる環境が揃います。

まずは全国ではなく居住区域の静岡県のみで作成して、テストして動くようなら適用範囲を拡大していく

まずは静岡県のみで最小セットを作って、表示・経路探索・レイヤ結合が全部通るか確認

「海岸線(KSJ)+津波浸水想定(KSJ)+道路(OSM)+避難所(自治体オープンデータがあれば)」をGeoPackage 1ファイルにまとめて、テスト用に「自宅→最寄り避難所の最短経路(浸水域は回避)」まで動かす

使うデータ(静岡県のみ)
* 海岸線:国土数値情報 C23 海岸線(Shapefile/GML)
※ページ:国土数値情報ダウンロード → C23(都道府県別を選択して「静岡県」をDL)
* 津波浸水想定区域:国土数値情報(都道府県別配布の「津波浸水想定」)
※静岡県のファイルをDL
* 道路:OSMnxでOSMから静岡県全域の道路グラフを取得(最新更新が反映されやすい)
* 避難所:静岡県/市町のオープンデータ(CSV/GeoJSON、なければ後で追加)
ポイント
* 最初は海岸線+津波浸水+道路の3つでOK。
* 経路探索やバッファ計算など**距離単位(m)**が必要になるので、計算時は **WGS84/UTM Zone 54N(EPSG:32654)**に投影して処理します(静岡は概ね54N)。

このために必要なファイルのダウンロード
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C23.html?utm_source=chatgpt.com
から海岸線データのダウンロード
今回は静岡県のみ

次に
津波浸水想定区域データのダウンロード
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A40-2024.html
からダウンロード
これもまずは静岡県のみ
ただし更新が2016年になっている
他の地域だと2024年まで更新されている

難所:静岡県/市町のオープンデータ(CSV/GeoJSON、なければ後で追加)
をやろうとしたけど
居住の市にCSVファイルがない…

https://www.city.fukuroi.shizuoka.jp/opendata/opendata/data/1490706118930.html

ふじのくにオープンデータカタログ
にあるらしいがリンクが違っているので検索する

https://opendata.pref.shizuoka.jp/dataset/

ふじのくにオープンデータカタログ

https://opendata.pref.shizuoka.jp/dataset/search/?s%5Barea_id%5D=36
に袋井市があった

指定緊急避難場所一覧
https://opendata.pref.shizuoka.jp/dataset/10246.html
最新版は
20250801_指定緊急避難場所一覧 (CSV 17.8KB)
なので
これをダウンロードする

次にプロジェクトフォルダの作成

mkdir FamilySafetyMap

名前を決めるのは苦手なので gemini を使った

源地から津波が到達するまでの時間と津波の高さを計算し GIS情報を使い避難するエリアを算出しルート案内できるようにしたい 家族の避難経路・地域防災マップをGraphDBで整理 避難所や重要連絡先も紐付け 使い方例: 「地震が来たら最寄りの避難所はどこ?」 「非常食の賞味期限が近いものは?」 など この用途で使うプロジェクトファイル名を決めたい。
使用するのは家族がメイン

実用・機能重視FamilySafetyMap家族の安全と地図(避難経路、GIS情報)という機能が明確。シンプルで覚えやすい。

が答え

cd FamilySafetyMap 
mkdir -p data/coastline_shizuoka/
mkdir -p data/tsunami_shizuoka/  

ファイルを設置しようとしたけど
どれを入れるかわからなくなったので

どっちが海岸線?どっちが津波浸水想定?
* C23-06_22_GML → 中に *_Coastline.* があり、これは 海岸線(C23)
例:C23-06_22-g_Coastline.shp ← これを 海岸線 として使う
* A40-16_22_GML → A40-16_22.* は 津波浸水想定区域(A40) のシェープ
つまり:
海岸線 = C23、津波浸水 = A40 です(認識を入れ替えてください)。

置くべきファイル(Shapefileの必須構成)
各フォルダに、少なくとも下記が同じ名前で揃っていればOK:
* *.shp(本体)
* *.shx(インデックス)
* *.dbf(属性)
* *.prj(座標系:あると良い・ほぼ入ってます)
*.xml(メタデータ)はあって良い(無くても読み込みは可能)。

あなたの手元の例
* 海岸線フォルダ(C23-06_22_GML)
* C23-06_22-g_Coastline.shp / .shx / .dbf / .prj … 使用するのはこれ
* C23-06_22-g.xml, KS-META-*.xml … メタデータ(任意)
* 津波浸水フォルダ(A40-16_22_GML)
* A40-16_22.shp / .shx / .dbf / .prj … これを使う
* A40-16_22.xml, KS-META-*.xml … メタデータ(任意)

スクリプトに設定するパス(例)

COAST_SHP   = "data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp"
TSUNAMI_SHP = "data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp"

あと
20250801_指定緊急避難場所一覧.CSV
の文字コードを調べる

file -I 20250801_指定緊急避難場所一覧.CSV
20250801_指定緊急避難場所一覧.CSV: text/csv; charset=unknown-8bit

となるので

brew install nkf uchardet

としたが

Error: Cannot install under Rosetta 2 in ARM default prefix (/opt/homebrew)!
To rerun under ARM use:
    arch -arm64 brew install ...
To install under x86_64, install Homebrew into /usr/local.

となる
久々で忘れていたけど
M1Mac の場合は

arch -arm64 brew install nkf uchardet

となる

uchardet 20250801_指定緊急避難場所一覧.CSV
SHIFT_JIS
 nkf -g 20250801_指定緊急避難場所一覧.CSV
CP932

このため

import pandas as pd
df = pd.read_csv("20250801_指定緊急避難場所一覧.CSV", encoding="cp932")

の文字コードはそのままにする

これらを元に

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

# =========================
# 0) 入力(あなたの実ファイルに合わせて修正)
# =========================
COAST_SHP   = "data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp"  # C23=海岸線
TSUNAMI_SHP = "data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp"                # A40=津波浸水想定
SHELTER_CSV = "data/20250801_指定緊急避難場所一覧.CSV"                            # cp932想定

OUT_GPKG = "outputs/shizuoka_hazard.gpkg"
os.makedirs(os.path.dirname(OUT_GPKG), exist_ok=True)

# 自宅の位置(WGS84)※必ずあなたの座標に置換(lat, lon)
home_wgs = (34.728, 137.978)  # 例: 袋井市周辺

# OSM道路グラフ取得範囲(自宅中心 半径[m])
GRAPH_RADIUS_M = 10000  # 10km程度で様子見。広げたければ 15000~20000


# =========================
# 1) データ読み込み・CRS整備
# =========================
# 海岸線(C23)/ 津波浸水想定(A40)をWGS84に
coast = gpd.read_file(COAST_SHP)
if coast.crs is None:
    raise ValueError("海岸線ShapefileにCRSがありません(.prjが同フォルダにあるか確認)")
coast = coast.to_crs(4326)

tsunami = gpd.read_file(TSUNAMI_SHP)
if tsunami.crs is None:
    raise ValueError("津波浸水想定ShapefileにCRSがありません(.prjが同フォルダにあるか確認)")
tsunami = tsunami.to_crs(4326)
tsunami_union = tsunami.unary_union  # 高速化用

# 避難所CSV(cp932想定)→ 可能な列名を自動検出して点に
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 = lat_candidates[0]
lon_col = 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")
)

# 浸水域外の避難所のみを「安全」とする
safe_shelters = shelters.loc[~shelters.intersects(tsunami_union)].copy()


# =========================
# 2) 道路グラフ(OSMnx)
# =========================
# 自宅中心の徒歩グラフ(避難行動向け)。範囲広すぎると重いので半径で制限
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)
# 以降WGS84で統一
if edges.crs is None:
    edges.set_crs(4326, inplace=True)
elif edges.crs.to_epsg() != 4326:
    edges = edges.to_crs(4326)

# 浸水域と交差する道路エッジを除去(形状は切らず、そのまま抜く)
bad = edges.intersects(tsunami_union)
edges_keep = edges.loc[~bad].copy()

# 使われるノードだけに絞る
used_node_ids = pd.unique(edges_keep[["u", "v"]].values.ravel())
nodes_keep = nodes.loc[nodes.index.isin(used_node_ids)].copy()

# サブグラフ再構築
G2 = ox.utils_graph.graph_from_gdfs(nodes_keep, edges_keep, graph_attrs=G.graph)

# =========================
# 3) 出発・目的地ノード
# =========================
home_point = Point(home_wgs[1], home_wgs[0])

# 目的地:安全な避難所のうち「自宅に最も近い」1つ(直線距離)
if len(safe_shelters) == 0:
    # もし全て浸水域なら、とりあえず最も近い避難所を使う(要検討)
    target_shelters = shelters.copy()
else:
    target_shelters = safe_shelters.copy()

# 距離計算はメートル系CRSで
METRIC_CRS = "EPSG:32654"  # 静岡周辺
tp = target_shelters.to_crs(METRIC_CRS)
hp = gpd.GeoSeries([home_point], crs=4326).to_crs(METRIC_CRS).iloc[0]
tp["__dist_m__"] = tp.geometry.distance(hp)
dest_pt = target_shelters.loc[tp["__dist_m__"].idxmin(), "geometry"]

# 近傍ノード
orig = ox.distance.nearest_nodes(G2, X=home_point.x, Y=home_point.y)
dest = ox.distance.nearest_nodes(G2, X=dest_pt.x, Y=dest_pt.y)

# =========================
# 4) 最短経路(浸水域エッジを除去済みグラフ)
# =========================
try:
    route = nx.shortest_path(G2, orig, dest, weight="length")
except nx.NetworkXNoPath:
    # 到達不能の場合は半径を増やす or 車道(network_type="drive")に切替など検討
    raise RuntimeError("安全な避難所へ到達するパスが見つかりません。半径を増やす/目的地を変更してください。")

# ルート形状(LineString)を作成
route_coords = [(G2.nodes[n]["x"], G2.nodes[n]["y"]) for n in route]
# GeoDataFrame化(WGS84)
route_gdf = gpd.GeoDataFrame(geometry=[gpd.GeoSeries([Point(xy) for xy in route_coords]).unary_union.convex_hull], crs=4326)
# ↑簡易化:見た目重視ならPolylineをFoliumへ直接渡す(下でやります)

# =========================
# 5) 出力(GeoPackage)& Folium 可視化
# =========================
# GPKGへ保存(必要なら既存を削除)
if os.path.exists(OUT_GPKG):
    os.remove(OUT_GPKG)

coast.to_file(OUT_GPKG, layer="coastline", driver="GPKG")
tsunami.to_file(OUT_GPKG, layer="tsunami_inundation", driver="GPKG")
shelters.to_file(OUT_GPKG, layer="shelters_all", driver="GPKG")
safe_shelters.to_file(OUT_GPKG, layer="shelters_safe", driver="GPKG")

# 走行可能エッジ(浸水交差なし)だけ保存
edges_keep.to_file(OUT_GPKG, layer="roads_safe", driver="GPKG")

# ルートはFolium用に直接線で描くので、GPKGは省略してもOK
# 保存したいなら下のようにLineStringを作って保存
from shapely.geometry import LineString
route_line = LineString(route_coords)
gpd.GeoDataFrame(geometry=[route_line], crs=4326).to_file(OUT_GPKG, layer="route", driver="GPKG")

print(f"✓ 出力: {OUT_GPKG}")

# Folium地図
m = folium.Map(location=[home_wgs[0], home_wgs[1]], zoom_start=12, control_scale=True)

# 津波浸水想定
folium.GeoJson(tsunami.__geo_interface__, name="津波浸水想定(A40)", tooltip=None).add_to(m)

# 海岸線
folium.GeoJson(coast.__geo_interface__, name="海岸線(C23)", tooltip=None).add_to(m)

# 避難所(安全/全体)
folium.GeoJson(safe_shelters.__geo_interface__, name="安全な避難所", 
               marker=folium.CircleMarker(radius=4),
               tooltip=folium.GeoJsonTooltip(fields=[:5])
               ).add_to(m)
folium.GeoJson(shelters.__geo_interface__, name="避難所(全体・灰)", 
               marker=folium.CircleMarker(radius=3, color="gray", fill=True),
               tooltip=folium.GeoJsonTooltip(fields=[:3])
               ).add_to(m)

# ルート線(G2ノード列から)
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(route_latlon, weight=5, opacity=0.9, tooltip="避難経路").add_to(m)

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

folium.LayerControl().add_to(m)
m  # Jupyterならこれで表示。スクリプトなら m.save("evac.html")

としたが

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[2], line 33
     31 coast = gpd.read_file(COAST_SHP)
     32 if coast.crs is None:
---> 33     raise ValueError("海岸線ShapefileにCRSがありません(.prjが同フォルダにあるか確認)")
     34 coast = coast.to_crs(4326)
     36 tsunami = gpd.read_file(TSUNAMI_SHP)

ValueError: 海岸線ShapefileにCRSがありません(.prjが同フォルダにあるか確認)

対処としては
海岸線(C23)の .prj が無いため CRS が未設定になっています。
C23 の公式仕様は JGD2000(経緯度) です。つまり EPSG:4612 を与えればOK。

対処方法は
A40(津波浸水想定)は配布物に .prj が入っていることが多いので、そのまま読める想定です。データ仕様の説明はここ(CRSは配布年度によりJGD2000/2011ですが、多くは .prj に記載あり)

coast = gpd.read_file(COAST_SHP)
if coast.crs is None:
    # C23 海岸線は JGD2000 (経緯度) → EPSG:4612
    coast = coast.set_crs(epsg=4612)
coast = coast.to_crs(4326)  # Folium等で扱いやすいWGS84へ

もしくは
コマンドで WKT を出して .prj として保存

# EPSG:4612 のWKTを出力(どちらか)
gdalsrsinfo EPSG:4612 -o wkt > C23-06_22-g_Coastline.prj
# または
projinfo EPSG:4612 --wkt2=EPSG > C23-06_22-g_Coastline.prj

念のための確認ポイント
* coast.total_bounds を見て、経度 ~122–154、緯度 ~20–46 の範囲に収まっていれば経緯度(度)の座標と判断可。
* メタデータXML(C23-06_22-g.xml など)に「座標系: JGD2000 / (B,L)」の記述があります。国土数値情報
これで ValueError は解消できます。次に進めて、浸水域クリップや経路探索を回してみてください

GPTで
対処(最短)と代替(.prj を作ってしまう) どちらが良い? その根拠は?
と調べてみる

結論:.prj を作って“直す”方(恒久対処)がオススメ。
理由は下記のとおりです。
どちらが良い?(要点)
* .prj を作成(恒久対処)
* 一度直せば QGIS/GeoPandas/同僚PC すべてで自動認識。再発しない。
* 将来の再利用・共有・配布に強い(ミスの温床を潰せる)。
* 前提:正しいCRSが確定していること(C23=JGD2000)。国土数値情報
* コード内で set_crs(暫定対処)
* いま手元のスクリプトをすぐ動かすなら最短。
* ただし、毎回コード側で設定を忘れるリスク、他ツール(QGIS等)を使うたびに同じエラーが再発。
このC23「海岸線」は仕様上 JGD2000(緯度経度)= EPSG:4612。
一方、A40「津波浸水想定」は JGD2011(緯度経度) が明記されています(年度によるが少なくとも2020版でJGD2011)。
→ どのみち最終的に同一CRS(例:WGS84=EPSG:4326)へ統一してから解析・可視化します。

これを実行するためのコマンドをインストール

arch -arm64 brew install gdal proj

🧭 1. PROJ(座標変換ライブラリ)
概要
PROJ は、地理座標系(緯度・経度など)と投影座標系(メートル単位の平面座標など)を相互変換するためのオープンソースライブラリです。
* URL: https://proj.org/
* 開発元: OSGeo Foundation
* 最新バージョン例: PROJ 9.x(2025年現在)
主な機能
機能 説明
座標変換 緯度経度(WGS84など)↔ 平面直角座標(UTM、EPSG:2451など)の変換
測地計算 2点間距離、方位角などを地球楕円体上で計算
EPSG対応 EPSGコード(例: 4326, 6674, 2451など)を使って簡単に座標系を指定可能
datum変換 世界測地系(WGS84)↔ 日本測地系(Tokyo Datum)などの変換
WKT/WKT2対応 座標系定義をWKT形式で入出力できる
利用例(コマンドライン)

projinfo EPSG:4326

→ WGS84の詳細を表示します。

🌍 2. GDAL(Geospatial Data Abstraction Library)
概要
GDAL は、地理空間データ(ラスタ・ベクタ形式)を読み書き・変換するためのライブラリ群です。
* URL: https://gdal.org/
* 開発元: OSGeo Foundation
* Pythonから利用するときは import gdal または from osgeo import gdal として使用可能。
主な機能
分類 説明
ファイル形式サポート GeoTIFF, Shapefile, GPKG, KML, GML, CSV, NetCDF など300以上
座標変換 内部的にPROJを使用して座標系を変換
ベクタ操作 属性情報付きの地理データ(例: Shapefile)を読み書き
ラスタ操作 画像や地形データ(DEMなど)の切り出し、リサンプリング、投影変換
コマンドツール ogr2ogr, gdalinfo, gdalwarp, gdalsrsinfo など多機能CLIツールを提供
利用例

gdalinfo sample.tif

→ GeoTIFFファイルのメタデータ(座標系・範囲など)を表示。

ogr2ogr -f "GPKG" output.gpkg input.shp

→ ShapefileをGeoPackage形式に変換。

🧠 まとめ
項目 内容
PROJ 座標変換ライブラリ(EPSGコード、測地変換)
GDAL 地理空間データ入出力ライブラリ(ShapefileやGeoTIFFの変換)
arch -arm64 Apple SiliconでARMネイティブ実行を強制
主な用途 GeoPandas, Rasterio, QGIS, OSMnx, PyProjなどのPython地理ライブラリの基盤

インストール完了後

==> Running `brew cleanup gdal`...
Disable this behaviour by setting `HOMEBREW_NO_INSTALL_CLEANUP=1`.
Hide these hints with `HOMEBREW_NO_ENV_HINTS=1` (see `man brew`).
==> Running `brew cleanup proj`...
==> No outdated dependents to upgrade!

となるので
バージョンを確認

 gdalinfo --version

を実行するとGPTで書いてあったが

GDAL 3.11.4 "Eganville", released 2025/09/04

 projinfo --version
Unrecognized option: --version
usage: projinfo [-o formats] [-k crs|operation|datum|ensemble|ellipsoid] [--summary] [-q]
                ([--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]) 
                [--spatial-test contains|intersects]
                [--crs-extent-use none|both|intersection|smallest]
                [--grid-check none|discard_missing|sort|known_available] 
                [--pivot-crs always|if_no_direct_transformation|never|{auth:code[,auth:code]*}]
                [--show-superseded] [--hide-ballpark] [--accuracy {accuracy}]
                [--allow-ellipsoidal-height-as-vertical-crs]
                [--boundcrs-to-wgs84]
                [--authority name]
                [--main-db-path path] [--aux-db-path path]*
                [--identify] [--3d]
                [--output-id AUTH:CODE]
                [--c-ify] [--single-line]
                --searchpaths | --remote-data |
                --list-crs [list-crs-filter] |
                --dump-db-structure [{object_definition} | {object_reference}] |
                {object_definition} | {object_reference} |
                (-s {srs_def} [--s_epoch {epoch}] -t {srs_def} [--t_epoch {epoch}]) |
                ({srs_def} {srs_def})

-o: formats is a comma separated combination of: all,default,PROJ,WKT_ALL,WKT2:2015,WKT2:2019,WKT1:GDAL,WKT1:ESRI,PROJJSON,SQL
    Except 'all' and 'default', other format can be preceded by '-' to disable them

list-crs-filter is a comma separated combination of: allow_deprecated,geodetic,geocentric,
geographic,geographic_2d,geographic_3d,vertical,projected,compound

{object_definition} might be a PROJ string, a WKT string, a AUTHORITY:CODE, or urn:ogc:def:OBJECT_TYPE:AUTHORITY::CODE

となる

正確には

proj          

を実行

Rel. 9.7.0, September 1st, 2025
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]

これで準備ができたので
.prj を作成
海岸線の .shp と同じフォルダ(= data/coastline_shizuoka/C23-06_22_GML)で、同じベース名の .prj を作ればOKです。
(津波浸水想定 A40 側は既に A40-16_22.prj があるので不要)

とのことなので

cd data/coastline_shizuoka/C23-06_22_GML

次にコマンドの実行だが
提案されたコマンドは2つ

# どちらか一方でOK(ESRI形式のほうがShapefileと相性が良い)
gdalsrsinfo EPSG:4612 -o wkt_esri > C23-06_22-g_Coastline.prj

# または

projinfo EPSG:4612 --wkt2=EPSG > C23-06_22-g_Coastline.prj

どちらが良い?

結論:gdalsrsinfo -o wkt_esri を推奨です。
理由(短く)
* Shapefile は ESRI 由来なので、.prj は ESRI 方言(wkt_esri)が最も相性◎。
* 互換性が広い(ArcGIS/QGIS/古いツールでもズレにくい)。
* EPSG WKT2 でも動きますが、古いツールで軸順やデータム解釈の差がまれに出ます。

ということなので

gdalsrsinfo EPSG:4612 -o wkt_esri > C23-06_22-g_Coastline.prj

を実行

結果は

head -n 3 data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.prj

で確認

GEOGCS["GCS_JGD_2000",
DATUM["D_JGD_2000",

しかし

python - <<'PY'
import geopandas as gpd
gdf = gpd.read_file("data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp")
print(gdf.crs) # → EPSG:4612 と出ればOK
PY
None

となる

ベース名が1文字でも違う(拡張子以外が完全一致でない)
GDAL/FionaがESRI方言(wkt_esri)をうまく読めない環境
.prjは正しいが、別の場所の同名SHPを読んでいる/パス違い
まれに .prj はあるが内容が空/壊れている
が主な原因らしい

ls data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.*
data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.dbf
data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.prj
data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp
data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shx

次に

ogrinfo -so data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp C23-06_22-g_Coastline


OGRがSRSを認識しているか

ogrinfo -so data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp C23-06_22-g_Coastline
INFO: Open of `data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp'
      using driver `ESRI Shapefile' successful.

Layer name: C23-06_22-g_Coastline
Metadata:
  DBF_DATE_LAST_UPDATE=1941-11-24
Geometry: Line String
Feature Count: 1141
Extent: (137.484966, 34.572136) - (139.176559, 35.144928)
Layer SRS WKT:
(unknown)
C23_001: String (5.0)
C23_002: String (1.0)
C23_003: String (4.0)
C23_004: String (46.0)
C23_005: String (1.0)
C23_006: String (10.0)
C23_007: String (5.0)

これらにより
.prj を置いたのに GDAL が Layer SRS WKT: (unknown)。
これは その環境の GDAL/Fiona が WKT をうまく解釈できていない ケース

ogr2ogr -a_srs で“CRS付きの新ファイル”を作るのが一番おすすめです。
理由:
* その場しのぎではなく恒久対処(QGIS/GeoPandas/他PCでも確実に認識)。
* .prj の方言差や環境差に左右されにくい(今回まさに .prj を置いても読まれなかった)。
* 元データを壊さず安全に複製できる。

ということで
CRSを付与した“新しい”SHPを作成(EPSG:4612 = JGD2000)

cd data/coastline_shizuoka/C23-06_22_GML

で移動

ogr2ogr -a_srs EPSG:4612 C23_4612.shp C23-06_22-g_Coastline.shp

しかし

Warning 1: One or several characters couldn't be converted correctly from UTF-8 to ISO-8859-1.  This warning will not be emitted anymore.

となる

この警告は属性の文字エンコード(日本語)が
Shapefile の既定(ISO-8859-1)に合わず、
変換できない文字があったという意味で
幾何やCRS付与は完了しているが、属性が文字化けする可能性がある

解決策は
1) Shapefileのエンコーディングを指定して作り直す(手早い)
2) いっそ GeoPackage にする(おすすめ:今後楽)

どの方法が良いかを相談

コンセプトが

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

✅ あなたの用途の要件
* 精度:震源→沿岸までの到達時間・波高推定の入力に使うので、誤差が大きすぎると危険。
* 範囲:主に日本周辺で十分。全球は必須ではない。
* 更新頻度:将来的に 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でノード/リレーションを整理して「質問に答えられる防災知識ベース」を構築
なら、どの方法を推奨する?

結論:公的データ優先の“ハイブリッド型”を強く推奨します。
「浸水範囲は公式シミュ結果(A40)」「到達時刻・高さは公的速報を取り込む」「経路はOSM+浸水域回避」で組むのが、安全・現実的・拡張容易です。
推奨アーキテクチャ(段階導入)
フェーズ1:静岡・オフライン基盤(最小成功)
* 海岸線:KSJ C23(CRS整備済み)
* 津波浸水想定:KSJ A40(静岡)を絶対基準に採用
* 道路:OSMnx(walk/drive)取得 → A40と交差するエッジを除去
* 避難所:自治体CSV(cp932)→ GeoDataFrame(WGS84)
* 結果:
* 「自宅 → 最寄り“安全な”避難所」の最短経路
* 到達圏(isochrone):歩行/車速で10分・20分…を事前計算
* 保存:1つの GeoPackage(coastline / tsunami / roads_safe / shelters_safe / routes / isochrones)
ここまでで、「地震が来たら最寄りの避難所は?」に即答できます。
高度な波高計算はまだ入れず、安全サイドの浸水域回避を徹底。
フェーズ2:リアルタイム連携(到達“時刻”の提示)
* 到達時刻・津波高の“情報源”は公的速報(気象庁の津波警報・到達予想)。
* ワークフロー:
1. 速報の沿岸予報区ごとの到達見込み時刻/高さを取得
2. 県内の該当海岸セグメントへ付与(C23を海岸線ノードに分割)
3. 自宅の最短海岸ノードまでの**“海→陸の到達時刻”を速報値**で表示
4. ルート案内UIに「残り時間(ETA)」を明示(保守的に短めに表示)
物理モデルでの即時計算は難度・不確実性が高いので、公式速報を真に受ける方が安全かつ信用性◎。
どうしても“目安”を補うなら、深海平均水深での線形速波
c
=
g
h
c=gh による暫定下限ETAを併記(保守的に短く)。
フェーズ3:GraphDB連携(質問に答える防災知識ベース)
* ノード:Home / Shelter / RoadNode|RoadEdge / TsunamiPolygon(A40) / CoastNode / Item(備蓄) / Contact
* リレーション:ROUTE_TO(time,length) / INTERSECTS / AVOIDS / NEAREST / EXPIRING_SOON
* QA例:
* 「津波到達まで何分?」→ 最近接 CoastNode の公式ETA(フェーズ2)
* 「最短でどの避難所?」→ ROUTE_TO の最小 time を返す
* 「非常食の賞味期限近い?」→ EXPIRING_SOON <= 14d の Item 列挙 この方式を推す理由 * 安全性:浸水範囲は公式A40で担保。到達時刻・高さも公的速報準拠。独自計算の誤差リスクを避ける。 * 実装容易 & 拡張性:いまのPython+OSMnxコードの発展で実装可能。将来は別県や別ハザード(洪水・高潮)も流用可。 * 運用:GeoPackageに集約→QGIS/GeoPandas両対応。GraphDBで家族向けの質問応答が作りやすい。 実装メモ(最短でやるなら) * まずはフェーズ1を静岡で完成:A40で道路遮断 → 最短経路 & isochrone 出力。 * フェーズ2は速報取り込みのI/Oだけ追加し、海岸ノードにETA属性を付与。UIで「残り時間」を表示。 * すべて保守的(安全側)に倒す:ETAは短め、通行可能性に不確実があれば通れない扱い。 要約:
あなたの要求(精度・安全・拡張性)に最も合うのは、A40×OSM×公的速報のハイブリッド。
先に静岡でフェーズ1を完成させ、問題なければ**速報連携(フェーズ2)→GraphDB化(フェーズ3)**へ進めましょう。 となる つまり、今回なら 1) Shapefileのエンコーディングを指定して作り直す(手早い) 2) いっそ GeoPackage にする(おすすめ:今後楽) はどちらにする? それとも別の方法をとる? 結論:② GeoPackage(GPKG)に統一するのがベストです。 理由:文字化け・.prj問題・属性名10文字制限・2GB制限など Shapefileの弱点を一掃でき、今後の拡張(複数レイヤ管理、QGIS/GeoPandas/GraphDB連携)が安定します やることは 出力は一つの outputs/shizuoka_hazard.gpkg にレイヤ分けして格納 らしいが # 1) 海岸線(C23):CRS未付与なので付けてから WGS84 へ投影して GPKG へ を実行する前に気になるので確認しておく 現在 # 1) 海岸線(C23):CRS未付与なので付けてから WGS84 へ投影して GPKG へ[shell] ogr2ogr -a_srs EPSG:4612 -t_srs EPSG:4326 \ -f GPKG outputs/shizuoka_hazard.gpkg \ data/coastline_shizuoka/C23-06_22_GML/C23-06_22-g_Coastline.shp \ -nln coastline [/shell] # 2) 津波浸水想定(A40):.prj あり → そのCRSから WGS84 へ投影して同じ GPKG に追加[shell] ogr2ogr -t_srs EPSG:4326 \ -f GPKG -update outputs/shizuoka_hazard.gpkg \ data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp \ -nln tsunami_inundation[/shell] を実行する前のフォルダ構成として [shell] data/coastline_shizuoka/C23-06_22_GML C23-06_22-g.xml C23_4612.dbf C23-06_22-g_Coastline.dbf C23_4612.prj C23-06_22-g_Coastline.prj C23_4612.shp C23-06_22-g_Coastline.shp C23_4612.shx C23-06_22-g_Coastline.shx KS-META-C23-06_22-g.xml[/shell] となっているが .prj ファイルは一度削除した方が良いか? 結論:消さなくてOKです。 大事なのは「どのSHPを入力に使うか」と「source SRS を明示して出力する」こと。GPKGを作るときに正しく指定すれば、フォルダ内に .prj が残っていても問題ありません 既にCRS付きコピーを作ってあるので、それを入力にして WGS84 で GPKG に出力すればOK 先に outputsフォルダの作成が必要なので [shell] cd -[/shell] で移動しておく ~/aw10s/FamilySafetyMap へ移動して [shell] mkdir -p outputs[/shell] でフォルダ作成し[shell] ogr2ogr -t_srs EPSG:4326 \ -f GPKG outputs/shizuoka_hazard.gpkg \ data/coastline_shizuoka/C23-06_22_GML/C23_4612.shp \ -nln coastline[/shell] で 既にCRS付きコピーを作ってあるので、それを入力にして WGS84 で GPKG に出力。 動作確認[shell] ogrinfo -so outputs/shizuoka_hazard.gpkg coastline INFO: Open of `outputs/shizuoka_hazard.gpkg' using driver `GPKG' successful. Layer name: coastline Metadata: DBF_DATE_LAST_UPDATE=2025-10-10 Geometry: Line String Feature Count: 1141 Extent: (137.484966, 34.572136) - (139.176559, 35.144928) Layer SRS WKT: GEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)"], MEMBER["World Geodetic System 1984 (G730)"], MEMBER["World Geodetic System 1984 (G873)"], MEMBER["World Geodetic System 1984 (G1150)"], MEMBER["World Geodetic System 1984 (G1674)"], MEMBER["World Geodetic System 1984 (G1762)"], MEMBER["World Geodetic System 1984 (G2139)"], MEMBER["World Geodetic System 1984 (G2296)"], ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], USAGE[ SCOPE["Horizontal component of 3D system."], AREA["World."], BBOX[-90,-180,90,180]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 FID Column = fid Geometry Column = geom C23_001: String (5.0) C23_002: String (1.0) C23_003: String (4.0) C23_004: String (46.0) C23_005: String (1.0) C23_006: String (10.0) C23_007: String (5.0) ogrinfo -so outputs/shizuoka_hazard.gpkg tsunami_inundation INFO: Open of `outputs/shizuoka_hazard.gpkg' using driver `GPKG' successful. ERROR 1: Couldn't fetch requested layer tsunami_inundation. [/shell] となった tsunami_inundation レイヤが GPKGに存在していない いま入っているレイヤ一覧を確認[shell] ogrinfo outputs/shizuoka_hazard.gpkg INFO: Open of `outputs/shizuoka_hazard.gpkg' using driver `GPKG' successful. 1: coastline (Line String) [/shell] A40 のソースCRSを確認[shell] ogrinfo -so data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp A40-16_22 INFO: Open of `data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp' using driver `ESRI Shapefile' successful. Layer name: A40-16_22 Metadata: DBF_DATE_LAST_UPDATE=2016-12-08 Geometry: Polygon Feature Count: 113272 Extent: (137.483132, 34.594813) - (139.176599, 35.163526) Layer SRS WKT: GEOGCRS["JGD2011", DATUM["Japanese Geodetic Datum 2011", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], USAGE[ SCOPE["Horizontal component of 3D system."], AREA["Japan - onshore and offshore."], BBOX[17.09,122.38,46.05,157.65]], ID["EPSG",6668]] Data axis to CRS axis mapping: 2,1 A40_001: String (10.0) A40_002: String (2.0) A40_003: String (40.0) [/shell] この結果から A40 はちゃんと EPSG:6668 (JGD2011) で読めています。あとは GPKG に “追記” できればOK 進捗表示・無効ジオメトリを自動修復・マルチ化・空間インデックス作成する [shell] ogr2ogr -progress -t_srs EPSG:4326 \ -makevalid -nlt PROMOTE_TO_MULTI \ -f GPKG -update -lco SPATIAL_INDEX=YES -lco GEOMETRY_NAME=geom \ outputs/shizuoka_hazard.gpkg \ data/tsunami_shizuoka/A40-16_22_GML/A40-16_22.shp \ -nln tsunami_inundation [/shell] を実行 補足(なぜこのオプション?) * -update:既存 GPKG に 追記 * -nln tsunami_inundation:レイヤ名を固定 * -t_srs EPSG:4326:WGS84 に統一(Folium/QGISで扱いやすい) * -makevalid:自己交差などの 無効ポリゴンを自動修復(A40は件数多く、稀に必要) * -nlt PROMOTE_TO_MULTI:MultiPolygon に正規化(後段の処理が安定) * -lco SPATIAL_INDEX=YES:高速化のための 空間インデックス 付与 次に動作確認 [shell] ogrinfo -so outputs/shizuoka_hazard.gpkg tsunami_inundation INFO: Open of `outputs/shizuoka_hazard.gpkg' using driver `GPKG' successful. Layer name: tsunami_inundation Metadata: DBF_DATE_LAST_UPDATE=2016-12-08 Geometry: Multi Polygon Feature Count: 113272 Extent: (137.483133, 34.594813) - (139.176599, 35.163526) Layer SRS WKT: GEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)"], MEMBER["World Geodetic System 1984 (G730)"], MEMBER["World Geodetic System 1984 (G873)"], MEMBER["World Geodetic System 1984 (G1150)"], MEMBER["World Geodetic System 1984 (G1674)"], MEMBER["World Geodetic System 1984 (G1762)"], MEMBER["World Geodetic System 1984 (G2139)"], MEMBER["World Geodetic System 1984 (G2296)"], ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ENSEMBLEACCURACY[2.0]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], USAGE[ SCOPE["Horizontal component of 3D system."], AREA["World."], BBOX[-90,-180,90,180]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 FID Column = fid Geometry Column = geom A40_001: String (10.0) A40_002: String (2.0) A40_003: String (40.0) [/shell] これだけだとわからないのでGPTで解説を依頼 これは `ogrinfo -so` コマンド(GDALのツール)で、 GeoPackage ファイル(`outputs/shizuoka_hazard.gpkg`)のレイヤー `tsunami_inundation`(津波浸水想定区域)に関する**メタデータと構造情報**を表示した結果です。 以下に各項目の意味を解説します👇 --- ## 🗂 基本情報 | 項目 | 内容 | | ------------------------------------------------------------- | ------------------------------------------------------------------------------------ | | **Open of ... using driver 'GPKG' successful** | GeoPackage 形式(`.gpkg`)のファイルを正常に開けたことを示します。GDALドライバ `GPKG`(GeoPackageフォーマット)を使用しています。 | | **Layer name: tsunami_inundation** | レイヤー名(地理データの層の名前)。ここでは津波浸水区域を示すレイヤーです。 | | **Geometry: Multi Polygon** | 地物のジオメトリ(形状)が「複数ポリゴン」で構成されていることを示します。津波の浸水区域など、面(エリア)を表現するために使われます。 | | **Feature Count: 113272** | 含まれる地物(フィーチャー)の総数。ここでは 113,272 個の津波浸水ポリゴンが含まれています。 | | **Extent: (137.483133, 34.594813) - (139.176599, 35.163526)** | データがカバーしている範囲の外接矩形(最小緯度経度と最大緯度経度)。静岡県の東西南北の境界座標です。 | | **Layer SRS WKT** | 空間参照系(Spatial Reference System)の詳細情報。ここでは世界測地系(WGS84 / EPSG:4326)です。 | --- ## 🌍 空間参照系(SRS: Spatial Reference System) ``` GEOGCRS["WGS 84", ... ID["EPSG",4326]] ``` これは座標系の定義で、簡単に言うと: * **EPSG:4326** → 世界測地系(WGS 84)。 GPSなどで使われる緯度経度座標系。 * **CS[ellipsoidal,2]** → 2次元の楕円体座標系(緯度・経度)。 * **PRIMEM["Greenwich",0]** → 基準子午線はグリニッジ子午線。 * **ANGLEUNIT["degree",...]** → 度単位で表されている。 つまり、データは**緯度経度(度単位)**で表現されており、地図投影(メートル単位の平面座標変換)はまだ行われていません。 --- ## 📊 属性カラム(テーブル構造) | カラム名 | データ型 | 意味(推定) | | ----------- | ------------- | -------------------- | | **A40_001** | String (10.0) | コードやID(市区町村や地域コードなど) | | **A40_002** | String (2.0) | 区分や種別(例:津波の想定区分など) | | **A40_003** | String (40.0) | 名称(例:「静岡県○○市△△地区」など) | これらのフィールド名は国土数値情報(国交省が公開する津波浸水想定データ)と一致しています。 具体的な意味はデータ仕様書により異なりますが、概ね以下の通りです: | フィールド名 | 意味(国交省仕様例) | | ------- | --------------------- | | A40_001 | 行政コード(市区町村コード) | | A40_002 | 津波区分(1:想定, 2:最大クラスなど) | | A40_003 | 地域名・想定名などのラベル | --- ## 🧭 Data axis to CRS axis mapping: 2,1 これは軸の並び順を示します。 * 通常は **(緯度, 経度)** の順ですが、 * このデータでは **(経度, 緯度)** の順に格納されていることを意味します。 GDALで処理する際、軸の順序に注意が必要です(例:FoliumやGeoPandasで描画する際に位置がズレる場合があります)。 --- ## 🔍 FID / Geometry Column | 項目 | 内容 | | -------------------------- | -------------------------- | | **FID Column = fid** | 各地物を識別する一意のID(Feature ID)。 | | **Geometry Column = geom** | ポリゴン形状が格納されているカラム名。 | --- ## 🧩 まとめ このGeoPackageは、 **「静岡県内の津波浸水想定区域」**をWGS84座標で格納したもので、 113,272件のポリゴンが `geom` フィールドに保存されています。 主な使い方の例: [shell] import geopandas as gpd gdf = gpd.read_file("outputs/shizuoka_hazard.gpkg", layer="tsunami_inundation") print(gdf.head()) ``` 出力例: | fid | A40_001 | A40_002 | A40_003 | geometry | | --- | ------- | ------- | ------- | ----------------- | | 1 | 22210 | 1 | 静岡市駿河区 | MULTIPOLYGON(...) | [/shell] ご希望であれば、このデータを **Folium地図上に可視化するコード例** や **GeoPandasでの属性抽出・座標変換(EPSG:4326 → EPSG:6677など)** の実例も示せます

海岸線データの取得

海岸線データの取得

import numpy as np
import geopandas as gpd
import rasterio
from rasterio import features
from shapely.geometry import Point
from skimage.graph import route_through_array  # 1対1のとき
import osmnx as ox
import networkx as nx
import folium

# --------------------
# 0) 入力
# --------------------
epicenter_wgs = (35.5, 142.0)  # 例: 緯度, 経度(仮)
gebco_tif = "data/gebco_bathymetry.tif"   # 海底地形(m, 海は正の深さ or 負の標高のことも)
dem_tif   = "data/dem_coast.tif"          # 陸域標高(m)
coastline = gpd.read_file("data/coastline.gpkg")  # 海岸線
shelters  = gpd.read_file("data/shelters.geojson")# 避難所点データ(WGS84想定)

# 便利関数:WGS84→ラスタ座標系へ
def to_crs_point(lat, lon, crs_epsg):
    g = gpd.GeoSeries([Point(lon, lat)], crs=4326).to_crs(crs_epsg)
    p = g.iloc[0]
    return p.y, p.x  # (Y, X) in projected CRS

# --------------------
# 1) 海上到達時間ラスタ(簡易)
# --------------------
with rasterio.open(gebco_tif) as src:
    bathy = src.read(1).astype("float32")
    prof  = src.profile
    # 海域の水深(m)に統一(負の標高なら反転)
    depth = np.where(bathy < 0, -bathy, bathy)  # 例:負値を深さに
    depth = np.where(depth<=0, np.nan, depth)

    g = 9.81
    v = np.sqrt(g * depth)  # m/s
    cost = 1.0 / v          # s/m(セル通過時間はこれ×セルサイズ)

    # 震源セル
    epic_y, epic_x = to_crs_point(epicenter_wgs[0], epicenter_wgs[1], src.crs.to_epsg())
    row, col = src.index(epic_x, epic_y)

    # coast マスク(海→陸境界の近傍)を後で使うためにベクタ→ラスタ化しておくのも可

# ここではデモとして「震源→任意の海岸点」1対1のrouteで travel time を出す例。
# 実際は multi-source 最短時間(fast marching/多点Dijkstra)が望ましい。
# skimage.graphの route_through_array はコスト配列にセルサイズを掛けて合計。
# セルサイズ(m):
pix_m = (prof["transform"][0])  # x解像度(m/px)。投影ラスタ前提
# 到達時間ラスタの厳密生成はここでは割愛(本格的には多源最短経路/コスト距離を実装)

# --------------------
# 2) 浸水域(閾値で簡易)
# --------------------
with rasterio.open(dem_tif) as demsrc:
    dem = demsrc.read(1).astype("float32")
    # 仮の沿岸波高 H(m)を設定(現実は断層モデルや公表想定を用いる)
    H0 = 1.0      # 沖合初期波高 1 m(例)
    h0 = 3000.0   # 沖合水深 3000 m(例)
    h_coast = 10.0# 沿岸水深(例)
    H_coast = H0*((h0/h_coast)**0.25)

    tide = 0.5    # 満潮加算(例)
    safety_margin = 0.5
    flood_threshold = H_coast + tide + safety_margin  # m

    # 海から連結する低地 = 浸水候補
    flood_mask = (dem <= flood_threshold).astype(np.uint8)

# flood_mask をポリゴン化(例)
shapes = features.shapes(flood_mask, transform=demsrc.transform)
flood_polys = [sh for sh, val in shapes if val == 1]
flood = gpd.GeoDataFrame(geometry=gpd.GeoSeries.from_iterable(flood_polys), crs=demsrc.crs)

# 海から連結していない内陸窪地は除外(簡易には海ポリゴンとの接触でフィルタ)
# coastlineを同一CRSへ
coastline_proj = coastline.to_crs(flood.crs)
flood = flood[flood.intersects(coastline_proj.unary_union)].buffer(0)

# --------------------
# 3) 安全エリア&避難所(浸水域の外)
# --------------------
shelters_proj = shelters.to_crs(flood.crs)
safe_shelters = shelters_proj[~shelters_proj.intersects(flood.unary_union)]

# --------------------
# 4) 避難ルート(浸水域の道路を除外)
# --------------------
place = "Sumida-ku, Tokyo, Japan"
G = ox.graph_from_place(place, network_type="walk")  # 徒歩網の方が避難向き

nodes, edges = ox.graph_to_gdfs(G)
edges_proj = edges.to_crs(flood.crs)

# 浸水と交差するエッジを除外
bad = edges_proj.intersects(flood.unary_union)
edges_keep = edges_proj.loc[~bad]

# サブグラフ再構築
G2 = ox.utils_graph.graph_from_gdfs(nodes.to_crs(edges_keep.crs), edges_keep)

# 出発地(例:現在地)と最寄り安全避難所
current_wgs = (35.71, 139.81)  # 例
cur_pt = gpd.GeoSeries([Point(current_wgs[1], current_wgs[0])], crs=4326).to_crs(edges_keep.crs).iloc[0]
dest_pt = safe_shelters.geometry.iloc[0] if len(safe_shelters) else shelters_proj.geometry.iloc[0]

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

route = nx.shortest_path(G2, orig, dest, weight="length")

# --------------------
# 5) Foliumで可視化
# --------------------
center = (current_wgs[0], current_wgs[1])
m = folium.Map(location=center, zoom_start=13)

# 浸水域
flood_wgs = flood.to_crs(4326)
folium.GeoJson(flood_wgs.__geo_interface__, name="Flood (simple)").add_to(m)

# ルート
route_latlon = [(G2.nodes[n]["y"], G2.nodes[n]["x"]) for n in route]
folium.PolyLine(route_latlon, color="red", weight=5, opacity=0.9, tooltip="Evacuation Route").add_to(m)

# 出発・目的地
folium.Marker(location=[cur_pt.y, cur_pt.x], popup="Start", icon=folium.Icon(color="green")).add_to(m)
folium.Marker(location=[dest_pt.to_crs(4326).y, dest_pt.to_crs(4326).x], popup="Shelter", icon=folium.Icon(color="blue")).add_to(m)

# 避難所レイヤ
folium.GeoJson(safe_shelters.to_crs(4326).__geo_interface__, name="Safe shelters",
               tooltip=folium.GeoJsonTooltip(fields=[:4])
               ).add_to(m)

folium.LayerControl().add_to(m)

# Jupyterなら m だけで表示
m  # もしくは m.save("evac.html")

を実行しても

DataSourceError: data/coastline.gpkg: No such file or directory 

となるため
データの取得先を調べる

data/coastline.gpkg は自分で用意する必要があります。用途(全球/日本のみ、精度、更新頻度)に応じて、下のどれかをダウンロード→GeoPackage に変換すれば OK です。
入手先(おすすめ順)
1. Natural Earth – Coastline(手軽・軽量)
1:10m/50m/110m の海岸線 Shapefile。全球で数MB、試作には十分です。Natural Earth Data+2Natural Earth Data+2
2. OSM Coastlines / Land polygons(高頻度更新・実用)
OpenStreetMap から抽出された最新の海岸線。WGS84 / Web Mercator の Shapefile がそのまま落とせます。OpenStreetMap Data+1
3. GSHHG(高解像度・全球)
解像度別(c/l/i/h/f)で海岸線を配布。Shapefile もあり。厳密な可視化に向く一方、古い部分もあるので用途に注意。ngdc.noaa.gov+2ngdc.noaa.gov+2
4. 国土数値情報「海岸線データ」(日本限定・行政用途向け属性つき)
都道府県ごとに海岸線(線)を提供。JPGIS/GML と Shapefile が選べます。nlftp.mlit.go.jp+2geospatial.jp+2

使うときの注意
* CRS を統一:あなたの dem_tif / gebco_tif / 道路網(OSMnx)と同じ座標系に coastline を揃えてください。
例:coastline = coastline.to_crs(flood.crs)(すでにコードでやっている通り)
* 粒度の選択:
* 試作:Natural Earth(軽い)
* 実運用・最新形状:OSM Coastlines
* 日本限定・行政属性重視:国土数値情報
* ライン or ポリゴン:浸水域との連結判定には、**海岸線(ライン)**だけでなく 陸/水ポリゴンがあると便利です(OSM の Land/Water polygons 参照)。

種類が多いためGPTで候補を割り出す

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

✅ あなたの用途の要件
* 精度:震源→沿岸までの到達時間・波高推定の入力に使うので、誤差が大きすぎると危険。
* 範囲:主に日本周辺で十分。全球は必須ではない。
* 更新頻度:将来的に 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でノード/リレーションを整理して「質問に答えられる防災知識ベース」を構築