I need to do some processing with groups of polygons. If any of the polygons in the group overlaps other polygon the result is unpredictable and don't throw any error. So, I need to detect if there is any overlapping polygon inside the group and which polygons are overlapping.

The groups are made of about 15 polygons, so checking each polygon against the others would be costly. I wonder if there is an efficient way to do that.

The code is in Python and can use either the OGR bindings or Shapely.

A theoretical answer is also appreciated.

  • 2
    In principle, any GIS platform that forms the (set-theoretic) union of each group of polygons will do so efficiently (using, say, a line-sweep algorithm); that is to say, it should require O(nlog(n)) time for a total number of *n vertices. Computing the sum of the original areas and computing the area of the union are both O(n) operations. The latter area is less than the former area (up to floating point error) if and only if the interiors of two (or more) of the polygons overlap. Therefore, this test should asymptotically be O(n*log(n)), which for general polygons is optimal. – whuber Jun 8 '13 at 19:08
  • 1
    @whuber, nice, that give me an idea of a good start: First get all individual areas (size) and store them. Second check if there is any intersection in the whole group by checking the sum of individual areas against the union area. Only if true, then check each one, since the individual area size is already stored it's 2 area calculations less for each comparison. I'll develop the script a little more, as soon as I get something consistent I'll post the results. – Pablo Jun 12 '13 at 13:40

There are 105 combinations for testing the intersections of pairs of fifteen polygons, or mathematically (using combinations) expressed as 15 choose 2, or symbolically (15 over 2).

To borrow some similar logic with Shapely, you can do this:

from itertools import combinations

# gather the 15 Shapely polygons in something iterable, like a list
shapes = [poly1, poly2, ..., poly15]

# test the intersection on the combinations of pairs
intersections = [pair[0].intersects(pair[1]) for pair in combinations(shapes, 2)]
# intersections is a list with 105 elements of True or False

if any(intersections):
    print("yes, %d of %d combinations intersect" % (sum(intersections), len(intersections)))
    print("no intersections")
| improve this answer | |

I recommend using OGR because it will leverage a Shapefile Spatial Index, either .sbn/.sbx or .qix. Iterating thru each group, you can first perform a faster Overlap test between BoundingBox geometries (lyr.SetSpatialFilterRect). When true, test for overlap on the real geometries. Another option, Spatialite is easy to setup and offers PostGIS like functions and performance. It may be worth investigating. There is a Python module pyspatialite.

| improve this answer | |

You can take advantage of matrix operations for this problem and achieve an O(n) time complexity.

The idea is to have an ones-and-zeros matrix e.g. mask per polygon. Then all you need to do is to add these matrices together(linear time and very efficiently if you use numpy sparse matrix). If there is a number greater than one it means there is overlap.

If you have two decimal points you can multiply by 100 to still be able to convert to a matrix.

You can extend the same solution to find out which ones are overlapping if you choose exponentially increasing numbers. For instance, the first polygon is a matrix of ones-and-zeros, the second is tens-and-zeros, third is hundreds and zeros, etc. this way if you have a 1000101 in a cell at the end it means polygons number 1, 3, and 7 overlap.

| improve this answer | |

It's not a scripting answer but you could put your data into a feature data set and create some topology rules. This all requires a geodatabase.


That will give you information on the basics, creating a topology, validating it and correcting errors. The result is basically something that you can add to your ArcMap view and it will highlight areas that are violating the topology rules you have set.

Here are topology rules: http://resources.arcgis.com/en/help/main/10.1/index.html#//01mm0000000m000000 (note that other links in this sub menu on the left side are pertinent to this subject and will be helpful to you). You'll probably want the polygon rule "Must Not Overlap"

| improve this answer | |

You could do a spatial join on the features that intersect, then select the input polygons that were joined and collect the ID numbers. In your code add a condition to check your lest of IDs and skip (or do a different operation on) the polygons that are in the list of overlapping IDs. Might be a little slow upfront, but once you had the list it shouldn't slow your script much.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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