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

Question

我有一个大文件的点，我试图找到这些点和另一组点之间的距离。最初我使用地质标准的^{cd1>}函数来转换crs，这样当我做^{cd2>}时，我就可以得到一个以米为单位的精确距离测量。但是，由于文件很大，所以转换文件的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并进行分析可能会更快。对于这组点：

^{pr2}$

我将这个EPSG 4326坐标转换成局部投影EPSG 4808，得到了：

^{pr3}$

这给出了0.09760780527657992的距离测量值。Google地图给了我一个距离测量，坐标为10.79km，坐标为^{cd3>}和^{cd4>}。看起来，我的代码的距离度量给出的答案是实际距离的10^-3倍。为什么这样？我的代码正确吗？

Answer 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

然后您的代码将返回正确的结果：

^{pr2}$