Full gravity corrections for Bouguer Anomaly

[craigmillernz]

I'm happy to code up the remaining gravity corrections for the Bouguer Anomaly as per Hinze 2005 http://library.seg.org/doi/10.1190/1.1988183

including : spherical cap correction for Bouguer slab atmospheric correction free air correction

Could either just stop at the correction calculation or calculate the anomalies, including free air, simple Bouguer (no terrain correction) anomaly and complete Bouguer anomaly (including terrain). As terrain correction is a work in progress perhaps just the free air and simple Bouguer? This would require the normal gravity as calculated in Boule.

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible.

You might also want to take a look at our Contributing Guide and Code of Conduct.

[leouieda]

Hi @craigmillernz good to see you here :slightly_smiling_face: It would be great to have a more complete pipeline in Harmonica since that is missing from most open-source software for gravity.

The spherical cap and atmospheric corrections would be most useful I think. Would be willing to code those up? For a general layout, have a look at the Bouguer correction function.

The free-air correction is unnecessary in this case since we use the closed-from normal gravity formula in Boule and can calculate at any height. That's much better than calculating it on the ellipsoid and applying the first- or second-order correction.

Could either just stop at the correction calculation or calculate the anomalies, including free air, simple Bouguer (no terrain correction) anomaly and complete Bouguer anomaly (including terrain).

A particular design decision we made with Harmonica is to only implement the corrections. So a user has to know what they are using and can't simply call bg = hm.bouguer_anomaly(raw_gravity). The reason is that there are a lot of assumptions in the corrections (which ellipsoid, Bouguer density, etc) that are crucial to understanding the final product. We can explain in the form of tutorials how to go about getting things like the Bouguer anomaly (which is not even really an anomaly but a disturbance in our case). That has the added bonus of helping spread accurate and updated information about gravimetry, which is sadly absent from most textbooks about general geophysics.

[craigmillernz]

Thanks @Leo, will do. Agreed about only the corrections, and leaving the user to assemble the parts, so they know what has gone in. I'd like to include the free air correction just for completeness, to allow people to use either approach. I'll include a note to say it's not necessary if the normal gravity is calculated atthe station elevation. A tutorial would be good to show how to assemble the Bouguer anomaly.

[craigmillernz]

Can you please give me permissions to make a pull request?

[leouieda]

I'd like to include the free air correction just for completeness, to allow people to use either approach. I'll include a note to say it's not necessary if the normal gravity is calculated at the station elevation.

I understand the desire for completeness but I also feel strange about implementing something and putting a warning telling people not to use it. Should we still be doing free-air corrections in 2020 when the analytical expression exists?

@santisoler @birocoles any opinions on this?

Can you please give me permissions to make a pull request?

Go right ahead :slightly_smiling_face: You don't need any special permissions for making PRs. Are you getting an error?

[santisoler]

Thanks @craigmillernz for opening this issue. We are glad to get help to take Harmonica to the next level.

I'd like to include the free air correction just for completeness, to allow people to use either approach. I'll include a note to say it's not necessary if the normal gravity is calculated at the station elevation.

I understand the desire for completeness but I also feel strange about implementing something and putting a warning telling people not to use it. Should we still be doing free-air corrections in 2020 when the analytical expression exists?

@santisoler @birocoles any opinions on this?

I always felt annoyed with the concept of free-air correction. Firstly, I think it's a misleading concept when trying to model Earth's interior (specially when you are starting to learn gravity observations): it's usually misunderstood as "bringing the observed point down to the ellipsoid". Secondly, the free-air anomaly is a first-order approximation of the normal gravity function (Blakely, 1995, p. 140), therefore by using it we are introducing errors in our gravity data process. Computing the normal gravity through its analytical form it's done very quickly nowadays, so we would be introducing errors without any benefit.

Moreover, free-air anomaly (as any approximation we make) should be accompanied by an error estimation. That's why, for example, tesseroids forward models include adaptative discretization algorithms that guarantee a certain accuracy threshold by comparing their effect with analytical solutions.

So, I think we should avoid adding a function for it if we can. Maybe would be nice to have it for educational purposes, but I would only use it to test and show how far this approximation is valid. In that case I would create a gallery example or tutorial including a simple free-air correction function for those purposes.

I'm glad this discussion was raised. Would be nice to hear some benefit of the free-air correction I ignore.

[craigmillernz]

I can think of one use case for a stand alone free air correction and that is for microgravity (timelapse) correcting for station elevation change (ground deformation). In which case you are entering the delta_h (station elevation change), not the "topography". Topography is a bit of an odd term to use here I think? Would "elevation" be better?

For Bouguer corrections i'm happy to leave it out but perhaps Boule documentation needs to be explicit that free air isn't needed as the normal gravity is calculated at station elevation, not on ellipsoid. It certainly had me puzzled the first time I used it (probably my own ignorance), when I was comparing to my own derivation which uses the ellipsoid surface. I'm not sure many text books are using the analytic expression, but still refer to the free air correction (even a 2nd order free air that is latitude dependent).

