Shapely中查找MultiPolygon内各Polygon最长内嵌线段的实现求助

问题描述 #

使用Shapely处理空间几何数据时，需查找MultiPolygon中每个Polygon的最长内嵌线段，已明确实现逻辑但代码运行异常，需技术修复。

实现逻辑 #

• 最长内嵌线段必经过多边形的两个顶点
• 遍历MultiPolygon中的每个Polygon
• 对每个Polygon，遍历其所有顶点对
• 为每对顶点创建无限延伸的直线
• 遍历Polygon的所有边，检查直线与边是否相交
• 将所有交点存入intersection_points数组
• 遍历intersection_points中的点对，生成线段并检查是否完全在多边形内，更新最大长度

伪代码 #

for polygon in multipolygon {
    for point1 in polygon {
        for point2 in polygon {
            if point1 == point2 {
                continue
            }
            line = makeLine(point1,point2)
            intersection_points = []
            for edge in polygon {
                if line_intersects_edge(line,edge) {
                    intersection_points.insert( intersection_point(line,edge) )
                }
            }
        }
    }
}

max_len = 0
for intersection_pt1 in intersection_points {
    for intersection_pt2 in intersection_points {
        if intersection_pt1 != intersection_pt2 {
            line_segment = make_line_segment(intersection_pt1,intersection_pt2)
            if line_segment.lies_within(polygon) {
                max_len = max( max_len, line_segment.length() )
            }
        }
    }
}

待修复Python代码 #

from shapely.geometry import MultiPolygon, Polygon, LineString, Point

def line_intersects_edge(line, edge):
    # Implement a function to check if a line intersects with an edge
    return line.intersects(edge)

def intersection_point(line, edge):
    # Implement a function to find the intersection point between a line and an edge
    return line.intersection(edge)

def make_line_segment(point1, point2):
    # Implement a function to create a line segment from two points
    return LineString([point1, point2])

multipolygon = MultiPolygon(
    [
        Polygon([(7, 10), (8, 11), (9, 11), (8, 10), (7, 9.5), (7, 10)]),
        Polygon([(9.5, 8.5), (10, 9), (10, 10), (11, 9), (9.5, 8.5)]),
    ]
)

max_len = 0

# Iterate over each polygon in the multipolygon
for polygon in multipolygon.geoms:
    # Iterate over each point in the polygon
    for point1 in polygon.exterior.coords:
        # Iterate over each point in the polygon again
        for point2 in polygon.exterior.coords:
            # Skip if points are the same
            if point1 == point2:
                continue
        
            # Create a line from the two points
            line = LineString([Point(point1), Point(point2)])
        
            # Find intersection points
            intersection_points = []
            for edge in polygon.exterior.coords[:-1]:
                if line_intersects_edge(line, LineString([edge, polygon.exterior.coords[polygon.exterior.coords.index(edge) + 1]])):
                   
                  intersection_points.append(intersection_point(line, LineString([edge, polygon.exterior.coords[polygon.exterior.coords.index(edge) + 1]])))
        
            # Iterate over intersection points
            for intersection_pt1 in intersection_points:
                for intersection_pt2 in intersection_points:
                    if intersection_pt1 != intersection_pt2:
                        # Create a line segment
                        line_segment = make_line_segment(intersection_pt1, intersection_pt2)
                        # Check if line segment lies within the polygon
                        if line_segment.within(polygon):
                            max_len = max(max_len, line_segment.length)

print("Max length:", max_len)

问题分析 #

1. 无限直线创建错误：原代码用LineString([Point(point1), Point(point2)])创建的是线段，并非无限直线，无法正确计算与多边形边的所有交点。
2. 边遍历方式错误：使用polygon.exterior.coords.index(edge)获取下一个顶点，若多边形存在重复顶点会导致索引错误，且效率低下。
3. 交点处理不严谨：未过滤非Point类型的交点（如直线与边重合时返回LineString），也未去重，导致后续计算重复或出错。
4. 顶点对重复遍历：遍历所有point1和point2的组合（包括顺序相反的对），造成不必要的重复计算。
5. 线段存在性判断问题：within方法不包含边界，若最长线段沿多边形边界会被忽略，需根据需求调整为包含边界的判断。

