diff --git a/geomagio/edge/FDSNFactory.py b/geomagio/edge/FDSNFactory.py index f682ff4868d454c85014960a32dccbe39eb90b53..fc27f1052d760073915826947d498833964527a2 100644 --- a/geomagio/edge/FDSNFactory.py +++ b/geomagio/edge/FDSNFactory.py @@ -244,10 +244,7 @@ class FDSNFactory(TimeseriesFactory): ) except FDSNNoDataException: data = Stream() - print("The data count is", data.count) - # make sure data is 32bit int - for trace in data: - trace.data = trace.data.astype("i4") + data.merge() if data.count() == 0 and add_empty_channels: data += self._get_empty_trace( @@ -264,29 +261,12 @@ class FDSNFactory(TimeseriesFactory): TimeseriesUtility.pad_and_trim_trace( trace=data[0], starttime=starttime, endtime=endtime ) - # Beneath is necessary code to check the reported azimuth - # for the LF1 (E or Y) and LF2 (H or X) channels and determine if intstrument axes have been reversed. - azi, dip = FDSNFactory._get_orientations(self, data[0], starttime, sncl) - - # Checks Azimuth for LF2 channel - if 270 > azi and azi > 90.0 and data[0].stats.channel[-2:] == "F2": - data[0].data *= -1 - - # Checks Azimuth for LF1 channel - elif azi > 180 and azi < 360 and data[0].stats.channel[-2:] == "F1": - data[0].data *= -1 - - # Checks Dip for LFZ channel - elif dip > 0 and dip < 180 and data[0].stats.channel[-2:] == "FZ": - data[0].data *= -1 - - # Remove channel Response: - # The function "remove_response" appears to be doing what we want; - # i.e. convert from counts to NT, but this may be a placeholder - # at least until we see how this function behaves if - # a station has a frequency response. - if data.count != 0: - data.remove_response(output="DEF", zero_mean=False, taper=False) + + # Rotate the trace into a right handed coordinate frame. + # This will worrk assuming the metadata is reported correctly. + data = self._rotate_and_return_requested_channel( + sncl, starttime, endtime, channel + ) self._set_metadata(data, observatory, channel, type, interval) @@ -358,30 +338,63 @@ class FDSNFactory(TimeseriesFactory): trace.stats, observatory, channel, type, interval ) - def _get_orientations( - self, trace: Trace, starttime: UTCDateTime, sncl - ) -> Tuple[float, float]: - # Retrieve station orientation information using FDSN for each trace - - FDSN = FDSNClient("IRIS") - inv = FDSN.get_stations( - network=sncl.network, - station=sncl.station, - channel=sncl.channel, - level="channel", - ) - - # Construct the channel code using the current trace's information - channel_code = f"{sncl.network}.{sncl.station}.{sncl.location}.{sncl.channel}" + def _rotate_and_return_requested_channel( + self, sncl, starttime: UTCDateTime, endtime: UTCDateTime, requested_channel: str + ) -> Stream: + """ + Retrieves the necessary LF channels, rotates them to ZNE, and returns the requested channel. + Channels are mapped as follows: 'X' -> 'LF2', 'Y' -> 'LF1', 'Z' -> 'LFZ'. + """ + # Mapping X, Y, Z to LF channels + channel_map = {"X": "LF2", "Y": "LF1", "Z": "LFZ"} + # Initialize FDSN client and get the necessary data + FDSN = self.Client + # Determine if the requested channel is X, Y, or Z + if requested_channel in channel_map: + # Pull the LF* channels needed for rotation + lf_channels = ["LF1", "LF2", "LFZ"] + + inv = FDSN.get_stations( + network=sncl.network, + station=sncl.station, + location=sncl.location, # Use location if necessary + channel=",".join(lf_channels), # Request LF1, LF2, and LFZ + starttime=starttime, + endtime=endtime, + level="response", + ) - # Get orientation for the constructed channel code and time - orient = inv.get_orientation(channel_code, starttime) + # # Get the waveform data for LF* channels + st = FDSN.get_waveforms( + network=sncl.network, + station=sncl.station, + location=sncl.location, + channel="LF*", + starttime=starttime, + endtime=endtime, + ) - azimuth = orient["azimuth"] - dip = orient["dip"] + # Rotate the stream to ZNE + st.rotate(method="->ZNE", inventory=inv) + print(st) - # Set azimuth in trace.stats.azimuth - trace.stats.azimuth = azimuth + # Now return only the requested channel (mapped to Z, N, or E) + if requested_channel == "X": + return st.select(channel="LFN") # N after rotation + elif requested_channel == "Y": + return st.select(channel="LFE") # E after rotation + elif requested_channel == "Z": + return st.select(channel="LFZ") # Z remains Z + print(st) - # Return both azimuth and dip - return azimuth, dip + else: + # If the requested channel is not X, Y, or Z, just retrieve and return that specific channel + st = FDSN.get_waveforms( + network=sncl.network, + station=sncl.station, + location=sncl.location, + channel=sncl.channel, + starttime=starttime, + endtime=endtime, + ) + return st