Spatial index to find points within polygon, if points and polygon have same minimum bounding box
<p>I have a shapely <em>polygon</em> representing the boundaries of the city of Los Angeles. I also have a set of ~1 million lat-long <em>points</em> in a geopandas GeoDataFrame, all of which fall within that polygon's minimum bounding box. Some of these points lie within the polygon itself, but others do not. I want to retain only those points within Los Angeles's boundaries, and due to Los Angeles's irregular shape, only approximately 1/3 of the points within its minimum bounding box are within the polygon itself.</p>

<p><strong>Using Python, what is the fastest way to identify which of these points lie within the polygon, given that the points and the polygon have the same minimum bounding box?</strong></p>

<p>I tried using geopandas and its r-tree spatial index:</p>

<pre><code>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)]
</code></pre>

<p>This uses the GeoDataFrame's r-tree spatial index to quickly find the <em>possible</em> matches, and then finds the exact intersection of the polygon and those possible matches. However, because the polygon's minimum bounding box is the same as the set of points' minimum bounding box, r-tree considers <em>every point</em> to be a possible match. Thus, using an r-tree spatial index makes the intersection run no faster than it would without the spatial index. This method is very slow: it takes ~30 minutes to complete.</p>

<p>I also tried dividing my polygon into small sub-polygons, then using the spatial index to find which points possibly intersect with each of these sub-polygons. This method successfully finds fewer possible matches because each of the sub-polygons' minimum bounding boxes is much smaller than the set of points minimum bounding box. However, intersecting this set of possible matches with my polygon still only shaves off ~25% of my computation time, so it's still a brutally slow process.</p>

<p>Is there a better spatial index method I should use? And what is the fastest way to find which points are within the polygon, if the points and polygon have the same minimum bounding box?</p>

<p>A little example to duplicate the problem a bit</p>

<pre><code>import pandas as pd
import shapely
import matplotlib.pyplot as plt

from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon
from shapely.geometry import Point
import seaborn as sns
import numpy as np

# some lon/lat points in a DataFrame
n = 1000000
data = {'lat':np.random.uniform(low=0.0, high=3.0, size=(n,)), 'lon':np.random.uniform(low=0.0, high=3.0, size=(n,))}
df = pd.DataFrame(data)

# the 'bounding' polygon
poly1 = shapely.geometry.Polygon([(1,1), (1.5,1.2), (2,.7), (2.1,1.2), (1.8,2.3), (1.6,1.8), (1.2,3)])
# poly2 = shapely.geometry.Polygon([(1,1), (1.3,1.6), (1.4,1.55), (1.5,1.2), (2,.7), (2.1,1.2), (1.8,2.3), (1.6,1.8), (1.2,3), (.8,1.5),(.91,1.3)])
# poly3 = shapely.geometry.Polygon([(1,1), (1.3,1.6), (1.4,1.55), (1.5,1.2), (2,.7), (2.1,1.2), (1.8,2.3), (1.6,1.8), (1.5,2), (1.4,2.5),(1.3,2.4), (1.2,3), (.8,2.8),(1,2.8),(1.3,2.2),(.7,1.5),(.66,1.4)])

# limit DataFrame to interior points
mask = [poly1.intersects(shapely.geometry.Point(lat,lon)) for lat,lon in zip(df.lat,df.lon)]
df = df[mask]

# plot bounding polygon
fig1, ax1 = sns.plt.subplots(1, figsize=(4,4))
patches  = PatchCollection([Polygon(poly1.exterior)], facecolor='red', linewidth=.5, alpha=.5)
ax1.add_collection(patches, autolim=True)

# plot the lat/lon points
df.plot(x='lat',y='lon', kind='scatter',ax=ax1)
plt.show()
</code></pre>

<p>Calling intersects() with a million points on a simple polygon doesn't take very much time. Using poly1, I get the following image. Finding the lat/lon points inside the polygon takes less than 10 seconds. Plotting only the interior points on top of the bounding polygon looks like this:</p>

<p><a href="http://i.stack.imgur.com/8us3h.png" rel="nofollow"><img src="http://i.stack.imgur.com/8us3h.png" alt="enter image description here"></a></p>

