diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py index 642b63e74556f0a21fa304268caf8e0621da440a..c6bf3b3f4242acecadac62359a4fb59520968b96 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 27985a83aee4e672c2ac345e583989b94133a00b..5a4cb4966cf1e8da1a3207ebceee51704d47ed7e 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.