海岸線データの取得

海岸線データの取得

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でノード/リレーションを整理して「質問に答えられる防災知識ベース」を構築