Create and Interpolate Free-air anomaly for comparing to Bouguer Anomlay + Bouguer Disturbance

[geojoeyphillips]

Thank you for these amazing gravity analysis tools. I'm stuck on a relative gravity dataset problem where I want to create a grid for the free-air correction and compare that with a grid of Bouguer anomaly vs. Bouguer disturbance.

To get the free-air on a sphere use 0 instead of height above ellipsoid:
data = data.assign(fa_normal_grav_ellipsoid = bl.WGS84.normal_gravity (data.Latitude, 0))

Then for Bouguer anomaly rather than disturbance:
data = data.assign(gravity_anomaly_mgal = data.Gravity - (data.fa_normal_grav_ellipsoid - 0.3086*(data.Elevation)))

I think for the gravity_anomaly_mgal I can use the equivalent sources method in the tutorial, but I don't think I can to get a free-air grid? I've tried interpolating the data then converting to xr.DataArray but having trouble plotting in PyGMT.

interpolator = scipy.interpolate.NearestNDInterpolator(data[['Longitude', 'Latitude']].values, data['fa_normal_grav_ellipsoid'].values)
y = np.linspace(lat[0],lat[1],100)
x = np.linspace(long[0],long[1],200)
xx, yy = np.meshgrid(x,y)
values = interplator(xx.flat, yy.flat).reshape(xx.shape)
filled = xr.DataArray(values, dims = ['Latitude','Longitude'], coords = [y,x])

Is there a better way to interpolate for the grid to plot in PyGMT? The above attempt doesn't work too well. I know using the ellipsoid to calculate gravity disturbance is the better method, but I would still like to compare the results from each method. I don't think I can use https://www.fatiando.org/boule/v0.4.1/user_guide/gravity_disturbance.html in the first case to look at the free-air correction itself?

Sorry if this is too much to ask I have just been stuck on this problem. Thank you!

[santisoler]

Hi @geojoeyphillips! Thanks for opening this question. There are a few things going on here that I'll try to sort out.

The first one is more conceptual. You mention that:

To get the free-air on a sphere use 0 instead of height above ellipsoid:

data = data.assign(fa_normal_grav_ellipsoid = bl.WGS84.normal_gravity (data.Latitude, 0))

But the output of bl.WGS84.normal_gravity() at zero height doesn't contain any free air contribution: it's just the gravity acceleration of the WGS84 ellipsoid computed on its surface. This is equivalent to using the Somigliana equation in case you are familiar with it.
The gravity anomaly $\Delta g$ is defined as

$\Delta g = g_\text{obs} - \gamma_q$

Where $g_\text{obs}$ is the observed gravity and $\gamma_q$ is the normal gravity on a point located at the same latitude as the observation point but located on the surface of the ellipsoid. So what you are computing in that line is this $\gamma_q$, not the free-air correction nor the free-air anomaly.

After that you compute the Bouguer anomaly as:

Then for Bouguer anomaly rather than disturbance:

data = data.assign(gravity_anomaly_mgal = data.Gravity - (data.fa_normal_grav_ellipsoid - 0.3086*(data.Elevation)))

The free-air anomaly $\Delta g_\text{fa}$ is defined as:

$\Delta g_\text{fa} = g_\text{obs} - (\gamma_q - 0.3086 h)$

Where $h$ is the height of the observation point and the gravity is given in mGal. So, what you are actually computing is the free-air anomaly and not the Bouguer anomaly.

Nevertheless, since you are exploring the differences in using the Bouguer anomaly and the Bouguer disturbance, the last equation provides a good insight: in it you can see that the free-air anomaly is equal to the observed gravity minus a first order approximation of the normal gravity on each observation point. Instead of using any approximation, Boule allow us to compute the analytical solution for it through the normal_gravity method.

---

The second one is more related to the code. You can grid with scipy interpolators, with Verde and if you are working with harmonic fields you could also use equivalent sources through Harmonica. Both with Verde and Harmonica you can use the grid method and get a DataArray without the hassle to build it yourself. In any case you could also use the verde.make_xarray_grid function to create a DataArray from the coordinates arrays and the 2d array of values.

The best interpolation method to use heavily depends on the data you are working with. If you already have a regular grid and you want to register in another grid, neirest neighbouts is enough. If you have gravity disturbances observed at scattered points at significantly different observation heights, then I would recommend you equivalent sources. So I wouldn't recommend you one in favour of the other without having more details about your data.

[santisoler]

Hi @geojoeyphillips! The make_xarray_grid function that not perform any interpolation or gridding. It's just a helper function to easily transform a 2d array with its coordinates into a DataArray. But to use it you already need a 2D grid. So you need to interpolate first and then you can use that function to convert the 2D grid into a DataArray.

Now that I'm reading the docstring of make_xarray_grid I think it's kind of confusing and might lead people to think that it actually performs some kind of interpolation. I'll try to change it to make it clearer.

[geojoeyphillips]

Hi @santisoler thank you for the response! That cleared up my confusion.
grav_anom_grd = vd.ScipyGridder(method = 'cubic').fit(coordinates = coordinates, data = data.grav_anom_mgal,weights = None)
region = vd.get_region((data.Longitude,data.Latitude))
gravity_anomaly_grid = grav_anom_grd.grid(region = region,spacing = 1/6000,projection = projection, dims = ["latitude","longitude"], data_names = "gravity_anomaly")
Worked out perfectly! If I can ask a follow-up question, is there a way to remove the NaN/gray space on the edges of the map region and plot a topography grdimage below? I can ask the PyGMT community as well since I've asked so much already. Thank you so much again!

[santisoler]

Glad to be helpful!

If you are using pygmt.Figure.grdimage for plotting the grids, I think you can make the nans transparent by setting nan_transparent=True.

[geojoeyphillips]

@santisoler Wow can't believe I missed that in the documentation. Thank you very much again! You're a life saver.

[leouieda]

Hi @geojoeyphillips @santisoler just to add a small correction: the classic gravity anomaly is defined as $\Delta g(\lambda, \phi) = g_{obs}(\lambda, \phi, h=geoid) - \gamma_{ell}(\lambda, \phi, h=ellipsoid)$. So the actual observed data aren't used where they were measured and need to be downward continue to the geoid (using the free-air correction) and the Bouguer correction tries to account for the extra mass between the observation point and the geoid. Hence why geoid/orthometric heights are used instead of ellipsoidal/geometric heights. Geodesy has moved on from this concept a bit already but it's still heavily used in geophysics. Notice that $\Delta g$ does not depend on the height since it's a difference between quantities at different heights.