Расчет координат по азимуту и ​​расстоянию - PullRequest
4 голосов
/ 18 мая 2009

У меня проблемы с реализацией описанной здесь функции здесь .

Это моя реализация Java:

private static double[] pointRadialDistance(double lat1, double lon1, 
        double radianBearing, double radialDistance) {
     double lat = Math.asin(Math.sin(lat1)*Math.cos(radialDistance)+Math.cos(lat1)
             *Math.sin(radialDistance)*Math.cos(radianBearing));
     double lon;
     if(Math.cos(lat) == 0) {  // Endpoint a pole
        lon=lon1;      
     }
     else {
        lon = ((lon1-Math.asin(Math.sin(radianBearing)*Math.sin(radialDistance)/Math.cos(lat))
                +Math.PI) % (2*Math.PI)) - Math.PI;
     }
    return (new double[]{lat, lon});
}

Я преобразую градусный градус в радианы и преобразую расстояние (км) в радиан, прежде чем вызывать функцию - так что это не проблема.

Однако при вводе таких координат, как: широта = 49,25705; lon = -123,140259; с азимутом 225 (юго-запад) и расстоянием 1 км

Я получил это возвращено: широта: -1.0085434360125864 длинный: -3,7595299668539504

Это явно не правильно, кто-нибудь может увидеть, что я делаю неправильно?

Спасибо

Ответы [ 6 ]

14 голосов
/ 19 мая 2009

Кажется, что это проблемы в вашем коде:

  1. Вам нужно преобразовать lat1 и lon1 в радианы перед вызовом вашей функции.
  2. Возможно, вы неправильно масштабируете radialDistance.
  3. Проверка числа с плавающей точкой на равенство опасна. Два числа, которые равны после точной арифметики, могут не быть точно равными после арифметики с плавающей точкой. Таким образом, abs(x-y) < threshold безопаснее, чем x == y, для проверки двух чисел с плавающей точкой x и y на равенство.
  4. Я думаю, вы хотите преобразовать lat и lon из радианов в градусы.

Вот моя реализация вашего кода на Python:

#!/usr/bin/env python

from math import asin,cos,pi,sin

rEarth = 6371.01 # Earth's average radius in km
epsilon = 0.000001 # threshold for floating-point equality


def deg2rad(angle):
    return angle*pi/180


def rad2deg(angle):
    return angle*180/pi


def pointRadialDistance(lat1, lon1, bearing, distance):
    """
    Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
    (lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in km]
    """
    rlat1 = deg2rad(lat1)
    rlon1 = deg2rad(lon1)
    rbearing = deg2rad(bearing)
    rdistance = distance / rEarth # normalize linear distance to radian angle

    rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )

    if cos(rlat) == 0 or abs(cos(rlat)) < epsilon: # Endpoint a pole
        rlon=rlon1
    else:
        rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi

    lat = rad2deg(rlat)
    lon = rad2deg(rlon)
    return (lat, lon)


def main():
    print "lat1 \t lon1 \t\t bear \t dist \t\t lat2 \t\t lon2"
    testcases = []
    testcases.append((0,0,0,1))
    testcases.append((0,0,90,1))
    testcases.append((0,0,0,100))
    testcases.append((0,0,90,100))
    testcases.append((49.25705,-123.140259,225,1))
    testcases.append((49.25705,-123.140259,225,100))
    testcases.append((49.25705,-123.140259,225,1000))
    for lat1, lon1, bear, dist in testcases:
        (lat,lon) = pointRadialDistance(lat1,lon1,bear,dist)
        print "%6.2f \t %6.2f \t %4.1f \t %6.1f \t %6.2f \t %6.2f" % (lat1,lon1,bear,dist,lat,lon)


if __name__ == "__main__":
    main()

Вот вывод:

lat1     lon1        bear    dist        lat2        lon2
  0.00     0.00       0.0       1.0        0.01        0.00
  0.00     0.00      90.0       1.0        0.00       -0.01
  0.00     0.00       0.0     100.0        0.90        0.00
  0.00     0.00      90.0     100.0        0.00       -0.90
 49.26   -123.14     225.0      1.0       49.25      -123.13
 49.26   -123.14     225.0    100.0       48.62      -122.18
 49.26   -123.14     225.0   1000.0       42.55      -114.51
