From bea12abe84d61cf8e2004706122c6dfdf18fef5b Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 5 Apr 2024 15:43:33 +1100 Subject: [PATCH 1/5] Update waveforms.py Implementation of wavefield normalization from PR #252. Adding new feature (median-based normalization), and fixing station_amplitude. --- mtuq/graphics/waveforms.py | 114 +++++++++++++++++++++++++++++-------- 1 file changed, 91 insertions(+), 23 deletions(-) diff --git a/mtuq/graphics/waveforms.py b/mtuq/graphics/waveforms.py index 27d22f26a..a85d5cdb6 100644 --- a/mtuq/graphics/waveforms.py +++ b/mtuq/graphics/waveforms.py @@ -62,6 +62,17 @@ def plot_waveforms1(filename, 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 + else: + raise ValueError("Invalid normalization method specified.") + # # loop over stations # @@ -100,7 +111,7 @@ def plot_waveforms1(filename, continue _plot_ZRT(axes[ir], 1, dat, syn, component, - normalize, trace_label_writer, max_amplitude, total_misfit) + normalize, trace_label_writer, max_amplitudes[_i], total_misfit) ir += 1 @@ -153,6 +164,23 @@ def plot_waveforms2(filename, 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))]) + else: + raise ValueError("Invalid normalization method specified.") + # # loop over stations # @@ -191,7 +219,7 @@ def plot_waveforms2(filename, continue _plot_ZR(axes[ir], 1, dat, syn, component, - normalize, trace_label_writer, max_amplitude_bw, total_misfit_bw) + normalize, trace_label_writer, max_amplitudes_bw[_i], total_misfit_bw) # @@ -216,7 +244,7 @@ def plot_waveforms2(filename, continue _plot_ZRT(axes[ir], 3, dat, syn, component, - normalize, trace_label_writer, max_amplitude_sw, total_misfit_sw) + normalize, trace_label_writer, max_amplitudes_sw[_i], total_misfit_sw) ir += 1 @@ -373,7 +401,7 @@ def _initialize(nrows=None, ncolumns=None, column_width_ratios=None, def _plot_ZRT(axes, ic, dat, syn, component, normalize='maximum_amplitude', trace_label_writer=None, - max_amplitude=1., total_misfit=1.): + normalization_amplitude=1., total_misfit=1.): # plot traces if component=='Z': @@ -387,17 +415,13 @@ def _plot_ZRT(axes, ic, dat, syn, component, _plot(axis, dat, syn) - # normalize amplitude + # 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': - max_stream = _max(stream_dat, stream_syn) - ylim = [-1.5*max_stream, +1.5*max_stream] - axis.set_ylim(*ylim) - elif normalize=='maximum_amplitude': - ylim = [-0.75*max_amplitude, +0.75*max_amplitude] + 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: @@ -406,7 +430,7 @@ def _plot_ZRT(axes, ic, dat, syn, component, def _plot_ZR(axes, ic, dat, syn, component, normalize='maximum_amplitude', trace_label_writer=None, - max_amplitude=1., total_misfit=1.): + normalization_amplitude=1., total_misfit=1.): # plot traces if component=='Z': @@ -418,20 +442,15 @@ def _plot_ZR(axes, ic, dat, syn, component, _plot(axis, dat, syn) - # normalize amplitude + # 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': - max_stream = _max(stream_dat, stream_syn) - ylim = [-1.5*max_stream, +1.5*max_stream] - axis.set_ylim(*ylim) - elif normalize=='maximum_amplitude': - ylim = [-0.75*max_amplitude, +0.75*max_amplitude] + 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) @@ -450,9 +469,9 @@ def _plot(axis, dat, syn, label=None): s = syn.data axis.plot(t, d, 'k', linewidth=1.5, - clip_on=False, zorder=10) + clip_on=True, zorder=10) axis.plot(t, s[start:stop], 'r', linewidth=1.25, - clip_on=False, zorder=10) + clip_on=True, zorder=10) def _add_component_labels1(axes, body_wave_labels=True, surface_wave_labels=True): @@ -534,6 +553,19 @@ 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: + 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()), @@ -552,6 +584,43 @@ 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. + + 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 def _hide_axes(axes): @@ -597,4 +666,3 @@ def _get_tag(tags, pattern): else: return None - From 40d1bfe19cbb3053610a22ec87782d3f612cd930 Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 5 Apr 2024 15:44:23 +1100 Subject: [PATCH 2/5] Update annotations.py Improved degree's distance so that it uses latitude/longitude instead of m_to_deg. --- mtuq/graphics/annotations.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mtuq/graphics/annotations.py b/mtuq/graphics/annotations.py index 9e8221324..aff599e05 100644 --- a/mtuq/graphics/annotations.py +++ b/mtuq/graphics/annotations.py @@ -3,6 +3,8 @@ from matplotlib import pyplot from obspy.geodetics import gps2dist_azimuth +# location to degree distance with obspy +from obspy.geodetics import locations2degrees def station_label_writer(ax, station, origin, units='km'): @@ -31,7 +33,8 @@ def station_label_writer(ax, station, origin, units='km'): label = '%d km' % round(distance_in_m/1000.) elif units=='deg': - label = '%d%s' % (round(m_to_deg(distance_in_m)), u'\N{DEGREE SIGN}') + label = '%d%s' % (round(locations2degrees(origin.latitude, origin.longitude, + station.latitude, station.longitude)), u'\N{DEGREE SIGN}') pyplot.text(0.2,0.35, label, fontsize=11, transform=ax.transAxes) @@ -89,4 +92,3 @@ def _getattr(trace, name, *args): else: raise TypeError("Wrong number of arguments") - From 36e859d167d7c2e5ecc7765200cd26a593d7e687 Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 5 Apr 2024 15:44:49 +1100 Subject: [PATCH 3/5] Update header.py Fixed typo. --- mtuq/graphics/header.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mtuq/graphics/header.py b/mtuq/graphics/header.py index 2f3398fa8..616b36544 100644 --- a/mtuq/graphics/header.py +++ b/mtuq/graphics/header.py @@ -114,7 +114,7 @@ def parse_data_processing(self): if not self.process_bw: pass if not self.process_sw: - raise Excpetion() + raise Exception() if self.process_sw.freq_max > 1.: units = 'Hz' From a363b70a4d28418d1aba0366740aa6b79cce144a Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 5 Apr 2024 18:29:42 +1100 Subject: [PATCH 4/5] Update unpack.bash Removed extra character in one of the file's name. Would prevent unpacking otherwise. --- data/examples/unpack.bash | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data/examples/unpack.bash b/data/examples/unpack.bash index da3897878..6c798a50d 100755 --- a/data/examples/unpack.bash +++ b/data/examples/unpack.bash @@ -9,7 +9,7 @@ cd $(dirname ${BASH_SOURCE[0]}) wd=$PWD for filename in \ - 20090407201255351.tgz 20210809074550.tgz 20SPECFEM3D_SGT.tgz SPECFEM3D_SAC.tgz; + 20090407201255351.tgz 20210809074550.tgz SPECFEM3D_SGT.tgz SPECFEM3D_SAC.tgz; do cd $wd cd $(dirname $filename) From da8ea0e691f4e7b6cf0a5dabc2040919802a9032 Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 5 Apr 2024 23:14:17 +1100 Subject: [PATCH 5/5] Update setup.py Scipy version fix to avoid conflict with obspy. (scipy.signal.hann is deprecated since scipy 1.13.0 and called in obspy). --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 6e4543725..46e789d58 100644 --- a/setup.py +++ b/setup.py @@ -100,7 +100,7 @@ def run_tests(self): # (consider using a conda based installation instead) install_requires=[ "numpy", - "scipy", + "scipy<1.13.0", "pandas", "xarray", "netCDF4",