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

The advantage of using an equivalent layer is that it takes into account the 3D
nature of the observations, not just their horizontal positions. It also allows
data uncertainty to be taken into account and noise to be suppressed though the
least-squares fitting process. The main disadvantage is the increased
computational load (both in terms of time and memory).
"""
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import verde as vd
import harmonica as hm


# Fetch the sample total-field magnetic anomaly data from Great Britain
data = hm.datasets.fetch_britain_magnetic()

# Slice a smaller portion of the survey data to speed-up calculations for this
# example
region = [-5.5, -4.7, 57.8, 58.5]
inside = vd.inside((data.longitude, data.latitude), region)
data = data[inside]
print("Number of data points:", data.shape[0])
print("Mean height of observations:", data.altitude_m.mean())

# Since this is a small area, we'll project our data and use Cartesian
# coordinates
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
easting, northing = projection(data.longitude.values, data.latitude.values)
coordinates = (easting, northing, data.altitude_m)

# Create the equivalent layer. We'll use the default point source configuration

(epsg:27700) coordinate system to WGS84 (epsg:4326) using to_crs function in
GeoPandas.

See the original data for more processing information.

If the file isn't already in your data directory, it will be downloaded
automatically.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import verde as vd
import harmonica as hm
import numpy as np

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

# Plot the observations in a Mercator map using Cartopy
fig = plt.figure(figsize=(7.5, 10))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Magnetic data from Great Britain", pad=25)
maxabs = np.percentile(data.total_field_anomaly_nt, 99)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.total_field_anomaly_nt,
    s=0.001,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
    transform=ccrs.PlateCarree(),