Geographic Information Systems Asked by Barnard87 on August 9, 2020
I have a Polygon of transects that contain overlapping features where they converge. I need to erase where they overlap, but still keep one of the polygons there (preferably the preceding one).
Essentially, I want to iterate through each feature, have the current feature (f1) see where it overlaps with all other features (for loop, f2) erase that overlap from f2 and assign its new geometry as that.
I’ve used this code within the QGIS Python Console which works, but I’m developing a notebook to automate the generation of transcripts, and this is the final task. Here is the code I use within Q:
for f1 in layer.getFeatures():
for f2 in layer2.getFeatures():
if f1.id() < f2.id():
geom1 = f1.geometry()
geom2 = f2.geometry()
new_geom = geom2.difference(geom1)
layer2.dataProvider().changeGeometryValues({f2.id(): new_geom})
I found a chunk of code which gets me a shapefile of every overlapping area, now I need to erase that geometry, but only from the succeeding feature (f2) not f1. Here is the code.
g1 = gpd.GeoDataFrame.from_file(linesBufferFolder + "Lines_Buffer")
g1.shape
import itertools
geoms = g1['geometry'].tolist()
intersection_iter = gpd.GeoDataFrame(gpd.GeoSeries([poly[0].intersection(poly[1]) for poly in itertools.combinations(geoms, 2) if poly[0].intersects(poly[1])]), columns=['geometry'])
intersection_iter.to_file(linesBufferFolder + "Lines_int_iter.shp")
I’m assuming I need something in that for loop to tell it to erase while it is there?
Above is the end result. The highlighted polygon comes after the one overlapping it on its right, therefore the polygon to the right remains the same, and the highlighted on is now the shape you see. If I highlight the one to the left of the currently highlighted, it would also show no overlap. All of these polygons are square 1m buffers.
You need to store the intermediary results as well, there must be some feedback between the iteration steps.
from itertools import combinations
import geopandas as gpd
from shapely.geometry import LineString
# Some example data
poly = gpd.GeoSeries([
LineString([(4, 4), (9, 9)]).buffer(1, cap_style=3),
LineString([(0, 0), (3, 3)]).buffer(1, cap_style=3),
LineString([(0, 0), (-3, 1)]).buffer(1, cap_style=3),
])
for p1_idx, p2_idx in combinations(poly.index, 2):
if poly.loc[p1_idx].intersects(poly.loc[p2_idx]):
# Store intermediary results back to poly
poly.loc[p2_idx] -= poly.loc[p1_idx]
poly.plot(alpha=0.75, cmap="tab10")
Before ⇒ After:
Correct answer by mwil.me on August 9, 2020
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP