1

I have a MultiPolygon representing the city boundary of Houston and it has an extremely complicated boundary. I also have a set of ~900,000 points (that has the same minimum bounding box as that Houston polygon). About ~400,000 of these points are within the polygon but the others lie outside it. Using python, geopandas, and shapely I tried intersecting this polygon with my points using r-tree. But because the points and polygon have the same minimum bounding box, r-tree offers no speed-up. The process currently takes 30+ minutes.

Which (if any) type of spatial index can I use to accelerate my intersection query when the polygon and points have the same minimum bounding box?

Edit to add code snippet here:

sindex = gdf['geometry'].sindex
possible_matches_index = list(sindex.intersection(polygon.bounds))
possible_matches = gdf.iloc[possible_matches_index]
points_in_polygon = possible_matches[possible_matches.intersects(polygon)]
Improve this question
12
  • I do not quite understand. Bounding box of a point is a point, isn't it?
    – user30184
    Sep 26, 2016 at 15:59
  • The two bounding boxes in question are 1) the minimum bounding box of the polygon, and 2) the minimum bounding box of a set of 900,000 points.
    – eos
    Sep 26, 2016 at 16:02
  • But your query compares bbox of each point with bbox of the polygon one by one. At least I hope so. It may still be not selective if the points outside the polygon are not outside the bbox of the polygon. Is that your case?
    – user30184
    Sep 26, 2016 at 16:03
  • I believe you want to do this gaia-gis.it/spatialite-3.0.0-BETA1/WorldBorders.pdf. Split your massive multipolygon to small polygons with few vertices and your spatial index will be selective and rock.
    – user30184
    Sep 26, 2016 at 16:15
  • Yes, basically the spatial index doesn't help because none of the points are outside the bounding box of the polygon. I previously tried what you suggested and divided my polygon into 1,000 sub-polygons then did a fast r-tree intersect to get possible matches. Then I intersected the possible matches with the full polygon geometry to get the actual precise matches. This was the fastest method I could come up with, but it really only yields speed-ups of 2x-3x (almost entirely because that intersection of possible matches with with full geom is still slow). I was hoping for like 100x speed ups.
    – eos
    Sep 26, 2016 at 17:12

1 Answer 1

Reset to default
1

I try to explain with an image. If you do not agree that this strategy works I will delete my answer. However, I will continue to use it in my own work.

I took Texas as an example and I split it into 1 by 1 degree pieces. The grid that shows behind corresponds with the R-Tree index despite at the boundaries where the real R-Tree BBOXes are smaller.

In the Texas case a point never hits more than one R-Tree box. In case of adjacent polygons there may be more hits, but never more than a handful. Next thing to do is to make an Intersects test with a true geometry of the corresponding piece. If point is inside that piece it is inside Texas as well, and in the USA, North America etc.

enter image description here

Improve this answer
1
  • Yes that image helps explain the logic. I will try to implement something like this in my python project, seems like it should offer considerable speed-up.
    – eos
    Sep 27, 2016 at 22:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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