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

Input:

lat = 43.21009
lon = -78.120123
h = 124
print(geodetic_to_geocentric(wgs84, lat, lon, h))

Output:

(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 ?

I can't explain the discrepancy with the online conversion tool without seeing the source code for the latter.

1. Review

1. The comment for the ellipsoid parameters should give the units for the semi-major axis. Also, I think "reciprocal flattening" would be a better description of the second parameter, given that you use the name rf later on in the code.




2. The docstring for geodetic_to_geocentric should give the units (degrees) for the latitude and longitude parameters.




3. "Height above ellipsoid" would be a clearer description for the h parameter.




4. The parameters to geodetic_to_geocentric could be spelled out in full, which would help understand their meaning. Parameter names are documentation too.




5. Comments would help understand the meaning of values like N.




6. The expressions (N + h) * cos(lat), sin(lat), and (1 - 1 / rf) ** 2 are evaluated twice. These could be stored in local variables to avoid duplicated work.

2. Revised code

from math import cos, radians, sin, sqrt

# Ellipsoid parameters: semi major axis in metres, reciprocal flattening.
GRS80 = 6378137, 298.257222100882711
WGS84 = 6378137, 298.257223563

def geodetic_to_geocentric(ellipsoid, latitude, longitude, height):
 """Return geocentric (Cartesian) Coordinates x, y, z corresponding to
 the geodetic coordinates given by latitude and longitude (in
 degrees) and height above ellipsoid. The ellipsoid must be
 specified by a pair (semi-major axis, reciprocal flattening).

 """
 φ = radians(latitude)
 λ = radians(longitude)
 sin_φ = sin(φ)
 a, rf = ellipsoid # semi-major axis, reciprocal flattening
 e2 = 1 - (1 - 1 / rf) ** 2 # eccentricity squared
 n = a / sqrt(1 - e2 * sin_φ ** 2) # prime vertical radius
 r = (n + height) * cos(φ) # perpendicular distance from z axis
 x = r * cos(λ)
 y = r * sin(λ)
 z = (n * (1 - e2) + height) * sin_φ
 return x, y, z

Precision in numeric calculus is usually handled by the decimal module. Compare:

>>> 1.001 + 0.0001
1.0010999999999999
>>> Decimal('1.001') + Decimal('0.0001')
Decimal('1.0011')

However, trigonometric functions would convert a decimal.Decimal back to a float, thus defeating the purpose. However, the decimal module provide some recipe to help fill the void. Using the recipe for pi() you can reimplement math.radians easily; cos() and sin() are covered by the recipes; and math.sqrt is available directly on Decimal objects. So you could write:

def geodetic_to_geocentric(ellps, lat, lon, h):
 a, rf = ellps
 lat_rad = radians(lat)
 lon_rad = radians(lon)
 N = a / (1 - (1 - (1 - 1 / rf) ** 2) * (sin(lat_rad)) ** 2).sqrt()
 X = (N + h) * cos(lat_rad) * cos(lon_rad)
 Y = (N + h) * cos(lat_rad) * sin(lon_rad)
 Z = ((1 - 1 / rf) ** 2 * N + h) * sin(lat_rad)

 return X, Y, Z

Which would yield:

>>> lat = Decimal('43.21009')
>>> lon = Decimal('-78.120123')
>>> h = Decimal('124')
>>> print(geodetic_to_geocentric(wgs84, lat, lon, h))
(Decimal('958506.0102730405668418845812'), Decimal('-4556367.372558064424670955239'), Decimal('4344627.161663229280843240044'))

Which is closer to your result than the one you expected. So I wouldn't worry in this case.

Note however that computation using Decimals are usually slower

>>> timeit.timeit("geodetic_to_geocentric((Decimal('6378137'), Decimal('298.257223563')), Decimal('43.21009'), Decimal('-78.120123'), Decimal('124'))", setup='from __main__ import geodetic_to_geocentric_decimal as geodetic_to_geocentric; from decimal import Decimal')
86.26855880199582
>>> timeit.timeit("geodetic_to_geocentric((6378137, 298.257223563), 43.21009, -78.120123, 124)", setup='from __main__ import geodetic_to_geocentric')
1.40149550899514

Lastly, since you cache values of angles in radian, you could as well store some other values that you compute twice in variable. This includes:

• sin(lat_rand)


• cos(lat_rand)


• (1 - 1 / rf) ** 2

And talking about variables, you could make an effort to improve naming a bit. Some variable names may make sense to you as they may be usual in your domain; but rf, ellps or N doesn't convey much meaning to me, for instance.