Shapely中完全位于多边形内的最长线段

Question

我有一个多边形组合（Multipolygon），我想找到每个多边形中最长的线段，使用shapely这个库。请问如何找到可以完全位于多边形内的最长线段呢？

补充说明：我已经理解了逻辑，但在实现上遇到了一些困难。

1. 最长的线段必须经过多边形的两个顶点（参考 - ref1, ref2, ref3）
2. 我们要遍历每一个polygon（多边形）在multipolygon（多边形组合）中
3. 对于每个polygon，我们要遍历所有的vertices（顶点）对
4. 对于每一对vertices，我们找到一条连接这两个顶点的line（线），这条线在两边无限延伸
5. 对于每条形成的line，我们遍历所有的edges（边），检查这条线是否与边相交
6. 我们把所有的交点存储在一个数组intersection_points[]中
7. 对于intersection_points[]中的每一对points（点），我们形成一段line，并检查它是否完全位于多边形内。如果是，并且它的长度大于max_length，那么我们就更新max_length

参考链接中有图片说明 - ref1, ref2, ref3

伪代码 -

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)

Answer 1

我觉得你的脚本和你描述的算法相比，有两个错误：

• 候选线没有被延伸到“无穷远”，这样就找不到多边形的外边缘，所以找到的交点并不是所有能形成最长线的交点。
• 对找到的交点进行循环的层级放错了地方。

虽然没有经过彻底测试，但我觉得这个可以解决这两个问题：

from typing import Tuple
import shapely
from shapely import LineString, box, Point

def line_interpolate_to_bbox(
    p1: Tuple[float, float],
    p2: Tuple[float, float],
    bbox: Tuple[float, float, float, float],
) -> LineString:
    """
    Interpolates a line so both points are onto the boundaries of a given bounding box.

    Args:
        p1 (Tuple[float, float]): The first point of the line.
        p2 (Tuple[float, float]): The second point of the line.
        bbox (Tuple[float, float, float, float]): A tuple representing the boundary of
            the bounding box to interpolate to in the format (minx, miny, maxx, maxy).

    Returns:
        LineString: The interpolated line on the boundary of the bbox.
    """
    minx, miny, maxx, maxy = bbox
    if p1[0] == p2[0]:  # vertical line
        return LineString([(p1[0], miny), (p1[0], maxy)])
    elif p1[1] == p2[1]:  # horizontal line
        return LineString([(minx, p1[1]), (maxx, p1[1])])
    else:
        # linear equation: y = k*x + m
        k = (p2[1] - p1[1]) / (p2[0] - p1[0])
        m = p1[1] - k * p1[0]
        y0 = k * minx + m
        y1 = k * maxx + m
        x0 = (miny - m) / k
        x1 = (maxy - m) / k
        points_on_boundary_lines = [
            Point(minx, y0),
            Point(maxx, y1),
            Point(x0, miny),
            Point(x1, maxy),
        ]
        bbox = box(minx, miny, maxx, maxy)
        points_sorted_by_distance = sorted(points_on_boundary_lines, key=bbox.distance)
        return LineString(points_sorted_by_distance[:2])

def longest_line(poly):
    max_len = 0

    # Iterate over each polygon in the multipolygon
    intersection_points = []
    max_len = 0
    for polygon in poly.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

                # Extend the line to cover the mbr of the polygon
                line = pygeoops.line_interpolate_to_bbox(point1, point2, polygon.bounds)

                # Find intersection points between the line and the polygon
                intersection = LineString(polygon.boundary).intersection(line)
                intersection_points.extend(shapely.get_coordinates(intersection))

        # Iterate over intersection points
        for intersection_pt1 in intersection_points:
            for intersection_pt2 in intersection_points:
                if all(intersection_pt1 != intersection_pt2):
                    # Create a line segment
                    line_segment = LineString([intersection_pt1, intersection_pt2])
                    # Check if line segment lies within the polygon
                    if line_segment.within(polygon):
                        length = line_segment.length
                        if length > max_len:
                            max_len = length
                            max_line = line_segment

    return max_line

Answer 2

这段代码能完成任务。我用几个复杂的测试案例试过，结果都很有效。

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


def get_infinite_line(p1: Point, p2: Point, scale=10000) -> LineString:
  """
  Creates a LineString extending indefinitely in the direction of p2 from p1.

  Args:
      p1: The first Point object defining the line (Shapely Point).
      p2: The second Point object defining the line (Shapely Point).
      scale: A factor to scale the extension length (default: 10).

  Returns:
      A LineString object representing the infinite line.
  """


  # Get the vector representing the direction using x and y components
  direction = (p2[0] - p1[0], p2[1] - p1[1])

  # Extend the line in both directions by scaling the direction vector
  extended_p1 = Point(p1[0] - direction[0] * scale, p1[1] - direction[1] * scale)
  extended_p2 = Point(p2[0] + direction[0] * scale, p2[1] + direction[1] * scale)

  # Create the LineString object
  return LineString([extended_p1, p2, extended_p2])


# 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)]),
#     ]
# )

multipolygon = MultiPolygon(
    [
        Polygon([(9, 10.5),(8.5, 10.5),(8.5, 11),(8,10.5),(7.5, 9),(8, 9.5),(8.5, 9),(9, 9.5),(8.5, 10),(7.8, 9.8),(9, 10.5)])
    
    ]
)



max_len = 0
line_ans = None
# Assuming multipolygon is a Shapely MultiPolygon object
for polygon in multipolygon.geoms:
    for point1 in polygon.exterior.coords:
        for point2 in polygon.exterior.coords:
            if point1 == point2:
                continue

            line = get_infinite_line(point1, point2)
            intersection_points = [point1,point2]

            for i, edge in enumerate(polygon.exterior.coords[:-1]):
                next_edge_index = (i + 1) % len(polygon.exterior.coords)
                edge_line = sg.LineString([edge, polygon.exterior.coords[next_edge_index]])

                if line.intersects(edge_line):
                    intersection_point = line.intersection(edge_line)
                    intersection_points.append(intersection_point.coords[0])  # Access coordinates of intersection

            for intersection_pt1 in intersection_points:
                for intersection_pt2 in intersection_points:
                    if intersection_pt1 != intersection_pt2:
                        line_segment = sg.LineString([intersection_pt1, intersection_pt2])
                        if polygon.contains(line_segment) and max_len < line_segment.length:
                            max_len = max(max_len, line_segment.length)
                            line_ans = line_segment

print(max_len)
print(line_ans)