Checking if a point is contained in a polygon/multipolygon - for many points

user_1016428 0
Viewed 1745

I have a 2D map divided in a grid of rectangles - about 45,000 of them - and a number of polygons/multipolygons derived from shapefiles (I currently read them with the shapefile library pyshp). Unfortunately, a few of these polygons are rather complex and made up by a large number of points (one of them has 640,000 points) and can have holes in them.

What I am trying to do is to check, for each of these polygons, which cell centers (the cells of my grid) fall inside that specific polygon. However, having about 45,000 cell centers and 150 polygons is taking quite a while to check everything using shapely. This is what I am doing, more or less:

# nx and ny are the number of cells in the x and y direction respectively
# x, y are the coordinates of the cell centers in the grid - numpy arrays
# shapes is a list of shapely shapes - either Polygon or MultiPolygon
# Point is a shapely.geometry.Point class
for j in range(ny):
    for i in range(nx):
        p = Point([x[i], y[j]])
        for s in shapes:
            if s.contains(p):
                # The shape s contains the point p - remove the cell

Depending on the polygon in question, it takes between 200 microseconds to 80 milliseconds per check, and with more than 6 million checks to do the runtime goes easily into the hours.

I guess there must be a more intelligent way to handle this problem - I'd rather stay with shapely if possible, but any suggestion is most appreciated.

Thank you in advance.

shapely gis python
2 Answers

In the accepted answer, the use of rtrees in shapely only supports the extend - or bounding box - of a polygon. As stated in the shapely documentation:

https://shapely.readthedocs.io/en/stable/manual.html#str-packed-r-tree

as of 1.7.1 the rtree.query needs to be followed up by a check whether the polygon actually contains the point. This means additional checks but could still speed the algorithm up if the polygons are not in a very jumbled configuration

That is, the code should read something like

from shapely.strtree import STRtree
from shapely.geometry import Polygon, Point    

# this is the grid. It includes a point for each rectangle center. 
points = [Point(i, j) for i in range(10) for j in range(10)]
tree = STRtree(points) # create a 'database' of points

# this is one of the polygons. 
p = Polygon([[0.5, 0], [2, 0.2], [3, 4], [0.8, 0.5], [0.5, 0]])

res = [o for o in tree.query(p) if p.contains(o)] # run the query (a single shot) - and test if the found points are actually inside the polygon. 

# do something with the results
print([o.wkt for o in res])

Indeed, the result is now

['POINT (2 1)', 'POINT (2 2)']

Which can be seen to be the points in the interior of the polygon here: enter image description here

-- bonus edit ---

I have observed better speed improvements by localizing the problem in the following sense:

A square grid is made over the area, and STRtree structures are made of the points and the polygons. The points and polygons are queried on each individual grid-polygon - within each grid-polygon the queried points are checked for containement inside the queried polygons. Thereby reducing the amount of checks to individual polygons.

from shapely.strtree import STRtree
from shapely.geometry import Polygon, Point
import random    

# build a grid spread over the area. 
grid = [Polygon([[i, j],[i+1,j],[i,j+1],[i+1,j+1]]) for i in range(10) for j in range(10)]

pointtree = STRtree([Point(random.uniform(1,10),random.uniform(1,10)) for q in range(50)]) # create a 'database' of 50 random points - and build a STRtree structure over it

# take the same polygon as above, but put in an STRtree 
polytree = STRtree([Polygon([[0.5, 0], [2, 0.2], [3, 4], [0.8, 0.5], [0.5, 0]])])
res = []
#do nested queries inside the grid
for g in grid:
    polyquery = polytree.query(g)
    pointquery = pointtree.query(g)
    for point in pointquery:
        for poly in polyquery:
            if poly.contains(point):
                 res.append(point)

# do something with the results
print([o.wkt for o in res])
user_58431790

If you'd like to perform less comparison operations, you can try to use the shapely str-tree feature. Consider the following code: 

from shapely.strtree import STRtree
from shapely.geometry import Polygon, Point

# this is the grid. It includes a point for each rectangle center. 
points = [Point(i, j) for i in range(10) for j in range(10)]
tree = STRtree(points). # create a 'database' of points

# this is one of the polygons. 
p = Polygon([[0.5, 0], [2, 0.2], [3, 4], [0.8, 0.5], [0.5, 0]])

res = tree.query(p). # run the query (a single shot). 

# do something with the results
print([o.wkt for o in res])

The result of this code is: 

['POINT (3 0)', 'POINT (2 0)', 'POINT (1 0)', 'POINT (1 1)', 'POINT (3 1)',
 'POINT (2 1)', 'POINT (2 2)', 'POINT (3 2)', 'POINT (1 2)', 'POINT (2 3)',
 'POINT (3 3)', 'POINT (1 3)', 'POINT (2 4)', 'POINT (1 4)', 'POINT (3 4)']
user_11055600
Related