计算Cython中两个地理代码之间的距离

Question

我的网络服务当前从 google API 检索地址的地理编码。比较两个地理代码之间的距离时，我尝试使用此处描述的算法 http://andrew.hedges.name/experiments/haversine/ 不幸的是，它 returns 值不准确。

我有两个版本

cdef double M_PI = 3.141592653589793;
cdef double rads = (180/M_PI)

cdef double lon1 = a.geoCode[1];
cdef double lon2 = b.geoCode[1];

cdef double lat1 = a.geoCode[0];
cdef double lat2 = b.geoCode[0];

cdef double dlat = lon2 - lon1
cdef double dlon = lat2 - lat1

cdef double x  = (pow(sin(dlat/2),2)) + cos(lat1)*cos(lat2)*pow(sin(dlon/2),2);

cdef double c = 2.0 * atan2(sqrt(x), sqrt(1-x));

cdef double d = 6370.0 * c;

return d

处理的值的示例是

lat lon a 38.898556 -77.037852
lat lon b 38.897147 -77.043934
Difference in lat -0.00608199999999
Difference in long -0.001409
x 9.31324375062e-06
sqrt(x) 0.00305176076235
sqrt(1-x) 0.999995343367
atan2(sqrt(x),sqrt(1-x)) 1.5677445613
c 0.00610353099867
distance in kiloemters 38.8794924615

和

cdef double M_PI = 3.141592653589793;
cdef double rads = (180/M_PI)

cdef double lon1 = a.geoCode[1]*rads;
cdef double lon2 = b.geoCode[1]*rads;

cdef double lat1 = a.geoCode[0]*rads;
cdef double lat2 = b.geoCode[0]*rads;

cdef double dlat = lon2 - lon1
cdef double dlon = lat2 - lat1

cdef double x  = (pow(sin(dlat/2),2)) + cos(lat1)*cos(lat2)*pow(sin(dlon/2),2);

cdef double c = 2.0 * atan2(sqrt(x), sqrt(1-x));

cdef double d = 6370.0 * c;
return d

哪些进程

lat lon a 2228.72308795 -4413.94378235
lat lon b 2228.6423582 -4414.29225528
Difference in lat -0.348472930998
Difference in long -0.0807297533338
x 0.0301717372754
sqrt(x) 0.173700136083
sqrt(1-x) 0.984798589928
atan2(sqrt(x),sqrt(1-x)) 1.39621064139
c 0.349171370807
distance in kiloemters 2224.22163204

问题是实际答案应该是0.549公里

我不确定我所做的与此公式有何不同

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c (where R is the radius of the Earth)

Answer 1

我认为弧度转换是倒转的（180/π 而不是 π/180）。其他好像没问题。

我在自己的一个应用中使用了下面的代码，效果不错（是Objective-C，但是类似）：

// Credit:
// http://www.movable-type.co.uk/scripts/latlong.html

double R = 6371000; // metres
double phi1 = lat1 * M_PI / 180.0;
double phi2 = lat2 * M_PI / 180.0;
double deltaphi = (lat2-lat1) * M_PI / 180.0;
double deltalambda = (long2-long1) * M_PI / 180.0;

double a = sin(deltaphi/2) * sin(deltaphi/2) +
    cos(phi1) * cos(phi2) *
    sin(deltalambda/2) * sin(deltalambda/2);
double c = 2 * atan2(sqrt(a), sqrt(1.0 - a));
double d = R * c;
return d;

我的代码没有使用 pow() 但你的代码在这方面是正确的，因为它避免了多余的 sin() 计算。

计算Cython中两个地理代码之间的距离

Calculating distance between two geoCodes in Cython

c

python

google-maps

geocoding

cython