当从另一个点只知道 dx dy 时,获取 pos 的经纬度
Get lat lon for a pos when only dx dy is known from an other point
我遇到的问题是,当我知道一个原点位置和到另一点的偏移量(以米为单位)时,获取某个位置的纬度和经度
我想使用以下函数(或任何其他能给我正确的 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。在功能等方面有什么建议吗?
您可以通过两个步骤获得接近预期值的结果:
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
但是这个计算仍然在概念上是错误的。
如果我们假设地球是一个球体(而不是椭圆体),那么第二步 (Go North/South) 就可以了,因为您在子午线上移动,这是一个大圆。沿着大圆的任何方式都是它的起点和终点之间的最短(表面上)路径,除非它绕球体的距离超过一半。测地线是椭圆体上的最短路径,因此适合。
然而,对于第一步,我们想要向东 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°(直线向西)改变了!
因此,测地线计算可能不是解决问题的正确方法。另请注意,顺序或步骤(无论是先 North/South 然后 East/West,还是先 East/West 然后 North/South)都会显着影响这样计算的结果.
'correct' 的顺序(如果有的话)由 dx
和 dy
.
引用的投影/笛卡尔坐标参考系统决定
感谢试用!
不幸的是,我使用了一些旧算法,所以 "correct answer" 是一个更微不足道的函数。
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)
我遇到的问题是,当我知道一个原点位置和到另一点的偏移量(以米为单位)时,获取某个位置的纬度和经度
我想使用以下函数(或任何其他能给我正确的 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。在功能等方面有什么建议吗?
您可以通过两个步骤获得接近预期值的结果:
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
但是这个计算仍然在概念上是错误的。
如果我们假设地球是一个球体(而不是椭圆体),那么第二步 (Go North/South) 就可以了,因为您在子午线上移动,这是一个大圆。沿着大圆的任何方式都是它的起点和终点之间的最短(表面上)路径,除非它绕球体的距离超过一半。测地线是椭圆体上的最短路径,因此适合。
然而,对于第一步,我们想要向东 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°(直线向西)改变了!
因此,测地线计算可能不是解决问题的正确方法。另请注意,顺序或步骤(无论是先 North/South 然后 East/West,还是先 East/West 然后 North/South)都会显着影响这样计算的结果.
'correct' 的顺序(如果有的话)由 dx
和 dy
.
感谢试用!
不幸的是,我使用了一些旧算法,所以 "correct answer" 是一个更微不足道的函数。
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)