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

:class:`harmonica.EQLHarmonicSpherical` implements the equivalent layer
technique in spherical coordinates. It has the same advantages as the Cartesian
equivalent layer (:class:`harmonica.EQLHarmonic`) while taking into account the
curvature of the Earth.
"""
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import boule as bl
import verde as vd
import harmonica as hm


# Fetch the sample gravity data from South Africa
data = hm.datasets.fetch_south_africa_gravity()

# Downsample the data using a blocked mean to speed-up the computations
# for this example. This is preferred over simply discarding points to avoid
# aliasing effects.
blocked_mean = vd.BlockReduce(np.mean, spacing=0.2, drop_coords=False)
(longitude, latitude, elevation), gravity_data = blocked_mean.filter(
    (data.longitude, data.latitude, data.elevation), data.gravity,
)

# Compute gravity disturbance by removing the gravity of normal Earth
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(latitude, height=elevation)
gravity_disturbance = gravity_data - gamma

# Convert data coordinates from geodetic (longitude, latitude, height) to
# spherical (longitude, spherical_latitude, radius).

Republic of South Africa. The data was made available by the `National Centers
for Environmental Information (NCEI) `__ (formerly
NGDC) and are in the `public domain
`__. The
entire dataset is stored in a :class:`pandas.DataFrame` with columns:
longitude, latitude, elevation (above sea level) and gravity(mGal). See the
documentation for :func:`harmonica.datasets.fetch_south_africa_gravity` for
more information.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import verde as vd
import harmonica as hm

# Fetch the data in a pandas.DataFrame
data = hm.datasets.fetch_south_africa_gravity()
print(data)

# Plot the observations in a Mercator map using Cartopy
fig = plt.figure(figsize=(6.5, 5))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Observed gravity data from South Africa", pad=25)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.gravity,
    s=0.8,
    cmap="viridis",
    transform=ccrs.PlateCarree(),
)
plt.colorbar(
    tmp, ax=ax, label="observed gravity [mGal]", aspect=50, pad=0.1, shrink=0.92