为什么我的代码计算转换坐标之间的距离不正确？答案

【问题标题】：Why is my code calculating an incorrect distance between transformed coordinates?为什么我的代码计算转换坐标之间的距离不正确？
【发布时间】：2019-11-01 17:53:56
【问题描述】：

我有一个大的点文件，我正在尝试找出这些点与另一组点之间的距离。最初，我使用 geopandas 的 to_crs 函数来转换 crs，以便在我执行 df.distance(point) 时可以获得以米为单位的准确距离测量值。但是，由于文件很大，仅仅转换文件的crs就花了很长时间。该代码运行了 2 个小时，但仍未完成转换。因此，我改用了这段代码。

inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')

for index, row in demand_indo_gdf.iterrows():
    o = Point(row['origin_longitude'], row['origin_latitude'])
    o_proj = Point(transform(inProj, outProj, o.x, o.y))

    for i, r in bus_indo_gdf.iterrows():
        stop = r['geometry']
        stop_proj = Point(transform(inProj, outProj, stop.x, stop.y))
        print ('distance:', o_proj.distance(stop_proj), '\n\n')

我认为单独转换 crs 并执行我的分析可能会更快。对于这组点：

o = (106.901024 -6.229162)
stop = (106.804 -6.21861)

我将这个 EPSG 4326 坐标转换为本地投影 EPSG 4808，得到了这个：

o_proj = (0.09183386384156803 -6.229330112968891)
stop_proj = (-0.005201753272169649 -6.218776788266844)

这给出了 0.09760780527657992 的距离测量值。谷歌地图给了我一个距离测量，坐标o和stop，10.79公里。看起来我的代码的距离测量给出的答案比实际距离小 10^-3 倍。为什么呢？我的代码正确吗？

【问题讨论】：

标签： python python-3.x distance geopandas pyproj

【解决方案1】：

您使用的坐标以度为单位，而不是米。

您正在使用的Points 类可能并不关心这一点，而是计算它们之间的笛卡尔距离。

请使用Haversine 等式，例如来自here的那个：

from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

那么你的代码返回正确的结果：

from pyproj import Proj, transform

inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')

o = (106.901024, -6.229162)
stop = (106.804, -6.21861)

o_proj = transform(inProj, outProj, o[0], o[1])
stop_proj = transform(inProj, outProj, stop[0], stop[1])

print(haversine(o_proj[0], o_proj[1], stop_proj[0], stop_proj[1]))

【讨论】：