Convert data of the format [longitude, latitude, altitude] into the format [x, y, z]

Question

I have a dataset from flights. The airplane communicates with a ground station and through GPS I have the distance between the current airplane location and the ground station location, alongside with the respective longitude, latitude and altitude. What I need to do now is convert the long, lat, alt position data into a x,y,z coordinate system.

Since the ground station is at a static position I want to have it in the origin (0,0,0) of the coordinate system. I'm in search of a python module that can handle this transformation.

As an example I provide the actual data from the ground station alongside 1 position of the airplane which needs to be converted:

Ground Station [should be (0, 0, 0)] Latitude: 41.947.694
Longitude: 3.209.083 Altitude(m): 379.41

Airplane Position: Lat.: 419.763.015.750 Long.: 34.890.851.772 Alt.[m]: 971.32

Moving the ground station to the origin I would have to subtract its position from the airplane position if I understand correctly which would mean:

Ground Station: Latitude: 0 Longitude: 0 Altitude(m): 0

Airplane Position: Lat.: 419.721.068.056
Long.: 34.887.642.689 Alt.[m]: 591.91

Since I'm new to working with such data and tried some simple solutions and did not find anything that was converting this data I hope you can help me find a solution.

I need the converted data to simulate ellipsoids in 3D to find possible intersections. This could also be done with a different approach than python. So matlab would be fine too eve if I have never worked with that.

Answer 1

Transforming latitude, longitude and altitude GPS coordinates, which are geodetic coordinates (the vertical direction is the normal to the ellipsoid and not the direction to the center of the Earth) based on WGS84 reference ellipsoid, to "xyz" coordinates, geocentric, in the so-called Earth-Centered Earth-Fixed frame is a function implemented in several libraries, in several languages. Since Matlab was mentioned, it is implemented as geodetic2ecef.

It is also described in several texts and web sites, even on wikipedia, so we can safely implement it from scratch.

If you transform both the station and the plane coordinates to this ECEF x,y,z coordinate, you may compute their difference, and you already have an XYZ coordinate sistem that is (0, 0, 0) at the station.

Note, that, indeed, as mentioned in the comments, you cannot perform the differences of the latitudes or longitudes themselves with any meaningful interpretation. To find the distance between two points we can only subtract sizes, as vectors or cartesian coordinates; latitudes and longitudes don't depend linearly on those.

Now, the xyz differences computed as described before are not very significant in themselves, because their reference directions are global - the z-axis is oriented parallel to the axis of rotation of the Earth and the xz plane is parallel to the prime meridian.

For that, I propose to "localize" these differences by adding a (double) rotation of those global directions to bring them to point to directions relevant for the reference station: the new z-axis to be the geodetic vertical (towards Zenith), the new x-axis towards North from the station and the y-axis towards East.

Again, these calculations are well documented in the literature, and online and it is safe to implement them from scratch here.

Here's an example, a possible implementation in python, not very efficient since it repeats the computation of the trigonometric functions for the base, but my intention was to make it as clear as possible rather than efficient.

import math

a = 6378137  # semi major axis, WGS84
# f = 1 / 298.257223563 # WGS84
e2 = 6.69437999014e-3  # eccentricity, WGS84 e2 = 1 - (1-f)*(1-f)


def gps_to_geocentric_ecef(lat, lon, alt):
    # converts geodetic latitude, longitude, altitude based on WGS84 (GPS coordinates)
    # to geocentric x, y, z in meters (ECEF), sqrt(x^2+y^2+z^2) = geocentric radius
    # |z| distance to equatorial plane, |x| distance to the plane of the prime meridian
    # https://en.wikipedia.org/wiki/Geodetic_coordinates#Conversion
    # https://github.com/proj4js/proj4js/blob/master/lib/datumUtils.js#L68C5-L70C44
    # verify: https://tool-online.com/en/coordinate-converter.php (right select WORLD and XYZ (geocentric))
    # or: https://www.oc.nps.edu/oc2902w/coord/llhxyz.htm - paste lat, lon, height, then press "LLH to ECEF"

    lat = math.radians(lat)
    lon = math.radians(lon)
    sin_lat = math.sin(lat)
    cos_lat = math.cos(lat)
    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)

    rn = a / (math.sqrt(1 - e2 * sin_lat * sin_lat))
    x = (rn + alt) * cos_lat * cos_lon
    y = (rn + alt) * cos_lat * sin_lon
    z = ((rn * (1 - e2)) + alt) * sin_lat

    return [x, y, z]


