从经纬度点到小弧段的距离

Question

从经纬度点到小弧段的距离

matlabgpsgisdistancelatitude-longitude

23

我需要计算从一个经纬度GPS点P到由另外两个经纬度GPS点A和B描述的线段的最短距离。

"越轨距离"可以帮助我计算P与由A和B描述的大圆之间的最短距离。

然而，这不是我想要的。我需要的是P和A-B的线段之间的距离，而不是整个大圆。

我已经使用了来自http://www.movable-type.co.uk/scripts/latlong.html的以下实现。

Formula:    dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
where:
δ13 is (angular) distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius

以下图片有助于说明我正在尝试解决的问题：

在第一张图片中，由绿色线表示的交叉距离是正确的，确实是到线段AB的最短距离。

在第二张图片中，显示了交叉距离的问题。在这种情况下，我希望最短距离是简单距离AP，但交叉距离给出的是红色线指示的距离。

如何更改我的算法以考虑这一点，或者检查点X是否位于AB内？是否可能进行计算？还是迭代是唯一可能的（昂贵）解决方案？（沿AB取N个点，然后计算从P到所有这些点的最小距离）

为了简化起见，图像中所有线都是直线。实际上，这些是大圆上的小弧。

- ChrisDekker

在该应用程序中，最大距离是多少？几百米还是高达几千公里？ - AlexWien

线段是（驾驶）路线中的航点。我会说它们之间的距离大约在100-1000米左右。 - ChrisDekker

你已经将交叉轨道距离翻译成Matlab了吗？如果是的话，能否将代码放入您的问题中？ - Daniel

在平面上，您只需检查ABP三角形中的所有角度是否<= 90度。是否有类似的东西可以从球面三角学中应用？ - Chad Kennedy

当然 Chad，一开始也是这么想的，但我可以想到几种情况，如果使用这种方法将不会产生准确的结果，包括跨越不同半球的点和弧段很小的情况。 - ChrisDekker

显示剩余2条评论

7个回答

3

并添加了Sga实现的Python翻译：

    def bear(latA, lonA, latB, lonB):
        # BEAR Finds the bearing from one lat / lon point to another.
        return math.atan2(
            math.sin(lonB - lonA) * math.cos(latB),
            math.cos(latA) * math.sin(latB) - math.sin(latA) * math.cos(latB) * math.cos(lonB - lonA)
        )


    def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
        lat1 = math.radians(lat1)
        lat2 = math.radians(lat2)
        lat3 = math.radians(lat3)
        lon1 = math.radians(lon1)
        lon2 = math.radians(lon2)
        lon3 = math.radians(lon3)
        R = 6378137

        bear12 = bear(lat1, lon1, lat2, lon2)
        bear13 = bear(lat1, lon1, lat3, lon3)
        dis13 = distance( (lat1, lon1), (lat3, lon3)).meters

        # Is relative bearing obtuse?
        if math.fabs(bear13 - bear12) > (math.pi / 2):
            return dis13

        # Find the cross-track distance.
        dxt = math.asin(math.sin(dis13 / R) * math.sin(bear13 - bear12)) * R

        # Is p4 beyond the arc?
        dis12 = distance((lat1, lon1), (lat2, lon2)).meters
        dis14 = math.acos(math.cos(dis13 / R) / math.cos(dxt / R)) * R
        if dis14 > dis12:
            return distance((lat2, lon2), (lat3, lon3)).meters
        return math.fabs(dxt)

- user13626872

2

将Java版本添加到wdickerson的回答中：

最初的回答：

public static double pointToLineDistance(double lon1, double lat1, double lon2, double lat2, double lon3, double lat3) {
    lat1 = Math.toRadians(lat1);
    lat2 = Math.toRadians(lat2);
    lat3 = Math.toRadians(lat3);
    lon1 = Math.toRadians(lon1);
    lon2 = Math.toRadians(lon2);
    lon3 = Math.toRadians(lon3);

    // Earth's radius in meters
    double R = 6371000;

    // Prerequisites for the formulas
    double bear12 = bear(lat1, lon1, lat2, lon2);
    double bear13 = bear(lat1, lon1, lat3, lon3);
    double dis13 = dis(lat1, lon1, lat3, lon3);

    // Is relative bearing obtuse?
    if (Math.abs(bear13 - bear12) > (Math.PI / 2))
        return dis13;

    // Find the cross-track distance.
    double dxt = Math.asin(Math.sin(dis13 / R) * Math.sin(bear13 - bear12)) * R;

    // Is p4 beyond the arc?
    double dis12 = dis(lat1, lon1, lat2, lon2);
    double dis14 = Math.acos(Math.cos(dis13 / R) / Math.cos(dxt / R)) * R;
    if (dis14 > dis12)
        return dis(lat2, lon2, lat3, lon3);
    return Math.abs(dxt);
}

