震源地から津波が到達するまでの時間と津波の高さを計算

震源地から津波が到達するまでの時間と津波の高さを計算

✅ 津波の到達時間の基本式
津波は浅水波近似で速度が次の式で表されます：
v=√g⋅h

v：津波の速度 (m/s)
g：重力加速度 (9.81 m/s²)
h：海の水深 (m)

例
水深 4000 m → 速度 ≈ 200 m/s (720 km/h)
水深 100 m → 速度 ≈ 31 m/s (110 km/h)

✅ 津波の高さ（沿岸での増幅）
震源での初期波高
H
0
H0 は断層変位や地殻変動から推定されます。
沿岸での津波高さ
H
H は「グリーンの法則」によって水深比に応じて増幅します：

h0 ：沖合の水深
h：沿岸の水深

つまり、水深が浅くなるにつれて波高は増大します。

✅ Pythonでの計算イメージ

import math

def tsunami_travel_time(distance_km, depth_m):
    g = 9.81
    v = math.sqrt(g * depth_m)  # m/s
    t_sec = (distance_km*1000) / v
    return t_sec / 60  # 分

def tsunami_height(H0, h0, h):
    return H0 * ((h0 / h) ** 0.25)

# 例: 震源から200 km, 平均水深 3000 m
t = tsunami_travel_time(200, 3000)
print(f"到達時間 ≈ {t:.1f} 分")

# 例: 初期波高 1m, 沖合水深 3000 m, 沿岸水深 10 m
H = tsunami_height(1.0, 3000, 10)
print(f"沿岸波高 ≈ {H:.2f} m")

✅ データが必要になる部分
* 震源位置と震源から沿岸までの距離 → ジオデータ（緯度経度）を使って測る
* 海底地形（バソメトリー）データ → GEBCO や NOAA の水深データが必要
* 初期波高 → 地震断層モデルから計算（研究所や気象庁の速報値を利用するのが実務的）

✅ まとめ
* 津波高さは「水深比の1/4乗則（グリーンの法則）」で沿岸増幅を近似
* 実務では GEBCO の海底地形データを組み合わせて地図上で計算する


震源地から津波が到達するまでの時間と津波の高さを計算し GIS情報を使い避難するエリアを算出しルート案内できるようにしたい

精密な津波シミュレータではなく、GIS＋簡易物理モデルで「到達時間」「想定浸水域」「安全エリア」「避難ルート」を一気通貫で出す手順と最小コード


※簡易近似。自治体ハザードマップがあればそれを優先
1. 想定浸水域の作成（簡易版）
陸域DEMから標高閾値（例：満潮＋推定波高＋α）以下の連結領域を抽出し、海から連結する部分を浸水ポリゴンとする。
（精密な遡上計算は専門コードが必要。ここでは避難計画用の予防的ゾーニング）
2. 安全エリアの抽出
浸水ポリゴン外かつ標高>閾値（例：5–10m以上）＆避難所候補のある場所を安全エリアとする。
3. 避難ルート
OSMnxで道路グラフを取得し、浸水ポリゴン内のエッジを除外して、現在地→最寄り安全避難所までの最短路（距離 or 時間）を算出。Foliumで地図化。

必要なパッケージ

pip install geopandas rasterio shapely numpy scikit-image osmnx folium contextily networkx

でインストール

ミニマル実装（雛形）
ファイルパスはあなたの環境に合わせて置き換え。座標系は**投影座標（m単位）**に揃えるのがコツ（例：対象地域のUTM帯、または EPSG:3857）。日本なら JGD2011 の平面直角系でもOK。

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")

実装のコツ
* CRSを徹底：ラスタ・ベクタ・OSMのCRSを揃える（距離・面積は投影座標）
* 閾値は保守的に：不確実性を考慮し安全側に（潮位・護岸・地形・内水氾濫は別論点）
* 道路グラフのカット：浸水ポリゴンだけでなく河川・橋梁も個別判定できると尚良
* 到達時間：本気でやるなら海上は1/vのコスト面で多源最短時間、陸上は時々刻々の浸水到達時刻 vs. 歩行速度で**到達可能圏（isochrone）**を比較

重要な注意（精度・安全）
* ここで示したのは教育・訓練用の簡易推定です。実運用や避難計画は、自治体のハザードマップ・気象庁/JMA・国研のモデル結果を最優先してください。
* グリーンの法則や低地抽出は大まかな上限見積に有効ですが、局所地形（湾形状、護岸、河口、運河）で大きな差が出ます。

次の拡張
海上到達時間：ラスタ全域での多源Dijkstra / Fast Marching（skfmm等）で全海岸の到達時刻マップを作る
公式ハザードデータ（津波浸水想定、指定緊急避難場所）との連携


とりあえずコードを jupyter で実行してみる