我有一个 Multipolygon，我想在 shapely 中找到 Multipolygon 中每个多边形中最长的线。 如何在 Multipolygon 中找到位于多边形内的最大直线？
编辑： 我明白了逻辑，但执行起来遇到困难。
polygon中的每个
multipolygonpolygonvertices的所有对
points(
polygonverticeslinevertexlineedgeslineedgepointsintersection_points[]pointsintersection_points[]linemax_length图像插图可在参考链接中找到 - ref1、ref2、ref3
伪代码-
max_len = 0
for polygon in multipolygon{
for point1 in polygon{
for point2 in polygon{
if point1 == point2{
continue
}
else{
line = makeLine(point1,point2)
intersection_points = []
for edge in polygon{
if line_intersects_edge(line,edge) == true{
intersection_points.insert( intersection_point(line,edge))
}
}
}
}
}
}
for intersection_pt1 in intersection_points{
for intersection_pt2 in intersection_points{
if intersection_pt1 != intersection_pt2{
line_segement = 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)
我认为与您描述的算法相比，您的脚本中有两个错误：
尚未彻底测试，但这解决了我认为的这两个问题：
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