diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py index 0d1cec2156a4c0cc5d6bcd60b5f353ef71227ff5..ddc108a4daeecd4e9007bb1d49de9cc63e040ca6 100644 --- a/geomagio/residual/Calculation.py +++ b/geomagio/residual/Calculation.py @@ -11,6 +11,7 @@ from .MeasurementType import ( MARK_TYPES, ) from .Measurement import AverageMeasurement, Measurement, average_measurement +from .Diagnostics import Diagnostics from .Reading import Reading @@ -29,7 +30,7 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading: # reference measurement, used to adjust absolutes reference = reading[mt.WEST_DOWN][0] # calculate inclination - inclination, f, mean = calculate_I( + inclination, f, i_mean = calculate_I( hemisphere=reading.hemisphere, measurements=reading.measurements ) corrected_f = f + reading.pier_correction @@ -37,15 +38,20 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading: absoluteH, absoluteZ = calculate_HZ_absolutes( corrected_f=corrected_f, inclination=inclination, - mean=mean, + mean=i_mean, reference=adjust_reference and reference or None, ) - absoluteD = calculate_D_absolute( + absoluteD, meridian = calculate_D_absolute( azimuth=reading.azimuth, h_baseline=absoluteH.baseline, measurements=reading.measurements, reference=adjust_reference and reference or None, ) + # populate diagnostics object with averaged measurements + diagnostics = Diagnostics( + inclination=inclination, + meridian=meridian, + ) # calculate scale scale_value = None scale_measurements = reading[mt.NORTH_DOWN_SCALE] @@ -59,8 +65,9 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading: calculated = Reading( absolutes=[absoluteD, absoluteH, absoluteZ], scale_value=scale_value, + diagnostics=diagnostics, # copy other attributes - **reading.dict(exclude={"absolutes", "scale_value"}), + **reading.dict(exclude={"absolutes", "scale_value", "diagnostics"}), ) return calculated @@ -70,7 +77,7 @@ def calculate_D_absolute( azimuth: float, h_baseline: float, reference: Measurement, -) -> Absolute: +) -> Tuple[Absolute, Diagnostics]: """Calculate D absolute. Parameters @@ -87,14 +94,14 @@ def calculate_D_absolute( mean = average_measurement(measurements, DECLINATION_TYPES) reference = reference or mean # average mark - average_mark = average_measurement(measurements, MARK_TYPES).angle + average_mark = average_measurement(measurements, MARK_TYPES) # adjust based on which is larger mark_up = average_measurement(measurements, [mt.FIRST_MARK_UP]).angle mark_down = average_measurement(measurements, [mt.FIRST_MARK_DOWN]).angle if mark_up < mark_down: - average_mark += 90 + average_mark.angle += 90 else: - average_mark -= 90 + average_mark.angle -= 90 # declination measurements declination_measurements = [ average_measurement(measurements, [t]) for t in DECLINATION_TYPES @@ -116,16 +123,20 @@ def calculate_D_absolute( shift = -180 # add subtract average mark angle from average meridian angle and add # azimuth to get the declination baseline - d_b = (meridian - average_mark) + azimuth + shift + d_b = (meridian - average_mark.angle) + azimuth + shift # calculate absolute d_abs = d_b + np.degrees(np.arctan(reference.e / (reference.h + h_baseline))) - return Absolute( - element="D", - absolute=d_abs, - baseline=d_b, - shift=shift, - starttime=mean.time, - endtime=mean.endtime, + + return ( + Absolute( + element="D", + absolute=d_abs, + baseline=d_b, + shift=shift, + starttime=mean.time, + endtime=mean.endtime, + ), + meridian, ) diff --git a/geomagio/residual/Diagnostics.py b/geomagio/residual/Diagnostics.py new file mode 100644 index 0000000000000000000000000000000000000000..fd9dac1c9e515e1e1033f5fb8735272ee15e1d57 --- /dev/null +++ b/geomagio/residual/Diagnostics.py @@ -0,0 +1,17 @@ +from typing import Optional + +from .Measurement import Measurement +from pydantic import BaseModel + + +class Diagnostics(BaseModel): + """Computed diagnostics during calculation. + + Attributes + ---------- + inclination: Average of inclination measurements + meridian: Calculated meridian value + """ + + inclination: float + meridian: float diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py index c91fcd6cc91e3c365bacdb1f65ae52cbe37a8c79..06c8efb975310d9590e7c422787153614baae5c9 100644 --- a/geomagio/residual/Reading.py +++ b/geomagio/residual/Reading.py @@ -9,6 +9,7 @@ from .. import TimeseriesUtility from ..TimeseriesFactory import TimeseriesFactory from .Absolute import Absolute from .Measurement import Measurement, average_measurement +from .Diagnostics import Diagnostics from .MeasurementType import MeasurementType @@ -33,6 +34,7 @@ class Reading(BaseModel): metadata: Dict = {} pier_correction: float = 0 scale_value: float = None + diagnostics: Diagnostics = None def __getitem__(self, measurement_type: MeasurementType): """Provide access to measurements by type. diff --git a/geomagio/residual/SpreadsheetAbsolutesFactory.py b/geomagio/residual/SpreadsheetAbsolutesFactory.py index db8283787320657e5ffd045aca904393f4ef41dd..1f50f794faff107c7663f9c7db9ba947749a1c1b 100644 --- a/geomagio/residual/SpreadsheetAbsolutesFactory.py +++ b/geomagio/residual/SpreadsheetAbsolutesFactory.py @@ -6,6 +6,13 @@ from obspy.core import UTCDateTime import openpyxl from .Absolute import Absolute +from .Calculation import ( + DECLINATION_TYPES, + MARK_TYPES, + INCLINATION_TYPES, + average_measurement, +) +from .Diagnostics import Diagnostics from .Measurement import Measurement from .MeasurementType import MeasurementType as mt from .Reading import Reading @@ -299,6 +306,7 @@ class SpreadsheetAbsolutesFactory(object): metadata=metadata, pier_correction=metadata["pier_correction"], scale_value=numpy.degrees(metadata["scale_value"]), + diagnostics=self._parse_diagnostics(calculation_sheet), ) def _parse_absolutes( @@ -402,6 +410,18 @@ class SpreadsheetAbsolutesFactory(object): "precision": measurement_sheet["H8"].value, } + def _parse_diagnostics( + self, + sheet: openpyxl.worksheet, + ) -> Diagnostics: + """ + Gather diagnostics from list of measurements + """ + return Diagnostics( + inclination=sheet["H40"].value, + meridian=sheet["E36"].value, + ) + def convert_precision(angle, precision="DMS"): """