Yelper 之查找离用户最近的商店

Notebook: Calc_Nearest_Businesses

---

当给出用户的地理坐标（纬度，精度）时，如何快速获取到离该用户最近的商店？

KNN（K Nearest Neighbors）算法跟我们的需求有点像，都涉及到求离某一点最近的点，只不过它一般是用来基于这些最近的点，投票决定该点的类别。咋看起来，KNN 不能满足我们的需求，但仔细一想，KNN 跟我们的需求本质上是一样的，都需要计算距离（多种定义，常见的是欧氏距离）。当我们查找到离用户最远距离为 D 的所有商店时，K 这个值也就确定了；当我们找到 K 个最近的点时，离用户最远的距离 D 也就确定了。所以，接下来我们就用 KNN 算法来解决我们的问题。

不过，这里还有一个问题就是用户和商店（Business）的地理坐标都是经纬度，直接计算直线距离是不行的，而且这些经纬度能够传入到 KNN 算法中计算。这些处理我们在后文中进行说明。

KNN 算法简要说明

KNN 算法算是一个比较简单算法，它没有显示的学习过程，也很容易理解。例如，对于下图：
链接

中央的绿色点是待分类点（分为蓝色方形或红色三角形），当：

• K = 3 时，距离绿色点最近的 3 个点中有 2 个是红色三角形，1 个是蓝色方形，所以该绿色点被分为红色三角形类别
• K = 5 时，距离绿色点最近的 5 个点中有 2 个是红色三角形，3 个是蓝色方形，所以该绿色点被分为蓝色方形类别
• 其余以此类推

虽然该算法简单，但高效计算就不是一件容易的事了。简单拿一个点跟所有点都计算一次距离在大数据情况下是不大可行的。对此，KNN 使用 KD 树来解决这个问题。关于 KD 树详细可参考博客 KD 树算法之思路篇、KD 树算法之详细篇。下文也将使用 KD 树来查找离用户最近的商店。

查找离用户最近的商店

这部分主要是参考了 KD-Tree 2 Real Distance、How to Convert KD-Tree Distance from Query to Kilometers 两篇资料。前篇资料是在后者基础上写成的。这部分的代码我们就直接摘录自前篇文档。

坐标转换

因为我们不能直接计算地球上两点的直线距离，所以需要先将经纬度转换为地心坐标系下的坐标。下边是维基关于地心坐标系的描述：

地心地固坐标系（Earth-Centered, Earth-Fixed，简称ECEF）简称地心坐标系，是一种以地心为原点的地固坐标系（也称地球坐标系），是一种笛卡儿坐标系。原点 O (0,0,0)为地球质心，z 轴与地轴平行指向北极点，x 轴指向本初子午线与赤道的交点，y 轴垂直于xOz平面(即东经90度与赤道的交点)构成右手坐标系。

如何转换经纬度为地心坐标系下坐标看下图就比较清晰了：
图片链接

写成代码如下：

def to_Cartesian(lat, lng):
    '''
    function to convert latitude and longitude to 3D cartesian coordinates
    '''
    R = 6371 # radius of the Earth in kilometers

    x = R * cos(lat) * cos(lng)
    y = R * cos(lat) * sin(lng)
    z = R * sin(lat)
    return x, y, z

距离转换

地理坐标转换后，还需要把地表距离转换为地心坐标系下的距离，这样才可以利用 KD 树来计算指定范围内的商店，代码如下：

def rad2deg(rad):
    '''
    function to convert radian to degree
    '''
    degree = rad/2/np.pi * 360
    return(degree)

def kmToDIST(x):
    '''
    function to convert real distance in km to cartesian distance 
    '''
    R = 6371 # earth radius
    gamma = 2*np.arcsin(x/2./R) 
    
    dist = 2*R*rad2deg(sin(gamma / 2.))
    return(dist)

进行反向转换的代码如下：

def deg2rad(degree):
    '''
    function to convert degree to radian
    '''
    rad = degree * 2*np.pi / 360
    return(rad)

def distToKM(x):
    '''
    function to convert cartesian distance to real distance in km
    '''
    R = 6371 # earth radius
    gamma = 2*np.arcsin(deg2rad(x/(2*R))) # compute the angle of the isosceles triangle
    dist = 2*R*sin(gamma/2) # compute the side of the triangle
    return(dist)

查找离用户最近的商店

完成以上两步后，就可以用 KD 树来计算指定范围内的商店了，代码如下：

# convert the grid points and reference points to cartesian coordinates
x, y, z = zip(*map(to_Cartesian, lats_1d, lons_1d))
x_ref, y_ref, z_ref = to_Cartesian(lat_0, lon_0)

# convert the 30 km to cartesian coordinates distance
dist = kmToDIST(30)	# 原始代码这里指定了范围为 30km，可以根据情况修改

# create the KD-tree using the 3D cartesian coordinates
coordinates = list(zip(x, y, z))
tree = spatial.cKDTree(coordinates)

# get all the points within 30 km from the reference point
ix = tree.query_ball_point((x_ref, y_ref, z_ref), dist)

得到的结果 ix 就是经纬度列表的下标。

详细处理过程请参考 Notebook Calc_Nearest_Businesses