Remove the intersecting parts in a set of polygons [closed]

Question

I have a set of (~200k) polygons showing deforestation in the Amazon. Each polygon is associated with a date when the deforestation happened.

Some of them are intersecting with each other, ie. one location is covered by multiple polygons. This is an error of the dataset - one place cannot be deforested twice (within a few years' time).

Since it is not feasible to manually check and resolve each of these situations, I want to get rid of these areas. So, for every polygon of the dataset. I want to remove/cut all of its parts where it intersects with other datasets.

Example:

On the picture, the purple and yellow rectangles represent two polygons that overlap. I'd like to cut out this overlapping part from both of the polygons.

I am preferably looking for a solution in Python. I have tried but failed to come up with a reasonable solution. Since this seems like a common spatial problem I imagine there should be some alternate solutions available.

Answer 1

Using GeoPandas (as it was mentioned in a comment). Inspired by this question.

import geopandas as gpd
import itertools

# load data
gpkg = r'D:\OneDrive\geodata\misc.gpkg'
layer = r'overlapping_polygons'
gdf1 = gpd.read_file(gpkg, layer=layer)

# get list of geometries
geoms = gdf1['geometry'].tolist()

# iterate over all combinations of polygons and get the intersections (overlaps)
overlaps = gpd.GeoDataFrame(gpd.GeoSeries([poly[0].intersection(poly[1]) for poly in itertools.combinations(geoms, 2) if poly[0].intersects(poly[1])]), columns=['geometry'])

# set the crs
overlaps.crs = gdf1.crs

# erase the overlaps from the original geodataframe
gdf2 = gdf1.overlay(overlaps, how='difference')

# plot before and after
gdf1.plot(alpha=0.75, cmap="tab10")
gdf2.plot(alpha=0.75, cmap="tab10")

Answer 2

A sort of brute-force way to do this is to put together a list of all groups of intersecting polygons, then for each group compute the difference between every polygon and every other polygon in the group. Obviously this is imaginary code and it's O(N^2) but I think it would work.

intsxns_all = []
test_polys = list(polygons)
for poly in polygons:
    intsxns = [poly]
    if poly not in test_polys:
        continue
    else:
        test_polys.remove(poly)

    tested = []
    for test_poly in test_polys:
        if test_poly.intersects(poly):
            intsxns.append(test_poly)
            tested.append(test_poly)

    for test in tested:
        test_polys.remove(test)

    intsxns_all.append(intsxns)

# let's use numbers as stand-ins for geometries:
#   intsxns_all maybe looks like [[1], [2, 3], [4, 5, 6]]

output = []
for intsxn in intsxns_all:
    for poly in intsxn:
        diff = poly

        subtract_polys = [p for p in intsxn if not p.equals(poly)]
        for subtract_poly in subtract_polys:
            diff = diff.difference(subtract_poly)

        output.append(diff)

# based on example in previous comment:
#   output looks like [1, 2-3, 3-2, 4-5-6, 5-4-6, 6-4-5]

Answer 3

To find all the intersections in a GeoDataFRame, you can use libpysal, see geopandas self intersection grouping (you can create fuzzy contiguity weights matrix based on intersection and get labels of connected components (ibpysal.weights.fuzzy_contiguity))

I use here a shapefile with some polygons that intersects and others not, and I want to preserve the original attributes.

import geopandas as gpd
gdf = gpd.read_file("poly_inters.shp")

import libpysal
W = libpysal.weights.fuzzy_contiguity(gdf,predicate='intersects')
print( W.component_labels)
[0 0 1 2 3 4 3]

There are 2 groups with polygon intersections, 0,1 and 4,6

gdf['inter']= W.component_labels #labels of connected components:
print(gdf)
    id                   geometry                          inter
0   0  POLYGON ((273.036 -149.055, 483.364 -149.317, ...      0
1   1  POLYGON ((274.142 -258.862, 275.036 -347.706, ...      0
2   2  POLYGON ((-51.161 16.699, 113.915 22.717, 82.9...      1
3   3  POLYGON ((409.675 -427.801, 579.049 -308.293, ...      2
4   4  POLYGON ((-147.455 -279.061, -110.316 -291.816...      3
5   5  POLYGON ((469.859 -21.991, 597.964 29.595, 616...      4
6   6  POLYGON ((-59.758 -235.213, 45.993 -221.457, -...      3


groups = gdf.groupby(['inter'])
def eliminter(g):
   gdf.loc[g.index[0],"geometry"] = g.geometry[g.index[0]]- g.geometry[g.index[1]]
   gdf.loc[g.index[1],"geometry"]=  g.geometry[g.index[1]]- g.geometry[g.index[0]]

for name, group in groups:
   if len(group) > 1: # intersection 
       eliminter(group)

You can also use the GeoPandas index (or the PyGEOS index, see GeoPandas issue: Add example or implement dissolve of contiguous or otherwise related geometries and performance example.ipynb)

heads, tails = gdf.sindex.query_bulk(gdf.geometry, predicate="intersects")
ix = heads != tails
heads = heads[ix]
tails = tails[ix]
print(list(zip(heads, tails))) 
[(0, 1), (1, 0), (4, 6), (6, 4)]

And you find the same groups.