From e91adbc791b005beb50e0422b6585924aef6cb09 Mon Sep 17 00:00:00 2001 From: Andrew Medlin Date: Fri, 28 Feb 2020 10:18:18 +1100 Subject: [PATCH 1/2] Changed to high accuracy interpolater to prevent quantization noise. Cut input dataset to OA region to prevent running out of memory on VDI instances. Changed gravity interpolated resampling along transect to use uniform spatial sampling rate. --- seismic/receiver_fn/plot_ccp_batch.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/seismic/receiver_fn/plot_ccp_batch.py b/seismic/receiver_fn/plot_ccp_batch.py index c22a4d7c..2287d57d 100644 --- a/seismic/receiver_fn/plot_ccp_batch.py +++ b/seismic/receiver_fn/plot_ccp_batch.py @@ -8,8 +8,10 @@ import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec +# import matplotlib.tri as tri import pandas as pd from scipy import interpolate +from obspy.geodetics.base import gps2dist_azimuth import seismic.receiver_fn.rf_util as rf_util from seismic.ASDFdatabase import FederatedASDFDataSet @@ -215,7 +217,8 @@ def gravity_subplot(hf, metadata, grav_map): pos_cb.y0 += 0.2 hf.axes[1].set_position(pos_cb) - # Add gravity plot + # Add gravity plot. Ignores curvature of earth in transect projection to distance + # from starting point, just uses Euclidean geometry. start = metadata['transect_start'] end = metadata['transect_end'] dirn = metadata['transect_dirn'] @@ -227,7 +230,10 @@ def gravity_subplot(hf, metadata, grav_map): ax_grav.set_position(pos_grav) plt.sca(ax_grav) plt.title("Gravity survey", fontsize=8, y=0.80) - grav_pos = np.linspace(start, end, 1000) + RESAMPLES_PER_KM = 2.0 + transect_len_km = gps2dist_azimuth(*start, *end)[0]/1000.0 + num_samples = int(np.ceil(transect_len_km*RESAMPLES_PER_KM)) + grav_pos = np.linspace(start, end, num_samples) grav_vals = grav_map(grav_pos) grav_dist = np.dot((grav_pos - start), dirn)*KM_PER_DEG plt.plot(grav_dist, grav_vals) @@ -277,12 +283,16 @@ def main(transect_file, output_folder, rf_file, waveform_database, stack_scale, assert len(channel) == 1, "Batch stack only on one channel at a time" # Custom plot modifiers. Leave commented for now until refactoring in ticket PST-479 - # print("Loading gravity data...") - # grav = np.load('post_analysis/GravityGrid.xyz.npy') - # print("Creating interpolator...") - # grav_map = interpolate.NearestNDInterpolator(grav[:, 0:2], grav[:, 2]) + print("Loading gravity data...") + grav = np.load('post_analysis/GravityGrid.xyz.npy') + # Cut gravity map to region of interest + OA_region_mask = (grav[:, 1] >= -23) & (grav[:, 1] <= -17) & (grav[:, 0] >= 132) & (grav[:, 0] <= 142) + grav = grav[OA_region_mask, :] + print("Creating interpolator...") + grav_map = interpolate.CloughTocher2DInterpolator(grav[:, 0:2], grav[:, 2]) # annotators = [moho_annotator, lambda hf, md: gravity_subplot(hf, md, grav_map)] - annotators = None + annotators = [lambda hf, md: gravity_subplot(hf, md, grav_map)] + # annotators = None print("Producing plot...") run_batch(transect_file, rf_file, waveform_database, stack_scale=stack_scale, width=width, spacing=spacing, max_depth=depth, channel=channel, output_folder=output_folder, colormap=colormap, From 241f998d7019080f622b719bec841d1bed16320e Mon Sep 17 00:00:00 2001 From: Andrew Medlin Date: Fri, 28 Feb 2020 13:20:09 +1100 Subject: [PATCH 2/2] Added low pass filtering to gravity signal. --- seismic/receiver_fn/plot_ccp_batch.py | 42 ++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/seismic/receiver_fn/plot_ccp_batch.py b/seismic/receiver_fn/plot_ccp_batch.py index 2287d57d..28bc5837 100644 --- a/seismic/receiver_fn/plot_ccp_batch.py +++ b/seismic/receiver_fn/plot_ccp_batch.py @@ -8,9 +8,10 @@ import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec -# import matplotlib.tri as tri import pandas as pd from scipy import interpolate +from scipy.signal.windows import tukey +from scipy.signal import firwin, tf2sos, sosfiltfilt from obspy.geodetics.base import gps2dist_azimuth import seismic.receiver_fn.rf_util as rf_util @@ -208,6 +209,31 @@ def gravity_subplot(hf, metadata, grav_map): return # end if + def tukey_lowpass_filter(x, filter_len_km): + filter_len = int(np.ceil(filter_len_km * RESAMPLES_PER_KM)) + if (filter_len % 2) == 0: + filter_len += 1 + # end if + kern = firwin(filter_len, 1.0 / filter_len_km, window=('tukey', 0.3), fs=RESAMPLES_PER_KM) + sos_kern = tf2sos(kern, 1) + sig_filt = sosfiltfilt(sos_kern, x) + return sig_filt + # end func + + def tukey_rolling_mean(x, filter_len_km): + filter_len = int(np.ceil(filter_len_km * RESAMPLES_PER_KM)) + if (filter_len % 2) == 0: + filter_len += 1 + # end if + kern = tukey(filter_len, 0.3) + kern = kern/np.sum(kern) + sos_kern = tf2sos(kern, 1) + sig_filt = sosfiltfilt(sos_kern, x) + return sig_filt + # end func + + RESAMPLES_PER_KM = 2.0 + # Move the bottom of the main axes bounding box up to make space to gravity plot beneath pos = hf.axes[0].get_position() pos.y0 += 0.2 @@ -230,14 +256,22 @@ def gravity_subplot(hf, metadata, grav_map): ax_grav.set_position(pos_grav) plt.sca(ax_grav) plt.title("Gravity survey", fontsize=8, y=0.80) - RESAMPLES_PER_KM = 2.0 - transect_len_km = gps2dist_azimuth(*start, *end)[0]/1000.0 + transect_len_km = gps2dist_azimuth(start[1], start[0], end[1], end[0])[0]/1000.0 num_samples = int(np.ceil(transect_len_km*RESAMPLES_PER_KM)) grav_pos = np.linspace(start, end, num_samples) grav_vals = grav_map(grav_pos) grav_dist = np.dot((grav_pos - start), dirn)*KM_PER_DEG - plt.plot(grav_dist, grav_vals) + grav_vals_lp_filt = tukey_lowpass_filter(grav_vals, 15.0) + # grav_vals_tk_filt = tukey_rolling_mean(grav_vals, 15.0) + grav_vals_lp_filt_long = tukey_lowpass_filter(grav_vals, 70.0) + # grav_vals_tk_filt_long = tukey_rolling_mean(grav_vals, 70.0) + plt.plot(grav_dist, grav_vals, label='Original Bouger') + plt.plot(grav_dist, grav_vals_lp_filt, '--', alpha=0.6, label='Designed Tukey filter (15 km)') + # plt.plot(grav_dist, grav_vals_tk_filt, '-.', alpha=0.6, label='Tukey rolling window (15 km)') + plt.plot(grav_dist, grav_vals_lp_filt_long, '--', alpha=0.6, label='Designed Tukey filter (70 km)') + # plt.plot(grav_dist, grav_vals_tk_filt_long, '-.', alpha=0.6, label='Tukey rolling window (70 km)') plt.grid("#80808080", linestyle=':') + plt.legend(loc=3, fontsize=8) tickstep_x = 50.0 xlim = hf.axes[0].get_xlim() plt.xticks(np.arange(0.0, xlim[1], tickstep_x), fontsize=12)