From c20a839ee28b62f6a7c0b09101297277ce841249 Mon Sep 17 00:00:00 2001 From: spencer <swilbur@usgs.gov> Date: Fri, 22 Dec 2023 09:42:36 -0700 Subject: [PATCH] Created a new branch to contain changes to IRISFactory.py and IRISSNCL.py, the channels reported by ASL and the expected channel coversions have been changed so that LF1 and LF2 correspond to V and U. The reported azimuth and dip is also checked and changed if the polarity is reversed. --- geomagio/edge/IRISFactory.py | 95 ++++++++++++++++++++++----------- geomagio/edge/IRISSNCL.py | 8 +-- test/edge_test/IRISSNCL_test.py | 16 +++--- 3 files changed, 77 insertions(+), 42 deletions(-) diff --git a/geomagio/edge/IRISFactory.py b/geomagio/edge/IRISFactory.py index 7f7889c31..a9b01b1b9 100644 --- a/geomagio/edge/IRISFactory.py +++ b/geomagio/edge/IRISFactory.py @@ -2,7 +2,7 @@ from typing import List, Optional from obspy import Trace, UTCDateTime, Stream from obspy.clients.iris.client import Client - +from obspy.clients.fdsn import Client as FDSNClient from ..geomag_types import DataInterval, DataType from geomagio.VariometerMetadata import VariometerMetadata from .. import TimeseriesUtility @@ -79,6 +79,7 @@ class IRISFactory(MiniSeedFactory): self.base_url = base_url self.network = network self.client = Client(base_url=base_url) + self.metaClient = FDSNClient def _get_timeseries( self, @@ -150,6 +151,25 @@ class IRISFactory(MiniSeedFactory): TimeseriesUtility.pad_and_trim_trace( trace=data[0], starttime=starttime, endtime=endtime ) + # Beneath is necessary code to correct the azimuth reported for the LF1(H) and LF2(E) Channels for a Stream Object. + azimuth = self._get_azimuth_from_orientations(data[0], starttime) + + # Checks Azimuth for LF2 channel + if 270 > azimuth[0] and azimuth[0] > 90.0 and data[0].stats.channel == "LF2": + data[0].data *= -1 + + # Checks Azimuth for LF1 channel + elif azimuth[0] > 180 and azimuth[0] < 360 and data[0].stats.channel == "LF1": + data[0].data *= -1 + + # Checks Azimuth for LFZ channel + elif azimuth[0] != 0 and data[0].stats.channel == "LFZ": + data[0].data *= 0 + + # Checks Dip for LFZ channel + elif azimuth[1] > 0 and azimuth[1] < 180 and data[0].stats.channel == "LFZ": + data[0].data *= -1 + self._set_metadata(data, observatory, channel, type, interval) return data @@ -208,32 +228,47 @@ class IRISFactory(MiniSeedFactory): ################################### - # def _get_azimuth_from_orientations(self, trace: Trace, starttime: UTCDateTime) -> float: - # # Retrieve station orientation information using FDSN for each trace - # print(trace.stats.location) - # sncl = IRISSNCL.get_sncl( - # data_type=trace.stats.location, - # element=trace.stats.channel, - # interval=self.interval, - # station=trace.stats.station, - # network=trace.stats.network, - # location=trace.stats.location, - # ) - # print(sncl) - # # station = trace.stats.station - # inv = FDSN.get_stations(self.metaClient, network=sncl.network, station=sncl.station, channel=sncl.channel, level="channel") - - # if inv: - # # Construct the channel code using the current trace's information - # channel_code = f"{sncl.network}.{sncl.station}.{sncl.location}.{sncl.channel}" - - # # Get orientation for the constructed channel code and time - # print(channel_code) - # orient = inv.get_orientation(channel_code, starttime) - - # if orient: - # # Return the azimuth from orientation - # trace.stats.orientation = orient - # return orient["azimuth"] - - # return 0.0 # Default azimuth if metadata is not available + def _get_azimuth_from_orientations( + self, trace: Trace, starttime: UTCDateTime + ) -> tuple[float, float]: + # Retrieve station orientation information using FDSN for each trace + # print(trace.stats.location) + sncl = IRISSNCL.get_sncl( + data_type=trace.stats.location, + element=trace.stats.channel, + interval=self.interval, + station=trace.stats.station, + network=trace.stats.network, + location=trace.stats.location, + ) + print(sncl) + # station = trace.stats.station + + FDSN = FDSNClient("IRIS") + inv = FDSN.get_stations( + network=sncl.network, + station=sncl.station, + channel=sncl.channel, + level="channel", + ) + + if inv: + # Construct the channel code using the current trace's information + channel_code = ( + f"{sncl.network}.{sncl.station}.{sncl.location}.{sncl.channel}" + ) + + # Get orientation for the constructed channel code and time + print(channel_code) + orient = inv.get_orientation(channel_code, starttime) + + if orient: + # Return the azimuth from orientation + + azimuth = orient["azimuth"] + dip = orient["dip"] + + # Return both azimuth and dip + return azimuth, dip + + return 0.0 # Default azimuth if metadata is not available diff --git a/geomagio/edge/IRISSNCL.py b/geomagio/edge/IRISSNCL.py index 0fdcc4170..1993c755f 100644 --- a/geomagio/edge/IRISSNCL.py +++ b/geomagio/edge/IRISSNCL.py @@ -51,9 +51,9 @@ class IRISSNCL(SNCL): def _get_iris_element(channel: str) -> str: channel_end = channel[1:] if channel_end == "F1": - return "U" - elif channel_end == "F2": return "V" + elif channel_end == "F2": + return "U" elif channel_end == "FZ": return "W" raise ValueError(f"Unsupported channel: {channel}") @@ -73,9 +73,9 @@ def get_iris_channel( def _get_channel_end(element: str) -> str: if element in ["H", "U"]: - return "F1" - elif element in ["E", "V"]: return "F2" + elif element in ["E", "V"]: + return "F1" elif element in ["Z", "W"]: return "FZ" elif element[0:2] == "LF": # predefined element diff --git a/test/edge_test/IRISSNCL_test.py b/test/edge_test/IRISSNCL_test.py index b1bfec0cf..06021877a 100644 --- a/test/edge_test/IRISSNCL_test.py +++ b/test/edge_test/IRISSNCL_test.py @@ -16,11 +16,11 @@ def test_element(): """edge_test.IRISSNCL_test.test_element()""" assert ( IRISSNCL(station="ANMO", network="IU", channel="LF1", location="40").element - == "U" + == "V" ) assert ( IRISSNCL(station="ANMO", network="IU", channel="LF2", location="40").element - == "V" + == "U" ) assert ( IRISSNCL(station="ANMO", network="IU", channel="LFZ", location="40").element @@ -35,7 +35,7 @@ def test_get_iris_channel(): """edge_test.IRISSNCL_test.test_get_iris_channel()""" assert ( get_iris_channel( - element="H", + element="E", interval="second", data_type="variation", network="IU", @@ -45,7 +45,7 @@ def test_get_iris_channel(): ) assert ( get_iris_channel( - element="U", + element="V", interval="second", data_type="variation", network="IU", @@ -55,7 +55,7 @@ def test_get_iris_channel(): ) assert ( get_iris_channel( - element="E", + element="H", interval="second", data_type="variation", network="IU", @@ -65,7 +65,7 @@ def test_get_iris_channel(): ) assert ( get_iris_channel( - element="V", + element="U", interval="second", data_type="variation", network="IU", @@ -114,7 +114,7 @@ def test_get_sncl(): station="ANMO", network="IU", location="40", - ) == IRISSNCL(station="ANMO", network="IU", channel="LF1", location="40") + ) == IRISSNCL(station="ANMO", network="IU", channel="LF2", location="40") with pytest.raises(ValueError) as error: IRISSNCL.get_sncl( data_type="adjusted", @@ -154,6 +154,6 @@ def test_parse_sncl(): "station": "ANMO", "network": "IU", "data_type": "variation", - "element": "U", + "element": "V", "interval": "second", } -- GitLab