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

Question

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

我有数百万个点(timestep，long，lat，pointID(和数百个不规则的非重叠多边形(vertex_long，vertex_lat，polygonID(。点和多边形格式示例

我想知道每个点包含的多边形是什么。

我可以这样做：

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格式包含与多边形一样多的行和与点一样多的列，其中，在点相对于多边形索引的行(给定的第一个多边形ID->第一行，依此类推(，每一列都为true＝1。

代码的工作速度很慢。当在规则网格中或点半径内定位点时，我成功地尝试实现了KDtree，这大大提高了速度，但我不能对不规则的非重叠多边形做同样或更快的事情。

我看到了一些相关的问题，但与其问点是什么多边形，不如问点是否在多边形内。

你知道吗？

Answer 1

您尝试过Geopandas Spatial加入吗？

使用pip安装软件包pip install geopandas或condaconda 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')

操作中的选项是：

• 相交
• 包含
• 在

在您的情况下，当您连接两个类型为Polygon和Point 的几何体列时，所有列都应该执行相同的操作