Longest line lying completely inside a polygon in shapely

Question

I have a Multipolygon and I want to find the longest line in each polygon in Multipolygon in shapely. How do I find the largest line that can lie within a polygon in Multipolygon?

edit: I have understood the logic but facing difficulty in implementation.

1. The longest line must pass through two vertices of the polygon( reference - ref1, ref2, ref3
2. We iterate though each polygon in multipolygon
3. for each polygon we iterate through all pair of vertices (points) of the polygon
4. for each pair of vertices, we find the line( extending both sides indefinitely ) joining the vertex pairs
5. for each line formed, we traverse all edges( line segments ) and check if the line intersects the edge
6. we store all intersection points in an array intersection_points[]
7. for each pair of points in intersection_points[] we form a line segment and check if it lies entirely within the polygon, if it does and its length is grater than max_length then we update max_length

Image illustrations is available in the reference links - ref1, ref2, ref3

Pseudocode -

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() )
            }
        }
    }
}

The python code is not working though, any help would be highly appreciated!

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

This code does the job. I have tested it with a few complex test cases and it is giving valid answer

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)

Answer 2

I think there are 2 errors in your script compared to the algorithm you describe:

• the candidate lines are not elongated "to infinity" to find the outer edges of the polygon, so the intersection points found aren't all the ones that can lead to the longest line
• the loop over the intersection points found is at the wrong level of the loops

Not thoroughly tested, but this solves both issues I think:

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