找出每个点在哪个多边形中

Question

我是 Python 新手，所以我为基本的编程技能道歉，我知道我使用了太多的“循环 for”（来自 Matlab，它拖累了我）。

我有数百万个点（timestep、long、lat、pointID）和数百个不规则的非重叠多边形（vertex_long、vertex_lat、polygonID）。points and polygons format sample

我想知道每个点包含哪个多边形。

我可以这样做：

from matplotlib import path
def inpolygon(lon_point, lat_point, lon_poly, lat_poly):
   shape = lon_point.shape
   lon_point = lon_point.reshape(-1)
   lat_point = lat_point.reshape(-1)
   lon_poly = lon_poly.values.reshape(-1)
   lat_poly = lat_poly.values.reshape(-1)
   points = [(lon_point[i], lat_point[i]) for i in range(lon_point.shape[0])]
   polys = path.Path([(lon_poly[i], lat_poly[i]) for i in range(lon_poly.shape[0])])
   return polys.contains_points(points).reshape(shape)

然后

import numpy as np
import pandas as pd
Areas_Lon = Areas.iloc[:,0]
Areas_Lat = Areas.iloc[:,1]
Areas_ID  = Areas.iloc[:,2]
Unique_Areas = np.unique(Areas_ID)

Areas_true=np.zeros((Areas_ID.shape[0],Unique_Areas.shape[0]))
for i in range(Areas_ID.shape[0]):
    for ii in range(Unique_Areas.shape[0]):
        Areas_true[i,ii]=(Areas_ID[i]==Unique_Areas[ii])

Areas_Lon_Vertex=np.zeros(Unique_Areas.shape[0],dtype=object)
Areas_Lat_Vertex=np.zeros(Unique_Areas.shape[0],dtype=object)
for i in range(Unique_Areas.shape[0]):
    Areas_Lon_Vertex[i]=(Areas_Lon[(Areas_true[:,i]==1)])
    Areas_Lat_Vertex[i]=(Areas_Lat[(Areas_true[:,i]==1)])

import f_inpolygon as inpolygon
Areas_in=np.zeros((Unique_Areas.shape[0],Points.shape[0]))
for i in range (Unique_Areas.shape[0]):
    for ii in range (PT.shape[0]):
        Areas_in[i,ii]=(inpolygon.inpolygon(Points[ii,2], Points[ii,3], Areas_Lon_Vertex[i], Areas_Lat_Vertex[i]))
        

这样，最终结果 Areas_in Areas_in format 包含与多边形一样多的行和与点一样多的列，其中每一列在点相对于多边形索引的行处为 true=1（第一个给定的多边形 ID -- > 第一行，依此类推）。

代码可以正常工作，但是对于它应该做的事情来说非常缓慢。当在规则网格或点半径内定位点时，我已成功尝试实现 KDtree，这显着提高了速度，但对于不规则的非重叠多边形，我无法做到相同或更快。

我看到了一些相关的问题，但不是询问一个点是什么多边形，而是关于一个点是否在多边形内。

有什么想法吗？

Answer 1

您尝试过 Geopandas 空间连接吗？

使用 pip 安装包 pip install geopandas 或康达 conda install -c conda-forge geopandas

那么您应该能够将数据读取为 GeoDataframe

import geopandas 

df = geopandas.read_file("file_name1.csv") # you can read shp files too.
right_df = geopandas.read_file("file_name2.csv") # you can read shp files too.

# Convert into geometry column 
geometry = [Point(xy) for xy in zip(df['longitude'], df['latitude'])] # Coordinate reference system : WGS84

crs = {'init': 'epsg:4326'}
# Creating a Geographic data frame 
left_df = geopandas.GeoDataFrame(df, crs=crs, geometry=geometry)

然后你就可以申请sjoin了

jdf = geopandas.sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuffix='right')

op 中的选项是：

• 相交
• 包含
• 内

当你加入两个多边形和点类型的几何列时，所有的都应该做同样的事情