有效地将点与几何体（多边形中的点）匹配，以获取点和多边形的大量集合

问题描述 投票：0回答：3

这里有很多关于有效匹配多边形中的点的问题（示例：Here和Here）。这些中感兴趣的主要变量是大量的点 N 和多边形顶点的数量 V。这些都很好且有用，但我正在查看大量的点 N 和多边形 G。这也意味着我的输出将有所不同（我主要看到的输出由落在多边形内的点组成，但在这里我想知道附加到点的多边形）。

我有一个包含大量多边形（数十万）的形状文件。多边形可以接触，但它们之间几乎没有重叠（内部的任何重叠都可能是错误的结果 - 想想人口普查区块组）。我还有一个包含点（数百万）的 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 导入的形状文件）：

# 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

我应该注意到，实际的形状文件将具有比该网格更复杂的形状。我认为有几个地方可以加快速度：

  1. 不必循环遍历每个多边形（随着 P 的增加，这会变得难以处理）
  2. 对多边形内的点使用不同的算法。我觉得应该有一种方法可以使用 STRtree 来做到这一点，但目前我只能返回点（而不是索引）。
  3. 向量化数据帧操作（避免应用）。
  4. 可能还有其他东西不在我的雷达范围内（并行化或其他东西）

开始基准测试： 网格大小为 10，点密度为 10（1440 点）：花费约 180ms 网格大小为 20，点密度为 10（4840 点）：耗时约 2.8s 网格大小为 30，点密度为 10（10240 点）：耗时约 12.8s 网格大小为 50，点密度为 10（27040 点）：大约需要 1.5 分钟 所以我们可以看到它的扩展性很差。

python geopandas shapely
3个回答
3
投票

geopandas 没有将其视为多边形中的质量点，而是有一种在这里有用的空间连接方法。它实际上相当快，至少在这个玩具示例中似乎并没有受到多边形数量的影响（不过我不能排除这可能是由于这些多边形的简单性所致）。

空间连接采用两个地理数据帧并将它们合并在一起。在本例中，我希望将多边形的属性附加到位于其中的点。所以我的代码如下所示：

joined=gpd.sjoin(gdf_points,gdf_polys,how='left',op='within')

退货：

    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 万个多边形，其中大多数多边形比我的玩具示例中使用的几何形状复杂得多。这似乎是一种易于管理的方法，但我有兴趣学习进一步提高效率的方法。

1
投票

听起来您可以通过使用最近的 STRtree 算法来避免迭代所有多边形，如 文档中所写（以及上面关于恢复多边形索引的注释） - 并检查该点是否位于最近的多边形内。 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)
0
投票

利用此处和其他线程上其他答案的启发，对于非常大的点集（数十亿）和可能跨越大面积的多边形来说，最适合我的解决方案是 geopandas.sjoin 和切片数据的组合分成几个部分，借用 osmnx 库的实现。

这个想法是将每个多边形切成规则形状，以便 sjoin 可以使用边界框更有效地进行查找。

import math
import numpy as np
from shapely.geometry import MultiPolygon, LineString
from shapely.ops import split
import geopandas as gp

def quadrat_cut_geometry(geometry, quadrat_width, min_num=1):
    """
    Split a Polygon or MultiPolygon up into sub-polygons of a specified size.
    Parameters
    ----------
    geometry : shapely.geometry.Polygon or shapely.geometry.MultiPolygon
        the geometry to split up into smaller sub-polygons
    quadrat_width : numeric
        the linear width of the quadrats with which to cut up the geometry (in
        the units the geometry is in)
    min_num : int
        the minimum number of linear quadrat lines (e.g., min_num=3 would
        produce a quadrat grid of 4 squares)
    Returns
    -------
    geometry : shapely.geometry.MultiPolygon
    """
    # create n evenly spaced points between the min and max x and y bounds
    west, south, east, north = geometry.bounds
    x_num = math.ceil((east - west) / quadrat_width) + 1
    y_num = math.ceil((north - south) / quadrat_width) + 1
    x_points = np.linspace(west, east, num=max(x_num, min_num))
    y_points = np.linspace(south, north, num=max(y_num, min_num))

    # create a quadrat grid of lines at each of the evenly spaced points
    vertical_lines = [LineString([(x, y_points[0]), (x, y_points[-1])]) for x in x_points]
    horizont_lines = [LineString([(x_points[0], y), (x_points[-1], y)]) for y in y_points]
    lines = vertical_lines + horizont_lines

    # recursively split the geometry by each quadrat line
    for line in lines:
        geometry = MultiPolygon(split(geometry, line))

    return geometry


def assign_region(points, polygons, region_name):
    """
    Assigns region definitions to a GeoDataFrame of points.
    Parameters
    ----------
    points : geopandas.GeoDataFrame
        Input data that contains a geometry field with the point locations
    polygons : geopandas.GeoDataFrame
        Polygon data for each region to be matched
    region_name : str
        Name of the column in the polygons GeoDataFrame that contains the 
        region name which will be assigned to the points data
    Returns
    -------
    geopandas.GeoDataFrame
    """
    cut_polygons = []
    for region, poly in polygons.set_index(region_name)['geometry'].iteritems():
        cut_geom = quadrat_cut_geometry(poly, 0.1)
        cut_geom = [{region_name: region, 'geometry': geom.buffer(1e-14).buffer(0)} for geom in cut_geom]
        cut_polygons.extend(cut_geom)
    cut_polygons = gp.GeoDataFrame(cut_polygons)
    merged = gp.sjoin(points, cut_polygons, how='inner', op='within')
    return merged

使用时，只需将包含多边形的 GeoDataFrame 和另一个包含点的 GeoDataFrame 传递给 

assign_region
函数即可。这将返回点 DataFrame 和一个新列，其中包含多边形 DataFrame 中的标识符列。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.