private static double dis(double latA, double lonA, double latB, double lonB) {
    double R = 6371000;
    return Math.acos(Math.sin(latA) * Math.sin(latB) + Math.cos(latA) * Math.cos(latB) * Math.cos(lonB - lonA)) * R;
}

private static double bear(double latA, double lonA, double latB, double lonB) {
    // BEAR Finds the bearing from one lat / lon point to another.
    return Math.atan2(Math.sin(lonB - lonA) * Math.cos(latB), Math.cos(latA) * Math.sin(latB) - Math.sin(latA) * Math.cos(latB) * Math.cos(lonB - lonA));
}

- Sga

不确定我是否实现有误，但是无论我尝试了什么，结果都只有应有结果的三分之一。 - Forseth11

1

对于100-1000米范围的球形问题，可以通过等距矩形投影将其转换为笛卡尔空间。然后继续使用学校数学：使用易于找到已实现的“线段上点到线段距离”的函数。该函数使用（有时返回）线段AB上的投影点X的相对前向/后向位置。如果投影点在线段内部，则值在[0,1]之间。如果X在A之前外部，则为负；如果在B之后外部，则大于1。如果相对位置在0和1之间，则取正常距离；如果在外部，则取起点和线段端点A、B中较短的距离。这样的笛卡尔实现示例或非常类似的示例可参考点与线段之间的最短距离。

- AlexWien

当距离增加时，准确度下降是多少？考虑到我们将球面上的距离视为二维平面上的距离。 - ChrisDekker

是的，它基于简单的等经纬投影，将地球视为圆柱体，而不是球体。对于小距离，平面与球体表面相同（只要坐标系被转换为在两个轴上具有相等单位，这对于经度来说并不是这种情况）。有更好的投影方式，如UTM，用于徒步旅行地图等，但它们更加复杂。 - AlexWien

0

/**
 * Calculates the euclidean distance from a point to a line segment.
 *
 * @param v     the point
 * @param a     start of line segment
 * @param b     end of line segment 
 * @return      an array of 2 doubles:
 *              [0] distance from v to the closest point of line segment [a,b],
 *              [1] segment coeficient of the closest point of the segment.
 *              Coeficient values < 0 mean the closest point is a.
 *              Coeficient values > 1 mean the closest point is b.
 *              Coeficient values between 0 and 1 mean how far along the segment the closest point is.
 *
 * @author         Afonso Santos
 */
public static
double[]
distanceToSegment( final R3 v, final R3 a, final R3 b )
{
    double[] results    = new double[2] ;

    final R3     ab_    = b.sub( a ) ;
    final double ab     = ab_.modulus( ) ;

    final R3     av_    = v.sub( a ) ;
    final double av     = av_.modulus( ) ;

    if (ab == 0.0)                       // a and b coincide
    {
        results[0] = av ;                // Distance
        results[1] = 0.0 ;               // Segment coeficient.
    }
    else
    {
        final double avScaProjAb  = av_.dot(ab_) / ab ;
        final double abCoeficient = results[1] = avScaProjAb / ab ;

        if (abCoeficient <= 0.0)                 // Point is before start of the segment ?
            results[0] = av ;                    // Use distance to start of segment.
        else if (abCoeficient >= 1.0)            // Point is past the end of the segment ?
            results[0] = v.sub( b ).modulus() ;    // Use distance to end of segment.
        else                                       // Point is within the segment's start/end perpendicular boundaries.
        {
            if (avScaProjAb >= av)                    // Test to avoid machine float representation epsilon rounding errors that would result in expection on sqrt.
                results[0] = 0.0 ;                    // a, b and v are colinear.
            else
                results[0] = Math.sqrt( av * av - avScaProjAb * avScaProjAb ) ;        // Perpendicular distance from point to segment.
        }
    }

    return results ;
}

上述方法需要笛卡尔三维空间参数，而您要求使用纬度/经度参数。要进行转换，请使用

/**
 * Calculate 3D vector (from center of earth).
 * 
 * @param latDeg    latitude (degrees)
 * @param lonDeg    longitude (degrees)
 * @param eleMtr    elevation (meters)
 * @return          3D cartesian vector (from center of earth).
 * 
 * @author          Afonso Santos
 */
public static
R3
cartesian( final double latDeg, final double lonDeg, final double eleMtr )
{
    return versor( latDeg, lonDeg ).scalar( EARTHMEANRADIUS_MTR + eleMtr ) ;
}

如需了解其余3D/R3代码或如何计算到路径/路线/轨迹的距离，请查看https://sourceforge.net/projects/geokarambola/

- Afonso Santos

