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()
The removal of intersections between many polygons is most likely a complex problem. Moreover, I used your method as the solver in my solution.
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.
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
.
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.
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).
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
})