gis - RTree : Count points in the neighbourhoods within each point of another set of points

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

我试过理查德的例子，但没有用

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 树非常适合在两个或三个欧几里得维度中进行边界框查询。

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

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

如何正确做

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

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

对于您的半径搜索，从 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)