用 OSM 打造城市空間指標：從歷史快照抓取到 500m 緩衝區七大指標計算

在項目開始之初爲什麽選 OSM 而不是政府開放資料

做城市量化分析時，第一個問題永遠是「資料從哪來」。選項通常有三類：

項目選擇 OSM 的理由很簡單：做跨年比較、加上跨城市複用，OSM 是唯一能在「歷史一致性 × schema 統一」兩個維度同時滿足的來源。

只是代價是要面對「社群貢獻偏差」（後面 A3 篇會專門講怎麼修正），但這個問題是可量化、可修正的。

從 Overpass 拉歷史 OSM 的關鍵：date 過濾器

OSM 整個資料庫是時間序列，而每個 element 都有 timestamp，Overpass 支援用 date 過濾器查「某個時間點的快照」：

# dataloader.py:56-64
date_str = f"{year}-01-01T00:00:00Z"
ox.settings.overpass_settings = f'[out:json][timeout:180][date:"{date_str}"]'

gdf = ox.features_from_place(city, tags)
file_path = f"{SAVE_DIR}/osm_{year}.geojson"
gdf.to_file(file_path, driver='GeoJSON')

NOTE：

1. [date:"2018-01-01T00:00:00Z"]：設定 Overpass 獲取 2018/1/1 OSM 快照」。
2. ox.features_from_place(city, tags) ：osmnx 將 city name 解析成 polygon，再用該 polygon 去 Overpass 抓 elements。
3. 獲取后存成 .geojson 本地快取，方便後續繼續處理。這邊需要注意以下幾點

• Overpass 公開端點有 query quota，重複跑很容易被擋。
• 一次抓 8 年 × 5 種 tag 的 query，沒快取每次重跑要 2 小時起跳。
• 後續所有指標計算都從本地 geojson 讀，跟網路完全脫鉤。

設計本地 geojson 快取的雙路徑架構

項目中 OSM 模組有兩層 cache：

# osm/fetch_osm.py:44-69
# 1. 優先從本地 osm_{year}.geojson 讀取並過濾
main_geojson_path = save_dir / f"osm_{year}.geojson"
if main_geojson_path.exists() and main_geojson_path.stat().st_size > 0:
    full_gdf = gpd.read_file(main_geojson_path)
    filtered_gdf = full_gdf.copy()
    for key, value in tags.items():
        if value is True:
            filtered_gdf = filtered_gdf[filtered_gdf[key].notna()]
        else:
            filtered_gdf = filtered_gdf[filtered_gdf[key] == value]
    return filtered_gdf

# 2. 從快取文件讀取
cache_path = save_dir / f"osm_{label}_{year}.geojson"
if cache_path.exists():
    return gpd.read_file(cache_path)

# 3. 本地沒有，返回 None（不打網路）

設計雙層緩存 cache 主要基於檔案處理時間跟空間考量，如下：：

• 主檔 osm_<year>.geojson 一次下完所有 tag（建築 + 設施 + 商店 + 休閒），體積大但只下載一次。後續按 tag 過濾在記憶體裡做，幾毫秒就完成。
• 快取檔 osm_<label>_<year>.geojson 是針對單一 tag 的精細下載，給某些主檔太大、想精細控制的場景用。
• 第三條路徑直接回 None 而不打網路 —— 這是工程上很重要的決定。讓「沒資料」變成可預期的訊號，而不是讓 pipeline 卡在 retry。

OSM 數據處理 - 500m 緩衝區 + 空間索引

抓到資料後，每個採樣點要算它周圍 500m 範圍的指標：

# osm/fetch_osm.py:29
BUFFER_DEG = 0.005  # ~500m

# osm/fetch_osm.py:123-132
def _query_nearby(gdf, area):
    """用空間索引快速查詢 buffer 範圍內的要素"""
    if gdf is None or len(gdf) == 0:
        return gpd.GeoDataFrame()
    # 用 sindex 先粗篩 bounding box，再精確 intersects
    candidates_idx = list(gdf.sindex.intersection(area.bounds))
    if not candidates_idx:
        return gpd.GeoDataFrame()
    candidates = gdf.iloc[candidates_idx]
    return candidates[candidates.geometry.intersects(area)]

其中有兩個細節：

1. BUFFER_DEG = 0.005

OSM 預設座標系是 EPSG:4326（經緯度）。在台中緯度 24°，經度 0.005° ≈ 510m，緯度 0.005° ≈ 555m，平均 ~500m。

其中為什麼不直接用公尺？ 則是因查詢階段要跟原始 EPSG:4326 的 geometry 比對，先 buffer 在度單位、再精算才不會在每個點都做一次投影轉換，每個點的投影轉化成本說大不大，但如果 500m 換成小單位則數量又會變大。

2. gdf.sindex.intersection() 的兩階段過濾

geopandas 的 sindex 是 R-tree。空間查詢的標準套路：

• 先用 R-tree 找 bounding box 相交的候選（O(log n)）
• 再對候選做 geometry.intersects() 精算（候選數量已經很少）

直接 gdf[gdf.geometry.intersects(area)] 是 O(n)，台中 2025 年 5 萬棟建築一次掃完要好幾秒；所以在這邊利用空間索引將處理時間降到毫秒級。

程式七指標計算

把上面的零件組起來，單一採樣點 500m 範圍的七項指標：

# osm/fetch_osm.py:135-184 (節錄)
def compute_osm_metrics_for_point(lon, lat, buildings_gdf, amenity_gdf, shop_gdf, leisure_gdf, roads_gdf):
    area = Point(lon, lat).buffer(BUFFER_DEG)
    result = {"lon": lon, "lat": lat}

    # 建築物
    nearby = _query_nearby(buildings_gdf, area)
    result["building_count"] = len(nearby)
    if len(nearby) > 0:
        projected = nearby.to_crs(epsg=3826)
        result["building_area_mean"] = float(projected.geometry.area.mean())
    else:
        result["building_area_mean"] = 0.0

    # Amenity / Shop / Leisure
    result["amenity_count"] = len(_query_nearby(amenity_gdf, area))
    result["shop_count"]    = len(_query_nearby(shop_gdf, area))
    result["leisure_count"] = len(_query_nearby(leisure_gdf, area))

    # POI 多樣性 (amenity + shop + leisure 的 Shannon entropy)
    type_counts = [result["amenity_count"], result["shop_count"], result["leisure_count"]]
    result["poi_diversity"] = shannon_entropy(type_counts)

    # 路網總長度
    nearby_r = _query_nearby(roads_gdf, area)
    if len(nearby_r) > 0:
        projected_r = nearby_r.to_crs(epsg=3826)
        result["road_length_total"] = float(projected_r.geometry.length.sum())
    else:
        result["road_length_total"] = 0.0

    return result

poi_diversity 的 Shannon entropy 邏輯後續會單獨展開，這裡先當黑盒：值愈大 = 功能愈混合 = 單一功能區的程度愈低。