From 082fa73720f6ff95efb4d4baa7ef8866063b9be0 Mon Sep 17 00:00:00 2001 From: "Karl N. Kappler" Date: Fri, 8 Nov 2024 16:50:53 -0800 Subject: [PATCH 01/11] housekeeping --- .../parkfield/make_parkfield_mth5.py | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/aurora/test_utils/parkfield/make_parkfield_mth5.py b/aurora/test_utils/parkfield/make_parkfield_mth5.py index 843e1d98..4d4e95cd 100644 --- a/aurora/test_utils/parkfield/make_parkfield_mth5.py +++ b/aurora/test_utils/parkfield/make_parkfield_mth5.py @@ -6,21 +6,17 @@ import pathlib from aurora.test_utils.dataset_definitions import TEST_DATA_SET_CONFIGS -from mth5.utils.helpers import read_back_data -from mth5.helpers import close_open_files -from aurora.sandbox.io_helpers.fdsn_dataset import FDSNDataset from aurora.sandbox.io_helpers.make_mth5_helpers import create_from_server_multistation from aurora.test_utils.parkfield.path_helpers import PARKFIELD_PATHS from loguru import logger -from typing import Union - +from mth5.utils.helpers import read_back_data +from mth5.helpers import close_open_files +from typing import Optional, Union DATA_SOURCES = ["NCEDC", "https://service.ncedc.org/"] DATASET_ID = "pkd_sao_test_00" FDSN_DATASET = TEST_DATA_SET_CONFIGS[DATASET_ID] -# - def select_data_source() -> None: """ @@ -44,16 +40,17 @@ def select_data_source() -> None: except: logger.warning(f"Data source {data_source} not initializing") if not ok: - logger.error("No data sources for Parkfield / Hollister initializing") - logger.error("NCEDC probably down") - raise ValueError + msg = "No data sources for Parkfield / Hollister initializing\n" + msg += "NCEDC probably down" + logger.error(msg) + raise ValueError(msg) else: return data_source def make_pkdsao_mth5( fdsn_dataset: FDSN_DATASET, - target_folder: Union[str, pathlib.Path, None] = PARKFIELD_PATHS["data"], + target_folder: Optional[Union[str, pathlib.Path]] = PARKFIELD_PATHS["data"], ) -> pathlib.Path: """ Makes MTH5 file with data from Parkfield and Hollister stations to use for testing. @@ -83,12 +80,13 @@ def ensure_h5_exists( Parameters ---------- h5_path: Union[pathlib.Path, None] - + The target path to build to mth5. Returns ------- h5_path: pathlib.Path The path to the PKD SAO mth5 file to be used for testing. + """ h5_path = target_folder.joinpath(FDSN_DATASET.h5_filebase) if h5_path.exists(): From f0c41d03a0b59ef3a9b02d191fdcf7d627756a8e Mon Sep 17 00:00:00 2001 From: "Karl N. Kappler" Date: Sat, 30 Nov 2024 16:14:01 -0800 Subject: [PATCH 02/11] remove uneeded ipynb form tests - also, pytest was complaining about it locally, and it doesn't seem worth looking into --- .github/workflows/tests.yml | 1 - tests/test_run_on_commit.ipynb | 35 ---------------------------------- 2 files changed, 36 deletions(-) delete mode 100644 tests/test_run_on_commit.ipynb diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index b12d023e..59ad13e8 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -63,7 +63,6 @@ jobs: jupyter nbconvert --to notebook --execute docs/tutorials/processing_configuration.ipynb jupyter nbconvert --to notebook --execute docs/tutorials/process_cas04_multiple_station.ipynb jupyter nbconvert --to notebook --execute docs/tutorials/synthetic_data_processing.ipynb - jupyter nbconvert --to notebook --execute tests/test_run_on_commit.ipynb # Replace "notebook.ipynb" with your notebook's filename # - name: Commit changes (if any) diff --git a/tests/test_run_on_commit.ipynb b/tests/test_run_on_commit.ipynb deleted file mode 100644 index 6ac8e9bc..00000000 --- a/tests/test_run_on_commit.ipynb +++ /dev/null @@ -1,35 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "99fe1b05-7951-4666-8c24-7bb28426611d", - "metadata": {}, - "outputs": [], - "source": [ - "assert(True)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "aurora-test", - "language": "python", - "name": "aurora-test" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.10" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From 89e000dbbd2a63cf487f8be60d955854ab25d0a5 Mon Sep 17 00:00:00 2001 From: "Karl N. Kappler" Date: Tue, 17 Dec 2024 10:10:48 -0800 Subject: [PATCH 03/11] Merge minor changes from fix_mtpy_issue_67 - these don't hurt and make the flow slightly more robust --- aurora/pipelines/process_mth5.py | 64 ++++++++----------- .../transfer_function/weights/edf_weights.py | 12 ++-- 2 files changed, 32 insertions(+), 44 deletions(-) diff --git a/aurora/pipelines/process_mth5.py b/aurora/pipelines/process_mth5.py index 94250a3a..4ad1e15e 100644 --- a/aurora/pipelines/process_mth5.py +++ b/aurora/pipelines/process_mth5.py @@ -63,9 +63,7 @@ # ============================================================================= -def make_stft_objects( - processing_config, i_dec_level, run_obj, run_xrds, units="MT" -): +def make_stft_objects(processing_config, i_dec_level, run_obj, run_xrds, units="MT"): """ Operates on a "per-run" basis. Applies STFT to all time series in the input run. @@ -103,9 +101,7 @@ def make_stft_objects( ].channel_scale_factors elif run_obj.station_metadata.id == processing_config.stations.remote[0].id: scale_factors = ( - processing_config.stations.remote[0] - .run_dict[run_id] - .channel_scale_factors + processing_config.stations.remote[0].run_dict[run_id].channel_scale_factors ) stft_obj = calibrate_stft_obj( @@ -152,9 +148,7 @@ def process_tf_decimation_level( The transfer function values packed into an object """ frequency_bands = config.decimations[i_dec_level].frequency_bands_obj() - transfer_function_obj = TTFZ( - i_dec_level, frequency_bands, processing_config=config - ) + transfer_function_obj = TTFZ(i_dec_level, frequency_bands, processing_config=config) dec_level_config = config.decimations[i_dec_level] # segment_weights = coherence_weights(dec_level_config, local_stft_obj, remote_stft_obj) transfer_function_obj = process_transfer_functions( @@ -183,9 +177,7 @@ def triage_issue_289(local_stfts: list, remote_stfts: list): for i_chunk in range(n_chunks): ok = local_stfts[i_chunk].time.shape == remote_stfts[i_chunk].time.shape if not ok: - logger.warning( - "Mismatch in FC array lengths detected -- Issue #289" - ) + logger.warning("Mismatch in FC array lengths detected -- Issue #289") glb = max( local_stfts[i_chunk].time.min(), remote_stfts[i_chunk].time.min(), @@ -196,18 +188,13 @@ def triage_issue_289(local_stfts: list, remote_stfts: list): ) cond1 = local_stfts[i_chunk].time >= glb cond2 = local_stfts[i_chunk].time <= lub - local_stfts[i_chunk] = local_stfts[i_chunk].where( - cond1 & cond2, drop=True - ) + local_stfts[i_chunk] = local_stfts[i_chunk].where(cond1 & cond2, drop=True) cond1 = remote_stfts[i_chunk].time >= glb cond2 = remote_stfts[i_chunk].time <= lub remote_stfts[i_chunk] = remote_stfts[i_chunk].where( cond1 & cond2, drop=True ) - assert ( - local_stfts[i_chunk].time.shape - == remote_stfts[i_chunk].time.shape - ) + assert local_stfts[i_chunk].time.shape == remote_stfts[i_chunk].time.shape return local_stfts, remote_stfts @@ -306,9 +293,7 @@ def load_stft_obj_from_mth5( An STFT from mth5. """ station_obj = station_obj_from_row(row) - fc_group = station_obj.fourier_coefficients_group.get_fc_group( - run_obj.metadata.id - ) + fc_group = station_obj.fourier_coefficients_group.get_fc_group(run_obj.metadata.id) fc_decimation_level = fc_group.get_decimation_level(f"{i_dec_level}") stft_obj = fc_decimation_level.to_xarray(channels=channels) @@ -369,10 +354,7 @@ def save_fourier_coefficients(dec_level_config, row, run_obj, stft_obj) -> None: raise NotImplementedError(msg) # Get FC group (create if needed) - if ( - run_obj.metadata.id - in station_obj.fourier_coefficients_group.groups_list - ): + if run_obj.metadata.id in station_obj.fourier_coefficients_group.groups_list: fc_group = station_obj.fourier_coefficients_group.get_fc_group( run_obj.metadata.id ) @@ -393,9 +375,7 @@ def save_fourier_coefficients(dec_level_config, row, run_obj, stft_obj) -> None: dec_level_name, decimation_level_metadata=decimation_level_metadata, ) - fc_decimation_level.from_xarray( - stft_obj, decimation_level_metadata.sample_rate - ) + fc_decimation_level.from_xarray(stft_obj, decimation_level_metadata.sample_rate) fc_decimation_level.update_metadata() fc_group.update_metadata() else: @@ -535,9 +515,7 @@ def process_mth5_legacy( local_merged_stft_obj, remote_merged_stft_obj, ) - ttfz_obj.apparent_resistivity( - tfk.config.channel_nomenclature, units=units - ) + ttfz_obj.apparent_resistivity(tfk.config.channel_nomenclature, units=units) tf_dict[i_dec_level] = ttfz_obj if show_plot: @@ -549,10 +527,20 @@ def process_mth5_legacy( tf_dict=tf_dict, processing_config=tfk.config ) - tf_cls = tfk.export_tf_collection(tf_collection) - - if z_file_path: - tf_cls.write(z_file_path) + try: + tf_cls = tfk.export_tf_collection(tf_collection) + if z_file_path: + tf_cls.write(z_file_path) + except Exception as e: + msg = "TF collection could not export to mt_metadata TransferFunction\n" + msg += f"Failed with exception {e}\n" + msg += "Perhaps an unconventional mixture of input/output channels was used\n" + msg += f"Input channels were {tfk.config.decimations[0].input_channels}\n" + msg += f"Output channels were {tfk.config.decimations[0].output_channels}\n" + msg += "No z-file will be written in this case\n" + msg += "Will return a legacy TransferFunctionCollection object, not mt_metadata object." + logger.error(msg) + return_collection = True tfk.dataset.close_mth5s() if return_collection: @@ -602,9 +590,7 @@ def process_mth5( The transfer function object """ if processing_type not in SUPPORTED_PROCESSINGS: - raise NotImplementedError( - f"Processing type {processing_type} not supported" - ) + raise NotImplementedError(f"Processing type {processing_type} not supported") if processing_type == "legacy": try: diff --git a/aurora/transfer_function/weights/edf_weights.py b/aurora/transfer_function/weights/edf_weights.py index 590cfa12..4c8e104c 100644 --- a/aurora/transfer_function/weights/edf_weights.py +++ b/aurora/transfer_function/weights/edf_weights.py @@ -198,18 +198,20 @@ def effective_degrees_of_freedom_weights( Weights for reducing leverage points. """ + # Initialize the weights + n_observations_initial = len(X.observation) + weights = np.ones(n_observations_initial) + # validate num channels num_channels = len(X.data_vars) if num_channels != 2: - logger.error("edfwts only works for 2 input channels") - raise Exception + logger.error(f"edfwts only works for 2 input channels, not {num_channels}") + return weights + X = X.to_array(dim="channel") if R is not None: R = R.to_array(dim="channel") - n_observations_initial = len(X.observation) - weights = np.ones(n_observations_initial) - # reduce the data to only valid (non-nan) observations if R is not None: keep_x_indices = ~np.isnan(X.data).any(axis=0) From e4df7e68fa15659fab9416a05dda535c7e27427c Mon Sep 17 00:00:00 2001 From: kkappler Date: Thu, 19 Dec 2024 21:40:04 -0800 Subject: [PATCH 04/11] housekeeping / doc strings --- aurora/pipelines/fourier_coefficients.py | 43 +++++++++++++----------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/aurora/pipelines/fourier_coefficients.py b/aurora/pipelines/fourier_coefficients.py index 3b1c93c8..97a0a420 100644 --- a/aurora/pipelines/fourier_coefficients.py +++ b/aurora/pipelines/fourier_coefficients.py @@ -61,7 +61,6 @@ # Imports # ============================================================================= -import mt_metadata.timeseries.time_period import mth5.mth5 import pathlib @@ -71,6 +70,7 @@ from loguru import logger from mth5.mth5 import MTH5 from mth5.utils.helpers import path_or_mth5_object +from mt_metadata.timeseries.time_period import TimePeriod from mt_metadata.transfer_functions.processing.fourier_coefficients import ( Decimation as FCDecimation, ) @@ -82,10 +82,10 @@ def fc_decimations_creator( initial_sample_rate: float, - decimation_factors: Optional[Union[list, None]] = None, + decimation_factors: Optional[list] = None, max_levels: Optional[int] = 6, - time_period: mt_metadata.timeseries.TimePeriod = None, -) -> list: + time_period: Optional[TimePeriod] = None, +) -> list[FCDecimation]: """ Creates mt_metadata FCDecimation objects that parameterize Fourier coefficient decimation levels. @@ -97,11 +97,12 @@ def fc_decimations_creator( ---------- initial_sample_rate: float Sample rate of the "level0" data -- usually the sample rate during field acquisition. - decimation_factors: list (or other iterable) + decimation_factors: Optional[list] The decimation factors that will be applied at each FC decimation level - max_levels: int + max_levels: Optional[int] The maximum number of decimation levels to allow - time_period: + time_period: Optional[TimePeriod] + Provides the start and end times Returns ------- @@ -137,7 +138,7 @@ def fc_decimations_creator( fc_dec.sample_rate_decimation = current_sample_rate if time_period: - if isinstance(time_period, mt_metadata.timeseries.time_period.TimePeriod): + if isinstance(time_period, TimePeriod): fc_dec.time_period = time_period else: msg = ( @@ -152,15 +153,14 @@ def fc_decimations_creator( @path_or_mth5_object -def add_fcs_to_mth5( - m: MTH5, fc_decimations: Optional[Union[list, None]] = None -) -> None: +def add_fcs_to_mth5(m: MTH5, fc_decimations: Optional[Union[str, list]] = None) -> None: """ Add Fourier Coefficient Levels ot an existing MTH5. **Notes:** - - This module computes the FCs differently than the legacy aurora pipeline. It uses scipy.signal.spectrogram. There is a test in Aurora to confirm that there are equivalent if we are not using fancy pre-whitening. + - This module computes the FCs differently than the legacy aurora pipeline. It uses scipy.signal.spectrogram. + There is a test in Aurora to confirm that there are equivalent if we are not using fancy pre-whitening. - Nomenclature: "usssr_grouper" is the output of a group-by on unique {survey, station, sample_rate} tuples. @@ -168,10 +168,11 @@ def add_fcs_to_mth5( ---------- m: MTH5 object The mth5 file, open in append mode. - fc_decimations: Union[str, None, List] + fc_decimations: Optional[Union[str, list]] This specifies the scheme to use for decimating the time series when building the FC layer. - None: Just use default (something like four decimation levels, decimated by 4 each time say. - String: Controlled Vocabulary, values are a work in progress, that will allow custom definition of the fc_decimations for some common cases. For example, say you have stored already decimated time + None: Just use default (something like four decimation levels, decimated by 4 each time say.) + String: Controlled Vocabulary, values are a work in progress, that will allow custom definition of + the fc_decimations for some common cases. For example, say you have stored already decimated time series, then you want simply the zeroth decimation for each run, because the decimated time series live under another run container, and that will get its own FCs. This is experimental. List: (**UNTESTED**) -- This means that the user thought about the decimations that they want to create and is @@ -287,17 +288,19 @@ def get_degenerate_fc_decimation(sample_rate: float) -> list: @path_or_mth5_object -def read_back_fcs(m: Union[MTH5, pathlib.Path, str], mode="r"): +def read_back_fcs(m: Union[MTH5, pathlib.Path, str], mode: str = "r") -> None: """ - This is mostly a helper function for tests. It was used as a sanity check while debugging the FC files, and + + This is a helper function for tests. It was used as a sanity check while debugging the FC files, and also is a good example for how to access the data at each level for each channel. The Time axis of the FC array will change from level to level, but the frequency axis will stay the same shape (for now -- storing all fcs by default) - Args: - m: pathlib.Path, str or an MTH5 object - The path to an h5 file that we will scan the fcs from + Parameters + ---------- + m: Union[MTH5, pathlib.Path, str] + Either a path to an mth5, or an MTH5 object that the FCs will be read back from. """ From 49b1f33d03dc073bcefdb036aa9b3d6eec5f25fa Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 08:17:02 -0800 Subject: [PATCH 05/11] fix typehint --- aurora/pipelines/fourier_coefficients.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aurora/pipelines/fourier_coefficients.py b/aurora/pipelines/fourier_coefficients.py index 97a0a420..d67acb7a 100644 --- a/aurora/pipelines/fourier_coefficients.py +++ b/aurora/pipelines/fourier_coefficients.py @@ -85,7 +85,7 @@ def fc_decimations_creator( decimation_factors: Optional[list] = None, max_levels: Optional[int] = 6, time_period: Optional[TimePeriod] = None, -) -> list[FCDecimation]: +) -> List[FCDecimation]: """ Creates mt_metadata FCDecimation objects that parameterize Fourier coefficient decimation levels. From a7406dd5ab86aac70f7796a9aaca018ab0765020 Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 16:39:01 -0800 Subject: [PATCH 06/11] housekeeping & typehinting --- aurora/pipelines/time_series_helpers.py | 66 ++++++++++++++------ tests/synthetic/test_fourier_coefficients.py | 9 +++ 2 files changed, 57 insertions(+), 18 deletions(-) diff --git a/aurora/pipelines/time_series_helpers.py b/aurora/pipelines/time_series_helpers.py index 5102d7c2..877c10ac 100644 --- a/aurora/pipelines/time_series_helpers.py +++ b/aurora/pipelines/time_series_helpers.py @@ -11,11 +11,23 @@ from aurora.time_series.windowed_time_series import WindowedTimeSeries from aurora.time_series.windowing_scheme import window_scheme_from_decimation - - -def validate_sample_rate(run_ts, expected_sample_rate, tol=1e-4): +from mt_metadata.transfer_functions.processing.aurora.decimation_level import ( + DecimationLevel as AuroraDecimationLevel, +) +from mt_metadata.transfer_functions.processing.fourier_coefficients import ( + Decimation as FCDecimation, +) +from mth5.groups import RunGroup +from mth5.timeseries import RunTS +from typing import Optional, Union + + +def validate_sample_rate( + run_ts: RunTS, expected_sample_rate: float, tol: Optional[float] = 1e-4 +) -> None: """ Check that the sample rate of a run_ts is the expected value, and warn if not. + :raise ValueError if outside tolderance Parameters ---------- @@ -24,6 +36,8 @@ def validate_sample_rate(run_ts, expected_sample_rate, tol=1e-4): expected_sample_rate: float The sample rate the time series is expected to have. Normally taken from the processing config + tol: float + This is the tolerance for the difference in sample rates allowed before an exception is raised. """ if run_ts.sample_rate != expected_sample_rate: @@ -34,22 +48,25 @@ def validate_sample_rate(run_ts, expected_sample_rate, tol=1e-4): logger.warning(msg) delta = run_ts.sample_rate - expected_sample_rate if np.abs(delta) > tol: - msg = f"Delta sample rate {delta} > {tol} tolerance" - msg += "TOL should be a percentage" - raise Exception(msg) + msg = f"Difference between expected sample rate and run_ts sample rate {delta} > {tol} tolerance" + raise ValueError(msg) -def apply_prewhitening(decimation_obj, run_xrds_input): +def apply_prewhitening( + decimation_obj: Union[AuroraDecimationLevel, FCDecimation], + run_xrds_input: xr.Dataset, +) -> xr.Dataset: """ Applies pre-whitening to time series to avoid spectral leakage when FFT is applied. - If "first difference", may want to consider clipping first and last sample from - the differentiated time series. + TODO: If "first difference", consider clipping first and last sample from the + differentiated time series. Parameters ---------- - decimation_obj : mt_metadata.transfer_functions.processing.aurora.DecimationLevel + decimation_obj : Union[AuroraDecimationLevel, FCDecimation] Information about how the decimation level is to be processed + run_xrds_input : xarray.core.dataset.Dataset Time series to be pre-whitened @@ -64,7 +81,6 @@ def apply_prewhitening(decimation_obj, run_xrds_input): if decimation_obj.prewhitening_type == "first difference": run_xrds = run_xrds_input.differentiate("time") - else: msg = f"{decimation_obj.prewhitening_type} pre-whitening not implemented" logger.exception(msg) @@ -72,7 +88,10 @@ def apply_prewhitening(decimation_obj, run_xrds_input): return run_xrds -def apply_recoloring(decimation_obj, stft_obj): +def apply_recoloring( + decimation_obj: Union[AuroraDecimationLevel, FCDecimation], + stft_obj: xr.Dataset, +) -> xr.Dataset: """ Inverts the pre-whitening operation in frequency domain. @@ -118,7 +137,7 @@ def apply_recoloring(decimation_obj, stft_obj): return stft_obj -def run_ts_to_stft_scipy(decimation_obj, run_xrds_orig): +def run_ts_to_stft_scipy(decimation_obj: FCDecimation, run_xrds_orig: xr.Dataset): """ Converts a runts object into a time series of Fourier coefficients. This method uses scipy.signal.spectrogram. @@ -173,7 +192,10 @@ def run_ts_to_stft_scipy(decimation_obj, run_xrds_orig): return stft_obj -def truncate_to_clock_zero(decimation_obj, run_xrds): +def truncate_to_clock_zero( + decimation_obj: Union[AuroraDecimationLevel, FCDecimation], + run_xrds: RunGroup, +): """ Compute the time interval between the first data sample and the clock zero Identify the first sample in the xarray time series that corresponds to a @@ -242,7 +264,9 @@ def nan_to_mean(xrds: xr.Dataset) -> xr.Dataset: return xrds -def run_ts_to_stft(decimation_obj, run_xrds_orig): +def run_ts_to_stft( + decimation_obj: AuroraDecimationLevel, run_xrds_orig: xr.Dataset +) -> xr.Dataset: """ Converts a runts object into a time series of Fourier coefficients. Similar to run_ts_to_stft_scipy, but in this implementation operations on individual @@ -291,7 +315,9 @@ def run_ts_to_stft(decimation_obj, run_xrds_orig): return stft_obj -def calibrate_stft_obj(stft_obj, run_obj, units="MT", channel_scale_factors=None): +def calibrate_stft_obj( + stft_obj: xr.Dataset, run_obj, units="MT", channel_scale_factors=None +) -> xr.Dataset: """ Calibrates frequency domain data into MT units. @@ -326,7 +352,7 @@ def calibrate_stft_obj(stft_obj, run_obj, units="MT", channel_scale_factors=None logger.warning(msg) if channel_id == "hy": msg = "Channel hy has no filters, try using filters from hx" - logger.warning("Channel HY has no filters, try using filters from HX") + logger.warning(msg) channel_response = run_obj.get_channel("hx").channel_response indices_to_flip = channel_response.get_indices_of_filters_to_remove( @@ -338,18 +364,22 @@ def calibrate_stft_obj(stft_obj, run_obj, units="MT", channel_scale_factors=None filters_to_remove = [channel_response.filters_list[i] for i in indices_to_flip] if not filters_to_remove: logger.warning("No filters to remove") + calibration_response = channel_response.complex_response( stft_obj.frequency.data, filters_list=filters_to_remove ) + if channel_scale_factors: try: channel_scale_factor = channel_scale_factors[channel_id] except KeyError: channel_scale_factor = 1.0 calibration_response /= channel_scale_factor + if units == "SI": logger.warning("Warning: SI Units are not robustly supported issue #36") - # TODO: This often raises a runtime warning due to DC term in calibration response=0 + + # TODO: FIXME Sometimes raises a runtime warning due to DC term in calibration response = 0 stft_obj[channel_id].data /= calibration_response return stft_obj diff --git a/tests/synthetic/test_fourier_coefficients.py b/tests/synthetic/test_fourier_coefficients.py index 0a454b33..c8b8199a 100644 --- a/tests/synthetic/test_fourier_coefficients.py +++ b/tests/synthetic/test_fourier_coefficients.py @@ -136,6 +136,15 @@ def test_fc_decimations_creator(self): fc_decimations_creator(1.0, time_period=time_period) return cfgs + def test_spectrogram(self): + """ + Place holder method. TODO: Move this into MTH5 + + Development Notes: + Currently mth5 does not have any STFT methods. Once that + :return: + """ + def test_create_then_use_stored_fcs_for_processing(self): """""" from test_processing import process_synthetic_2 From c81729c27f089c8adf6fe4708c7fe7cbb0e05900 Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 16:51:10 -0800 Subject: [PATCH 07/11] doc and typehints --- aurora/pipelines/time_series_helpers.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/aurora/pipelines/time_series_helpers.py b/aurora/pipelines/time_series_helpers.py index 877c10ac..a2c0efb9 100644 --- a/aurora/pipelines/time_series_helpers.py +++ b/aurora/pipelines/time_series_helpers.py @@ -19,7 +19,7 @@ ) from mth5.groups import RunGroup from mth5.timeseries import RunTS -from typing import Optional, Union +from typing import Literal, Optional, Union def validate_sample_rate( @@ -77,6 +77,8 @@ def apply_prewhitening( """ if not decimation_obj.prewhitening_type: + msg = "No prewhitening specified - skipping this step" + logger.info(msg) return run_xrds_input if decimation_obj.prewhitening_type == "first difference": @@ -316,7 +318,10 @@ def run_ts_to_stft( def calibrate_stft_obj( - stft_obj: xr.Dataset, run_obj, units="MT", channel_scale_factors=None + stft_obj: xr.Dataset, + run_obj: RunGroup, + units: Literal["MT", "SI"] = "MT", + channel_scale_factors: Optional[dict] = None, ) -> xr.Dataset: """ Calibrates frequency domain data into MT units. @@ -385,7 +390,7 @@ def calibrate_stft_obj( def prototype_decimate( - config: mt_metadata.transfer_functions.processing.aurora.decimation.Decimation, + config: AuroraDecimationLevel, run_xrds: xr.Dataset, ) -> xr.Dataset: """ From b7cb76ccd881b30df420ec2b9519bf4e09c0eeb3 Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 17:12:16 -0800 Subject: [PATCH 08/11] Minor change - make prewhitening correction depend only on keyword, not all decimation params --- aurora/pipelines/time_series_helpers.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/aurora/pipelines/time_series_helpers.py b/aurora/pipelines/time_series_helpers.py index a2c0efb9..b1fd1356 100644 --- a/aurora/pipelines/time_series_helpers.py +++ b/aurora/pipelines/time_series_helpers.py @@ -117,13 +117,13 @@ def apply_recoloring( return stft_obj if decimation_obj.prewhitening_type == "first difference": - freqs = decimation_obj.fft_frequencies - prewhitening_correction = 1.0j * 2 * np.pi * freqs # jw - - stft_obj /= prewhitening_correction + # first difference prewhitening correction is to divide by jw + freqs = stft_obj.frequency.data # was freqs = decimation_obj.fft_frequencies + jw = 1.0j * 2 * np.pi * freqs + stft_obj /= jw # suppress nan and inf to mute later warnings - if prewhitening_correction[0] == 0.0: + if jw[0] == 0.0: cond = stft_obj.frequency != 0.0 stft_obj = stft_obj.where(cond, complex(0.0)) # elif decimation_obj.prewhitening_type == "ARMA": From 080709cd81f218e878a1011afdeeebc4da8e1f88 Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 17:38:16 -0800 Subject: [PATCH 09/11] doc / type hints --- aurora/pipelines/time_series_helpers.py | 8 +++++++- tests/synthetic/test_fourier_coefficients.py | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/aurora/pipelines/time_series_helpers.py b/aurora/pipelines/time_series_helpers.py index b1fd1356..8c3627e5 100644 --- a/aurora/pipelines/time_series_helpers.py +++ b/aurora/pipelines/time_series_helpers.py @@ -139,15 +139,21 @@ def apply_recoloring( return stft_obj -def run_ts_to_stft_scipy(decimation_obj: FCDecimation, run_xrds_orig: xr.Dataset): +def run_ts_to_stft_scipy( + decimation_obj: Union[AuroraDecimationLevel, FCDecimation], + run_xrds_orig: xr.Dataset, +) -> xr.Dataset: """ Converts a runts object into a time series of Fourier coefficients. This method uses scipy.signal.spectrogram. + Parameters ---------- decimation_obj : mt_metadata.transfer_functions.processing.aurora.DecimationLevel Information about how the decimation level is to be processed + Note: This works with FCdecimation and AuroraDecimationLevel becuase test_fourier_coefficients + and test_stft_methods_agree both use them) run_xrds_orig : : xarray.core.dataset.Dataset Time series to be processed diff --git a/tests/synthetic/test_fourier_coefficients.py b/tests/synthetic/test_fourier_coefficients.py index c8b8199a..8b89ff4d 100644 --- a/tests/synthetic/test_fourier_coefficients.py +++ b/tests/synthetic/test_fourier_coefficients.py @@ -115,6 +115,7 @@ def test_123(self): x.to_fc_decimation() for x in processing_config.decimations ] # For code coverage, have a case where fc_decimations is None + # This also (indirectly) tests a different FCDeecimation object. if mth5_path.stem == "test1": fc_decimations = None From 021ad304ab13ab2bf675ac665319c6ba35fa3c4c Mon Sep 17 00:00:00 2001 From: kkappler Date: Fri, 20 Dec 2024 17:47:03 -0800 Subject: [PATCH 10/11] Move recoloring boolean check above call to apply_recoloring --- aurora/pipelines/time_series_helpers.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/aurora/pipelines/time_series_helpers.py b/aurora/pipelines/time_series_helpers.py index 8c3627e5..7d120ead 100644 --- a/aurora/pipelines/time_series_helpers.py +++ b/aurora/pipelines/time_series_helpers.py @@ -113,8 +113,9 @@ def apply_recoloring( # No recoloring needed if prewhitening not appiled, or recoloring set to False if not decimation_obj.prewhitening_type: return stft_obj - if not decimation_obj.recoloring: - return stft_obj + # TODO Delete after tests (20241220) -- this check has been moved above the call to this function + # if not decimation_obj.recoloring: + # return stft_obj if decimation_obj.prewhitening_type == "first difference": # first difference prewhitening correction is to divide by jw @@ -194,8 +195,8 @@ def run_ts_to_stft_scipy( coords={"frequency": ff, "time": time_axis}, ) stft_obj.update({channel_id: xrd}) - - stft_obj = apply_recoloring(decimation_obj, stft_obj) + if decimation_obj.recoloring: + stft_obj = apply_recoloring(decimation_obj, stft_obj) return stft_obj @@ -318,7 +319,8 @@ def run_ts_to_stft( detrend_type=decimation_obj.extra_pre_fft_detrend_type, ) - stft_obj = apply_recoloring(decimation_obj, stft_obj) + if decimation_obj.recoloring: + stft_obj = apply_recoloring(decimation_obj, stft_obj) return stft_obj From b4af28becbcaec1f7c967f9283caa44b13b2461a Mon Sep 17 00:00:00 2001 From: kkappler Date: Thu, 26 Dec 2024 09:50:09 -0800 Subject: [PATCH 11/11] tidy doc --- aurora/pipelines/fourier_coefficients.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/aurora/pipelines/fourier_coefficients.py b/aurora/pipelines/fourier_coefficients.py index d67acb7a..2eb31a2a 100644 --- a/aurora/pipelines/fourier_coefficients.py +++ b/aurora/pipelines/fourier_coefficients.py @@ -264,7 +264,7 @@ def get_degenerate_fc_decimation(sample_rate: float) -> list: Makes a default fc_decimation list. WIP This "degnerate" config will only operate on the first decimation level. - This is useful for testing but could be used in future if an MTH5 stored time series in decimation + This is useful for testing. It could also be used in future on an MTH5 stored time series in decimation levels already as separate runs. Parameters @@ -275,7 +275,8 @@ def get_degenerate_fc_decimation(sample_rate: float) -> list: Returns ------- output: list - List has only one element which is of type mt_metadata.transfer_functions.processing.fourier_coefficients.Decimation. + List has only one element which is of type FCDecimation, aka. + mt_metadata.transfer_functions.processing.fourier_coefficients.Decimation. """ output = fc_decimations_creator( sample_rate,