From 966828e1285cb2f23cb00550e85f0cdec23c6bea Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Nov 2024 22:35:52 -0700 Subject: [PATCH 1/9] Refactored graphics/waveforms.py --- mtuq/graphics/waveforms.py | 458 ++++++++++--------------------------- 1 file changed, 123 insertions(+), 335 deletions(-) diff --git a/mtuq/graphics/waveforms.py b/mtuq/graphics/waveforms.py index 86bfe2e9..a44a053f 100644 --- a/mtuq/graphics/waveforms.py +++ b/mtuq/graphics/waveforms.py @@ -6,7 +6,6 @@ import numpy as np import matplotlib.pyplot as pyplot -from collections import defaultdict from matplotlib.font_manager import FontProperties from mtuq.dataset import Dataset from mtuq.event import MomentTensor, Force @@ -22,7 +21,8 @@ # high-level plotting functions # -def plot_waveforms1(filename, +def plot_waveforms1( + filename, data, synthetics, stations, @@ -60,18 +60,14 @@ def plot_waveforms1(filename, _add_component_labels1(axes) - max_amplitude = _max(data, synthetics) - - if normalize == 'median_amplitude': - # Using the updated _median_amplitude function to calculate the median of non-zero maximum amplitudes - max_amplitude_median = _median_amplitude(data, synthetics) - max_amplitudes = np.array([max_amplitude_median if len(data[i]) > 0 and len(synthetics[i]) > 0 else 0.0 for i in range(len(data))]) - elif normalize == 'maximum_amplitude': - max_amplitudes = np.array([max_amplitude if len(data[i]) > 0 and len(synthetics[i]) > 0 else 0.0 for i in range(len(data))]) - elif normalize == 'station_amplitude' or normalize == 'trace_amplitude': - pass + # scaling factor applied to all traces + if normalize=='maximum_amplitude': + factor = _max(data, synthetics) + elif normalize=='median_amplitude': + factor = 2.*_median(data, synthetics) else: - raise ValueError("Invalid normalization method specified.") + factor = None + # # loop over stations @@ -89,29 +85,10 @@ def plot_waveforms1(filename, if station_label_writer is not None: station_label_writer(axes[ir][0], stations[_i], origin) - # # plot traces - # - - stream_dat = data[_i] - stream_syn = synthetics[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZRT(axes[ir], 1, dat, syn, component, - normalize, trace_label_writer, max_amplitudes[_i], total_misfit) + _plot_stream(axes[ir], [1,2,3], ['Z','R','T'], + data[_i], synthetics[_i], + normalize, factor, trace_label_writer, total_misfit) ir += 1 @@ -120,7 +97,8 @@ def plot_waveforms1(filename, -def plot_waveforms2(filename, +def plot_waveforms2( + filename, data_bw, data_sw, synthetics_bw, @@ -159,28 +137,16 @@ def plot_waveforms2(filename, _add_component_labels2(axes) - # determine maximum trace amplitudes - max_amplitude_bw = _max(data_bw, synthetics_bw) - max_amplitude_sw = _max(data_sw, synthetics_sw) - - - if normalize == 'median_amplitude': - # For body wave data and synthetics - bw_median = _median_amplitude(data_bw, synthetics_bw) - max_amplitudes_bw = np.array([bw_median if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - - # For surface wave data and synthetics - sw_median = _median_amplitude(data_sw, synthetics_sw) - max_amplitudes_sw = np.array([sw_median if len(data_sw[i]) > 0 and len(synthetics_sw[i]) > 0 else 0.0 for i in range(len(data_sw))]) - elif normalize == 'maximum_amplitude': - max_amplitudes_bw = np.array([max_amplitude_bw if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - max_amplitudes_sw = np.array([max_amplitude_sw if len(data_sw[i]) > 0 and len(synthetics_sw[i]) > 0 else 0.0 for i in range(len(data_sw))]) - elif normalize == 'station_amplitude' or normalize == 'trace_amplitude': - max_amplitudes_bw = np.array([_max(data_bw[i], synthetics_bw[i]) if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - max_amplitudes_sw = np.array([_max(data_sw[i], synthetics_sw[i]) if len(data_sw[i]) > 0 and len(synthetics_sw[i]) > 0 else 0.0 for i in range(len(data_sw))]) + if normalize=='maximum_amplitude': + factor_bw = _max(data_bw, synthetics_bw) + factor_sw = _max(data_sw, synthetics_sw) + elif normalize=='median_amplitude': + factor_bw = 2.*_median(data_bw, synthetics_bw) + factor_sw = 2.*_median(data_sw, synthetics_sw) else: - raise ValueError("Invalid normalization method specified.") - + factor_bw = None + factor_sw = None + # # loop over stations # @@ -194,58 +160,20 @@ def plot_waveforms2(filename, continue # add station labels - if station_label_writer is not None: + try: station_label_writer(axes[ir][0], stations[_i], origin) + except: + pass - # # plot body wave traces - # - - stream_dat = data_bw[_i] - stream_syn = synthetics_bw[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZR(axes[ir], 1, dat, syn, component, - normalize, trace_label_writer, max_amplitudes_bw[_i], total_misfit_bw) - + _plot_stream(axes[ir], [1,2], ['Z','R'], + data_bw[_i], synthetics_bw[_i], + normalize, factor_bw, trace_label_writer, total_misfit_bw) - # # plot surface wave traces - # - - stream_dat = data_sw[_i] - stream_syn = synthetics_sw[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZRT(axes[ir], 3, dat, syn, component, - normalize, trace_label_writer, max_amplitudes_sw[_i], total_misfit_sw) - + _plot_stream(axes[ir], [3,4,5], ['Z','R','T'], + data_sw[_i], synthetics_sw[_i], + normalize, factor_sw, trace_label_writer, total_misfit_sw) ir += 1 @@ -253,7 +181,8 @@ def plot_waveforms2(filename, pyplot.close() -def plot_waveforms3(filename, +def plot_waveforms3( + filename, data_bw, data_rayleigh, data_love, @@ -295,35 +224,21 @@ def plot_waveforms3(filename, _add_component_labels2(axes) - # determine maximum trace amplitudes - max_amplitude_bw = _max(data_bw, synthetics_bw) - max_amplitude_rayleigh = _max(data_rayleigh, synthetics_rayleigh) - max_amplitude_love = _max(data_love, synthetics_love) - - if normalize == 'median_amplitude': - # For body wave data and synthetics - bw_median = _median_amplitude(data_bw, synthetics_bw) - max_amplitudes_bw = np.array([bw_median if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - - # For Rayleigh wave data and synthetics - rayleigh_median = _median_amplitude(data_rayleigh, synthetics_rayleigh) - max_amplitudes_rayleigh = np.array([rayleigh_median if len(data_rayleigh[i]) > 0 and len(synthetics_rayleigh[i]) > 0 else 0.0 for i in range(len(data_rayleigh))]) - - # For Love wave data and synthetics - love_median = _median_amplitude(data_love, synthetics_love) - max_amplitudes_love = np.array([love_median if len(data_love[i]) > 0 and len(synthetics_love[i]) > 0 else 0.0 for i in range(len(data_love))]) - - elif normalize == 'maximum_amplitude': - max_amplitudes_bw = np.array([max_amplitude_bw if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - max_amplitudes_rayleigh = np.array([max_amplitude_rayleigh if len(data_rayleigh[i]) > 0 and len(synthetics_rayleigh[i]) > 0 else 0.0 for i in range(len(data_rayleigh))]) - max_amplitudes_love = np.array([max_amplitude_love if len(data_love[i]) > 0 and len(synthetics_love[i]) > 0 else 0.0 for i in range(len(data_love))]) - - elif normalize == 'station_amplitude' or normalize == 'trace_amplitude': - max_amplitudes_bw = np.array([_max(data_bw[i], synthetics_bw[i]) if len(data_bw[i]) > 0 and len(synthetics_bw[i]) > 0 else 0.0 for i in range(len(data_bw))]) - max_amplitudes_rayleigh = np.array([_max(data_rayleigh[i], synthetics_rayleigh[i]) if len(data_rayleigh[i]) > 0 and len(synthetics_rayleigh[i]) > 0 else 0.0 for i in range(len(data_rayleigh))]) - max_amplitudes_love = np.array([_max(data_love[i], synthetics_love[i]) if len(data_love[i]) > 0 and len(synthetics_love[i]) > 0 else 0.0 for i in range(len(data_love))]) + if normalize=='maximum_amplitude': + factor_bw = _max(data_bw, synthetics_bw) + factor_rayleigh = _max(data_rayleigh, synthetics_rayleigh) + factor_love = _max(data_love, synthetics_love) + + elif normalize=='median_amplitude': + factor_bw = 2.*_median(data_bw, synthetics_bw) + factor_rayleigh = 2.*_median(data_rayleigh, synthetics_rayleigh) + factor_love = 2.*_median(data_love, synthetics_love) + else: - raise ValueError("Invalid normalization method specified.") + factor_bw = None + factor_rayeligh = None + factor_love = None + # # loop over stations @@ -338,81 +253,25 @@ def plot_waveforms3(filename, continue # add station labels - if station_label_writer is not None: + try: station_label_writer(axes[ir][0], stations[_i], origin) + except: + pass - # - # plot body wave traces - # - - stream_dat = data_bw[_i] - stream_syn = synthetics_bw[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZR(axes[ir], 1, dat, syn, component, - normalize, trace_label_writer, max_amplitudes_bw[_i], total_misfit_bw) - - # - # plot Rayleigh wave traces - # - - stream_dat = data_rayleigh[_i] - stream_syn = synthetics_rayleigh[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZRT(axes[ir], 3, dat, syn, component, - normalize, trace_label_writer, max_amplitudes_rayleigh[_i], total_misfit_rayleigh) - - # - # plot Love wave traces - # - - stream_dat = data_love[_i] - stream_syn = synthetics_love[_i] - - for dat in stream_dat: - component = dat.stats.channel[-1].upper() - weight = _getattr(dat, 'weight', 1.) - - if not weight: - continue - - # skip missing components - try: - syn = stream_syn.select(component=component)[0] - except: - warn('Missing component, skipping...') - continue - - _plot_ZRT(axes[ir], 3, dat, syn, component, - normalize, trace_label_writer, max_amplitudes_love[_i], total_misfit_love) + # plot body waves + _plot_stream(axes[ir], [1,2], ['Z','R'], + data_bw[_i], synthetics_bw[_i], + normalize, factor_bw, trace_label_writer, total_misfit_bw) + + # plot Rayleigh waves + _plot_stream(axes[ir], [3,4], ['Z','R'], + data_rayleigh[_i], synthetics_rayeleigh[_i], + normalize, factor_rayleigh, trace_label_writer, total_misfit_rayleigh) + # plot Love waves + _plot_stream(axes[ir], [5], ['T'], + data_love[_i], synthetics_love[_i], + normalize, factor_love, trace_label_writer, total_misfit_love) ir += 1 @@ -420,7 +279,8 @@ def plot_waveforms3(filename, pyplot.close() -def plot_data_greens1(filename, +def plot_data_greens1( + filename, data, greens, process_data, @@ -439,9 +299,6 @@ def plot_data_greens1(filename, # collect synthetic waveforms with misfit attributes attached synthetics = misfit.collect_synthetics(data, greens.select(origin), source) - # Prune synthetics to only include components that are present in the misfit time_shift_groups - _prune_synthetics(misfit, synthetics) - # calculate total misfit for display in figure header total_misfit = misfit(data, greens.select(origin), source, optimization_level=0) @@ -490,9 +347,6 @@ def plot_data_greens2(filename, synthetics_sw = misfit_sw.collect_synthetics( data_sw, greens_sw.select(origin), source) - # Prune synthetics to only include components that are present in the misfit time_shift_groups - _prune_synthetics(misfit_bw, synthetics_bw) - _prune_synthetics(misfit_sw, synthetics_sw) # calculate total misfit for display in figure header total_misfit_bw = misfit_bw( @@ -522,7 +376,8 @@ def plot_data_greens2(filename, header=header, **kwargs) -def plot_data_greens3(filename, +def plot_data_greens3( + filename, data_bw, data_rayleigh, data_love, @@ -544,7 +399,7 @@ def plot_data_greens3(filename, """ Creates data/synthetics comparison figure with 5 columns (Pn Z, Pn R, Rayleigh Z, Rayleigh R, Love T) - Different input arguments, same result as plot_waveforms3. Considers three data groups. + Different input arguments, same result as plot_waveforms3 """ # collect synthetic waveforms with misfit attributes attached @@ -557,11 +412,6 @@ def plot_data_greens3(filename, synthetics_love = misfit_love.collect_synthetics( data_love, greens_love.select(origin), source) - # Prune synthetics to only include components that are present in the misfit time_shift_groups - _prune_synthetics(misfit_bw, synthetics_bw) - _prune_synthetics(misfit_rayleigh, synthetics_rayleigh) - _prune_synthetics(misfit_love, synthetics_love) - # calculate total misfit for display in figure header total_misfit_bw = misfit_bw( data_bw, greens_bw.select(origin), source, optimization_level=0) @@ -570,10 +420,7 @@ def plot_data_greens3(filename, data_rayleigh, greens_rayleigh.select(origin), source, optimization_level=0) total_misfit_love = misfit_love( - data_love, greens_love.select(origin), source, optimization_level=0) - - # Print the length of the data and synthetics for each group - # Hav + data_love, greens_love.select(origin), source, optimization_level=0) # prepare figure header if 'header' in kwargs: @@ -593,7 +440,9 @@ def plot_data_greens3(filename, data_bw, data_rayleigh, data_love, synthetics_bw, synthetics_rayleigh, synthetics_love, stations, origin, - total_misfit_bw=total_misfit_bw, total_misfit_rayleigh=total_misfit_rayleigh, total_misfit_love=total_misfit_love, + total_misfit_bw=total_misfit_bw, + total_misfit_rayleigh=total_misfit_rayleigh, + total_misfit_love=total_misfit_love, header=header, **kwargs) @@ -647,60 +496,50 @@ def _initialize(nrows=None, ncolumns=None, column_width_ratios=None, return fig, axes -def _plot_ZRT(axes, ic, dat, syn, component, - normalize='maximum_amplitude', trace_label_writer=None, - normalization_amplitude=1., total_misfit=1.): - - # plot traces - if component=='Z': - axis = axes[ic+0] - elif component=='R': - axis = axes[ic+1] - elif component=='T': - axis = axes[ic+2] - else: - return - - _plot(axis, dat, syn) - - # normalize amplitude -- logic for station_amplitude, median_amplitude, and maximum_amplitude is done at higher level - if normalize=='trace_amplitude': - max_trace = _max(dat, syn) - ylim = [-1.5*max_trace, +1.5*max_trace] - axis.set_ylim(*ylim) - elif normalize=='station_amplitude' or normalize=='median_amplitude' or normalize=='maximum_amplitude': - ylim = [-1.25*normalization_amplitude, +1.25*normalization_amplitude] - axis.set_ylim(*ylim) - - if trace_label_writer is not None: - trace_label_writer(axis, dat, syn, total_misfit) +def _plot_stream( + axes, + column_indices, + components, + stream_dat, + stream_syn, + normalize='maximum_amplitude', + amplitude_factor=None, + trace_label_writer=None, + total_misfit=1. + ): + + for _i, component, in enumerate(components): + axis = axes[column_indices[_i]] + + try: + dat = stream_dat.select(component=component)[0] + syn = stream_syn.select(component=component)[0] + except: + continue + weight = _getattr(dat, 'weight', 1.) + if not weight: + continue -def _plot_ZR(axes, ic, dat, syn, component, - normalize='maximum_amplitude', trace_label_writer=None, - normalization_amplitude=1., total_misfit=1.): + _plot(axis, dat, syn) - # plot traces - if component=='Z': - axis = axes[ic+0] - elif component=='R': - axis = axes[ic+1] - else: - return + if normalize=='trace_amplitude': + max_trace = _max(dat, syn) + ylim = [-1.5*max_trace, +1.5*max_trace] + axis.set_ylim(*ylim) - _plot(axis, dat, syn) + elif normalize=='station_amplitude': + max_station = _max(stream_dat, stream_syn) + ylim = [-1.25*max_station, +1.25*max_station] - # normalize amplitude -- logic for station_amplitude, median_amplitude, and maximum_amplitude is done at higher level - if normalize=='trace_amplitude': - max_trace = _max(dat, syn) - ylim = [-1.5*max_trace, +1.5*max_trace] - axis.set_ylim(*ylim) - elif normalize=='station_amplitude' or normalize=='median_amplitude' or normalize=='maximum_amplitude': - ylim = [-1.25*normalization_amplitude, +1.25*normalization_amplitude] - axis.set_ylim(*ylim) + elif amplitude_factor: + ylim = [-amplitude_factor, +amplitude_factor] + axis.set_ylim(*ylim) - if trace_label_writer is not None: - trace_label_writer(axis, dat, syn, total_misfit) + try: + trace_label_writer(axis, dat, syn, total_misfit) + except: + pass def _plot(axis, dat, syn, label=None): @@ -801,19 +640,8 @@ def _isempty(dataset): def _max(dat, syn): - """ - Computes the maximum value from a set of two input data objects (observed and synthetics). - - Parameters: - dat (Trace, Stream, or Dataset): observed data. - syn (Trace, Stream, or Dataset): synthetics. + # returns maximum amplitude over traces, streams, or datasets - Returns: - float: The maximum value for normalization purposes. - - Raises: - TypeError: If the input objects are not of the same type (Trace, Stream, or Dataset). - """ if type(dat)==type(syn)==Trace: return max( abs(dat.max()), @@ -832,43 +660,18 @@ def _max(dat, syn): else: raise TypeError -def _median_amplitude(data, synthetics): - """ - Computes the median of the maximum non-zero amplitudes for pairs of data and synthetic traces. - - Args: - data: A list of of observed data (can be Trace, Stream, or Dataset objects). - synthetics: A list of synthetic traces corresponding to the observed data. +def _median(*datasets): + # returns median Linf amplitude over traces in dataset - Returns: - The median of the non-zero maximum amplitudes computed across all pairs. - - Raises: - ValueError: If the lengths of data and synthetics lists differ. - """ - # Validate input lengths - # If Trace directly input, make it a list - data = [data] if isinstance(data, Trace) else data - synthetics = [synthetics] if isinstance(synthetics, Trace) else synthetics - - # Validate lengths - if len(data) != len(synthetics): - raise ValueError("Data and synthetics lists must have the same length.") - - max_amplitudes = [] - - # Iterate over pairs and handle empty traces - This gets a list of maximum amplitudes for each pair of data and synthetics - for dat, syn in zip(data, synthetics): - if not dat or not syn: - max_amplitudes.append(0) - else: - max_amplitudes.append(_max(dat, syn)) - - # Convert to NumPy array for efficient filtering - max_amplitudes = np.array(max_amplitudes) - - # Compute median of non-zero values or return 0 if none exist - return np.median(max_amplitudes[max_amplitudes > 0]) if np.any(max_amplitudes > 0) else 0.0 + _list = [] + for ds in datasets: + assert type(ds)==Dataset + for stream in ds: + for trace in stream: + trace_max = abs(trace.max()) + if trace_max > 0.: + _list.append(trace_max) + return np.median(np.array(_list)) def _hide_axes(axes): @@ -914,18 +717,3 @@ def _get_tag(tags, pattern): else: return None -def _prune_synthetics(misfit, synthetics): - """ Prunes synthetics to only include components that are present in the misfit time_shift_groups - - note: Ideally, this would be done internally when collecting synthetics in the misfit object - but it is done here to avoid modifying important moving parts of the codebase. - """ - components = [] - for group in misfit.time_shift_groups: - for component in group: - components.append(component) - - for sta in synthetics: - for tr in sta: - if tr.stats.channel not in components: - sta.remove(tr) \ No newline at end of file From 7109e97531df59fd1dfa588c5d377043ad91d523 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Nov 2024 22:37:40 -0700 Subject: [PATCH 2/9] collect_synthetics() options for all stations and components, or certain subsets --- mtuq/misfit/waveform/__init__.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/mtuq/misfit/waveform/__init__.py b/mtuq/misfit/waveform/__init__.py index 4070ad68..6f585e54 100644 --- a/mtuq/misfit/waveform/__init__.py +++ b/mtuq/misfit/waveform/__init__.py @@ -5,7 +5,7 @@ from mtuq.misfit.waveform import level0, level1, level2 from mtuq.misfit.waveform._stats import estimate_sigma, calculate_norm_data from mtuq.util import Null, iterable, warn -from mtuq.util.math import isclose, list_intersect_with_indices +from mtuq.util.math import isclose, list_intersect, list_intersect_with_indices from mtuq.util.signal import check_padding, get_components, isempty @@ -247,9 +247,32 @@ def collect_attributes(self, data, greens, source, normalize=False): return deepcopy(attrs) - def collect_synthetics(self, data, greens, source, normalize=False): + def collect_synthetics(self, data, greens, source, normalize=False, mode=2): """ Collects synthetics with misfit, time shifts and other attributes attached """ + + if mode==1: + # returns synthetics from all stations and components + components = [['Z','R','T']]*len(data) + + elif mode==2: + # returns synthetics only from stations and components which exist in + # observed data + components = data.get_components() + + elif mode==3: + # returns synthetics only from stations and components which exist + # and are included in misfit function evaluation + _active = [] + for group in self.time_shift_groups: + for component in group: + _active.append(component) + + components = [] + for _exist in data.get_components(): + components.append(list_intersect(_active, _exist)) + + # checks that dataset is nonempty if isempty(data): warn("Empty data set. No attributes will be returned") @@ -259,8 +282,9 @@ def collect_synthetics(self, data, greens, source, normalize=False): # shift bounds check_padding(greens, self.time_shift_min, self.time_shift_max) + synthetics = greens.get_synthetics( - source, components=data.get_components(), stats=data.get_stats(), + source, components=components, stats=data.get_stats(), mode='map', inplace=True) # attaches attributes to synthetics @@ -269,6 +293,7 @@ def collect_synthetics(self, data, greens, source, normalize=False): self.time_shift_min, self.time_shift_max, msg_handle=Null(), normalize=False, set_attributes=True) + return deepcopy(synthetics) From 8aefdf5a09adf2717a324ced3a32cfd28a20ea5f Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Nov 2024 23:22:07 -0700 Subject: [PATCH 3/9] Less code duplication --- mtuq/graphics/waveforms.py | 101 +++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 54 deletions(-) diff --git a/mtuq/graphics/waveforms.py b/mtuq/graphics/waveforms.py index a44a053f..84da822c 100644 --- a/mtuq/graphics/waveforms.py +++ b/mtuq/graphics/waveforms.py @@ -40,27 +40,14 @@ def plot_waveforms1( raise Exception # how many stations have at least one trace? - nstations = _count([data]) + nrows = _count([data]) - # - # initialize figure - # - - fig, axes = _initialize( - nrows=nstations, - ncolumns=4, - column_width_ratios=[1.,1.,1.], - height=1.25*nstations, - width=8.5, - margin_right=0.5, - header=header, - header_height=1.5, - station_labels=bool(station_label_writer), - ) + # intialize figure + fig, axes = _initialize1(nrows, header) _add_component_labels1(axes) - # scaling factor applied to all traces + # optional normalization if normalize=='maximum_amplitude': factor = _max(data, synthetics) elif normalize=='median_amplitude': @@ -117,26 +104,16 @@ def plot_waveforms2( """ Creates data/synthetics comparison figure with 5 columns (Pn Z, Pn R, Rayleigh Z, Rayleigh R, Love T) """ - # how many stations have at least one trace? - nstations = _count([data_bw, data_sw]) - # - # initialize figure - # + # how many stations have at least one trace? + nrows = _count([data_bw, data_sw]) - fig, axes = _initialize( - nrows=nstations, - ncolumns=6, - column_width_ratios=[0.5,0.5,1.,1.,1.], - height=1.25*nstations, - width=10., - header=header, - header_height=2., - station_labels=bool(station_label_writer), - ) + # intialize figure + fig, axes = _initialize2(nrows, header) _add_component_labels2(axes) + # optional normalization if normalize=='maximum_amplitude': factor_bw = _max(data_bw, synthetics_bw) factor_sw = _max(data_sw, synthetics_sw) @@ -201,29 +178,19 @@ def plot_waveforms3( ): """ Creates data/synthetics comparison figure with 5 columns - (Pn Z, Pn R, Rayleigh Z, Rayleigh R, Love T), for three data groups + (Pn Z, Pn R, Rayleigh Z, Rayleigh R, Love T) """ # how many stations have at least one trace? - nstations = _count([data_bw, data_rayleigh, data_love]) + nrows = _count([data_bw, data_rayleigh, data_love]) - # - # initialize figure - # - - fig, axes = _initialize( - nrows=nstations, - ncolumns=6, - column_width_ratios=[0.5,0.5,1.,1.,1.], - height=1.25*nstations, - width=10., - header=header, - header_height=2., - station_labels=bool(station_label_writer), - ) + # intialize figure + fig, axes = _initialize2(nrows, header) _add_component_labels2(axes) + + # optional normalization if normalize=='maximum_amplitude': factor_bw = _max(data_bw, synthetics_bw) factor_rayleigh = _max(data_rayleigh, synthetics_rayleigh) @@ -450,8 +417,35 @@ def plot_data_greens3( # low-level plotting utilities # +def _initialize1(nrows, header): + return _initialize( + nrows=nrows, + ncolumns=4, + column_width_ratios=[1.,1.,1.], + height=1.25*nrows, + width=8.5, + margin_right=0.5, + header=header, + header_height=1.5, + station_labels=True, + ) + + +def _initialize2(nrows, header): + return _initialize( + nrows=nrows, + ncolumns=6, + column_width_ratios=[0.5,0.5,1.,1.,1.], + height=1.25*nrows, + width=10., + header=header, + header_height=2., + station_labels=True, + ) + -def _initialize(nrows=None, ncolumns=None, column_width_ratios=None, +def _initialize( + nrows=None, ncolumns=None, column_width_ratios=None, header=None, height=None, width=None, margin_top=0.25, margin_bottom=0.25, margin_left=0.25, margin_right=0.25, header_height=1.5, station_labels=True, station_label_width=0.5): @@ -521,12 +515,9 @@ def _plot_stream( if not weight: continue - _plot(axis, dat, syn) - if normalize=='trace_amplitude': max_trace = _max(dat, syn) ylim = [-1.5*max_trace, +1.5*max_trace] - axis.set_ylim(*ylim) elif normalize=='station_amplitude': max_station = _max(stream_dat, stream_syn) @@ -534,7 +525,9 @@ def _plot_stream( elif amplitude_factor: ylim = [-amplitude_factor, +amplitude_factor] - axis.set_ylim(*ylim) + + _plot_trace(axis, dat, syn) + axis.set_ylim(*ylim) try: trace_label_writer(axis, dat, syn, total_misfit) @@ -542,7 +535,7 @@ def _plot_stream( pass -def _plot(axis, dat, syn, label=None): +def _plot_trace(axis, dat, syn, label=None): """ Plots data and synthetics time series on current axes """ t1,t2,nt,dt = _time_stats(dat) From 282cea55a7c8522cc4e1ab42548c243076599a3a Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 18 Nov 2024 13:48:35 -0700 Subject: [PATCH 4/9] Use raw strings in code generator --- setup/code_generator.py | 46 ++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/setup/code_generator.py b/setup/code_generator.py index b81eb03d..59f3b88f 100755 --- a/setup/code_generator.py +++ b/setup/code_generator.py @@ -197,7 +197,7 @@ # import matplotlib - matplotlib.use('Agg', warn=False, force=True) + matplotlib.use('Agg', force=True) import matplotlib """ @@ -909,7 +909,7 @@ print('Reading data...\\n') data = read(path_data, format='sac', event_id=event_id, - station_id_list=station_id_list, + #station_id_list=station_id_list, tags=['units:m', 'type:velocity']) @@ -1833,27 +1833,27 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): file.write( replace( Main_GridSearch, - 'origin', - 'origins', - 'Reading Greens functions...\\\\n', + r'origin', + r'origins', + r'Reading Greens functions...\\\\n', ( - 'Reading Greens functions...\\\\n\\\\n'+ - ' Downloads can sometimes take as long as a few hours!\\\\n' + r'Reading Greens functions...\\\\n\\\\n'+ + r' Downloads can sometimes take as long as a few hours!\\\\n' ), - 'download_greens_tensors\(stations, origin, model\)', - 'download_greens_tensors(stations, origin, model, verbose=True)', + r'download_greens_tensors\(stations, origin, model\)', + r'download_greens_tensors(stations, origin, model, verbose=True)', )) file.write( replace( WrapUp_GridSearch_DoubleCoupleMagnitudeDepth, - 'DC\+Z', - 'DC+XY', - 'misfit_depth', - 'misfit_latlon', - "title=event_id", - "title=event_id, colorbar_label='L2 misfit'", - 'show_magnitudes=True, ', - '', + r'DC\+Z', + r'DC+XY', + r'misfit_depth', + r'misfit_latlon', + r"title=event_id", + r"title=event_id, colorbar_label='L2 misfit'", + r'show_magnitudes=True, ', + r'', )) @@ -2005,9 +2005,9 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): file.write( replace( Main1_SerialGridSearch_DoubleCouple, - 'greens = download_greens_tensors\(stations, origin, model\)', - 'db = open_db(path_greens, format=\'FK\', model=model)\n ' - +'greens = db.get_greens_tensors(stations, origin)', + r'greens = download_greens_tensors\(stations, origin, model\)', + r'db = open_db(path_greens, format=\'FK\', model=model)\n ' + +r'greens = db.get_greens_tensors(stations, origin)', )) file.write( replace( @@ -2071,9 +2071,9 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): file.write( replace( Main1_SerialGridSearch_DoubleCouple, - 'greens = download_greens_tensors\(stations, origin, model\)', - 'db = open_db(path_greens, format=\'FK\', model=model)\n ' - +'greens = db.get_greens_tensors(stations, origin)', + r'greens = download_greens_tensors\(stations, origin, model\)', + r'db = open_db(path_greens, format=\'FK\', model=model)\n ' + +r'greens = db.get_greens_tensors(stations, origin)', )) file.write(Main_TestMisfit) From 79fae332cde0e3ffa97fe2471a25b65f80d6486b Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 18 Nov 2024 15:54:40 -0700 Subject: [PATCH 5/9] Updated test_graphics.py --- setup/code_generator.py | 70 +++++++++++++++++++++++++++-------- tests/test_graphics.py | 81 ++++++++++++++++++++++++++++++++--------- 2 files changed, 119 insertions(+), 32 deletions(-) diff --git a/setup/code_generator.py b/setup/code_generator.py index 59f3b88f..7b98172b 100755 --- a/setup/code_generator.py +++ b/setup/code_generator.py @@ -191,9 +191,10 @@ # # Tests data, synthetics and beachball plotting utilities # - # Note that in the figures created by this script, the data and synthetics - # are not expected to fit epsecially well; currently, the only requirement - # is that the script runs without errors + + # + # The idea is for a test that runs very quickly, suitable for CI testing; + # eventually we may more detailed tests to tests/graphics # import matplotlib @@ -703,8 +704,13 @@ Grid_TestGraphics=""" - mt = np.sqrt(1./3.)*np.array([1., 1., 1., 0., 0., 0.]) # explosion - mt *= 1.e16 + from mtuq import MomentTensor + + mt = MomentTensor( + 1.e16 * np.sqrt(1./3.)*np.array([1., 1., 1., 0., 0., 0.])) # explosion + + mt_dict = { + 'rho':1.,'v':0.,'w':3/8*np.pi,'kappa':0.,'sigma':0.,'h':0.} wavelet = Trapezoid( magnitude=4.5) @@ -935,23 +941,56 @@ # Generate figures # - print('Figure 1 of 3\\n') + print('Plot data (1 of 6)\\n') - plot_data_greens2('graphics_test_1.png', - data_bw, data_sw, greens_bw, greens_sw, - process_bw, process_sw, misfit_bw, misfit_sw, - stations, origin, mt, header=False) + from mtuq.graphics import plot_waveforms2 + from mtuq.util import Null + + plot_waveforms2('graphics_test_1.png', + data_bw, data_sw, Null(), Null(), + stations, origin, header=False) + + + print('Plot synthetics (2 of 6)\\n') + + synthetics_bw = greens_bw.get_synthetics(mt, components=['Z','R']) + synthetics_sw = greens_sw.get_synthetics(mt, components=['Z','R','T']) + + + plot_waveforms2('graphics_test_2.png', + synthetics_bw, synthetics_sw, Null(), Null(), + stations, origin, header=False) - print('Figure 2 of 3\\n') - plot_data_greens2('graphics_test_2.png', + print('Plot synthetics (3 of 6)\\n') + + synthetics_bw = misfit_bw.collect_synthetics(data_bw, greens_bw, mt) + synthetics_sw = misfit_sw.collect_synthetics(data_sw, greens_sw, mt) + + + plot_waveforms2('graphics_test_3.png', + synthetics_bw, synthetics_sw, Null(), Null(), + stations, origin, header=False) + + + print('Plot data and synthetics without header (4 of 6)\\n') + + plot_waveforms2('graphics_test_4.png', + data_bw, data_sw, synthetics_bw, synthetics_sw, + stations, origin, header=False) + + + print('Plot data and synthetics without header (5 of 6)\\n') + + plot_data_greens2('graphics_test_5.png', data_bw, data_sw, greens_bw, greens_sw, process_bw, process_sw, misfit_bw, misfit_sw, - stations, origin, mt, header=False) + stations, origin, mt, mt_dict) - print('Figure 3 of 3\\n') - plot_beachball('graphics_test_3.png', + print('Plot explosion bechball (6 of 6)\\n') + + plot_beachball('graphics_test_6.png', mt, None, None) print('\\nFinished\\n') @@ -2129,6 +2168,7 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): 'FK_database=path_greens,', )) file.write(MisfitDefinitions) + file.write(OriginDefinitions) file.write(Grid_TestGraphics) file.write(Main_TestGraphics) diff --git a/tests/test_graphics.py b/tests/test_graphics.py index dfacd2c2..a6abb81e 100644 --- a/tests/test_graphics.py +++ b/tests/test_graphics.py @@ -18,13 +18,14 @@ # # Tests data, synthetics and beachball plotting utilities # - # Note that in the figures created by this script, the data and synthetics - # are not expected to fit epsecially well; currently, the only requirement - # is that the script runs without errors + + # + # The idea is for a test that runs very quickly, suitable for CI testing; + # eventually we may more detailed tests to tests/graphics # import matplotlib - matplotlib.use('Agg', warn=False, force=True) + matplotlib.use('Agg', force=True) import matplotlib path_greens= fullpath('data/tests/benchmark_cap/greens/scak') @@ -72,8 +73,21 @@ ) - mt = np.sqrt(1./3.)*np.array([1., 1., 1., 0., 0., 0.]) # explosion - mt *= 1.e16 + origin = Origin({ + 'time': '2009-04-07T20:12:55.000000Z', + 'latitude': 61.454200744628906, + 'longitude': -149.7427978515625, + 'depth_in_m': 33033.599853515625, + }) + + + from mtuq import MomentTensor + + mt = MomentTensor( + 1.e16 * np.sqrt(1./3.)*np.array([1., 1., 1., 0., 0., 0.])) # explosion + + mt_dict = { + 'rho':1.,'v':0.,'w':3/8*np.pi,'kappa':0.,'sigma':0.,'h':0.} wavelet = Trapezoid( magnitude=4.5) @@ -82,7 +96,7 @@ print('Reading data...\n') data = read(path_data, format='sac', event_id=event_id, - station_id_list=station_id_list, + #station_id_list=station_id_list, tags=['units:m', 'type:velocity']) @@ -108,23 +122,56 @@ # Generate figures # - print('Figure 1 of 3\n') + print('Plot data (1 of 6)\n') - plot_data_greens2('graphics_test_1.png', - data_bw, data_sw, greens_bw, greens_sw, - process_bw, process_sw, misfit_bw, misfit_sw, - stations, origin, mt, header=False) + from mtuq.graphics import plot_waveforms2 + from mtuq.util import Null + + plot_waveforms2('graphics_test_1.png', + data_bw, data_sw, Null(), Null(), + stations, origin, header=False) + + + print('Plot synthetics (2 of 6)\n') + + synthetics_bw = greens_bw.get_synthetics(mt, components=['Z','R']) + synthetics_sw = greens_sw.get_synthetics(mt, components=['Z','R','T']) - print('Figure 2 of 3\n') - plot_data_greens2('graphics_test_2.png', + plot_waveforms2('graphics_test_2.png', + synthetics_bw, synthetics_sw, Null(), Null(), + stations, origin, header=False) + + + print('Plot synthetics (3 of 6)\n') + + synthetics_bw = misfit_bw.collect_synthetics(data_bw, greens_bw, mt) + synthetics_sw = misfit_sw.collect_synthetics(data_sw, greens_sw, mt) + + + plot_waveforms2('graphics_test_3.png', + synthetics_bw, synthetics_sw, Null(), Null(), + stations, origin, header=False) + + + print('Plot data and synthetics without header (4 of 6)\n') + + plot_waveforms2('graphics_test_4.png', + data_bw, data_sw, synthetics_bw, synthetics_sw, + stations, origin, header=False) + + + print('Plot data and synthetics without header (5 of 6)\n') + + plot_data_greens2('graphics_test_5.png', data_bw, data_sw, greens_bw, greens_sw, process_bw, process_sw, misfit_bw, misfit_sw, - stations, origin, mt, header=False) + stations, origin, mt, mt_dict) + - print('Figure 3 of 3\n') + print('Plot explosion bechball (6 of 6)\n') - plot_beachball('graphics_test_3.png', + plot_beachball('graphics_test_6.png', mt, None, None) print('\nFinished\n') From fe7a8aaa96270ce089c4d6ab1e0b77c2d8f9723d Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 18 Nov 2024 16:08:29 -0700 Subject: [PATCH 6/9] Better handle empty datasets --- mtuq/graphics/waveforms.py | 156 +++++++++++++++++++++---------------- 1 file changed, 90 insertions(+), 66 deletions(-) diff --git a/mtuq/graphics/waveforms.py b/mtuq/graphics/waveforms.py index 84da822c..923c165a 100644 --- a/mtuq/graphics/waveforms.py +++ b/mtuq/graphics/waveforms.py @@ -12,7 +12,7 @@ from mtuq.graphics.annotations import trace_label_writer, station_label_writer,\ _getattr from mtuq.graphics.header import MomentTensorHeader, ForceHeader -from mtuq.util import warn +from mtuq.util import Null, warn from mtuq.util.signal import get_components, m_to_deg from obspy import Stream, Trace @@ -40,7 +40,7 @@ def plot_waveforms1( raise Exception # how many stations have at least one trace? - nrows = _count([data]) + nrows = _count(data, synthetics) # intialize figure fig, axes = _initialize1(nrows, header) @@ -106,7 +106,7 @@ def plot_waveforms2( """ # how many stations have at least one trace? - nrows = _count([data_bw, data_sw]) + nrows = _count(data_bw, data_sw, synthetics_bw, synthetics_bw) # intialize figure fig, axes = _initialize2(nrows, header) @@ -161,16 +161,16 @@ def plot_waveforms2( def plot_waveforms3( filename, data_bw, - data_rayleigh, + data_rayl, data_love, synthetics_bw, - synthetics_rayleigh, + synthetics_rayl, synthetics_love, stations, origin, header=None, total_misfit_bw=1., - total_misfit_rayleigh=1., + total_misfit_rayl=1., total_misfit_love=1., normalize='maximum_amplitude', trace_label_writer=trace_label_writer, @@ -182,7 +182,7 @@ def plot_waveforms3( """ # how many stations have at least one trace? - nrows = _count([data_bw, data_rayleigh, data_love]) + nrows = _count(data_bw, data_rayl, data_love) # intialize figure fig, axes = _initialize2(nrows, header) @@ -193,12 +193,12 @@ def plot_waveforms3( # optional normalization if normalize=='maximum_amplitude': factor_bw = _max(data_bw, synthetics_bw) - factor_rayleigh = _max(data_rayleigh, synthetics_rayleigh) + factor_rayl = _max(data_rayl, synthetics_rayl) factor_love = _max(data_love, synthetics_love) elif normalize=='median_amplitude': factor_bw = 2.*_median(data_bw, synthetics_bw) - factor_rayleigh = 2.*_median(data_rayleigh, synthetics_rayleigh) + factor_rayl = 2.*_median(data_rayl, synthetics_rayl) factor_love = 2.*_median(data_love, synthetics_love) else: @@ -216,7 +216,7 @@ def plot_waveforms3( for _i in range(len(stations)): # skip empty stations - if len(data_bw[_i]) == len(data_rayleigh[_i]) == len(data_love[_i]) == 0: + if len(data_bw[_i]) == len(data_rayl[_i]) == len(data_love[_i]) == 0: continue # add station labels @@ -232,8 +232,8 @@ def plot_waveforms3( # plot Rayleigh waves _plot_stream(axes[ir], [3,4], ['Z','R'], - data_rayleigh[_i], synthetics_rayeleigh[_i], - normalize, factor_rayleigh, trace_label_writer, total_misfit_rayleigh) + data_rayl[_i], synthetics_rayl[_i], + normalize, factor_rayl, trace_label_writer, total_misfit_rayl) # plot Love waves _plot_stream(axes[ir], [5], ['T'], @@ -346,16 +346,16 @@ def plot_data_greens2(filename, def plot_data_greens3( filename, data_bw, - data_rayleigh, + data_rayl, data_love, greens_bw, - greens_rayleigh, + greens_rayl, greens_love, process_data_bw, - process_data_rayleigh, + process_data_rayl, process_data_love, misfit_bw, - misfit_rayleigh, + misfit_rayl, misfit_love, stations, origin, @@ -373,8 +373,8 @@ def plot_data_greens3( synthetics_bw = misfit_bw.collect_synthetics( data_bw, greens_bw.select(origin), source) - synthetics_rayleigh = misfit_rayleigh.collect_synthetics( - data_rayleigh, greens_rayleigh.select(origin), source) + synthetics_rayl = misfit_rayl.collect_synthetics( + data_rayl, greens_rayl.select(origin), source) synthetics_love = misfit_love.collect_synthetics( data_love, greens_love.select(origin), source) @@ -383,8 +383,8 @@ def plot_data_greens3( total_misfit_bw = misfit_bw( data_bw, greens_bw.select(origin), source, optimization_level=0) - total_misfit_rayleigh = misfit_rayleigh( - data_rayleigh, greens_rayleigh.select(origin), source, optimization_level=0) + total_misfit_rayl = misfit_rayl( + data_rayl, greens_rayl.select(origin), source, optimization_level=0) total_misfit_love = misfit_love( data_love, greens_love.select(origin), source, optimization_level=0) @@ -398,17 +398,17 @@ def plot_data_greens3( header = _prepare_header( model, solver, source, source_dict, origin, - process_data_bw, process_data_rayleigh, misfit_bw, misfit_rayleigh, - total_misfit_bw, total_misfit_rayleigh, best_misfit_sw_supp=total_misfit_love, - misfit_sw_supp = misfit_love, data_bw=data_bw, data_sw=data_rayleigh, + process_data_bw, process_data_rayl, misfit_bw, misfit_rayl, + total_misfit_bw, total_misfit_rayl, best_misfit_sw_supp=total_misfit_love, + misfit_sw_supp = misfit_love, data_bw=data_bw, data_sw=data_rayl, data_sw_supp=data_love, process_sw_supp=process_data_love) plot_waveforms3(filename, - data_bw, data_rayleigh, data_love, - synthetics_bw, synthetics_rayleigh, synthetics_love, + data_bw, data_rayl, data_love, + synthetics_bw, synthetics_rayl, synthetics_love, stations, origin, total_misfit_bw=total_misfit_bw, - total_misfit_rayleigh=total_misfit_rayleigh, + total_misfit_rayl=total_misfit_rayl, total_misfit_love=total_misfit_love, header=header, **kwargs) @@ -450,12 +450,14 @@ def _initialize( margin_left=0.25, margin_right=0.25, header_height=1.5, station_labels=True, station_label_width=0.5): - if header: - height += header_height + if not header: + header_height = 0. if not station_labels: station_label_width = 0. + height += header_height + height += margin_bottom height += margin_top width += margin_left @@ -507,9 +509,13 @@ def _plot_stream( try: dat = stream_dat.select(component=component)[0] + except: + dat = None + + try: syn = stream_syn.select(component=component)[0] except: - continue + syn = None weight = _getattr(dat, 'weight', 1.) if not weight: @@ -527,7 +533,11 @@ def _plot_stream( ylim = [-amplitude_factor, +amplitude_factor] _plot_trace(axis, dat, syn) - axis.set_ylim(*ylim) + + try: + axis.set_ylim(*ylim) + except ValueError: + pass try: trace_label_writer(axis, dat, syn, total_misfit) @@ -538,20 +548,26 @@ def _plot_stream( def _plot_trace(axis, dat, syn, label=None): """ Plots data and synthetics time series on current axes """ - t1,t2,nt,dt = _time_stats(dat) + if dat is None and syn is None: + # nothing to plot + return + + if dat: + t,d = _time_series(dat) + + axis.plot(t, d, 'k', linewidth=1.5, + clip_on=True, zorder=10) - # which start and stop indices will correctly align synthetics? - start = _getattr(syn, 'idx_start', 0) - stop = _getattr(syn, 'idx_stop', len(syn.data)) + if syn: + # which start and stop indices will correctly align synthetics? + start = _getattr(syn, 'idx_start', 0) + stop = _getattr(syn, 'idx_stop', len(syn.data)) - t = np.linspace(0,t2-t1,nt,dt) - d = dat.data - s = syn.data + t,s = _time_series(syn) + + axis.plot(t[start:stop], s[start:stop], 'r', linewidth=1.25, + clip_on=True, zorder=10) - axis.plot(t, d, 'k', linewidth=1.5, - clip_on=True, zorder=10) - axis.plot(t, s[start:stop], 'r', linewidth=1.25, - clip_on=True, zorder=10) def _add_component_labels1(axes, body_wave_labels=True, surface_wave_labels=True): @@ -604,21 +620,22 @@ def _add_component_labels2(axes, body_wave_labels=True, surface_wave_labels=True # utility functions # -def _time_stats(trace): - # returns time scheme - return ( - float(trace.stats.starttime), - float(trace.stats.endtime), - trace.stats.npts, - trace.stats.delta, - ) +def _time_series(trace): + t1 = float(trace.stats.starttime) + t2 = float(trace.stats.endtime) + nt = trace.stats.npts + dt = trace.stats.delta + return np.linspace(0,t2-t1,nt,dt), trace.data + -def _count(datasets): +def _count(*datasets): # counts number of nonempty streams in dataset(s) count = 0 for streams in zip(*datasets): for stream in streams: + if stream is None: + continue if len(stream) > 0: count += 1 break @@ -629,29 +646,36 @@ def _isempty(dataset): if not dataset: return True else: - return bool(_count([dataset])==0) + return bool(_count(dataset)==0) -def _max(dat, syn): +def _max(*datasets): # returns maximum amplitude over traces, streams, or datasets - if type(dat)==type(syn)==Trace: - return max( - abs(dat.max()), - abs(syn.max())) + maxall = -np.Inf - elif type(dat)==type(syn)==Stream: - return max( - max(map(abs, dat.max())), - max(map(abs, syn.max()))) + for ds in datasets: + if type(ds) not in [Dataset, Stream, Trace, Null]: + print('Expected a type Dataset, Stream, Trace, Null ' + 'but received a type f{type(ds)}.\nSkipping...') - elif type(dat)==type(syn)==Dataset: - return max( - abs(dat.max()), - abs(syn.max())) + if type(ds)==Trace: + maxval = abs(ds.max()) + + elif type(ds)==Stream: + maxval = map(abs, ds.max()) + + elif type(ds)==dsaset: + maxval = abs(ds.max()) + + elif type(ds)==Null: + continue + + if maxval > maxall: + maxall = maxval + + return maxall - else: - raise TypeError def _median(*datasets): # returns median Linf amplitude over traces in dataset From 7ef57b44dc16b6df4264eff66301a1e463a616e7 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 18 Nov 2024 16:09:37 -0700 Subject: [PATCH 7/9] Try out CI graphics test --- .github/workflows/python-app.yaml | 3 +++ setup/code_generator.py | 42 +++++++++++++++---------------- tests/test_graphics.py | 4 +-- 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/.github/workflows/python-app.yaml b/.github/workflows/python-app.yaml index 5bcefdd4..b8991cd6 100644 --- a/.github/workflows/python-app.yaml +++ b/.github/workflows/python-app.yaml @@ -42,6 +42,9 @@ jobs: python tests/test_grid_search_mt_depth.py --no_figures python tests/test_greens_SPECFEM3D_SAC.py --no_figures python tests/test_time_shifts.py --no_figures + - name: Graphics + run: | + python tests/test_graphics.py # unfortunately, these Conda installation tests exceed the resource limits # for GitHub workflows diff --git a/setup/code_generator.py b/setup/code_generator.py index 7b98172b..546b5d94 100755 --- a/setup/code_generator.py +++ b/setup/code_generator.py @@ -977,10 +977,10 @@ plot_waveforms2('graphics_test_4.png', data_bw, data_sw, synthetics_bw, synthetics_sw, - stations, origin, header=False) + stations, origin, header=None) - print('Plot data and synthetics without header (5 of 6)\\n') + print('Plot data and synthetics with header (5 of 6)\\n') plot_data_greens2('graphics_test_5.png', data_bw, data_sw, greens_bw, greens_sw, @@ -1874,25 +1874,25 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): Main_GridSearch, r'origin', r'origins', - r'Reading Greens functions...\\\\n', + 'Reading Greens functions...\\\\n', ( - r'Reading Greens functions...\\\\n\\\\n'+ - r' Downloads can sometimes take as long as a few hours!\\\\n' + 'Reading Greens functions...\\\\n\\\\n'+ + ' Downloads can sometimes take as long as a few hours!\\\\n' ), - r'download_greens_tensors\(stations, origin, model\)', - r'download_greens_tensors(stations, origin, model, verbose=True)', + 'download_greens_tensors\(stations, origin, model\)', + 'download_greens_tensors(stations, origin, model, verbose=True)', )) file.write( replace( WrapUp_GridSearch_DoubleCoupleMagnitudeDepth, - r'DC\+Z', - r'DC+XY', - r'misfit_depth', - r'misfit_latlon', - r"title=event_id", - r"title=event_id, colorbar_label='L2 misfit'", - r'show_magnitudes=True, ', - r'', + 'DC\+Z', + 'DC+XY', + 'misfit_depth', + 'misfit_latlon', + "title=event_id", + "title=event_id, colorbar_label='L2 misfit'", + 'show_magnitudes=True, ', + '', )) @@ -2044,9 +2044,9 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): file.write( replace( Main1_SerialGridSearch_DoubleCouple, - r'greens = download_greens_tensors\(stations, origin, model\)', - r'db = open_db(path_greens, format=\'FK\', model=model)\n ' - +r'greens = db.get_greens_tensors(stations, origin)', + 'greens = download_greens_tensors\(stations, origin, model\)', + 'db = open_db(path_greens, format=\'FK\', model=model)\n ' + +'greens = db.get_greens_tensors(stations, origin)', )) file.write( replace( @@ -2110,9 +2110,9 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): file.write( replace( Main1_SerialGridSearch_DoubleCouple, - r'greens = download_greens_tensors\(stations, origin, model\)', - r'db = open_db(path_greens, format=\'FK\', model=model)\n ' - +r'greens = db.get_greens_tensors(stations, origin)', + 'greens = download_greens_tensors\(stations, origin, model\)', + 'db = open_db(path_greens, format=\'FK\', model=model)\n ' + +'greens = db.get_greens_tensors(stations, origin)', )) file.write(Main_TestMisfit) diff --git a/tests/test_graphics.py b/tests/test_graphics.py index a6abb81e..ec0ef508 100644 --- a/tests/test_graphics.py +++ b/tests/test_graphics.py @@ -158,10 +158,10 @@ plot_waveforms2('graphics_test_4.png', data_bw, data_sw, synthetics_bw, synthetics_sw, - stations, origin, header=False) + stations, origin, header=None) - print('Plot data and synthetics without header (5 of 6)\n') + print('Plot data and synthetics with header (5 of 6)\n') plot_data_greens2('graphics_test_5.png', data_bw, data_sw, greens_bw, greens_sw, From 3ecf41eec0403dbd09d36fb8a91cc10186882e31 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Wed, 20 Nov 2024 15:59:50 -0700 Subject: [PATCH 8/9] Single station plotting fix --- mtuq/graphics/waveforms.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mtuq/graphics/waveforms.py b/mtuq/graphics/waveforms.py index 923c165a..d1dd1174 100644 --- a/mtuq/graphics/waveforms.py +++ b/mtuq/graphics/waveforms.py @@ -478,8 +478,6 @@ def _initialize( hspace=0., ) - _hide_axes(axes) - if header: header.write( header_height, width, @@ -489,6 +487,8 @@ def _initialize( if nrows==1: axes = [axes] + _hide_axes(axes) + return fig, axes From 0a3b8768dc1e0de188d542e2830719e749507a6b Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 16 Jan 2025 20:56:16 -0700 Subject: [PATCH 9/9] Updated README links --- README.md | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 350a3ca9..3bf3846c 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # mtuq -[![Build Status](https://github.com/uafgeotools/mtuq/actions/workflows/python-app.yaml/badge.svg)](https://github.com/rmodrak/mtuq/blob/master/.github/workflows/python-app.yaml) +[![Build Status](https://github.com/mtuqorg/mtuq/actions/workflows/python-app.yaml/badge.svg)](https://github.com/rmodrak/mtuq/blob/master/.github/workflows/python-app.yaml) [![SCOPED](https://img.shields.io/endpoint?url=https://runkit.io/wangyinz/scoped/branches/master/MTUQ)](https://github.com/SeisSCOPED/container/pkgs/container/MTUQ) MTUQ provides *m*oment *t*ensor estimates and *u*ncertainty *q*uantification from broadband seismic data. @@ -8,33 +8,33 @@ MTUQ provides *m*oment *t*ensor estimates and *u*ncertainty *q*uantification fro ## Getting started -[Installation](https://uafgeotools.github.io/mtuq/install/index.html) +[Installation](https://mtuqorg.github.io/mtuq/install/index.html) -[Quick start](https://uafgeotools.github.io/mtuq/quick_start.html) +[Quick start](https://mtuqorg.github.io/mtuq/quick_start.html) ## Documentation -[Acquiring seismic data](https://uafgeotools.github.io/mtuq/user_guide/02.html) +[Acquiring seismic data](https://mtuqorg.github.io/mtuq/user_guide/02.html) -[Acquiring Green's functions](https://uafgeotools.github.io/mtuq/user_guide/03.html) +[Acquiring Green's functions](https://mtuqorg.github.io/mtuq/user_guide/03.html) -[Data processing](https://uafgeotools.github.io/mtuq/user_guide/04.html) +[Data processing](https://mtuqorg.github.io/mtuq/user_guide/04.html) -[Visualization galleries](https://uafgeotools.github.io/mtuq/user_guide/05.html) +[Visualization galleries](https://mtuqorg.github.io/mtuq/user_guide/05.html) -[Library reference](https://uafgeotools.github.io/mtuq/library/index.html) +[Library reference](https://mtuqorg.github.io/mtuq/library/index.html) ## Highlights -Common use cases include [double couple moment tensor](https://github.com/uafgeotools/mtuq/blob/master/examples/SerialGridSearch.DoubleCouple.py), [full moment tensor](https://github.com/uafgeotools/mtuq/blob/master/examples/GridSearch.FullMomentTensor.py), [depth](https://github.com/rmodrak/mtuq/blob/master/examples/GridSearch.DoubleCouple%2BMagnitude%2BDepth.py) and [hypocenter](https://github.com/rmodrak/mtuq/blob/master/examples/GridSearch.DoubleCouple%2BMagnitude%2BHypocenter.py) uncertainty analysis. Applications involving composite sources, force sources, constrained moment tensor sources, source-time functions, and other source parameters are also possible. +Common use cases include [double couple moment tensor](https://github.com/mtuqorg/mtuq/blob/master/examples/SerialGridSearch.DoubleCouple.py), [full moment tensor](https://github.com/mtuqorg/mtuq/blob/master/examples/GridSearch.FullMomentTensor.py), [depth](https://github.com/rmodrak/mtuq/blob/master/examples/GridSearch.DoubleCouple%2BMagnitude%2BDepth.py) and [hypocenter](https://github.com/rmodrak/mtuq/blob/master/examples/GridSearch.DoubleCouple%2BMagnitude%2BHypocenter.py) uncertainty analysis. Applications involving composite sources, force sources, constrained moment tensor sources, source-time functions, and other source parameters are also possible. ### Solver interfaces -[I/O functions](https://uafgeotools.github.io/mtuq/library/index.html#data-i-o) +[I/O functions](https://mtuqorg.github.io/mtuq/library/index.html#data-i-o) are included for reading AxiSEM, SPECFEM3D, and FK Green's functions as well as downloading Green's functions from remote [syngine](http://ds.iris.edu/ds/products/syngine/) databases. @@ -42,22 +42,22 @@ downloading Green's functions from remote [syngine](http://ds.iris.edu/ds/produc ### Misfit evaluation -Waveform difference and cross-correlation time-shift [misfit evaluation](https://uafgeotools.github.io/mtuq/library/index.html#data-processing-and-inversion) +Waveform difference and cross-correlation time-shift [misfit evaluation](https://mtuqorg.github.io/mtuq/library/index.html#data-processing-and-inversion) on body-wave and surface-wave windows is implemented in C-accelerated Python. -These misfit functions can be used with [mtuq.grid_search](https://uafgeotools.github.io/mtuq/library/generated/mtuq.grid_search.grid_search.html), which automatically partitions the grid over multiple MPI processes if invoked from an MPI environment. For efficient and unbiased uncertainty quantification, [uniform grids](https://uafgeotools.github.io/mtuq/library/index.html#moment-tensor-and-force-grids) can be used for the grid search, drawing from [Tape2015](https://academic.oup.com/gji/article/202/3/2074/613765). +These misfit functions can be used with [mtuq.grid_search](https://mtuqorg.github.io/mtuq/library/generated/mtuq.grid_search.grid_search.html), which automatically partitions the grid over multiple MPI processes if invoked from an MPI environment. For efficient and unbiased uncertainty quantification, [uniform grids](https://mtuqorg.github.io/mtuq/library/index.html#moment-tensor-and-force-grids) can be used for the grid search, drawing from [Tape2015](https://academic.oup.com/gji/article/202/3/2074/613765). Alternatively, MTUQ misfit functions can be used as a starting point for Bayesian uncertainty quantification using pymc or other MCMC libraries. ### Visualization -[Visualization utilities](https://uafgeotools.github.io/mtuq/user_guide/05/gallery_mt.html) are included for both the [eigenvalue lune](https://onlinelibrary.wiley.com/doi/10.1111/j.1365-246X.2012.05491.x) and [v,w rectangle](https://academic.oup.com/gji/article/202/3/2074/613765), with matplotlib and Generic Mapping Tools graphics backends. +[Visualization utilities](https://mtuqorg.github.io/mtuq/user_guide/05/gallery_mt.html) are included for both the [eigenvalue lune](https://onlinelibrary.wiley.com/doi/10.1111/j.1365-246X.2012.05491.x) and [v,w rectangle](https://academic.oup.com/gji/article/202/3/2074/613765), with matplotlib and Generic Mapping Tools graphics backends. ### Testing -The package has been tested against [legacy Perl/C codes](https://github.com/uafgeotools/mtuq/blob/master/tests/benchmark_cap_vs_mtuq.py) as well as [published studies](https://github.com/rmodrak/mtbench). +The package has been tested against [legacy Perl/C codes](https://github.com/mtuqorg/mtuq/blob/master/tests/benchmark_cap_vs_mtuq.py) as well as [published studies](https://github.com/rmodrak/mtbench).