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

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

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

我不会直接回答你的问题，而是认为你做错了。经过争论后，我会给出更好的答案。

\n\n

为什么你做错了

\n\n

r 树非常适合两个或三个欧几里得维度的边界框查询。

\n\n

您正在三维空间中弯曲的二维表面上查找经纬度点。结果是您的坐标系将产生奇点和不连续性： 180\xc2\xb0W 与 180\xc2\xb0E 相同，2\xc2\xb0E by 90\xc2\xb0N 接近 2\xc2\xb0W by 90\ xc2\xb0N。r 树不捕获这些类型的东西！

\n\n

但是，即使它们是一个很好的解决方案，您采用lat\xc2\xb1r和lon\xc2\xb1r的想法也会产生一个正方形区域；相反，您可能需要围绕您的点有一个圆形区域。

\n\n

如何正确做事

\n\n

\n

1. 不要将点保留为 lon-lat 格式，而是使用球坐标转换将它们转换为 xyz 格式。现在它们处于 3D 欧几里得空间中，并且不存在奇点或不连续性。

\n

2. 将点放置在三维kd 树中。这使您能够在O(log n)时间内快速提出诸如“到此点的 k 个最近邻是什么？”之类的问题。和“这些点的半径r内有哪些点？” SciPy 附带了一个实现。

\n

3. 对于半径搜索，请从大圆半径转换为弦：这使得 3 空间中的搜索相当于对包裹在球体（在本例中为地球）表面的圆进行的半径搜索。

\n

\n\n

正确执行的代码

\n\n

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

\n\n

#/usr/bin/env python3\n\nimport numpy as np\nimport scipy as sp\nimport scipy.spatial\n\nRearth = 6371\n\n#Generate uniformly-distributed lon-lat points on a sphere\n#See: http://mathworld.wolfram.com/SpherePointPicking.html\ndef GenerateUniformSpherical(num):\n  #Generate random variates\n  pts      = np.random.uniform(low=0, high=1, size=(num,2))\n  #Convert to sphere space\n  pts[:,0] = 2*np.pi*pts[:,0]          #0-360 degrees\n  pts[:,1] = np.arccos(2*pts[:,1]-1)   #0-180 degrees\n  #Convert to degrees\n  pts = np.degrees(pts)\n  #Shift ranges to lon-lat\n  pts[:,0] -= 180\n  pts[:,1] -= 90\n  return pts\n\ndef ConvertToXYZ(lonlat):\n  theta  = np.radians(lonlat[:,0])+np.pi\n  phi    = np.radians(lonlat[:,1])+np.pi/2\n  x      = Rearth*np.cos(theta)*np.sin(phi)\n  y      = Rearth*np.sin(theta)*np.sin(phi)\n  z      = Rearth*np.cos(phi)\n  return np.transpose(np.vstack((x,y,z)))\n\n#Get all points which lie with `r_km` Great Circle kilometres of the query\n#points `qpts`.\ndef GetNeighboursWithinR(qpts,kdtree,r_km):\n  #We need to convert Great Circle kilometres into chord length kilometres in\n  #order to use the kd-tree\n  #See: http://mathworld.wolfram.com/CircularSegment.html\n  angle        = r_km/Rearth\n  chord_length = 2*Rearth*np.sin(angle/2)\n  pts3d        = ConvertToXYZ(qpts)\n  #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\n  #p=2 implies Euclidean distance, eps=0 implies no approximation (slower)\n  return kdtree.query_ball_point(pts3d,chord_length,p=2,eps=0) \n\n\n##############################################################################\n#WARNING! Do NOT alter pts3d or kdtree will malfunction and need to be rebuilt\n##############################################################################\n\n##############################\n#Correctness tests on the North, South, East, and West poles, along with Kolkata\nptsll = np.array([[0,90],[0,-90],[0,0],[-180,0],[88.3639,22.5726]])\npts3d = ConvertToXYZ(ptsll)\nkdtree = sp.spatial.KDTree(pts3d, leafsize=10) #Stick points in kd-tree for fast look-up\n\nqptsll = np.array([[-3,88],[5,-85],[10,10],[-178,3],[175,4]])\nGetNeighboursWithinR(qptsll, kdtree, 2000)\n\n##############################\n#Stress tests\nptsll = GenerateUniformSpherical(100000)    #Generate uniformly-distributed lon-lat points on a sphere\npts3d = ConvertToXYZ(ptsll)                 #Convert points to 3d\n#See: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.html\nkdtree = sp.spatial.KDTree(pts3d, leafsize=10) #Stick points in kd-tree for fast look-up\n\nqptsll = GenerateUniformSpherical(100)      #We\'ll find neighbours near these points\nGetNeighboursWithinR(qptsll, kdtree, 500)\n

\n