RTree：计算另一组点的每个点内邻域中的点

Question

为什么这不返回每个邻域（边界框）中的点数？

import geopandas as gpd

def radius(points_neighbour, points_center, new_field_name, r):
    """
    :param points_neighbour:
    :param points_center:
    :param new_field_name: new field_name attached to points_center
    :param r: radius around points_center
    :return:
    """
    sindex = points_neighbour.sindex
    pts_in_neighbour = []
    for i, pt_center in points_center.iterrows():
        nearest_index = list(sindex.intersection((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r)))
        pts_in_this_neighbour = points_neighbour[nearest_index]
        pts_in_neighbour.append(len(pts_in_this_neighbour))
    points_center[new_field_name] = gpd.GeoSeries(pts_in_neighbour)

每个循环都给出相同的结果。

第二个问题，如何找到第 k 个最近的邻居？

有关问题本身的更多信息：

• 我们正在以非常小的规模进行，例如美国华盛顿州或加拿大不列颠哥伦比亚省

• 我们希望尽可能多地利用 geopandas，因为它与 pandas 类似，并且支持空间索引：RTree

• 比如这里的sindex有方法最近，交集等

如果您需要更多信息，请发表评论。这是 GeoPandasBase 类中的代码

@property
def sindex(self):
    if not self._sindex_generated:
        self._generate_sindex()
    return self._sindex

我尝试了 Richard 的示例，但没有成功

def radius(points_neighbour, points_center, new_field_name, r):
    """
    :param points_neighbour:
    :param points_center:
    :param new_field_name: new field_name attached to points_center
    :param r: radius around points_center
    :return:
    """
    sindex = points_neighbour.sindex
    pts_in_neighbour = []
    for i, pt_center in points_center.iterrows():
        pts_in_this_neighbour = 0
        for n in sindex.intersection(((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r))):
            dist = pt_center.distance(points_neighbour['geometry'][n])
            if dist < radius:
                pts_in_this_neighbour = pts_in_this_neighbour + 1
        pts_in_neighbour.append(pts_in_this_neighbour)
    points_center[new_field_name] = gpd.GeoSeries(pts_in_neighbour)

要下载形状文件，请转到 https://catalogue.data.gov.bc.ca/dataset/hellobc-activities-and-attractions-listing 并选择 ArcView 下载

Answer 1

与其直接回答你的问题，我认为你做错了。争论完这个，我会给出一个更好的答案。

为什么你做错了

r-tree 非常适合在两个或三个欧几里得维度中进行边界框查询。

您正在查找在三维空间中弯曲的二维曲面上的经纬度点。结果是您的坐标系将产生奇点和不连续性：180°W 与 180°E 相同，2°E x 90°N 接近 2°W x 90°N。 r-tree 不会捕获这些东西！

但是，即使它们是一个好的解决方案，您采用 lat±r 和 lon±r 的想法也会产生一个正方形区域；相反，您可能希望在您的点周围有一个圆形区域。

如何正确操作

1. 不要将点保留为 lon-lat 格式，而是使用 spherical coordinate conversion 将它们转换为 xyz 格式。现在它们处于 3D 欧几里得空间中，没有奇点或不连续性。

2. 将点放置在三维kd-tree 中。这使您可以在 O(log n) 时间内快速提出诸如“到目前为止的 k 最近邻是什么？”之类的问题。和“这些点的半径 r 内的所有点是什么？” SciPy 附带an implementation。

3. 对于您的半径搜索，将Great Circle radius 转换为chord：这使得在 3 空间中的搜索相当于在包裹到球体表面的圆上进行半径搜索（在这种情况下，地球）。

正确操作的代码

我已经在 Python 中实现了上述内容作为演示。请注意，所有球面点都使用 lon=[-180,180], lat=[-90,90] 方案以 (longitude,latitude)/(x-y) 格式存储。所有 3D 点都以 (x,y,z) 格式存储。

#/usr/bin/env python3

import numpy as np
import scipy as sp
import scipy.spatial

Rearth = 6371

#Generate uniformly-distributed lon-lat points on a sphere
#See: http://mathworld.wolfram.com/SpherePointPicking.html
def GenerateUniformSpherical(num):
  #Generate random variates
  pts      = np.random.uniform(low=0, high=1, size=(num,2))
  #Convert to sphere space
  pts[:,0] = 2*np.pi*pts[:,0]          #0-360 degrees
  pts[:,1] = np.arccos(2*pts[:,1]-1)   #0-180 degrees
  #Convert to degrees
  pts = np.degrees(pts)
  #Shift ranges to lon-lat
  pts[:,0] -= 180
  pts[:,1] -= 90
  return pts

def ConvertToXYZ(lonlat):
  theta  = np.radians(lonlat[:,0])+np.pi
  phi    = np.radians(lonlat[:,1])+np.pi/2
  x      = Rearth*np.cos(theta)*np.sin(phi)
  y      = Rearth*np.sin(theta)*np.sin(phi)
  z      = Rearth*np.cos(phi)
  return np.transpose(np.vstack((x,y,z)))

#Get all points which lie with `r_km` Great Circle kilometres of the query
#points `qpts`.
def GetNeighboursWithinR(qpts,kdtree,r_km):
  #We need to convert Great Circle kilometres into chord length kilometres in
  #order to use the kd-tree
  #See: http://mathworld.wolfram.com/CircularSegment.html
  angle        = r_km/Rearth
  chord_length = 2*Rearth*np.sin(angle/2)
  pts3d        = ConvertToXYZ(qpts)
  #See: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.query_ball_point.html#scipy.spatial.KDTree.query_ball_point
  #p=2 implies Euclidean distance, eps=0 implies no approximation (slower)
  return kdtree.query_ball_point(pts3d,chord_length,p=2,eps=0) 


