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.
polygon in multipolygonpolygon we iterate through all pair of vertices (points) of the polygonvertices, we find the line( extending both sides indefinitely ) joining the vertex pairsline formed, we traverse all edges( line segments ) and check if the line intersects the edgepoints in an array intersection_points[]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_lengthImage 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)
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)