diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py
index f8b8999ef87531ddb8811e09e3287cda8e5b1d85..ea24041dc1b597a1e33165468f3dec6c720b36c1 100644
--- a/geomagio/adjusted/AdjustedMatrix.py
+++ b/geomagio/adjusted/AdjustedMatrix.py
@@ -1,11 +1,12 @@
import numpy as np
from obspy import UTCDateTime
from pydantic import BaseModel
-from typing import Any, List, Optional, Tuple
+from typing import Any, List, Optional
from .. import ChannelConverter
from .. import pydantic_utcdatetime
from .Metric import Metric
+from ..residual.Reading import Reading, get_absolutes_xyz, get_ordinates
class AdjustedMatrix(BaseModel):
@@ -39,26 +40,22 @@ class AdjustedMatrix(BaseModel):
adjusted = adjusted[0 : len(outchannels)]
return adjusted
- # TODO: MOVE GET ABSOLUTES/ORDINATES INTO READING
- def get_metrics(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- ) -> List[Metric]:
+ def get_metrics(self, readings: List[Reading]) -> List[Metric]:
"""Computes mean absolute error and standard deviation for X, Y, Z, and dF between expected and predicted values.
Attributes
----------
- absolutes: X, Y and Z absolutes
- ordinates: H, E and Z ordinates
+ readings: list of Readings
matrix: composed matrix
Outputs
-------
metrics: list of Metric objects
"""
- ordinates = np.vstack((ordinates, np.ones_like(ordinates[0])))
- predicted = self.matrix @ ordinates
+ absolutes = get_absolutes_xyz(readings=readings)
+ ordinates = get_ordinates(readings=readings)
+ stacked_ordinates = np.vstack((ordinates, np.ones_like(ordinates[0])))
+ predicted = self.matrix @ stacked_ordinates
metrics = []
elements = ["X", "Y", "Z", "dF"]
expected = absolutes + tuple(
diff --git a/geomagio/adjusted/Affine.py b/geomagio/adjusted/Affine.py
index fc3252deaf8e6e6fbb0cf21c2ea2411d59957db1..20100eff92b6fe9b99d1538b8483614f95fc9dd8 100644
--- a/geomagio/adjusted/Affine.py
+++ b/geomagio/adjusted/Affine.py
@@ -5,7 +5,12 @@ from pydantic import BaseModel
from typing import List, Optional, Tuple
from .. import pydantic_utcdatetime
-from ..residual import Reading
+from ..residual.Reading import (
+ Reading,
+ get_absolutes_xyz,
+ get_baselines,
+ get_ordinates,
+)
from .AdjustedMatrix import AdjustedMatrix
from .transform import RotationTranslationXY, TranslateOrigins, Transform
@@ -95,7 +100,7 @@ class Affine(BaseModel):
"""
absolutes = get_absolutes_xyz(readings)
baselines = get_baselines(readings)
- ordinates = get_ordinates(readings, self.observatory)
+ ordinates = get_ordinates(readings)
times = get_times(readings)
Ms = []
weights = []
@@ -131,7 +136,7 @@ class Affine(BaseModel):
matrix=M_composed,
pier_correction=pier_correction,
)
- matrix.metrics = matrix.get_metrics(absolutes=absolutes, ordinates=ordinates)
+ matrix.metrics = matrix.get_metrics(readings=readings)
return AdjustedMatrix(
matrix=M_composed,
pier_correction=pier_correction,
@@ -214,39 +219,6 @@ def filter_iqrs(
return weights * good
-def get_absolutes(
- readings: List[Reading],
-) -> Tuple[List[float], List[float], List[float]]:
- """Get H, D and Z absolutes"""
- h_abs = np.array([reading.get_absolute("H").absolute for reading in readings])
- d_abs = np.array([reading.get_absolute("D").absolute for reading in readings])
- z_abs = np.array([reading.get_absolute("Z").absolute for reading in readings])
-
- return (h_abs, d_abs, z_abs)
-
-
-def get_absolutes_xyz(
- readings: List[Reading],
-) -> Tuple[List[float], List[float], List[float]]:
- """Get X, Y and Z absolutes from H, D and Z baselines"""
- h_abs, d_abs, z_abs = get_absolutes(readings)
- # convert from cylindrical to Cartesian coordinates
- x_a = h_abs * np.cos(np.radians(d_abs))
- y_a = h_abs * np.sin(np.radians(d_abs))
- z_a = z_abs
- return (x_a, y_a, z_a)
-
-
-def get_baselines(
- readings: List[Reading],
-) -> Tuple[List[float], List[float], List[float]]:
- """Get H, D and Z baselines"""
- h_bas = np.array([reading.get_absolute("H").baseline for reading in readings])
- d_bas = np.array([reading.get_absolute("D").baseline for reading in readings])
- z_bas = np.array([reading.get_absolute("Z").baseline for reading in readings])
- return (h_bas, d_bas, z_bas)
-
-
def get_epochs(
epochs: List[float],
time: UTCDateTime,
@@ -277,31 +249,6 @@ def get_epochs(
return epoch_start, epoch_end
-def get_ordinates(
- readings: List[Reading], observatory: str
-) -> Tuple[List[float], List[float], List[float]]:
- """Calculates ordinates from absolutes and baselines"""
- h_abs, d_abs, z_abs = get_absolutes(readings)
- h_bas, d_bas, z_bas = get_baselines(readings)
-
- # recreate ordinate variometer measurements from absolutes and baselines
- h_ord = h_abs - h_bas
- d_ord = d_abs - d_bas
- z_ord = z_abs - z_bas
-
- # WebAbsolutes defines/generates h differently than USGS residual
- # method spreadsheets. The following should ensure that ordinate
- # values are converted back to their original raw measurements:
- # TODO: REMOVE OR FIX IN RESIDUALS
- e_o = h_abs * d_ord * 60 / 3437.7468
- if observatory in ["DED", "CMO"]:
- h_o = np.sqrt(h_ord ** 2 - e_o ** 2)
- else:
- h_o = h_ord
- z_o = z_ord
- return (h_o, e_o, z_o)
-
-
def get_times(readings: List[UTCDateTime]):
return np.array([reading.get_absolute("H").endtime for reading in readings])
diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py
index c6ede9ebf28862ba409e7461964400dd9f5211ea..20a423bf6660f4507ec11ab380401470e18b54c7 100644
--- a/geomagio/residual/Reading.py
+++ b/geomagio/residual/Reading.py
@@ -1,6 +1,7 @@
-from typing import Dict, List, Optional
+from typing import Dict, List, Optional, Tuple
from typing_extensions import Literal
+import numpy as np
from obspy import Stream, UTCDateTime
from pydantic import BaseModel
@@ -40,6 +41,15 @@ class Reading(BaseModel):
"""
return [m for m in self.measurements if m.measurement_type == measurement_type]
+ def get_absolute(
+ self,
+ element: str,
+ ) -> Optional[Absolute]:
+ for absolute in self.absolutes:
+ if absolute.element == element:
+ return absolute
+ return None
+
def load_ordinates(
self,
observatory: str,
@@ -119,11 +129,60 @@ class Reading(BaseModel):
):
return True
- def get_absolute(
- self,
- element: str,
- ) -> Optional[Absolute]:
- for absolute in self.absolutes:
- if absolute.element == element:
- return absolute
- return None
+
+def get_absolutes(
+ readings: List[Reading],
+) -> Tuple[List[float], List[float], List[float]]:
+ """Get H, D and Z absolutes"""
+ h_abs = np.array([reading.get_absolute("H").absolute for reading in readings])
+ d_abs = np.array([reading.get_absolute("D").absolute for reading in readings])
+ z_abs = np.array([reading.get_absolute("Z").absolute for reading in readings])
+
+ return (h_abs, d_abs, z_abs)
+
+
+def get_absolutes_xyz(
+ readings: List[Reading],
+) -> Tuple[List[float], List[float], List[float]]:
+ """Get X, Y and Z absolutes from H, D and Z baselines"""
+ h_abs, d_abs, z_abs = get_absolutes(readings)
+ # convert from cylindrical to Cartesian coordinates
+ x_a = h_abs * np.cos(np.radians(d_abs))
+ y_a = h_abs * np.sin(np.radians(d_abs))
+ z_a = z_abs
+ return (x_a, y_a, z_a)
+
+
+def get_baselines(
+ readings: List[Reading],
+) -> Tuple[List[float], List[float], List[float]]:
+ """Get H, D and Z baselines"""
+ h_bas = np.array([reading.get_absolute("H").baseline for reading in readings])
+ d_bas = np.array([reading.get_absolute("D").baseline for reading in readings])
+ z_bas = np.array([reading.get_absolute("Z").baseline for reading in readings])
+ return (h_bas, d_bas, z_bas)
+
+
+def get_ordinates(
+ readings: List[Reading],
+) -> Tuple[List[float], List[float], List[float]]:
+ """Calculates ordinates from absolutes and baselines"""
+ h_abs, d_abs, z_abs = get_absolutes(readings)
+ h_bas, d_bas, z_bas = get_baselines(readings)
+
+ # recreate ordinate variometer measurements from absolutes and baselines
+ h_ord = h_abs - h_bas
+ d_ord = d_abs - d_bas
+ z_ord = z_abs - z_bas
+
+ # WebAbsolutes defines/generates h differently than USGS residual
+ # method spreadsheets. The following should ensure that ordinate
+ # values are converted back to their original raw measurements:
+ # TODO: REMOVE OR FIX IN RESIDUALS
+ e_o = h_abs * d_ord * 60 / 3437.7468
+ if len(readings) > 0 and readings[0].metadata["station"] in ["DED", "CMO"]:
+ h_o = np.sqrt(h_ord ** 2 - e_o ** 2)
+ else:
+ h_o = h_ord
+ z_o = z_ord
+ return (h_o, e_o, z_o)
diff --git a/geomagio/residual/SpreadsheetSummaryFactory.py b/geomagio/residual/SpreadsheetSummaryFactory.py
index 9d1e1f718a8444e4872ee9a2c608ea0b3564d0ab..7a1f5082c11e88c3fd26f5637eb3251e25f049e9 100644
--- a/geomagio/residual/SpreadsheetSummaryFactory.py
+++ b/geomagio/residual/SpreadsheetSummaryFactory.py
@@ -61,7 +61,7 @@ class SpreadsheetSummaryFactory(object):
date = sheet["I1"].value
date = f"{date.year}{date.month:02}{date.day:02}"
return {
- "observatory": sheet["D49"].value[0:3],
+ "station": sheet["D49"].value[0:3],
"pier_correction": sheet["C5"].value,
"instrument": sheet["B3"].value,
"date": date,