##############################################################################
#WARNING! Do NOT alter pts3d or kdtree will malfunction and need to be rebuilt
##############################################################################

##############################
#Correctness tests on the North, South, East, and West poles, along with Kolkata
ptsll = np.array([[0,90],[0,-90],[0,0],[-180,0],[88.3639,22.5726]])
pts3d = ConvertToXYZ(ptsll)
kdtree = sp.spatial.KDTree(pts3d, leafsize=10) #Stick points in kd-tree for fast look-up

qptsll = np.array([[-3,88],[5,-85],[10,10],[-178,3],[175,4]])
GetNeighboursWithinR(qptsll, kdtree, 2000)

##############################
#Stress tests
ptsll = GenerateUniformSpherical(100000)    #Generate uniformly-distributed lon-lat points on a sphere
pts3d = ConvertToXYZ(ptsll)                 #Convert points to 3d
#See: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.html
kdtree = sp.spatial.KDTree(pts3d, leafsize=10) #Stick points in kd-tree for fast look-up

qptsll = GenerateUniformSpherical(100)      #We'll find neighbours near these points
GetNeighboursWithinR(qptsll, kdtree, 500)

Answer 2

我附上了代码，经过一些小的修改，应该可以做你想做的事情。

我认为您的问题是出于以下两个原因之一：

1. 您没有正确构建空间索引。您对我的 cmets 的回复表明您并不完全了解空间索引是如何制作的。

2. 您的空间查询的边界框未正确构建。

我将在下面讨论这两种可能性。

构建空间索引

事实证明，空间索引只需键入：

sindex = gpd_df.sindex

魔法。

但是gpd_df.sindex 从哪里得到它的数据呢？它假定数据以shapely 格式存储在名为geometry 的列中。如果您尚未向此类列添加数据，则会引发警告。

数据框的正确初始化如下所示：

#Generate random points throughout Oregon
x = np.random.uniform(low=oregon_xmin, high=oregon_xmax, size=10000)
y = np.random.uniform(low=oregon_ymin, high=oregon_ymax, size=10000)

#Turn the lat-long points into a geodataframe
gpd_df = gpd.GeoDataFrame(data={'x':x, 'y':y})
#Set up point geometries so that we can index the data frame
#Note that I am using x-y points!
gpd_df['geometry'] = gpd_df.apply(lambda row: shapely.geometry.Point((row['x'], row['y'])), axis=1)

#Automagically constructs a spatial index from the `geometry` column
gpd_df.sindex 

在您的问题中看到上述类型的示例代码将有助于诊断您的问题并继续解决它。

由于您没有收到在几何列丢失时引发的极其明显的警告geopandas：

AttributeError: 还没有几何数据集（预计在“几何”列中。

我认为您可能已经正确地完成了这部分。

构造边界框

在你的问题中，你形成了一个像这样的边界框：

nearest_index = list(sindex.intersection((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r)))

事实证明，边界框具有以下形式：

(West, South, East, North)

至少，它们适用于 X-Y 样式点，例如shapely.geometry.Point(Lon,Lat)

在我的代码中，我使用以下内容：

bbox = (cpt.x-radius, cpt.y-radius, cpt.x+radius, cpt.y+radius)

工作示例

将以上内容放在一起可以让我看到这个工作示例。请注意，我还演示了如何按距离对点进行排序，回答您的第二个问题。

#!/usr/bin/env python3

import numpy as np
import numpy.random
import geopandas as gpd
import shapely.geometry
import operator

oregon_xmin = -124.5664
oregon_xmax = -116.4633
oregon_ymin = 41.9920
oregon_ymax = 46.2938

def radius(gpd_df, cpt, radius):
  """
  :param gpd_df: Geopandas dataframe in which to search for points
  :param cpt:    Point about which to search for neighbouring points
  :param radius: Radius about which to search for neighbours
  :return:       List of point indices around the central point, sorted by
                 distance in ascending order
  """
  #Spatial index
  sindex = gpd_df.sindex
  #Bounding box of rtree search (West, South, East, North)
  bbox = (cpt.x-radius, cpt.y-radius, cpt.x+radius, cpt.y+radius)
  #Potential neighbours
  good = []
  for n in sindex.intersection(bbox):
    dist = cpt.distance(gpd_df['geometry'][n])
    if dist<radius:
      good.append((dist,n))
  #Sort list in ascending order by `dist`, then `n`
  good.sort() 
  #Return only the neighbour indices, sorted by distance in ascending order
  return [x[1] for x in good]

#Generate random points throughout Oregon
x = np.random.uniform(low=oregon_xmin, high=oregon_xmax, size=10000)
y = np.random.uniform(low=oregon_ymin, high=oregon_ymax, size=10000)

#Turn the lat-long points into a geodataframe
gpd_df = gpd.GeoDataFrame(data={'x':x, 'y':y})
#Set up point geometries so that we can index the data frame
gpd_df['geometry'] = gpd_df.apply(lambda row: shapely.geometry.Point((row['x'], row['y'])), axis=1)

#The 'x' and 'y' columns are now stored as part of the geometry, so we remove
#their columns in order to save space
del gpd_df['x']
del gpd_df['y']

for i, row in gpd_df.iterrows():
  neighbours = radius(gpd_df,row['geometry'],0.5)
  print(neighbours)
  #Use len(neighbours) here to construct a new row for the data frame

（我在 cmets 中要求的是与上述类似的代码，但它可以说明您的问题。注意使用 random 简洁地生成用于实验的数据集。）