2

I have multiple files that essentially have anywhere from 2000-5000 points each.

I have a polygon file that has about 10000 polygon shapes.

For each file containing the points, I want to get all the polygons that the file intersected.

Checking each polygon for each point and using Contains() takes too long. 2000 x 10000 iterations for each file essentially.

My current approach is to create a polygon around the points then get the intersections of that one bounding polygon to the 10000. This is fast but leads to me picking up intersecting polygons that don't actually contain the points but intersect the bounding polygon I had made around the points. So I have to go back then and do another check for intersections that don't actually contain points.

#go through files containing points and create a bounding polygon around each
for file in point_files:
        pts_ds = ogr.Open(file)
        pts_lyr = pts_ds.GetLayer()
        geom_collection = ogr.Geometry(ogr.wkbGeometryCollection)
        for feat in pts_lyr:
            geom_collection.AddGeometry(feature.GetGeometryRef())
        bounding_polygon = geom_collection.ConvexHull()
        #then save bounding polygon to file

#open up each of the bounding polygons I had made and intersect it against the comparison polygons
comparison_polygons_ds = ogr.Open(comparison_polygons.shp)
comparison_polygons_lyr = comparison_polygons.GetLayer()
matches = {}
for file in bounding_polygon_files:
    b_polygon_ds = ogr.Open(file)
    b_polygon_lyr = b_polygon_ds.GetLayer()
    for b_polygon_feat in b_polygon_lyr:
        b_polygon_geom = b_polygon_feat.GetGeometryRef()
        for comparison_polygons_feat in comparison_polygons_lyr:
            comparison_polygon_geom = comparison_polygons_feat.GetGeometryRef() 
            if b_polygon_geom.Intersects(comparison_polygon_geom):
                matches.setdefault(comparison_polygon_geom.GetField(“ID”), []).append(file)

The end result I want is essentially: file1 intersects x,y,z polygons. file2 intersects a,b,c polygons.

How do I achieve that?

I am unfamiliar with GIS and using Python with the OGR package.

Improve this question
7
  • 1
    Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works
    – Ian Turton ♦
    Apr 24, 2021 at 15:01
  • @IanTurton Hi Ian, I edited the post with code I'm currently using. Looking for a better approach since it leads to an unnecessary step of having to do a second cleanup. Is this good?
    – BlueCoop
    Apr 24, 2021 at 15:18
  • It's possible that contains() does a bounding box search with each polygon as well, for the same reasons that you have. Have you looked at shapely, the topological library? It may allow you to customise your methodology a little better than OGR. I think OGR uses shapely internally. Also, if the points in each file are geographically sorted, maybe you can subset them before determining the MBR. Or compute an MBR for each polygon and test that against the MBR of the points. How long is it taking to run a file against the polygons?
    – wingnut
    Apr 24, 2021 at 15:57
  • Comparing each point against each polygon with Contains(), it ended up taking hours before I manually stopped it. Making a bounding box around each file's points, then getting the intersecting polygons of the bounding polygon doesn't take long but the intersections obviously aren't accurate to the actual points. But it filters it down at least then I do another second round based on the actual points, at which point the whole thing takes about 1 minute for all files. Not bad but it feels like I'm doing an unnecessary intermediary step. I'll check out Shapely!
    – BlueCoop
    Apr 24, 2021 at 16:15
  • I guess I'm assuming there's a more efficient way because I can imagine that rather than checking if each polygon contains the point, it might be possible somehow to just get the coordinate of the point I'm checking, place it in the polygon file, then just see what polygon that coordinate is in more geographically/mathematically than iteratively checking every polygon.
    – BlueCoop
    Apr 24, 2021 at 16:18

1 Answer 1

Reset to default
1

I aggregated each files points into a MultiPoint and ran Intersects() on that for every comparison polygon and it was a lot faster (immediate). I'm not sure why since I imagine it still has to test every coordinate within the MultiPoint and see which polygons intersect with it. Maybe a MultiPoint object allows OGR to incorporate better spatial indexing algorithms as suggested by @bugmenot123.

Not completely sure how its working under the hood but it gives me the same results as I got with the bounding polygon method after second pass through cleaning, so results appear accurate.

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.