gpt4 book ai didi

python - 使用 Python、Shapely 和 Fiona 编辑多边形坐标

转载 作者:行者123 更新时间:2023-12-04 15:02:09 25 4
gpt4 key购买 nike

我需要编辑相交多边形的几何形状,但我不知道如何将修改后的几何形状保存到 shapefile。有可能吗?

from shapely.geometry import Polygon, shape
import matplotlib.pyplot as plt
import fiona


c = fiona.open('polygon23.shp', 'r')
d = fiona.open('polygon23.shp', 'r')

for poly in c.values():
for poly2 in d.values():
p_poly = shape(poly['geometry'])
p_poly2 = shape(poly2['geometry'])
intersect_polygons = p_poly.intersection(p_poly2)
if type(intersect_polygons) == Polygon:
intersect_polygons = p_poly.intersection(p_poly2).exterior.coords
if p_poly.exterior.xy != p_poly2.exterior.xy:

y_difference = abs(intersect_polygons[0][1]) - abs(intersect_polygons[2][1])

coords_polygonB = p_poly2.exterior.coords[:]

coords_polygonB[0] = (coords_polygonB[0][0], coords_polygonB[0][1] + (y_difference))
coords_polygonB[1] = (coords_polygonB[1][0], coords_polygonB[1][1] + (y_difference))
coords_polygonB[2] = (coords_polygonB[2][0], coords_polygonB[2][1] + (y_difference))
coords_polygonB[3] = (coords_polygonB[3][0], coords_polygonB[3][1] + (y_difference))
coords_polygonB[4] = (coords_polygonB[4][0], coords_polygonB[4][1] + (y_difference))

p_poly2 = Polygon(coords_polygonB)

x,y = p_poly.exterior.xy
plt.plot(x,y)
x,y = p_poly2.exterior.xy
plt.plot(x,y)
plt.show()

最佳答案

移除许多多边形之间的交点很可能是一个复杂的问题。此外,我在我的解决方案中使用了您的方法作为求解器。

回答

您的问题的答案是肯定的。您可以纠正 shp 文件中多边形之间的交点;但是,您需要创建新的多边形对象,不能只更改现有多边形的外部坐标。

存储原始 shp 文件的元数据和光盘

下面的解决方案将生成的多边形集写入/创建到一个新的 shp 文件中。这需要我们存储原始 shp 文件中的元数据,并将其传递给新文件。我们还需要存储每个多边形的属性,我将它们存储在一个单独的列表中,set_of_properties

不需要两个for循环

您不需要 for 循环,只需使用 itertools 标准库中的组合来遍历所有可能的多边形组合。我使用索引组合来替换与新多边形相交的多边形。

外层 do...while 循环

在非常棘手的情况下,使用您的方法进行的整改实际上可能会引入新的交叉点。我们可以捕获这些并通过循环求解器来纠正它们,直到没有交叉点为止。这需要一个 do...while 循环,但是 Python 中没有 do...while 循环。此外,它可以用 while-loops 实现(参见实现方案)。

解决方案

from itertools import combinations

from shapely.geometry import Polygon, Point, shape, mapping
import matplotlib.pyplot as plt
import fiona

SHOW_NEW_POLYGONS = False

polygons, set_of_properties = [], []
with fiona.open("polygon23.shp", "r") as source:
for line in source:
polygons.append(shape(line["geometry"]))
set_of_properties.append(line["properties"])
meta = source.meta

poly_index_combinations = combinations(tuple(range(len(polygons))), 2)
while True:
intersection_record = []
for i_poly_a, i_poly_b in poly_index_combinations:
poly_a, poly_b = polygons[i_poly_a], polygons[i_poly_b]
if poly_a.exterior.xy == poly_b.exterior.xy:
# print(f"The polygons have identical exterior coordinates:\n{poly_a} and {poly_b}\n")
continue

intersecting = poly_a.intersection(poly_b)
if type(intersecting) != Polygon:
continue
intersecting_polygons = intersecting.exterior.coords
if not intersecting_polygons:
# print(f"No intersections between\n{poly_a} and {poly_b}\n")
continue

print("Rectifying intersection")
y_diff = abs(intersecting_polygons[0][1]) - abs(intersecting_polygons[2][1])

new_poly_b = Polygon((
Point(float(poly_b.exterior.coords[0][0]), float(poly_b.exterior.coords[0][1] + y_diff)),
Point(float(poly_b.exterior.coords[1][0]), float(poly_b.exterior.coords[1][1] + y_diff)),
Point(float(poly_b.exterior.coords[2][0]), float(poly_b.exterior.coords[2][1] + y_diff)),
Point(float(poly_b.exterior.coords[3][0]), float(poly_b.exterior.coords[3][1] + y_diff)),
Point(float(poly_b.exterior.coords[4][0]), float(poly_b.exterior.coords[4][1] + y_diff))
))

if SHOW_NEW_POLYGONS:
x, y = poly_a.exterior.xy
plt.plot(x, y)
x, y = new_poly_b.exterior.xy
plt.plot(x, y)
plt.show()

polygons[i_poly_b] = new_poly_b
intersection_record.append(True)

if not intersection_record:
break

with fiona.open("new.shp", "w", **meta) as sink:
for poly, properties in zip(polygons, set_of_properties):
sink.write({
"geometry": mapping(poly),
"properties": properties
})

关于python - 使用 Python、Shapely 和 Fiona 编辑多边形坐标,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/66784929/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com