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 有 nearest、intersection 等方法。

如果需要更多信息，请留言。这是 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进行下载。

- ZHU

---

你能否发布一下生成rtree的代码？ - Richard

@Richard points_neighbour.sindex 这是你想要的吗？ - ZHU

是的，那应该就是了。 - Richard

请问您是否有“points_neighbour.sindex”代码？ - Richard

你需要编辑你的问题，包括你正在使用的代码；否则，你会浪费别人的时间，让他们试图解密你发布的链接。 - Richard

Answer 1

不是直接回答你的问题，我认为你做错了。在提出这个观点后，我会给出更好的答案。

为什么你做错了

r-tree非常适用于在二维或三维欧几里得空间中进行边界框查询。

你正在查找位于三维空间中弯曲的二维表面上的经纬度点。结果是，你的坐标系将产生奇点和不连续性：180°W与180°E相同，2°E by 90°N接近于2°W by 90°N。r-tree无法捕获这些类型的问题！

但即使它们是一个好的解决方案，你的想法是采取lat±r和lon±r，这将产生一个正方形区域；相反，你可能需要一个围绕你的点的圆形区域。

如何正确地做

1. 将点从经纬度格式转换为xyz格式，使用球坐标转换。现在它们在三维欧几里得空间中，没有奇点或不连续性。

2. 将点放置在三维kd-tree中。这允许您快速地在O(log n)时间内询问问题，例如“这个点的k个最近邻是什么？”和“距离这个点小于半径r的所有点是什么？”SciPy带有一个实现。

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)

- Richard

---

我实际上是在一个较小的范围内进行，比如说华盛顿州。所以我认为奇异点并不重要？ - ZHU

我想你知道我在问什么。你能给另一个实现吗？尝试在geopandas中利用RTree索引。 - ZHU

已经编辑过了。还有其他需要了解的吗？ - ZHU

你最后一次编辑这个答案是什么时候？它看起来和我上次看到它的时候不同。 - ZHU

我已经苦恼了几周如何对球体进行空间索引的问题（我已经阅读了几篇关于工作但复杂的系统的论文，这些系统最终会比我所有其他代码加起来还要重），当我突然想起上周看到过这个建议时，感觉非常惊喜。优雅、简单、实用，使用现成的库，最重要的是没有边缘情况。谢谢。 - dmckee --- ex-moderator kitten

Answer 2

我已经附上了一些代码，稍加修改就可以完成你想要的功能。

我认为你的问题可能出现在以下两个方面之一：

1. 你没有正确构建空间索引。你对我的评论回复表明你并不完全了解空间索引是如何生成的。 2. 空间查询的边界框没有正确构建。

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

构建空间索引

事实上，只需输入以下命令即可构建空间索引：

sindex = gpd_df.sindex

神奇。

但是gpd_df.sindex从哪里获取数据呢？它假设数据存储在一个名为geometry的列中，格式为shapely。如果您没有向这样的列添加数据，它将引发警告。

正确初始化数据框应该如下：

#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: No geometry data set yet (expected in column 'geometry'.

我想你可能已经正确完成了这部分内容。

构建边界框

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

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

在评论中我一直在要求的是像上述代码那样的代码，它展示了你的问题。请注意使用random来简洁地生成一个数据集用于实验。

- Richard

---

@ZHU：这并不是特别具体的描述。我提供的示例是否有效？如果无效，可能是您的设置出了问题：在我的端上它是可以工作的。如果有效，那么您需要弄清楚您正在做什么不同。 - Richard

看到我的编辑后的问题，其中包括了你的实现，但我得到的全部计数都是0。而当我使用我的原始代码时，总共计数了41231个点。 - ZHU

在这里，您可以下载数据集 https://catalogue.data.gov.bc.ca/dataset/hellobc-activities-and-attractions-listing 除了geopandas库之外，没有使用其他额外的库。 - ZHU

我遇到的问题仅限于这个数据集。我认为随机生成点是没有意义的。 - ZHU

@ZHU：问题可能出在很多地方，包括你的数据。我的示例展示了如何使用随机生成的点与你提出的算法。请构建一个使用随机生成的点复制你的问题的示例。如果你无法这样做，那么问题可能出在你的数据上。如果你能够这样做，那么问题就在你的代码中，因为我们将同时查看同一件事情，所以很容易找到。 - Richard