4
\$\begingroup\$

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 ?

\$\endgroup\$
  • \$\begingroup\$ What are the units of the output? How far off is it, in meters? \$\endgroup\$ – 200_success Jun 6 '18 at 6:41
  • \$\begingroup\$ Units are m yes. So it's actually very close but I'm just wondering why the Z computation is off after the 4th decimal when X and Y are pretty consistent. \$\endgroup\$ – YeO Jun 6 '18 at 6:43
  • 2
    \$\begingroup\$ So, the heights are in agreement to within a tenth of a millimeter, but that's not good enough? What application is this code for? \$\endgroup\$ – 200_success Jun 6 '18 at 6:48
  • 1
    \$\begingroup\$ Because I really want to ? :-) More seriously, I want to do this within QGIS where sympy, bigfloat and such aren't available by default. numpy is actually acceptable as it is available. There won't be a huge amount of points, less than 100 and usually actually around 10~20 max. \$\endgroup\$ – YeO Jun 6 '18 at 8:47
  • 1
    \$\begingroup\$ Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. \$\endgroup\$ – Mathias Ettinger Jun 6 '18 at 9:22
7
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ Thanks. Technically, the value e2 is not the eccentricity squared, that would be 1 - (1 - 1/rf) ** 2. Also, n is the prime vertical radius of curvature, not the distance from centre (φ is the Geodetic Latitude, not the Geocentric one). \$\endgroup\$ – YeO Jun 6 '18 at 10:59
  • \$\begingroup\$ @YeO: Thanks for the corrections, now fixed. \$\endgroup\$ – Gareth Rees Jun 6 '18 at 11:06
  • \$\begingroup\$ There's a reasonable numerical analytic case to be made for evaluating e2 as (2 - 1 / rf) / rf: if using a different ellipsoid which has rf very close to 1 or very large, you eliminate potential loss of significance due to subtracting a number very close to 1 from 1. \$\endgroup\$ – Peter Taylor Jun 6 '18 at 14:48
  • 1
    \$\begingroup\$ Upvoted just for using greek letters in variable names. This is a severely under-used feature of Python that could help dramatically simplify the correlation between code and traditional math or academic papers. \$\endgroup\$ – scnerd Jun 6 '18 at 18:45
  • \$\begingroup\$ @PeterTaylor Thanks, that is actually the kind of advice I was looking for! In this case, rf will usually be ~300 as I'm interested in actual coordinates on earth. \$\endgroup\$ – YeO Jun 6 '18 at 19:19
5
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ thanks for the advice. I always thought caching a value did not retain as much precision as doing the computation each time. With regards to variable naming, fair enough remark, but I think a for semi major axis and rf for reciprocal flattening are pretty standard notations (they are the proj4 parameters for instance) for ellipsoid parameters. \$\endgroup\$ – YeO Jun 6 '18 at 9:06
  • \$\begingroup\$ @YeO If they are standard notation for your domain, then it might be fine, it's just that I know very few of it and these names don't "speak" to me. \$\endgroup\$ – Mathias Ettinger Jun 6 '18 at 9:07
  • 2
    \$\begingroup\$ @YeO: In Python it makes no difference whether you assign an intermediate value to a variable or not; the value is the same whether or not you give it a name. \$\endgroup\$ – Gareth Rees Jun 6 '18 at 10:41
  • \$\begingroup\$ @GarethRees thanks for this clarification! It's very good to know. \$\endgroup\$ – YeO Jun 6 '18 at 11:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.