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

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

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保存 → 追加の堅牢化まで整う