修复方案 #

1. 模拟无限直线：基于顶点对创建超出多边形包围盒的线段，模拟无限直线（因为直线与多边形边的交点仅会出现在多边形范围内）。
2. 正确遍历多边形边：用zip遍历连续坐标对，直接获取每条边的两个端点。
3. 过滤并去重交点：仅保留Point类型的交点，通过坐标元组去重。
4. 优化顶点对遍历：仅遍历索引i<j的顶点对，避免重复计算。
5. 调整线段判断逻辑：使用covered_by方法包含边界情况，或根据需求选择within。

修复后的代码 #

from shapely.geometry import MultiPolygon, Polygon, LineString, Point

def create_infinite_line(point1, point2, polygon):
    # 基于多边形包围盒创建足够长的线段，模拟无限直线
    min_x, min_y, max_x, max_y = polygon.bounds
    dx = point2[0] - point1[0]
    dy = point2[1] - point1[1]
    
    # 延长线段至包围盒外10倍距离，确保覆盖所有可能交点
    scale = max(max_x - min_x, max_y - min_y) * 10
    if dx == 0 and dy == 0:
        return LineString([point1, point2])
    
    ext_point1 = (point1[0] - dx * scale, point1[1] - dy * scale)
    ext_point2 = (point1[0] + dx * scale, point1[1] + dy * scale)
    return LineString([ext_point1, ext_point2])

def get_intersection_points(infinite_line, polygon):
    intersections = []
    coords = list(polygon.exterior.coords)
    # 遍历所有边
    for i in range(len(coords)-1):
        edge = LineString([coords[i], coords[i+1]])
        if infinite_line.intersects(edge):
            inter = infinite_line.intersection(edge)
            # 仅保留Point类型的交点
            if isinstance(inter, Point):
                # 保留6位小数避免浮点精度问题，转换为元组去重
                inter_coord = (round(inter.x, 6), round(inter.y, 6))
                intersections.append(inter_coord)
    # 去重后返回
    return list(set(intersections))

def make_line_segment(coord1, coord2):
    return LineString([coord1, coord2])

multipolygon = MultiPolygon(
    [
        Polygon([(7, 10), (8, 11), (9, 11), (8, 10), (7, 9.5), (7, 10)]),
        Polygon([(9.5, 8.5), (10, 9), (10, 10), (11, 9), (9.5, 8.5)]),
    ]
)

max_len = 0

for polygon in multipolygon.geoms:
    coords = list(polygon.exterior.coords[:-1])  # 去掉最后一个重复的闭合顶点
    num_points = len(coords)
    # 遍历i<j的顶点对，避免重复计算
    for i in range(num_points):
        point1 = coords[i]
        for j in range(i+1, num_points):
            point2 = coords[j]
            # 创建模拟无限直线
            infinite_line = create_infinite_line(point1, point2, polygon)
            # 获取所有有效交点
            intersection_coords = get_intersection_points(infinite_line, polygon)
            # 遍历交点对计算最长内嵌线段
            num_inter = len(intersection_coords)
            for m in range(num_inter):
                coord_m = intersection_coords[m]
                for n in range(m+1, num_inter):
                    coord_n = intersection_coords[n]
                    segment = make_line_segment(coord_m, coord_n)
                    # 检查线段是否在多边形内（包含边界）
                    if segment.covered_by(polygon):
                        current_len = segment.length
                        if current_len > max_len:
                            max_len = current_len

print("Max length:", round(max_len, 6))

验证结果 #

运行修复后的代码，输出结果为Max length: 2.5，对应第一个多边形中(7, 9.5)到(9, 11)的线段长度，符合预期。

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