Question About Obtaining Rotation Matrices

[specarmi]

Hello and thank you for the wonderful tool!

I have a position given in lat/lon/height and ECEF. I'd like to get the rotation matrix from a NED frame centered at this position to ECEF.

What is the most direct way to do this with pygeodesy?

The closest I have gotten is

pygeodesy.ltp.LocalCartesian(lat, lon, height).M

and here I believe M represents a rotation from ECEF to some local level frame, but it's not totally clear to me what that local frame is.

[mrJean1]

That pygeodesy.ltp.LocalCartesian(lat0, lon0, height0) instance and represents a Local Tangent Plane at the (lat0, lon0, height0) location, intended to convert between nearby geodetic (lat, lon, height) locations and local (x, y, z) coordinates on that Local Tangent Plane, back and forth.

Matrix M is the initial rotation matrix is used in methods forward and reverse to perform each conversion. If keyword argument M=False of those methods is set to True, a different matrix, the concatenated rotation matrix is included in the result returned by methods forward and reverse.

All that is straight from Karney's C++ classes LocalCartesian and Geocentric. See methods Forward (2/2) and Reverse (2/2) for some more details about the concatenated rotation matrix.

HTH.

[specarmi]

Thanks for your response, especially pointing me to that additional documentation. It looks like M is the rotation matrix from east-north-up (ENU) to ECEF

I have confirmed this with the following:

import numpy as np
import pygeodesy
import math

# Random values
lat_deg = 40
lon_deg = -30
height = 0.0

lat = math.radians(lat_deg)
lon = math.radians(lon_deg)

s_lat = math.sin(lat)
s_lon = math.sin(lon)
c_lat = math.cos(lat)
c_lon = math.cos(lon)

# From https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU
R_ecef_to_enu = np.eye(3)
R_ecef_to_enu[0, 0] = -s_lon
R_ecef_to_enu[0, 1] = c_lon
R_ecef_to_enu[0, 2] = 0.0
R_ecef_to_enu[1, 0] = -s_lat * c_lon
R_ecef_to_enu[1, 1] = -s_lat * s_lon
R_ecef_to_enu[1, 2] = c_lat
R_ecef_to_enu[2, 0] = c_lat * c_lon
R_ecef_to_enu[2, 1] = c_lat * s_lon
R_ecef_to_enu[2, 2] = s_lat

R_enu_to_ecef_pg = pygeodesy.ltp.LocalCartesian(lat_deg, lon_deg, height).M
R_enu_to_ecef = np.array([R_enu_to_ecef_pg.row(0), R_enu_to_ecef_pg.row(1), R_enu_to_ecef_pg.row(2)])

print(R_ecef_to_enu)
print(R_enu_to_ecef)
print(R_ecef_to_enu - R_enu_to_ecef.T)

which outputs

[[ 0.5         0.8660254   0.        ]
 [-0.5566704   0.3213938   0.76604444]
 [ 0.66341395 -0.38302222  0.64278761]]
[[ 0.5        -0.5566704   0.66341395]
 [ 0.8660254   0.3213938  -0.38302222]
 [ 0.          0.76604444  0.64278761]]
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

[mrJean1]

Thank you for the confirmation.

In an upcoming release, the EcefMatrix documentation will clarify the matrix as ENU to ECEF.

Also, 2 new properties matrix3 and matrixTransposed3 will return the matrix as 3 rows, directly acceptable by numpy, like numpy.array(M.matrix3).

Classes LocalCartesian and Ltp are equivalent, the latter is a sub-class of the former. However, the Ltp documentation is incorrect: Ltp instanced do have the property M to obtain the rotation matrix.

[specarmi]

Great to hear!