Получить lat lon для позиции, когда из другой точки известно только dx dy

У меня проблема с получением широты и долготы для позиции, когда я знаю одну исходную позицию и смещение к другой точке (в метрах)

Я хочу использовать следующую функцию (или любую другую функцию, которая дает мне правильный lat2 и lon2):

import geographiclib
geographiclib.geodesic.Geodesic.Direct(self, lat1, lon1, azi1, s12)

У меня есть некоторые примеры данных:

lon1 = 11.62113333
lat1 = 55.9862

dx = -51659.25 #meter
dy = -33702.33 #meter

Это результат, которого я надеюсь достичь:

azi1 = -120.95109978727244
s12 = 61691.57175978693         
lat2 = 55.69834 
lon2 = 10.77969

Когда я использую конвертер UTM2latlon, я получаю положение, которое далеко. Я думаю, что система координат, которая используется для расчета dx и dy, - это ESPG: 5596.

Так как расстояния велики и расстояния учитывались для земли, теорема Пифагора неприменима для вычисления s12 и azi1. Любые предложения по функциям и т.д.?


python utm geographic-distance
person axel_ande    schedule 18.08.2015    source источник
comment
Вы имеете в виду EPSG:5596? Возможно, GDAL или pyproj могут выполнить преобразование. См. gis.stackexchange.com/questions/78838/.   -  person    schedule 18.08.2015

Ответы (2)


arrow_upward
0
arrow_downward

Вы можете получить результат, близкий к вашим ожидаемым значениям, за два шага:

from geographiclib.geodesic import Geodesic

lon1 = 11.62113333
lat1 = 55.9862

dx = -51659.25 #meter
dy = -33702.33 #meter

tmp = Geodesic.WGS84.Direct(
    lat1, lon1,
    90, # Go East ...
    dx) #         ... (negatively, so actually West)

destination = Geodesic.WGS84.Direct(
    tmp['lat2'], tmp['lon2'],
    0,  # Go North ...
    dy) #          ... (negatively, so actually South)

lat2, lon2 = destination['lat2'], destination['lon2']
# 55.680721562111955, 10.793499275609594

Но этот расчет по-прежнему концептуально неверен.

Если мы предположим, что Земля является сферой (а не эллипсоидом), второй шаг (Движение на север/юг) будет в порядке, поскольку вы движетесь по меридиану, то есть по большому кругу. Любой путь вдоль большого круга является кратчайшим (на поверхности) путем между его начальной и конечной точками, если только он не проходит более половины пути вокруг сферы. Геодезическая — это кратчайший путь на эллипсоиде, так что это подходит.

Однако для первого шага нам нужно пройти dx метра на восток (или -dx метра на запад). Следуя геодезической, направление на запад изначально не будет делать этого, потому что параллели (за исключением экватора) не являются большими кругами / геодезическими. Таким образом, рассчитанный путь будет отклоняться от параллели, что легко увидеть, взглянув на промежуточное местоположение:

tmp
# {'a12': -0.4645520407433718,
#  'azi1': 90.0,
#  'azi2': 89.31397924368152,
#  'lat1': 55.9862,
#  'lat2': 55.98342231490044,
#  'lon1': 11.62113333,
#  'lon2': 10.793499275609594,
#  's12': -51659.25}

Широта там (tmp['lat2']) не является нашей начальной широтой, и азимут (tmp['azi2']) тоже изменился с 90° (прямо на запад)!

Таким образом, геодезические расчеты, вероятно, не являются правильным способом решения вашей проблемы. Также обратите внимание, что порядок или шаги (идете ли вы сначала на север/юг, а затем на восток/запад или сначала на восток/запад, а затем на север/юг) существенно повлияют на результат при таком расчете.

Какой порядок будет «правильным» (если есть) определяется проекционной/декартовой системой координат, на которую ссылаются dx и dy.

person das-g    schedule 18.08.2015

arrow_upward
0
arrow_downward

Спасибо за попытку!

К сожалению, я работал с каким-то старым алгоритмом, поэтому «правильный ответ» был более тривиальной функцией.

def melat(lat):
    return  degrees( log ( tan ( radians( lat / 2 + 45 ) ) ))

def latit(mlt):
    return ( degrees(atan ( exp ( radians(mlt) ) )) - 45 ) * 2

def XY2LatLong(x, y):
    EQUATOR_MINUTE_LENGTH = 1851.8518519
    GeoPoint = namedtuple('GeoPoint', ('lat', 'lon'))
    Point = namedtuple('Point', ('x', 'y'))
    base = GeoPoint(55.98619572, 11.62113936)
    factor = cos( radians(base.lat) ) * EQUATOR_MINUTE_LENGTH * 60.
    return GeoPoint(
        latit(melat(base.lat) + y / factor),
        base.lon + x / factor)
person axel_ande    schedule 14.09.2015