圆与街道多边形交集结果缺失问题排查求助

我使用多边形代表街道（蓝色部分），还有一个基于点生成的圆形缓冲区（橙色部分），希望找出圆内的所有街道片段。

当前实现代码如下：

from myproject import geometry as geo
...
# 加载城市街道数据并构建空间索引树
# 城市街道由多个多边形组成
city = geo.load_geojson("my_city.geojson")
tree = shapely.STRtree(city)

# 生成圆形缓冲区
center = shapely.Point(x, y)
circle = geo.generate_circle_around_point(center, radius)

# 查询与圆形相交的街道并计算交集
shapes = [city[idx].intersection(circle) for idx in tree.query(circle)]

但这段代码得到的交集结果质量不佳，存在部分本该在圆内的街道片段缺失的情况：

补充说明：以下是生成橙色圆形的代码（因使用GPS坐标，需将半径从米转换为度）：

def generate_circle_around_point(point: shapely.Point, radius: int) -> shapely.Polygon:
    # 定义本地投影坐标系（UTM）
    local_projection = pyproj.Proj(proj="utm", zone=18, ellps="WGS84")

    # 定义WGS84转UTM的转换函数
    to_local_projection = partial(
        pyproj.transform, pyproj.CRS("EPSG:4326"), local_projection
    )

    # 将中心点转换到UTM坐标系
    center_point_local = shapely.ops.transform(to_local_projection, point)

    # 创建指定半径的圆形缓冲区
    circle_buffer = center_point_local.buffer(radius, resolution=200)

    # 定义UTM转WGS84的转换函数
    to_gps_projection = partial(
        pyproj.transform, local_projection, pyproj.CRS("EPSG:4326")
    )
    circle_buffer_gps = shapely.ops.transform(to_gps_projection, circle_buffer)

    return circle_buffer_gps

可能的原因及解决办法 #

• 空间索引的候选集漏判
  STRtree的query方法依赖空间索引做快速筛选，存在极小概率误判：比如完全被圆包含的街道多边形，可能因索引边界计算误差被排除在候选集外。
  解决方式：对索引返回的候选集做二次精确校验，只保留真正与圆形相交或被包含的街道：

candidate_indices = tree.query(circle)
valid_shapes = []
for idx in candidate_indices:
    street = city[idx]
    if street.intersects(circle) or circle.contains(street):
        valid_shapes.append(street.intersection(circle))

• 投影转换的精度与匹配问题

  1. 固定UTM带不匹配：代码中硬编码了zone=18，如果你的城市不在UTM 18带的经度范围（约-78°到-72°）内，投影误差会导致圆形形状严重失真，进而影响交集计算。解决：根据中心点经度动态计算UTM带：

     def get_utm_zone(lon):
         return int((lon + 180) / 6) + 1
     
     utm_zone = get_utm_zone(center.x)
     local_projection = pyproj.Proj(proj="utm", zone=utm_zone, ellps="WGS84")
     

  2. 多次投影转换的精度损失：先转UTM生成缓冲区再转回WGS84，两次转换会累积精度误差。解决：将街道数据也转换到UTM坐标系下完成所有空间计算，最后再统一转回WGS84：

     # 将所有街道转换到UTM
     city_utm = [shapely.ops.transform(to_local_projection, geom) for geom in city]
     tree_utm = shapely.STRtree(city_utm)
     # 直接在UTM下计算交集
     shapes_utm = [city_utm[idx].intersection(circle_buffer) for idx in tree_utm.query(circle_buffer)]
     # 转回WGS84
     shapes = [shapely.ops.transform(to_gps_projection, geom) for geom in shapes_utm]
     

• 缓冲区的折线近似误差
  圆形缓冲区在UTM下由大量折线顶点模拟，投影回WGS84后，折线间隙可能导致部分街道边缘无法被正确检测到交集。解决：适当提高buffer方法的resolution参数（比如设为300），或者对生成的圆形做微小缓冲消除间隙：

circle_buffer_gps = shapely.ops.transform(to_gps_projection, circle_buffer).buffer(0.00001)

• 街道多边形的拓扑错误
  如果原始街道多边形存在自相交、无效拓扑等问题，Shapely的交集计算会出现异常结果。解决：加载数据时先校验并修复几何体：

city = []
for geom in geo.load_geojson("my_city.geojson"):
    if not geom.is_valid:
        # 通过buffer(0)修复无效几何体
        geom = geom.buffer(0)
    city.append(geom)

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