<pre><code>In [45]: %timeit mask = [Point(lat,lon).intersects(poly1) for lat,lon in zip(df.lat,df.lon)]
1 loops, best of 3: 9.23 s per loop
</code></pre>

<p>Poly3 is a bit bigger and funnier. The new image looks like this and it takes about a minute to get through the bottle-neck intersects() line. </p>

<p><a href="http://i.stack.imgur.com/HEunh.png" rel="nofollow"><img src="http://i.stack.imgur.com/HEunh.png" alt="enter image description here"></a></p>

<pre><code>In [2]: %timeit mask = [poly3.intersects(shapely.geometry.Point(lat,lon)) for lat,lon in zip(df.lat,df.lon)]
1 loops, best of 3: 51.4 s per loop
</code></pre>

<p>So the criminal isn't necessarily the number of lat/lon points. Just as bad is the complexity of the bounding polygon. First, I would recommend <code>poly.simplify()</code>, or anything you can do to reduce the number of points in your bounding polygon (without changing it drastically, obviously).</p>

<p>Next, I would recommend thinking of some probabilistic approach. If a point <code>p</code> is surrounded by points which are all inside the bounding polygon, there's a pretty good chance <code>p</code> is also in the bounding polygon. Generally, there's a bit of a trade-off between speed and accuracy, but maybe it can be a way of reducing the number of points you need to check. Here's my attempt using <a href="http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsClassifier.html#sklearn.neighbors.KNeighborsClassifier" rel="nofollow">k-nearest neighbors classifier</a>:</p>

<pre><code>from sklearn.neighbors import KNeighborsClassifier

# make a knn object, feed it some training data
neigh = KNeighborsClassifier(n_neighbors=4)
df_short = df.sample(n=40000)
df_short['labels'] = np.array([poly3.intersects(shapely.geometry.Point(lat,lon)) for lat,lon in zip(df_short.lat,df_short.lon)])*1
neigh.fit(df_short[['lat','lon']], df_short['labels'])

# now use the training data to guess whether a point is in polygon or not
df['predict'] = neigh.predict(df[['lat','lon']])
</code></pre>

<p>Gives me this image. Not perfect, but %timeit for this block takes only 3.62 seconds (4.39 seconds for n=50000), compared to about 50 seconds to check every single point.</p>

<p><a href="http://i.stack.imgur.com/xeXnX.png" rel="nofollow"><img src="http://i.stack.imgur.com/xeXnX.png" alt="enter image description here"></a></p>

<p>If instead, I just want to drop points that have, say, a 30% chance of being in the polygon (just throwing out the obvious offenders and checking the rest by hand). I can use a <a href="http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsRegressor.html#sklearn.neighbors.KNeighborsRegressor" rel="nofollow">knn regression</a>:</p>

<pre><code>from sklearn.neighbors import KNeighborsRegressor
neigh = KNeighborsRegressor(n_neighbors=3, weights='distance')
#everything else using 'neigh' is the same as before

# only keep points with more than 30\% chance of being inside
df = df[df.predict&gt;.30]
</code></pre>

<p>Now I only have about 138000 points to check, which would go pretty quickly if I wanted to check each one using <code>intersects()</code>. </p>

<p>Of course if I increase the number of neighbors, or the size of the training set, I can get sharper images still. Some nice things about this probabilistic approach are (1) it's algorithmic, so you can throw it at any funky bounding polygon, (2) you can easily tune its accuracy up/down, (3) it's a lot faster and scales pretty well (better to brute force at least).</p>

<p>Like many things in machine learning, there can be 100 ways to do it. Hopefully this helps you figure something out that works. Here's one more picture with the following settings (using a classifier, not a regression). You can see it's getting better.</p>

<pre><code>neigh = KNeighborsClassifier(n_neighbors=3, weights='distance')
df_short = df.sample(n=80000)
</code></pre>

<p><a href="http://i.stack.imgur.com/GhHbt.png" rel="nofollow"><img src="http://i.stack.imgur.com/GhHbt.png" alt="enter image description here"></a></p>

python gis geospatial shapely geopandas