[craigmillernz]

I've added an eotvos corrections as well for moving platforms.

[craigmillernz]

Go right ahead 🙂 You don't need any special permissions for making PRs. Are you getting an error?

Sorted now, pilot error... Once we decide about the free air correction or not I can make a pull request.

[santisoler]

I can think of one use case for a stand alone free air correction and that is for microgravity (timelapse) correcting for station elevation change (ground deformation). In which case you are entering the delta_h (station elevation change), not the "topography".

In this example you have two gravity measurements on the same location (longitude, latitude) at different times and with slightly different observation heights, right? In that case, isn't free-air correction easily replaceable by the analytic normal gravity? I would process each measurement independently, so the difference on the elevation would be accounted in the normal gravity analytical form. For example:

import boule as bl

# Both observation points are located on the same longitude, latitude point,
# but with different heights
coordinates_1 = (longitude, latitude, height_1)
coordinates_2 = (longitude, latitude, height_1)

# The observed gravity on these points is (in mGal)
gravity_1 = 978001.3  # some random values
gravity_2 = 977999.5

# Compute gravity disturbances by removing the normal gravity
# on each observation point
ellipsoid = bl.WGS84
gravity_disturbance_1 = gravity_1 - ellipsoid.normal_gravity(latitude, height_1)
gravity_disturbance_2 = gravity_2 - ellipsoid.normal_gravity(latitude, height_1)

This way we are obtaining the gravity disturbance on both locations, although they are not comparable because of the height difference. But for that case I think we should up or down continue the gravity disturbance from either one of the two measurements, right?

So, I'm missing the spot where free-air correction becomes necessary. What data process would you use?

Topography is a bit of an odd term to use here I think? Would "elevation" be better?

Yes, for sure, we are not limiting our computations to ground surveys. Height of measurement can either be on top of the relief, at airborne heights or at satellite orbits (or even inside the Earth). Some authors refer to elevation when they are talking about height above the geoid and height when talking about ellipsoidal height. In Harmonica we defined three coordinate systems: ellipsoidal coordinates, spherical geocentric coordinates and Cartesian projected coordinates. The elevation is called the upward coordinate both in ellipsoidal and Cartesian coordinates, while the "vertical" coordinate in spherical is radius.

For Bouguer corrections i'm happy to leave it out but perhaps Boule documentation needs to be explicit that free air isn't needed as the normal gravity is calculated at station elevation, not on ellipsoid. It certainly had me puzzled the first time I used it (probably my own ignorance), when I was comparing to my own derivation which uses the ellipsoid surface.

This kind of feedback is really needed in Fatiando, thanks a lot! We maintainers are so into Fatiando's code that some things we consider to be absolutely clear may not be for other people. Maybe we could point to the Harmonica's Gravity Disturbances gallery example on the Boule's Normal Gravity documentation.

I'm not sure many text books are using the analytic expression, but still refer to the free air correction (even a 2nd order free air that is latitude dependent).

You are right about this. The normal gravity analytical form vs free-air anomaly is not a subject that many text books mention, specially the old ones. In my path of learning this stuff, I found very interesting this paper from Li and Götze (2001). They even mention that the second order approximation is accurate enough,. But again, if we can use an analytical form, why approximate it?

Sorted now, pilot error... Once we decide about the free air correction or not I can make a pull request.

Maybe would be nice to have a single pull-request per feature, so feel free to open one for any feature you want to develop, so you don't need to wait until we decide what to do with free-air correction. Maybe we don't mention this very often, but we usually start opening PRs including WIP on their title just to show they are work in progress. This way we can get feedback before we consider it to be ready. In summary, don't wait for the code to be absolutely right to open the PR.

[craigmillernz]

Alright first pull request made for spherical cap and atmospheric corrections.

I'll need to put some thought into your suggestion to use the normal gravity change in lieu of free air, it's an interesting idea, but i'll need to think if there are any fishhooks.

[leouieda]

perhaps Boule documentation needs to be explicit that free air isn't needed as the normal gravity is calculated at station elevation, not on ellipsoid.

Isn't that always the case :slightly_smiling_face:

I'm not sure many text books are using the analytic expression, but still refer to the free air correction (even a 2nd order free air that is latitude dependent).

None of them do. Which is a big problem since almost none of the available software do either. Actually, most people don't even use any "software" for this, just Excel. Since "gravity is easy" there is no need to think about. If only...

I'll need to put some thought into your suggestion to use the normal gravity change in lieu of free air, it's an interesting idea, but i'll need to think if there are any fishhooks.

That is something to think about. If you find that there really is a use case for a pure free-air correction, we could add a second order one and explain in the docs the use cases.