Convert latitude, longitude & altitude to local ENU coordinates in Python

Question

How can I convert geodetic (latitude, longitude, altitude) coordinates into local tangent plane ENU (East, North, Up) coordinates with Python?

pyproj package does not seem to have the right functionality ...

Answer 1

Use pyproj without 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

Ref 2 https://www.nsstc.uah.edu/users/phillip.bitzer/python_doc/pyltg/_modules/pyltg/utilities/latlon.html

Ref 3 https://gist.github.com/sbarratt/a72bede917b482826192bf34f9ff5d0b

Answer 2

You could use the pymap3d package:

Installation

pip install pymap3d

Simple example

Let's take the example values shown on 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)

which produces

(-7134.757195979863, -4556.321513844541, 2852.3904239436915)

which are the east, north and up components, respectively.

The used ellipsoid is WGS84 by default. All available ellipsoid models (which can be fed as an argument to geodetic2enu) can be seen here. Here is how you would calculate the same ENU coordinates using WGS72 reference ellipsoid:

pm.geodetic2enu(lat, lon, h, lat0, lon0, h0, ell=pm.utils.Ellipsoid('wgs72'))
# (-7134.754845247729, -4556.320150825548, 2852.3904257449926)

Answer 3

Pymap3d module, https://scivision.github.io/pymap3d/ provides coordinate transforms and geodesy functions including ENU <--> (long, lat, alt).

Here are some examples.

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)

Suppose we have an ENU coordinate system defined with its origin at (lat0, lon0, h0 = 5.0, 48.0, 10.0). And let a point (point_1) with coordinate ENU: (0,0,0) be a test point, this point_1 will be used to do transformation, both direct and inverse.

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)

Now, using the obtained (lat1, lon1, h1), we compute geodetic to ENU conversion.

e2, n2, u2 = pymap3d.geodetic2enu(lat1, lon1, h1, \
                           lat0, lon0, h0, \
                           ell=ell_clrk66, deg=True)
print(e2, n2, u2)

The output should agree with (e1, n1, u1). Small discrepancies is normal for this kind of computation.

In the above computation, the option ell has WGS84 ellipsoid as the default value.