对于点和多边形的大量集合，有效地将点与几何（多边形中的点）匹配
Efficiently match point to geometry (point in poly) for large collections of both points and polygons
这里有很多关于有效匹配多边形点的问题（示例： and ）。这些中感兴趣的主要变量是大量的点 N 和多边形顶点的数量 V。这些都很好而且有用，但我正在查看大量的点 N 和多边形 G。这也意味着我的输出将不同（我主要看到输出由落在多边形内的点组成，但在这里我想知道附加到一个点的多边形）。
我有一个包含大量多边形（数十万个）的 shapefile。多边形可以接触，但它们之间几乎没有重叠（任何内部重叠都是错误的结果——想想人口普查区块组）。我还有一个带有点（数百万）的 csv，我想根据点落入的多边形对这些点进行分类（如果有的话）。有些可能不会落入多边形（继续我的例子，想想海洋上的点）。下面我搭建一个toy example来看看这个问题。
设置：
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)
Returns:
我现在想用多边形匹配点。因此，所需的输出将是 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)
其中 return 个：
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 来做到这一点，但目前我只能 return 点（而不是索引）。
- 矢量化数据帧操作（避免应用）。
- 可能是我没有注意到的其他东西（并行化或其他东西）
开始基准测试：
网格大小为 10，点密度为 10（1440 点）：大约需要 180 毫秒
网格大小为 20，点密度为 10（4840 点）：大约需要 2.8 秒
网格大小为 30，点密度为 10（10240 点）：大约需要 12.8 秒
网格大小为 50，点密度为 10（27040 点）：大约需要 1.5 分钟
所以我们可以看到这个比例很差。
听起来您可以通过使用最近的 STRtree 算法避免遍历所有多边形，如 the documentation 中所写（以及上面关于恢复多边形索引的注释）- 并检查该点是否位于最近的多边形。 IE。像
from shapely.strtree import STRtree
#... coords is the list of shapely points and poly_list is the list of polygons ...
#to recover the polygon id, use their unique python id.
poly_id = dict((id(poly), i) for i, poly in enumerate(poly_list))
#form stretree of polygons
poly_tree = STRtree(poly_list)
pt_to_id = []
#fill pt_to_id with the nearest polygon if it contains the given point. If the point is within no polygon write None.
for c in coords:
near = poly_tree.nearest(c)
if near.contains(c):
pt_to_id.append(poly_id[id(near)])
else:
pt_to_id.append(None)
geopandas 没有将其视为多边形中的大量点，而是采用了一种在这里很有用的空间连接方法。它实际上非常快，至少在这个玩具示例中似乎并没有受到多边形数量的影响（我不能排除这可能是由于这些多边形的简单性）。
Spatial join 获取两个地理数据框并将它们合并在一起。在这种情况下，我想要附加到位于其中的点的多边形的属性。所以我的代码看起来像这样：
joined=gpd.sjoin(gdf_points,gdf_polys,how='left',op='within')
Returns:
x y geometry poly index_right id numeric string included
0 18.651358 26.920261 POINT (18.65136 26.92026) 908 908.0 908.0 0.0 foo True
1 38.577101 1.505424 POINT (38.57710 1.50542) 1863 1863.0 1863.0 0.0 foo True
2 15.430436 15.543219 POINT (15.43044 15.54322) 750 750.0 750.0 0.0 foo True
3 44.928141 7.726345 POINT (44.92814 7.72635) 2163 2163.0 2163.0 0.0 foo True
4 34.259632 5.373809 POINT (34.25963 5.37381) 1671 1671.0 1671.0 0.0 foo True
... ... ... ... ... ... ... ... ... ...
27035 32.386086 23.440186 POINT (32.38609 23.44019) 1591 1591.0 1591.0 0.0 foo True
27036 7.569414 1.836633 POINT (7.56941 1.83663) 344 344.0 344.0 0.0 foo True
27037 1.141440 34.739388 POINT (1.14144 34.73939) 83 83.0 83.0 0.0 foo True
27038 -0.770784 14.027607 POINT (-0.77078 14.02761) 0 NaN NaN NaN NaN NaN
27039 12.695803 33.405048 POINT (12.69580 33.40505) 621 621.0 621.0 0.0 foo True
与我的初始代码相比，这非常快。 运行 我测试的最大尺寸（27k 点）耗时不到 60 毫秒（与之前代码的 1.5 分钟相比）。放大到我的一些实际工作，100 万个点只用了 13 秒多一点的时间来匹配不到 20 万个多边形，其中大部分比我的玩具示例中使用的几何体复杂得多。这似乎是一种易于管理的方法，但我有兴趣学习进一步提高效率的方法。
这里有很多关于有效匹配多边形点的问题（示例：
我有一个包含大量多边形（数十万个）的 shapefile。多边形可以接触，但它们之间几乎没有重叠（任何内部重叠都是错误的结果——想想人口普查区块组）。我还有一个带有点（数百万）的 csv，我想根据点落入的多边形对这些点进行分类（如果有的话）。有些可能不会落入多边形（继续我的例子，想想海洋上的点）。下面我搭建一个toy example来看看这个问题。
设置：
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)
Returns:
我现在想用多边形匹配点。因此，所需的输出将是 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)
其中 return 个：
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 来做到这一点，但目前我只能 return 点（而不是索引）。
- 矢量化数据帧操作（避免应用）。
- 可能是我没有注意到的其他东西（并行化或其他东西）
开始基准测试： 网格大小为 10，点密度为 10（1440 点）：大约需要 180 毫秒 网格大小为 20，点密度为 10（4840 点）：大约需要 2.8 秒 网格大小为 30，点密度为 10（10240 点）：大约需要 12.8 秒 网格大小为 50，点密度为 10（27040 点）：大约需要 1.5 分钟 所以我们可以看到这个比例很差。
听起来您可以通过使用最近的 STRtree 算法避免遍历所有多边形，如 the documentation 中所写（以及上面关于恢复多边形索引的注释）- 并检查该点是否位于最近的多边形。 IE。像
from shapely.strtree import STRtree
#... coords is the list of shapely points and poly_list is the list of polygons ...
#to recover the polygon id, use their unique python id.
poly_id = dict((id(poly), i) for i, poly in enumerate(poly_list))
#form stretree of polygons
poly_tree = STRtree(poly_list)
pt_to_id = []
#fill pt_to_id with the nearest polygon if it contains the given point. If the point is within no polygon write None.
for c in coords:
near = poly_tree.nearest(c)
if near.contains(c):
pt_to_id.append(poly_id[id(near)])
else:
pt_to_id.append(None)
geopandas 没有将其视为多边形中的大量点，而是采用了一种在这里很有用的空间连接方法。它实际上非常快，至少在这个玩具示例中似乎并没有受到多边形数量的影响（我不能排除这可能是由于这些多边形的简单性）。
Spatial join 获取两个地理数据框并将它们合并在一起。在这种情况下，我想要附加到位于其中的点的多边形的属性。所以我的代码看起来像这样：
joined=gpd.sjoin(gdf_points,gdf_polys,how='left',op='within')
Returns:
x y geometry poly index_right id numeric string included
0 18.651358 26.920261 POINT (18.65136 26.92026) 908 908.0 908.0 0.0 foo True
1 38.577101 1.505424 POINT (38.57710 1.50542) 1863 1863.0 1863.0 0.0 foo True
2 15.430436 15.543219 POINT (15.43044 15.54322) 750 750.0 750.0 0.0 foo True
3 44.928141 7.726345 POINT (44.92814 7.72635) 2163 2163.0 2163.0 0.0 foo True
4 34.259632 5.373809 POINT (34.25963 5.37381) 1671 1671.0 1671.0 0.0 foo True
... ... ... ... ... ... ... ... ... ...
27035 32.386086 23.440186 POINT (32.38609 23.44019) 1591 1591.0 1591.0 0.0 foo True
27036 7.569414 1.836633 POINT (7.56941 1.83663) 344 344.0 344.0 0.0 foo True
27037 1.141440 34.739388 POINT (1.14144 34.73939) 83 83.0 83.0 0.0 foo True
27038 -0.770784 14.027607 POINT (-0.77078 14.02761) 0 NaN NaN NaN NaN NaN
27039 12.695803 33.405048 POINT (12.69580 33.40505) 621 621.0 621.0 0.0 foo True
与我的初始代码相比，这非常快。 运行 我测试的最大尺寸（27k 点）耗时不到 60 毫秒（与之前代码的 1.5 分钟相比）。放大到我的一些实际工作，100 万个点只用了 13 秒多一点的时间来匹配不到 20 万个多边形，其中大部分比我的玩具示例中使用的几何体复杂得多。这似乎是一种易于管理的方法，但我有兴趣学习进一步提高效率的方法。