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

===========

The geoid is the equipotential surface of the Earth's gravity potential that
coincides with mean sea level. It's often represented by "geoid heights", which
indicate the height of the geoid relative to the reference ellipsoid (WGS84 in
this case). Negative values indicate that the geoid is below the ellipsoid
surface and positive values that it is above. The data are on a regular grid
with 0.5 degree spacing and was generated from the spherical harmonic model
EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the geoid grid
data = hm.datasets.fetch_geoid_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100))
pc = data.geoid.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)
plt.colorbar(
    pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Geoid heights (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()

thing to note is that the ETOPO1 topography is referenced to the geoid, not the
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)