WIP atmospheric and spherical cap corrections

[craigmillernz]

Created functions for the Bouguer spherical cap correction and atmospheric correction.

[welcome[bot]]

💖 Thanks for opening your first pull request! 💖

Please make sure you read the following:

• Authorship Guidelines: Our rules for giving you credit for your contributions, including authorship on publications and Zenodo archives.
• Contributing Guide: What the review process is like and our infrastructure for testing and documentation.
• Code of Conduct: How we expect people to interact in our projects.

A few things to keep in mind:

• Remember to run make format to make sure your code follows our style guide.
• If you need help writing tests, take a look at the existing ones for inspiration. If you don't know where to start, let us know and we'll walk you through it.
• All new features should be documented. It helps to write the docstrings for your functions/classes before writing the code. This will help you think about your code design and results in better code.
• No matter what, we are really grateful that you put in the effort to do this! ⭐

[leouieda]

@craigmillernz sorry this has lapsed. I'll try to find some time to look at the code but this week and next will be a bit busy as we're finishing the teaching semester here. Thanks for all your work!

[leouieda]

@craigmillernz I've been reading the Hinze et al paper and have a few thoughts:

1. They adopt the geoid as a vertical datum and then correct for the "geophysical indirect effect". This is entirely unnecessary if we adopt the ellipsoid as the vertical datum (which we do when calculating disturbances instead of anomalies). So this is something we can skip.
2. The atmospheric correction in there is an evaluation of an expression for a given atmospheric model. If possible, it would be better to have the actual expression and particular atmospheric model separated in the code so that we could potentially swap the atmosphere if we desire (also to make it clear where those numbers come from). I haven't dug into the references yet so I'm not sure how difficult this would be.
3. The paper mentions "In the revised procedure, to account for the effect of the curvature of the earth, the horizontal slab equation is replaced by the closed-form formula for a spherical cap of radius 166.7km (LaFehr, 1991b)" which sounds to me like they completely replace the infinite plate with a spherical cap. In this case, there is analytical expression for the effect of a cap in Heck and Seitz (2007). This is a general expression that can be used at any observation height, earth radius, and even cap size. It's similar to LaFehr 1991 but doesn't assume that observations are on the topography.

Many of these older papers take these short cuts because of the computing limitations of the time. But in 2020 we don't need to do that any more. So we could potential set new trends by ditching some of these old approximations in favour of analytical expressions.

Given that, what I would propose is:

1. Leave the Bouguer correction as a simple Bouguer planar correction
2. Add a spherical_cap_correction that replaces the Bouguer correction (instead of complementing it). This would implement the full equation from Heck and Seitz.
3. Have the analytical expression for the atmospheric correction with defaults (or maybe a separate variable) for the Earth-related coefficients.

The cap correction is a bit tricky because of differences between spherical and geodetic (ellipsoid) coordinates (the latitude to spherical latitude conversion is height dependent). But we can probably ignore that and assume a spherical approximation (ignoring the flattening and using a mean radius).

The inputs would be: topography (height and lat, lon), heights of observations, reference_radius, density_above (above the reference = topography), density_below (below the reference = water - crust if Earth).

The approximation would be that latitude and spherical latitude are the same. The reference radius can be a constant to assume that Earth is spherical. But it could also be an array to have variable radius with latitude and take into account the flattening (with boule.WGS84.geodetic_to_spherical(longitude, latitude, height=0) for example).

[craigmillernz]

@leouieda Thanks for your reply and thoughts. Just to correct your first point, Hinze 2005 don't recommend using the geoid plus "indirect effect". They recommend using the ellipsoid. Refer to Table 1 where they list the "conventional use" and "recommended revisions". As we are already using the ellipsoid it's all good. Hinze et al do recommend replacing the planar bouguer slab with the spherical cap equivalent. This is what i've coded up at the moment.

I'll take a look at the Heck and Steiz paper for the spherical cap and see what's involved. Is there a particular equation I should look at in that paper?

Can I just clarify what you are proposing, these two points seem a bit contradictory

1. Keep the bouguer planar slab correction as is
2. Implement the spherical cap that replaces the bouguer slab (instead of complementing it) which contradicts point 1.

Do you mean implement both the planar slab correction AND the spherical cap correction then the user can choose which version they would like to use?

I think for "backwards compatibility" with older data it would be wise to keep the planar bouguer slab available for those who want it, but to recommend for new data to use the spherical cap bouguer correction.

Finally for atmospheric correction, do you mean so that say Mars atmosphere could be used instead of Earth?