在Python中将纬度、经度和高度转换为本地ENU坐标
如何使用Python将大地坐标(纬度、经度、高度)转换为局部切面ENU (东、北、上)坐标？
pyproj包似乎没有正确的功能...
回答 3
Stack Overflow用户
发布于 2020-06-23 03:14:32
您可以使用pymap3d包：
安装
pip install pymap3d简单的例子
让我们以this page上显示的示例值为例。
import pymap3d as pm
# The local coordinate origin (Zermatt, Switzerland)
lat0 = 46.017 # deg
lon0 = 7.750 # deg
h0 = 1673 # meters
# The point of interest
lat = 45.976 # deg
lon = 7.658 # deg
h = 4531 # meters
pm.geodetic2enu(lat, lon, h, lat0, lon0, h0)它会产生
(-7134.757195979863, -4556.321513844541, 2852.3904239436915)它们分别是东、北和上分量。
默认情况下，使用的椭球体为WGS84。所有可用的椭球模型(可以作为geodetic2enu的参数)都可以看到here。下面是如何使用WGS72参考椭球计算相同的ENU坐标：
pm.geodetic2enu(lat, lon, h, lat0, lon0, h0, ell=pm.utils.Ellipsoid('wgs72'))
# (-7134.754845247729, -4556.320150825548, 2852.3904257449926)Stack Overflow用户
发布于 2020-11-28 18:16:10
使用不带pymap3d的pyproj
import numpy as np
import pyproj
import scipy.spatial.transform
def geodetic2enu(lat, lon, alt, lat_org, lon_org, alt_org):
transformer = pyproj.Transformer.from_crs(
{"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
{"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
)
x, y, z = transformer.transform( lon,lat, alt,radians=False)
x_org, y_org, z_org = transformer.transform( lon_org,lat_org, alt_org,radians=False)
vec=np.array([[ x-x_org, y-y_org, z-z_org]]).T
rot1 = scipy.spatial.transform.Rotation.from_euler('x', -(90-lat_org), degrees=True).as_matrix()#angle*-1 : left handed *-1
rot3 = scipy.spatial.transform.Rotation.from_euler('z', -(90+lon_org), degrees=True).as_matrix()#angle*-1 : left handed *-1
rotMatrix = rot1.dot(rot3)
enu = rotMatrix.dot(vec).T.ravel()
return enu.T
def enu2geodetic(x,y,z, lat_org, lon_org, alt_org):
transformer1 = pyproj.Transformer.from_crs(
{"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
{"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
)
transformer2 = pyproj.Transformer.from_crs(
{"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
{"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
)
x_org, y_org, z_org = transformer1.transform( lon_org,lat_org, alt_org,radians=False)
ecef_org=np.array([[x_org,y_org,z_org]]).T
rot1 = scipy.spatial.transform.Rotation.from_euler('x', -(90-lat_org), degrees=True).as_matrix()#angle*-1 : left handed *-1
rot3 = scipy.spatial.transform.Rotation.from_euler('z', -(90+lon_org), degrees=True).as_matrix()#angle*-1 : left handed *-1
rotMatrix = rot1.dot(rot3)
ecefDelta = rotMatrix.T.dot( np.array([[x,y,z]]).T )
ecef = ecefDelta+ecef_org
lon, lat, alt = transformer2.transform( ecef[0,0],ecef[1,0],ecef[2,0],radians=False)
return [lat,lon,alt]
if __name__ == '__main__':
# The local coordinate origin (Zermatt, Switzerland)
lat_org = 46.017 # deg
lon_org = 7.750 # deg
alt_org = 1673 # meters
# The point of interest
lat = 45.976 # deg
lon = 7.658 # deg
alt = 4531 # meters
res1 = geodetic2enu(lat, lon, alt, lat_org, lon_org, alt_org)
print (res1)
#[-7134.75719598 -4556.32151385 2852.39042395]
x=res1[0]
y=res1[1]
z=res1[2]
res2 = enu2geodetic(x,y,z, lat_org, lon_org, alt_org)
print (res2)
#[45.97600000000164, 7.658000000000001, 4531.0000001890585]参考1个https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates
参考文献2 https://www.nsstc.uah.edu/users/phillip.bitzer/python_doc/pyltg/_modules/pyltg/utilities/latlon.html
参考文献3 https://gist.github.com/sbarratt/a72bede917b482826192bf34f9ff5d0b
Stack Overflow用户
发布于 2019-06-01 18:05:28
Pymap3d模块，https://scivision.github.io/pymap3d/提供坐标转换和大地测量功能，包括ENU <--> (long，lat，alt)。
这里有一些例子。
import pymap3d
# create an ellipsoid object
ell_clrk66 = pymap3d.Ellipsoid('clrk66')
# print ellipsoid's properties
ell_clrk66.a, ell_clrk66.b, ell_clrk66.f
# output
(6378206.4, 6356583.8, 0.0033900753039287634)假设我们有一个原点为(lat0，lon0，h0 = 5.0，48.0，10.0)的ENU坐标系。让一个坐标为ENU：(0,0,0)的点(point_1)作为一个测试点，这个point_1将被用来做正反变换。
lat0, lon0, h0 = 5.0, 48.0, 10.0 # origin of ENU, (h is height above ellipsoid)
e1, n1, u1 = 0.0, 0.0, 0.0 # ENU coordinates of test point, `point_1`
# From ENU to geodetic computation
lat1, lon1, h1 = pymap3d.enu2geodetic(e1, n1, u1, \
lat0, lon0, h0, \
ell=ell_clrk66, deg=True) # use clark66 ellisoid
print(lat1, lon1, h1)
# display: (5.000000000000001, 48.0, 10.000000000097717)现在，使用获得的(lat1，lon1，h1)，我们计算大地坐标到ENU的转换。
e2, n2, u2 = pymap3d.geodetic2enu(lat1, lon1, h1, \
lat0, lon0, h0, \
ell=ell_clrk66, deg=True)
print(e2, n2, u2)输出应与(e1，n1，u1)一致。对于这种计算来说，小的差异是正常的。
在上面的计算中，选项ell使用WGS84椭球体作为默认值。
https://stackoverflow.com/questions/53408780
复制相似问题
腾讯云开发者
