WIP atmospheric and spherical cap corrections

[craigmillernz]

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

[welcome]

💖 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! ⭐

[stickler-ci]

W291 trailing whitespace

[stickler-ci]

E302 expected 2 blank lines, found 1

[stickler-ci]

W293 blank line contains whitespace

[stickler-ci]

W291 trailing whitespace

[stickler-ci]

W293 blank line contains whitespace

[stickler-ci]

E222 multiple spaces after operator

[stickler-ci]

W293 blank line contains whitespace

[stickler-ci]

W293 blank line contains whitespace

[stickler-ci]

W293 blank line contains whitespace

[stickler-ci]

W391 blank line at end of file
W293 blank line contains whitespace
E303 too many blank lines (5)

[leouieda]

@craigmillernz any reason for this? The r is usually needed when writing math to avoid having to use \\ instead of \.

[craigmillernz]

opps my mistake, I thought it a typo

[leouieda]

Suggested change

	"""
	r"""

[leouieda]

Is this a correction applied to the Bouguer correction or the effect of a spherical cap? In other words, does this require a Bouguer correction first or can be it directly applied to the gravity disturbance?

[craigmillernz]

It's my understanding that the spherical cap correction is applied to the Bouguer slab correction. The Bouguer correction has 3 parts, the infinite slab (2piGph, also known as Bullard A), the spherical cap correction (also known as Bullard B) and the terrain correction (Bullard C).

However we can reformulate it to include the spherical cap correction in the Bouguer slab correction, so that there is only one "Bouguer correction" step to apply.

[leouieda]

Could we use the analytical formula for the cap instead of this one? That way, the size of the cap can be a parameter instead of restricted to 167km (which might not apply for other planetary bodies or for satellite data).

[craigmillernz]

Yes that would be better, you mean as defined by LaFehr 1991 (Geophysics 56)? That allows for density, cap radius and earth radius to be specified.
I've coded this up and compares as follows:

1. Using the approximation (present method)

spherical_bouguer_cap_correction(np.array(500))
Out[137]: 0.643743325

2. Using the full analytic

spherical_bouguer_cap_correction_full(np.array(500))
Out[138]: 0.6440463304229793

The approximation assumes a constant density of 2670 kg/m3, so when a different reduction density is used the affect is more apparent. shown here for 2400 kg/m3 density

spherical_bouguer_cap_correction_full(np.array(500), 2400)
Out[143]: 0.5789180498184083

when allowing for different ellipsoids and cap extents.

spherical_bouguer_cap_correction_full(np.array(500), 2400, 1040, planet_radius=bl.MARS.mean_radius, cap_extent=120000)
Out[150]: 0.7782434040018343

[craigmillernz]

reformulated to include the spherical cap correction in the "Bouguer slab formulation" to make it one step

 bouguer_correction(np.array(500))
Out[16]: 55.97179595519349

The spherical_bouguer_cap_correction now returns both the correction and the spherical corrected infinite slab

spherical_bouguer_cap_correction(np.array(500))
Out[17]: (0.6440463304229793, 56.61589317141049)

and the logic check between the infinite slab value (55.971) plus the spherical cap correction (0.644) equals the spherical bouguer correction (56.615) checks out.

[craigmillernz]

Let me know your preference for a 1 step or 2 step correction.

[leouieda]

Right, I makes more sense now. I'm leaning towards a 1 step process since this is only ever used with the Bouguer correction anyway. We could include an argument like spherical=False to turn the correction on or off.

[craigmillernz]

I agree. What I proposed is a single function called bouger_correction() where if spherical=True is passed you get the full spherical bouguer correction (the 56.6 mGal value in the above example), otherwise you get the regular infinite slab correction (the 55.9 mGal value). We could make spherical=True the default and if spherical=False is passed you get the inifinite slab bouguer_correction. As you say, no need to ever return just the spherical cap correction (0.64 mGal value) as you'd only ever need that to apply to the infinite slab bouguer correction.

[leouieda]

This strikes me as a bit odd since the effect of a cap depends on 2 things: the height of the observation point and the thickness of the cap (i.e., height of topography). Is this equation assuming that they are the same?

