Расстояние от Lat/Lng указывает на сегмент Minor Arc

Мне нужно рассчитать кратчайшее расстояние от точки GPS lat/lng GPS до отрезка линии, описываемого двумя другими точками GPS и A точек BAT/lng.

"Расстояние между дорожками" помогает мне рассчитать кратчайшее расстояние между 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

Следующие изображения, надеюсь, демонстрируют проблему, которую я пытаюсь решить:

На первом изображении расстояние Cross-Track, обозначенное строкой зеленая, является правильным и действительно самым коротким расстоянием до отрезка AB.

Во втором изображении показана проблема с расстоянием между дорожками. В этом случае я хотел бы, чтобы кратчайшее расстояние было простым расстоянием AP, но расстояние между дорожками дает мне расстояние, обозначенное значком красный.

Как мне изменить свой алгоритм, чтобы принять это во внимание, или проверить, находится ли точка X внутри AB. Можно ли это сделать вычислительно? Или это итеративное единственное возможное (дорогостоящее) решение? (возьмите N точек вдоль AB и вычислите минимальное расстояние от P до всех этих точек)

Для целей простоты все линии на изображениях прямые. На самом деле это небольшие дуги на большом круге

Ответ 1

Во-первых, некоторые номенклатуры:
Наша дуга рисуется от p1 до p2.
Наш третий пункт - p3.
Мнимая точка, пересекающая большой круг, равна р4.
p1 определяется lat1, lon1; p2 - lat2, lon2; и т.д.
dis12 - расстояние от p1 до p2; и т.д.
bear12 - это подшипник от p1 до p2; и т.д.
dxt - расстояние между дорожками.
dxa - это поперечное расстояние, наша цель!

Обратите внимание, что формула кросс-трека опирается на относительный подшипник, bear13-bear12

У нас есть 3 случая.

Случай 1: Относительный подшипник тупой. Итак, dxa = dis13.

Случай 2.1: Относительный подшипник острый, AND 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 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);

    %// Is relative bearing obtuse?
    if abs(bear13-bear12)>(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
для генерации изображений карты.

Ответ 2

Для сферических задач 100-1000 м легко перейти к декартово пространство, используя равноугольную проекцию.
Затем он продолжается со школьной математикой:
Используйте функцию "расстояние от сегмента линии", которое легко найти в готовом исполнении. Это fucntion использует (и иногда возвращает) относительную позицию вперед/назад для проецируемой точки X на линии A, B. Значение

в интервале [0,1], если проецируемая точка находится внутри сегмента линии.
отрицательно, если X вне до А,
это > 1, если снаружи после B.
Если относительное положение находится между 0,1, то выполняется нормальное расстояние, если за пределами более короткого расстояния между начальными и конечными точками A, B.

Примером такой/или очень похожей декартовой реализации является кратчайшее расстояние между точкой и сегментом линии

Ответ 3

/**
 * 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 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 ;
}

приведенный выше метод требует декартовых пространственных аргументов, и вы попросили использовать аргументы lat/lon. Чтобы сделать преобразование, используйте

/**
 * 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/