Coverage for C: \ Users \ peaco \ OneDrive \ Documents \ GitHub \ mth5 \ mth5 \ timeseries \ channel_ts.py: 78%
753 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-27 20:09 -0800
1# -*- coding: utf-8 -*-
2"""
3Channel time series module for MT data.
5This module provides the `ChannelTS` class for handling magnetotelluric (MT)
6time series data with comprehensive metadata management, calibration,
7and signal processing capabilities.
9Notes
10-----
11- Time series are stored in `xarray.DataArray` for efficient operations.
12- Metadata follows the mt_metadata standard with Survey/Station/Run/Channel hierarchy.
13- Supports instrument response removal, resampling, merging, and Obspy integration.
15"""
17from __future__ import annotations
19# ==============================================================================
20# Imports
21# ==============================================================================
22import inspect
23from typing import Any
25import mt_metadata.timeseries as metadata
26import numpy as np
27import pandas as pd
28import scipy
29import xarray as xr
30from loguru import logger
31from mt_metadata.common.list_dict import ListDict
32from mt_metadata.common.mttime import MTime
33from mt_metadata.common.units import get_unit_object
34from mt_metadata.timeseries.filters import ChannelResponse
35from obspy.core import Trace
36from scipy import signal
38from mth5.timeseries.ts_filters import RemoveInstrumentResponse
39from mth5.timeseries.ts_helpers import get_decimation_sample_rates, make_dt_coordinates
40from mth5.utils import fdsn_tools
43# =============================================================================
44# make a dictionary of available metadata classes
45# =============================================================================
46meta_classes = dict(inspect.getmembers(metadata, inspect.isclass))
49# ==============================================================================
50# Channel Time Series Object
51# ==============================================================================
52class ChannelTS:
53 """
54 Time series container for a single MT channel with full metadata.
56 Stores equally-spaced time series data in an `xarray.DataArray` with
57 a time coordinate index. Integrates comprehensive metadata from
58 Survey/Station/Run/Channel hierarchy and supports calibration,
59 resampling, merging, and format conversions.
61 Parameters
62 ----------
63 channel_type : {'electric', 'magnetic', 'auxiliary'}, default 'auxiliary'
64 Type of the channel.
65 data : array-like, optional
66 Time series data (numpy array, pandas DataFrame/Series, xarray.DataArray).
67 channel_metadata : mt_metadata.timeseries.Electric | Magnetic | Auxiliary | dict, optional
68 Channel-specific metadata.
69 station_metadata : mt_metadata.timeseries.Station | dict, optional
70 Station metadata.
71 run_metadata : mt_metadata.timeseries.Run | dict, optional
72 Run metadata.
73 survey_metadata : mt_metadata.timeseries.Survey | dict, optional
74 Survey metadata.
75 **kwargs
76 Additional attributes to set on the object.
78 Attributes
79 ----------
80 ts : numpy.ndarray
81 The time series data array.
82 sample_rate : float
83 Sample rate in samples per second.
84 start : MTime
85 Start time (UTC).
86 end : MTime
87 End time (UTC), derived from start + duration.
88 n_samples : int
89 Number of samples.
90 component : str
91 Component name (e.g., 'ex', 'hy', 'temperature').
92 channel_response : ChannelResponse
93 Full instrument response filter chain.
95 Notes
96 -----
97 - End time is a derived property and cannot be set directly.
98 - Leverages xarray for efficient interpolation, resampling, and groupby operations.
99 - Metadata follows mt_metadata standards with automatic time period updates.
101 Examples
102 --------
103 Create an auxiliary channel with synthetic data::
105 >>> from mth5.timeseries import ChannelTS
106 >>> import numpy as np
107 >>> ts_obj = ChannelTS('auxiliary')
108 >>> ts_obj.sample_rate = 8
109 >>> ts_obj.start = '2020-01-01T12:00:00+00:00'
110 >>> ts_obj.ts = np.random.randn(4096)
111 >>> ts_obj.station_metadata.id = 'MT001'
112 >>> ts_obj.run_metadata.id = 'MT001a'
113 >>> ts_obj.component = 'temperature'
114 >>> print(ts_obj)
116 Calibrate and remove instrument response::
118 >>> calibrated = ts_obj.remove_instrument_response()
119 >>> calibrated.channel_metadata.units
120 """
122 def __init__(
123 self,
124 channel_type: str = "auxiliary",
125 data: (
126 np.ndarray | pd.DataFrame | pd.Series | xr.DataArray | list | tuple | None
127 ) = None,
128 channel_metadata: (
129 metadata.Electric | metadata.Magnetic | metadata.Auxiliary | dict | None
130 ) = None,
131 station_metadata: metadata.Station | dict | None = None,
132 run_metadata: metadata.Run | dict | None = None,
133 survey_metadata: metadata.Survey | dict | None = None,
134 **kwargs: Any,
135 ) -> None:
136 self.logger = logger
138 self._channel_type = self._validate_channel_type(channel_type)
139 self._survey_metadata = self._initialize_metadata()
141 self.data_array = xr.DataArray([1], coords=[("time", [1])], name="ts")
142 self._channel_response = ChannelResponse() # type: ignore
144 self.survey_metadata = survey_metadata
145 self.station_metadata = station_metadata
146 self.run_metadata = run_metadata
147 self.channel_metadata = channel_metadata
148 self._sample_rate = self.get_sample_rate_supplied_at_init(channel_metadata)
149 # input data
150 if data is not None:
151 self.ts = data
152 else:
153 self._update_xarray_metadata()
155 for key in list(kwargs.keys()):
156 setattr(self, key, kwargs[key])
158 def get_sample_rate_supplied_at_init(
159 self,
160 channel_metadata: (
161 metadata.Electric | metadata.Magnetic | metadata.Auxiliary | dict | None
162 ),
163 ) -> float | None:
164 """
165 Extract sample_rate from channel_metadata if available.
167 Parameters
168 ----------
169 channel_metadata : mt_metadata.timeseries.Electric | Magnetic | Auxiliary | dict | None
170 Metadata that may contain a sample_rate field.
172 Returns
173 -------
174 float | None
175 Sample rate if found, otherwise None.
177 Notes
178 -----
179 Supports nested dict structures like ``{"electric": {"sample_rate": 8.0}}``.
182 """
183 sr = None
184 if channel_metadata is None:
185 sr = None
186 elif isinstance(channel_metadata, dict):
187 # check first two layers for sample_rate key
188 if "sample_rate" in channel_metadata.keys():
189 sr = channel_metadata["sample_rate"]
190 else:
191 for k, v in channel_metadata.items():
192 if isinstance(v, dict):
193 if "sample_rate" in v.keys():
194 sr = v["sample_rate"]
195 else:
196 try:
197 # if an mt_metadata.timeseries access attr
198 sr = channel_metadata.sample_rate
199 except AttributeError:
200 sr = None
201 return sr
203 def __str__(self) -> str:
204 """
205 Return a summary string representation of the channel.
207 Returns
208 -------
209 str
210 Multi-line summary including survey, station, run, component,
211 sample rate, time range, and sample count.
212 """
213 lines = [
214 f"Survey: {self.survey_metadata.id}",
215 f"Station: {self.station_metadata.id}",
216 f"Run: {self.run_metadata.id}",
217 f"Channel Type: {self.channel_type}",
218 f"Component: {self.component}",
219 f"Sample Rate: {self.sample_rate}",
220 f"Start: {self.start}",
221 f"End: {self.end}",
222 f"N Samples: {self.n_samples}",
223 ]
225 return "\n\t".join(["Channel Summary:"] + lines)
227 def __repr__(self) -> str:
228 return self.__str__()
230 def __eq__(self, other: object) -> bool:
231 """
232 Test equality with another ChannelTS.
234 Parameters
235 ----------
236 other : object
237 Object to compare with.
239 Returns
240 -------
241 bool
242 True if metadata and data arrays are equal.
244 Raises
245 ------
246 TypeError
247 If `other` is not a ChannelTS instance.
248 """
249 if not isinstance(other, ChannelTS):
250 raise TypeError(f"Cannot compare ChannelTS with {type(other)}.")
251 if not other.channel_metadata == self.channel_metadata:
252 return False
253 if self.data_array.equals(other.data_array) is False:
254 msg = "timeseries are not equal"
255 self.logger.info(msg)
256 return False
257 return True
259 def __ne__(self, other: object) -> bool:
260 return not self.__eq__(other)
262 def __lt__(self, other: ChannelTS) -> bool:
263 """
264 Compare start times of two channels.
266 Parameters
267 ----------
268 other : ChannelTS
269 Channel to compare with.
271 Returns
272 -------
273 bool
274 True if self.start < other.start and sample rates match.
276 Raises
277 ------
278 TypeError
279 If `other` is not a ChannelTS instance.
280 """
281 if not isinstance(other, ChannelTS):
282 raise TypeError(f"Cannot compare ChannelTS with {type(other)}")
283 self.logger.info("Only testing start time")
284 return self.start < other.start
286 def __gt__(self, other: ChannelTS) -> bool:
287 if not isinstance(other, ChannelTS):
288 raise TypeError(f"Cannot compare ChannelTS with {type(other)}")
289 return self.start > other.start
291 def __add__(self, other: ChannelTS) -> ChannelTS:
292 """
293 Combine two channels with the same component.
295 Combines using `xr.combine_by_coords`, computes a monotonic time index,
296 and reindexes with linear interpolation.
298 Parameters
299 ----------
300 other : ChannelTS
301 Channel to combine with this one.
303 Returns
304 -------
305 ChannelTS
306 Combined channel with monotonic time index.
308 Raises
309 ------
310 TypeError
311 If `other` is not a ChannelTS.
312 ValueError
313 If components differ.
315 Examples
316 --------
317 Merge two sequential segments::
319 >>> combined = ch1 + ch2
320 """
321 if not isinstance(other, ChannelTS):
322 raise TypeError(f"Cannot combine {type(other)} with ChannelTS.")
323 if self.component != other.component:
324 raise ValueError(
325 "Cannot combine channels with different components. "
326 f"{self.component} != {other.component}"
327 )
328 if self.data_array.name != self.component:
329 self.data_array.name = self.component
330 if other.data_array.name != self.component:
331 other.data_array.name = self.component
332 # combine into a data set use override to keep attrs from original
333 combined_ds = xr.combine_by_coords(
334 [self.data_array, other.data_array], combine_attrs="override"
335 )
337 # Handle datetime.timedelta for Python 3.12+ compatibility
338 duration = combined_ds.time.max().values - combined_ds.time.min().values
339 if hasattr(duration, "total_seconds"):
340 # Python datetime.timedelta
341 duration_ns = duration.total_seconds() * 1e9
342 elif hasattr(duration, "view"):
343 # numpy timedelta64 - convert to nanoseconds
344 duration_ns = float(duration / np.timedelta64(1, "ns"))
345 else:
346 # Already numeric
347 duration_ns = float(duration)
349 n_samples = (self.sample_rate * duration_ns / 1e9) + 1
351 new_dt_index = make_dt_coordinates(
352 combined_ds.time.min().values, self.sample_rate, n_samples
353 )
355 new_channel = ChannelTS(
356 channel_type=self.channel_metadata.type,
357 channel_metadata=self.channel_metadata,
358 run_metadata=self.run_metadata,
359 station_metadata=self.station_metadata,
360 survey_metadata=self.survey_metadata,
361 channel_response=self.channel_response,
362 )
364 new_channel.data_array = combined_ds.interp(
365 time=new_dt_index, method="slinear"
366 ).to_array()
368 new_channel.channel_metadata.time_period.start = new_channel.start
369 new_channel.channel_metadata.time_period.end = new_channel.end
371 new_channel.run_metadata.update_time_period()
372 new_channel.station_metadata.update_time_period()
373 new_channel.survey_metadata.update_time_period()
375 new_channel._update_xarray_metadata()
377 return new_channel
379 def _initialize_metadata(self) -> metadata.Survey:
380 """
381 Create a Survey metadata hierarchy with default Station/Run/Channel.
383 Returns
384 -------
385 mt_metadata.timeseries.Survey
386 Initialized survey metadata with default IDs.
387 """
389 survey_metadata = metadata.Survey(id="0")
390 survey_metadata.stations.append(metadata.Station(id="0"))
391 survey_metadata.stations[0].runs.append(metadata.Run(id="0"))
393 # Create temporary channel metadata with valid default components
394 channel_type_lower = self.channel_type.lower()
395 if channel_type_lower == "electric":
396 ch_metadata = meta_classes[self.channel_type](component="ex")
397 elif channel_type_lower == "magnetic":
398 ch_metadata = meta_classes[self.channel_type](component="hx")
399 elif channel_type_lower == "auxiliary":
400 ch_metadata = meta_classes[self.channel_type](component="temperature")
401 else:
402 # Fallback for unknown types - try with a generic component
403 ch_metadata = meta_classes[self.channel_type](component="temp")
405 ch_metadata.type = self.channel_type.lower()
406 survey_metadata.stations[0].runs[0].channels.append(ch_metadata)
408 return survey_metadata
410 def _validate_channel_type(self, channel_type: str | None) -> str:
411 """
412 Validate and normalize channel type.
414 Parameters
415 ----------
416 channel_type : str | None
417 Channel type string.
419 Returns
420 -------
421 str
422 Capitalized valid channel type: 'Electric', 'Magnetic', or 'Auxiliary'.
424 Raises
425 ------
426 ValueError
427 If channel type is not recognized.
428 """
430 if channel_type is None:
431 channel_type = "auxiliary"
432 if channel_type.lower() not in ["electric", "magnetic"]:
433 channel_type = "auxiliary"
434 if not channel_type.capitalize() in meta_classes.keys():
435 msg = (
436 "Channel type is undefined, must be [ electric | "
437 "magnetic | auxiliary ]"
438 )
439 self.logger.error(msg)
440 raise ValueError(msg)
441 return channel_type.capitalize()
443 def _validate_channel_metadata(
444 self,
445 channel_metadata: (
446 metadata.Electric | metadata.Magnetic | metadata.Auxiliary | dict
447 ),
448 ) -> metadata.Electric | metadata.Magnetic | metadata.Auxiliary:
449 """
450 Validate and normalize channel metadata input.
452 Parameters
453 ----------
454 channel_metadata : mt_metadata.timeseries.Electric | Magnetic | Auxiliary | dict
455 Metadata to validate.
457 Returns
458 -------
459 mt_metadata.timeseries.Electric | Magnetic | Auxiliary
460 Validated metadata object.
462 Raises
463 ------
464 TypeError
465 If input is not an expected type.
466 """
467 expected_types = (
468 metadata.Electric,
469 metadata.Magnetic,
470 metadata.Auxiliary,
471 )
472 if isinstance(channel_metadata, expected_types):
473 return channel_metadata.copy()
475 if not isinstance(channel_metadata, dict):
476 msg = (
477 f"input metadata must be type {type(self.channel_metadata)}"
478 f" or dict, not {type(channel_metadata)}"
479 )
480 self.logger.error(msg)
481 raise TypeError(msg)
483 channel_metadata_lower_keys = [x.lower() for x in channel_metadata.keys()]
484 if self.channel_type.lower() not in channel_metadata_lower_keys:
485 try:
486 self.channel_type = channel_metadata["type"]
487 except KeyError:
488 pass
489 channel_metadata = {self.channel_type: channel_metadata}
490 self.channel_type = list(channel_metadata.keys())[0]
491 # Create channel metadata with proper default component
492 channel_type_lower = self.channel_type.lower()
493 if channel_type_lower == "electric":
494 ch_metadata = meta_classes[self.channel_type]()
495 elif channel_type_lower == "magnetic":
496 ch_metadata = meta_classes[self.channel_type]()
497 elif channel_type_lower == "auxiliary":
498 ch_metadata = meta_classes[self.channel_type]()
499 else:
500 ch_metadata = meta_classes[self.channel_type]()
501 self.logger.debug("Loading from metadata dict")
502 ch_metadata.from_dict(channel_metadata)
503 return ch_metadata
505 def _validate_run_metadata(self, run_metadata: metadata.Run | dict) -> metadata.Run:
506 """
507 Validate and normalize run metadata input.
509 Parameters
510 ----------
511 run_metadata : mt_metadata.timeseries.Run | dict
512 Run metadata to validate.
514 Returns
515 -------
516 mt_metadata.timeseries.Run
517 Validated run metadata object.
519 Raises
520 ------
521 TypeError
522 If input is not a Run object or dict.
523 """
525 if not isinstance(run_metadata, metadata.Run):
526 if isinstance(run_metadata, dict):
527 if "run" not in [cc.lower() for cc in run_metadata.keys()]:
528 run_metadata = {"Run": run_metadata}
529 r_metadata = metadata.Run()
530 r_metadata.from_dict(run_metadata)
531 self.logger.debug("Loading from metadata dict")
532 return r_metadata
533 else:
534 msg = (
535 f"input metadata must be type {type(self.run_metadata)} "
536 f"or dict, not {type(run_metadata)}"
537 )
538 self.logger.error(msg)
539 raise TypeError(msg)
540 return run_metadata.copy()
542 def _validate_station_metadata(
543 self, station_metadata: metadata.Station | dict
544 ) -> metadata.Station:
545 """
546 Validate and normalize station metadata input.
548 Parameters
549 ----------
550 station_metadata : mt_metadata.timeseries.Station | dict
551 Station metadata to validate.
553 Returns
554 -------
555 mt_metadata.timeseries.Station
556 Validated station metadata object.
558 Raises
559 ------
560 TypeError
561 If input is not a Station object or dict.
562 """
563 if not isinstance(station_metadata, metadata.Station):
564 if isinstance(station_metadata, dict):
565 if "station" not in [cc.lower() for cc in station_metadata.keys()]:
566 station_metadata = {"Station": station_metadata}
567 st_metadata = metadata.Station()
568 st_metadata.from_dict(station_metadata)
569 self.logger.debug("Loading from metadata dict")
570 return st_metadata
571 else:
572 msg = (
573 f"input metadata must be type {type(self.station_metadata)}"
574 f" or dict, not {type(station_metadata)}"
575 )
576 self.logger.error(msg)
577 raise TypeError(msg)
578 return station_metadata.copy()
580 def _validate_survey_metadata(
581 self, survey_metadata: metadata.Survey | dict
582 ) -> metadata.Survey:
583 """
584 Validate and normalize survey metadata input.
586 Parameters
587 ----------
588 survey_metadata : mt_metadata.timeseries.Survey | dict
589 Survey metadata to validate.
591 Returns
592 -------
593 mt_metadata.timeseries.Survey
594 Validated survey metadata object.
596 Raises
597 ------
598 TypeError
599 If input is not a Survey object or dict.
600 """
601 if not isinstance(survey_metadata, metadata.Survey):
602 if isinstance(survey_metadata, dict):
603 if "survey" not in [cc.lower() for cc in survey_metadata.keys()]:
604 survey_metadata = {"Survey": survey_metadata}
605 sv_metadata = metadata.Survey()
606 sv_metadata.from_dict(survey_metadata)
607 self.logger.debug("Loading from metadata dict")
608 return sv_metadata
609 else:
610 msg = (
611 f"input metadata must be type {type(self.survey_metadata)}"
612 f" or dict, not {type(survey_metadata)}"
613 )
614 self.logger.error(msg)
615 raise TypeError(msg)
616 return survey_metadata.copy()
618 def copy(self, data: bool = True) -> ChannelTS:
619 """
620 Create a copy of the ChannelTS object.
622 Parameters
623 ----------
624 data : bool, default True
625 Include data in the copy (True) or only metadata (False).
627 Returns
628 -------
629 ChannelTS
630 Copy of the channel.
632 Examples
633 --------
634 Copy metadata structure without data::
636 >>> ch_copy = ts_obj.copy(data=False)
637 """
639 if not data:
640 return ChannelTS(
641 channel_type=self.channel_metadata.type,
642 channel_metadata=self.channel_metadata.copy(),
643 run_metadata=self.run_metadata.copy(),
644 station_metadata=self.station_metadata.copy(),
645 survey_metadata=self.survey_metadata.copy(),
646 channel_response=self.channel_response.copy(),
647 )
648 else:
649 return ChannelTS(
650 channel_type=self.channel_metadata.type,
651 data=self.ts,
652 channel_metadata=self.channel_metadata.copy(),
653 run_metadata=self.run_metadata.copy(),
654 station_metadata=self.station_metadata.copy(),
655 survey_metadata=self.survey_metadata.copy(),
656 channel_response=self.channel_response.copy(),
657 )
659 ### Properties ------------------------------------------------------------
660 @property
661 def survey_metadata(self) -> metadata.Survey:
662 """
663 Survey metadata.
665 Returns
666 -------
667 mt_metadata.timeseries.Survey
668 Survey metadata with updated keys.
669 """
670 self._survey_metadata.stations[0].runs.update_keys()
671 self._survey_metadata.stations.update_keys()
672 return self._survey_metadata
674 @survey_metadata.setter
675 def survey_metadata(self, survey_metadata: metadata.Survey | dict | None) -> None:
676 """
677 Set survey metadata.
679 Parameters
680 ----------
681 survey_metadata : mt_metadata.timeseries.Survey | dict | None
682 Survey metadata object or dictionary.
683 """
685 if survey_metadata is not None:
686 survey_metadata = self._validate_survey_metadata(survey_metadata)
687 self._survey_metadata.update(survey_metadata)
689 @property
690 def station_metadata(self) -> metadata.Station:
691 """
692 Station metadata.
694 Returns
695 -------
696 mt_metadata.timeseries.Station
697 Station metadata from the first station in the survey.
698 """
699 self._survey_metadata.stations.update_keys()
700 return self.survey_metadata.stations[0]
702 @station_metadata.setter
703 def station_metadata(
704 self, station_metadata: metadata.Station | dict | None
705 ) -> None:
706 """
707 Set station metadata.
709 Parameters
710 ----------
711 station_metadata : mt_metadata.timeseries.Station | dict | None
712 Station metadata to set.
713 """
715 if station_metadata is not None:
716 station_metadata = self._validate_station_metadata(station_metadata)
718 runs = ListDict()
719 if self.run_metadata.id not in ["0", 0, None]:
720 runs.append(self.run_metadata.copy())
721 runs.extend(station_metadata.runs)
722 if len(runs) == 0:
723 runs[0] = metadata.Run(id="0")
724 # be sure there is a level below
725 if len(runs[0].channels) == 0:
726 # Create channel metadata with proper default component
727 channel_type_lower = self.channel_type.lower()
728 if channel_type_lower == "electric":
729 ch_metadata = meta_classes[self.channel_type](component="ex")
730 elif channel_type_lower == "magnetic":
731 ch_metadata = meta_classes[self.channel_type](component="hx")
732 elif channel_type_lower == "auxiliary":
733 ch_metadata = meta_classes[self.channel_type](
734 component="temperature"
735 )
736 else:
737 ch_metadata = meta_classes[self.channel_type]()
738 ch_metadata.type = self.channel_type.lower()
739 runs[0].channels.append(ch_metadata)
740 stations = ListDict()
741 stations.append(station_metadata)
742 stations[0].runs = runs
744 self.survey_metadata.stations = stations
746 @property
747 def run_metadata(self) -> metadata.Run:
748 """
749 Run metadata.
751 Returns
752 -------
753 mt_metadata.timeseries.Run
754 Run metadata from the first run in the station.
755 """
757 self._survey_metadata.stations[0].runs.update_keys()
758 self._survey_metadata.stations[0].runs[0].channels.update_keys()
759 return self.survey_metadata.stations[0].runs[0]
761 @run_metadata.setter
762 def run_metadata(self, run_metadata: metadata.Run | dict | None) -> None:
763 """
764 Set run metadata.
766 Parameters
767 ----------
768 run_metadata : mt_metadata.timeseries.Run | dict | None
769 Run metadata to set.
770 """
772 # need to make sure the first index is the desired channel
773 if run_metadata is not None:
774 run_metadata = self._validate_run_metadata(run_metadata)
776 runs = ListDict()
777 runs.append(run_metadata)
778 channels = ListDict()
779 if self.component is not None:
780 key = str(self.component)
782 channels.append(self.station_metadata.runs[0].channels[key])
783 # add existing channels
784 channels.extend(self.run_metadata.channels, skip_keys=[key, "0"])
785 # add channels from input metadata
786 channels.extend(run_metadata.channels)
788 runs[0].channels = channels
789 runs.extend(self.station_metadata.runs, skip_keys=[run_metadata.id, "0"])
791 self._survey_metadata.stations[0].runs = runs
793 @property
794 def channel_metadata(
795 self,
796 ) -> metadata.Electric | metadata.Magnetic | metadata.Auxiliary:
797 """
798 Channel metadata.
800 Returns
801 -------
802 mt_metadata.timeseries.Electric | Magnetic | Auxiliary
803 Channel metadata from the first channel in the run.
804 """
805 ch_metadata = self._survey_metadata.stations[0].runs[0].channels[0]
806 if self.has_data():
807 ch_metadata.sample_rate = self.sample_rate
808 return ch_metadata
810 @channel_metadata.setter
811 def channel_metadata(
812 self,
813 channel_metadata: (
814 metadata.Electric | metadata.Magnetic | metadata.Auxiliary | dict | None
815 ),
816 ) -> None:
817 """
818 Set channel metadata.
820 Parameters
821 ----------
822 channel_metadata : mt_metadata.timeseries.Electric | Magnetic | Auxiliary | dict | None
823 Channel metadata to set.
825 Raises
826 ------
827 ValueError
828 If the channel component is None.
829 """
831 if channel_metadata is not None:
832 channel_metadata = self._validate_channel_metadata(channel_metadata)
833 if channel_metadata.component is not None:
834 channels = ListDict()
835 if channel_metadata.component in self.run_metadata.channels.keys():
836 channels.append(
837 self.run_metadata.channels[channel_metadata.component]
838 )
839 channels[0].update(channel_metadata)
840 else:
841 channels.append(channel_metadata)
842 channels.extend(
843 self.run_metadata.channels,
844 skip_keys=[channel_metadata.component, None],
845 )
847 self.run_metadata.channels = channels
848 self.channel_type = self.run_metadata.channels[0].type
849 else:
850 raise ValueError("Channel 'component' cannot be None")
852 def _check_pd_index(self, ts_arr: pd.DataFrame | pd.Series) -> pd.DatetimeIndex:
853 """
854 Check and return the time index from a pandas DataFrame or Series.
856 Parameters
857 ----------
858 ts_arr : pandas.DataFrame | pandas.Series
859 Time series data.
861 Returns
862 -------
863 pandas.DatetimeIndex
864 Time index (existing or reconstructed from start/sample_rate).
865 """
866 if isinstance(ts_arr.index[0], pd._libs.tslibs.timestamps.Timestamp):
867 return ts_arr.index
868 else:
869 return make_dt_coordinates(self.start, self.sample_rate, ts_arr.shape[0])
871 def _validate_dataframe_input(
872 self, ts_arr: pd.DataFrame
873 ) -> tuple[pd.DataFrame, pd.DatetimeIndex]:
874 """
875 Validate pandas DataFrame input.
877 Parameters
878 ----------
879 ts_arr : pandas.DataFrame
880 DataFrame containing a 'data' column.
882 Returns
883 -------
884 tuple[pandas.DataFrame, pandas.DatetimeIndex]
885 Validated DataFrame and time index.
887 Raises
888 ------
889 ValueError
890 If 'data' column is missing or has object dtype that can't convert.
891 """
892 if "data" not in ts_arr.columns:
893 msg = (
894 "Data frame needs to have a column named `data` "
895 "where the time series data is stored"
896 )
897 self.logger.error(msg)
898 raise ValueError(msg)
900 if isinstance(type(ts_arr.data.dtype), type(np.object_)):
901 try:
902 ts_arr["data"] = ts_arr.data.astype(float)
903 except ValueError:
904 raise ValueError(
905 "DataFrame dtype is 'object' and cannot convert "
906 "data to float values, check data dtype."
907 )
908 dt = self._check_pd_index(ts_arr)
909 return ts_arr, dt
911 def _validate_series_input(
912 self, ts_arr: pd.Series
913 ) -> tuple[pd.Series, pd.DatetimeIndex]:
914 """
915 Validate pandas Series input.
917 Parameters
918 ----------
919 ts_arr : pandas.Series
920 Series containing time series data.
922 Returns
923 -------
924 tuple[pandas.Series, pandas.DatetimeIndex]
925 Validated Series and time index.
927 Raises
928 ------
929 ValueError
930 If Series has object dtype that can't convert to float.
931 """
933 if isinstance(type(ts_arr.dtype), type(np.object_)):
934 try:
935 ts_arr = ts_arr.astype(float)
936 except ValueError:
937 raise ValueError(
938 "Series dtype is 'object' and cannot convert "
939 "data to float values, check data dtype."
940 )
942 dt = self._check_pd_index(ts_arr)
943 return ts_arr, dt
945 @property
946 def ts(self) -> np.ndarray:
947 """
948 Time series data as a numpy array.
950 Returns
951 -------
952 numpy.ndarray
953 The time series data.
954 """
955 return self.data_array.data
957 @ts.setter
958 def ts(
959 self,
960 ts_arr: np.ndarray | list | tuple | pd.DataFrame | pd.Series | xr.DataArray,
961 ) -> None:
962 """
963 Set the time series data.
965 Parameters
966 ----------
967 ts_arr : numpy.ndarray | list | tuple | pandas.DataFrame | pandas.Series | xarray.DataArray
968 Time series data. DataFrames must have a 'data' column.
970 Raises
971 ------
972 TypeError
973 If data type is not supported.
975 Notes
976 -----
977 - For pandas DataFrames/Series, time index is extracted or reconstructed.
978 - For xarray.DataArray, metadata is extracted from attrs.
979 """
981 if isinstance(ts_arr, (np.ndarray, list, tuple)):
982 if not isinstance(ts_arr, np.ndarray):
983 ts_arr = np.array(ts_arr)
984 # Validate an input array to make sure its 1D
985 if len(ts_arr.shape) == 2:
986 if 1 in ts_arr.shape:
987 ts_arr = ts_arr.reshape(ts_arr.size)
988 else:
989 msg = f"Input array must be 1-D array not {ts_arr.shape}"
990 self.logger.error(msg)
991 raise ValueError(msg)
992 dt = make_dt_coordinates(self.start, self.sample_rate, ts_arr.size)
993 self.data_array = xr.DataArray(
994 ts_arr, coords=[("time", dt)], name=self.component
995 )
996 self._update_xarray_metadata()
997 elif isinstance(ts_arr, pd.core.frame.DataFrame):
998 ts_arr, dt = self._validate_dataframe_input(ts_arr)
999 self.data_array = xr.DataArray(
1000 ts_arr["data"], coords=[("time", dt)], name=self.component
1001 )
1002 self._update_xarray_metadata()
1004 elif isinstance(ts_arr, pd.core.series.Series):
1005 ts_arr, dt = self._validate_series_input(ts_arr)
1006 self.data_array = xr.DataArray(
1007 ts_arr.values, coords=[("time", dt)], name=self.component
1008 )
1009 self._update_xarray_metadata()
1010 elif isinstance(ts_arr, xr.DataArray):
1011 # TODO: need to validate the input xarray
1012 self.data_array = ts_arr
1013 # need to pull out the metadata as a separate dictionary
1014 meta_dict = dict([(k, v) for k, v in ts_arr.attrs.items()])
1016 # need to get station and run metadata out
1017 survey_dict = {}
1018 station_dict = {}
1019 run_dict = {}
1021 for key in [k for k in meta_dict.keys() if "survey." in k]:
1022 survey_dict[key.split("station.")[-1]] = meta_dict.pop(key)
1023 for key in [k for k in meta_dict.keys() if "station." in k]:
1024 station_dict[key.split("station.")[-1]] = meta_dict.pop(key)
1025 for key in [k for k in meta_dict.keys() if "run." in k]:
1026 run_dict[key.split("run.")[-1]] = meta_dict.pop(key)
1027 self.channel_type = meta_dict["type"]
1028 # Create channel metadata with proper default component
1029 channel_type_lower = self.channel_type.lower()
1030 if channel_type_lower == "electric":
1031 ch_metadata = meta_classes[self.channel_type](component="ex")
1032 elif channel_type_lower == "magnetic":
1033 ch_metadata = meta_classes[self.channel_type](component="hx")
1034 elif channel_type_lower == "auxiliary":
1035 ch_metadata = meta_classes[self.channel_type](component="temperature")
1036 else:
1037 ch_metadata = meta_classes[self.channel_type]()
1038 ch_metadata.from_dict({self.channel_type: meta_dict})
1040 self.survey_metadata.from_dict({"survey": survey_dict})
1041 self.station_metadata.from_dict({"station": station_dict})
1042 self.run_metadata.from_dict({"run": run_dict})
1043 self.channel_metadata = ch_metadata
1044 # need to run this incase things are different.
1045 self._update_xarray_metadata()
1046 else:
1047 msg = (
1048 "Data type {0} not supported".format(type(ts_arr))
1049 + ", ts needs to be a numpy.ndarray, pandas DataFrame, "
1050 + "or xarray.DataArray."
1051 )
1052 raise TypeError(msg)
1054 @property
1055 def time_index(self) -> np.ndarray:
1056 """
1057 Time index as a numpy array.
1059 Returns
1060 -------
1061 numpy.ndarray
1062 Array of datetime64[ns] timestamps.
1063 """
1065 try:
1066 return self.data_array.time.to_numpy()
1067 except AttributeError:
1068 return self.data_array.time.values
1070 @property
1071 def channel_type(self) -> str:
1072 """
1073 Channel type.
1075 Returns
1076 -------
1077 str
1078 Channel type: 'Electric', 'Magnetic', or 'Auxiliary'.
1079 """
1080 return self._channel_type
1082 @channel_type.setter
1083 def channel_type(self, value: str) -> None:
1084 """change channel type means changing the metadata type"""
1086 value = self._validate_channel_type(value)
1087 if value != self._channel_type:
1088 m_dict = self.channel_metadata.to_dict(single=True)
1090 msg = (
1091 f"Changing metadata from {self.channel_type} to {value}, "
1092 "will translate any similar attributes."
1093 )
1094 # Create channel metadata with proper default component
1095 if value.lower() == "electric":
1096 channel_metadata = meta_classes[value](component="ex")
1097 elif value.lower() == "magnetic":
1098 channel_metadata = meta_classes[value](component="hx")
1099 elif value.lower() == "auxiliary":
1100 channel_metadata = meta_classes[value](component="temperature")
1101 else:
1102 channel_metadata = meta_classes[value]()
1103 self.logger.debug(msg)
1104 for key in channel_metadata.to_dict(single=True).keys():
1105 # need to skip type otherwise it keeps the same type
1106 if key in ["type"]:
1107 continue
1108 # Skip component when changing channel types to avoid validation conflicts
1109 # The new metadata already has appropriate default component for the type
1110 if key in ["component"] and self._channel_type != value:
1111 continue
1112 try:
1113 channel_metadata.update_attribute(key, m_dict[key])
1114 except KeyError:
1115 pass
1116 self._channel_type = value
1117 self.run_metadata.channels[0] = channel_metadata
1119 def _update_xarray_metadata(self) -> None:
1120 """
1121 Update xarray attrs with current metadata.
1123 Notes
1124 -----
1125 Synchronizes channel_metadata fields into data_array.attrs and adds
1126 station/run IDs for convenient access.
1127 """
1128 self.logger.debug("Updating xarray attributes")
1130 self.channel_metadata.time_period.start = self.start.iso_no_tz
1131 self.channel_metadata.time_period.end = self.end.iso_no_tz
1132 self.channel_metadata.sample_rate = self.sample_rate
1134 self.data_array.attrs.update(
1135 self.channel_metadata.to_dict()[self.channel_metadata._class_name]
1136 )
1137 # add station and run id's here, for now this is all we need but may need
1138 # more metadata down the road.
1139 self.data_array.attrs["station.id"] = self.station_metadata.id
1140 self.data_array.attrs["run.id"] = self.run_metadata.id
1141 self.data_array.name = self.component
1143 @property
1144 def component(self):
1145 """component"""
1146 return self.channel_metadata.component
1148 @component.setter
1149 def component(self, comp):
1150 """set component in metadata and carry through"""
1151 if self.channel_metadata.type == "electric":
1152 if comp[0].lower() != "e":
1153 msg = (
1154 "The current timeseries is an electric channel. "
1155 "Cannot change channel type, create a new ChannelTS object."
1156 )
1157 self.logger.error(msg)
1158 raise ValueError(msg)
1159 elif self.channel_metadata.type == "magnetic":
1160 if comp[0].lower() not in ["h", "b"]:
1161 msg = (
1162 "The current timeseries is a magnetic channel. "
1163 "Cannot change channel type, create a new ChannelTS object."
1164 )
1165 self.logger.error(msg)
1166 raise ValueError(msg)
1167 if self.channel_metadata.type == "auxiliary":
1168 if comp[0].lower() in ["e", "h", "b"]:
1169 msg = (
1170 "The current timeseries is an auxiliary channel. "
1171 "Cannot change channel type, create a new ChannelTS object."
1172 )
1173 self.logger.error(msg)
1174 raise ValueError(msg)
1175 self.channel_metadata.component = comp
1177 # need to update the keys in the list dict
1178 channels = ListDict()
1179 channels.append(self.channel_metadata)
1180 if len(self.run_metadata.channels) > 1:
1181 for ch in self.run_metadata.channels[1:]:
1182 channels.append(ch)
1183 self.run_metadata.channels = channels
1185 self._update_xarray_metadata()
1187 # --> number of samples just to make sure there is consistency
1188 @property
1189 def n_samples(self):
1190 """number of samples"""
1191 return int(self.ts.size)
1193 @n_samples.setter
1194 def n_samples(self, n_samples):
1195 """number of samples (int)"""
1196 self.logger.warning(
1197 "Cannot set the number of samples. Use `ChannelTS.resample` or `get_slice`"
1198 )
1200 def has_data(self):
1201 """
1202 check to see if there is an index in the time series
1203 """
1204 if self.data_array.data.size > 1:
1205 if isinstance(
1206 self.data_array.indexes["time"][0],
1207 pd._libs.tslibs.timestamps.Timestamp,
1208 ):
1209 return True
1210 return False
1211 else:
1212 return False
1214 def is_high_frequency(self, threshold_dt=1e-4):
1215 """
1216 Quasi hard-coded condition to check if data are logged at more than 10kHz
1217 can be parameterized in future
1218 """
1219 if (
1220 self.data_array.coords.indexes["time"][1]
1221 - self.data_array.coords.indexes["time"][0]
1222 ).total_seconds() < threshold_dt:
1223 return True
1224 else:
1225 return False
1227 def compute_sample_rate(self):
1228 """
1229 Two cases, high_frequency (HF) data and not HF data.
1231 # Original comment about the HF case:
1232 Taking the median(diff(timestamps)) is more accurate for high sample rates, the way pandas.date_range
1233 rounds nanoseconds is not consistent between samples, therefore taking the median provides better results
1234 if the time series is long this can be inefficient so test first
1236 """
1237 if self.is_high_frequency():
1238 dt_array = np.diff(self.data_array.coords.indexes["time"])
1239 best_dt, counts = scipy.stats.mode(dt_array)
1241 # Calculate total seconds of the best dt and calculate sample rate
1242 best_dt_seconds = float(best_dt) / 1e9
1243 sr = 1 / best_dt_seconds
1244 else:
1245 t_diff = (
1246 self.data_array.coords.indexes["time"][-1]
1247 - self.data_array.coords.indexes["time"][0]
1248 )
1249 sr = self.data_array.size / t_diff.total_seconds()
1250 return np.round(sr, 0)
1252 # --> sample rate
1253 @property
1254 def sample_rate(self):
1255 """sample rate in samples/second"""
1256 if self.has_data():
1257 if self._sample_rate is None:
1258 self._sample_rate = self.compute_sample_rate()
1259 return self._sample_rate
1260 else:
1261 self.logger.debug("Data has not been set yet, sample rate is from metadata")
1262 sr = self.channel_metadata.sample_rate
1263 if sr is None:
1264 sr = 0.0
1266 if sr >= 1:
1267 return np.round(sr, 0)
1268 else:
1269 return np.round(sr, 6)
1271 @sample_rate.setter
1272 def sample_rate(self, sample_rate):
1273 """
1274 sample rate in samples/second
1276 :param sample_rate: sample rate in samples per second
1277 :type sample_rate: float
1278 """
1279 if self.has_data():
1280 self.logger.warning(
1281 "Resetting sample_rate assumes same start time and "
1282 "same number of samples, resulting in new end time. "
1283 "If you want to downsample existing time series "
1284 "use the method channelTS.resample()"
1285 )
1286 self.logger.debug(
1287 f"Resetting sample rate from {self.sample_rate} to {sample_rate}"
1288 )
1289 new_dt = make_dt_coordinates(self.start, sample_rate, self.n_samples)
1290 self.data_array.coords["time"] = new_dt
1291 else:
1292 if self.channel_metadata.sample_rate not in [0.0, None]:
1293 self.logger.warning(
1294 f"Resetting ChannelTS.channel_metadata.sample_rate to {sample_rate}. "
1295 )
1296 self.channel_metadata.sample_rate = sample_rate
1297 self._sample_rate = sample_rate
1298 self._update_xarray_metadata()
1300 @property
1301 def sample_interval(self):
1302 """
1303 Sample interval = 1 / sample_rate
1305 :return: sample interval as time distance between time samples
1306 :rtype: float
1308 """
1310 if self.sample_rate != 0:
1311 return 1.0 / self.sample_rate
1312 return 0.0
1314 ## set time and set index
1315 @property
1316 def start(self):
1317 """MTime object"""
1318 if self.has_data():
1319 return MTime(
1320 time_stamp=self.data_array.coords.indexes["time"][0].isoformat()
1321 )
1322 else:
1323 self.logger.debug(
1324 "Data not set yet, pulling start time from "
1325 "metadata.time_period.start"
1326 )
1327 return MTime(time_stamp=self.channel_metadata.time_period.start)
1329 @start.setter
1330 def start(self, start_time):
1331 """
1332 start time of time series in UTC given in some format or a datetime
1333 object.
1335 Resets epoch seconds if the new value is not equivalent to previous
1336 value.
1338 Resets how the ts data frame is indexed, setting the starting time to
1339 the new start time.
1341 :param start_time: start time of time series, can be string or epoch seconds
1343 """
1345 if not isinstance(start_time, MTime):
1346 start_time = MTime(time_stamp=start_time)
1347 self.channel_metadata.time_period.start = start_time.isoformat()
1348 if self.has_data():
1349 if start_time == MTime(
1350 time_stamp=self.data_array.coords.indexes["time"][0].isoformat()
1351 ):
1352 return
1353 else:
1354 new_dt = make_dt_coordinates(
1355 start_time, self.sample_rate, self.n_samples
1356 )
1357 self.data_array.coords["time"] = new_dt
1358 # make a time series that the data can be indexed by
1359 else:
1360 self.logger.debug("No data, just updating metadata start")
1361 self._survey_metadata.stations[0].runs[0].update_time_period()
1362 self._survey_metadata.stations[0].update_time_period()
1363 self._survey_metadata.update_time_period()
1365 self._update_xarray_metadata()
1367 @property
1368 def end(self):
1369 """MTime object"""
1370 if self.has_data():
1371 return MTime(
1372 time_stamp=self.data_array.coords.indexes["time"][-1].isoformat()
1373 )
1374 else:
1375 self.logger.debug(
1376 "Data not set yet, pulling end time from metadata.time_period.end"
1377 )
1378 return MTime(time_stamp=self.channel_metadata.time_period.end)
1380 @end.setter
1381 def end(self, end_time):
1382 """
1383 end time of time series in UTC given in some format or a datetime
1384 object.
1386 Resets epoch seconds if the new value is not equivalent to previous
1387 value.
1389 Resets how the ts data frame is indexed, setting the starting time to
1390 the new start time.
1391 """
1392 self.logger.warning(
1393 "Cannot set `end`. If you want a slice, then use get_slice method"
1394 )
1396 @property
1397 def channel_response(self):
1398 """
1399 Full channel response filter
1401 :return: full channel response filter
1402 :rtype: :class:`mt_metadata.timeseries.filters.ChannelResponse`
1404 """
1406 return self._channel_response
1408 @channel_response.setter
1409 def channel_response(self, value):
1410 """
1412 :param value: channel response filter
1413 :type value: :class:`mt_metadata.timeseries.filters.`
1416 """
1417 if value is None:
1418 return
1419 if not isinstance(value, ChannelResponse):
1420 msg = (
1421 "channel response must be a "
1422 "mt_metadata.timeseries.filters.ChannelResponse object "
1423 f"not {type(value)}."
1424 )
1425 self.logger.error(msg)
1426 raise TypeError(msg)
1427 self._channel_response = value
1429 # update channel metadata
1430 if self.channel_metadata.filter_names != value.names:
1431 for ch_filter in self._channel_response.filters_list:
1432 if ch_filter.name in self.channel_metadata.filter_names:
1433 # update existing filter info
1434 existing_filter = self.channel_metadata.get_filter(ch_filter.name)
1435 existing_filter.applied = False
1436 existing_filter.stage = ch_filter.sequence_number
1437 existing_filter.comments = ch_filter.comments
1438 else:
1439 self.channel_metadata.add_filter(
1440 name=ch_filter.name,
1441 applied=False,
1442 stage=ch_filter.sequence_number,
1443 comments=ch_filter.comments,
1444 )
1446 def get_calibration_operation(self):
1447 if self.channel_response.units_out == self.channel_metadata.unit_object.name:
1448 calibration_operation = "divide"
1449 elif self.channel_response.units_in == self.channel_metadata.unit_object.name:
1450 calibration_operation = "multiply"
1451 self.logger.warning(
1452 "Unexpected Inverse Filter is being corrected -- something maybe wrong here "
1453 )
1454 else:
1455 msg = "cannot determine multiply or divide via units -- setting to divide"
1456 self.logger.warning(msg)
1457 calibration_operation = "divide"
1458 return calibration_operation
1460 def get_calibrated_units(self):
1461 """
1462 Follows the FDSN standard which has the filter stages starting with physical units to digital counts.
1464 The channel_response is expected to have a list of filter "stages" of which the first stage
1465 has input units corresponding to the the physical quantity that the instrument measures, and the last is
1466 normally counts.
1468 channel_response can be viewed as the chaining together of all of these filters.
1470 Thus it is normal for channel_response.units_out will be in the same units as the archived raw
1471 time series, and for the units after the response is corrected for will be the units_in of
1473 The units of the channel metadata are compared to the input and output units of the channel_response.
1476 :return: tuple, calibration_operation, either "mulitply" or divide", and a string for calibrated units
1477 :rtype: tuple (of two strings_
1478 """
1480 if self.channel_response.units_out == self.channel_metadata.unit_object.name:
1481 calibrated_units = self.channel_response.units_in
1482 elif (
1483 self.channel_response.units_in == None
1484 and self.channel_response.units_out == None
1485 ):
1486 msg = "No Units are associated with the channel_response"
1487 self.logger.warning(msg)
1488 msg = "cannot determine multiply or divide via units -- setting to divide:/"
1489 self.logger.warning(msg)
1490 calibrated_units = self.channel_metadata.units
1491 else:
1492 logger.critical(
1493 "channel response filter units are likely corrupt or channel_ts has no units"
1494 )
1495 calibrated_units = self.channel_response.units_in
1496 unit_object = get_unit_object(calibrated_units)
1497 calibrated_units = unit_object.name
1498 return calibrated_units
1500 def remove_instrument_response(
1501 self, include_decimation=False, include_delay=False, **kwargs
1502 ):
1503 """
1504 Remove instrument response from the given channel response filter
1506 The order of operations is important (if applied):
1508 1) detrend
1509 2) zero mean
1510 3) zero pad
1511 4) time window
1512 5) frequency window
1513 6) remove response
1514 7) undo time window
1515 8) bandpass
1517 :param include_decimation: Include decimation in response,
1518 defaults to True
1519 :type include_decimation: bool, optional
1520 :param include_delay: include delay in complex response,
1521 defaults to False
1522 :type include_delay: bool, optional
1524 **kwargs**
1526 :param plot: to plot the calibration process [ False | True ]
1527 :type plot: boolean, default True
1528 :param detrend: Remove linar trend of the time series
1529 :type detrend: boolean, default True
1530 :param zero_mean: Remove the mean of the time series
1531 :type zero_mean: boolean, default True
1532 :param zero_pad: pad the time series to the next power of 2 for efficiency
1533 :type zero_pad: boolean, default True
1534 :param t_window: Time domain windown name see `scipy.signal.windows` for options
1535 :type t_window: string, default None
1536 :param t_window_params: Time domain window parameters, parameters can be
1537 found in `scipy.signal.windows`
1538 :type t_window_params: dictionary
1539 :param f_window: Frequency domain windown name see `scipy.signal.windows` for options
1540 :type f_window: string, defualt None
1541 :param f_window_params: Frequency window parameters, parameters can be
1542 found in `scipy.signal.windows`
1543 :type f_window_params: dictionary
1544 :param bandpass: bandpass freequency and order {"low":, "high":, "order":,}
1545 :type bandpass: dictionary
1547 """
1549 def bool_flip(x):
1550 return bool(int(x) - 1)
1552 if hasattr(self.channel_metadata, "filter"):
1553 if self.channel_metadata.filter.applied is []:
1554 self.logger.warning("No filters to apply to calibrate time series data")
1555 return self.copy()
1556 elif self.channel_metadata.filters is []:
1557 self.logger.warning("No filters to apply to calibrate time series data")
1558 return self.copy()
1560 calibrated_ts = self.copy(data=False)
1562 # Make a list of the filters whose response will be removed.
1563 # We make the list here so that we have access to the indices to flip
1564 indices_to_flip = self.channel_response.get_indices_of_filters_to_remove(
1565 include_decimation=include_decimation,
1566 include_delay=include_delay,
1567 )
1568 filters_to_remove = [
1569 self.channel_response.filters_list[i] for i in indices_to_flip
1570 ]
1572 remover = RemoveInstrumentResponse(
1573 self.ts,
1574 self.time_index,
1575 self.sample_interval,
1576 self.channel_response,
1577 **kwargs,
1578 )
1580 calibration_operation = self.get_calibration_operation()
1581 calibrated_ts.ts = remover.remove_instrument_response(
1582 filters_to_remove=filters_to_remove,
1583 operation=calibration_operation,
1584 )
1586 # update "applied" booleans
1587 if hasattr(calibrated_ts.channel_metadata, "filter"):
1588 applied_filters = calibrated_ts.channel_metadata.filter.applied
1589 for idx in indices_to_flip:
1590 applied_filters[idx] = bool_flip(applied_filters[idx])
1591 calibrated_ts.channel_metadata.filter.applied = applied_filters
1592 else:
1593 for idx in indices_to_flip:
1594 calibrated_ts.channel_metadata.filters[idx].applied = bool_flip(
1595 calibrated_ts.channel_metadata.filters[idx].applied
1596 )
1598 # update units
1599 calibrated_units = self.get_calibrated_units()
1600 calibrated_ts.data_array.attrs["units"] = calibrated_units
1601 calibrated_ts.channel_metadata.units = calibrated_units
1602 calibrated_ts._update_xarray_metadata()
1604 return calibrated_ts
1606 def get_slice(self, start, end=None, n_samples=None):
1607 """
1608 Get a slice from the time series given a start and end time.
1610 Looks for >= start & <= end
1612 Uses loc to be exact with milliseconds
1614 :param start: start time of the slice
1615 :type start: string, MTime
1616 :param end: end time of the slice
1617 :type end: string, MTime
1618 :param n_samples: number of sample to get after start time
1619 :type n_samples: integer
1620 :return: slice of the channel requested
1621 :rtype: ChannelTS
1623 """
1625 if n_samples is None and end is None:
1626 msg = "Must input either end_time or n_samples."
1627 self.logger.error(msg)
1628 raise ValueError(msg)
1629 if n_samples is not None and end is not None:
1630 msg = "Must input either end_time or n_samples, not both."
1631 self.logger.error(msg)
1632 raise ValueError(msg)
1633 if not isinstance(start, MTime):
1634 start = MTime(time_stamp=start)
1635 if n_samples is not None:
1636 n_samples = int(n_samples)
1637 end = start + ((n_samples - 1) / self.sample_rate)
1638 if end is not None:
1639 if not isinstance(end, MTime):
1640 end = MTime(time_stamp=end)
1641 chunk = self.data_array.indexes["time"].slice_indexer(
1642 start=np.datetime64(start.iso_no_tz),
1643 end=np.datetime64(end.iso_no_tz),
1644 )
1645 new_ts = self.data_array.isel(indexers={"time": chunk})
1647 new_ch_ts = ChannelTS(
1648 channel_type=self.channel_type,
1649 data=new_ts,
1650 survey_metadata=self.survey_metadata,
1651 channel_response=self.channel_response,
1652 )
1654 return new_ch_ts
1656 # decimate data
1657 def decimate(self, new_sample_rate, inplace=False, max_decimation=8):
1658 """
1659 decimate the data by using scipy.signal.decimate
1661 :param dec_factor: decimation factor
1662 :type dec_factor: int
1664 * refills ts.data with decimated data and replaces sample_rate
1666 """
1668 sr_list = get_decimation_sample_rates(
1669 self.sample_rate, new_sample_rate, max_decimation
1670 )
1672 # need to fill nans with 0 otherwise they wipeout the decimation values
1673 # and all becomes nan.
1674 new_ts = self.data_array.fillna(0)
1675 for step_sr in sr_list:
1676 new_ts = new_ts.sps_filters.decimate(step_sr)
1677 new_ts.attrs["sample_rate"] = new_sample_rate
1678 self.channel_metadata.sample_rate = new_ts.attrs["sample_rate"]
1680 if inplace:
1681 self.ts = new_ts
1682 else:
1683 new_ts.attrs.update(
1684 self.channel_metadata.to_dict()[self.channel_metadata._class_name]
1685 )
1686 # return new_ts
1687 return ChannelTS(
1688 self.channel_metadata.type,
1689 data=new_ts,
1690 metadata=self.channel_metadata,
1691 )
1693 def resample_poly(self, new_sample_rate, pad_type="mean", inplace=False):
1694 """
1695 Use scipy.signal.resample_poly to resample data while using an FIR
1696 filter to remove aliasing.
1698 :param new_sample_rate: DESCRIPTION
1699 :type new_sample_rate: TYPE
1700 :param pad_type: DESCRIPTION, defaults to "mean"
1701 :type pad_type: TYPE, optional
1702 :return: DESCRIPTION
1703 :rtype: TYPE
1705 """
1707 # need to fill nans with 0 otherwise they wipeout the decimation values
1708 # and all becomes nan.
1709 new_ts = self.data_array.fillna(0)
1710 new_ts = new_ts.sps_filters.resample_poly(new_sample_rate, pad_type=pad_type)
1712 new_ts.attrs["sample_rate"] = new_sample_rate
1713 self.channel_metadata.sample_rate = new_ts.attrs["sample_rate"]
1715 if inplace:
1716 self.ts = new_ts
1717 else:
1718 new_ts.attrs.update(
1719 self.channel_metadata.to_dict()[self.channel_metadata._class_name]
1720 )
1721 # return new_ts
1722 return ChannelTS(
1723 self.channel_metadata.type,
1724 data=new_ts,
1725 metadata=self.channel_metadata,
1726 )
1728 def merge(
1729 self,
1730 other,
1731 gap_method="slinear",
1732 new_sample_rate=None,
1733 resample_method="poly",
1734 ):
1735 """
1736 merg two channels or list of channels together in the following steps
1738 1. xr.combine_by_coords([original, other])
1739 2. compute monotonic time index
1740 3. reindex(new_time_index, method=gap_method)
1742 If you want a different method or more control use merge
1744 :param other: Another channel
1745 :type other: :class:`mth5.timeseries.ChannelTS`
1746 :raises TypeError: If input is not a ChannelTS
1747 :raises ValueError: if the components are different
1748 :return: Combined channel with monotonic time index and same metadata
1749 :rtype: :class:`mth5.timeseries.ChannelTS`
1751 """
1752 if new_sample_rate is not None:
1753 merge_sample_rate = new_sample_rate
1754 if resample_method == "decimate":
1755 combine_list = [self.decimate(new_sample_rate).data_array]
1756 elif resample_method == "poly":
1757 combine_list = [self.resample_poly(new_sample_rate).data_array]
1758 else:
1759 merge_sample_rate = self.sample_rate
1760 combine_list = [self.data_array]
1761 if isinstance(other, (list, tuple)):
1762 for ch in other:
1763 if not isinstance(ch, ChannelTS):
1764 raise TypeError(f"Cannot combine {type(ch)} with ChannelTS.")
1765 if self.component != ch.component:
1766 raise ValueError(
1767 "Cannot combine channels with different components. "
1768 f"{self.component} != {ch.component}"
1769 )
1770 if new_sample_rate is not None:
1771 if resample_method == "decimate":
1772 ch = ch.decimate(new_sample_rate)
1773 elif resample_method == "poly":
1774 ch = ch.resample_poly(new_sample_rate)
1775 combine_list.append(ch.data_array)
1776 else:
1777 if not isinstance(other, ChannelTS):
1778 raise TypeError(f"Cannot combine {type(other)} with ChannelTS.")
1779 if self.component != other.component:
1780 raise ValueError(
1781 "Cannot combine channels with different components. "
1782 f"{self.component} != {other.component}"
1783 )
1784 if new_sample_rate is not None:
1785 if resample_method == "decimate":
1786 other = other.decimate(new_sample_rate)
1787 elif resample_method == "poly":
1788 other = other.resample_poly(new_sample_rate)
1789 combine_list.append(other.data_array)
1790 # combine into a data set use override to keep attrs from original
1792 combined_ds = xr.combine_by_coords(combine_list, combine_attrs="override")
1794 # compute duration between max and min times robustly (handles
1795 # numpy.timedelta64 and datetime.timedelta across Python versions)
1796 time_max = combined_ds.time.max().values
1797 time_min = combined_ds.time.min().values
1798 delta = time_max - time_min
1800 try:
1801 # numpy.timedelta64 -> get seconds via division
1802 seconds = delta / np.timedelta64(1, "s")
1803 except Exception:
1804 # fallback for datetime.timedelta or other types
1805 import datetime
1807 if isinstance(delta, datetime.timedelta):
1808 seconds = delta.total_seconds()
1809 else:
1810 # as a last resort, use pandas to_timedelta which handles
1811 # various timedelta-like representations
1812 seconds = pd.to_timedelta(delta).total_seconds()
1814 n_samples = int(merge_sample_rate * float(seconds)) + 1
1816 new_dt_index = make_dt_coordinates(
1817 combined_ds.time.min().values, merge_sample_rate, n_samples
1818 )
1820 channel_metadata = self.channel_metadata.copy()
1821 channel_metadata.sample_rate = merge_sample_rate
1822 run_metadata = self.run_metadata.copy()
1823 run_metadata.sample_rate = merge_sample_rate
1825 new_channel = ChannelTS(
1826 channel_type=self.channel_metadata.type,
1827 channel_metadata=channel_metadata,
1828 run_metadata=self.run_metadata,
1829 station_metadata=self.station_metadata,
1830 survey_metadata=self.survey_metadata,
1831 channel_response=self.channel_response,
1832 )
1834 new_channel.data_array = combined_ds.interp(
1835 time=new_dt_index, method=gap_method
1836 ).to_array()
1838 new_channel.channel_metadata.time_period.start = new_channel.start
1839 new_channel.channel_metadata.time_period.end = new_channel.end
1841 new_channel.run_metadata.update_time_period()
1842 new_channel.station_metadata.update_time_period()
1843 new_channel.survey_metadata.update_time_period()
1845 new_channel._update_xarray_metadata()
1847 return new_channel
1849 def to_xarray(self):
1850 """
1851 Returns a :class:`xarray.DataArray` object of the channel timeseries
1852 this way metadata from the metadata class is updated upon return.
1854 :return: Returns a :class:`xarray.DataArray` object of the channel timeseries
1855 this way metadata from the metadata class is updated upon return.
1856 :rtype: :class:`xarray.DataArray`
1859 >>> import numpy as np
1860 >>> from mth5.timeseries import ChannelTS
1861 >>> ex = ChannelTS("electric")
1862 >>> ex.start = "2020-01-01T12:00:00"
1863 >>> ex.sample_rate = 16
1864 >>> ex.ts = np.random.rand(4096)
1867 """
1868 self._update_xarray_metadata()
1869 return self.data_array
1871 def to_obspy_trace(self, network_code=None, encoding=None):
1872 """
1873 Convert the time series to an :class:`obspy.core.trace.Trace` object. This
1874 will be helpful for converting between data pulled from IRIS and data going
1875 into IRIS.
1877 :param network_code: two letter code provided by FDSN DMC
1878 :type network_code: string
1879 :return: DESCRIPTION
1880 :rtype: TYPE
1882 """
1883 encoding_dict = {
1884 "INT16": np.int16,
1885 "INT32": np.int32,
1886 "INT64": np.int32,
1887 "FLOAT32": np.float32,
1888 "FLOAT64": np.float64,
1889 }
1890 if self.ts.dtype.type in [np.int64]:
1891 obspy_trace = Trace(self.ts.astype(np.int32))
1893 if encoding:
1894 try:
1895 obspy_trace = Trace(self.ts.astype(encoding_dict[encoding]))
1896 except KeyError:
1897 raise KeyError(
1898 f"{encoding} is not understood. Acceptable values are {list(encoding_dict.keys())}"
1899 )
1900 else:
1901 obspy_trace = Trace(self.ts)
1903 # add metadata
1904 obspy_trace.stats.channel = fdsn_tools.make_channel_code(self.channel_metadata)
1905 obspy_trace.stats.starttime = self.start.isoformat()
1906 obspy_trace.stats.sampling_rate = self.sample_rate
1907 if self.station_metadata.fdsn.id is None:
1908 self.station_metadata.fdsn.id = self.station_metadata.id
1909 obspy_trace.stats.station = self.station_metadata.fdsn.id.upper()
1910 obspy_trace.stats.network = network_code
1912 return obspy_trace
1914 def from_obspy_trace(self, obspy_trace):
1915 """
1916 Fill data from an :class:`obspy.core.Trace`
1918 :param obspy.core.trace obspy_trace: Obspy trace object
1920 """
1922 if not isinstance(obspy_trace, Trace):
1923 msg = f"Input must be obspy.core.Trace, not {type(obspy_trace)}"
1924 self.logger.error(msg)
1925 raise TypeError(msg)
1926 if obspy_trace.stats.channel[1].lower() in ["e", "q"]:
1927 self.channel_type = "electric"
1928 measurement = "electric"
1929 elif obspy_trace.stats.channel[1].lower() in ["h", "b", "f"]:
1930 self.channel_type = "magnetic"
1931 measurement = "magnetic"
1932 else:
1933 try:
1934 measurement = fdsn_tools.measurement_code_dict_reverse[
1935 obspy_trace.stats.channel[1]
1936 ]
1937 except KeyError:
1938 measurement = "auxiliary"
1939 self.channel_type = "auxiliary"
1940 mt_code = fdsn_tools.make_mt_channel(
1941 fdsn_tools.read_channel_code(obspy_trace.stats.channel)
1942 )
1944 self.channel_metadata.component = mt_code
1945 self.channel_metadata.type = measurement
1946 self.sample_rate = obspy_trace.stats.sampling_rate
1947 self.start = obspy_trace.stats.starttime.isoformat()
1948 self.station_metadata.fdsn.id = obspy_trace.stats.station
1949 # Handle None network values
1950 if (
1951 obspy_trace.stats.network is not None
1952 and obspy_trace.stats.network != "None"
1953 ):
1954 self.station_metadata.fdsn.network = obspy_trace.stats.network
1955 self.station_metadata.id = obspy_trace.stats.station
1956 self.channel_metadata.units = "counts"
1957 self.ts = obspy_trace.data
1958 self.run_metadata.id = f"sr{int(self.sample_rate)}_001"
1960 def plot(self):
1961 """
1962 Simple plot of the data
1964 :return: figure object
1965 :rtype: matplotlib.figure
1967 """
1969 return self.data_array.plot()
1971 def welch_spectra(self, window_length=2**12, **kwargs):
1972 """
1973 get welch spectra
1975 :param window_length: DESCRIPTION
1976 :type window_length: TYPE
1977 :param **kwargs: DESCRIPTION
1978 :type **kwargs: TYPE
1979 :return: DESCRIPTION
1980 :rtype: TYPE
1982 """
1984 plot_frequency, power = signal.welch(
1985 self.ts, fs=self.sample_rate, nperseg=window_length, **kwargs
1986 )
1988 return plot_frequency, power
1990 def plot_spectra(self, spectra_type="welch", window_length=2**12, **kwargs):
1991 """
1993 :param spectra_type: spectra type, defaults to "welch"
1994 :type spectra_type: string, optional
1995 :param window_length: window length of the welch method should be a
1996 power of 2, defaults to 2 ** 12
1997 :type window_length: int, optional
1998 :param **kwargs: DESCRIPTION
1999 :type **kwargs: TYPE
2001 """
2002 from matplotlib import pyplot as plt
2004 if spectra_type == "welch":
2005 plot_frequency, power = self.welch_spectra(
2006 window_length=window_length, **kwargs
2007 )
2008 fig = plt.figure()
2009 ax = fig.add_subplot(1, 1, 1)
2010 ax.loglog(1.0 / plot_frequency, power, lw=1.5)
2011 ax.set_xlabel("Period (s)", fontdict={"size": 10, "weight": "bold"})
2012 ax.set_ylabel("Power (dB)", fontdict={"size": 10, "weight": "bold"})
2013 ax.axis("tight")
2014 ax.grid(which="both")
2015 ax2 = ax.twiny()
2016 ax2.loglog(plot_frequency, power, lw=0)
2017 ax2.set_xlabel("Frequency (Hz)", fontdict={"size": 10, "weight": "bold"})
2018 ax2.set_xlim([1 / cc for cc in ax.get_xlim()])
2019 plt.show()