pythonmatplotlibshapefileshapelyfiona

Edit polygon coords using Python, Shapely and Fiona


I need to edit the geometry of intersecting polygons and I don't know how I can save modified geometry to a shapefile. Is it even possible?

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

Solution

  • The removal of intersections between many polygons is most likely a complex problem. Moreover, I used your method as the solver in my solution.

    Answer

    The answer to your question, is yes. You can rectify the intersections between the polygons in your shp file; however, you need to create new Polygon objects, you can't just change the exterior coordinates of an existing Polygon.

    Store metadata and disc from original shp file

    The solution below writes/creates the resulting polygon set to a new shp file. This requires us to store the metadata from the original shp file, and pass it to the new one. We also need to store the properties of each polygon, I store these in a separate list, set_of_properties.

    No need for two for loops

    You don't need to for loops, just use combinations from the itertools standard library to loop through all possible polygon combinations. I use index combinations to replace polygons that are intersecting with new ones.

    Outer do...while-loop

    In very cringe caes, a rectification using your method may actually introduce new intersections. We can catch these and rectify them by looping through your solver until there are no intersections left. This requires a do... while loop, but there is no do...while loop in Python. Moreover, it can be implemented with while-loops (see Solution for implementation).

    Solution

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