Конвертировать лат длинные координаты в локальные координаты ENU в Python

(Как) можно преобразовать координаты длинной широты в локальные координаты ENU (восток, север, вверх) с помощью Python? Я, кажется, не нашел нужную функциональность в пакете pyproj...

Источник

0 ответов

Вы можете использовать pymap3d пакет:

Установка

pip install pymap3d

Простой пример

Возьмем примеры значений, показанные на этой странице.

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) можно увидеть здесь. Вот как можно вычислить те же координаты ENU, используя опорный эллипсоид WGS72:

pm.geodetic2enu(lat, lon, h, lat0, lon0, h0, ell=pm.utils.Ellipsoid('wgs72'))
# (-7134.754845247729, -4556.320150825548, 2852.3904257449926)
Источник

Используйте pyproj без pymap3d

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]

Ref. 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

Источник

Модуль Pymap3d, https://scivision.github.io/pymap3d/, предоставляет преобразования координат и функции геодезии, включая ENU <--> (long, lat, h).

Вот несколько примеров.

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)

Предположим, у нас есть система координат ENU, определенная с ее началом в (lat0, lon0, h0 = 5.0, 48.0, 10.0). И пусть точка (point_1) с координатой ENU: (0,0,0) будет контрольной точкой, это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 эллипсоид в качестве значения по умолчанию.

Источник
Другие вопросы по тегам