0

添加一个Python+Numpy实现（现在您可以将经纬度作为数组传递并同时计算所有距离，无需循环）。

def _angularSeparation(long1, lat1, long2, lat2):
    """All radians
    """
    t1 = np.sin(lat2/2.0 - lat1/2.0)**2
    t2 = np.cos(lat1)*np.cos(lat2)*np.sin(long2/2.0 - long1/2.0)**2
    _sum = t1 + t2

    if np.size(_sum) == 1:
        if _sum < 0.0:
            _sum = 0.0
    else:
        _sum = np.where(_sum < 0.0, 0.0, _sum)

    return 2.0*np.arcsin(np.sqrt(_sum))


def bear(latA, lonA, latB, lonB):
    """All radians
    """
    # BEAR Finds the bearing from one lat / lon point to another.
    result = np.arctan2(np.sin(lonB - lonA) * np.cos(latB),
                        np.cos(latA) * np.sin(latB) - np.sin(latA) * np.cos(latB) * np.cos(lonB - lonA)
                        )

    return result


def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
    """All radians
    points 1 and 2 define an arc segment,
    this finds the distance of point 3 to the arc segment. 
    """

    result = lon1*0
    needed = np.ones(result.size, dtype=bool)

    bear12 = bear(lat1, lon1, lat2, lon2)
    bear13 = bear(lat1, lon1, lat3, lon3)
    dis13 = _angularSeparation(lon1, lat1, lon3, lat3)

    # Is relative bearing obtuse?
    diff = np.abs(bear13 - bear12)
    if np.size(diff) == 1:
        if diff > np.pi:
            diff = 2*np.pi - diff
        if diff > (np.pi / 2):
            return dis13
    else:
        solved = np.where(diff > (np.pi / 2))[0]
        result[solved] = dis13[solved]
        needed[solved] = 0
    
    # Find the cross-track distance.
    dxt = np.arcsin(np.sin(dis13) * np.sin(bear13 - bear12))

    # Is p4 beyond the arc?
    dis12 = _angularSeparation(lon1, lat1, lon2, lat2)
    dis14 = np.arccos(np.cos(dis13) / np.cos(dxt))
    if np.size(dis14) == 1:
        if dis14 > dis12:
            return _angularSeparation(lon2, lat2, lon3, lat3)
    else:
        solved = np.where(dis14 > dis12)[0]
        result[solved] = _angularSeparation(lon2[solved], lat2[solved], lon3[solved], lat3[solved])

    if np.size(lon1) == 1:
        return np.abs(dxt)
    else:
        result[needed] = np.abs(dxt[needed])
        return result

- I.P. Freeley

0

添加一个wdickerson实现的ObjectiveC翻译：

#define DEGREES_RADIANS(angle) ((angle) / 180.0 * M_PI)
#define RADIANS_DEGREES(angle) ((angle) / M_PI * 180)

- (double)crossArcFromCoord:(CLLocationCoordinate2D)fromCoord usingArcFromCoord:(CLLocationCoordinate2D)arcCoord1 toArcCoord:(CLLocationCoordinate2D)arcCoord2 {

        fromCoord.latitude = DEGREES_RADIANS(fromCoord.latitude);
        fromCoord.longitude = DEGREES_RADIANS(fromCoord.longitude);

        arcCoord1.latitude = DEGREES_RADIANS(arcCoord1.latitude);
        arcCoord1.longitude = DEGREES_RADIANS(arcCoord1.longitude);

        arcCoord2.latitude = DEGREES_RADIANS(arcCoord2.latitude);
        arcCoord2.longitude = DEGREES_RADIANS(arcCoord2.longitude);

        double R = 6371000; // Earth's radius in meters

        // Prerequisites for the formulas
        double bear12 = [self bearFromCoord:arcCoord1 toCoord:arcCoord2];
        double bear13 = [self bearFromCoord:arcCoord1 toCoord:fromCoord];

        double dis13 = [self distFromCoord:arcCoord1 toCoord:fromCoord];

        double diff = fabs(bear13 - bear12);

        if (diff > M_PI) {
            diff = 2 * M_PI - diff;
        }

        // Is relative bearing obtuse?
        if (diff > (M_PI/2)) {
            return dis13;
        }

        // Find the cross-track distance
        double dxt = asin(sin(dis13 / R) * sin(bear13 - bear12)) * R;

        // Is p4 beyond the arc?
        double dis12 = [self distFromCoord:arcCoord1 toCoord:arcCoord2];
        double dis14 = acos(cos(dis13 / R) / cos(dxt / R)) * R;

        if (dis14 > dis12) {
            return [self distFromCoord:arcCoord2 toCoord:fromCoord];
        }

        return fabs(dxt);
    }

    - (double)distFromCoord:(CLLocationCoordinate2D)coord1 toCoord:(CLLocationCoordinate2D)coord2 {

        double R = 6371000;

        return acos(sin(coord1.latitude) * sin(coord2.latitude) + cos(coord1.latitude) * cos(coord2.latitude) * cos(coord2.longitude - coord2.longitude)) * R;
    }

    - (double)bearFromCoord:(CLLocationCoordinate2D)fromCoord toCoord:(CLLocationCoordinate2D)toCoord {

        return atan2(sin(toCoord.longitude - fromCoord.longitude) * cos(toCoord.latitude),
                     cos(fromCoord.latitude) * sin(toCoord.latitude) - (sin(fromCoord.latitude) * cos(toCoord.latitude) * cos(toCoord.longitude - fromCoord.longitude)));
     }

