【发布时间】：2021-07-14 03:09:07
【问题描述】：
这里有很多关于有效匹配多边形中的点的问题（例如：Here 和 Here）。这些中感兴趣的主要变量是大量的点 N 和多边形顶点的数量 V。这些都很好且有用，但我正在查看大量的点 N 和多边形 G。这也意味着我的输出将有所不同（我主要看到的输出由落在多边形内的点组成，但在这里我想知道附加到一个点的多边形）。
我有一个包含大量多边形（数十万个）的 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 来执行此操作，但目前我只能返回点（而不是索引）。
- 矢量化数据帧操作（避免应用）。
- 可能是我不在意的其他事情（并行化或其他事情）
开始基准测试： 网格大小为 10，点密度为 10（1440 个点）：耗时约 180ms 网格大小为 20，点密度为 10（4840 个点）：大约需要 2.8 秒 网格大小为 30，点密度为 10（10240 个点）：耗时约 12.8s 网格大小为 50，点密度为 10（27040 个点）：大约需要 1.5 分钟 所以我们可以看到这个规模很差。
【问题讨论】：