如果点和多边形具有相同的最小边界框，则在多边形内查找点的空间索引
Spatial index to find points within polygon, if points and polygon have same minimum bounding box
我有一个形状匀称的 多边形 代表洛杉矶市的边界。我在 geopandas GeoDataFrame 中也有一组约 100 万经纬度 points，所有这些都落在该多边形的最小边界框内。其中一些点位于多边形本身内，但其他点不在多边形内。我只想保留洛杉矶边界内的那些点，并且由于洛杉矶的形状不规则，其最小边界框内的点中只有大约 1/3 在多边形本身内。
使用 Python，如果这些点和多边形具有相同的最小边界框，那么识别这些点中的哪些点位于多边形内的最快方法是什么？
我尝试使用 geopandas 及其 r-tree 空间索引：
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)]
这使用 GeoDataFrame 的 r-tree 空间索引快速找到 可能的 匹配，然后找到多边形与那些可能匹配的精确交集。但是，由于多边形的最小边界框与点集的最小边界框相同，r-tree 认为 每个点 都是可能的匹配项。因此，使用 r-tree 空间索引使得交集 运行 不会比没有空间索引时更快。此方法非常慢：大约需要 30 分钟才能完成。
我还尝试将我的多边形划分为小的子多边形，然后使用空间索引来查找哪些点可能与这些子多边形中的每一个相交。这种方法成功地找到了更少的可能匹配项，因为每个子多边形的最小边界框都比点集的最小边界框小得多。然而，将这组可能的匹配项与我的多边形相交仍然只节省了大约 25% 的计算时间，因此它仍然是一个极其缓慢的过程。
我应该使用更好的空间索引方法吗？如果点和多边形具有相同的最小边界框，找到哪些点在多边形内的最快方法是什么？
一个小例子来复制问题
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()
在一个简单的多边形上用一百万个点调用 intersects() 不会花费太多时间。使用 poly1，我得到以下图像。在多边形内找到 lat/lon 个点只需不到 10 秒。仅绘制边界多边形顶部的内部点如下所示：
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
Poly3 更大更有趣。新图像看起来像这样，它需要大约一分钟才能通过 bottle-neck intersects() 行。
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
所以罪犯不一定是lat/lon点数。同样糟糕的是边界多边形的复杂性。首先，我会推荐 poly.simplify()，或者您可以做的任何事情来减少边界多边形中的点数（显然不会大幅改变它）。
接下来，我建议考虑一些概率方法。如果点 p 周围的点都在边界多边形内，那么 p 很有可能也在边界多边形中。通常，速度和准确性之间存在一点 trade-off，但也许它可以减少您需要检查的点数。这是我使用 k-nearest neighbors classifier:
的尝试
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']])
给我这张图片。不完美，但是这个块的 %timeit 只需要 3.62 秒（n=50000 是 4.39 秒），相比之下检查每个点大约需要 50 秒。
如果相反，我只想删除有 30% 机会位于多边形中的点（只是扔掉明显的违规者并手动检查其余部分）。我可以使用 knn regression:
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>.30]
现在我只有大约 138000 个点要检查，如果我想使用 intersects() 检查每个点，这会很快。
当然，如果我增加邻居的数量或训练集的大小，我仍然可以获得更清晰的图像。这种概率方法的一些好处是 (1) 它是算法的，所以你可以把它放在任何时髦的边界多边形上，(2) 你可以轻松地调整它的准确性 up/down，(3) 它更快并且可以很好地缩放好吧（至少比蛮力更好）。
就像机器学习中的许多事情一样，有 100 种方法可以做到。希望这可以帮助您找到可行的方法。这是具有以下设置的另一张图片（使用分类器，而不是回归）。你可以看到它越来越好。
neigh = KNeighborsClassifier(n_neighbors=3, weights='distance')
df_short = df.sample(n=80000)
总结问题：当多边形的边界框与点集相同时，r-tree 将每个点都识别为可能的匹配项，因此没有提供加速。当结合大量点和具有大量顶点的多边形时，相交过程非常缓慢。
解决方案：从此geopandas r-tree spatial index tutorial，使用quadrat例程将多边形划分为子多边形。然后，对于每个子多边形，首先将其与点的 r-tree 索引相交以获得一小组可能的匹配，然后将这些可能的匹配与子多边形相交以获得精确匹配的集合。这提供了大约 100 倍的加速。
我有一个形状匀称的 多边形 代表洛杉矶市的边界。我在 geopandas GeoDataFrame 中也有一组约 100 万经纬度 points，所有这些都落在该多边形的最小边界框内。其中一些点位于多边形本身内，但其他点不在多边形内。我只想保留洛杉矶边界内的那些点，并且由于洛杉矶的形状不规则，其最小边界框内的点中只有大约 1/3 在多边形本身内。
使用 Python，如果这些点和多边形具有相同的最小边界框，那么识别这些点中的哪些点位于多边形内的最快方法是什么？
我尝试使用 geopandas 及其 r-tree 空间索引：
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)]
这使用 GeoDataFrame 的 r-tree 空间索引快速找到 可能的 匹配，然后找到多边形与那些可能匹配的精确交集。但是，由于多边形的最小边界框与点集的最小边界框相同，r-tree 认为 每个点 都是可能的匹配项。因此，使用 r-tree 空间索引使得交集 运行 不会比没有空间索引时更快。此方法非常慢：大约需要 30 分钟才能完成。
我还尝试将我的多边形划分为小的子多边形，然后使用空间索引来查找哪些点可能与这些子多边形中的每一个相交。这种方法成功地找到了更少的可能匹配项，因为每个子多边形的最小边界框都比点集的最小边界框小得多。然而，将这组可能的匹配项与我的多边形相交仍然只节省了大约 25% 的计算时间，因此它仍然是一个极其缓慢的过程。
我应该使用更好的空间索引方法吗？如果点和多边形具有相同的最小边界框，找到哪些点在多边形内的最快方法是什么？
一个小例子来复制问题
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()
在一个简单的多边形上用一百万个点调用 intersects() 不会花费太多时间。使用 poly1，我得到以下图像。在多边形内找到 lat/lon 个点只需不到 10 秒。仅绘制边界多边形顶部的内部点如下所示：
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
Poly3 更大更有趣。新图像看起来像这样，它需要大约一分钟才能通过 bottle-neck intersects() 行。
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
所以罪犯不一定是lat/lon点数。同样糟糕的是边界多边形的复杂性。首先，我会推荐 poly.simplify()，或者您可以做的任何事情来减少边界多边形中的点数（显然不会大幅改变它）。
接下来，我建议考虑一些概率方法。如果点 p 周围的点都在边界多边形内，那么 p 很有可能也在边界多边形中。通常，速度和准确性之间存在一点 trade-off，但也许它可以减少您需要检查的点数。这是我使用 k-nearest neighbors classifier:
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']])
给我这张图片。不完美，但是这个块的 %timeit 只需要 3.62 秒（n=50000 是 4.39 秒），相比之下检查每个点大约需要 50 秒。
如果相反，我只想删除有 30% 机会位于多边形中的点（只是扔掉明显的违规者并手动检查其余部分）。我可以使用 knn regression:
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>.30]
现在我只有大约 138000 个点要检查，如果我想使用 intersects() 检查每个点，这会很快。
当然，如果我增加邻居的数量或训练集的大小，我仍然可以获得更清晰的图像。这种概率方法的一些好处是 (1) 它是算法的，所以你可以把它放在任何时髦的边界多边形上，(2) 你可以轻松地调整它的准确性 up/down，(3) 它更快并且可以很好地缩放好吧（至少比蛮力更好）。
就像机器学习中的许多事情一样，有 100 种方法可以做到。希望这可以帮助您找到可行的方法。这是具有以下设置的另一张图片（使用分类器，而不是回归）。你可以看到它越来越好。
neigh = KNeighborsClassifier(n_neighbors=3, weights='distance')
df_short = df.sample(n=80000)
总结问题：当多边形的边界框与点集相同时，r-tree 将每个点都识别为可能的匹配项，因此没有提供加速。当结合大量点和具有大量顶点的多边形时，相交过程非常缓慢。
解决方案：从此geopandas r-tree spatial index tutorial，使用quadrat例程将多边形划分为子多边形。然后，对于每个子多边形，首先将其与点的 r-tree 索引相交以获得一小组可能的匹配，然后将这些可能的匹配与子多边形相交以获得精确匹配的集合。这提供了大约 100 倍的加速。