基于OSMnx批量计算海量Quadkey 14瓦片街道密度的技术咨询

一、先下载Quadkey 9大图再提取子瓦片的方案完全可行 #

这种"一次性下载大区域路网+本地切割子图"的思路，能直接解决你原方案中500万次请求的效率问题，同时完全契合你想用官方basic_stats()的需求。

二、具体操作步骤 #

1. 把Quadkey转换为地理边界 #

用quadkey库将Quadkey编码转成经纬度范围，9级或14级瓦片都适用：

import quadkey

# 转换单个Quadkey 9为边界（示例）
qk9 = quadkey.QuadKey("你的Quadkey9编码")
bbox9 = qk9.to_geo()  # 返回格式：(min_lon, min_lat, max_lon, max_lat)

# 批量转换所有Quadkey 14的边界
all_qk14_bboxes = []
for qk_str in 你的Quadkey14列表:
    qk = quadkey.QuadKey(qk_str)
    all_qk14_bboxes.append((qk_str, qk.to_geo()))

2. 下载Quadkey 9范围的路网并提前投影 #

一次性下载大区域路网，然后投影到墨卡托CRS（和basic_stats()默认使用的坐标系一致），避免后续重复投影浪费时间：

import osmnx as ox

# 下载路网，按需设置network_type（比如drive/walk/bike）
G = ox.graph_from_bbox(bbox9[3], bbox9[1], bbox9[2], bbox9[0], network_type="drive")
# 投影到墨卡托CRS
G_proj = ox.project_graph(G)

3. 批量裁剪子图并计算密度 #

遍历所有Quadkey 14的边界，从已投影的大图中裁剪出对应子图，直接调用basic_stats()：

stats_results = []
for qk_str, bbox14 in all_qk14_bboxes:
    # 裁剪子图
    G_sub = ox.truncate.truncate_graph_bbox(G_proj, bbox14[3], bbox14[1], bbox14[2], bbox14[0])
    # 计算统计值，空图会返回默认值，可按需过滤
    sub_stats = ox.basic_stats(G_sub)
    sub_stats["quadkey14"] = qk_str
    stats_results.append(sub_stats)

三、关于OSMnx多请求的问题 #

这个方案只需要发起1次（或少量几次，如果你需要覆盖多个Quadkey9区域）下载请求，后续全是本地操作，完全不需要处理多请求的频率限制。如果确实需要下载多个大区域，提前配置OSMnx的请求参数即可：

ox.config(requests_timeout=120, retry_delay=5, max_retries=10)

四、解释GeoDataFrame方法的差异问题 #

你之前遇到的length列和墨卡托长度不符的问题，是因为OSMnx边的length字段是WGS84球面距离，而墨卡托下的几何长度是平面距离，两者本身就有差异。basic_stats()里的street_density_km是用投影后的平面长度计算的，所以如果想用GeoDataFrame方法对齐结果，必须统一用墨卡托的平面长度和面积：

import geopandas as gpd
from shapely.geometry import box

# 把大图转成边的GeoDataFrame（已投影）
edges_proj = ox.graph_to_gdfs(G_proj, nodes=False, edges=True)

for qk_str, bbox14 in all_qk14_bboxes:
    # 把子瓦片边界转成GeoDataFrame并投影到同一CRS
    bbox_poly = box(bbox14[0], bbox14[1], bbox14[2], bbox14[3])
    bbox_gdf = gpd.GeoDataFrame({"geometry": [bbox_poly]}, crs="EPSG:4326").to_crs(G_proj.graph["crs"])
    # 筛选和子瓦片相交的边
    edges_sub = edges_proj[edges_proj.intersects(bbox_gdf.iloc[0]["geometry"])]
    # 计算总长度（转千米）和面积（转平方千米）
    total_length_km = edges_sub["geometry"].length.sum() / 1000
    area_km2 = bbox_gdf["geometry"].area.iloc[0] / 10**6
    # 密度计算
    density = total_length_km / area_km2

不过还是推荐你用子图裁剪+basic_stats()的方式，完全贴合官方实现，不用手动处理细节，误差更小。

内容的提问来源于stack exchange，提问作者Luca Ceribelli