CODEWITHSUNDEEP
pythonmatplotlibshapefileshapelyfiona
anyryg
anyryg

Reputation: 313

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

Upvotes: 3

Views: 2247

Answers (1)

Can H. Tartanoglu
Can H. Tartanoglu

Reputation: 378

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

Upvotes: 2

Related Questions