diff --git a/etc/adjusted/BOU201601vmin.min b/etc/adjusted/BOU201601vmin.min
index 57c20c44233ae76dfd9dae307ce91c14086b4651..83ef32f4830f78ffc089e167dc6230832db4e2ee 100644
--- a/etc/adjusted/BOU201601vmin.min
+++ b/etc/adjusted/BOU201601vmin.min
@@ -5,7 +5,7 @@
Geodetic Latitude 40.137 |
Geodetic Longitude 254.764 |
Elevation 1682 |
- Reported HEZF |
+ Reported UVWF |
Sensor Orientation HDZF |
Digital Sampling 100.0 second |
Data Interval Type filtered 1-minute (00:15-01:45) |
@@ -19,7 +19,7 @@
# CONDITIONS OF USE: The Conditions of Use for data provided |
# through INTERMAGNET and acknowledgement templates can be found at |
# www.intermagnet.org |
-DATE TIME DOY BOUH BOUE BOUZ BOUF |
+DATE TIME DOY BOUU BOUV BOUW BOUF |
2016-01-01 00:00:00.000 001 20735.93 -99.77 47370.21 52248.63
2016-01-01 00:01:00.000 001 20735.48 -95.30 47370.53 52248.72
2016-01-01 00:02:00.000 001 20736.30 -91.25 47371.05 52249.49
diff --git a/etc/adjusted/BOU202005adj.min b/etc/adjusted/BOU202005adj.min
index 55f8fcc64580668ad815ae8706c1628cbbea0285..551fd3ed7ec4fabb3ecf4b5a360b782ce4f45b3c 100644
--- a/etc/adjusted/BOU202005adj.min
+++ b/etc/adjusted/BOU202005adj.min
@@ -5,7 +5,7 @@
Geodetic Latitude 40.137 |
Geodetic Longitude 254.763 |
Elevation 1682 |
- Reported HENULNUL |
+ Reported UVNULNUL |
Sensor Orientation HDZF |
Digital Sampling 0.01 second |
Data Interval Type filtered 1-minute (00:15-01:45) |
@@ -17,7 +17,7 @@
# CONDITIONS OF USE: The Conditions of Use for data provided |
# through INTERMAGNET and acknowledgement templates can be found at |
# www.intermagnet.org |
-DATE TIME DOY BOUH BOUE BOUNUL BOUNUL |
+DATE TIME DOY BOUU BOUV BOUNUL BOUNUL |
2020-05-20 16:15:00.000 141 -20799.04 80.92 99999.00 99999.00
2020-05-20 16:16:00.000 141 -20798.26 79.97 99999.00 99999.00
2020-05-20 16:17:00.000 141 -20799.43 81.22 99999.00 99999.00
diff --git a/etc/adjusted/BOU202005vmin.min b/etc/adjusted/BOU202005vmin.min
index 545671c9b844f2f36fa64467d195b52c2780c762..f3d0f15bc47bdb0d77d374942c155d20e95955ee 100644
--- a/etc/adjusted/BOU202005vmin.min
+++ b/etc/adjusted/BOU202005vmin.min
@@ -5,7 +5,7 @@
Geodetic Latitude 40.137 |
Geodetic Longitude 254.763 |
Elevation 1682 |
- Reported HENULNUL |
+ Reported UVNULNUL |
Sensor Orientation HDZF |
Digital Sampling 0.01 second |
Data Interval Type filtered 1-minute (00:15-01:45) |
@@ -19,7 +19,7 @@
# CONDITIONS OF USE: The Conditions of Use for data provided |
# through INTERMAGNET and acknowledgement templates can be found at |
# www.intermagnet.org |
-DATE TIME DOY BOUH BOUE BOUNUL BOUNUL |
+DATE TIME DOY BOUU BOUV BOUNUL BOUNUL |
2020-05-20 16:15:00.000 141 20799.04 -80.92 99999.00 99999.00
2020-05-20 16:16:00.000 141 20798.26 -79.97 99999.00 99999.00
2020-05-20 16:17:00.000 141 20799.43 -81.22 99999.00 99999.00
diff --git a/geomagio/ChannelConverter.py b/geomagio/ChannelConverter.py
index a7c3bfc8b255f77521bce41916ab0b7484a02693..3884af3d286050fba482a6eee352432d87c94c9f 100644
--- a/geomagio/ChannelConverter.py
+++ b/geomagio/ChannelConverter.py
@@ -5,7 +5,7 @@ the geomagnetic community.
We use three coordinate systems.
Geo: Based on Geographic North. X, Y, Z
X is north, Y is east
-Obs: Based on the observatories orientaion. H, E, Z [d]
+Obs: Based on the observatories orientaion. U, V, W [d]
Mag: Based on Magnetic North. H, D, Z [E]
d0: Declination baseline in radians
@@ -29,15 +29,15 @@ R2M = 180.0 / numpy.pi * 60 # Radians to Minutes
# ###
-def get_geo_from_obs(h, e, d0=0):
+def get_geo_from_obs(u, v, d0=0):
"""gets the geographical components given the observatory components.
Parameters
__________
- h: array_like
- the h component from the observatory
- e: array_like
- the e component from the observatory
+ u: array_like
+ the u component from the observatory
+ v: array_like
+ the v component from the observatory
d0: float
the declination baseline angle in radians
@@ -47,7 +47,7 @@ def get_geo_from_obs(h, e, d0=0):
[0]: x component as a float
[1]: y component as a float
"""
- mag_h, mag_d = get_mag_from_obs(h, e, d0)
+ mag_h, mag_d = get_mag_from_obs(u, v, d0)
return get_geo_from_mag(mag_h, mag_d)
@@ -112,15 +112,15 @@ def get_geo_y_from_mag(h, d):
# ###
# get magnetic north coordinates from....
# ###
-def get_mag_from_obs(h, e, d0=0):
+def get_mag_from_obs(u, v, d0=0):
"""gets the magnetic north components given the observatory components.
Parameters
__________
- h: array_like
- the h component from the observatory
- e: array_like
- the e component from the observatory
+ u: array_like
+ the u component from the observatory
+ v: array_like
+ the v component from the observatory
d0: float
the declination baseline angle in radians
@@ -130,8 +130,8 @@ def get_mag_from_obs(h, e, d0=0):
[0]: total h component as a float
[1]: total d declination as a float
"""
- mag_h = get_mag_h_from_obs(h, e)
- mag_d = get_mag_d_from_obs(h, e, d0)
+ mag_h = get_mag_h_from_obs(u, v)
+ mag_d = get_mag_d_from_obs(u, v, d0)
return (mag_h, mag_d)
@@ -156,15 +156,15 @@ def get_mag_from_geo(x, y):
return (mag_h, mag_d)
-def get_mag_d_from_obs(h, e, d0=0):
+def get_mag_d_from_obs(u, v, d0=0):
"""gets the magnetic d component given the observatory components.
Parameters
__________
- h: array_like
- the h component from the observatory
- e: array_like
- the e component from the observatory
+ u: array_like
+ the u component from the observatory
+ v: array_like
+ the v component from the observatory
d0: float
the declination baseline angle in radians
@@ -173,7 +173,7 @@ def get_mag_d_from_obs(h, e, d0=0):
array_like
the total magnetic declination
"""
- return numpy.add(d0, get_obs_d_from_obs(h, e))
+ return numpy.add(d0, get_obs_d_from_obs(u, v))
def get_mag_d_from_geo(x, y):
@@ -194,22 +194,22 @@ def get_mag_d_from_geo(x, y):
return numpy.arctan2(y, x)
-def get_mag_h_from_obs(h, e):
+def get_mag_h_from_obs(u, v):
"""gets the magnetic h component given the observatory components.
Parameters
__________
- h: array_like
- the h component from the observatory
- e: array_like
- the e component from the observatory
+ u: array_like
+ the u component from the observatory
+ v: array_like
+ the v component from the observatory
Returns
_______
array_like
the total magnetic h component
"""
- return numpy.hypot(h, e)
+ return numpy.hypot(u, v)
def get_mag_h_from_geo(x, y):
@@ -271,32 +271,32 @@ def get_obs_from_mag(h, d, d0=0):
Returns
_______
tuple of array_like
- [0]: observatory h component
- [1]: observatory e component
+ [0]: observatory u component
+ [1]: observatory v component
[2]: observatory d declination
"""
- obs_h = get_obs_h_from_mag(h, d, d0)
- obs_e = get_obs_e_from_mag(h, d, d0)
- return (obs_h, obs_e)
+ obs_u = get_obs_u_from_mag(h, d, d0)
+ obs_v = get_obs_v_from_mag(h, d, d0)
+ return (obs_u, obs_v)
# inividual get obs from calls
-def get_obs_d_from_obs(h, e):
+def get_obs_d_from_obs(u, v):
"""gets the observatory d declination given the observatory components.
Parameters
__________
- h: array_like
- the h component from the observatory
- e: array_like
- the e component from the observatory
+ u: array_like
+ the u component from the observatory
+ v: array_like
+ the v component from the observatory
Returns
_______
array_like
the observatory d declination
"""
- return numpy.arctan2(e, h)
+ return numpy.arctan2(v, u)
def get_obs_d_from_mag_d(d, d0=0):
@@ -318,8 +318,8 @@ def get_obs_d_from_mag_d(d, d0=0):
return numpy.subtract(d, d0)
-def get_obs_e_from_mag(h, d, d0=0):
- """gets the observatory e component given the magnetic components.
+def get_obs_v_from_mag(h, d, d0=0):
+ """gets the observatory v component given the magnetic components.
Parameters
__________
@@ -333,32 +333,32 @@ def get_obs_e_from_mag(h, d, d0=0):
Returns
_______
array_like
- the observatory e component
+ the observatory v component
"""
obs_d = get_obs_d_from_mag_d(d, d0)
return numpy.multiply(h, numpy.sin(obs_d))
-def get_obs_e_from_obs(h, d):
- """gets the observatory e component given the observatory components.
+def get_obs_v_from_obs(h, d):
+ """gets the observatory v component given the observatory components.
Parameters
__________
- h: array_like
- the observatory h component.
+ u: array_like
+ the observatory u component.
d: array_like
the observatory d declination.
Returns
_______
array_like
- the observatory e component
+ the observatory v component
"""
return numpy.multiply(h, numpy.tan(d))
-def get_obs_h_from_mag(h, d, d0=0):
- """gets the observatory h component given the magnetic north components
+def get_obs_u_from_mag(h, d, d0=0):
+ """gets the observatory u component given the magnetic north components
Parameters
__________
@@ -372,7 +372,7 @@ def get_obs_h_from_mag(h, d, d0=0):
Returns
_______
array_like
- the observatory h component
+ the observatory u component
"""
obs_d = get_obs_d_from_mag_d(d, d0)
return numpy.multiply(h, numpy.cos(obs_d))
@@ -405,7 +405,7 @@ def get_computed_f_using_squares(x, y, z):
Notes
-----
This works for geographic coordinates, or observatory coordinates.
- ie x, y, z or h, e, z
+ ie x, y, z or u, v, w
We're using variables x,y,z to represent generic cartisian coordinates.
"""
x2 = numpy.square(x)
diff --git a/geomagio/StreamConverter.py b/geomagio/StreamConverter.py
index 4ef3232c5db1ef5b0ec7cf4eb62357a9d5550413..4677acee27626dc968216793895ba4102668b302 100644
--- a/geomagio/StreamConverter.py
+++ b/geomagio/StreamConverter.py
@@ -3,7 +3,7 @@
We use three coordinate systems.
Geo: Based on Geographic North. X, Y, Z, F
X is north, Y is east, Z is down
-Obs: Based on the observatories orientaion. H, E, Z, F, d0
+Obs: Based on the observatories orientaion. U, V, W, F, d0
Mag: Based on Magnetic North. H, D, Z, F
"""
from __future__ import absolute_import
@@ -48,7 +48,7 @@ def get_geo_from_obs(obs):
Parameters
----------
stream : obspy.core.Stream
- stream containing observatory components H, D or E, Z, and F.
+ stream containing observatory components U, D or V, W, and F.
Returns
-------
@@ -86,20 +86,20 @@ def get_deltaf_from_obs(obs):
Parameters
----------
obs: obspy.core.Stream
- stream containing the observatory components H, D or E, Z, and F.
+ stream containing the observatory components U, D or V, W, and F.
Returns
-------
obspy.core.Stream
stream object containing delta f values
"""
- h = obs.select(channel="H")[0]
- z = obs.select(channel="Z")[0]
+ u = obs.select(channel="U")[0]
+ w = obs.select(channel="W")[0]
fs = obs.select(channel="F")[0]
- e = __get_obs_e_from_obs(obs)
- fv = ChannelConverter.get_computed_f_using_squares(h, e, z)
+ v = __get_obs_v_from_obs(obs)
+ fv = ChannelConverter.get_computed_f_using_squares(u, v, w)
G = ChannelConverter.get_deltaf(fv, fs)
- return obspy.core.Stream((__get_trace("G", h.stats, G),))
+ return obspy.core.Stream((__get_trace("G", u.stats, G),))
def get_mag_from_geo(geo):
@@ -137,28 +137,28 @@ def get_mag_from_obs(obs):
Parameters
----------
obs : obspy.core.Stream
- stream containing observatory components H, D or E, Z, and F.
+ stream containing observatory components U, D or V, W, and F.
Returns
-------
obspy.core.Stream
new stream object containing magnetic components H, D, Z, and F.
"""
- h = obs.select(channel="H")[0]
- e = __get_obs_e_from_obs(obs)
- z = obs.select(channel="Z")
+ u = obs.select(channel="U")[0]
+ v = __get_obs_v_from_obs(obs)
+ w = obs.select(channel="W")
f = obs.select(channel="F")
- obs_h = h.data
- obs_e = e.data
+ obs_u = u.data
+ obs_v = v.data
d0 = ChannelConverter.get_radians_from_minutes(
- numpy.float64(e.stats.declination_base) / 10
+ numpy.float64(v.stats.declination_base) / 10
)
- (mag_h, mag_d) = ChannelConverter.get_mag_from_obs(obs_h, obs_e, d0)
+ (mag_h, mag_d) = ChannelConverter.get_mag_from_obs(obs_u, obs_v, d0)
return (
obspy.core.Stream(
- (__get_trace("H", h.stats, mag_h), __get_trace("D", e.stats, mag_d))
+ (__get_trace("H", u.stats, mag_h), __get_trace("D", v.stats, mag_d))
)
- + z
+ + w
+ f
)
@@ -193,7 +193,7 @@ def get_obs_from_mag(mag, include_d=False):
Returns
-------
obspy.core.Stream
- new stream object containing observatory components H, D, E, Z, and F
+ new stream object containing observatory components U, D, V, W, and F
"""
h = mag.select(channel="H")[0]
d = mag.select(channel="D")[0]
@@ -205,43 +205,43 @@ def get_obs_from_mag(mag, include_d=False):
d0 = ChannelConverter.get_radians_from_minutes(
numpy.float64(d.stats.declination_base) / 10
)
- (obs_h, obs_e) = ChannelConverter.get_obs_from_mag(mag_h, mag_d, d0)
+ (obs_u, obs_v) = ChannelConverter.get_obs_from_mag(mag_h, mag_d, d0)
- traces = (__get_trace("H", h.stats, obs_h), __get_trace("E", d.stats, obs_e))
+ traces = (__get_trace("U", h.stats, obs_u), __get_trace("V", d.stats, obs_v))
if include_d:
- obs_d = ChannelConverter.get_obs_d_from_obs(obs_h, obs_e)
+ obs_d = ChannelConverter.get_obs_d_from_obs(obs_u, obs_v)
traces = traces + (__get_trace("D", d.stats, obs_d),)
return obspy.core.Stream(traces) + z + f
-def get_obs_from_obs(obs, include_e=False, include_d=False):
+def get_obs_from_obs(obs, include_v=False, include_d=False):
"""Fill in the observatory parameters as requested
Parameters
----------
stream: obspy.core.Stream
- stream containing the observatory components H, D or E, Z, and F.
+ stream containing the observatory components U, D or V, W, and F.
include_e: boolean
- whether to include the e component
+ whether to include the v component
include_d: boolean
whether to include the d component
Returns
-------
obspy.core.Stream
- new stream object containing observatory components H, D, E, Z, and F
+ new stream object containing observatory components U, D, V, W, and F
"""
- h = obs.select(channel="H")[0]
- z = obs.select(channel="Z")
+ u = obs.select(channel="U")[0]
+ w = obs.select(channel="W")
f = obs.select(channel="F")
- traces = (h,)
+ traces = (u,)
if include_d:
d = __get_obs_d_from_obs(obs)
traces = traces + (d,)
- if include_e:
- e = __get_obs_e_from_obs(obs)
- traces = traces + (e,)
- return obspy.core.Stream(traces) + z + f
+ if include_v:
+ v = __get_obs_v_from_obs(obs)
+ traces = traces + (v,)
+ return obspy.core.Stream(traces) + w + f
def __get_trace(channel, stats, data):
@@ -269,12 +269,12 @@ def __get_trace(channel, stats, data):
def __get_obs_d_from_obs(obs):
"""Get trace containing observatory D component.
- Returns D if found, otherwise computes D from H, E, D0.
+ Returns D if found, otherwise computes D from U, V, D0.
Parameters
----------
obs : obspy.core.Stream
- observatory components (D) or (H, E).
+ observatory components (D) or (U, V).
Returns
-------
@@ -284,35 +284,35 @@ def __get_obs_d_from_obs(obs):
try:
d = obs.select(channel="D")[0]
except IndexError:
- h = obs.select(channel="H")[0]
- e = obs.select(channel="E")[0]
+ u = obs.select(channel="U")[0]
+ v = obs.select(channel="V")[0]
d = __get_trace(
- "D", e.stats, ChannelConverter.get_obs_d_from_obs(h.data, e.data)
+ "D", v.stats, ChannelConverter.get_obs_d_from_obs(u.data, v.data)
)
return d
-def __get_obs_e_from_obs(obs):
- """Get trace containing observatory E component.
+def __get_obs_v_from_obs(obs):
+ """Get trace containing observatory V component.
- Returns E if found, otherwise computes E from H,D.
+ Returns V if found, otherwise computes V from U,D.
Parameters
----------
obs : obspy.core.Stream
- observatory components (E) or (H, D).
+ observatory components (V) or (U, D).
Returns
-------
obspy.core.Trace
- observatory component E.
+ observatory component V.
"""
try:
- e = obs.select(channel="E")[0]
+ v = obs.select(channel="V")[0]
except IndexError:
- h = obs.select(channel="H")[0]
+ u = obs.select(channel="U")[0]
d = obs.select(channel="D")[0]
- e = __get_trace(
- "E", d.stats, ChannelConverter.get_obs_e_from_obs(h.data, d.data)
+ v = __get_trace(
+ "V", d.stats, ChannelConverter.get_obs_v_from_obs(u.data, d.data)
)
- return e
+ return v
diff --git a/geomagio/TimeseriesFactory.py b/geomagio/TimeseriesFactory.py
index 7446f81ee6cc4877dde540044c53c742e9a3564f..c8294c82d67ff4a7992f0b1f1b05c2406354add1 100644
--- a/geomagio/TimeseriesFactory.py
+++ b/geomagio/TimeseriesFactory.py
@@ -1,11 +1,13 @@
"""Abstract Timeseries Factory Interface."""
from __future__ import absolute_import, print_function
+from typing import List
import numpy
-import obspy.core
+from obspy import Stream
import os
import sys
from .TimeseriesFactoryException import TimeseriesFactoryException
+from . import derived
from . import TimeseriesUtility
from . import Util
@@ -59,6 +61,43 @@ class TimeseriesFactory(object):
self.urlTemplate = urlTemplate
self.urlInterval = urlInterval
+ def get_derived_timeseries(self, timeseries: Stream) -> Stream:
+ """compute missing derived channels"""
+ missing_channels = derived.get_missing_channels(timeseries=timeseries)
+ if not missing_channels:
+ return timeseries
+ input_timeseries = self.add_required_channels(
+ timeseries=timeseries, channels=missing_channels
+ )
+ output_channels = TimeseriesUtility.get_channels(stream=timeseries)
+ return derived.get_timeseries(
+ timeseries=input_timeseries, channels=output_channels
+ )
+
+ def add_required_channels(self, timeseries: Stream, channels: List[str]) -> Stream:
+ """adds required channels into timeseries for derived products"""
+ required_channels = derived.get_required_channels(
+ channels=channels, data_type=timeseries[0].stats.data_type
+ )
+ output_stream = Stream()
+ stats = timeseries[0].stats
+ defined_channels = TimeseriesUtility.get_channels(
+ stream=timeseries, contains_data=True
+ )
+ for channel in required_channels:
+ if channel in defined_channels:
+ output_stream += timeseries.select(channel=channel)[0]
+ output_stream += self.get_timeseries(
+ starttime=stats.starttime,
+ endtime=stats.endtime,
+ observatory=stats.station,
+ channels=(channel,),
+ type=stats.data_type,
+ interval=TimeseriesUtility.get_interval_from_delta(stats.delta),
+ derive_missing=False,
+ )
+ return output_stream
+
def get_timeseries(
self,
starttime,
@@ -67,7 +106,8 @@ class TimeseriesFactory(object):
channels=None,
type=None,
interval=None,
- ):
+ derive_missing: bool = True,
+ ) -> Stream:
"""Get timeseries data.
Support for specific channels, types, and intervals varies
@@ -93,11 +133,8 @@ class TimeseriesFactory(object):
interval : {'day', 'hour', 'minute', 'month', 'second'}
data interval, optional.
uses default if unspecified.
-
- Returns
- -------
- obspy.core.Stream
- stream containing traces for requested timeseries.
+ derive_missing
+ if True, calculates trivial derived products missing from initial request
Raises
------
@@ -109,7 +146,7 @@ class TimeseriesFactory(object):
type = type or self.type
interval = interval or self.interval
- timeseries = obspy.core.Stream()
+ timeseries = Stream()
urlIntervals = Util.get_intervals(
starttime=starttime, endtime=endtime, size=self.urlInterval
)
@@ -140,7 +177,7 @@ class TimeseriesFactory(object):
print("Error parsing data: " + str(e), file=sys.stderr)
print(data, file=sys.stderr)
if channels is not None:
- filtered = obspy.core.Stream()
+ filtered = Stream()
for channel in channels:
filtered += timeseries.select(channel=channel)
timeseries = filtered
@@ -152,6 +189,8 @@ class TimeseriesFactory(object):
pad=True,
fill_value=numpy.nan,
)
+ if derive_missing:
+ timeseries = self.get_derived_timeseries(timeseries=timeseries)
return timeseries
def parse_string(self, data, **kwargs):
diff --git a/geomagio/TimeseriesUtility.py b/geomagio/TimeseriesUtility.py
index e0729c21b81796707a06cff1184765aeab5be29c..a2e152b682384b9fc1b71abf82ba64cd25fed442 100644
--- a/geomagio/TimeseriesUtility.py
+++ b/geomagio/TimeseriesUtility.py
@@ -3,6 +3,8 @@ from builtins import range
from datetime import datetime
import math
import sys
+from typing import List
+
import numpy
from obspy.core import Stats, Stream, Trace, UTCDateTime
@@ -288,21 +290,22 @@ def get_merged_gaps(gaps):
return merged_gaps
-def get_channels(stream):
+def get_channels(stream: Stream, contains_data: bool = False) -> List[str]:
"""Get a list of channels in a stream.
-
- Parameters
- ----------
- stream : Stream
-
- Returns
- -------
- channels : array_like
+ If contains_data is True, only channels containing data are returned.
"""
channels = {}
for trace in stream:
- channel = trace.stats.channel
+ stats = trace.stats
+ channel = stats.channel
if channel:
+ if contains_data and not has_all_channels(
+ stream=stream,
+ channels=(channel,),
+ starttime=stats.starttime,
+ endtime=stats.endtime,
+ ):
+ continue
channels[channel] = True
return [ch for ch in channels]
@@ -610,6 +613,9 @@ def split_trace(trace: Trace, size: int = 86400) -> Stream:
size=size,
trim=True,
):
+ if interval["start"] == interval["end"]:
+ stream += out_trace
+ continue
stream += out_trace.slice(
starttime=interval["start"],
endtime=interval["end"] - out_trace.stats.delta,
diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py
index 88e7169bff4200f65eb1899cfa89077afbafab5b..f02671ae2161d2e751b758b13cead236f97e73e1 100644
--- a/geomagio/adjusted/AdjustedMatrix.py
+++ b/geomagio/adjusted/AdjustedMatrix.py
@@ -32,7 +32,7 @@ class AdjustedMatrix(BaseModel):
def process(
self,
stream: Stream,
- inchannels=["H", "E", "Z", "F"],
+ inchannels=["U", "V", "W", "F"],
outchannels=["X", "Y", "Z", "F"],
):
""" Apply matrix to raw data. Apply pier correction to F when necessary """
diff --git a/geomagio/algorithm/AdjustedAlgorithm.py b/geomagio/algorithm/AdjustedAlgorithm.py
index fcbe2dcebb8e66e2d95e9235a7134e905f52eca6..ef92e2ed471e8b7e767f28b0a23a6377961b2065 100644
--- a/geomagio/algorithm/AdjustedAlgorithm.py
+++ b/geomagio/algorithm/AdjustedAlgorithm.py
@@ -23,7 +23,7 @@ class AdjustedAlgorithm(Algorithm):
inchannels=None,
outchannels=None,
):
- inchannels = inchannels or ["H", "E", "Z", "F"]
+ inchannels = inchannels or ["U", "V", "W", "F"]
outchannels = outchannels or ["X", "Y", "Z", "F"]
Algorithm.__init__(
self,
diff --git a/geomagio/algorithm/DeltaFAlgorithm.py b/geomagio/algorithm/DeltaFAlgorithm.py
index 22ebab00988a260a6fc986048ef076caa5b2ad2a..d2e1735f28f1b9a02c959864016f5685c94d0707 100644
--- a/geomagio/algorithm/DeltaFAlgorithm.py
+++ b/geomagio/algorithm/DeltaFAlgorithm.py
@@ -10,11 +10,11 @@ from .. import StreamConverter
# List of channels by geomagnetic observatory orientation.
# geo represents a geographic north/south orientation
# obs represents the sensor orientation aligned close to the mag orientation
-# obsd is the same as obs, but with D(declination) instead of E (e/w vector)
+# obsd is the same as obs, but with D(declination) instead of V (e/w vector)
CHANNELS = {
"geo": ["X", "Y", "Z", "F"],
- "obs": ["H", "E", "Z", "F"],
- "obsd": ["H", "D", "Z", "F"],
+ "obs": ["U", "V", "W", "F"],
+ "obsd": ["U", "D", "W", "F"],
}
diff --git a/geomagio/algorithm/FilterAlgorithm.py b/geomagio/algorithm/FilterAlgorithm.py
index da0f2df18f115f7be5d9e3d25c86147917a65e4a..304992771bd8eaeaa64abfcbb008dbb9d5549d75 100644
--- a/geomagio/algorithm/FilterAlgorithm.py
+++ b/geomagio/algorithm/FilterAlgorithm.py
@@ -411,7 +411,7 @@ class FilterAlgorithm(Algorithm):
steps = self.get_filter_steps()
# calculate start/end from inverted step array
for step in reversed(steps):
- start_interval = get_nearest_time(step=step, output_time=start, left=False)
+ start_interval = get_nearest_time(step=step, output_time=start, left=True)
end_interval = get_nearest_time(step=step, output_time=end, left=True)
start, end = start_interval["data_start"], end_interval["data_end"]
return (start, end)
diff --git a/geomagio/algorithm/XYZAlgorithm.py b/geomagio/algorithm/XYZAlgorithm.py
index 1f43c87853c8da21d7a4636ad1bf7560f35e5471..766545224a6152b2c1c7a3e228560a1b97da9ef3 100644
--- a/geomagio/algorithm/XYZAlgorithm.py
+++ b/geomagio/algorithm/XYZAlgorithm.py
@@ -12,12 +12,12 @@ from .. import StreamConverter
# geo represents a geographic north/south orientation
# mag represents the (calculated)instantaneous mangnetic north orientation
# obs represents the sensor orientation aligned close to the mag orientation
-# obsd is the same as obs, but with D(declination) instead of E (e/w vector)
+# obsd is the same as obs, but with D(declination) instead of V (e/w vector)
CHANNELS = {
"geo": ["X", "Y", "Z", "F"],
"mag": ["H", "D", "Z", "F"],
- "obs": ["H", "E", "Z", "F"],
- "obsd": ["H", "D", "Z", "F"],
+ "obs": ["U", "V", "W", "F"],
+ "obsd": ["U", "D", "W", "F"],
}
diff --git a/geomagio/derived.py b/geomagio/derived.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f0807ecae2fe5d75d8270388d7d19bdb255654a
--- /dev/null
+++ b/geomagio/derived.py
@@ -0,0 +1,98 @@
+from typing import List
+
+from obspy import Stream, Trace
+
+from .algorithm.DeltaFAlgorithm import DeltaFAlgorithm
+from .algorithm.XYZAlgorithm import XYZAlgorithm
+from .TimeseriesUtility import get_channels, has_all_channels
+
+
+def compute_trace(timeseries: Stream, channel: str) -> Trace:
+ """calculates derived product from input timeseries"""
+ data_type = timeseries[0].stats.data_type
+ input_channels = get_derived_input_channels(channel=channel, data_type=data_type)
+ input_timeseries = Stream([timeseries.select(channel=c)[0] for c in input_channels])
+ input_format = _get_input_format(data_type=data_type)
+ output_format = _get_output_format(data_type=data_type, channel=channel)
+ if channel == "G":
+ return DeltaFAlgorithm(informat=input_format).process(
+ timeseries=input_timeseries
+ )[0]
+ elif channel in ["D", "X", "Y"] and data_type == "variation":
+ output_timeseries = XYZAlgorithm(
+ informat=input_format, outformat=output_format
+ ).process(timeseries=input_timeseries)
+ elif channel in ["D", "H"] and data_type == "adjusted":
+ output_timeseries = XYZAlgorithm(
+ informat=input_format, outformat=output_format
+ ).process(timeseries=input_timeseries)
+ else:
+ output_timeseries = input_timeseries.select(channel=channel)
+ return output_timeseries.select(channel=channel)[0]
+
+
+def get_timeseries(timeseries: Stream, channels: List[str]) -> Stream:
+ """calculates requested derived products from input timeseries"""
+ output_timeseries = Stream()
+ for channel in channels:
+ output_timeseries += compute_trace(timeseries=timeseries, channel=channel)
+ return output_timeseries
+
+
+def get_missing_channels(timeseries: Stream) -> List[str]:
+ """retrieves empty channels from stream"""
+ missing_channels = []
+ stats = timeseries[0].stats
+ for channel in get_channels(stream=timeseries):
+ if not has_all_channels(
+ stream=timeseries,
+ channels=(channel,),
+ starttime=stats.starttime,
+ endtime=stats.endtime,
+ ):
+ missing_channels.append(channel)
+ return missing_channels
+
+
+def get_derived_input_channels(channel: str, data_type: str) -> List[str]:
+ """returns channels needed to compute derived channel"""
+ if channel == "G":
+ if data_type == "variation":
+ return ["U", "V", "W", "F"]
+ else:
+ return ["X", "Y", "Z", "F"]
+ elif channel in ["X", "Y", "H", "D"] and data_type == "variation":
+ return ["U", "V"]
+ elif channel in ["H", "D"] and data_type != "variation":
+ return ["X", "Y"]
+ else:
+ # not a derived channel
+ return [channel]
+
+
+def get_required_channels(channels: List[str], data_type: str) -> List[str]:
+ """returns channels needed to compute derived channels"""
+ required_channels = []
+ for channel in channels:
+ input_channels = get_derived_input_channels(
+ channel=channel, data_type=data_type
+ )
+ required_channels.extend(
+ [c for c in input_channels if c not in required_channels]
+ )
+ return required_channels
+
+
+def _get_input_format(data_type: str) -> str:
+ if data_type == "variation":
+ return "obs"
+ return "geo"
+
+
+def _get_output_format(data_type: str, channel: str) -> str:
+ if data_type == "variation":
+ if channel in ["X", "Y"]:
+ return "geo"
+ elif channel == "D":
+ return "obsd"
+ return "mag"
diff --git a/geomagio/edge/EdgeFactory.py b/geomagio/edge/EdgeFactory.py
index 682d03f733ed2b24fa04269925859f3279f179ed..90a68111306d3eb14faa374d39ad17c49206ea2b 100644
--- a/geomagio/edge/EdgeFactory.py
+++ b/geomagio/edge/EdgeFactory.py
@@ -111,6 +111,7 @@ class EdgeFactory(TimeseriesFactory):
channels=None,
type=None,
interval=None,
+ derive_missing: bool = True,
):
"""Get timeseries data
@@ -128,6 +129,8 @@ class EdgeFactory(TimeseriesFactory):
data type.
interval: {'day', 'hour', 'minute', 'second', 'tenhertz'}
data interval.
+ dervie_missing
+ if True, calculates trivial derived products missing from initial request
Returns
-------
@@ -166,7 +169,8 @@ class EdgeFactory(TimeseriesFactory):
# restore stdout
sys.stdout = original_stdout
self._post_process(timeseries, starttime, endtime, channels)
-
+ if derive_missing:
+ timeseries = self.get_derived_timeseries(timeseries=timeseries)
return timeseries
def put_timeseries(
@@ -200,6 +204,11 @@ class EdgeFactory(TimeseriesFactory):
channel and that trace should have an ndarray, with nan's
representing gaps.
"""
+ if not TimeseriesUtility.has_any_channels(
+ stream=timeseries, channels=channels, starttime=starttime, endtime=endtime
+ ):
+ sys.stderr.write("No data to write\n")
+ return
stats = timeseries[0].stats
observatory = observatory or stats.station or self.observatory
channels = channels or self.channels
@@ -443,6 +452,7 @@ class EdgeFactory(TimeseriesFactory):
trace_send.data
)
trace_send = self._convert_trace_to_int(trace_send)
+ trace.stats.channel = sncl.element
ric.send_trace(interval, trace_send)
if self.forceout:
ric.forceout()
diff --git a/geomagio/edge/LegacySNCL.py b/geomagio/edge/LegacySNCL.py
index 1b8cd2169c2ac5b4e45c207bd116e343f2032f44..87f47d5f8be9f178ccc5a22e7f68712100260bda 100644
--- a/geomagio/edge/LegacySNCL.py
+++ b/geomagio/edge/LegacySNCL.py
@@ -31,7 +31,9 @@ class LegacySNCL(SNCL):
return LegacySNCL(
station=station,
network=network,
- channel=get_channel(element=element, interval=interval),
+ channel=get_channel(
+ data_type=data_type, element=element, interval=interval
+ ),
location=location or get_location(element=element, data_type=data_type),
)
@@ -55,9 +57,10 @@ class LegacySNCL(SNCL):
raise ValueError(f"Unexcepted interval code: {channel_start}")
-def get_channel(element: str, interval: str) -> str:
+def get_channel(data_type: str, element: str, interval: str) -> str:
return _check_predefined_channel(element=element, interval=interval) or (
- _get_channel_start(interval=interval) + _get_channel_end(element=element)
+ _get_channel_start(interval=interval)
+ + _get_channel_end(data_type=data_type, element=element)
)
@@ -87,6 +90,12 @@ def _get_channel_start(interval: str) -> str:
def _get_element(channel: str, location: str) -> str:
"""Translates channel/location to element"""
element_start = channel[2]
+ if element_start == "U":
+ element_start = "H"
+ elif element_start == "V":
+ element_start = "E"
+ elif element_start == "W":
+ element_start = "Z"
channel = channel
channel_middle = channel[1]
location_end = location[1]
@@ -116,7 +125,7 @@ def _check_predefined_channel(element: str, interval: str) -> Optional[str]:
return None
-def _get_channel_end(element: str) -> str:
+def _get_channel_end(data_type: str, element: str) -> str:
channel_middle = "V"
if "_Volt" in element:
channel_middle = "E"
@@ -127,6 +136,12 @@ def _get_channel_end(element: str) -> str:
elif element in ["F", "G"]:
channel_middle = "S"
channel_end = element.split("_")[0]
+ if channel_end == "U" and data_type == "variation":
+ channel_end = "H"
+ elif channel_end == "V" and data_type == "variation":
+ channel_end = "E"
+ elif channel_end == "W" and data_type == "variation":
+ channel_end = "Z"
return channel_middle + channel_end
diff --git a/geomagio/edge/MiniSeedFactory.py b/geomagio/edge/MiniSeedFactory.py
index d56bab7713c0684c199fa227a8d815d3621c9cb1..d94437114f6cd716ce123c94ecf3dac2d75abfd8 100644
--- a/geomagio/edge/MiniSeedFactory.py
+++ b/geomagio/edge/MiniSeedFactory.py
@@ -68,7 +68,7 @@ class MiniSeedFactory(TimeseriesFactory):
self,
host="cwbpub.cr.usgs.gov",
port=2061,
- write_port=7981,
+ write_port=7974,
observatory=None,
channels=None,
type=None,
@@ -97,6 +97,7 @@ class MiniSeedFactory(TimeseriesFactory):
channels=None,
type=None,
interval=None,
+ derive_missing: bool = True,
):
"""Get timeseries data
@@ -114,6 +115,8 @@ class MiniSeedFactory(TimeseriesFactory):
data type.
interval: {'day', 'hour', 'minute', 'second', 'tenhertz'}
data interval.
+ derive_missing
+ if True, calculates trivial derived products missing from initial request
Returns
-------
@@ -158,6 +161,8 @@ class MiniSeedFactory(TimeseriesFactory):
sys.stdout = original_stdout
self._post_process(timeseries, starttime, endtime, channels)
+ if derive_missing:
+ self.get_derived_timeseries(timeseries=timeseries)
return timeseries
def put_timeseries(
@@ -191,6 +196,11 @@ class MiniSeedFactory(TimeseriesFactory):
channel and that trace should have an ndarray, with nan's
representing gaps.
"""
+ if not TimeseriesUtility.has_any_channels(
+ stream=timeseries, channels=channels, starttime=starttime, endtime=endtime
+ ):
+ sys.stderr.write("No data to write\n")
+ return
stats = timeseries[0].stats
observatory = observatory or stats.station or self.observatory
channels = channels or self.channels
@@ -329,6 +339,10 @@ class MiniSeedFactory(TimeseriesFactory):
data = self.client.get_waveforms(
sncl.network, sncl.station, sncl.location, sncl.channel, starttime, endtime
)
+ if len(data) > 1:
+ data = obspy.core.Stream(
+ [d for d in data if d.data.dtype != "int32"]
+ ) # for shared channels names between earthworm and miniseed(T1-4)
data.merge()
if data.count() == 0:
data += TimeseriesUtility.create_empty_trace(
diff --git a/geomagio/edge/MiniSeedInputClient.py b/geomagio/edge/MiniSeedInputClient.py
index a25d31ce89741d8483aa9a1bf55220b1d2c3358e..d2dc2e148e6bb8351b846bc2ccca4d268a544a35 100644
--- a/geomagio/edge/MiniSeedInputClient.py
+++ b/geomagio/edge/MiniSeedInputClient.py
@@ -25,7 +25,7 @@ class MiniSeedInputClient(object):
Floating point precision for output data
"""
- def __init__(self, host, port=2061, encoding="float32"):
+ def __init__(self, host, port=2061, encoding="float64"):
self.host = host
self.port = port
self.encoding = encoding
diff --git a/geomagio/processing/__init__.py b/geomagio/processing/__init__.py
index be97cfbd5e32485df54fc814fd646c595685fed7..a429346ca2b1f86fe0aa54564cf1c56321be0a1d 100644
--- a/geomagio/processing/__init__.py
+++ b/geomagio/processing/__init__.py
@@ -3,21 +3,26 @@
Note that these implementations are subject to change,
and should be considered less stable than other packages in the library.
"""
+from .derived import adjusted, average, sqdist
from .factory import get_edge_factory, get_miniseed_factory
-from .observatory import adjusted, average, deltaf, rotate, sqdist_minute
-from .obsrio import obsrio_minute, obsrio_second, obsrio_temperatures, obsrio_tenhertz
+from .obsrio import (
+ filter_days,
+ filter_hours,
+ filter_minutes,
+ filter_temperatures,
+ filter_tenhertz,
+)
__all__ = [
"adjusted",
"average",
- "deltaf",
+ "filter_days",
+ "filter_hours",
+ "filter_minutes",
+ "filter_temperatures",
+ "filter_tenhertz",
"get_edge_factory",
"get_miniseed_factory",
- "obsrio_minute",
- "obsrio_second",
- "obsrio_temperatures",
- "obsrio_tenhertz",
- "rotate",
- "sqdist_minute",
+ "sqdist",
]
diff --git a/geomagio/processing/derived.py b/geomagio/processing/derived.py
new file mode 100644
index 0000000000000000000000000000000000000000..e8f4a565fc11f902313bd4dce139d5082f32064f
--- /dev/null
+++ b/geomagio/processing/derived.py
@@ -0,0 +1,108 @@
+from typing import Optional, List
+
+from obspy import UTCDateTime
+
+from ..algorithm import (
+ AdjustedAlgorithm,
+ AverageAlgorithm,
+ SqDistAlgorithm,
+)
+from ..Controller import Controller, get_realtime_interval
+from ..edge import MiniSeedFactory
+from .factory import get_miniseed_factory
+
+
+def adjusted(
+ observatory: str,
+ statefile: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
+):
+ """ Transform 1Hz UVWF to 1HZ XYZF """
+ controller = Controller(
+ inputFactory=input_factory or get_miniseed_factory(data_type="variation"),
+ outputFactory=output_factory or get_miniseed_factory(data_type="adjusted"),
+ algorithm=AdjustedAlgorithm(
+ statefile=statefile,
+ data_type="adjusted",
+ location="A0",
+ ),
+ inputInterval="second",
+ outputInterval="second",
+ )
+ controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
+ starttime=starttime,
+ endtime=endtime,
+ input_channels=("U", "V", "W", "F"),
+ output_channels=("X", "Y", "Z", "F"),
+ realtime=endtime - starttime,
+ update_limit=update_limit,
+ )
+
+
+def average(
+ observatories: List[str],
+ input_channel: str,
+ output_channel: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
+):
+ """ Computes single channel average across multiple observatories """
+ controller = Controller(
+ inputFactory=input_factory,
+ outputFactory=output_factory,
+ algorithm=AverageAlgorithm(observatories=observatories, channel=output_channel),
+ inputInterval="minute",
+ outputInterval="minute",
+ )
+ controller.run_as_update(
+ observatory=observatories,
+ output_observatory="USGS",
+ starttime=starttime,
+ endtime=endtime,
+ input_channels=(input_channel,),
+ output_channels=(output_channel,),
+ realtime=endtime - starttime,
+ update_limit=update_limit,
+ )
+
+
+def sqdist(
+ observatory: str,
+ statefile: str,
+ realtime_interval: int = 1800,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
+):
+ """ Computes H_SQ, H_SV, and H_Dist """
+ starttime, endtime = get_realtime_interval(interval_seconds=realtime_interval)
+ controller = Controller(
+ inputFactory=input_factory or get_miniseed_factory(interval="minute"),
+ outputFactory=output_factory or get_miniseed_factory(interval="minute"),
+ algorithm=SqDistAlgorithm(
+ alpha=2.3148e-5,
+ gamma=3.3333e-2,
+ beta=0,
+ m=1440,
+ mag=True,
+ smooth=180,
+ statefile=statefile,
+ ),
+ inputInterval="minute",
+ outputInterval="minute",
+ )
+ controller.run(
+ observatory=(observatory,),
+ starttime=starttime,
+ endtime=endtime,
+ input_channels=("X", "Y", "Z", "F"),
+ output_channels=("H_Dist", "H_SQ", "H_SV"),
+ )
diff --git a/geomagio/processing/factory.py b/geomagio/processing/factory.py
index d9a065497d9cf6b0a38c391355e11a9d698795d2..a12bd93e6c5ceab6df5dc5bc11c3316d08383ddd 100644
--- a/geomagio/processing/factory.py
+++ b/geomagio/processing/factory.py
@@ -1,27 +1,31 @@
import os
-from ..TimeseriesFactory import TimeseriesFactory
from ..edge import EdgeFactory, MiniSeedFactory
def get_edge_factory(
- data_type="variation", interval="second", **kwargs
-) -> TimeseriesFactory:
+ host=os.getenv("DATA_HOST", "127.0.0.1"),
+ data_type="variation",
+ interval="second",
+ **kwargs,
+):
return EdgeFactory(
- host=os.getenv("EDGE_HOST", "127.0.0.1"),
- interval=interval,
+ host=host,
type=data_type,
- **kwargs
+ interval=interval,
+ **kwargs,
)
def get_miniseed_factory(
- data_type="variation", interval="second", **kwargs
-) -> TimeseriesFactory:
+ host=os.getenv("DATA_HOST", "127.0.0.1"),
+ data_type="variation",
+ interval="second",
+ **kwargs,
+):
return MiniSeedFactory(
- convert_channels=("U", "V", "W"),
- host=os.getenv("EDGE_HOST", "127.0.0.1"),
- interval=interval,
+ host=host,
type=data_type,
- **kwargs
+ interval=interval,
+ **kwargs,
)
diff --git a/geomagio/processing/observatory.py b/geomagio/processing/observatory.py
index 85a9cfd7f033c4e72bb3401456f939c4e154a8c1..d31b31517c468d4d4c8e492f35de16c14cd43af1 100644
--- a/geomagio/processing/observatory.py
+++ b/geomagio/processing/observatory.py
@@ -1,243 +1,372 @@
-from typing import List, Optional
+from enum import Enum
+import json
+from typing import List, Optional, Tuple
-import numpy
+from obspy import UTCDateTime
+import typer
-from ..algorithm import (
- AdjustedAlgorithm,
- AverageAlgorithm,
- DeltaFAlgorithm,
- SqDistAlgorithm,
- XYZAlgorithm,
-)
-from ..Controller import Controller, get_realtime_interval
-from ..TimeseriesFactory import TimeseriesFactory
-from .factory import get_edge_factory
+from ..Controller import get_realtime_interval
+from .factory import get_edge_factory, get_miniseed_factory
+from . import obsrio, derived, pcdcp
-def adjusted(
- observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- interval: str = "second",
- output_factory: Optional[TimeseriesFactory] = None,
- matrix: Optional[numpy.ndarray] = None,
- pier_correction: Optional[float] = None,
- statefile: Optional[str] = None,
+app = typer.Typer()
+
+
+def main():
+ app()
+
+
+class DataFormat(str, Enum):
+ OBSRIO = "obsrio"
+ PCDCP = "pcdcp"
+
+
+class DataType(str, Enum):
+ ADJUSTED = "adjusted"
+ VARIATION = "variation"
+
+
+@app.command(
+ name="average",
+ short_help="computes Dst-3 and Dst-4 from 1-minute variation data",
+ help="""
+ Dst-3: Computed from HON, GUA and SJG(H_Dist)
+
+ Dst-4: Computed from HON, SJG, HER, KAK(H_Dist)
+ """,
+)
+def compute_average(
+ data_type: DataType = DataType.VARIATION,
realtime_interval: int = 600,
+ update_limit: int = 10,
+ input_host: str = "127.0.0.1",
+ output_host: str = "127.0.0.1",
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
):
- """Run Adjusted algorithm.
-
- Parameters
- ----------
- observatory: observatory to calculate
- input_factory: where to read, should be configured with data_type
- interval: data interval
- output_factory: where to write, should be configured with data_type
- matrix: adjusted matrix
- pier_correction: adjusted pier correction
- statefile: adjusted statefile
- realtime_interval: window in seconds
-
- Uses update_limit=10.
- """
- if not statefile and (not matrix or not pier_correction):
- raise ValueError("Either statefile or matrix and pier_correction are required.")
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=AdjustedAlgorithm(
- matrix=matrix,
- pier_correction=pier_correction,
- statefile=statefile,
- data_type="adjusted",
- location="A0",
- ),
- inputFactory=input_factory or get_edge_factory(data_type="variation"),
- inputInterval=interval,
- outputFactory=output_factory or get_edge_factory(data_type="adjusted"),
- outputInterval=interval,
+ starttime, endtime = get_processing_interval(
+ realtime_interval=realtime_interval,
+ starttime=starttime,
+ endtime=endtime,
+ )
+ input_factory = get_miniseed_factory(host=input_host, data_type=data_type)
+ output_factory = get_miniseed_factory(host=output_host, data_type=data_type)
+ derived.average(
+ observatories=["HON", "SJG", "HER", "KAK"],
+ input_channel="H_Dist",
+ output_channel="UX4",
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
)
- controller.run_as_update(
- observatory=(observatory,),
- output_observatory=(observatory,),
+ derived.average(
+ observatories=["HON", "GUA", "SJG"],
+ input_channel="H_Dist",
+ output_channel="UX3",
starttime=starttime,
endtime=endtime,
- input_channels=("H", "E", "Z", "F"),
- output_channels=("X", "Y", "Z", "F"),
- realtime=realtime_interval,
- update_limit=10,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
)
-def average(
- observatories: List[str],
- input_channel: str,
- input_factory: Optional[TimeseriesFactory] = None,
- interval: str = "second",
- output_channel: str = None,
- output_factory: Optional[TimeseriesFactory] = None,
- output_observatory: str = "USGS",
+@app.command(
+ name="adjusted",
+ short_help="Computes adjusted XYZF and derived channels at 1Hz and 1 minute",
+ help="""
+ Inputs: U, V, W, F(1Hz)
+
+ Outputs: X, Y, Z, F, H_Dist, H_SQ, H_SV(1Hz/1 min)
+ """,
+)
+def process_adjusted(
+ observatory: str,
+ adjusted_statefile: str,
+ sqdist_statefile: str,
realtime_interval: int = 600,
+ update_limit: int = 10,
+ input_host: str = "127.0.0.1",
+ output_host: str = "127.0.0.1",
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
):
- """Run Average algorithm.
-
- Parameters
- ----------
- observatories: input observatories to calculate
- input_channel: channel from multiple observatories to average
- input_factory: where to read, should be configured with data_type and interval
- interval: data interval
- output_channel: channel to write (defaults to input_channel).
- output_factory: where to write, should be configured with data_type and interval
- output_observatory: observatory where output is written
- realtime_interval: window in seconds
-
- Uses update_limit=10.
- """
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=AverageAlgorithm(observatories=observatories, channel=output_channel),
- inputFactory=input_factory or get_edge_factory(),
- inputInterval=interval,
- outputFactory=output_factory or get_edge_factory(),
- outputInterval=interval,
+ process_sqdist = False
+ if starttime is None and endtime is None:
+ process_sqdist = True
+ starttime, endtime = get_processing_interval(
+ realtime_interval=realtime_interval,
+ starttime=starttime,
+ endtime=endtime,
+ )
+ input_factory = get_miniseed_factory(host=input_host, data_type="adjusted")
+ output_factory = get_miniseed_factory(host=output_host, data_type="adjusted")
+ derived.adjusted(
+ observatory=observatory,
+ statefile=adjusted_statefile,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=get_miniseed_factory(host=input_host, data_type="variation"),
+ output_factory=output_factory,
)
- controller.run_as_update(
- observatory=observatories,
- output_observatory=(output_observatory,),
+ obsrio.filter_minutes(
+ observatory=observatory,
starttime=starttime,
endtime=endtime,
- output_channels=(output_channel or input_channel,),
- realtime=realtime_interval,
- update_limit=10,
+ update_limit=update_limit,
+ channels=("X", "Y", "Z", "F"),
+ input_factory=input_factory,
+ output_factory=output_factory,
)
+ if process_sqdist and valid_sqdist_statefile(
+ statefile=sqdist_statefile, starttime=starttime
+ ):
+ derived.sqdist(
+ observatory=observatory,
+ statefile=sqdist_statefile,
+ realtime_interval=realtime_interval,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ pcdcp.legacy_adjusted(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=get_edge_factory(host=output_host, data_type="adjusted"),
+ )
+
+
+@app.command(
+ name="day",
+ short_help="Computes daily UVWF and derived channels(variation/adjusted)",
+ help="""
+ Variation:
+
+ Inputs: U, V, W, F(1 min)
+
+ Outputs: U, V, W, F(1 day)
+
+ Adjusted:
+ Inputs: X, Y, Z, F(1 min)
-def deltaf(
+ Outputs: X, Y, Z, F(1 day)
+ """,
+)
+def process_day(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- interval: str = "second",
- output_factory: Optional[TimeseriesFactory] = None,
- deltaf_from="obs",
- realtime_interval: int = 600,
+ data_type: DataType = DataType.VARIATION,
+ realtime_interval: int = 86400,
+ update_limit: int = 7,
+ input_host: str = "127.0.0.1",
+ output_host: str = "127.0.0.1",
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
):
- """Run Delta-F algorithm.
-
- Parameters
- ----------
- observatory: observatory to calculate
- input_factory: where to read, should be configured with data_type and interval
- output_factory: where to write, should be configured with data_type and interval
- deltaf_from: one of {"obs", "mag", "geo"}
- realtime_interval: window in seconds
-
- Uses update_limit=10.
- """
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=DeltaFAlgorithm(informat=deltaf_from),
- inputFactory=input_factory or get_edge_factory(),
- inputInterval=interval,
- outputFactory=output_factory or get_edge_factory(),
- outputInterval=interval,
+ starttime, endtime = get_processing_interval(
+ realtime_interval=realtime_interval,
+ starttime=starttime,
+ endtime=endtime,
)
- controller.run_as_update(
- observatory=(observatory,),
- output_observatory=(observatory,),
+ input_factory = get_miniseed_factory(host=input_host, data_type=data_type)
+ output_factory = get_miniseed_factory(host=output_host, data_type=data_type)
+ obsrio.filter_days(
+ observatory=observatory,
starttime=starttime,
endtime=endtime,
- output_channels=("G",),
- realtime=realtime_interval,
- update_limit=10,
+ channels=("U", "V", "W", "F")
+ if data_type == DataType.VARIATION
+ else ("X", "Y", "Z", "F"),
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
)
-def rotate(
+@app.command(
+ name="hour",
+ short_help="Computes hourly UVWF and derived channels(variation/adjusted)",
+ help="""
+ Variation:
+
+ Inputs: U, V, W, F(1 min)
+
+ Outputs: U, V, W, F(1 hour)
+
+ Adjusted:
+
+ Inputs: X, Y, Z, F(1 min)
+
+ Outputs: X, Y, Z, F(1 hour)
+ """,
+)
+def process_hour(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- interval: str = "second",
- output_channels=("X", "Y"),
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
- xyz_from="obs",
- xyz_to="geo",
+ data_type: DataType = DataType.VARIATION,
+ realtime_interval: int = 3600,
+ update_limit: int = 7,
+ input_host: str = "127.0.0.1",
+ output_host: str = "127.0.0.1",
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
):
- """Run XYZ rotation algorithm.
-
- Parameters
- ----------
- observatory: observatory to calculate
- input_factory: where to read, should be configured with data_type and interval
- output_channels: which channels to write
- output_factory: where to write, should be configured with data_type and interval
- realtime_interval: window in seconds
- xyz_from: one of {"obs", "mag", "geo"}
- xyz_to: one of {"obs", "obsd", "mag", "geo"}
-
- Uses update_limit=10.
- """
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=XYZAlgorithm(informat=xyz_from, outformat=xyz_to),
- inputFactory=input_factory or get_edge_factory(),
- inputInterval=interval,
- outputFactory=output_factory or get_edge_factory(),
- outputInterval=interval,
+ starttime, endtime = get_processing_interval(
+ realtime_interval=realtime_interval,
+ starttime=starttime,
+ endtime=endtime,
)
- controller.run_as_update(
- observatory=(observatory,),
- output_observatory=(observatory,),
+ input_factory = get_miniseed_factory(host=input_host, data_type=data_type)
+ output_factory = get_miniseed_factory(host=output_host, data_type=data_type)
+ obsrio.filter_hours(
+ observatory=observatory,
starttime=starttime,
endtime=endtime,
- output_channels=output_channels,
- realtime=realtime_interval,
- update_limit=10,
+ channels=("U", "V", "W", "F")
+ if data_type == DataType.VARIATION
+ else ("X", "Y", "Z", "F"),
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
)
-def sqdist_minute(
+@app.command(
+ name="variation",
+ short_help="Computes variation UVWF and derived channels at 1Hz and 1 minute",
+ help="""
+ Inputs:
+
+ PCDCP:
+
+ H, E, Z, F(1Hz)
+
+ T1-4(1 min)
+
+ OBSRIO:
+
+ U, V ,W(10Hz)
+
+ F, T1-4(1Hz)
+
+ Outputs:
+
+ U, V, W, F(1 Hz)
+
+ U, V, W, F, T1-4, H_SQ, H_SV, H_Dist(1 min)
+ """,
+)
+def process_variation(
observatory: str,
- statefile: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 1800,
+ sqdist_statefile: str,
+ data_format: DataFormat = DataFormat.PCDCP,
+ realtime_interval: int = 600,
+ update_limit: int = 10,
+ input_host: str = "127.0.0.1",
+ output_host: str = "127.0.0.1",
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
):
- """Run SqDist algorithm.
-
- Only supports "minute" interval.
-
- Parameters
- ----------
- observatory: observatory to calculate
- statefile: sqdist statefile must already exist
- input_factory: where to read, should be configured with data_type and interval
- output_factory: where to write, should be configured with data_type and interval
- realtime_interval: window in seconds
- """
- if not statefile:
- raise ValueError("Statefile is required.")
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=SqDistAlgorithm(
- alpha=2.3148e-5,
- gamma=3.3333e-2,
- m=1440,
- mag=True,
- smooth=180,
- statefile=statefile,
- ),
- inputFactory=input_factory or get_edge_factory(interval="minute"),
- inputInterval="minute",
- outputFactory=output_factory or get_edge_factory(interval="minute"),
- outputInterval="minute",
+ process_sqdist = False
+ if starttime is None and endtime is None:
+ process_sqdist = True
+ starttime, endtime = get_processing_interval(
+ realtime_interval=realtime_interval,
+ starttime=starttime,
+ endtime=endtime,
+ )
+ input_factory = get_miniseed_factory(host=input_host)
+ output_factory = get_miniseed_factory(host=output_host)
+ if data_format == DataFormat.PCDCP:
+ pcdcp.prepare(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=get_edge_factory(host=input_host),
+ output_factory=output_factory,
+ )
+ else:
+ obsrio.prepare(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=get_miniseed_factory(
+ host=input_host, convert_channels=("U", "V", "W")
+ ),
+ output_factory=output_factory,
+ )
+ obsrio.copy_channels(
+ observatory=observatory,
+ channels=[
+ ("UK1", "UK1"),
+ ("UK2", "UK2"),
+ ("UK3", "UK3"),
+ ("UK4", "UK4"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="minute",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=get_edge_factory(host=input_host),
+ )
+ obsrio.filter_minutes(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
)
- # sqdist is stateful, use run
- controller.run(
- observatory=(observatory,),
- output_observatory=(observatory,),
+ if process_sqdist and valid_sqdist_statefile(
+ statefile=sqdist_statefile, starttime=starttime
+ ):
+ derived.sqdist(
+ observatory=observatory,
+ statefile=sqdist_statefile,
+ realtime_interval=realtime_interval,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ pcdcp.legacy_variation(
+ observatory=observatory,
starttime=starttime,
endtime=endtime,
- input_channels=("X", "Y", "Z", "F"),
- output_channels=("MDT", "MSQ", "MSV"),
- realtime=realtime_interval,
- rename_output_channel=(("H_Dist", "MDT"), ("H_SQ", "MSQ"), ("H_SV", "MSV")),
- update_limit=10,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=get_edge_factory(host=output_host),
)
+
+
+def get_processing_interval(
+ realtime_interval: int,
+ starttime: Optional[str] = None,
+ endtime: Optional[str] = None,
+) -> Tuple[UTCDateTime, UTCDateTime]:
+ if starttime is not None and endtime is not None:
+ return (UTCDateTime(starttime), UTCDateTime(endtime))
+ elif starttime is not None and endtime is None:
+ raise ValueError("endtime not provided")
+ elif starttime is None and endtime is not None:
+ raise ValueError("starttime not provided")
+ else:
+ return get_realtime_interval(interval_seconds=realtime_interval)
+
+
+def valid_sqdist_statefile(statefile: str, starttime: UTCDateTime) -> bool:
+ try:
+ with open(statefile, "r") as f:
+ data = f.read()
+ data = json.loads(data)
+ if data["next_starttime"] > starttime:
+ return False
+ return True
+ except:
+ return False
diff --git a/geomagio/processing/obsrio.py b/geomagio/processing/obsrio.py
index bf0acf2b81687c6a3f590e6b4cd8b902f4730aa3..2db021e02b9b3ce52c90208ecba5f27d6369f94c 100644
--- a/geomagio/processing/obsrio.py
+++ b/geomagio/processing/obsrio.py
@@ -1,341 +1,262 @@
-from typing import Optional
+from typing import List, Optional, Tuple
-import typer
+from obspy import UTCDateTime
from ..algorithm import Algorithm, FilterAlgorithm
-from ..edge import EdgeFactory, MiniSeedFactory
-from ..Controller import (
- Controller,
- get_realtime_interval,
-)
-from ..TimeseriesFactory import TimeseriesFactory
+from ..Controller import Controller
+from ..edge import MiniSeedFactory
+from .. import TimeseriesFactory
from .factory import get_edge_factory, get_miniseed_factory
-def main():
- typer.run(obsrio_filter)
-
-
-def obsrio_filter(
- interval: str,
+def copy_channels(
observatory: str,
- input_factory: Optional[str] = None,
- host: str = "127.0.0.1",
- port: str = 2061,
- output_factory: Optional[str] = None,
- output_port: int = typer.Option(
- 2061, help="Port where output factory writes data."
- ),
- output_read_port: int = typer.Option(
- 2061, help="Port where output factory reads data"
- ),
- realtime_interval: int = 600,
+ channels: List[Tuple[str]],
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ data_interval: str = "second",
update_limit: int = 10,
+ input_factory: Optional[TimeseriesFactory] = None,
+ output_factory: Optional[TimeseriesFactory] = None,
):
- if interval == "realtime":
- filter_realtime(
- observatory=observatory,
- input_factory=input_factory,
- host=host,
- port=port,
- output_factory=output_factory,
- output_port=output_port,
- output_read_port=output_read_port,
- realtime_interval=realtime_interval,
+ """ Copy channels between legacy and miniseed formats """
+ controller = Controller(
+ inputFactory=input_factory or get_miniseed_factory(),
+ outputFactory=output_factory or get_edge_factory(),
+ inputInterval=data_interval,
+ outputInterval=data_interval,
+ )
+ for input_channel, output_channel in channels:
+ controller.run_as_update(
+ algorithm=Algorithm(
+ inchannels=input_channel,
+ outchannels=output_channel,
+ ),
+ observatory=(observatory,),
+ output_observatory=(observatory,),
+ starttime=starttime,
+ endtime=endtime,
+ input_channels=(input_channel,),
+ output_channels=(output_channel,),
+ rename_output_channel=((input_channel, output_channel),),
+ realtime=endtime - starttime,
update_limit=update_limit,
)
- elif interval in ["hour", "day"]:
- input_factory = EdgeFactory(host=host, port=port)
- output_factory = MiniSeedFactory(
- host=host, port=output_read_port, write_port=output_port
- )
- if interval == "hour":
- obsrio_hour(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
- elif interval == "day":
- obsrio_day(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
- else:
- raise ValueError("Invalid interval")
-def filter_realtime(
+def filter_days(
observatory: str,
- input_factory: Optional[str] = None,
- host: str = "127.0.0.1",
- port: str = 2061,
- output_factory: Optional[str] = None,
- output_port: int = typer.Option(
- 2061, help="Port where output factory writes data."
- ),
- output_read_port: int = typer.Option(
- 2061, help="Port where output factory reads data"
- ),
- realtime_interval: int = 600,
- update_limit: int = 10,
-):
- """Filter 10Hz miniseed, 1 second, one minute, and temperature data.
- Defaults set for realtime processing; can also be implemented to update legacy data"""
- if input_factory == "miniseed":
- input_factory = MiniSeedFactory(host=host, port=port)
- elif input_factory == "edge":
- input_factory = EdgeFactory(host=host, port=port)
- if output_factory == "miniseed":
- output_factory = MiniSeedFactory(
- host=host, port=output_read_port, write_port=output_port
- )
- elif output_factory == "edge":
- output_factory = EdgeFactory(
- host=host, port=output_read_port, write_port=output_port
- )
-
- obsrio_tenhertz(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
- obsrio_second(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
- obsrio_minute(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
- obsrio_temperatures(
- observatory=observatory,
- input_factory=input_factory,
- output_factory=output_factory,
- realtime_interval=realtime_interval,
- update_limit=update_limit,
- )
-
-
-def obsrio_day(
- observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 86400,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
update_limit: int = 7,
+ channels: List[str] = ("U", "V", "W", "F"),
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
):
- """Filter 1 second edge H,E,Z,F to 1 day miniseed U,V,W,F."""
- starttime, endtime = get_realtime_interval(realtime_interval)
+ """ Filter 1 min UVWF to 1 day UVWF """
controller = Controller(
- inputFactory=input_factory or get_edge_factory(),
+ inputFactory=input_factory or get_miniseed_factory(interval="minute"),
+ outputFactory=output_factory or get_miniseed_factory(interval="day"),
inputInterval="minute",
- outputFactory=output_factory or get_miniseed_factory(),
outputInterval="day",
)
- renames = {"H": "U", "E": "V", "Z": "W", "F": "F"}
- for input_channel in renames.keys():
- output_channel = renames[input_channel]
+ for channel in channels:
controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
algorithm=FilterAlgorithm(
input_sample_period=60.0,
- output_sample_period=86400.0,
- inchannels=(input_channel,),
- outchannels=(output_channel,),
+ output_sample_period=3600.0,
+ inchannels=(channel,),
+ outchannels=(channels,),
),
- observatory=(observatory,),
- output_observatory=(observatory,),
starttime=starttime,
endtime=endtime,
- input_channels=(input_channel,),
- output_channels=(output_channel,),
- realtime=realtime_interval,
- rename_output_channel=((input_channel, output_channel),),
+ input_channels=(channel,),
+ output_channels=(channel,),
+ realtime=endtime - starttime,
update_limit=update_limit,
)
-def obsrio_hour(
+def filter_hours(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
- update_limit: int = 10,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 24,
+ channels: List[str] = ("U", "V", "W", "F"),
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
):
- """Filter 1 second edge H,E,Z,F to 1 hour miniseed U,V,W,F."""
- starttime, endtime = get_realtime_interval(realtime_interval)
+ """ Filter 1 min UVWF to 1 hour UVWF """
controller = Controller(
- inputFactory=input_factory or get_edge_factory(),
+ inputFactory=input_factory or get_miniseed_factory(interval="minute"),
+ outputFactory=output_factory or get_miniseed_factory(interval="hour"),
inputInterval="minute",
- outputFactory=output_factory or get_miniseed_factory(),
outputInterval="hour",
)
- renames = {"H": "U", "E": "V", "Z": "W", "F": "F"}
- for input_channel in renames.keys():
- output_channel = renames[input_channel]
+ for channel in channels:
controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
algorithm=FilterAlgorithm(
input_sample_period=60.0,
output_sample_period=3600.0,
- inchannels=(input_channel,),
- outchannels=(output_channel,),
+ inchannels=(channel,),
+ outchannels=(channels,),
),
- observatory=(observatory,),
- output_observatory=(observatory,),
starttime=starttime,
endtime=endtime,
- input_channels=(input_channel,),
- output_channels=(output_channel,),
- realtime=realtime_interval,
- rename_output_channel=((input_channel, output_channel),),
+ input_channels=(channel,),
+ output_channels=(channel,),
+ realtime=endtime - starttime,
update_limit=update_limit,
)
-def obsrio_minute(
+def filter_minutes(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
update_limit: int = 10,
+ channels: List[str] = ("U", "V", "W", "F"),
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
):
- """Filter 1Hz legacy H,E,Z,F to 1 minute legacy.
-
- Should be called after obsrio_second() and obsrio_tenhertz(),
- which populate 1Hz legacy H,E,Z,F.
- """
- starttime, endtime = get_realtime_interval(realtime_interval)
+ """ Filter 1Hz UVWF to 1 min UVWF """
controller = Controller(
- inputFactory=input_factory or get_edge_factory(),
+ inputFactory=input_factory or get_miniseed_factory(),
+ outputFactory=output_factory or get_miniseed_factory(),
inputInterval="second",
- outputFactory=output_factory or get_edge_factory(),
outputInterval="minute",
)
- for channel in ["H", "E", "Z", "F"]:
+ for channel in channels:
controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
algorithm=FilterAlgorithm(
- input_sample_period=1,
- output_sample_period=60,
+ input_sample_period=1.0,
+ output_sample_period=60.0,
inchannels=(channel,),
outchannels=(channel,),
),
- observatory=(observatory,),
- output_observatory=(observatory,),
starttime=starttime,
endtime=endtime,
input_channels=(channel,),
output_channels=(channel,),
- realtime=realtime_interval,
+ realtime=endtime - starttime,
update_limit=update_limit,
)
-def obsrio_second(
+def filter_temperatures(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
):
- """Copy 1Hz miniseed F to 1Hz legacy F."""
- starttime, endtime = get_realtime_interval(realtime_interval)
- controller = Controller(
- algorithm=Algorithm(inchannels=("F",), outchannels=("F",)),
- inputFactory=input_factory or get_miniseed_factory(),
- outputFactory=output_factory or get_edge_factory(),
- )
- controller.run_as_update(
- observatory=(observatory,),
- output_observatory=(observatory,),
- starttime=starttime,
- endtime=endtime,
- input_channels=("F",),
- output_channels=("F",),
- realtime=realtime_interval,
- update_limit=update_limit,
- )
-
-
-def obsrio_temperatures(
- observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
- update_limit: int = 10,
-):
- """Filter temperatures 1Hz miniseed (LK1-4) to 1 minute legacy (UK1-4)."""
- starttime, endtime = get_realtime_interval(realtime_interval)
+ """ Filter 1Hz T1-4 to 1 min T1-4 """
controller = Controller(
inputFactory=input_factory or get_miniseed_factory(),
+ outputFactory=output_factory or get_miniseed_factory(),
inputInterval="second",
- outputFactory=output_factory or get_edge_factory(),
outputInterval="minute",
)
- renames = {"LK1": "UK1", "LK2": "UK2", "LK3": "UK3", "LK4": "UK4"}
- for input_channel in renames.keys():
- output_channel = renames[input_channel]
+ input_channels = ["LK1", "LK2", "LK3", "LK4"]
+ output_channels = ["UK1", "UK2", "UK3", "UK4"]
+ for input_channel, output_channel in zip(input_channels, output_channels):
controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
algorithm=FilterAlgorithm(
- input_sample_period=1,
- output_sample_period=60,
+ input_sample_period=1.0,
+ output_sample_period=60.0,
inchannels=(input_channel,),
outchannels=(output_channel,),
),
- observatory=(observatory,),
- output_observatory=(observatory,),
starttime=starttime,
endtime=endtime,
input_channels=(input_channel,),
output_channels=(output_channel,),
- realtime=realtime_interval,
+ realtime=endtime - starttime,
rename_output_channel=((input_channel, output_channel),),
update_limit=update_limit,
)
-def obsrio_tenhertz(
+def filter_tenhertz(
observatory: str,
- input_factory: Optional[TimeseriesFactory] = None,
- output_factory: Optional[TimeseriesFactory] = None,
- realtime_interval: int = 600,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
):
- """Filter 10Hz miniseed U,V,W to 1Hz legacy H,E,Z."""
- starttime, endtime = get_realtime_interval(realtime_interval)
+ """ Filter 10HZ UVW to 1Hz UVW """
controller = Controller(
- inputFactory=input_factory or get_miniseed_factory(),
+ inputFactory=input_factory
+ or get_miniseed_factory(convert_channels=("U", "V", "W")),
+ outputFactory=output_factory or get_miniseed_factory(),
inputInterval="tenhertz",
- outputFactory=output_factory or get_edge_factory(),
outputInterval="second",
)
- renames = {"U": "H", "V": "E", "W": "Z"}
- for input_channel in renames.keys():
- output_channel = renames[input_channel]
+ for channel in ["U", "V", "W"]:
controller.run_as_update(
+ observatory=(observatory,),
+ output_observatory=(observatory,),
algorithm=FilterAlgorithm(
input_sample_period=0.1,
- output_sample_period=1,
- inchannels=(input_channel,),
- outchannels=(output_channel,),
+ output_sample_period=1.0,
+ inchannels=(channel,),
+ outchannels=(channel,),
),
- observatory=(observatory,),
- output_observatory=(observatory,),
starttime=starttime,
endtime=endtime,
- input_channels=(input_channel,),
- output_channels=(output_channel,),
- realtime=realtime_interval,
- rename_output_channel=((input_channel, output_channel),),
+ input_channels=(channel,),
+ output_channels=(channel,),
+ realtime=endtime - starttime,
update_limit=update_limit,
)
+
+
+def prepare(
+ observatory: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
+):
+ """Filter 10Hz UVW to 1Hz UVW
+ Filter 1Hz T1-4 to 1 min T1-4
+ Duplicate 1Hz W to Z
+ """
+ filter_tenhertz(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ input_factory.convert_channels = []
+ filter_temperatures(
+ observatory=observatory,
+ starttime=starttime,
+ endtime=endtime,
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ copy_channels(
+ observatory=observatory,
+ channels=[("W", "Z")],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="second",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
diff --git a/geomagio/processing/pcdcp.py b/geomagio/processing/pcdcp.py
new file mode 100644
index 0000000000000000000000000000000000000000..f19670118e445b9c75ba2fa761e4011fbd56bb85
--- /dev/null
+++ b/geomagio/processing/pcdcp.py
@@ -0,0 +1,160 @@
+from typing import List, Optional, Tuple
+
+from obspy import UTCDateTime
+
+from ..algorithm import Algorithm
+from ..Controller import Controller
+from ..edge import EdgeFactory, MiniSeedFactory
+from .. import TimeseriesFactory
+from .factory import get_edge_factory, get_miniseed_factory
+
+
+def copy_channels(
+ observatory: str,
+ channels: List[Tuple[str]],
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ data_interval: str = "second",
+ update_limit: int = 10,
+ input_factory: Optional[TimeseriesFactory] = None,
+ output_factory: Optional[TimeseriesFactory] = None,
+):
+ """ Copy channels between legacy and miniseed formats """
+ controller = Controller(
+ inputFactory=input_factory or get_edge_factory(),
+ outputFactory=output_factory or get_miniseed_factory(),
+ inputInterval=data_interval,
+ outputInterval=data_interval,
+ )
+ for input_channel, output_channel in channels:
+ controller.run_as_update(
+ algorithm=Algorithm(
+ inchannels=input_channel,
+ outchannels=output_channel,
+ ),
+ observatory=(observatory,),
+ output_observatory=(observatory,),
+ starttime=starttime,
+ endtime=endtime,
+ input_channels=(input_channel,),
+ output_channels=(output_channel,),
+ rename_output_channel=((input_channel, output_channel),),
+ realtime=endtime - starttime,
+ update_limit=update_limit,
+ )
+
+
+def legacy_adjusted(
+ observatory: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[EdgeFactory] = None,
+):
+ """ Copy 1Hz/1 min derived adjusted channels from miniseed to legacy """
+ input_factory = input_factory or get_miniseed_factory(data_type="adjusted")
+ output_factory = output_factory or get_edge_factory(data_type="adjusted")
+ copy_channels(
+ observatory=observatory,
+ channels=[
+ ("X", "X"),
+ ("Y", "Y"),
+ ("Z", "Z"),
+ ("F", "F"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="second",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ copy_channels(
+ observatory=observatory,
+ channels=[
+ ("X", "X"),
+ ("Y", "Y"),
+ ("Z", "Z"),
+ ("F", "F"),
+ ("H_SQ", "MSQ"),
+ ("H_SV", "MSV"),
+ ("H_Dist", "MDT"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="minute",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+
+
+def legacy_variation(
+ observatory: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[MiniSeedFactory] = None,
+ output_factory: Optional[EdgeFactory] = None,
+):
+ """ Copy 1Hz/1 min derived variation channels from miniseed to legacy """
+ input_factory = input_factory or get_miniseed_factory(data_type="variation")
+ output_factory = output_factory or get_edge_factory(data_type="variation")
+ copy_channels(
+ observatory=observatory,
+ channels=[
+ ("H_SQ", "MSQ"),
+ ("H_SV", "MSV"),
+ ("H_Dist", "MDT"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="minute",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+
+
+def prepare(
+ observatory: str,
+ starttime: UTCDateTime,
+ endtime: UTCDateTime,
+ update_limit: int = 10,
+ input_factory: Optional[EdgeFactory] = None,
+ output_factory: Optional[MiniSeedFactory] = None,
+):
+ """Copy 1Hz legacy HEZF to 1Hz miniseed UVWF
+ Copy 1 min T1-4 from legacy to miniseed
+ """
+ copy_channels(
+ observatory=observatory,
+ channels=[
+ ("H", "U"),
+ ("E", "V"),
+ ("Z", "W"),
+ ("F", "F"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="second",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
+ copy_channels(
+ observatory=observatory,
+ channels=[
+ ("UK1", "UK1"),
+ ("UK2", "UK2"),
+ ("UK3", "UK3"),
+ ("UK4", "UK4"),
+ ],
+ starttime=starttime,
+ endtime=endtime,
+ data_interval="minute",
+ update_limit=update_limit,
+ input_factory=input_factory,
+ output_factory=output_factory,
+ )
diff --git a/setup.py b/setup.py
index 0d7801818eb84eaf21ae6c77c8b308409cf748ed..2413e901186984bf715fcfd399b0677291bb2b84 100644
--- a/setup.py
+++ b/setup.py
@@ -28,7 +28,7 @@ setuptools.setup(
"generate-matrix=geomagio.processing.affine_matrix:main",
"geomag-metadata=geomagio.metadata.main:main",
"magproc-prepfiles=geomagio.processing.magproc:main",
- "obsrio-filter=geomagio.processing.obsrio:main",
+ "observatory=geomagio.processing.observatory:main",
],
},
)
diff --git a/test/ChannelConverter_test.py b/test/ChannelConverter_test.py
index bf050888bd9229c4fd2e311f81463975af8a1591..79ae68a2b6add67f32dd134c3eedb99381e35e77 100644
--- a/test/ChannelConverter_test.py
+++ b/test/ChannelConverter_test.py
@@ -21,9 +21,9 @@ class ChannelConverterTest:
_____
Observatory frame of reference.
- h: the component corresponding to the field strength along the
+ u: the component corresponding to the field strength along the
observatories primary horizontal axis.
- e: the component corresponding to the field strength along the
+ v: the component corresponding to the field strength along the
observatories secondary horizonal axis.
d: the angle between observatory h and the horizontal field value
(aka magnetic north in the horizontal plane.)
@@ -50,59 +50,59 @@ class ChannelConverterTest:
Computed f is the combination of the 3 vector components from
either the geographic coordinate system (X,Y,Z), or the observatory
- coordinate system (h,e,z).
+ coordinate system (u,v,w).
"""
- h = 2
- e = 2
- z = 2
- fv = channel.get_computed_f_using_squares(h, e, z)
+ u = 2
+ v = 2
+ w = 2
+ fv = channel.get_computed_f_using_squares(u, v, w)
fs = math.sqrt(12)
assert_almost_equal(fv, fs, 8, "Expect fv to almost equal sqrt(12)", True)
def test_get_geo_from_obs(self):
"""ChannelConverter_test.test_get_geo_from_obs()
- The observatory component ``h`` and ``e`` combined with the declination
+ The observatory component ``u`` and ``v`` combined with the declination
basline angle ``d0`` converts to the geographic north component
``X`` and east vector ``Y``
"""
- # 1) Call get_geo_from_obs using equal h,e values with a d0 of 0
+ # 1) Call get_geo_from_obs using equal u,v values with a d0 of 0
# the geographic values X,Y will be the same. This proves the simple
# case where the observatory is aligned with geographic north.
- h = 1
- e = 1
- (X, Y) = channel.get_geo_from_obs(h, e)
+ u = 1
+ v = 1
+ (X, Y) = channel.get_geo_from_obs(u, v)
assert_almost_equal(X, 1, 8, "Expect X to almost equal 1.", True)
assert_almost_equal(Y, 1, 8, "Expect Y to almost equal 1.", True)
- # 2) Call get_geo_from_obs using h,e values of cos(15), sin(15)
+ # 2) Call get_geo_from_obs using u,v values of cos(15), sin(15)
# (to create a d of 15 degrees) and a d0 of 15 degrees.
# X and Y will be cos(30), sin(30)
- h = cos(15 * D2R)
- e = sin(15 * D2R)
+ u = cos(15 * D2R)
+ v = sin(15 * D2R)
d0 = 15 * D2R
- (X, Y) = channel.get_geo_from_obs(h, e, d0)
+ (X, Y) = channel.get_geo_from_obs(u, v, d0)
assert_almost_equal(X, cos(30 * D2R), 8, "Expect X to equal cos(30)", True)
assert_almost_equal(Y, sin(30 * D2R), 8, "Expect Y to equal sin(30)", True)
- # 3) Call get_geo_from_obs using h,e values of 1,0 with a d0 of 315
+ # 3) Call get_geo_from_obs using u,v values of 1,0 with a d0 of 315
# degrees. The geographic X,Y will be cos(45), sin(-45)
- h = 1
- e = 0
+ u = 1
+ v = 0
d0 = 315 * D2R
- (X, Y) = channel.get_geo_from_obs(h, e, d0)
+ (X, Y) = channel.get_geo_from_obs(u, v, d0)
assert_almost_equal(X, cos(45 * D2R), 8, "Expect X to equal cos(45).", True)
assert_almost_equal(Y, sin(-45 * D2R), 8, "Expect Y to equal sin(45).", True)
- # 4) Call get_geo_from_obs using h,e values of cos_30,sin_30 and d0 of
+ # 4) Call get_geo_from_obs using u,v values of cos_30,sin_30 and d0 of
# 30 degrees. The geographic X,Y will be cos(-30), sin(-30), due to
# combined angle of the observatory declination of 30 degrees, and
# the d0 of -60 degrees.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
d0 = -60 * D2R
- (X, Y) = channel.get_geo_from_obs(h, e, d0)
+ (X, Y) = channel.get_geo_from_obs(u, v, d0)
assert_almost_equal(X, cos(-30 * D2R), 8, "Expect X to equal cos(60).", True)
assert_almost_equal(Y, sin(-30 * D2R), 8, "Expect Y to equal sin(60).", True)
@@ -164,19 +164,19 @@ class ChannelConverterTest:
def test_get_mag_from_obs(self):
"""ChannelConverter_test.test_get_geo_y_from_obs()
- ``h``, ``e`` are the primary and secondary axis of the ``H``
+ ``u``, ``v`` are the primary and secondary axis of the ``H``
vector in the horizontal plane of the magnetic field. ``d0``
is the declination baseline of the observatory frame of reference.
``D`` comes from the combination of ``d0`` and the angle produced
- from the ``h`` and ``e`` components.
+ from the ``u`` and ``v`` components.
"""
- # Call get_mag_from_obs using h,d of cos(30), sin(30) and
+ # Call get_mag_from_obs using u,v of cos(30), sin(30) and
# d0 of 15 degrees. Expect H,D to equal 1, 45.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
d0 = 15 * D2R
- H, D = channel.get_mag_from_obs(h, e, d0)
+ H, D = channel.get_mag_from_obs(u, v, d0)
assert_almost_equal(H, 1, 8, "Expect H to be 1.", True)
assert_almost_equal(D, 45 * D2R, 8, "Expect D to be 45.", True)
@@ -199,51 +199,51 @@ class ChannelConverterTest:
def test_get_mag_d_from_obs(self):
"""ChannelConverter_test.test_get_mag_d_from_obs()
- The observatory components ``h`` and ``e`` form an angle (d) with
+ The observatory components ``u`` and ``v`` form an angle (d) with
the horizontal magnetic vector. Adding d to the observatory
declination baseline ``d0`` the magnetic declination angle ``D`` can
be found.
"""
- # 1) Call get_mag_d_from_obs using h,e equal to each other.
+ # 1) Call get_mag_d_from_obs using u,v equal to each other.
# Expect D to be 45 degrees.
- h = 2
- e = 2
- D = channel.get_mag_d_from_obs(h, e)
+ u = 2
+ v = 2
+ D = channel.get_mag_d_from_obs(u, v)
assert_almost_equal(D, 45 * D2R, 8, "Expect D to be 45 degrees.", True)
- # 2) Call get_mag_d_from_obs using h,e cos(30), sin(30).
+ # 2) Call get_mag_d_from_obs using u,v cos(30), sin(30).
# Expect d of 30 degress.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
- D = channel.get_mag_d_from_obs(h, e)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
+ D = channel.get_mag_d_from_obs(u, v)
assert_almost_equal(D, 30 * D2R, 8, "Expect D to equal 30 degrees", True)
- # 3) Call get_mag_d_from_obs using h,e cos(30), sin(30),
+ # 3) Call get_mag_d_from_obs using u,v cos(30), sin(30),
# d0 = 30 degrees Expect d to be 60 degress.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
d0 = 30 * D2R
- D = channel.get_mag_d_from_obs(h, e, d0)
+ D = channel.get_mag_d_from_obs(u, v, d0)
assert_almost_equal(D, 60 * D2R, 8, "Expect D to equal 60 degrees", True)
- # 4) Call get_mag_d_from_obs using h,e cos(30), sin(30),
+ # 4) Call get_mag_d_from_obs using u,v cos(30), sin(30),
# d0 = 330 degrees Expect d of 360 degress.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
d0 = 330 * D2R
- D = channel.get_mag_d_from_obs(h, e, d0)
+ D = channel.get_mag_d_from_obs(u, v, d0)
assert_almost_equal(D, 360 * D2R, 8, "Expect D to equal 360 degrees", True)
- # 5) Call get_mag_d_from_obs using h,e cos(30), sin(30),
+ # 5) Call get_mag_d_from_obs using u,v cos(30), sin(30),
# d0 = -30 degrees Expect d of 0 degress.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
d0 = -30 * D2R
- D = channel.get_mag_d_from_obs(h, e, d0)
+ D = channel.get_mag_d_from_obs(u, v, d0)
assert_almost_equal(D, 0, 8, "Expect D to equal 0 degrees", True)
- # 6) Call get_mag_d_from_obs using h,e cos(30), -sin(30),
+ # 6) Call get_mag_d_from_obs using u,v cos(30), -sin(30),
# d0 = -30 degrees. Expect d of -60 degress.
- h = cos(30 * D2R)
- e = sin(-30 * D2R)
+ u = cos(30 * D2R)
+ v = sin(-30 * D2R)
d0 = -30 * D2R
- D = channel.get_mag_d_from_obs(h, e, d0)
+ D = channel.get_mag_d_from_obs(u, v, d0)
assert_almost_equal(D, -60 * D2R, 8, "Expect D to equal -60 degrees", True)
def test_get_mag_d_from_geo(self):
@@ -303,17 +303,17 @@ class ChannelConverterTest:
The geographic north and east components ``X`` and ``Y`` of the
magnetic field vector H, combined with the declination baseline angle
- ``d0`` can be used to produce the observatory components ``h`` and
- ```e``
+ ``d0`` can be used to produce the observatory components ``u`` and
+ ```v``
"""
# 1) Call get_obs_from_geo using equal X,Y values with a d0 of 0
- # the observatory values h,e will be the same.
+ # the observatory values u,v will be the same.
X = 1
Y = 1
- (h, e) = channel.get_obs_from_geo(X, Y)
- assert_almost_equal(h, 1.0, 8, "Expect h to be 1.", True)
- assert_almost_equal(e, 1.0, 8, "Expect e to be 1.", True)
+ (u, v) = channel.get_obs_from_geo(X, Y)
+ assert_almost_equal(u, 1.0, 8, "Expect u to be 1.", True)
+ assert_almost_equal(v, 1.0, 8, "Expect v to be 1.", True)
# 2) Call get_obs_from_geo using equal X,Y values to create a 45
# degree angle (D), with a d0 of 45/2. The observatory declination
# (d) will be 45/2, the difference between the total field angle,
@@ -321,8 +321,8 @@ class ChannelConverterTest:
X = 1
Y = 1
d0 = 22.5 * D2R
- (h, e) = channel.get_obs_from_geo(X, Y, d0)
- d = channel.get_obs_d_from_obs(h, e)
+ (u, v) = channel.get_obs_from_geo(X, Y, d0)
+ d = channel.get_obs_d_from_obs(u, v)
assert_almost_equal(d, 22.5 * D2R, 8, "Expect d to be 22.5 degrees.", True)
# 3) Call get_obs_from_geo using equal X,Y values to create a 45
# degree angle (D), with a d0 of 315 degrees. The observatory
@@ -330,21 +330,21 @@ class ChannelConverterTest:
X = 1
Y = 1
d0 = 315 * D2R
- (h, e) = channel.get_obs_from_geo(X, Y, d0)
- d = channel.get_obs_d_from_obs(h, e)
+ (u, v) = channel.get_obs_from_geo(X, Y, d0)
+ d = channel.get_obs_d_from_obs(u, v)
assert_almost_equal(d, 90 * D2R, 8, "Expect d to be 90 degrees.", True)
# 4) Call get_obs_from_geo using X,Y values of cos(60), sin(60), and
- # d0 of 30 degrees. The observatory values h,e will be cos(30)
+ # d0 of 30 degrees. The observatory values u,v will be cos(30)
# and sin(30), and the observatory declination will be 30 degrees.
# The observatory angle of 30 degrees + the d0 of 30 degrees produces
# the total declination (D) of 60 degrees.
X = cos(60 * D2R)
Y = sin(60 * D2R)
d0 = 30 * D2R
- (h, e) = channel.get_obs_from_geo(X, Y, d0)
- assert_almost_equal(h, cos(30 * D2R), 8, "Expect h to be cos(30).", True)
- assert_almost_equal(e, sin(30 * D2R), 8, "Expect e to be sin(30).", True)
- d = channel.get_obs_d_from_obs(h, e)
+ (u, v) = channel.get_obs_from_geo(X, Y, d0)
+ assert_almost_equal(u, cos(30 * D2R), 8, "Expect u to be cos(30).", True)
+ assert_almost_equal(v, sin(30 * D2R), 8, "Expect v to be sin(30).", True)
+ d = channel.get_obs_d_from_obs(u, v)
assert_almost_equal(d, 30 * D2R, 8, "Expect d to be 30 degrees.", True)
def test_get_obs_from_mag(self):
@@ -357,28 +357,28 @@ class ChannelConverterTest:
H = 1
D = -22.5 * D2R
- (h, e) = channel.get_obs_from_mag(H, D, 22.5 * D2R)
- assert_almost_equal(h, cos(45 * D2R), 8, "Expect h to be cos(45)", True)
- assert_almost_equal(e, -cos(45 * D2R), 8, "Expect e to be -cos(45).", True)
+ (u, v) = channel.get_obs_from_mag(H, D, 22.5 * D2R)
+ assert_almost_equal(u, cos(45 * D2R), 8, "Expect u to be cos(45)", True)
+ assert_almost_equal(v, -cos(45 * D2R), 8, "Expect v to be -cos(45).", True)
def test_get_obs_d_from_obs(self):
"""ChannelConverter_test.test_get_obs_d_from_obs()
- ``d`` is the angle formed by the observatory components ``h`` and
- ``e`` the primary and secondary axis of the horizontal magnetic
+ ``d`` is the angle formed by the observatory components ``u`` and
+ ``v`` the primary and secondary axis of the horizontal magnetic
field vector in the observatories frame of reference.
"""
- # 1) Call get_obs_d_from_obs usine h,e equal to cos(30), sin(30).
+ # 1) Call get_obs_d_from_obs usine u,v equal to cos(30), sin(30).
# Expect d to be 30.
- h = cos(30 * D2R)
- e = sin(30 * D2R)
- d = channel.get_obs_d_from_obs(h, e)
+ u = cos(30 * D2R)
+ v = sin(30 * D2R)
+ d = channel.get_obs_d_from_obs(u, v)
assert_almost_equal(d, 30 * D2R, 8, "Expect d to be 30 degrees.", True)
- # 2) Call get_obs_d_from_obs using h,e cos(30), -sin(30). Expect
+ # 2) Call get_obs_d_from_obs using u,v cos(30), -sin(30). Expect
# d to be 30.
- h = cos(30 * D2R)
- e = sin(-30 * D2R)
- d = channel.get_obs_d_from_obs(h, e)
+ u = cos(30 * D2R)
+ v = sin(-30 * D2R)
+ d = channel.get_obs_d_from_obs(u, v)
assert_almost_equal(d, -30 * D2R, 8, "Expect d to be 30 degrees.", True)
def test_get_obs_d_from_mag_d(self):
@@ -416,69 +416,69 @@ class ChannelConverterTest:
def test_get_obs_e_from_mag(self):
"""ChannelConverter_test.test_get_obs_e_from_mag()
- ``e`` is the secondary axis or 'east' component of the observatory
+ ``v`` is the secondary axis or 'east' component of the observatory
reference frame. Using the difference between the magnetic declination
angle ``D`` and the observatory baseline declination ``d0`` to get the
- observatory declination angle d the ``e`` component can be found
+ observatory declination angle d the ``v`` component can be found
as ``H`` * sin(d)
"""
- # 1) Call get_obs_e_from mag using H,D of 1,45. Expect e to be sin(45)
+ # 1) Call get_obs_v_from_mag using H,D of 1,45. Expect v to be sin(45)
H = 1
D = 45 * D2R
- e = channel.get_obs_e_from_mag(H, D)
- assert_almost_equal(e, sin(45 * D2R), 8, "Expect e to be sin(45).", True)
- # 2) Call get_obs_e_from_mag using H,D of 1, 30. Expect e to be sin(30)
+ v = channel.get_obs_v_from_mag(H, D)
+ assert_almost_equal(v, sin(45 * D2R), 8, "Expect v to be sin(45).", True)
+ # 2) Call get_obs_v_from_mag using H,D of 1, 30. Expect v to be sin(30)
H = 1
D = 30 * D2R
- e = channel.get_obs_e_from_mag(H, D)
- assert_almost_equal(e, sin(30 * D2R), 8, "Expect e to be sin(30).", True)
- # 3) Call get_obs_e_from_mag using H,D,d0 of 1, 15, -15. Expect e to
+ v = channel.get_obs_v_from_mag(H, D)
+ assert_almost_equal(v, sin(30 * D2R), 8, "Expect v to be sin(30).", True)
+ # 3) Call get_obs_v_from_mag using H,D,d0 of 1, 15, -15. Expect v to
# be sin(30)
H = 1
D = 15 * D2R
d0 = -15 * D2R
- e = channel.get_obs_e_from_mag(H, D, d0)
- assert_almost_equal(e, sin(30 * D2R), 8, "Expect e to be sin(30)", True)
+ v = channel.get_obs_v_from_mag(H, D, d0)
+ assert_almost_equal(v, sin(30 * D2R), 8, "Expect v to be sin(30)", True)
def test_get_obs_e_from_obs(self):
"""ChannelConverter_test.test_get_obs_e_from_obs()
- ``e`` is the seconday (east) component of the observatory vector ``h``.
- ``e`` is calculated using ``h`` * tan(``d``) where ``d`` is the
+ ``v`` is the seconday (east) component of the observatory vector ``u``.
+ ``v`` is calculated using ``u`` * tan(``d``) where ``d`` is the
declination angle of the observatory.
"""
# Call get_obs_e_from_obs using h,d of 2, 30.
- # Expect e to be 2 * tan(30)
+ # Expect v to be 2 * tan(30)
h = 2
d = 30 * D2R
- e = channel.get_obs_e_from_obs(h, d)
+ v = channel.get_obs_v_from_obs(h, d)
assert_almost_equal(
- e, 2 * tan(30 * D2R), 8, "Expect e to be 2 * tan(30).", True
+ v, 2 * tan(30 * D2R), 8, "Expect v to be 2 * tan(30).", True
)
- def test_get_obs_h_from_mag(self):
- """ChannelConverter_test.test_get_obs_h_from_mag()
+ def test_get_obs_u_from_mag(self):
+ """ChannelConverter_test.test_get_obs_u_from_mag()
- The observatories horizontal magnetic vector ``h`` is caculated from
+ The observatories horizontal magnetic vector ``u`` is caculated from
the magnetic north vector ``H`` and the observatory declination angle
``d``. ``d`` is the difference of the magnetic declination ``D`` and
the observatory baseline declination ``d0``
"""
- # 1) Call get_obs_h_from_mag using H,D 1,45. Expect h to be cos(45)
+ # 1) Call get_obs_u_from_mag using H,D 1,45. Expect u to be cos(45)
H = 1
D = 45 * D2R
- h = channel.get_obs_h_from_mag(H, D)
- assert_almost_equal(h, cos(45 * D2R), 8, "Expect h to be cos(45).", True)
- # 2) Call get_obs_h_from_mag using H,D,d0 1,30,15.
- # Expect h to be cos(15)
+ u = channel.get_obs_u_from_mag(H, D)
+ assert_almost_equal(u, cos(45 * D2R), 8, "Expect u to be cos(45).", True)
+ # 2) Call get_obs_u_from_mag using H,D,d0 1,30,15.
+ # Expect u to be cos(15)
H = 1
D = 30 * D2R
d0 = 15 * D2R
- h = channel.get_obs_h_from_mag(H, D, d0)
- assert_almost_equal(h, cos(15 * D2R), 8, "Expect h to be cos(15)", True)
+ u = channel.get_obs_u_from_mag(H, D, d0)
+ assert_almost_equal(u, cos(15 * D2R), 8, "Expect u to be cos(15)", True)
def test_geo_to_obs_to_geo(self):
"""ChannelConverter_test.test_geo_to_obs_to_geo()
@@ -488,14 +488,14 @@ class ChannelConverterTest:
get_geo_from_obs. Expect the end values to be the same
as the start values.
"""
- h_in = 20840.15
- e_in = -74.16
+ u_in = 20840.15
+ v_in = -74.16
d0 = dec_bas_rad
- (X, Y) = channel.get_geo_from_obs(h_in, e_in, d0)
- (h, e) = channel.get_obs_from_geo(X, Y, d0)
+ (X, Y) = channel.get_geo_from_obs(u_in, v_in, d0)
+ (u, v) = channel.get_obs_from_geo(X, Y, d0)
- assert_almost_equal(h, 20840.15, 8, "Expect h to = 20840.15.", True)
- assert_almost_equal(e, -74.16, 8, "Expect e to = -74.16", True)
+ assert_almost_equal(u, 20840.15, 8, "Expect u to = 20840.15.", True)
+ assert_almost_equal(v, -74.16, 8, "Expect v to = -74.16", True)
def test_get_radians_from_minutes(self):
"""ChannelConverter_test.test_get_radian_from_decimal()
diff --git a/test/StreamConverter_test.py b/test/StreamConverter_test.py
index 0c8b4207ac380f9d27efb5c526d24f9703dcec78..4421411516715a3bf73662d0dba8c0048550b451 100644
--- a/test/StreamConverter_test.py
+++ b/test/StreamConverter_test.py
@@ -63,16 +63,16 @@ def test_get_geo_from_obs():
"""StreamConverter_test.test_get_geo_from_obs()
The observatory stream containing the observatory traces
- ''h'', ''d'' or ''e'', ''z'', and ''f'' converts to the geographic
+ ''u'', ''d'' or ''v'', ''w'', and ''f'' converts to the geographic
stream containing the traces ''x'', ''y'', ''z'', and ''f''
"""
obs = obspy.core.Stream()
- # 1) Call get_geo_from_obs using equal h, e streams with a decbas of 0
+ # 1) Call get_geo_from_obs using equal u, v streams with a decbas of 0
# the geographic stream values X, Y will be the same.
- obs += __create_trace("H", [1])
- obs += __create_trace("E", [1])
- obs += __create_trace("Z", [1])
+ obs += __create_trace("U", [1])
+ obs += __create_trace("V", [1])
+ obs += __create_trace("W", [1])
obs += __create_trace("F", [1])
geo = StreamConverter.get_geo_from_obs(obs)
X = geo.select(channel="X")[0].data
@@ -81,13 +81,13 @@ def test_get_geo_from_obs():
assert_almost_equal(Y[0], 1, 9, "Expect Y to almost equal 1", True)
# 2) Call get_geo_from_obs using a decbas of 15 degrees, and streams
- # with H = [cos(15), cos(30)], and E = [sin(15), sin(30)].
+ # with U = [cos(15), cos(30)], and V = [sin(15), sin(30)].
# Expect streams of X = [cos(30), cos(45)] and Y = sin(30), sin(45)
obs = obspy.core.Stream()
DECBAS = 15 * D2I
- obs += __create_trace("H", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
- obs += __create_trace("E", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
- obs += __create_trace("Z", [1, 1], DECBAS)
+ obs += __create_trace("U", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
+ obs += __create_trace("V", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
+ obs += __create_trace("W", [1, 1], DECBAS)
obs += __create_trace("F", [1, 1], DECBAS)
geo = StreamConverter.get_geo_from_obs(obs)
X = geo.select(channel="X")[0].data
@@ -137,18 +137,18 @@ def test_get_mag_from_geo():
def test_get_mag_from_obs():
"""StreamConverter_test.test_get_mag_from_obs()
- The observatory stream containing the traces ''h'', ''e'' or ''d'',
- ''z'' and ''f''
+ The observatory stream containing the traces ''u'', ''v'' or ''d'',
+ ''w'' and ''f''
"""
obs = obspy.core.Stream()
- # Call get_mag_from_obs using a DECBAS of 15 degrees, a H stream of
- # [cos(15), cos(30)] and a E stream of [sin(15), sin(30)].
+ # Call get_mag_from_obs using a DECBAS of 15 degrees, a U stream of
+ # [cos(15), cos(30)] and a V stream of [sin(15), sin(30)].
# Expect a H stream of [1, 1] and a D stream of [30 degrees, 45 degrees]
DECBAS = 15 * D2I
- obs += __create_trace("H", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
- obs += __create_trace("E", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
- obs += __create_trace("Z", [1, 1], DECBAS)
+ obs += __create_trace("U", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
+ obs += __create_trace("V", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
+ obs += __create_trace("W", [1, 1], DECBAS)
obs += __create_trace("F", [1, 1], DECBAS)
mag = StreamConverter.get_mag_from_obs(obs)
H = mag.select(channel="H")[0].data
@@ -164,35 +164,35 @@ def test_get_obs_from_geo():
The geographic stream containing the traces ''x'', ''y'', ''z'', and
''f'' converts to the observatory stream containing the traces
- ''h'', ''d'' or ''e'', ''z'', and ''f''.
+ ''u'', ''d'' or ''v'', ''w'', and ''f''.
"""
geo = obspy.core.Stream()
# Call get_geo_from_obs using a decbas of 15, a X stream of
# [cos(30), cos(45)], and a Y stream of [sin(30), sin(45)].
- # Expect a H stream of [cos(15), cos(30)] and a
- # E stream of [sin(15), sin(30)]
+ # Expect a U stream of [cos(15), cos(30)] and a
+ # V stream of [sin(15), sin(30)]
DECBAS = 15 * D2I
geo += __create_trace("X", [cos(30 * D2R), cos(45 * D2R)], DECBAS)
geo += __create_trace("Y", [sin(30 * D2R), sin(45 * D2R)], DECBAS)
geo += __create_trace("Z", [1, 1], DECBAS)
geo += __create_trace("F", [1, 1], DECBAS)
obs = StreamConverter.get_obs_from_geo(geo, True)
- H = obs.select(channel="H")[0].data
- E = obs.select(channel="E")[0].data
+ U = obs.select(channel="U")[0].data
+ V = obs.select(channel="V")[0].data
D = obs.select(channel="D")[0].data
assert_almost_equal(
- H,
+ U,
[cos(15 * D2R), cos(30 * D2R)],
9,
- "Expect H to equal [cos(15), cos(30)]",
+ "Expect U to equal [cos(15), cos(30)]",
True,
)
assert_almost_equal(
- E,
+ V,
[sin(15 * D2R), sin(30 * D2R)],
9,
- "Expect E to equal [sin(15), sin(30)",
+ "Expect V to equal [sin(15), sin(30)",
True,
)
assert_almost_equal(
@@ -205,38 +205,38 @@ def test_get_obs_from_mag():
The magnetic stream containing the traces ''h'', ''d'', ''z'', and ''f''
converts to the observatory stream containing the traces
- ''h'', ''e'' and/or ''d'', ''z'', and ''f''
+ ''u'', ''v'' and/or ''d'', ''w'', and ''f''
"""
mag = obspy.core.Stream()
# Call get_obs_from_mag using a decbas of 15, a H stream of [1,1],
- # and a D stream of [30 degrees, 45 degrees]. Expect a H stream
+ # and a D stream of [30 degrees, 45 degrees]. Expect a U stream
# of [cos(15), cos(30)], a D stream of [30 degrees, 45 degrees],
- # and a E stream of [sin(15), sin(30)]
+ # and a V stream of [sin(15), sin(30)]
DECBAS = 15 * D2I
mag += __create_trace("H", [1, 1], DECBAS)
mag += __create_trace("D", [30 * D2R, 45 * D2R], DECBAS)
mag += __create_trace("Z", [1, 1], DECBAS)
mag += __create_trace("F", [1, 1], DECBAS)
obs = StreamConverter.get_obs_from_mag(mag, True)
- H = obs.select(channel="H")[0].data
+ U = obs.select(channel="U")[0].data
D = obs.select(channel="D")[0].data
- E = obs.select(channel="E")[0].data
+ V = obs.select(channel="V")[0].data
assert_almost_equal(
- H,
+ U,
[cos(15 * D2R), cos(30 * D2R)],
9,
- "Expect H to equal [cos(15), cos(30)",
+ "Expect U to equal [cos(15), cos(30)",
True,
)
assert_almost_equal(
D, [15 * D2R, 30 * D2R], 9, "Expect D to equal [15 degrees, 30 degrees", True
)
assert_almost_equal(
- E,
+ V,
[sin(15 * D2R), sin(30 * D2R)],
9,
- "Expect E to equal [sin(15), sin(30)",
+ "Expect V to equal [sin(15), sin(30)",
True,
)
@@ -244,19 +244,19 @@ def test_get_obs_from_mag():
def test_get_obs_from_obs():
"""StreamConverter_test.test_get_obs_from_obs()
- The observatory stream can contain either ''d'' or ''e'' depending
+ The observatory stream can contain either ''d'' or ''v'' depending
on it's source. get_obs_from_obs will return either or both as part
of the obs Stream.
"""
- # 1) Call get_obs_from_obs using a decbas of 15, a H stream of
- # [cos(15), cos(30)], and a E stream of [sin(15), sin(30)].
+ # 1) Call get_obs_from_obs using a decbas of 15, a U stream of
+ # [cos(15), cos(30)], and a V stream of [sin(15), sin(30)].
# Expect a D stream of [15 degrees, 30 degrees]
obs_e = obspy.core.Stream()
DECBAS = 15 * D2I
- obs_e += __create_trace("H", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
- obs_e += __create_trace("E", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
- obs_e += __create_trace("Z", [1, 1], DECBAS)
+ obs_e += __create_trace("U", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
+ obs_e += __create_trace("V", [sin(15 * D2R), sin(30 * D2R)], DECBAS)
+ obs_e += __create_trace("W", [1, 1], DECBAS)
obs_e += __create_trace("F", [1, 1], DECBAS)
obs_D = StreamConverter.get_obs_from_obs(obs_e, False, True)
d = obs_D.select(channel="D")[0].data
@@ -264,21 +264,21 @@ def test_get_obs_from_obs():
d, [15 * D2R, 30 * D2R], 9, "Expect D to equal [15 degrees, 30 degrees]", True
)
- # 2) Call get_obs_from_obs using a decbase of 15 degrees, a H stream of
+ # 2) Call get_obs_from_obs using a decbase of 15 degrees, a U stream of
# [cos(15), cos(30)], and a D stream of [15, 30].
# Expect a D stream of [sin(15), sin(30)]
obs_d = obspy.core.Stream()
- obs_d += __create_trace("H", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
+ obs_d += __create_trace("U", [cos(15 * D2R), cos(30 * D2R)], DECBAS)
obs_d += __create_trace("D", [15 * D2R, 30 * D2R], DECBAS)
- obs_d += __create_trace("Z", [1, 1], DECBAS)
+ obs_d += __create_trace("W", [1, 1], DECBAS)
obs_d += __create_trace("F", [1, 1], DECBAS)
- obs_E = StreamConverter.get_obs_from_obs(obs_d, True, False)
- e = obs_E.select(channel="E")[0].data
+ obs_V = StreamConverter.get_obs_from_obs(obs_d, True, False)
+ v = obs_V.select(channel="V")[0].data
assert_almost_equal(
- e,
+ v,
[sin(15 * D2R), sin(30 * D2R)],
9,
- "Expect E to equal [sin(15), sin(30)",
+ "Expect V to equal [sin(15), sin(30)",
True,
)
@@ -298,13 +298,13 @@ def test_verification_data():
DECBAS = 552.7
obs_v = obspy.core.Stream()
obs_v += __create_trace(
- "H", [20889.55, 20889.57, 20889.74, 20889.86, 20889.91, 20889.81], DECBAS
+ "U", [20889.55, 20889.57, 20889.74, 20889.86, 20889.91, 20889.81], DECBAS
)
obs_v += __create_trace(
- "E", [-21.10, -20.89, -20.72, -20.57, -20.39, -20.12], DECBAS
+ "V", [-21.10, -20.89, -20.72, -20.57, -20.39, -20.12], DECBAS
)
obs_v += __create_trace(
- "Z", [47565.29, 47565.34, 47565.39, 47565.45, 47565.51, 47565.54], DECBAS
+ "W", [47565.29, 47565.34, 47565.39, 47565.45, 47565.51, 47565.54], DECBAS
)
obs_v += __create_trace(
"F", [52485.77, 52485.84, 52485.94, 52486.06, 52486.11, 52486.10], DECBAS
diff --git a/test/algorithm_test/AdjustedAlgorithm_test.py b/test/algorithm_test/AdjustedAlgorithm_test.py
index 220d498c6f1bf6a65e59243523cd73789b6dc95e..f06b4405181222ec7730a5ed02d8ec5fd1f226ef 100644
--- a/test/algorithm_test/AdjustedAlgorithm_test.py
+++ b/test/algorithm_test/AdjustedAlgorithm_test.py
@@ -87,8 +87,8 @@ def test_process_reverse_polarity_AdjustedMatrix():
],
pier_correction=-22,
),
- inchannels=["H", "E"],
- outchannels=["H", "E"],
+ inchannels=["U", "V"],
+ outchannels=["U", "V"],
)
# load boulder May 20 files from /etc/ directory
@@ -101,7 +101,7 @@ def test_process_reverse_polarity_AdjustedMatrix():
adjusted = a.process(raw)
assert_streams_almost_equal(
- adjusted=adjusted, expected=expected, channels=["H", "E"]
+ adjusted=adjusted, expected=expected, channels=["U", "V"]
)
@@ -141,8 +141,8 @@ def test_process_reverse_polarity_statefile():
# load adjusted data transform matrix and pier correction
a = AdjustedAlgorithm(
statefile="etc/adjusted/adjbou_state_HE_.json",
- inchannels=["H", "E"],
- outchannels=["H", "E"],
+ inchannels=["U", "V"],
+ outchannels=["U", "V"],
)
# load boulder May 20 files from /etc/ directory
@@ -155,7 +155,7 @@ def test_process_reverse_polarity_statefile():
adjusted = a.process(raw)
assert_streams_almost_equal(
- adjusted=adjusted, expected=expected, channels=["H", "E"]
+ adjusted=adjusted, expected=expected, channels=["U", "V"]
)
diff --git a/test/algorithm_test/FilterAlgorithm_test.py b/test/algorithm_test/FilterAlgorithm_test.py
index b4c21d45afc0581d78d340d7310cdf6cdeb1e30a..9e16a14dd315db2b78727825e603fa0680a8cfd9 100644
--- a/test/algorithm_test/FilterAlgorithm_test.py
+++ b/test/algorithm_test/FilterAlgorithm_test.py
@@ -305,6 +305,37 @@ def test_align_trace():
assert_equal(starttime, UTCDateTime("2020-08-31T02:29:30"))
+def test_get_input_interval():
+ """algorithm_test.FilterAlgorithm_test.test_get_input_interval()
+ validates trimming of processing interval at 1Hz, 1-minute, 1-hour and 1-day
+ """
+ original_starttime = original_endtime = UTCDateTime("2021-06-30T12:12:00Z")
+ # 1Hz
+ starttime, endtime = FilterAlgorithm(
+ input_sample_period=0.1, output_sample_period=1.0
+ ).get_input_interval(start=original_starttime, end=original_endtime)
+ assert starttime == UTCDateTime("2021-06-30T12:11:53.9Z")
+ assert endtime == UTCDateTime("2021-06-30T12:12:06.1Z")
+ # 1-minute
+ starttime, endtime = FilterAlgorithm(
+ input_sample_period=1.0, output_sample_period=60.0
+ ).get_input_interval(start=original_starttime, end=original_endtime)
+ assert starttime == UTCDateTime("2021-06-30T12:11:15Z")
+ assert endtime == UTCDateTime("2021-06-30T12:12:45Z")
+ # 1-hour
+ starttime, endtime = FilterAlgorithm(
+ input_sample_period=60.0, output_sample_period=3600.0
+ ).get_input_interval(start=original_starttime, end=original_endtime)
+ assert starttime == UTCDateTime("2021-06-30T12:00:00Z")
+ assert endtime == UTCDateTime("2021-06-30T12:59:00Z")
+ # 1-day
+ starttime, endtime = FilterAlgorithm(
+ input_sample_period=60.0, output_sample_period=86400.0
+ ).get_input_interval(start=original_starttime, end=original_endtime)
+ assert starttime == UTCDateTime("2021-06-30T00:00:00Z")
+ assert endtime == UTCDateTime("2021-06-30T23:59:00Z")
+
+
def test_get_nearest__oneday_average():
"""algorithm_test.FilterAlgorithm_test.test_get_nearest__oneday_average()
Tests get_nearest_time for minute to day
diff --git a/test/algorithm_test/XYZAlgorithm_test.py b/test/algorithm_test/XYZAlgorithm_test.py
index 67be143a880a3ed11b5de59c0d711e0773c205d7..bc1edb16a3160b12512fca5f302e375a6276055b 100644
--- a/test/algorithm_test/XYZAlgorithm_test.py
+++ b/test/algorithm_test/XYZAlgorithm_test.py
@@ -13,9 +13,9 @@ def test_xyzalgorithm_process():
"""
algorithm = XYZAlgorithm("obs", "geo")
timeseries = Stream()
- timeseries += __create_trace("H", [1, 1])
- timeseries += __create_trace("E", [1, 1])
- timeseries += __create_trace("Z", [1, 1])
+ timeseries += __create_trace("U", [1, 1])
+ timeseries += __create_trace("V", [1, 1])
+ timeseries += __create_trace("W", [1, 1])
timeseries += __create_trace("F", [1, 1])
outputstream = algorithm.process(timeseries)
assert_equal(outputstream[0].stats.channel, "X")
@@ -28,7 +28,7 @@ def test_xyzalgorithm_channels():
informat/outformat during instantiation.
"""
algorithm = XYZAlgorithm("obs", "geo")
- inchannels = ["H", "E", "Z", "F"]
+ inchannels = ["U", "V", "W", "F"]
outchannels = ["X", "Y", "Z", "F"]
assert_equal(algorithm.get_input_channels(), inchannels)
assert_equal(algorithm.get_output_channels(), outchannels)
@@ -38,13 +38,13 @@ def test_xyzalgorithm_limited_channels():
"""XYZAlgorithm_test.test_xyzalgorithm_limited_channels()
confirms that only the required channels are necessary for processing
- ie. 'H' and 'E' are only needed to get 'X' and 'Y' without 'Z' or 'F'
+ ie. 'U' and 'V' are only needed to get 'X' and 'Y' without 'Z' or 'F'
"""
algorithm = XYZAlgorithm("obs", "mag")
count = 5
timeseries = Stream()
- timeseries += __create_trace("H", [2] * count)
- timeseries += __create_trace("E", [3] * count)
+ timeseries += __create_trace("U", [2] * count)
+ timeseries += __create_trace("V", [3] * count)
outstream = algorithm.process(timeseries)
ds = outstream.select(channel="D")
# there is 1 trace
@@ -66,14 +66,14 @@ def test_xyzalgorithm_uneccesary_channel_empty():
"""
algorithm = XYZAlgorithm("obs", "mag")
timeseries = Stream()
- timeseries += __create_trace("H", [1, 1])
- timeseries += __create_trace("E", [1, 1])
- timeseries += __create_trace("Z", [1, np.NaN])
+ timeseries += __create_trace("U", [1, 1])
+ timeseries += __create_trace("V", [1, 1])
+ timeseries += __create_trace("W", [1, np.NaN])
timeseries += __create_trace("F", [np.NaN, np.NaN])
outstream = algorithm.process(timeseries)
assert_equal(
- outstream.select(channel="Z")[0].data.all(),
- timeseries.select(channel="Z")[0].data.all(),
+ outstream.select(channel="W")[0].data.all(),
+ timeseries.select(channel="W")[0].data.all(),
)
assert_equal(
outstream.select(channel="F")[0].data.all(),
diff --git a/test/derived_test.py b/test/derived_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..0dbf55055cf11b21dee724397301846f9b9035f4
--- /dev/null
+++ b/test/derived_test.py
@@ -0,0 +1,84 @@
+from obspy import Stream, UTCDateTime
+
+from geomagio import derived, TimeseriesUtility
+
+
+def test_get_missing_channels():
+ """test.derived_test.test_get_missing_channels()"""
+ starttime = UTCDateTime("2021-07-01T00:00:00Z")
+ endtime = UTCDateTime("2021-07-01T00:01:00Z")
+ timeseries = Stream()
+ for channel in ["U", "V", "W", "F"]:
+ timeseries += TimeseriesUtility.create_empty_trace(
+ starttime=starttime,
+ endtime=endtime,
+ observatory="TEST",
+ channel=channel,
+ type="variation",
+ interval="second",
+ network="NT",
+ station="TEST",
+ location="R0",
+ )
+ assert derived.get_missing_channels(timeseries=timeseries) == [
+ "U",
+ "V",
+ "W",
+ "F",
+ ]
+
+
+def test_get_derived_input_channels():
+ """test.derived_test.test_get_derived_input_channels()"""
+ assert derived.get_derived_input_channels(channel="G", data_type="variation") == [
+ "U",
+ "V",
+ "W",
+ "F",
+ ]
+ assert derived.get_derived_input_channels(channel="G", data_type="adjusted") == [
+ "X",
+ "Y",
+ "Z",
+ "F",
+ ]
+ assert derived.get_derived_input_channels(channel="D", data_type="variation") == [
+ "U",
+ "V",
+ ]
+ assert derived.get_derived_input_channels(channel="X", data_type="variation") == [
+ "U",
+ "V",
+ ]
+ assert derived.get_derived_input_channels(channel="Y", data_type="variation") == [
+ "U",
+ "V",
+ ]
+ assert derived.get_derived_input_channels(channel="H", data_type="adjusted") == [
+ "X",
+ "Y",
+ ]
+ assert derived.get_derived_input_channels(channel="D", data_type="adjusted") == [
+ "X",
+ "Y",
+ ]
+
+
+def test_get_required_channels():
+ """test.derived_test.test_get_required_channels()"""
+ assert derived.get_required_channels(
+ channels=["G", "D"], data_type="variation"
+ ) == ["U", "V", "W", "F"]
+ assert derived.get_required_channels(
+ channels=["X", "Y"], data_type="variation"
+ ) == ["U", "V"]
+ assert derived.get_required_channels(channels=["G", "D"], data_type="adjusted") == [
+ "X",
+ "Y",
+ "Z",
+ "F",
+ ]
+ assert derived.get_required_channels(channels=["H", "D"], data_type="adjusted") == [
+ "X",
+ "Y",
+ ]
diff --git a/test/edge_test/LegacySNCL_test.py b/test/edge_test/LegacySNCL_test.py
index adbcc8a8388ba85a543cd075d2e3bccc2c0e5202..7e456401af2983a459ea55e0f3845bbf80b3ce79 100644
--- a/test/edge_test/LegacySNCL_test.py
+++ b/test/edge_test/LegacySNCL_test.py
@@ -35,7 +35,7 @@ def test_element():
channel="MVU",
location="R0",
).element
- == "U"
+ == "H"
)
assert (
LegacySNCL(
@@ -113,17 +113,24 @@ def test_element():
def test_get_channel():
"""edge_test.LegacySNCL_test.test_get_channel()"""
- assert get_channel(element="D", interval="second") == "SVD"
- assert get_channel(element="F", interval="minute") == "MSF"
- assert get_channel(element="H", interval="hour") == "HVH"
- assert get_channel(element="E-E", interval="day") == "DQE"
- assert get_channel(element="E-N", interval="minute") == "MQN"
- assert get_channel(element="SQ", interval="minute") == "MSQ"
- assert get_channel(element="SV", interval="minute") == "MSV"
- assert get_channel(element="UK1", interval="minute") == "UK1"
- assert get_channel(element="DIST", interval="minute") == "MDT"
- assert get_channel(element="DST", interval="minute") == "MGD"
- assert get_channel(element="UK1.R0", interval="minute") == "UK1"
+ assert get_channel(data_type="variation", element="D", interval="second") == "SVD"
+ assert get_channel(data_type="variation", element="F", interval="minute") == "MSF"
+ assert get_channel(data_type="variation", element="H", interval="hour") == "HVH"
+ assert get_channel(data_type="variation", element="E-E", interval="day") == "DQE"
+ assert get_channel(data_type="variation", element="E-N", interval="minute") == "MQN"
+ assert get_channel(data_type="variation", element="SQ", interval="minute") == "MSQ"
+ assert get_channel(data_type="variation", element="SV", interval="minute") == "MSV"
+ assert get_channel(data_type="variation", element="UK1", interval="minute") == "UK1"
+ assert (
+ get_channel(data_type="variation", element="DIST", interval="minute") == "MDT"
+ )
+ assert get_channel(data_type="variation", element="DST", interval="minute") == "MGD"
+ assert (
+ get_channel(data_type="variation", element="UK1.R0", interval="minute") == "UK1"
+ )
+ assert get_channel(data_type="variation", element="U", interval="minute") == "MVH"
+ assert get_channel(data_type="variation", element="V", interval="minute") == "MVE"
+ assert get_channel(data_type="variation", element="W", interval="minute") == "MVZ"
def test_get_location():
diff --git a/test/edge_test/MiniSeedFactory_test.py b/test/edge_test/MiniSeedFactory_test.py
index bd1367d37ad3487522ddad7f5df879c3c37ad5f5..502300d264684c273f158bb4d80c70c28d6f2b5c 100644
--- a/test/edge_test/MiniSeedFactory_test.py
+++ b/test/edge_test/MiniSeedFactory_test.py
@@ -50,7 +50,7 @@ def test__pre_process():
processed = MiniSeedInputClient(host=None)._pre_process(stream=Stream(trace))
assert len(processed) == 2
for trace in processed:
- assert trace.data.dtype == "float32"
+ assert trace.data.dtype == "float64"
stats = trace.stats
assert stats.npts == 86400
assert stats.starttime.timestamp % 86400 == 0
@@ -62,7 +62,7 @@ def test__format_miniseed():
buf = io.BytesIO()
trace = __create_trace(numpy.arange((86400 * 2) + 1), channel="H")
MiniSeedInputClient(host=None)._format_miniseed(stream=Stream(trace), buf=buf)
- block_size = 512
+ block_size = 1024
data = buf.getvalue()
n_blocks = int(len(data) / block_size)
assert n_blocks == 1516