How to use the harmonica.EQLHarmonic function in harmonica

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
# at a constant relative depth beneath each observation point. The damping
# parameter helps smooth the predicted data and ensure stability.
eql = hm.EQLHarmonic(relative_depth=1000, damping=1)

# Fit the layer coefficients to the observed magnetic anomaly.
eql.fit(coordinates, data.total_field_anomaly_nt)

# Evaluate the data fit by calculating an R² score against the observed data.
# This is a measure of how well layer the fits the data NOT how good the
# interpolation will be.
print("R² score:", eql.score(coordinates, data.total_field_anomaly_nt))

# Interpolate data on a regular grid with 500 m spacing. The interpolation
# requires an extra coordinate (upward height). By passing in 1500 m, we're
# effectively upward-continuing the data (mean flight height is 500 m).
grid = eql.grid(spacing=500, data_names=["magnetic_anomaly"], extra_coords=1500)

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)

# 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)

eql = hm.EQLHarmonic(depth=1000, damping=10)
eql.fit(*train)
print("R² score on testing set:", eql.score(*test))

# Interpolate data into the regular grid at 200m above the sea level
grid = eql.grid(spacing=400, data_names=["magnetic_anomaly"], extra_coords=200)
# The grid is a xarray.Dataset with values, coordinates, and metadata
print(grid)

# Plot original magnetic anomaly
maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
tmp = ax1.scatter(data.easting, data.northing, c=data.total_field_anomaly_nt, s=20,
                  vmin=-maxabs, vmax=maxabs, cmap="seismic")
plt.colorbar(tmp, ax=ax1, label="nT")
ax1.set_title("Observed Anomaly Magnetic data from Rio de Janeiro")