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")
# 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)
# 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 を作成")
# # ルートの延長[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, 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)
の後に
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"]]
を実行すると
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[31], line 2
1 import geopandas as gpd
----> 2 gpd.read_file("outputs/shizuoka_hazard.gpkg", layer="route_latest")[["meters","eta_min_walk","eta_min_walk_safe","eta_min_drive"]]
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: "['eta_min_walk_safe'] not in index"
となった
GPTによれば
原因は2つのどちらか（または両方
途中の「一回目の書き出し」で eta_min_walk_safe を含まないスキーマが作られた
最初に meters / eta_min_walk / eta_min_drive だけを持つ route_gdf を
書き出す箇所があります（OVERWRITE なし）。
この時点で route_latest の列定義が確定 します。
その後に eta_min_walk_safe を追加して上書きしても、環境によっては列が残らない／読み出し側で見えないことがあります
route_gdf.to_file("outputs/shizuoka_hazard.gpkg", layer="route_latest", driver="GPKG")
がその部分
次の理由として
キャッシュヒット分岐では route_latest を書き直していない
ルート: キャッシュヒット になった場合、後段の
# 6) route_latest へ保存（…） route_gdf.to_file(... OVERWRITE="YES")
がが実行されません。結果、古いスキーマのまま
このため
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")
の中の
# ❌ これを削除：
route_gdf.to_file("outputs/shizuoka_hazard.gpkg", layer="route_latest", driver="GPKG")
をコメントアウト
どの分岐でも同じスキーマで必ず上書き保存するように
共通の保存関数を作り、キャッシュヒット時も呼ぶようにする
from datetime import datetime
def save_route_latest(out_gpkg, route_geom, meters, eta_walk_min, eta_drive_min, eta_walk_safe):
gdf = gpd.GeoDataFrame(
{
"meters":[float(meters)],
"eta_min_walk":[float(eta_walk_min)],
"eta_min_walk_safe":[float(eta_walk_safe)],
"eta_min_drive":[float(eta_drive_min)],
"created_at":[datetime.utcnow().isoformat(timespec="seconds")+"Z"],
},
geometry=[route_geom], crs="EPSG:4326"
)
gdf.to_file(
out_gpkg, layer="route_latest",
driver="GPKG", engine="pyogrio",
layer_options={"OVERWRITE": "YES"}
)