Coordinate conversions#

Geodetic to geocentric spherical#

The boule.Ellipsoid class implements coordinate conversions between geocentric geodetic coordinates and geocentric spherical coordinates. Both are common in geophysical applications when dealing with spherical harmonics or spherical modeling of topography.

The example below will show you how to convert geodetic latitude and height into geocentric spherical latitude and radius.

import boule as bl
import numpy as np

longitude = 40
latitude = np.linspace(-90, 90, 45)
height = 481_000  # ICESat-2 orbit height in meters

longitude_sph, latitude_sph, radius = bl.WGS84.geodetic_to_spherical(
    longitude, latitude, height,
)

print("Geodetic longitude:", longitude)
print("Spherical longitude:", longitude_sph)
print("Geodetic latitude:", latitude)
print("Spherical latitude:", latitude_sph)
print("Height (m):", height)
print("Radius (m):", radius)

Geodetic longitude: 40
Spherical longitude: 40
Geodetic latitude: [-90.         -85.90909091 -81.81818182 -77.72727273 -73.63636364
 -69.54545455 -65.45454545 -61.36363636 -57.27272727 -53.18181818
 -49.09090909 -45.         -40.90909091 -36.81818182 -32.72727273
 -28.63636364 -24.54545455 -20.45454545 -16.36363636 -12.27272727
  -8.18181818  -4.09090909   0.           4.09090909   8.18181818
  12.27272727  16.36363636  20.45454545  24.54545455  28.63636364
  32.72727273  36.81818182  40.90909091  45.          49.09090909
  53.18181818  57.27272727  61.36363636  65.45454545  69.54545455
  73.63636364  77.72727273  81.81818182  85.90909091  90.        ]
Spherical latitude: [-90.         -85.88354756 -81.76762011 -77.65273153 -73.53937377
 -69.42800654 -65.31904771 -61.21286467 -57.10976692 -53.00999983
 -48.91373991 -44.82109151 -40.73208506 -36.6466769  -32.56475052
 -28.4861193  -24.41053065 -20.33767131 -16.26717386 -12.19862422
  -8.13156991  -4.06552908   0.           4.06552908   8.13156991
  12.19862422  16.26717386  20.33767131  24.41053065  28.4861193
  32.56475052  36.6466769   40.73208506  44.82109151  48.91373991
  53.00999983  57.10976692  61.21286467  65.31904771  69.42800654
  73.53937377  77.65273153  81.76762011  85.88354756  90.        ]
Height (m): 481000
Radius (m): [6837752.31424518 6837862.00852276 6838188.80559622 6838725.89905453
 6839462.11303663 6840382.14963168 6841466.92578056 6842693.99069968
 6844038.01293238 6845471.32464188 6846964.50972863 6848487.02178921
 6850007.81781446 6851495.99381081 6852921.40916165 6854255.28745742
 6855470.78263883 6856543.50054635 6857451.9672834  6858178.03712615
 6858707.23400555 6859029.02182043 6859137.         6859029.02182043
 6858707.23400555 6858178.03712615 6857451.9672834  6856543.50054635
 6855470.78263883 6854255.28745742 6852921.40916165 6851495.99381081
 6850007.81781446 6848487.02178921 6846964.50972863 6845471.32464188
 6844038.01293238 6842693.99069968 6841466.92578056 6840382.14963169
 6839462.11303663 6838725.89905453 6838188.80559622 6837862.00852276
 6837752.31424518]

Notice that:

The longitude is the same in both coordinates systems.
The latitude is slightly different except for the poles and equator.
The radius (distance from the center of the ellipsoid) varies even though the height is constant.

Tip

We used the WGS84 ellipsoid here but the workflow is the same for any other oblate ellipsoid. Checkout Available ellipsoids for options.

Geodetic to geocentric spherical using pymap3d#

Boule’s Ellipsoid and Sphere classes can be used with pymap3d for converting between different coordinate systems. While pymap3d defines some ellipsoids internally, you may want to use one from Boule if:

You want to be certain that the parameters used for coordinate conversions and gravity calculations are consistent.
You need to define your own ellipsoid, either because you need different parameters than the built-in ones or they aren’t available in either Boule or pymap3d.

The example below converts between geodetic and geocentric spherical using pymap3d.geodetic2spherical instead of boule.Ellipsoid.geodetic_to_spherical to achieve the same outcome as in the previous example.

import pymap3d

