From 83c15d2997ca29dd3756b5fdbccc983e0d606a90 Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Tue, 17 Mar 2020 15:34:48 -0600
Subject: [PATCH] WIP: Fix data type/mathematical errors
---
geomagio/residual/Calculation.py | 42 ++++++++++++++++++++------------
geomagio/residual/Reading.py | 8 +++---
2 files changed, 31 insertions(+), 19 deletions(-)
diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py
index 642b63e74..c6bf3b3f4 100644
--- a/geomagio/residual/Calculation.py
+++ b/geomagio/residual/Calculation.py
@@ -23,11 +23,15 @@ def calculate_I(measurements, ordinates, metadata):
northdown_ordinate = average_ordinate(ordinates, "NorthDown")
northup_ordinate = average_ordinate(ordinates, "NorthUp")
southup_ordinate = average_ordinate(ordinates, "SouthUp")
- total_ordinate = Ordinate()
- total_ordinate.h, total_ordinate.z, total_ordinate.e, total_ordinate.f = np.average(
- [southdown_ordinate, southup_ordinate, northdown_ordinate, northup_ordinate],
- axis=1,
- )
+ # gather ordinates into array
+ ordinates = [
+ southdown_ordinate,
+ northdown_ordinate,
+ northup_ordinate,
+ southup_ordinate,
+ ]
+
+ total_ordinate = average_ordinate(ordinates, None)
# calculate f for each measurement type
southdown_f = calculate_f(southdown_ordinate, total_ordinate, Iprime)
southup_f = calculate_f(southup_ordinate, total_ordinate, Iprime)
@@ -90,7 +94,7 @@ def calculate_scale(f, measurements, I, pier_correction):
time_delta = np.diff([m.time for m in measurements])
- delta_b = delta_f + time_delta
+ delta_b = delta_f + (time_delta / 60.0)
scale_value = f * np.deg2rad(delta_b / detla_r)
@@ -100,26 +104,32 @@ def calculate_scale(f, measurements, I, pier_correction):
def average_angle(measurements, type):
if type == "NorthDown":
# exclude final measurement, which is only used for scaling
- measurements = measurements[:-1]
- return np.average([m.angle for m in measurements[type]])
+ measurements = measurements[type][:-1]
+ else:
+ measurements = measurements[type]
+ return np.average([m.angle for m in measurements])
def average_residual(measurements, type):
if type == "NorthDown":
# exclude final measurement, which is only used for scaling
- measurements = measurements[:-1]
- return np.average([m.residual for m in measurements[type]])
+ measurements = measurements[type][:-1]
+ else:
+ measurements = measurements[type]
+ return np.average([m.residual for m in measurements])
def average_ordinate(ordinates, type):
+ ordinates = ordinates
if type == "NorthDown":
# exclude final measurement, which is only used for scaling
- ordinates = ordinates[:-1]
- ordinate = Ordinate()
- ordinate.h, ordinate.e, ordinate.z, ordinate.f = np.average(
- [o for o in ordinates[type]], axis=1
- )
- return ordinate
+ ordinates = ordinates[type][:-1]
+ elif type is not None:
+ ordinates = ordinates[type]
+ o = Ordinate(measurement_type=type)
+ avgs = np.average([[o.h, o.e, o.z, o.f] for o in ordinates], axis=0)
+ o.h, o.e, o.z, o.f = avgs
+ return o
def calculate_f(ordinate, total_ordinate, I):
diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py
index 27985a83a..5a4cb4966 100644
--- a/geomagio/residual/Reading.py
+++ b/geomagio/residual/Reading.py
@@ -43,8 +43,8 @@ class Reading(BaseModel):
def calculate(self):
# gather class object to perform calculations
metadata = self.metadata
- ordinates = self.ordinate_index
- measurements = self.measurement_index
+ ordinates = self.ordinate_index()
+ measurements = self.measurement_index()
# calculate inclination
inclination, f, ordinate = calculate_I(measurements, ordinates, metadata)
# calculate absolutes
@@ -54,7 +54,9 @@ class Reading(BaseModel):
# calculate baselines
Hb, Zb = calculate_baselines(Habs, Zabs, ordinate)
# calculate scale value for declination
- calculate_scale(f, measurements, inclination, metadata["pier_correction"])
+ calculate_scale(
+ f, measurements["NorthDown"], inclination, metadata["pier_correction"]
+ )
def measurement_index(self) -> Dict[MeasurementType, List[Measurement]]:
"""Generate index of measurements keyed by MeasurementType.
--
GitLab