拆分自相交多边形仅在 Shapely 中返回一个多边形
Splitting self-intersecting polygon only returned one polygon in Shapely
我在 Windows 7 64 位中使用 Python 3.5 64 位,匀称版本 1.5.13。
我有以下代码返回了一个自相交的多边形:
import numpy as np
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt
x = np.array([ 0.38517325, 0.40859912, 0.43296919, 0.4583215 , 0.4583215 ,
0.43296919, 0.40859912, 0.38517325, 0.36265506, 0.34100929])
y = np.array([ 62.5 , 56.17977528, 39.39698492, 0. ,
0. , 17.34605377, 39.13341671, 60.4180932 ,
76.02574417, 85.47008547])
polygon = Polygon(np.c_[x, y])
plt.plot(*polygon.exterior.xy)
这是正确的。然后我尝试使用 buffer(0)
:
获取两个单独的多边形
split_polygon = polygon.buffer(0)
plt.plot(*polygon.exterior.xy)
print(type(split_polygon))
plt.fill(*split_polygon.exterior.xy)
不幸的是,它只返回了两个多边形:
有人可以帮忙吗?谢谢!
第一步是关闭 LineString 以制作 LinearRing,这就是多边形的组成部分。
from shapely.geometry import LineString, MultiPolygon
from shapely.ops import polygonize, unary_union
# original data
ls = LineString(np.c_[x, y])
# closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple # False
但是请注意,它并不简单,因为线条交叉形成领结。 (根据我的经验,广泛使用的 buffer(0) 技巧通常不适用于固定领结)。这不适合 LinearRing,因此需要进一步的工作。使用 unary_union
:
使其变得简单和 MultiLineString
mls = unary_union(lr)
mls.geom_type # MultiLineString'
然后使用polygonize
从线条中找到多边形:
for polygon in polygonize(mls):
print(polygon)
或者,如果您想要一个多边形几何体:
mp = MultiPolygon(list(polygonize(mls)))
我还在2020年纠结了一段时间,终于写了一个清理self intersections的方法。
这需要 Shapely v 1.2.1 explain_validity() 方法才能工作。
def clean_bowtie_geom(base_linearring):
base_polygon = Polygon(base_linearring)
invalidity = explain_validity(base_polygon)
invalid_regex = re.compile('^(Self-intersection)[[](.+)\s(.+)[]]$')
match = invalid_regex.match(invalidity)
if match:
groups = match.groups()
intersect_point = (float(groups[1]), float(groups[2]))
new_linring_coords1 = []
new_linring_coords2 = []
pop_new_linring = False
for i in range(0, len(base_linearring.coords)):
if i == len(base_linearring.coords) - 1:
end_point = base_linearring.coords[0]
else:
end_point = base_linearring.coords[i + 1]
start_point = base_linearring.coords[i]
if not pop_new_linring:
if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
new_linring_coords2.append(intersect_point)
new_linring_coords1.append(intersect_point)
pop_new_linring = True
else:
new_linring_coords1.append(start_point)
else:
new_linring_coords2.append(start_point)
if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
new_linring_coords2.append(intersect_point)
pop_new_linring = False
corrected_linear_ring1 = LinearRing(coordinates=new_linring_coords1)
corrected_linear_ring2 = LinearRing(coordinates=new_linring_coords2)
polygon1 = Polygon(corrected_linear_ring1)
polygon2 = Polygon(corrected_linear_ring2)
def is_point_on_line_and_between(start, end, pt, tol=0.0005):
"""
Checks to see if pt is directly in line and between start and end coords
:param start: list or tuple of x, y coordinates of start point of line
:param end: list or tuple of x, y coordinates of end point of line
:param pt: list or tuple of x, y coordinates of point to check if it is on the line
:param tol: Tolerance for checking if point on line
:return: True if on the line, False if not on the line
"""
v1 = (end[0] - start[0], end[1] - start[1])
v2 = (pt[0] - start[0], pt[1] - start[1])
cross = cross_product(v1, v2)
if cross <= tol:
# The point lays on the line, but need to check if in between
if ((start[0] <= pt[0] <= end[0]) or (start[0] >= pt[0] >= end[0])) and ((start[1] <= pt[1] <= end[1]) or (start[1] >= pt[1] >= end[1])):
return True
return False
这不是最干净的代码,但它帮我完成了工作。
输入是一个具有自相交几何形状的 LinearRing (is_simple=False),输出可以是 2 个 LinearRing 或两个多边形,无论你喜欢哪个(或者有条件选择一个,世界是你的牡蛎,真的)。
编辑
Shapely 1.8.0,新增功能。
shapely.validation.make_valid() 将采用一个自相交多边形和 return 一个多多边形,每个多边形都是通过在自相交点处分割创建的。
我在 Windows 7 64 位中使用 Python 3.5 64 位,匀称版本 1.5.13。
我有以下代码返回了一个自相交的多边形:
import numpy as np
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt
x = np.array([ 0.38517325, 0.40859912, 0.43296919, 0.4583215 , 0.4583215 ,
0.43296919, 0.40859912, 0.38517325, 0.36265506, 0.34100929])
y = np.array([ 62.5 , 56.17977528, 39.39698492, 0. ,
0. , 17.34605377, 39.13341671, 60.4180932 ,
76.02574417, 85.47008547])
polygon = Polygon(np.c_[x, y])
plt.plot(*polygon.exterior.xy)
这是正确的。然后我尝试使用 buffer(0)
:
split_polygon = polygon.buffer(0)
plt.plot(*polygon.exterior.xy)
print(type(split_polygon))
plt.fill(*split_polygon.exterior.xy)
不幸的是,它只返回了两个多边形:
有人可以帮忙吗?谢谢!
第一步是关闭 LineString 以制作 LinearRing,这就是多边形的组成部分。
from shapely.geometry import LineString, MultiPolygon
from shapely.ops import polygonize, unary_union
# original data
ls = LineString(np.c_[x, y])
# closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple # False
但是请注意,它并不简单,因为线条交叉形成领结。 (根据我的经验,广泛使用的 buffer(0) 技巧通常不适用于固定领结)。这不适合 LinearRing,因此需要进一步的工作。使用 unary_union
:
mls = unary_union(lr)
mls.geom_type # MultiLineString'
然后使用polygonize
从线条中找到多边形:
for polygon in polygonize(mls):
print(polygon)
或者,如果您想要一个多边形几何体:
mp = MultiPolygon(list(polygonize(mls)))
我还在2020年纠结了一段时间,终于写了一个清理self intersections的方法。
这需要 Shapely v 1.2.1 explain_validity() 方法才能工作。
def clean_bowtie_geom(base_linearring):
base_polygon = Polygon(base_linearring)
invalidity = explain_validity(base_polygon)
invalid_regex = re.compile('^(Self-intersection)[[](.+)\s(.+)[]]$')
match = invalid_regex.match(invalidity)
if match:
groups = match.groups()
intersect_point = (float(groups[1]), float(groups[2]))
new_linring_coords1 = []
new_linring_coords2 = []
pop_new_linring = False
for i in range(0, len(base_linearring.coords)):
if i == len(base_linearring.coords) - 1:
end_point = base_linearring.coords[0]
else:
end_point = base_linearring.coords[i + 1]
start_point = base_linearring.coords[i]
if not pop_new_linring:
if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
new_linring_coords2.append(intersect_point)
new_linring_coords1.append(intersect_point)
pop_new_linring = True
else:
new_linring_coords1.append(start_point)
else:
new_linring_coords2.append(start_point)
if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
new_linring_coords2.append(intersect_point)
pop_new_linring = False
corrected_linear_ring1 = LinearRing(coordinates=new_linring_coords1)
corrected_linear_ring2 = LinearRing(coordinates=new_linring_coords2)
polygon1 = Polygon(corrected_linear_ring1)
polygon2 = Polygon(corrected_linear_ring2)
def is_point_on_line_and_between(start, end, pt, tol=0.0005):
"""
Checks to see if pt is directly in line and between start and end coords
:param start: list or tuple of x, y coordinates of start point of line
:param end: list or tuple of x, y coordinates of end point of line
:param pt: list or tuple of x, y coordinates of point to check if it is on the line
:param tol: Tolerance for checking if point on line
:return: True if on the line, False if not on the line
"""
v1 = (end[0] - start[0], end[1] - start[1])
v2 = (pt[0] - start[0], pt[1] - start[1])
cross = cross_product(v1, v2)
if cross <= tol:
# The point lays on the line, but need to check if in between
if ((start[0] <= pt[0] <= end[0]) or (start[0] >= pt[0] >= end[0])) and ((start[1] <= pt[1] <= end[1]) or (start[1] >= pt[1] >= end[1])):
return True
return False
这不是最干净的代码,但它帮我完成了工作。
输入是一个具有自相交几何形状的 LinearRing (is_simple=False),输出可以是 2 个 LinearRing 或两个多边形,无论你喜欢哪个(或者有条件选择一个,世界是你的牡蛎,真的)。
编辑
Shapely 1.8.0,新增功能。 shapely.validation.make_valid() 将采用一个自相交多边形和 return 一个多多边形,每个多边形都是通过在自相交点处分割创建的。