4

I am trying to find the largest intersecting area for each Polygon in a layer1 compared to the geometries in a layer2 in python. I use geopandas.overlay for it.

However, the geometries in layer2 don't cover the entire area, and for my needs the area NOT covered is also a candidate for being the largest intersecting area. Hence, I added a single polygon with the difference of the total bbox of layer1 with all geometries in layer2.

Before adding the "difference polygon" to layer2 it took about 10 minutes to calculate the intersection for an area of about 100+ municipalities. After adding it, it took the same time for about five municipalities: so about 2 orders of magnitude slower.

I applied multiprocessing to speed up the calculation, but the loss in speed is still about one order of magnitude additional time consumed.

The question is why did it become so slow? I only added a single polygon to one of the GeoDataFrames. I did assume, the new polygon has a lot of holes, so I at maximum doubled the number of points in all the polygons. Why is it so much slower than half-speed? Does it lead to problems with spatial indexing?

2
  • 1
    Intersection is one of those O(NxM) algorithms, that's a function of segments in both the input features. Reducing complexity of one helps, but of both is better.
    – Vince
    Commented May 18 at 13:50
  • If you use Union instead of Intersection you dont need to create the large polygon
    – Bera
    Commented May 20 at 15:39

1 Answer 1

5

1. Why is it so much slower

  • The spatial index works as follows: it ~stores the bboxes of the geometries of a layer in an efficient tree structure so it can be very quickly determined if a given bbox intersects with one or more of the bboxes of the geometries in the layer. This way the "intersection candidates" between two layers can be determined very quickly, and the proper intersection calculation only needs to be done between geometries in layer1 and layer2 that most likely have an intersection instead of having to calculate it between all combinations. However, as the huge geometry added to layer2 has a huge bbox (the total bbox of layer1), the spatial index won't do anything for this particular geometry as each geometry in layer1 will intersect with that big bbox. Hence, none of the expensive intersection calculations can be avoided by the spatial index.
  • Based on my experiences, calculating an intersection with a geometry of up to 10.000 points is quite fast, but higher than that becomes quickly slower and slower. So many very slow calculations = very very slow.

2. Some possible solutions

The typical thing to do for your need is to use geopandas.overlay(..., how="identity"). The "identity" overlay type performs an "intersection", but also adds all the parts of the geometries in layer1 that don't intersect layer2 to the output. All columns in the output that originate from layer2 will be None for those features. This way you can avoid adding the huge "inversed" polygon and avoid slowing things down. The time taken will be a bit longer than double the time of just the basic "intersection", but not orders of magnitude longer.

If you need/want a faster solution, you can try geofileops.identity, it uses multiprocessing and some other tricks to speed up the processing. Disclaimer: I'm the developer of geofileops.

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