From 720c5490f3e0464b3742350b0be5c560ff70367d Mon Sep 17 00:00:00 2001
From: Alex Wernle <awernle@usgs.gov>
Date: Tue, 22 Aug 2023 09:40:40 -0600
Subject: [PATCH] New BIQParserFactory to parse Geomag BIQ QD files.
---
geomagio/biq/BIQParserFactory.py | 129 +++++++++++++++++++++++++++++++
1 file changed, 129 insertions(+)
create mode 100644 geomagio/biq/BIQParserFactory.py
diff --git a/geomagio/biq/BIQParserFactory.py b/geomagio/biq/BIQParserFactory.py
new file mode 100644
index 00000000..59d0117a
--- /dev/null
+++ b/geomagio/biq/BIQParserFactory.py
@@ -0,0 +1,129 @@
+"""Factory that parses BIQ files into arrays and/or ObsPy Stream."""
+import numpy as np
+import obspy
+
+from typing import BinaryIO
+
+
+class BIQParser:
+ """Geomag BIQ file fields."""
+
+ def __init__(self, filename: BinaryIO):
+ self.filename = filename
+
+ # define the data types for a day of data
+ self.day_block_dtype = np.dtype(
+ [
+ ("Station ID", "S4"),
+ ("Year & Day", np.int32),
+ ("Co-Latitude", np.int32),
+ ("Longitude", np.int32),
+ ("Elevation", np.int32),
+ ("Orientation", "S4"),
+ ("Origination", "S4"),
+ ("D Conversion", np.int32),
+ ("Quality", "S4"),
+ ("Instrumentation", "S4"),
+ ("Custom Field 1", np.int32),
+ ("Custom Field 2", np.int32),
+ ("Custom Field 3", "S4"),
+ ("Custom Field 4", "S4"),
+ ("Custom Field 5", np.int32),
+ ("Custom Field 6", "S4"),
+ ("X-minute", np.int32, (1440,)),
+ ("Y-minute", np.int32, (1440,)),
+ ("Z-minute", np.int32, (1440,)),
+ ("G-minute", np.int32, (1440,)),
+ ("X-hourly", np.int32, (24,)),
+ ("Y-hourly", np.int32, (24,)),
+ ("Z-hourly", np.int32, (24,)),
+ ("G-hourly", np.int32, (24,)),
+ ("X-daily", np.int32, (1,)),
+ ("Y-daily", np.int32, (1,)),
+ ("Z-daily", np.int32, (1,)),
+ ("G-daily", np.int32, (1,)),
+ ("K-values", np.int32, (8,)),
+ ("Leftovers", np.int32, (4,)),
+ ]
+ )
+
+ def read_data(self):
+ with open(self.filename, "rb") as f:
+ data = np.fromfile(f, dtype=self.day_block_dtype)
+
+ # call the check_station_id method to check that the Station ID for each day is the same
+ self.check_station_id(data)
+
+ return data
+
+ def check_station_id(self, data):
+ # get the unique station IDs
+ station_ids = set(day_data["Station ID"].decode() for day_data in data)
+
+ # raise an error if there is more than one unique station ID
+ if len(station_ids) > 1:
+ raise ValueError(
+ "One or more station IDs do not match. File may be corrupted or written incorrectly!"
+ )
+
+ def get_stream(self):
+ data = self.read_data()
+
+ # combine the minute data from all the days into an obspy stream
+ st = obspy.Stream()
+ tr = obspy.Trace(
+ header={
+ "network": "NT",
+ "location": "Q0", # THIS IS AN ASSUMPTION!
+ "station": data[0]["Station ID"].decode().strip(),
+ "delta": 60,
+ "starttime": obspy.UTCDateTime(str(data[0]["Year & Day"])),
+ }
+ )
+
+ for day in data:
+ # add stats for header information
+ tr.stats.starttime = obspy.UTCDateTime(str(day["Year & Day"]))
+
+ x_tr = tr.copy()
+ y_tr = tr.copy()
+ z_tr = tr.copy()
+ g_tr = tr.copy()
+ x_tr.stats.channel = "X"
+ y_tr.stats.channel = "Y"
+ z_tr.stats.channel = "Z"
+ g_tr.stats.channel = "G"
+ # convert to floats so that we're not stuck dealing with ints later
+ x_tr.data = day["X-minute"].astype(np.float64)
+ y_tr.data = day["Y-minute"].astype(np.float64)
+ z_tr.data = day["Z-minute"].astype(np.float64)
+ g_tr.data = day["G-minute"].astype(np.float64)
+
+ st += obspy.Stream(traces=[x_tr, y_tr, z_tr, g_tr])
+ st.merge()
+
+ # data was stored as scaled integers, unscale them and translate the 99999's to nan
+ for tr in st:
+ tr.data[tr.data == 999999] = np.nan
+ tr.data = tr.data / 10
+
+ st = self.append_f(st)
+
+ return st
+
+ def append_f(self, st):
+ # calculate F from XYZ and G ("delta-F")
+ f_xyz = np.zeros_like(st.select(channel="G")[0].data)
+ for tr in st.select(channel="[X|Y|Z]"):
+ f_xyz = f_xyz + np.square(tr.data)
+ f_xyz = np.sqrt(f_xyz)
+
+ f_tr = obspy.Trace(
+ f_xyz - st.select(channel="G")[0].data,
+ header=st.select(channel="G")[0].stats,
+ )
+ f_tr.stats.channel = "F"
+
+ st.append(f_tr)
+
+ return st
--
GitLab