Equivalent sources in spherical coordinates¶

When interpolating gravity or magnetic data over a very large region that covers full continents we need to take into account the curvature of the Earth. In these cases projecting the data to plain Cartesian coordinates may introduce errors due to the distorsions caused by it. Therefore using the `harmonica.EquivalentSources` class is not well suited for it.

Instead, we can make use of equivalent sources that are defined in spherical geocentric coordinates through the `harmonica.EquivalentSourcesSph` class.

Lets start by fetching some gravity data over Southern Africa:

```import ensaio
import pandas as pd

fname = ensaio.fetch_southern_africa_gravity(version=1)
data
```
longitude latitude height_sea_level_m gravity_mgal
0 18.34444 -34.12971 32.2 979656.12
1 18.36028 -34.08833 592.5 979508.21
2 18.37418 -34.19583 18.4 979666.46
3 18.40388 -34.23972 25.0 979671.03
4 18.41112 -34.16444 228.7 979616.11
... ... ... ... ...
14354 21.22500 -17.95833 1053.1 978182.09
14355 21.27500 -17.98333 1033.3 978183.09
14356 21.70833 -17.99166 1041.8 978182.69
14357 21.85000 -17.95833 1033.3 978193.18
14358 21.98333 -17.94166 1022.6 978211.38

14359 rows × 4 columns

To speed up the computations on this simple example we are going to downsample the data using a blocked mean to a resolution of 6 arcminutes.

```import numpy as np
import verde as vd

blocked_mean = vd.BlockReduce(np.mean, spacing=6 / 60, drop_coords=False)
(longitude, latitude, height), gravity_data = blocked_mean.filter(
(data.longitude, data.latitude, data.height_sea_level_m),
data.gravity_mgal,
)

data = pd.DataFrame(
{
"longitude": longitude,
"latitude": latitude,
"height_sea_level_m": height,
"gravity_mgal": gravity_data,
}
)
data
```
longitude latitude height_sea_level_m gravity_mgal
0 19.394500 -34.919505 0.00 979750.15
1 19.465000 -34.953000 0.00 979751.50
2 19.554000 -34.996000 0.00 979750.20
3 19.748000 -34.979000 0.00 979747.00
4 19.299585 -34.832660 0.00 979726.95
... ... ... ... ...
8072 15.896660 -17.395000 1104.30 978169.29
8073 15.973330 -17.390000 1106.40 978170.99
8074 16.068330 -17.390000 1110.25 978168.94
8075 16.140000 -17.390000 1112.50 978168.79
8076 18.408330 -17.391660 1115.90 978157.49

8077 rows × 4 columns

And compute gravity disturbance by subtracting normal gravity using `boule` (see Gravity Disturbance for more information):

```import boule as bl

ellipsoid = bl.WGS84
normal_gravity = ellipsoid.normal_gravity(data.latitude, data.height_sea_level_m)
gravity_disturbance = data.gravity_mgal - normal_gravity
```

Lets define some equivalent sources in spherical coordinates. We will choose some first guess values for the `damping` and `depth` parameters.

```import harmonica as hm

eqs = hm.EquivalentSourcesSph(damping=1e-3, relative_depth=10000)
```

Check how we can estimate the damping and depth parameters using a cross-validation.

Before we can fit the sources’ coefficients we need to convert the data given in geographical coordinates to spherical ones. We can do it through the `boule.Ellipsoid.geodetic_to_spherical` method of the WGS84 ellipsoid defined in `boule`.

```coordinates = ellipsoid.geodetic_to_spherical(
data.longitude, data.latitude, data.height_sea_level_m
)
```

And then use them to fit the sources:

```eqs.fit(coordinates, gravity_disturbance)
```
`EquivalentSourcesSph(damping=0.001, relative_depth=10000)`
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.

We can then use these sources to predict the gravity disturbance on a regular grid defined in geodetic coordinates. We will generate a regular grid of computation points located at the maximum height of the data and with a spacing of 6 arcminutes.

```# Get the bounding region of the data in geodetic coordinates
region = vd.get_region((data.longitude, data.latitude))

# Get the maximum height of the data coordinates
max_height = data.height_sea_level_m.max()

# Define a regular grid of points in geodetic coordinates
grid_coords = vd.grid_coordinates(
region=region, spacing=6 / 60, extra_coords=max_height
)
```

But before we can tell the equivalent sources to predict the field we need to convert the grid coordinates to spherical.

```grid_coords_sph = ellipsoid.geodetic_to_spherical(*grid_coords)
```

And then predict the gravity disturbance on the grid points:

```gridded_disturbance = eqs.predict(grid_coords_sph)
```

Lastly we can generate a `xarray.DataArray` using `verde.make_xarray_grid`:

```grid = vd.make_xarray_grid(
grid_coords,
gridded_disturbance,
data_names=["gravity_disturbance"],
extra_coords_names="upward",
)
grid
```
```<xarray.Dataset>
Dimensions:              (northing: 177, easting: 209)
Coordinates:
* easting              (easting) float64 11.91 12.01 12.11 ... 32.65 32.75
* northing             (northing) float64 -35.0 -34.9 -34.8 ... -17.49 -17.39
upward               (northing, easting) float64 2.622e+03 ... 2.622e+03
Data variables:
gravity_disturbance  (northing, easting) float64 5.957 5.986 ... 3.885 3.892```

Since the data points don’t cover the entire area, we might want to mask those grid points that are too far away from any data point:

```grid = vd.distance_mask(
data_coordinates=(data.longitude, data.latitude), maxdist=0.5, grid=grid
)
```

Lets plot it:

```import matplotlib.pyplot as plt
import cartopy.crs as ccrs

maxabs = vd.maxabs(gravity_disturbance, grid.gravity_disturbance.values)

fig, (ax1, ax2) = plt.subplots(
nrows=1,
ncols=2,
figsize=(12, 9),
sharey=True,
subplot_kw=dict(
projection=ccrs.Mercator(central_longitude=100)
)
)
tmp = ax1.scatter(
data.longitude,
data.latitude,
c=gravity_disturbance,
s=3,
vmin=-maxabs,
vmax=maxabs,
cmap="seismic",
transform=ccrs.PlateCarree(),
)
tmp = grid.gravity_disturbance.plot.pcolormesh(
ax=ax2,
vmin=-maxabs,
vmax=maxabs,
cmap="seismic",
transform=ccrs.PlateCarree(),
)

ax1.gridlines(draw_labels=["bottom", "left"], linewidth=0)
ax1.set_title("Block-median reduced gravity disturbance")
ax2.gridlines(draw_labels=["bottom"], linewidth=0)
ax2.set_title("Gridded gravity disturbance")

for ax in (ax1, ax2):
ax.set_extent(region, crs=ccrs.PlateCarree())
ax.coastlines()
plt.colorbar(
tmp, ax=ax, label="mGal", pad=0.05, aspect=40, orientation="horizontal"
)

plt.show()
```