投稿者: Snowpool

海岸線データの取得

海岸線データの取得

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

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

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

✅ 津波の到達時間の基本式
津波は浅水波近似で速度が次の式で表されます：
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 で実行してみる

OSMxを使い最短経路の算出

OSMxを使い最短経路の算出

東京スカイツリーから
江戸東京博物館への最短経路を算出する

まずはスカイツリーと江戸東京博物館の座標をセットする

skytree = (35.7100,139.8180)
museum = (35.6963,139.7967)

次に一番近い node つまり交差点を取得

始点はスカイツリーなので

start = ox.get_nearest_node(G,skytree)

としたがエラー

--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[15], line 1 ----> 1 start = ox.get_nearest_node(G,skytree) AttributeError: module 'osmnx' has no attribute 'get_nearest_node'

これもバージョンの違いが原因

OSMnx 2.x では廃止され、ox.distance サブモジュールに移動
引数は (X=経度, Y=緯度) の順番で渡す必要があります
普通の Point(lon, lat) と同じ並び

start = ox.distance.nearest_nodes(G, X=skytree[1], Y=skytree[0])

とすればOK

同様に終点も

end = ox.distance.nearest_nodes(G,museum[1],museum[0])

次に最短経路

shortest_path = ox.shortest_path(G,start,end)

Foliumで表示できないので調べた

plot_route_folium
https://github.com/gboeing/osmnx/issues/1252
によれば

release (currently 1.9.4), then run your code. The latest v1 release will remain compatible with the v1 API, but any deprecated functionality you're using will issue a FutureWarning with simple guidance to make it v2 compatible.
Yes, the folium module was deprecated with a warning for removal back in 2023, and finally removed in v2.0.0. The good news: interactive web mapping with folium remains easy to do with OSMnx in ~2 lines of code! See this example in the Examples Gallery for a demonstration.

となっていて
foliumモジュールは2023年に削除の警告とともに非推奨となり、最終的にv2.0.0で削除されました。

とのこと

https://github.com/gboeing/osmnx-examples/blob/main/notebooks/11-interactive-web-mapping.ipynb
のサンプルが参考になりそう

OSMnx 2.0.6 には Folium 専用関数は無い → 自分で PolyLine を追加すればOK

再度GPTでコード生成

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

というように
weight=”length”)
があるので調べた

✅ weight の意味
networkx.shortest_path の weight 引数は、**「どのエッジ属性を最短経路計算に使うか」**を指定するものです。
OSMnx で作成した G（道路グラフ）には、各エッジ（道路の区間）に様々な属性があります：
* “length” : 区間の長さ（メートル単位） ← デフォルトでよく使う
* “travel_time” : 区間の推定走行時間（秒単位）※ ox.add_edge_speeds / ox.add_edge_travel_times を実行した場合のみ
* その他、highway（道路種別）や lanes（車線数）など

✅ 例
1. 最短距離ルート（距離を最小化）

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

2. 最短時間ルート（所要時間を最小化）
まず速度と所要時間をエッジに付与してから：

import osmnx as ox

G = ox.add_edge_speeds(G)        # 各道路に推定速度 (km/h) を追加
G = ox.add_edge_travel_times(G)  # 各道路に travel_time (秒) を追加

route = nx.shortest_path(G, orig, dest, weight="travel_time")

✅ まとめ
* weight=”length” → 「最短距離ルート」
* weight=”travel_time” → 「最短時間ルート」
* 他のエッジ属性を指定すれば、それに基づいた経路探索ができる

シェープファイルへの書き出し

シェープファイルへの書き出し

今回は作成した world を書き出す

これは GeodataFlameになっていれば簡単にできる

保存するときにはto_file を使う

シェープファイルで保存するなら

world.to_file(driver='ESRI Shapefile',filename='test_shape_world.shp')

というように
driverへ ESRI Shapefileを指定すればいい

でも

