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

Question

我需要编辑相交多边形的几何形状，但不知道如何将修改后的几何形状保存到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()

- anyryg

---

你能同时上传 polygon23.shp 文件吗？ - SKPS

https://drive.google.com/file/d/1MoaXBPbXZOQUMYdAhsSb7QXZhAZudLip/view?usp=sharing - anyryg

你打开的文件从未关闭，这是不好的实践，使用with语句（最简单的方法）。for循环逻辑是多余的。你没有保存任何暂存在coords_polygonB中的交集矩形修正。你可以在我的答案中看到我对你的代码所做的重要更改。 - Can H. Tartanoglu

Answer 1

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

答案

回答您的问题，是的。您可以纠正shp文件中多边形之间的交集；但是，您需要创建新的Polygon对象，而不能仅更改现有Polygon的外部坐标。

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

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

不需要两个for循环

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

外部do...while循环

在非常糟糕的情况下，使用您的方法进行矫正实际上可能会引入新的交点。我们可以通过循环运行您的求解器来捕获这些并进行纠正，直到没有交点为止。这需要一个do...while循环，但是Python中没有do...while循环。此外，它可以使用while循环实现（请参见解决方案以获取实现方式）。

解决方案

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

- Can H. Tartanoglu