From ee113b352070d1ea763d04c4530616d81be026eb Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Thu, 4 Feb 2021 10:59:28 -0700
Subject: [PATCH] Support legacy statefiles, add matrix process method
---
etc/adjusted/adjbou_state_HE_.json | 16 +++++----
geomagio/adjusted/AdjustedMatrix.py | 12 +++++++
geomagio/algorithm/AdjustedAlgorithm.py | 45 +++++++++++++------------
3 files changed, 45 insertions(+), 28 deletions(-)
diff --git a/etc/adjusted/adjbou_state_HE_.json b/etc/adjusted/adjbou_state_HE_.json
index f9c0e3d9d..4aa26dd35 100644
--- a/etc/adjusted/adjbou_state_HE_.json
+++ b/etc/adjusted/adjbou_state_HE_.json
@@ -1,8 +1,12 @@
{
- "matrix": [
- [-1, 0, 0],
- [0, -1, 0],
- [0, 0, 1]
- ],
- "pier_correction": -22.0
+ "PC": -22,
+ "M11": -1,
+ "M12": 0,
+ "M13": 0,
+ "M21": 0,
+ "M22": -1,
+ "M23": 0,
+ "M31": 0,
+ "M32": 0,
+ "M33": 1
}
\ No newline at end of file
diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py
index b2647ff91..ff38990dd 100644
--- a/geomagio/adjusted/AdjustedMatrix.py
+++ b/geomagio/adjusted/AdjustedMatrix.py
@@ -1,3 +1,4 @@
+import numpy as np
from obspy import UTCDateTime
from pydantic import BaseModel
from typing import Any, List, Optional
@@ -24,3 +25,14 @@ class AdjustedMatrix(BaseModel):
metrics: Optional[List[Metric]] = None
starttime: Optional[UTCDateTime] = None
endtime: Optional[UTCDateTime] = None
+
+ def process(self, values: List[List[float]], outchannels=["X", "Y", "Z", "F"]):
+ """ Apply matrix to raw data. Apply pier correction to F when necessary """
+ data = np.vstack([values[0:3]] + [np.ones_like(values[0])])
+ adjusted = self.matrix @ data
+ if "F" in outchannels:
+ f = values[-1] + self.pier_correction
+ adjusted = np.vstack([adjusted[0 : len(outchannels) - 1]] + [f])
+ else:
+ adjusted = adjusted[0 : len(outchannels)]
+ return adjusted
diff --git a/geomagio/algorithm/AdjustedAlgorithm.py b/geomagio/algorithm/AdjustedAlgorithm.py
index adc4b1521..78a38ae25 100644
--- a/geomagio/algorithm/AdjustedAlgorithm.py
+++ b/geomagio/algorithm/AdjustedAlgorithm.py
@@ -43,24 +43,33 @@ class AdjustedAlgorithm(Algorithm):
"""Load algorithm state from a file.
File name is self.statefile.
"""
- pier_correction = 0
if self.statefile is None:
return
- # Adjusted matrix defaults to identity matrix
- matrix_size = len([c for c in self.get_input_channels() if c != "F"]) + 1
- matrix = np.eye(matrix_size)
- data = None
try:
with open(self.statefile, "r") as f:
data = f.read()
data = json.loads(data)
except IOError as err:
- sys.stderr.write("I/O error {0}".format(err))
- if data is None or data == "":
- return
- self.matrix = AdjustedMatrix(
- matrix=np.array(data["matrix"]), pier_correction=data["pier_correction"]
- )
+ raise FileNotFoundError("statefile not found")
+ pier_correction = 0
+ # Adjusted matrix defaults to identity matrix
+ matrix_size = len([c for c in self.get_input_channels() if c != "F"]) + 1
+ matrix = np.eye(matrix_size)
+ if self.statefile is not None:
+ try:
+ if "pier_correction" in data.keys():
+ # read data
+ matrix = np.array(data["matrix"])
+ pier_correction = data["pier_correction"]
+ if "PC" in data.keys():
+ # read data from legacy format
+ for row in range(matrix_size):
+ for col in range(matrix_size):
+ matrix[row, col] = np.float64(data[f"M{row+1}{col+1}"])
+ pier_correction = np.float64(data["PC"])
+ except:
+ raise KeyError("matrix and pier correction not found in statefile")
+ self.matrix = AdjustedMatrix(matrix=matrix, pier_correction=pier_correction)
def save_state(self):
"""Save algorithm state to a file.
@@ -123,14 +132,9 @@ class AdjustedAlgorithm(Algorithm):
inchannels = self.get_input_channels()
outchannels = self.get_output_channels()
raws = np.vstack(
- [
- stream.select(channel=channel)[0].data
- for channel in inchannels
- if channel != "F"
- ]
- + [np.ones_like(stream[0].data)]
+ [stream.select(channel=channel)[0].data for channel in inchannels]
)
- adjusted = np.matmul(self.matrix.matrix, raws)
+ adjusted = self.matrix.process(values=raws, outchannels=outchannels)
out = Stream(
[
self.create_trace(
@@ -138,12 +142,9 @@ class AdjustedAlgorithm(Algorithm):
stream.select(channel=inchannels[i])[0].stats,
adjusted[i],
)
- for i in range(len(adjusted) - 1)
+ for i in range(len(adjusted))
]
)
- if "F" in inchannels and "F" in outchannels:
- f = stream.select(channel="F")[0]
- out += self.create_trace("F", f.stats, f.data + self.matrix.pier_correction)
return out
def can_produce_data(self, starttime, endtime, stream):
--
GitLab