3 голосов
/ 18 октября 2011

Я думаю, что есть проблема в алгоритме, представленном в сообщении 5.

Это работает, но только для широты, для долготы есть проблема из-за знака.

Данные говорят сами за себя:

49,26 -123,14 225,0 1,0 49,25 -123,13

Если вы начнете с -123.14 ° и пойдете на ЗАПАД, у вас должно быть что-то ДАЛЕКО на ЗАПАДЕ. Здесь мы возвращаемся на ВОСТОК (-123,13)!

Формула должна включать где-то:

СтепеньBearing = ((360-градусныйBearing)% 360)

до преобразования в радианы.

3 голосов
/ 18 мая 2009

По сути, кажется, что ваша проблема в том, что вы передаете широту, долготу и направление в градусах, а не в радианах. Попробуйте убедиться, что вы всегда передаете радианы для своей функции, и посмотрите, что вы получите обратно.

PS: смотри похожие вопросы, обсуждаемые здесь и здесь .

0 голосов
/ 07 мая 2011

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

Быстрый трал вашего любимого поискового движка для "Vincenty Formula", надеюсь, окажется полезным.

0 голосов
/ 01 октября 2010

Спасибо за ваш код Python Я попытался настроить его в моем случае использования где я пытаюсь найти точку между двумя другими на заданном расстоянии от первой точки, так что это довольно похоже на ваш код Похоже, что моя азимут рассчитывается динамически

начальная точка (lat1) lon1 / lat1 = 55.625541, -21.142463

конечная точка (lat2) lon2 / lat2 = 55.625792, -22.142248

мой результат должен быть точкой между этими двумя в lon3 / lat3 к сожалению, я получаю lon3 / lat3 = 0.0267695450609,0.0223553243666

Я думал, что это может быть разницей в лат, но нет когда я добавляю или добавляю это не хорошо

любой совет был бы действительно хорош Спасибо

вот моя реализация

расстояние = 0,001 эпсилон = 0,000001

Расчет подшипника динамически

y = math.sin(distance) * math.cos(lat2);
x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(distance);
bearing = math.atan2(y, x)

динамический расчет lat3 lon3

rlat1 = (lat1 * 180) / math.pi
rlon1 = (lon1 * 180) / math.pi
rbearing = (bearing * 180) / math.pi
rdistance = distance / R # normalize linear distance to radian angle

rlat = math.asin( math.sin(rlat1) * math.cos(rdistance) + math.cos(rlat1) * math.sin(rdistance) * math.cos(rbearing) )
if math.cos(rlat) == 0 or abs(math.cos(rlat)) < epsilon: # Endpoint a pole
      rlon=rlon1
else:
    rlon = ( (rlon1 + math.asin( math.sin(rbearing)* math.sin(rdistance) / math.cos(rlat) ) + math.pi ) % (2*math.pi) ) - math.pi

lat3 = (rlat * math.pi)/ 180
lon3 = (rlon * math.pi)/ 180
0 голосов
/ 27 августа 2009

Когда я реализовал это, мои полученные широты были правильными, но долготы были неправильными. Например, отправная точка: 36,9460678N 9,434807E, подшипник 45,03334, расстояние 15,0083313 км. Результат был 37.0412865N 9.315302E Это дальше на запад, чем моя отправная точка, а не дальше на восток. На самом деле это как если бы подшипник был 315.03334 градусов.

Больше поиска в сети привело меня к: http://www.movable -type.co.uk / scripts / latlong.html Код долготы показан ниже (в C # со всем в радианах)

        if ((Math.Cos(rLat2) == 0) || (Math.Abs(Math.Cos(rLat2)) < EPSILON))
        {
            rLon2 = rLon1;
        }
        else
        {
            rLon2 = rLon1 + Math.Atan2(Math.Sin(rBearing) * Math.Sin(rDistance) * Math.Cos(rLat1), Math.Cos(rDistance) - Math.Sin(rLat1) * Math.Sin(rLat2));
        }

Мне кажется, это нормально работает. Надеюсь, это полезно.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...