From 559314808f79a61dc44b1d1f0720386a191272cb Mon Sep 17 00:00:00 2001
From: "E. Joshua Rigler" <erigler@usgs.gov>
Date: Fri, 21 Feb 2025 16:48:05 -0700
Subject: [PATCH] Remove auto-sift, allow manual shift in calculate_D_absolute
Turns out there are operational limitations that don't allow a
deterministic removal of the 180 degree ambiguity that exists
for calculated Declination. The one previously implemented only
worked for a subset of observatories. Now, we force all calculated
Declinations to point north(ish), which is valid for most of the
globe, and all of the USGS observatories. However, if we ever
take absolute measurements in a location where the local magnetic
field points south(ish), we can simply provide a "shift" of +/-
180 degrees.
---
geomagio/residual/Calculation.py | 25 +++++++++++++++++++++----
1 file changed, 21 insertions(+), 4 deletions(-)
diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py
index d912d581..55a40d31 100644
--- a/geomagio/residual/Calculation.py
+++ b/geomagio/residual/Calculation.py
@@ -59,6 +59,7 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading:
h_baseline=absoluteH.baseline,
measurements=reading.measurements,
reference=reference,
+ shift=reading.absolutes[0].shift,
)
d_computed = calculate_D_computed(
@@ -145,6 +146,7 @@ def calculate_D_absolute(
azimuth: float,
h_baseline: float,
reference: Measurement,
+ shift: float = 0,
) -> Tuple[Absolute, Diagnostics]:
"""Calculate D absolute.
@@ -154,6 +156,7 @@ def calculate_D_absolute(
azimuth: azimuth to mark in decimal degrees.
h_baseline: calculated H baseline value.
reference: reference measurement used to adjust absolute.
+ shift: usually 0 or +/-180 for Declination
Returns
-------
@@ -192,12 +195,26 @@ def calculate_D_absolute(
for m in declination_measurements
]
)
- shift = 0
- if azimuth > 180:
- shift = -180
+
# add subtract average mark angle from average meridian angle and add
# azimuth to get the declination baseline
- d_b = (meridian - average_mark.angle) + azimuth + shift
+ d_b = (meridian - average_mark.angle) + azimuth
+
+ # given operational limitations, there is always an 180-degree ambiguity
+ # in our calculated d_b; therefore, we have to assume the magnetic field
+ # points north(ish), and force d_b to fall between -90 and +90 degrees
+ while d_b < -90:
+ d_b += 180
+ while d_b > 90:
+ d_b -= 180
+
+ # the preceeding lead to a wrong answer if the magnetic field actually
+ # points south(ish);
+ d_b += shift
+
+ # (for example, there are places in Antarctica where the magnetic field
+ # actually points south-ish)
+
# calculate absolute
d_abs = d_b + np.degrees(np.arctan(reference.e / (reference.h + h_baseline)))
--
GitLab