The library#

All functions and classes in Boule are available from the base namespace of the boule package. This means that you can access all of them with a single import:

# Boule is usually imported as bl
import boule as bl


Boule comes with several built-in ellipsoids that are available as global variables in the boule module. The ellipsoids can be printed to view their defining attributes:

Ellipsoid(name='WGS84', semimajor_axis=6378137, flattening=0.0033528106647474805, geocentric_grav_const=398600441800000.0, angular_velocity=7.292115e-05, long_name='World Geodetic System (1984)', reference='Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien\u202f; New York: Springer.', comments=None)
Sphere(name='Moon2015', radius=1737151, geocentric_grav_const=4902800070000.0, angular_velocity=2.6617073e-06, long_name='Moon spheroid (2015)', reference='Wieczorek, MA (2015). 10.05 - Gravity and Topography of the Terrestrial Planets, Treatise of Geophysics (Second Edition); Elsevier. doi:10.1016/B978-0-444-53802-4.00169-X', comments=None)
Ellipsoid(name='Vesta2017', semimajor_axis=278556, flattening=0.17459684946653456, geocentric_grav_const=17288000000.0, angular_velocity=0.0003267, long_name='Vesta reference ellipsoid (2017)', reference='Karimi, R., Azmoudeh Ardalan, A., & Vasheghani Farahani, S. (2017). The size, shape and orientation of the asteroid Vesta based on data from the Dawn mission. Earth and Planetary Science Letters, 475, 71–82. https://doi.org/10.1016/j.epsl.2017.07.033', comments='Geocentric biaxial ellipsoid')

Ellipsoids define a name (short and long version) and reference for the origin of the numbers used:

print(f"{bl.Mars2009.name}: {bl.Mars2009.long_name}")
Mars2009: Mars ellipsoid (2009)
Ardalan, A. A., Karimi, R., & Grafarend, E. W. (2009). A New Reference Equipotential Surface, and Reference Ellipsoid for the Planet Mars. Earth, Moon, and Planets, 106(1), 1. doi:10.1007/s11038-009-9342-7

Other derived properties of ellipsoids are calculated on demand when accessed:



You may have noticed that there are 3 different types of ellipsoids: Ellipsoid, Sphere, and TriaxialEllipsoid. They each offer different attributes and capabilities. Be sure to check out the List of functions and classes (API) for a full list of what each class offers.

Normal gravity#

Ellipsoids can be used for computations generally encountered in geodetic and geophysical applications. A common one is calculating normal gravity. Here is an example of using Boule to calculate the normal gravity of the WGS84 ellipsoid on the Earth’s surface (topography in the continents, the geoid in the oceans).

First, we need to import a few other packages:

import ensaio        # For downloading sample data
import pygmt         # For plotting maps
import xarray as xr  # For manipulating grids

Now we can download and open co-located grids of topography and geoid using ensaio and xarray:

fname_topo = ensaio.fetch_earth_topography(version=1)
fname_geoid = ensaio.fetch_earth_geoid(version=1)
topography = xr.load_dataarray(fname_topo)
geoid = xr.load_dataarray(fname_geoid)
<xarray.DataArray 'geoid' (latitude: 1081, longitude: 2161)> Size: 19MB
array([[-29.5, -29.5, -29.5, ..., -29.5, -29.5, -29.5],
       [-29.5, -29.5, -29.5, ..., -29.5, -29.5, -29.5],
       [-29.6, -29.6, -29.6, ..., -29.6, -29.6, -29.6],
       [ 14.7,  14.7,  14.7, ...,  14.7,  14.7,  14.7],
       [ 15.2,  15.2,  15.2, ...,  15.2,  15.2,  15.2],
       [ 15.4,  15.4,  15.4, ...,  15.4,  15.4,  15.4]])
  * longitude  (longitude) float64 17kB -180.0 -179.8 -179.7 ... 179.8 180.0
  * latitude   (latitude) float64 9kB -90.0 -89.83 -89.67 ... 89.67 89.83 90.0
    Conventions:     CF-1.8
    title:           Geoid height (EIGEN-6C4) with respect to WGS84
    crs:             WGS84
    source:          Generated from the EIGEN-6C4 model by the ICGEM Calculat...
    license:         Creative Commons Attribution 4.0 International Licence
    references:      https://doi.org/10.5880/icgem.2015.1
    long_name:       geoid height
    standard_name:   geoid_height_above_reference_ellipsoid
    description:     height of the geoid with respect to the WGS84 ellipsoid
    units:           m
    actual_range:    [-106.5   86. ]
    icgem_metadata:  generating_institute: gfz-potsdam\ngenerating_date: 2021...

The computation height can be defined by combining topography and geoid:

height = xr.where(
    topography >= 0,
    topography + geoid,  # geometric height of topography in the continents
    geoid,  # geoid height in the oceans

Finally, we can calculate normal gravity using normal_gravity at the given heights and plot it on a map with pygmt:

gamma = bl.WGS84.normal_gravity(topography.latitude, height)

fig = pygmt.Figure()
fig.grdimage(gamma, projection="W20c", cmap="viridis", shading="+a45+nt0.3")
fig.basemap(frame=["af", "WEsn"])
fig.colorbar(position="JCB+w10c", frame=["af", 'y+l"mGal"'])

See also

Normal gravity provides a more detailed tutorial, including the different definitions of normal gravity for each ellipsoid type.