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

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 2

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