点到直线距离的大圆函数无法正常工作。

Question

我需要获取一个经纬度点到一条线的距离，当然需要遵循大圆航线。

我在http://www.movable-type.co.uk/scripts/latlong.html找到了一篇很棒的文章，但代码不能正常工作。要么是我做错了什么，要么是有什么东西漏掉了。以下是相关的函数。如有需要，请查看链接中的其他函数。

    var R = 3961.3
    LatLon.crossTrack = function(lat1, lon1, lat2, lon2, lat3, lon3) {
     var d13 = LatLon.distHaversine(lat1, lon1, lat3, lon3);
     var brng12 = LatLon.bearing(lat1, lon1, lat2, lon2);
     var brng13 = LatLon.bearing(lat1, lon1, lat3, lon3);
     var dXt = Math.asin(Math.sin(d13/R)*Math.sin(brng13-brng12)) * R;
     return dXt;
    } 

纬度/经度1 = -94.127592，41.81762

纬度/经度2 = -94.087257，41.848202

纬度/经度3 = -94.046875，41.791057

这里报告的距离是0.865英里。实际距离是4.29905英里。

有什么线索可以解决这个问题吗？我不是数学家，只是一个老程序员。

- Brad Mathews

---

我也在计算越航距离时遇到了一些麻烦。我想使用你的代码作为示例来检查我的结果，但我不理解你的坐标。当纬度值只定义在-90和+90之间时，你是如何拥有小于-90的纬度值的？ - oschrenk

Answer 1

大多数三角函数需要弧度制。您的角度测量单位是度吗？也许需要使用以下常用公式进行转换：

2 * π 弧度 = 360 度

如果您查看Haversine公式的公式，您会看到以下内容：

（请注意，为了传递给三角函数，角度必须以弧度为单位）。

- duffymo

---

确实是这样！我刚才一时脑抽，没想到bearing方法的返回值是角度。真是个低级错误。 - Brad Mathews

Answer 2

你的函数对于这些坐标是否返回相同的值：

crossTrack(0,0,0,1,0.1,0.5);
crossTrack(0,0,0,1,0.1,0.6);
crossTrack(0,0,0,1,0.1,0.4);

我认为它应该，但我的结果并不是这样。第三个点始终在赤道以北0.1的位置。只有经度变化了，这本不应对结果产生影响。但看起来确实有影响。

- user181345

Answer 3

我尝试了这个点线距离测试，发送了 aalatlon 等等。

private static final double _eQuatorialEarthRadius = 6378.1370D;
private static final double _d2r = (Math.PI / 180D);
private static double PRECISION = 1;





// Haversine Algorithm
// source: https://dev59.com/enRC5IYBdhLWcg3wP-dh

private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
    return  (1000D * HaversineInKM(lat1, long1, lat2, long2));



}

private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
    double dlong = (long2 - long1) * _d2r;
    double dlat = (lat2 - lat1) * _d2r;
    double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
            * Math.pow(Math.sin(dlong / 2D), 2D);
    double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
    double d = _eQuatorialEarthRadius * c;





    return d;
}

// Distance between a point and a line

public static double pointLineDistanceTest(double[] aalatlng,double[] bblatlng,double[]cclatlng){



    double [] a = aalatlng;
    double [] b = bblatlng;
    double [] c = cclatlng;




    double[] nearestNode = nearestPointGreatCircle(a, b, c);
    //        System.out.println("nearest node: " + Double.toString(nearestNode[0])
    + ","+Double.toString(nearestNode[1]));
    double result =  HaversineInM(c[0], c[1], nearestNode[0], nearestNode[1]);

       //        System.out.println("result: " + Double.toString(result));



          return (result);






}

// source: https://dev59.com/93M_5IYBdhLWcg3wn0lO
private static double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
    double[] a_ = toCartsian(a);
    double[] b_ = toCartsian(b);
    double[] c_ = toCartsian(c);

    double[] G = vectorProduct(a_, b_);
    double[] F = vectorProduct(c_, G);
    double[] t = vectorProduct(G, F);

    return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}

@SuppressWarnings("unused")
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
   double[] t= nearestPointGreatCircle(a,b,c);
   if (onSegment(a,b,t))
     return t;

   return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}

 private static boolean onSegment (double[] a, double[] b, double[] t)
   {
     // should be   return distance(a,t)+distance(b,t)==distance(a,b), 
     // but due to rounding errors, we use: 
     return Math.abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
   }


// source: https://dev59.com/LnM_5IYBdhLWcg3w4nfb
private static double[] toCartsian(double[] coord) {
    double[] result = new double[3];
    result[0] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.cos(Math.toRadians(coord[1]));
    result[1] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.sin(Math.toRadians(coord[1]));
    result[2] = _eQuatorialEarthRadius * Math.sin(Math.toRadians(coord[0]));


    return result;
}

private static double[] fromCartsian(double[] coord){
    double[] result = new double[2];
    result[0] = Math.toDegrees(Math.asin(coord[2] / _eQuatorialEarthRadius));
    result[1] = Math.toDegrees(Math.atan2(coord[1], coord[0]));

    return result;
}


// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
    double[] result = new double[3];
    result[0] = a[1] * b[2] - a[2] * b[1];
    result[1] = a[2] * b[0] - a[0] * b[2];
    result[2] = a[0] * b[1] - a[1] * b[0];

    return result;
}

private static double[] normalize(double[] t) {
    double length = Math.sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
    double[] result = new double[3];
    result[0] = t[0]/length;
    result[1] = t[1]/length;
    result[2] = t[2]/length;
    return result;
}

private static double[] multiplyByScalar(double[] normalize, double k) {
    double[] result = new double[3];
    result[0] = normalize[0]*k;
    result[1] = normalize[1]*k;
    result[2] = normalize[2]*k;
    return result;
}

- Tomsmith