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)]
  • I do not quite understand. Bounding box of a point is a point, isn't it? – user30184 Sep 26 '16 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 '16 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 '16 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 '16 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 '16 at 17:12
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

  • 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 '16 at 22:34

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.