- rsarto

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- wdickerson · Accepted Answer

首先，一些术语：
我们的弧线从p1到p2绘制。
我们的第三个点是p3。
相交于大圆的虚拟点是p4。
p1由lat1，lon1定义；p2由lat2，lon2定义；以此类推。
dis12是从p1到p2的距离；等等。
bear12是从p1到p2的方位角；等等。
dxt是横跨距离。
dxa是横跨弧线的距离，也就是我们的目标！

请注意，横跨式公式依赖于相对方位角bear13-bear12

我们有3种情况要处理。

情况1：相对方位角是钝角。因此，dxa=dis13。

情况2.1：相对方位角是锐角，且p4在我们的弧线上。因此，dxa=dxt。

情况2.2：相对方位角是锐角，且p4在我们的弧线之外。因此，dxa=dis23

算法：

步骤1：如果相对方位角是钝角，则dxa=dis13
完成！
步骤2：如果相对方位角是锐角：
2.1：找到dxt。
2.3：找到dis12。
2.4：找到dis14。
2.4：如果dis14>dis12，则dxa=dis23。
完成！
2.5：如果我们到达这里，则dxa=abs(dxt)

MATLAB代码：

function [ dxa ] = crossarc( lat1,lon1,lat2,lon2,lat3,lon3 )
%// CROSSARC Calculates the shortest distance in meters 
%// between an arc (defined by p1 and p2) and a third point, p3.
%// Input lat1,lon1,lat2,lon2,lat3,lon3 in degrees.
    lat1=deg2rad(lat1); lat2=deg2rad(lat2); lat3=deg2rad(lat3);
    lon1=deg2rad(lon1); lon2=deg2rad(lon2); lon3=deg2rad(lon3);

    R=6371000; %// Earth's radius in meters
    %// Prerequisites for the formulas
    bear12 = bear(lat1,lon1,lat2,lon2);
    bear13 = bear(lat1,lon1,lat3,lon3);
    dis13 = dis(lat1,lon1,lat3,lon3);

    diff = abs(bear13-bear12);
    if diff > pi
        diff = 2 * pi - diff;
    end
    %// Is relative bearing obtuse?
    if diff>(pi/2)
        dxa=dis13;
    else
        %// Find the cross-track distance.
        dxt = asin( sin(dis13/R)* sin(bear13 - bear12) ) * R;

        %// Is p4 beyond the arc?
        dis12 = dis(lat1,lon1,lat2,lon2);
        dis14 = acos( cos(dis13/R) / cos(dxt/R) ) * R;
        if dis14>dis12
            dxa=dis(lat2,lon2,lat3,lon3);
        else
            dxa=abs(dxt);
        end   
    end
end

function [ d ] = dis( latA, lonA, latB, lonB )
%DIS Finds the distance between two lat/lon points.
R=6371000;
d = acos( sin(latA)*sin(latB) + cos(latA)*cos(latB)*cos(lonB-lonA) ) * R;
end

function [ b ] = bear( latA,lonA,latB,lonB )
%BEAR Finds the bearing from one lat/lon point to another.
b=atan2( sin(lonB-lonA)*cos(latB) , ...
    cos(latA)*sin(latB) - sin(latA)*cos(latB)*cos(lonB-lonA) );
end

示例输出：演示所有情况。请参见下面的地图。

>> crossarc(-10.1,-55.5,-15.2,-45.1,-10.5,-62.5)
ans =
   7.6709e+05
>> crossarc(40.5,60.5,50.5,80.5,51,69)
ans =
   4.7961e+05
>> crossarc(21.72,35.61,23.65,40.7,25,42)
ans =
   1.9971e+05

地图上的同样输出：

演示了情况1：

演示了情况2.1：

演示了情况2.2：

http://www.movable-type.co.uk/scripts/latlong.html提供公式，
http://www.darrinward.com/lat-long/?id=1788764生成地图图像。