/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 328239523 of field POP_EST of feature 4 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 270625568 of field POP_EST of feature 8 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 144373535 of field POP_EST of feature 18 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 127575529 of field POP_EST of feature 27 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 211049527 of field POP_EST of feature 29 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 200963599 of field POP_EST of feature 56 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 1366417754 of field POP_EST of feature 98 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 163046161 of field POP_EST of feature 99 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 216565318 of field POP_EST of feature 102 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 1397715000 of field POP_EST of feature 139 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 108116615 of field POP_EST of feature 147 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 126264931 of field POP_EST of feature 155 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 100388073 of field POP_EST of feature 163 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(
/Users/snowpool/.pyenv/versions/3.11.0/lib/python3.11/site-packages/pyogrio/raw.py:723: RuntimeWarning: Value 112078730 of field POP_EST of feature 165 not successfully written. Possibly due to too larger number with respect to field width
  ogr_write(

となる

GPTによれば
併せて知っておくと良い Shapefile の制約
* フィールド名は最大 10 文字（長い列名は切られる）
* 文字列の最大長、数値の幅と小数点以下桁を持つ
* 64bit 整数を正式には持てない
* 複数のファイル（.shp/.shx/.dbf …）が必要

警告は Shapefile のフィールド幅不足。
ベストは GeoPackage/GeoJSON へ保存。
Shapefile が必須なら layer_options で幅指定を。

これは Shapefile（.shp）側の制約で数値フィールドの桁数が足りないために出る警告です。
Shapefile は属性を dBase(DBF) で持ちますが、各数値フィールドには「幅（桁数）と小数桁」を設定する必要があり、デフォルト推定だと POP_EST のような 10桁超え（例: 1,306,370,000） が収まりきらずに “not successfully written / too large for field width” になります。

Shapefile は古く、制限が多いです。GeoPackageやGeoJSONなら桁数制限に悩まされません。

このためGeoJSONでj保存する

world.to_file(driver='GeoJSON',filename='test_shape_world.geojson')

driverにGeoJSON　ファイル拡張子を .geojsonにすればOK

GeoDataFlame の集計

GeoDataFlame の集計

これは dissolve を使う

これを使い地理データの集計を行う

# world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowers'))
# cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
import geopandas as gpd
import requests, zipfile, io

# 世界地図（国境）
url_world = "https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_countries.zip"
r = requests.get(url_world)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall("naturalearth_data")
world = gpd.read_file("naturalearth_data/ne_110m_admin_0_countries.shp")

# 都市データ
url_cities = "https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_populated_places.zip"
r = requests.get(url_cities)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall("naturalearth_data")
cities = gpd.read_file("naturalearth_data/ne_110m_populated_places.shp")

print(world.head())
print(cities.head())

で世界地理情報の取得

国名はADMINになっているので
日本の情報を取得するのなら

world[world['ADMIN']=='Japan']

とすれば

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
155	Admin-0 country	1	2	Japan	JPN	0	2	Sovereign country	1	Japan	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((141.8846 39.18086, 140.95949 3...

world[world['ADMIN']=='Japan'].plot()

で日本地図の表示



この世界データは国ごとの情報になっている

これを国ではなく大陸ごとに変えたい
この時に位置情報などを取り扱うには dissolve を使う

その前に人口が出てないので調べた

「人口情報が出てこない」のは、表示されていないだけの可能性が高い
Natural Earth の属性は列数が多く、print(df.head()) だと途中で省略される

# 国ごとの推定人口
print(world[['ADMIN','ISO_A3','POP_EST']].head())

# 都市の人口（大きい順に上位10件）
print(cities[['NAME','ADM0NAME','POP_MAX','POP_MIN']].sort_values('POP_MAX', ascending=False).head(10))

とすると

                         ADMIN ISO_A3      POP_EST
0                         Fiji    FJI     889953.0
1  United Republic of Tanzania    TZA   58005463.0
2               Western Sahara    ESH     603253.0
3                       Canada    CAN   37589262.0
4     United States of America    USA  328239523.0
             NAME                  ADM0NAME   POP_MAX   POP_MIN
233         Tokyo                     Japan  35676000   8336599
218      New York  United States of America  19040000   8008278
224   Mexico City                    Mexico  19028000  10811002
234        Mumbai                     India  18978000  12691836
239     São Paulo                    Brazil  18845000  10021295
232      Shanghai                     China  14987000  14608512
237       Kolkata                     India  14787000   4631392
171         Dhaka                Bangladesh  12797394   7000940
210  Buenos Aires                 Argentina  12795000  10929146
216   Los Angeles  United States of America  12500000   3694820

というように表示された

なので

japan = world[world['ADMIN'] == 'Japan']
print(japan[['ADMIN','ISO_A3','POP_EST']])

とすれば

     ADMIN ISO_A3      POP_EST
155  Japan    JPN  126264931.0

というように表示された

Natural Earth の Admin-0 Countries データには continent というカラムがあり

continent の意味
* 各国・地域が属している 大陸名 を示す属性です。
* 値の例：
* “Asia”（日本、中国、インドなど）
* “Europe”（フランス、ドイツなど）
* “Africa”
* “North America”
* “South America”
* “Oceania”
* “Seven seas (open ocean)”（公海などの海洋エリア）
* “Antarctica”

ということで大陸名を調べようと

japan = world[world['ADMIN'] == 'Japan']
print(japan[['ADMIN','continent']])

としたが

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[53], line 2
      1 japan = world[world['ADMIN'] == 'Japan']
----> 2 print(japan[['ADMIN','continent']])

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: "['continent'] not in index"

となる
どうやら変更されている

print(world.columns.tolist())

の結果は

['featurecla', 'scalerank', 'LABELRANK', 'SOVEREIGNT', 'SOV_A3', 'ADM0_DIF', 'LEVEL', 'TYPE', 'TLC', 'ADMIN', 'ADM0_A3', 'GEOU_DIF', 'GEOUNIT', 'GU_A3', 'SU_DIF', 'SUBUNIT', 'SU_A3', 'BRK_DIFF', 'NAME', 'NAME_LONG', 'BRK_A3', 'BRK_NAME', 'BRK_GROUP', 'ABBREV', 'POSTAL', 'FORMAL_EN', 'FORMAL_FR', 'NAME_CIAWF', 'NOTE_ADM0', 'NOTE_BRK', 'NAME_SORT', 'NAME_ALT', 'MAPCOLOR7', 'MAPCOLOR8', 'MAPCOLOR9', 'MAPCOLOR13', 'POP_EST', 'POP_RANK', 'POP_YEAR', 'GDP_MD', 'GDP_YEAR', 'ECONOMY', 'INCOME_GRP', 'FIPS_10', 'ISO_A2', 'ISO_A2_EH', 'ISO_A3', 'ISO_A3_EH', 'ISO_N3', 'ISO_N3_EH', 'UN_A3', 'WB_A2', 'WB_A3', 'WOE_ID', 'WOE_ID_EH', 'WOE_NOTE', 'ADM0_ISO', 'ADM0_DIFF', 'ADM0_TLC', 'ADM0_A3_US', 'ADM0_A3_FR', 'ADM0_A3_RU', 'ADM0_A3_ES', 'ADM0_A3_CN', 'ADM0_A3_TW', 'ADM0_A3_IN', 'ADM0_A3_NP', 'ADM0_A3_PK', 'ADM0_A3_DE', 'ADM0_A3_GB', 'ADM0_A3_BR', 'ADM0_A3_IL', 'ADM0_A3_PS', 'ADM0_A3_SA', 'ADM0_A3_EG', 'ADM0_A3_MA', 'ADM0_A3_PT', 'ADM0_A3_AR', 'ADM0_A3_JP', 'ADM0_A3_KO', 'ADM0_A3_VN', 'ADM0_A3_TR', 'ADM0_A3_ID', 'ADM0_A3_PL', 'ADM0_A3_GR', 'ADM0_A3_IT', 'ADM0_A3_NL', 'ADM0_A3_SE', 'ADM0_A3_BD', 'ADM0_A3_UA', 'ADM0_A3_UN', 'ADM0_A3_WB', 'CONTINENT', 'REGION_UN', 'SUBREGION', 'REGION_WB', 'NAME_LEN', 'LONG_LEN', 'ABBREV_LEN', 'TINY', 'HOMEPART', 'MIN_ZOOM', 'MIN_LABEL', 'MAX_LABEL', 'LABEL_X', 'LABEL_Y', 'NE_ID', 'WIKIDATAID', 'NAME_AR', 'NAME_BN', 'NAME_DE', 'NAME_EN', 'NAME_ES', 'NAME_FA', 'NAME_FR', 'NAME_EL', 'NAME_HE', 'NAME_HI', 'NAME_HU', 'NAME_ID', 'NAME_IT', 'NAME_JA', 'NAME_KO', 'NAME_NL', 'NAME_PL', 'NAME_PT', 'NAME_RU', 'NAME_SV', 'NAME_TR', 'NAME_UK', 'NAME_UR', 'NAME_VI', 'NAME_ZH', 'NAME_ZHT', 'FCLASS_ISO', 'TLC_DIFF', 'FCLASS_TLC', 'FCLASS_US', 'FCLASS_FR', 'FCLASS_RU', 'FCLASS_ES', 'FCLASS_CN', 'FCLASS_TW', 'FCLASS_IN', 'FCLASS_NP', 'FCLASS_PK', 'FCLASS_DE', 'FCLASS_GB', 'FCLASS_BR', 'FCLASS_IL', 'FCLASS_PS', 'FCLASS_SA', 'FCLASS_EG', 'FCLASS_MA', 'FCLASS_PT', 'FCLASS_AR', 'FCLASS_JP', 'FCLASS_KO', 'FCLASS_VN', 'FCLASS_TR', 'FCLASS_ID', 'FCLASS_PL', 'FCLASS_GR', 'FCLASS_IT', 'FCLASS_NL', 'FCLASS_SE', 'FCLASS_BD', 'FCLASS_UA', 'geometry']
The history saving thread hit an unexpected error (OperationalError('unable to open database file')).History will not be written to the database.

カラム一覧を見ると、このデータでは
* 大陸は CONTINENT（全部大文字）
* 地域は REGION_UN、SUBREGION

これを元に

# 日本だけ絞り込み
japan = world[world["ADMIN"] == "Japan"]

# 大陸名を見る（大文字の CONTINENT）
print(japan[["ADMIN", "CONTINENT"]])

# 代替として地域区分を見る
print(japan[["ADMIN", "REGION_UN", "SUBREGION"]])

とすれば

  ADMIN CONTINENT
155  Japan      Asia
     ADMIN REGION_UN     SUBREGION
155  Japan      Asia  Eastern Asia

となる

これで世界の大陸ごとに分けることができそうなので
大陸列に対して dissolve で集計する

集計方法は
aggfuncで指定する
今回は sum 

continents = world.dissolve(by='CONTINENT',aggfunc='sum')
continents

とすると

geometry	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	ADM0_A3	GEOU_DIF	GEOUNIT	GU_A3	SU_DIF	SUBUNIT	SU_A3	BRK_DIFF	NAME	NAME_LONG	BRK_A3	BRK_NAME	BRK_GROUP	ABBREV	POSTAL	FORMAL_EN	FORMAL_FR	NAME_CIAWF	NOTE_ADM0	NOTE_BRK	NAME_SORT	NAME_ALT	MAPCOLOR7	MAPCOLOR8	MAPCOLOR9	MAPCOLOR13	POP_EST	POP_RANK	POP_YEAR	GDP_MD	GDP_YEAR	ECONOMY	INCOME_GRP	FIPS_10	ISO_A2	ISO_A2_EH	ISO_A3	ISO_A3_EH	ISO_N3	ISO_N3_EH	UN_A3	WB_A2	WB_A3	WOE_ID	WOE_ID_EH	WOE_NOTE	ADM0_ISO	ADM0_DIFF	ADM0_TLC	ADM0_A3_US	ADM0_A3_FR	ADM0_A3_RU	ADM0_A3_ES	ADM0_A3_CN	ADM0_A3_TW	ADM0_A3_IN	ADM0_A3_NP	ADM0_A3_PK	ADM0_A3_DE	ADM0_A3_GB	ADM0_A3_BR	ADM0_A3_IL	ADM0_A3_PS	ADM0_A3_SA	ADM0_A3_EG	ADM0_A3_MA	ADM0_A3_PT	ADM0_A3_AR	ADM0_A3_JP	ADM0_A3_KO	ADM0_A3_VN	ADM0_A3_TR	ADM0_A3_ID	ADM0_A3_PL	ADM0_A3_GR	ADM0_A3_IT	ADM0_A3_NL	ADM0_A3_SE	ADM0_A3_BD	ADM0_A3_UA	ADM0_A3_UN	ADM0_A3_WB	REGION_UN	SUBREGION	REGION_WB	NAME_LEN	LONG_LEN	ABBREV_LEN	TINY	HOMEPART	MIN_ZOOM	MIN_LABEL	MAX_LABEL	LABEL_X	LABEL_Y	NE_ID	WIKIDATAID	NAME_AR	NAME_BN	NAME_DE	NAME_EN	NAME_ES	NAME_FA	NAME_FR	NAME_EL	NAME_HE	NAME_HI	NAME_HU	NAME_ID	NAME_IT	NAME_JA	NAME_KO	NAME_NL	NAME_PL	NAME_PT	NAME_RU	NAME_SV	NAME_TR	NAME_UK	NAME_UR	NAME_VI	NAME_ZH	NAME_ZHT	FCLASS_ISO	TLC_DIFF	FCLASS_TLC	FCLASS_US	FCLASS_FR	FCLASS_RU	FCLASS_ES	FCLASS_CN	FCLASS_TW	FCLASS_IN	FCLASS_NP	FCLASS_PK	FCLASS_DE	FCLASS_GB	FCLASS_BR	FCLASS_IL	FCLASS_PS	FCLASS_SA	FCLASS_EG	FCLASS_MA	FCLASS_PT	FCLASS_AR	FCLASS_JP	FCLASS_KO	FCLASS_VN	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	
CONTINENT																																																																																																																																																																								
Africa	MULTIPOLYGON (((-11.43878 6.78592, -11.70819 6.8601, -12.4281 7.26294, -12.94905 7.79865, -13.12403 8.16395, -13.24655 8.90305, -13.68515 9.49474, -14.07404 9.88617, -14.33008 10.01572, -14.5797 10.21447, -14.69323 10.6563, -14.83955 10.87657, -15.13031 11.04041, -15.66418 11.45847, -16.08521 11.52459, -16.31479 11.80651, -16.30895 11.9587, -16.61384 12.17091, -16.67745 12.38485, -16.84152 13.15139, -16.71373 13.59496, -17.12611 14.37352, -17.62504 14.72954, -17.18517 14.91948, -16.70071 15.62153, -16.4631 16.13504, -16.54971 16.67389, -16.27055 17.16696, -16.14635 18.10848, -16.25688 19.09672, -16.37765 19.59382, -16.27784 20.09252, -16.53632 20.56787, -17.06342 20.99975, -17.02043 21.42231, -16.97325 21.88574, -16.58914 22.15823, -16.26192 22.67934, -16.32641 23.01777, -15.98261 23.72336, -15.426 24.35913, -15.08933 24.52026, -14.82465 25.10353, -14.80093 25.63626, -14.43994 26.25442, -13.7738 26.61889, -13.13994 27.64015, -13.12161 27.65415, -12.61884 28.03819, -11.68892 28.14864, -10.90096 28.83214, -10.39959 29.09859, -9.56481 29.93357, -9.81472 31.17774, -9.43479 32.0381, -9.30069 32.56468, -8.65748 33.24025, -7.65418 33.69706, -6.91254 34.11048, -6.24434 35.14587, -5.92999 35.75999, -5.19386 35.75518, -4.59101 35.33071, -3.64006 35.39986, -2.60431 35.17909, -2.16991 35.1684, -1.2086 35.71485, -0.12745 35.88866, 0.50388 36.30127, 1.46692 36.60565, 3.1617 36.7839, 4.81576 36.86504, 5.32012 36.71652, 6.26182 37.11066, 7.33038 37.11838, 7.73708 36.88571, 8.42096 36.94643, 9.50999 37.34999, 10.21 37.23, 10.18065 36.72404, 11.02887 37.0921, 11.10003 36.9, 10.6 36.41, 10.59329 35.94744, 10.93952 35.69898, 10.80785 34.83351, 10.14959 34.33077, 10.33966 33.78574, 10.85684 33.76874, 11.1085 33.29334, 11.48879 33.137, 12.66331 32.79278, 13.08326 32.87882, 13.91868 32.71196, 15.24563 32.26508, 15.71394 31.37626, 16.61162 31.18218, 18.02109 30.76357, 19.08641 30.26639, 19.57404 30.52582, 20.05335 30.98576, 19.82033 31.75179, 20.13397 32.2382, 20.85452 32.7068, 21.54298 32.8432, 22.89576 32.63858, 23.2368 32.19149, 23.60913 32.18726, 23.9275 32.01667, 24.92114 31.89936, 25.16482 31.56915, 26.49533 31.58568, 27.45762 31.32126, 28.45048 31.02577, 28.91353 30.87005, 29.68342 31.18686, 30.09503 31.4734, 30.97693 31.55586, 31.68796 31.4296, 31.96041 30.9336, 32.19247 31.26034, 32.99392 31.02407, 33.7734 30.96746, 34.26543 31.21936, 34.26544 31.21936, 34.82324 29.76108, 34.9226 29.50133, 34.64174 29.09942, 34.42655 28.34399, 34.15451 27.8233, 33.92136 27.6487, 33.58811 27.97136, 33.13676 28.41765, 32.42323 29.85108, 32.32046 29.76043, 32.73482 28.70523, 33.34876 27.69989, 34.10455 26.14227, 34.47387 25.59856, 34.79507 25.03375, 35.69241 23.92671, 35.49372 23.75237, 35.52598 23.10244, 36.69069 22.20485, 36.86623 22, 37.18872 21.01885, 36.96941 20.83744, 37.1147 19.80796, 37.48179 18.61409, 37.86276 18.36786, 38.41009 17.99831, 38.99062 16.84063, 39.26611 15.92272, 39.81429 15.43565, 41.17927 14.49108, 41.73495 13.92104, 42.27683 13.34399, 42.58958 13.00042, 43.08123 12.69964, 43.31785 12.39015, 43.28638 11.97493, 42.71587 11.73564, 43.1453 11.46204, 43.47066 11.27771, 43.66667 10.86417, 44.1178 10.44554, 44.61426 10.44221, 45.55694 10.69803, 46.6454 10.81655, 47.52566 11.12723, 48.0216 11.19306, 48.37878 11.37548, 48.94821 11.41062, 48.9482 11.41062, 49.26776 11.43033, 49.72862 11.5789, 50.25878 11.67957, 50.73202 12.0219, 51.1112 12.02464, 51.13387 11.74815, 51.04153 11.16651, 51.04531 10.6409, 50.83418 10.27972, 50.55239 9.19874, 50.07092 8.08173, 49.4527 6.80466, 48.59455 5.33911, 47.74079 4.2194, 46.56476 2.85529, 45.56399 2.04576, 44.06815 1.05283, 43.13597 0.2922, 42.04157 -0.91916, 41.81095 -1.44647, 41.58513 -1.68325, 40.88477 -2.08255, 40.63785 -2.49979, 40.26304 -2.57309, 40.12119 -3.27768, 39.80006 -3.68116, 39.60489 -4.34653, 39.20222 -4.67677, 38.74054 -5.90895, 38.79977 -6.47566, 39.44 -6.84, 39.47 -7.1, 39.19469 -7.7039, 39.25203 -8.00781, 39.18652 -8.48551, 39.53574 -9.11237, 39.9496 -10.0984, 40.31659 -10.3171, 40.31659 -10.3171, 40.47839 -10.76544, 40.43725 -11.76171, 40.56081 -12.63918, 40.59962 -14.20198, 40.77548 -14.69176, 40.47725 -15.40629, 40.08926 -16.10077, 39.45256 -16.72089, 38.53835 -17.10102, 37.41113 -17.58637, 36.28128 -18.65969, 35.8965 -18.84226, 35.1984 -19.55281, 34.78638 -19.78401, 34.70189 -20.49704, 35.17613 -21.25436, 35.37343 -21.84084, 35.38585 -22.14, 35.56255 -22.09, 35.53393 -23.07079, 35.37177 -23.53536, 35.60747 -23.70656, 35.45875 -24.12261, 35.04073 -24.47835, 34.21582 -24.81631, 33.01321 -25.35757, 32.57463 -25.72732, 32.66036 -26.14858, 32.91596 -26.21587, 32.83012 -26.74219, 32.58026 -27.47016, 32.46213 -28.30101, 32.20339 -28.7524, 31.521 -29.25739, 31.32556 -29.40198, 30.90176 -29.90996, 30.62281 -30.42378, 30.05572 -31.14027, 28.92555 -32.17204, 28.21976 -32.77195, 27.46461 -33.22696, 26.41945 -33.61495, 25.90966 -33.66704, 25.78063 -33.94465, 25.17286 -33.79685, 24.67785 -33.98718, 23.59404 -33.79447, 22.98819 -33.91643, 22.57416 -33.86408, 21.5428 -34.25884, 20.68905 -34.41718, 20.07126 -34.79514, 19.61641 -34.81917, 19.19328 -34.4626, 18.85531 -34.44431, 18.42464 -33.99787, 18.37741 -34.13652, 18.2445 -33.86775, 18.25008 -33.28143, 17.92519 -32.61129, 18.24791 -32.42913, 18.22176 -31.66163, 17.56692 -30.72572, 17.06442 -29.87864, 17.06292 -29.87595, 16.34498 -28.57671, 15.60182 -27.82125, 15.21047 -27.09096, 14.98971 -26.11737, 14.74321 -25.39292, 14.40814 -23.85301, 14.38572 -22.65665, 14.25771 -22.11121, 13.86864 -21.69904, 13.3525 -20.87283, 12.82685 -19.67317, 12.60856 -19.04535, 11.79492 -18.06913, 11.7342 -17.30189, 11.6401 -16.67314, 11.77854 -15.79382, 12.12358 -14.87832, 12.17562 -14.44914, 12.5001 -13.5477, 12.73848 -13.13791, 13.31291 -12.48363, 13.63372 -12.03864, 13.73873 -11.29786, 13.68638 -10.73108, 13.38733 -10.37358, 13.12099 -9.7669, 12.87537 -9.16693, 12.92906 -8.95909, 13.23643 -8.56263, 12.93304 -7.59654, 12.7283 -6.92712, 12.22735 -6.29445, 12.32243 -6.10009, 12.18234 -5.78993, 11.91496 -5.03799, 11.09377 -3.97883, 10.06614 -2.96948, 9.40525 -2.14431, 8.798 -1.1113, 8.83009 -0.77907, 9.04842 -0.45935, 9.29135 0.26867, 9.49289 1.01012, 9.30561 1.16091, 9.64916 2.28387, 9.7952 3.0734, 9.40437 3.73453, 8.94812 3.90413, 8.74492 4.35222, 8.48882 4.49562, 8.50029 4.77198, 7.46211 4.41211, 7.0826 4.46469, 6.69807 4.24059, 5.89817 4.26245, 5.3628 4.88797, 5.03357 5.6118, 4.32561 6.27065, 3.57418 6.2583, 2.6917 6.25882, 1.86524 6.14216, 1.06012 5.92884, -0.50764 5.34347, -1.06362 5.00055, -1.96471 4.71046, -2.85613 4.99448, -3.31108 4.9843, -4.00882 5.17981, -4.64992 5.16826, -5.8345 4.9937, -6.52877 4.70509, -7.51894 4.33829, -7.71216 4.36457, -7.97411 4.35576, -9.00479 4.83242, -9.91342 5.59356, -10.76538 6.14071, -11.43878 6.78592)), ((50.20027 -16.00026, 49.86061 -15.41425, 49.67261 -15.7102, 49.86334 -16.45104, 49.77456 -16.87504, 49.49861 -17.10604, 49.43562 -17.95306, 49.04179 -19.11878, 48.54854 -20.49689, 47.93075 -22.3915, 47.54772 -23.78196, 47.09576 -24.94163, 46.28248 -25.17846, 45.40951 -25.60143, 44.83357 -25.3461, 44.03972 -24.98835, 43.76377 -24.46068, 43.69778 -23.57412, 43.34565 -22.7769, 43.25419 -22.05741, 43.4333 -21.33648, 43.89368 -21.16331, 43.89637 -20.83046, 44.37433 -20.07237, 44.4644 -19.43545, 44.23242 -18.96199, 44.04298 -18.33139, 43.96308 -17.40994, 44.31247 -16.8505, 44.44652 -16.21622, 44.94494 -16.17937, 45.50273 -15.97437, 45.87299 -15.79345, 46.31224 -15.78002, 46.88218 -15.21018, 47.70513 -14.5943, 48.00521 -14.09123, 47.86905 -13.66387, 48.29383 -13.78407, 48.84506 -13.08917, 48.86351 -12.48787, 49.19465 -12.04056, 49.54352 -12.46983, 49.80898 -12.89528, 50.05651 -13.55576, 50.21743 -14.75879, 50.47654 -15.22651, 50.37711 -15.70607, 50.20027 -16.00026)))

というようになってしまうので方法を変える

continents = world.dissolve(
    by='CONTINENT',
    aggfunc={
        'POP_EST': 'sum',
        'GDP_MD': 'sum'
    }
).reset_index()

continents[['CONTINENT','POP_EST','GDP_MD']]

としたが

CONTINENT	POP_EST	GDP_MD	
0	Africa	1.306370e+09	2455514
1	Antarctica	4.490000e+03	898
2	Asia	4.550277e+09	32725478
3	Europe	7.454125e+08	21587850
4	North America	5.837560e+08	25075988
5	Oceania	4.120487e+07	1647113
6	Seven seas (open ocean)	1.400000e+02	16
7	South America	4.270667e+08	3852015

となって geometry が無くなっている

計後に列を選び直したときに geometry 列を外してしまい、ただの DataFrame になったからです。geometry が無いと GeoDataFrame ではなくなるので地理情報が消えたように見える

これだと使えないので別の方法にする

GPTで調べたら
pandas の表示設定で抑える（geometryはそのまま）
Shapelyオブジェクトは文字列化されると長い座標列が出ますが、
pd.set_option を調整すると途中で省略される

import pandas as pd
pd.set_option("display.max_colwidth", 80)  # 長さ制限
continents.head()

とすれば
これなら geometry は Shapely のまま保持され、
計算や可視化にもそのまま使える

CONTINENT	geometry	POP_EST	GDP_MD	
0	Africa	MULTIPOLYGON (((-11.43878 6.78592, -11.70819 6.8601, -12.4281 7.26294, -12.9...	1.306370e+09	2455514
1	Antarctica	MULTIPOLYGON (((-61.13898 -79.98137, -60.61012 -79.62868, -59.57209 -80.0401...	4.490000e+03	898
2	Asia	MULTIPOLYGON (((48.67923 14.0032, 48.23895 13.94809, 47.93891 14.00723, 47.3...	4.550277e+09	32725478
3	Europe	MULTIPOLYGON (((-53.55484 2.3349, -53.77852 2.3767, -54.08806 2.10556, -54.5...	7.454125e+08	21587850
4	North America	MULTIPOLYGON (((-155.22217 19.23972, -155.54211 19.08348, -155.68817 18.9161...	5.837560e+08	25075988

座標系の変換

座標系の変換

緯度経度は度になっているが
面積、距離などを求める時には変換が必要

https://sorabatake.jp/20670/
も参考にする

Name: WGS 84

– Lat[north]: Geodetic latitude (degree)
– Lon[east]: Geodetic longitude (degree)

がよく使われる
これはGPSで使われる

WGS 84
これを
EPSGにすると4326
https://epsg.io/4326
単位は度

これは世界

日本近辺に特化したものが
EPSG:6668
https://epsg.io/6668
単位は度

日本付近に絞るならこちらを使う方が良さそう

EPSG:6691
https://epsg.io/6691
これは
単位がメートルになっている
こちらも日本周辺がメイン

Area of use: Japan – between 138°E and 144°E, onshore and offshore.
となっているので
それ以外のエリアの精度が正しくないことがあるので注意

とりあえず品川の情報を表示

shinagawa
	N03_001	N03_002	N03_003	N03_004	N03_007	geometry
87	東京都	None	品川区	None	13109	POLYGON ((139.75501 35.61771, 139.75507 35.617...
88	東京都	None	品川区	None	13109	POLYGON ((139.77231 35.62188, 139.77302 35.621...
89	東京都	None	品川区	None	13109	POLYGON ((139.75953 35.625, 139.75966 35.6216,...
90	東京都	None	品川区	None	13109	POLYGON ((139.71943 35.64167, 139.71953 35.641...

座標情報を見るなら crs で表示できる

shinagawa.crs

結果は

<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

この中で注目するのは

EPSG:6668

と

- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)

から
単位は度であることがわかる

同様に路線図の crs も確認する

gdf.crs

結果は

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

こちらも

EPSG:4326

と

- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)

から単位は度であることがわかる

この２つはEPSGコードが異なるため
このままだと一緒に使えないので変換をして合わせる

今回は路線図のEPSGコードを6668に変換するので
to_crs を使う

gdf.to_crs(epsg=6668,inplace=True)
gdf.crs

これで

<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

となる

OSMxを使い最短経路の算出

OSMxを使い最短経路の算出

東京スカイツリーから
江戸東京博物館への最短経路を算出する

まずはスカイツリーと江戸東京博物館の座標をセットする

skytree = (35.7100,139.8180)
museum = (35.6963,139.7967)

次に一番近い node つまり交差点を取得

始点はスカイツリーなので

start = ox.get_nearest_node(G,skytree)

としたがエラー

--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[15], line 1 ----> 1 start = ox.get_nearest_node(G,skytree) AttributeError: module 'osmnx' has no attribute 'get_nearest_node'

これもバージョンの違いが原因

OSMnx 2.x では廃止され、ox.distance サブモジュールに移動
引数は (X=経度, Y=緯度) の順番で渡す必要があります
普通の Point(lon, lat) と同じ並び

start = ox.distance.nearest_nodes(G, X=skytree[1], Y=skytree[0])

とすればOK

同様に終点も

end = ox.distance.nearest_nodes(G,museum[1],museum[0])

次に最短経路

shortest_path = ox.shortest_path(G,start,end)

Foliumで表示できないので調べた

plot_route_folium
https://github.com/gboeing/osmnx/issues/1252
によれば

release (currently 1.9.4), then run your code. The latest v1 release will remain compatible with the v1 API, but any deprecated functionality you’re using will issue a FutureWarning with simple guidance to make it v2 compatible.
Yes, the folium module was deprecated with a warning for removal back in 2023, and finally removed in v2.0.0. The good news: interactive web mapping with folium remains easy to do with OSMnx in ~2 lines of code! See this example in the Examples Gallery for a demonstration.

となっていて
foliumモジュールは2023年に削除の警告とともに非推奨となり、最終的にv2.0.0で削除されました。

とのこと

https://github.com/gboeing/osmnx-examples/blob/main/notebooks/11-interactive-web-mapping.ipynb
のサンプルが参考になりそう

OSMnx 2.0.6 には Folium 専用関数は無い → 自分で PolyLine を追加すればOK

再度GPTでコード生成

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

というように

weight="length")

があるので調べた

✅ weight の意味
networkx.shortest_path の weight 引数は、**「どのエッジ属性を最短経路計算に使うか」**を指定するものです。
OSMnx で作成した G（道路グラフ）には、各エッジ（道路の区間）に様々な属性があります：
* “length” : 区間の長さ（メートル単位） ← デフォルトでよく使う
* “travel_time” : 区間の推定走行時間（秒単位）※ ox.add_edge_speeds / ox.add_edge_travel_times を実行した場合のみ
* その他、highway（道路種別）や lanes（車線数）など

✅ 例
1. 最短距離ルート（距離を最小化）

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

2. 最短時間ルート（所要時間を最小化）
まず速度と所要時間をエッジに付与してから：

import osmnx as ox

G = ox.add_edge_speeds(G)        # 各道路に推定速度 (km/h) を追加
G = ox.add_edge_travel_times(G)  # 各道路に travel_time (秒) を追加

route = nx.shortest_path(G, orig, dest, weight="travel_time")

✅ まとめ
* weight=”length” → 「最短距離ルート」
* weight=”travel_time” → 「最短時間ルート」
* 他のエッジ属性を指定すれば、それに基づいた経路探索ができる

OSMnxについて

OSMnxについて

OSMは
Open street map の略

pythonで道路ネットワーク分析

も参考にする

先にインストール

pip install osmnx

でインストールできる

import folium
import osmnx as ox
import matplotlib.pyplot as plt
from shapely.geometry import Point,LineString,Polygon

で必要なものをインポート

まずは道路データの取得
今回は取得するのは道路なので drive にする
歩道なら　walkなどを指定

Queryでどこの道路データを取るか指定
今回は
日本の東京都の墨田区にするので

query ="Sumidaku,Tokyo,Japan"
G = ox.graph_from_place(query,network_type="drive")

とする

print(G)で表示すると

MultiDiGraph with 3000 nodes and 7834 edges

となる

Nodes は交差点
Edges は道路の数
みたいな感じ

これを foliumで表示する

しかし

ox.plot_graph_folium(G)

だと

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[9], line 1
----> 1 ox.plot_graph_folium(G)

AttributeError: module 'osmnx' has no attribute 'plot_graph_folium'

✅ 背景
* 以前の OSMnx (v1.x くらいまで) には ox.plot_graph_folium(G) がありました。
* しかし OSMnx 2.0 以降では Folium 連携機能は 専用モジュールに移動しています。

from osmnx import folium as ox_folium

を追加
しかしエラーになる

ラッキーピエロ全店を最短でめぐりたい

では

ox.plot.plot_graph_folium

を使っている

ちなみにドキュメントは
https://osmnx.readthedocs.io/en/stable/index.html

調べてもバージョンは2.06で最新

GPTは2.1が最新と言っているが実際は2.06が最新

https://speakerdeck.com/mopinfish/osmnxniyorujie-lu-gou-zao-nofen-xi-toke-shi-hua?slide=20
を参考に行う

pip install contextily
でインストール

📌 まとめ
* contextily = 背景地図を GeoPandas/Matplotlib に簡単に追加するモジュール
* OpenStreetMap などのタイルを呼び出せる
* EPSG:3857（Web地図座標系）に自動変換してくれる

https://qiita.com/yutowac/items/b8e3d7653d3dcad16faf
の場合使っているバージョンが古い可能性がある

import osmnx as ox
import matplotlib.pyplot as plt
import contextily as cx

# 道路ネットワークの取得
place = "Higashi-ku, Sapporo, Hokkaido, Japan"
G = ox.graph_from_place(place, network_type="drive")

# グラフを GeoDataFrame に変換
nodes, edges = ox.graph_to_gdfs(G)

# 地図の作成
fig, ax = plt.subplots(figsize=(15, 15))

# エッジの描画
edges.plot(ax=ax, linewidth=0.5, edgecolor="dimgray")

# 背景地図の追加（OpenStreetMap）
cx.add_basemap(ax, crs=edges.crs.to_string(),
               source=cx.providers.OpenStreetMap.Mapnik)

# 地図の表示
plt.show()

で表示できたので
これを元にGPTでコード生成

import osmnx as ox
import matplotlib.pyplot as plt
import contextily as cx

# 墨田区の道路ネットワークを取得
query = "Sumida-ku, Tokyo, Japan"
G = ox.graph_from_place(query, network_type="drive")

# グラフを GeoDataFrame に変換
nodes, edges = ox.graph_to_gdfs(G)

# 描画
fig, ax = plt.subplots(figsize=(12, 12))

# 道路（edges）を描く
edges.plot(ax=ax, linewidth=0.7, edgecolor="blue")

# 背景地図を追加（OpenStreetMap）
cx.add_basemap(ax,
               crs=edges.crs.to_string(),
               source=cx.providers.OpenStreetMap.Mapnik)

# 表示
plt.show()

これで地図と情報が表示できた

あとは別の方法として
交差点を赤
道路を青
背景を白
で表示する

ox.plot_graph(G,bgcolor="white",node_color="red",edge_color="blue",node_size=2)

で表示はされる
ただし地図にはならない

座標系の変換

座標系の変換

緯度経度は度になっているが
面積、距離などを求める時には変換が必要

https://sorabatake.jp/20670/
も参考にする

Name: WGS 84

– Lat[north]: Geodetic latitude (degree)
– Lon[east]: Geodetic longitude (degree)

がよく使われる
これはGPSで使われる

WGS 84
これを
EPSGにすると4326
https://epsg.io/4326
単位は度

これは世界

日本近辺に特化したものが
EPSG:6668
https://epsg.io/6668
単位は度

日本付近に絞るならこちらを使う方が良さそう

EPSG:6691
https://epsg.io/6691
これは
単位がメートルになっている
こちらも日本周辺がメイン

Area of use: Japan – between 138°E and 144°E, onshore and offshore.
となっているので
それ以外のエリアの精度が正しくないことがあるので注意

とりあえず品川の情報を表示

shinagawa
	N03_001	N03_002	N03_003	N03_004	N03_007	geometry
87	東京都	None	品川区	None	13109	POLYGON ((139.75501 35.61771, 139.75507 35.617...
88	東京都	None	品川区	None	13109	POLYGON ((139.77231 35.62188, 139.77302 35.621...
89	東京都	None	品川区	None	13109	POLYGON ((139.75953 35.625, 139.75966 35.6216,...
90	東京都	None	品川区	None	13109	POLYGON ((139.71943 35.64167, 139.71953 35.641...

座標情報を見るなら crs で表示できる

shinagawa.crs

結果は

<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

この中で注目するのは

EPSG:6668

と

- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)

から
単位は度であることがわかる

同様に路線図の crs も確認する

gdf.crs

結果は

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

こちらも

EPSG:4326

と

- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)

から単位は度であることがわかる

この２つはEPSGコードが異なるため
このままだと一緒に使えないので変換をして合わせる

今回は路線図のEPSGコードを6668に変換するので
to_crs を使う

gdf.to_crs(epsg=6668,inplace=True)
gdf.crs

これで

<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

となる

行政データのプロット

行政データのプロット

今回は国土数値情報の行政区域データを使う

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html
これは市区町村の境界が表示できる

今回は東京都の平成２９年のものを使う
N03-170101_13_GML.zip
をダウンロード

ダブルクリックで解凍し

mv ~/Downloads/N03-170101_13_GML .

で移動

この中にあるシェープファイルを使う
N03-17_13_170101.shp
を使う

拡張子 .shpがシェープファイル

まずはこのファイルの読み込み

gyosei = gpd.read_file('./N03-170101_13_GML/N03-17_13_170101.shp')
gyosei

で

N03_001	N03_002	N03_003	N03_004	N03_007	geometry	
0	東京都	None	千代田区	None	13101	POLYGON ((139.77287 35.7037, 139.77279 35.7031...
1	東京都	None	中央区	None	13102	POLYGON ((139.78341 35.69645, 139.78459 35.696...
2	東京都	None	港区	None	13103	POLYGON ((139.77129 35.62841, 139.77128 35.628...
3	東京都	None	港区	None	13103	POLYGON ((139.76689 35.62774, 139.76718 35.627...
4	東京都	None	港区	None	13103	POLYGON ((139.77022 35.63199, 139.77046 35.631...
...	...	...	...	...	...	...
6205	東京都	None	None	所属未定地	None	POLYGON ((139.8413 35.64702, 139.84131 35.6471...
6206	東京都	None	None	所属未定地	None	POLYGON ((139.80438 35.60061, 139.80399 35.600...
6207	東京都	None	None	所属未定地	None	POLYGON ((139.81937 35.60899, 139.81923 35.608...
6208	東京都	None	None	所属未定地	None	POLYGON ((139.81009 35.61355, 139.81069 35.613...
6209	東京都	None	None	所属未定地	None	POLYGON ((139.82664 35.5976, 139.827 35.59703,...
6210 rows × 6 columns

となる

次にこの中から品川区のデータに絞るため
まずは区のみに絞る

そのためにまずはNoneのあるデータを除く

gyosei_tmp = gyosei.dropna(subset=['N03_003'])
gyosei_tmp

これで

N03_001	N03_002	N03_003	N03_004	N03_007	geometry	
0	東京都	None	千代田区	None	13101	POLYGON ((139.77287 35.7037, 139.77279 35.7031...
1	東京都	None	中央区	None	13102	POLYGON ((139.78341 35.69645, 139.78459 35.696...
2	東京都	None	港区	None	13103	POLYGON ((139.77129 35.62841, 139.77128 35.628...
3	東京都	None	港区	None	13103	POLYGON ((139.76689 35.62774, 139.76718 35.627...
4	東京都	None	港区	None	13103	POLYGON ((139.77022 35.63199, 139.77046 35.631...
...	...	...	...	...	...	...
140	東京都	None	西多摩郡	瑞穂町	13303	POLYGON ((139.35786 35.74464, 139.35776 35.744...
141	東京都	None	西多摩郡	瑞穂町	13303	POLYGON ((139.35 35.79414, 139.35016 35.79364,...
142	東京都	None	西多摩郡	日の出町	13305	POLYGON ((139.19302 35.78875, 139.19307 35.788...
143	東京都	None	西多摩郡	檜原村	13307	POLYGON ((139.10611 35.77867, 139.10644 35.778...
144	東京都	None	西多摩郡	奥多摩町	13308	POLYGON ((139.01912 35.89826, 139.01941 35.898...

となる

次にこのデータから区を含むものだけに絞り込む

gyosei_ku = gyosei_tmp['N03_003'].str.contains('区')
gyosei_ku

とすれば

0       True
1       True
2       True
3       True
4       True
       ...  
140    False
141    False
142    False
143    False
144    False
Name: N03_003, Length: 117, dtype: bool

となって判別できるので

この条件式を使って

gyosei_ku = gyosei_tmp[gyosei_tmp['N03_003'].str.contains('区')]
gyosei_ku

とすれば

N03_001	N03_002	N03_003	N03_004	N03_007	geometry	
0	東京都	None	千代田区	None	13101	POLYGON ((139.77287 35.7037, 139.77279 35.7031...
1	東京都	None	中央区	None	13102	POLYGON ((139.78341 35.69645, 139.78459 35.696...
2	東京都	None	港区	None	13103	POLYGON ((139.77129 35.62841, 139.77128 35.628...
3	東京都	None	港区	None	13103	POLYGON ((139.76689 35.62774, 139.76718 35.627...
4	東京都	None	港区	None	13103	POLYGON ((139.77022 35.63199, 139.77046 35.631...
...	...	...	...	...	...	...
107	東京都	None	葛飾区	None	13122	POLYGON ((139.87626 35.79479, 139.87661 35.793...
108	東京都	None	江戸川区	None	13123	POLYGON ((139.86285 35.63532, 139.86299 35.635...
109	東京都	None	江戸川区	None	13123	POLYGON ((139.8638 35.63722, 139.86391 35.6371...
110	東京都	None	江戸川区	None	13123	POLYGON ((139.8556 35.63856, 139.85563 35.6385...
111	東京都	None	江戸川区	None	13123	POLYGON ((139.89018 35.75055, 139.89044 35.750...
112 rows × 6 columns

というように区のみにできる

これで東京２３区の表示をする

gyosei_ku.plot()

これで区の区切りが白い線になっているのがわかる

次に品川区のみ色を変更する

ax = gyosei_ku.plot()
gdf.plot(ax=ax,color='orange',markersize=1)
ax = gyosei_ku.plot()
gdf.plot(ax=ax,color='orange',markersize=1)

次に品川のみに絞り込む

shinagawa = gyosei[gyosei['N03_003']=='品川区']
shinagawa

これで

N03_001	N03_002	N03_003	N03_004	N03_007	geometry	
87	東京都	None	品川区	None	13109	POLYGON ((139.75501 35.61771, 139.75507 35.617...
88	東京都	None	品川区	None	13109	POLYGON ((139.77231 35.62188, 139.77302 35.621...
89	東京都	None	品川区	None	13109	POLYGON ((139.75953 35.625, 139.75966 35.6216,...
90	東京都	None	品川区	None	13109	POLYGON ((139.71943 35.64167, 139.71953 35.641...

となる

次に品川区のみ表示する

ax = shinagawa.plot()
gdf.plot(ax=ax,color='orange',markersize=2)