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

different heights. For a great number of applications we may need to interpolate these
data points into a regular grid. This can be done through an Equivalent Layer
interpolator. We will use :class:`harmonica.EQLHarmonic` to generate a set of point
sources beneath the observation points that predicts the observed data. Then these point
sources will be used to interpolate the data values into a regular grid at a constant
height.
"""
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import verde as vd
import harmonica as hm


# Fetch magnetic anomaly data from Rio de Janeiro
data = hm.datasets.fetch_rio_magnetic()

# Reduce region of the survey to speed things up
region = [-42.35, -42.10, -22.35, -22.15]
inside = vd.inside((data.longitude.values, data.latitude.values), region)
data = data[inside]
print("Number of data points:", data.shape[0])

# Project coordinates
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
data["easting"], data["northing"] = projection(
    data.longitude.values, data.latitude.values
)
coordinates = (data.easting, data.northing, data.altitude_m)

train, test = vd.train_test_split(coordinates, data.total_field_anomaly_nt, random_state=0)

================================================

A subsection from an airborne survey of the state of Rio de Janeiro, Brazil,
conducted in 1978. The data are stored in a :class:`pandas.DataFrame` with
columns: longitude, latitude, total field anomaly (nanoTesla), observation
height above the WGS84 ellipsoid (meters), type of flight line (LINE or TIE),
and flight line number. See the documentation for
:func:`harmonica.datasets.fetch_rio_magnetic` 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_rio_magnetic()
print(data)

# Plot the observations in a Mercator map using Cartopy
fig = plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Magnetic data from Rio de Janeiro", pad=25)
maxabs = vd.maxabs(data.total_field_anomaly_nt)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.total_field_anomaly_nt,
    s=0.8,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
    transform=ccrs.PlateCarree(),