How to use the harmonica.datasets.fetch_topography_earth function in harmonica

reference thickness of the continental crust in order to convert the
root/anti-root thickness into Moho depth. The function contains common default
values for the reference thickness and crust, mantle, and water densities
[TurcotteSchubert2014]_.

We'll use our sample topography data
(:func:`harmonica.datasets.fetch_topography_earth`) to calculate the Airy
isostatic Moho depth of Africa.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the elevation model and cut out the portion of the data corresponding to
# Africa
data = hm.datasets.fetch_topography_earth()
region = (-20, 60, -40, 45)
data_africa = data.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))
print("Topography/bathymetry grid:")
print(data_africa)

# Calculate the isostatic Moho depth using the default values for densities and
# reference Moho
moho = hm.isostasy_airy(data_africa.topography)
print("\nMoho depth grid:")
print(moho)

# Draw the maps
plt.figure(figsize=(8, 9.5))
ax = plt.axes(projection=ccrs.LambertCylindrical(central_longitude=20))
pc = moho.plot.pcolormesh(
    ax=ax, cmap="viridis_r", add_colorbar=False, transform=ccrs.PlateCarree()

"""
Earth Topography
================

The topography and bathymetry of the Earth according to the ETOPO1 model
[AmanteEakins2009]_. The original model has 1 arc-minute grid spacing but here
we downsampled to 0.5 degree grid spacing to save space and download times.
Heights are referenced to sea level.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the topography grid
data = hm.datasets.fetch_topography_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-30))
pc = data.topography.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="terrain"
)
plt.colorbar(
    pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Topography of the Earth (ETOPO1)")
ax.coastlines()
plt.tight_layout()
plt.show()

ellipsoid. Since we want to remove the masses between the surface of the Earth
and ellipsoid, we need to add the geoid height to the topography before Bouguer
correction.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import boule as bl
import harmonica as hm

# Load the global gravity, topography, and geoid grids
data = xr.merge(
    [
        hm.datasets.fetch_gravity_earth(),
        hm.datasets.fetch_geoid_earth(),
        hm.datasets.fetch_topography_earth(),
    ]
)
print(data)

# Calculate normal gravity and the disturbance
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell)
disturbance = data.gravity - gamma

# Reference the topography to the ellipsoid
topography_ell = data.topography + data.geoid

# Calculate the Bouguer planar correction and the topography-free disturbance.
# Use the default densities for the crust and ocean water.
bouguer = hm.bouguer_correction(topography_ell)
disturbance_topofree = disturbance - bouguer