对两个GeoDataFrame执行Shapely difference操作的最快方法是什么？

你当前循环计算耗时极长的核心原因是：遍历17000个建筑逐个执行difference操作，每次差集运算都要对当前已有的大几何对象做全量计算，时间复杂度极高，完全可以通过批量合并的方式将计算次数从17000次降到1次，实现量级级的性能提升。

1. 核心优化逻辑 #

你现有代码中已经用到了合并几何的接口，完全不需要循环逐个做差，直接把所有建筑几何批量合并为单个对象，再和城市范围几何做一次差集即可，结果和循环计算的结果完全一致：

注意：旧版本的cascaded_union已被弃用，推荐用性能更优的unary_union替代

buildings = gi_osm_geo.overpass_result_to_geodataframe(
        gi_osm_geo.query_overpass_types_in_area(
            region_or_tags=area_tags,
            search_tags=[("building", "")],
            link_by_or=True,
            relations=True,
            ways=True,
            nodes=False,
        )
    )

in_crs = buildings.crs
if in_crs is None:
    in_crs = "EPSG:4326"  # If no crs is set, assume EPSG:4326
    buildings = buildings.set_crs(in_crs)
utm = buildings.estimate_utm_crs()
buildings = buildings.to_crs(utm)

# 调整顺序：先过滤有效几何再做缓冲，避免无效几何生成错误结果
buildings = buildings[buildings["geometry"].is_valid]
buffered_geo = buildings.geometry.buffer(buffer_in_m)
# 批量合并所有建筑的缓冲几何，只需执行一次
buildings_union = buffered_geo.unary_union

2. 单步计算可用面积 #

假设存储城市范围的GeoDataFrame名为city_boundary，先转成和建筑相同的UTM投影，单次差集即可得到结果：

# 城市范围转同投影
city_utm = city_boundary.to_crs(utm)
# 仅执行一次差集，替代17000次循环运算
available_geo = city_utm.geometry.iloc[0].difference(buildings_union)
# 直接获取可用面积
available_area = available_geo.area

额外性能提升建议 #

• geopandas 0.10+版本默认集成pygeos矢量运算后端，相比原生shapely运算性能可提升3~10倍，建议升级到最新版本
• 如果存在城市范围外的建筑数据，可以先通过空间连接gpd.sjoin过滤后再做合并，进一步减少计算量

内容的提问来源于stack exchange，提问作者Siddhant Gupta