Блог

Я новичок в 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 содержит столько строк, сколько полигонов, и столько столбцов, сколько точек, где каждый столбец равен true = 1 в строке, где точка находится относительно индекса полигона (1-й заданный идентификатор полигона -> 1-я строка и так далее).

Код работает, но очень медленно для того, что он должен делать. При нахождении точек в обычной сетке или в радиусе точки я успешно попытался реализовать KDTree, что значительно увеличивает скорость, но я не могу сделать то же самое или что-то еще быстрее для неправильных неперекрывающихся полигонов.

Я видел несколько связанных вопросов, но вместо того, чтобы спрашивать, из каких полигонов состоит точка, речь шла о том, находится ли точка внутри многоугольника или нет.

Любая идея, пожалуйста?

Вы пробовали пространственное объединение Geopandas?

установите пакет, используя pip pip install geopandas или conda conda install -c conda-forge geopandas

затем вы должны иметь возможность считывать данные как фрейм геоданных

 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 является:

• пересекает
• содержит
• внутри

Все должно делать то же самое в вашем случае, когда вы соединяете два столбца геометрии типа Polygon и Point