I'd like to retain as much float precision as possible without having to resort to specialized libraries so I need this code to work with only standard modules.
I'm doing this in QGIS where specialized libraries aren't available by default.
numpy would be OK to use.
Sources for the computations can be found on Wikipedia.
Here is my code:
import math ## Ellipsoid Parameters as tuples (semi major axis, inverse flattening) grs80 = (6378137, 298.257222100882711) wgs84 = (6378137, 298.257223563) def geodetic_to_geocentric(ellps, lat, lon, h): """ Compute the Geocentric (Cartesian) Coordinates X, Y, Z given the Geodetic Coordinates lat, lon + Ellipsoid Height h """ a, rf = ellps lat_rad = math.radians(lat) lon_rad = math.radians(lon) N = a / math.sqrt(1 - (1 - (1 - 1 / rf) ** 2) * (math.sin(lat_rad)) ** 2) X = (N + h) * math.cos(lat_rad) * math.cos(lon_rad) Y = (N + h) * math.cos(lat_rad) * math.sin(lon_rad) Z = ((1 - 1 / rf) ** 2 * N + h) * math.sin(lat_rad) return X, Y, Z
lat = 43.21009 lon = -78.120123 h = 124 print(geodetic_to_geocentric(wgs84, lat, lon, h))
(958506.0102730404, -4556367.372558064, 4344627.16166323)
This online converter gives the following output:
X: 958506.01027304 Y: -4556367.37255806 Z: 4344627.16160147
So, quite close for X and Y but slightly different for Z. Is this a rounding error ? If so, on which side ? Is there a better strategy to do the computations that may retain better precision ?