Convert GPS coordinates into 3-dimensional cartesian coordinates (x,y,z)

Question

In a nutshell:

• What is a valid EPSG code for an earth fixed 3-dimensional cartesian coordinate system?
• If it is 6500: How can I trick the python library pyproj into using this 3D coordinate system instead of calculating 2-dimensional results optimized for Minnesota?
• If I can't use pyproj to convert GPS-coordinates into 3-dimensional cartesian coordinates: How can I still do it in python?

Here is the problem I need to solve:

Our project partner operates microwave links. Each of them has a transmitter on one end and a receiver on the other end of the link. The GPS coordiants of these devices are known and I have to calculate the length of the radio links. (We need these lengths together with other data to compare the signal attenuation with the amount of rain that falls in this section.)

We need the cartesian distances (no geodesic distance)
Microwaves don't follow earth curvature. They spread out like beams of light, i.e. along geometric straight lines. So, I don't need any geodesic distances along the surface of the earth (like you need it for "How long is the path of an airplane that flies from New York to Beijing?") Instead I need cartesian distances. For the distances we handle in our project (median length is 2.8 km, average length is 3.8 km, the shortest link is only 47 m long, the longest 43,8 km) there will be not a big difference, but I still want to do it correctly.

How I want to solve the problem:

In our project we use python as a programming language, and I decided to use pyproj for anything that has to do with coordinates on our planet's surface. And in the last weeks I've learned, that there are thousands of different coordinate systems.

GPS (Global Positioning System) seems to be the same as WGS84 (World Geodetic System 1984), or to be more precise: WGS84 seems to be the coordinate system that GPS uses. And WGS84 is the same as EPSG:4326. (EPSG = European Petroleum Survey Group Geodesy; they created a system of key numbers for different coordinate reference systems)

So, what I need to do, is to convert coordinates of two points given in the EPSG:4326 coordinate system into any earth fixed 3-dimensional cartesian system with well known length units. This should give me 2 sets of 3D coordinates from which I easily can calculate the cartesian distance (and if necessary convert this length into meters).

This is what I've tried so far:

So, when there is a EPSG code for GPS coordinates, there also should be a similar code for 3D-cartesian coordinates, and I even found this code. It is EPSG:6500. So, I wrote this little python program:

#!/usr/bin/env python3
from pyproj import Transformer
# EPSG-codes:
# 4326: https://epsg.io/4326 World Geodetic System 1984, used in GPS
# 6500: https://epsg.io/6500-cs Cartesian 3D coordinate system
trans_GPS_to_XYZ = Transformer.from_crs(4326, 6500)
gps = (48.20841,16.37239) # center of Vienna in Austria
xyz = trans_GPS_to_XYZ.transform(gps[0], gps[1])
print(xyz)

But the resulting tuple does not contain 3 coordinates, but only 2. Later I found out, that EPSG:6500 seems to be the code of two very different coordinate systems:

• A 3-dimensional cartesian coordinate system
• A 2-dimensional coordinate system designed to suit best for Minnesota

And pyproj obviously calculates its results for the Minnesota map which is very useless for our project, because we operate in Austria in Europe.

So, here is my question:
In a python program: How can I easily convert convert coordinates given in WGS84 = EPSG:4326 into 3-dimensional earth fixed cartesian coordinates?

Answer 1

To do 3D transformations, you need 3D CRS:

• EPSG:4326 (2D) -> EPSG:4979 (3D)

<Geographic 3D CRS: EPSG:4979>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World: ...
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

• EPSG:4978 (3D Geocentric CRS)

<Geocentric CRS: EPSG:4978>
Name: WGS 84
Axis Info [cartesian]:
- X[geocentricX]: Geocentric X (metre)
- Y[geocentricY]: Geocentric Y (metre)
- Z[geocentricZ]: Geocentric Z (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Then, you need to provide the Z coordinate on the input for the transform:

from pyproj import Transformer
trans_GPS_to_XYZ = Transformer.from_crs(4979, 4978, always_xy=True)
trans_GPS_to_XYZ.transform(16.37239, 48.20841, 0)

Output:

(4085787.068269556, 1200373.402366255, 4732351.138300765)

For maximal correctness, you have to enter the actual elevation (in meters) of your transmitters as Z coordinate. Using 0 as in this example corresponds to WGS84 sea level. – Hagen von Eitzen (from comment below)