1

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?

enter image description here

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.

2
  • Can you add a screenshot showing what you want to do? (for example the before and after from QGIS)
    – BERA
    Aug 4, 2020 at 17:48
  • @BERA just added- thanks for tuning in. I should have added end result screenshots.
    – Barnard87
    Aug 4, 2020 at 19:11

1 Answer 1

2

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:

Polygons before transformation Polygons after transformation

2
  • Yes thank you! Everything seems to have erased perfect. I'll use this to erase one GeoSeries from the next as well, it will take a little playing around. Onelast thing, do you know of a way to retain the polygon's attributes? Or will I need to reappend those after?
    – Barnard87
    Aug 5, 2020 at 13:07
  • 1
    ^^ Edit to add onto the above content. Instead of reading in as a GeoSeries, I just read in the shapefile as a GeoDataFrame, looked into its geometries by calling poly.geometry where you had "poly", and re-wrote the file back out. Loading it in as a GeoSeries only reads in its geometries
    – Barnard87
    Aug 5, 2020 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.