def rotate_ecef_to_local_xyz(dx, dy, dz, lat, lon):
    # rotates the coordinate differences dx, dy, dz fromm the global frame (ECEF) to
    # a local frame that has the geodetic normal to (lat, lon) as z-axis, while
    # the new x-axis is towards North and the y-axis towards East

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

    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)
    cos_lat = math.cos(lat)
    sin_lat = math.sin(lat)
    dx, dy, dz = (- dx * cos_lon * sin_lat - dy * sin_lon * sin_lat + dz * cos_lat,
                  - dx * sin_lon + dy * cos_lon,
                  dx * cos_lon * cos_lat + dy * sin_lon * cos_lat + dz * sin_lat)

    return [dx, dy, dz]


def main():
    bgr_lat, bgr_lon, bgr_alt = [41.947694, 3.209083, 379.41]
    airplane_gps = [
        [41.9763015750,  3.4890851772, 971.32]
    ]

    # ecef coordinates for the station:
    bgr_xyz = gps_to_geocentric_ecef(bgr_lat, bgr_lon, bgr_alt)

    # ecef coordinates for the (first) airplane:
    plane_xyz = gps_to_geocentric_ecef(*airplane_gps[0])

    # xyz position of the airplane relative to the station, ecef orientation:
    diff_xyz = [plane_xyz[i] - bgr_xyz[i] for i in range(3)]

    # xyz position of the airplane relative to the station, station-local orientation
    local_xyz = rotate_ecef_to_local_xyz(*diff_xyz, bgr_lat, bgr_lon)

    print(local_xyz)  # final result, in meters (x-towards North, y-towards East, z-towards Zenith)


if __name__ == '__main__':
    main()

Since the altitudes are given in meters, I used the WGS84 semi major axis also in meters, and the results are in meters: this plane is 3215m to the North, 23210m to the East, and 549m high (measured from the plane tangent to Earth that passes through the reference station, BGR)

Although I thoroughly verified this code in several ways, I still hope nobody will fly a plane based on it :)

Answer 2

You could use pyproj, and define

import pyproj
import math

P = pyproj.Proj(proj='utm', zone=32, ellps='WGS84', preserve_units=True)
G = pyproj.Geod(ellps='WGS84')


base_x_offset = 19930.039590830565
base_y_offset = 4660215.067408236

def lonlat_to_xy(Lon, Lat):
    p = P(Lon, Lat)    
    return (p[0] - base_x_offset, p[1] - base_y_offset)    

def xy_to_lonlat(x,y):
    return P(x + base_x_offset, y + base_y_offset, inverse=True)    

def distance(Lat1, Lon1, Lat2, Lon2):
    return G.inv(Lon1, Lat1, Lon2, Lat2)[2]

And check with:

base_lon = 3.209083
base_lat = 41.947694

plane_lon = 3.4890851772
plane_lat = 41.9763015750

p_base = lonlat_to_xy(base_lon, base_lat)
print(p_base)

p_plane = lonlat_to_xy(plane_lon, plane_lat)
print(p_plane)

and

xy_to_lonlat( lonlat_to_xy(plane_lon, plane_lat)[0] , lonlat_to_xy(plane_lon, plane_lat)[1] )

xy_to_lonlat( lonlat_to_xy(base_lon, base_lat)[0] , lonlat_to_xy(base_lon, base_lat)[1] )

This doesnt handle the z coord. With all the things that of course can go wrong with transforming a sphere to a flat map.