longitude = 40
latitude = np.linspace(-90, 90, 45)
height = 481_000  # ICESat-2 orbit height in meters

latitude_sph, longitude_sph, radius = pymap3d.geodetic2spherical(
    latitude, longitude, height, ell=bl.WGS84,
)

print("Geodetic longitude:", longitude)
print("Spherical longitude:", longitude_sph)
print("Geodetic latitude:", latitude)
print("Spherical latitude:", latitude_sph)
print("Height (m):", height)
print("Radius (m):", radius)

Geodetic longitude: 40
Spherical longitude: 40.0
Geodetic latitude: [-90.         -85.90909091 -81.81818182 -77.72727273 -73.63636364
 -69.54545455 -65.45454545 -61.36363636 -57.27272727 -53.18181818
 -49.09090909 -45.         -40.90909091 -36.81818182 -32.72727273
 -28.63636364 -24.54545455 -20.45454545 -16.36363636 -12.27272727
  -8.18181818  -4.09090909   0.           4.09090909   8.18181818
  12.27272727  16.36363636  20.45454545  24.54545455  28.63636364
  32.72727273  36.81818182  40.90909091  45.          49.09090909
  53.18181818  57.27272727  61.36363636  65.45454545  69.54545455
  73.63636364  77.72727273  81.81818182  85.90909091  90.        ]
Spherical latitude: [-90.         -85.88354756 -81.76762011 -77.65273153 -73.53937377
 -69.42800654 -65.31904771 -61.21286467 -57.10976692 -53.00999983
 -48.91373991 -44.82109151 -40.73208506 -36.6466769  -32.56475052
 -28.4861193  -24.41053065 -20.33767131 -16.26717386 -12.19862422
  -8.13156991  -4.06552908   0.           4.06552908   8.13156991
  12.19862422  16.26717386  20.33767131  24.41053065  28.4861193
  32.56475052  36.6466769   40.73208506  44.82109151  48.91373991
  53.00999983  57.10976692  61.21286467  65.31904771  69.42800654
  73.53937377  77.65273153  81.76762011  85.88354756  90.        ]
Height (m): 481000
Radius (m): [6837752.31424518 6837862.00852276 6838188.80559622 6838725.89905453
 6839462.11303663 6840382.14963169 6841466.92578056 6842693.99069968
 6844038.01293239 6845471.32464188 6846964.50972864 6848487.02178921
 6850007.81781447 6851495.99381081 6852921.40916165 6854255.28745742
 6855470.78263883 6856543.50054635 6857451.9672834  6858178.03712615
 6858707.23400555 6859029.02182044 6859137.         6859029.02182044
 6858707.23400555 6858178.03712615 6857451.9672834  6856543.50054635
 6855470.78263883 6854255.28745742 6852921.40916165 6851495.99381081
 6850007.81781447 6848487.02178921 6846964.50972864 6845471.32464188
 6844038.01293239 6842693.99069968 6841466.92578056 6840382.14963169
 6839462.11303663 6838725.89905453 6838188.80559622 6837862.00852276
 6837752.31424518]

Geocentric spherical to geodetic#

Another common coordinate conversion used in global studies is from geocentric spherical to geodetic coordinates. The example below demonstrate this conversion using the Cartesian coordinates of the Insight lander on Mars from [LeMaistre2023] and the Martian ellipsoid defined in Boule.

import boule as bl
import numpy as np

xyz = [-2_417_504.5, 2_365_954.5, 266_266.7]  # InSight lander coordinates

# convert Cartesian to geocentric spherical
radius = np.linalg.norm(xyz)
latitude_sph = np.rad2deg(np.atan2(xyz[2], np.linalg.norm(xyz[0:2])))
longitude_sph = np.rad2deg(np.atan2(xyz[1], xyz[0]))

mars_ellipsoid = bl.Mars2009

longitude, latitude, height = mars_ellipsoid.spherical_to_geodetic(
  longitude_sph, latitude_sph, radius,
)

print(f"Geodetic longitude: {longitude}")
print(f"Geodetic latitude: {latitude}")
print(f"Ellipsoidal height (m): {height}")

Geodetic longitude: 135.6174366918431
Geodetic latitude: 4.548092430561262
Ellipsoidal height (m): -2241.5527252679644

Coordinate conversions

Contents

Coordinate conversions#

Geodetic to geocentric spherical#

Geodetic to geocentric spherical using pymap3d#

Geocentric spherical to geodetic#