1

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:

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.

10
  • This is a typical overlay problem, and if you have access to geopandas, you might find a solution there. geopandas.org/en/stable/docs/user_guide/set_operations.html
    – Mintx
    Commented Apr 20, 2022 at 17:15
  • 3
    What have you tried so far? Have you seen the .symmetric_difference() in the shapely package? Please, be aware that your question can be closed, because it does not contain any code and you used the python tag.
    – Taras
    Commented Apr 20, 2022 at 17:15
  • Thanks for the tips. I am aware of overlay but I am not sure how to apply it in this case. I do use geopandas
    – Jan Pisl
    Commented Apr 20, 2022 at 17:29
  • All the examples I've seen, incl. geopandas, work with two dataframes. I only have one and if I use it twice then every polygon intersects with itself for example so it's not quite the same
    – Jan Pisl
    Commented Apr 20, 2022 at 17:32
  • 1
    I agree with you that a simple, well written question is beneficial. However, such an open question has little value. By not mentioning what you have already tried, or what packages etc you are using, leads to people spending time making suggestions that are unsuitable. I apologise for coming across unfriendly, I was frankly baffled by your aversion for including a code sample that would allow readers to see where you get stuck and serve as a starting point for offering help.
    – Matt
    Commented Apr 22, 2022 at 23:10

3 Answers 3

3

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

enter image description here

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]
2

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

enter image description here

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)

enter image description here

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.

2
  • how and where do you visualize those polygons? matplotlib?
    – Taras
    Commented Jun 13, 2022 at 5:22
  • Yes, Matplotlib only
    – gene
    Commented Jun 13, 2022 at 6:48

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