这里有很多关于如何高效地在多边形中匹配点的问题（例如：这里和这里）。这些问题涉及到的主要变量是大量点N和多边形顶点数V。 这些都是很好很有用的，但我想处理点数N很多且多边形G很多的情况。这也意味着我的输出将不同（我主要看到的输出由落入多边形内的点组成，但我想知道与一个点相关联的多边形）。
我有一个包含大量多边形（数十万个）的shapefile文件。多边形可以相接，但它们之间几乎没有重叠（任何内部重叠都是错误的结果 - 比如人口普查区块组）。 我还有一个包含数百万个点的csv文件，并且如果可能的话，我想根据点位于哪个多边形中对这些点进行分类。有些点可能不属于多边形（继续上面的例子，考虑在海洋上方的点）。下面我设置了一个玩具示例来研究这个问题。
设置：
创建多边形的地理数据框（模拟使用geopandas导入的shapefile）：
返回：
我有一个包含大量多边形（数十万个）的shapefile文件。多边形可以相接，但它们之间几乎没有重叠（任何内部重叠都是错误的结果 - 比如人口普查区块组）。 我还有一个包含数百万个点的csv文件，并且如果可能的话，我想根据点位于哪个多边形中对这些点进行分类。有些点可能不属于多边形（继续上面的例子，考虑在海洋上方的点）。下面我设置了一个玩具示例来研究这个问题。
设置：
import numpy as np
from shapely.geometry import shape, MultiPolygon, Point, Polygon
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from shapely.strtree import STRtree
#setup:
np.random.seed(12345)
# shape gridsize:
gridsize=10
avgpointspergridspace=10 #point density
创建多边形的地理数据框（模拟使用geopandas导入的shapefile）：
# creating a geodataframe (shapefile imported via geopandas):
garr=np.empty((gridsize,gridsize),dtype=object)
for i in range(gridsize):
for j in range(gridsize):
garr[i,j]=Point(i,j)
# polygons:
poly_list=[]
for i in range(gridsize-1):
for j in range(gridsize-1):
temp_points=[garr[i,j],garr[i,j+1],garr[i+1,j+1],garr[i+1,j],garr[i,j]]
poly=Polygon([[p.x,p.y] for p in temp_points])
poly_list+=[poly]
# creating a geodataframe, including some additional numeric and string variables:
gdf=gpd.GeoDataFrame()
gdf['geometry']= poly_list
gdf['id']=list(range(len(gdf['geometry'])))
gdf['numeric']=0
gdf['string']='foo'
# creating some holes in the grid:
gdf['included']=[np.random.choice([True,False],p=[.95,.05]) for x in range(len(gdf))]
gdf_polys=gdf[gdf['included']]
生成点（模拟通过pandas导入的csv）：
# creating a pandas dataframe with points (csv of coordinates imported to pandas):
npoints=(gridsize+2)**2*10
fgridend=gridsize+1
fgridstart=-1
xlist=[]
ylist=[]
points=[]
for i in range(npoints):
x=fgridstart+np.random.random()*fgridend
y=fgridstart+np.random.random()*fgridend
xlist+=[x]
ylist+=[y]
df=pd.DataFrame(list(zip(xlist,ylist)),columns=['x','y'])
coords=[Point(xy) for xy in zip(df['x'],df['y'])]
gdf_points=gpd.GeoDataFrame(df,geometry=coords)
绘制结果：
fig, ax = plt.subplots(figsize=[10,10])
gdf_polys.plot(ax=ax,facecolor='gray',alpha=.2,edgecolor='black',lw=2)
gdf_points.plot(ax=ax)
返回：
返回：
现在我想把点和多边形进行匹配。所以期望的输出是在gdf_points中增加一个附加列，该列包含一个标识符，表示该点与哪个多边形相关联（使用 gdf_polys ['id'] 列）。我的代码速度非常慢，但可以得出正确结果，如下所示：
def id_gen(row):
point=row['geometry']
out=0
for i,poly in shapes_list:
if poly.contains(point):
out=i
break
return out
#shapes_list=gdf_polys['geometry']
shapes_list=[(gdf_polys['id'].iloc[i],gdf_polys['geometry'].iloc[i]) for i in range(len(gdf_polys['geometry']))]
point_list=[]
gdf_points['poly']=gdf_points.apply(id_gen,axis=1)
返回：
x y geometry poly
0 4.865555 1.777419 POINT (4.86555 1.77742) 37
1 6.929483 3.041826 POINT (6.92948 3.04183) 57
2 4.485133 1.492326 POINT (4.48513 1.49233) 37
3 2.889222 6.159370 POINT (2.88922 6.15937) 24
4 2.442262 7.456090 POINT (2.44226 7.45609) 25
... ... ... ... ...
1435 6.414556 5.254309 POINT (6.41456 5.25431) 59
1436 6.409027 4.454615 POINT (6.40903 4.45461) 58
1437 5.763154 2.770337 POINT (5.76315 2.77034) 47
1438 9.613874 1.371165 POINT (9.61387 1.37116) 0
1439 6.013953 3.622011 POINT (6.01395 3.62201) 57
1440 rows × 4 columns
需要注意的是，实际的shapefile形状比这个网格复杂得多。我认为有几个地方可以加速：
- 不需要遍历每一个多边形（随着P的增加，这将变得笨重）
- 使用不同的点-多边形算法。我感觉应该有一种方法可以使用STRtree来实现这个，但目前我只能够返回点而不是索引。
- 向量化数据框操作（避免使用apply函数）
- 可能还有其他未被我考虑到的东西（并行化或其他）
开始基准测试： 网格大小为10，点密度为10（共1440个点）：大约需要180毫秒 网格大小为20，点密度为10（共4840个点）：大约需要2.8秒 网格大小为30，点密度为10（共10240个点）：大约需要12.8秒 网格大小为50，点密度为10（共27040个点）：大约需要1.5分钟 因此我们可以看到这样的规模扩展性很差。