In fact, it should actually depend on the radial coordinate for both of these, not really the heights (which the authors probably assumed as some mean Earth radius). This will cause the equation to fail for lunar or martian applications, where the reference radius is much smaller.

This can quickly get very complicate if we assume that the heights are relative to the ellipsoid not a sphere. One compromise could be:

• We assume a spherical approximation for this corrections. Height should be reference to the ellipsoid but we'll approximate it for a sphere with a given radius. The differences should be relatively small.
• Require passing the radius of this sphere. For the Earth, we can pass boule.WGS84.mean_radius easily. It might be best to not default to this value realatively sas it would introduce errors if users forget to specify it when applying to other planets.
• Allow passing the extent of the cap, which could default to 167 km.

[craigmillernz]

See comments above about using the full analytic expression. I think if we use that then we can include both the planet radius and cap extent (though there is pretty good evidence that beyond 167 km it gets more imprecise as average crustal densities are likely to vary over large scales.) See Figure 4 in LaFehr 1991.

The full expression does assume a constant radius and is most valid at mid latitudes where radius changes are small. If we have a method of calculating the specific raidus at the station location, then that could be used instead, but Figure 6 from LaFehr shows that it's only at high altitudes (4000m +) that the effect becomes more than instrument error (i.e max of ~20 uGal difference between pole and equator at 4000m), so unless you have a survey that covers the hemisphere I don't think it matter.

[leouieda]

I'm still a bit confused here about the radius referred in the papers. Is that assuming that observations are located at the surface of topography? In that case, the radius of the observation point and the radius for the spherical topography would be the same. But if we have observations from airplanes or satellites, this wouldn't work. Right?

But then again, the Bouguer correction itself doesn't work in these cases anyway.

Either way, I like the full analytic expression! It's good to have these things be explicit. I just had a bug in some code doing a Bouguer correction for Mars and forgetting to set density_water=0 since there are no oceans on Mars. Having defaults can be very dangerous.

[craigmillernz]

I'm not too familiar with airborne data processing (and even less satellite) but i did some googling and it seems airborne gravity is processed similar to ground gravity (except it requires the EotVos correction of course, which we could include as a separate function). See report below for processing steps:

http://www.geosciencebc.com/i/project_data/QUESTdata/GBCReport2008-8/Gravity_Technical_Report.pdf

There is a "curvature" correction that requires the input of the ground surface (not the aircraft altitude). When I computed this it is approximately the same as the spherical_bouguer_cap_correction, (0.64368 mGal).

So i think that for airborne data you would calculate the "normal gravity" at the aircraft elevation (which is where the sensor is) and then calculate the spherical bouguer slab at the elevation of the ground below the aircraft (which can be dervied from the laser altimeter data etc (flight height - laser altimeter = ground height).

In summary, it seems we can use the full spherical bouguer for both ground and airborne gravity. The user needs to ensure they use the correct elevations (aircraft vs ground) for the various parts of the complete Bouguer Anomaly formulation.

[leouieda]

In summary, it seems we can use the full spherical bouguer for both ground and airborne gravity. The user needs to ensure they use the correct elevations (aircraft vs ground) for the various parts of the complete Bouguer Anomaly formulation.

Yes. The problem I mentioned is that gravity should decay with observation height but in all of these different Bouguer corrections the observation height isn't a parameter. That is because the Bouguer plate has constant gravity with height due to being infinite. This is an OK approximation if the observation is on top of the topography, but it's an over-correction when the observations are at greater height (airborne).

I was wondering if the spherical version took this into account. But apparently it doesn't.

The user needs to ensure they use the correct elevations (aircraft vs ground) for the various parts of the complete Bouguer Anomaly formulation.

Yes! This is why we're avoiding the word "height" when talking about Bouguer and using "topography" instead to make it clear that this is not the observation height but the topographic height.

[leouieda]

Does the density not enter into this equation?

[craigmillernz]

See comments above about using the full analytic expression. This approximation assumes 2670, so need to scale to other density.

[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?