Python Shapely – How to Remove Overlapping Polygons in One Dataset

iterationoverlapping-featurespolygonpythonshapely

I am looking for a method using python to remove a large amount of overlapping polygons from one collection. These polygons are grouped together based on imagery collection and have been stored as WKB in a database. I am able to pull the WKBs out per collection. My current approach is to create a list of geometries and iterate using shapely.intersects per polygon to an empty list. If polygon in listA does not intersect with polygon in listB the polygon in listA is appended to listB and removed from listA. My thinking is as this iterates I will be able to move polygons with no overlaps to one list and keep the overlapping polygons in another. When I run this code I get the error: "AttributeError 'list' object has no attribute '_geom'".

polygon_list_A = [<shapely.geometry.polygon.Polygon object at 0x000005B3A>, <shapely.geometry.polygon.Polygon object at 0x000006E7A>, ...]

polygon_list_B = []

for polygon_A in polygon_list_A:
    if polygon_A.intersects(polygon_list_B) == True:
        polygon_list_B.append(polygon_A)
        polygon_list_A.remove(polygon_A)
    else:
        pass

Best Answer

AttributeError 'list' object has no attribute '_geom'

Why ? Because polygon_list_B is an empty list and not a Shapely geometry in if polygon_A.intersects(polygon_list_B): (you don't need == True, it is implicit)

A possible solution is to iterate through all pairs of geometries in polygon_list_A to find the intersections (as in Finding the common borders between polygons in the same shapefile)

import itertools
polygon_list_B = []
for geom1,geom2 in  itertools.combinations(polygon_list_A, 2):
     if geom1.intersects(geom2):
        polygon_list_B.append(geom1)
        polygon_list_A.remove(geom1)

polygon_list_A -> dark green colour, polygon_list_B -> light green colour

If you only want the polygons without intersection

polygon_list_B = []          
for geom1,geom2 in  itertools.combinations(polygon_list_A, 2):
    if geom1.intersects(geom2):
         polygon_list_B.append(geom1)
         polygon_list_B.append(geom2)
 result = [geom for geom in polygon_list_A if geom not in polygon_list_B]

result -> dark green colour, polygon_list_B -> light green colour

Related Solutions

[GIS] Counting overlapping polygons from one Shapefile using ArcGIS Desktop

I would recommend using the Count Overlapping Features (Analysis) tool.

Generates planarized overlapping features from the input features. The count of overlapping features is written to the output features.

[GIS] Loop to check multiple polygons overlap in r

I played with assign(), ls(), and mget() to accomplish something that I believe will improve your workflow. First I use ls() to get a list of all environment variables that start with "p":

names_poly <- ls(pattern='^p.')

I used combn() to find all the unique combinations of polygons

combos <- combn(names_poly,2)

I looped over combos using get() to fetch the polygons i and j named in combos. I also used assign() to create a new object whose name is pasted together (e.g. int_p1_p2)

for(k in seq_along(combos[1,])){
    i <- combos[1,k]
    j <- combos[2,k]
    print(paste("intersecting",i,j))
    assign(paste("int",i,j,sep="_"),gIntersection(get(i), get(j), byid = FALSE))
}

You may not need this step below since you have all the new int_pi_pj SpatialPolygons as separate objects, but if you want them in a list:

int <- mget(ls(pattern='^int_.'))
rm(list = ls(pattern='^int_.')   #cleanup

Then you can access the items in your list with $:

plot(int$int_p1_p2)

EDIT: In the above, when there is no overlap, the variable will be NULL.If you don't want to even see the combinations where there is no overlap, you could add this if statement to the same for loop:

for(k in seq_along(combos[1,])){
    i <- combos[1,k]
    j <- combos[2,k]
    print(paste("intersecting",i,j))
    if(gIntersects(get(i), get(j), byid = FALSE)){
        assign(paste("int",i,j,sep="_"),gIntersection(get(i), get(j), byid = FALSE))
    }
}

Best Answer

Related Solutions

[GIS] Counting overlapping polygons from one Shapefile using ArcGIS Desktop

[GIS] Loop to check multiple polygons overlap in r

Related Question