0

当我知道一个原点位置和到另一点的偏移量(以米为单位)时,我遇到的问题是获取一个位置的纬度和经度

我想使用以下函数(或任何其他给我正确 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。关于功能等有什么建议吗?

4

2 回答 2

0

您可以通过两个步骤获得接近预期值的结果:

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

于 2015-08-18T22:26:59.237 回答
0

感谢您的尝试!

不幸的是,我使用了一些旧算法,所以“正确答案”是一个更微不足道的函数。

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) 
于 2015-09-14T06:16:07.317 回答