diff --git a/.gitignore b/.gitignore index 027ebd1..81310e5 100644 --- a/.gitignore +++ b/.gitignore @@ -130,3 +130,6 @@ dmypy.json # And bespoke ignores: ./**/*.old + +# Mac OS specific +*.DS_store diff --git a/README.md b/README.md index 02e1bde..8f37424 100644 --- a/README.md +++ b/README.md @@ -4,3 +4,5 @@ SeisSeeker is a package for performing beamforming for earthquake detection. It A decription of how the package works can be found here: Thomas S. Hudson, Alex M. Brisbourne, Sofia-Katerina Kufner, J-Michael Kendall, and Andy M. Smith. (in review). "Array processing in cryoseismology". Submitted to: The Cryosphere. + +This version of SeisSeeker has been modified for application to monitoring of North Sea Microseismicity diff --git a/SeisSeeker/processing/detection.py b/SeisSeeker/processing/detection.py index f88564c..b754657 100755 --- a/SeisSeeker/processing/detection.py +++ b/SeisSeeker/processing/detection.py @@ -1,141 +1,202 @@ #!/usr/bin/python -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Script Description: # Script to perform earthquake detection using array processing methods. # Created by Tom Hudson, 10th August 2022 -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Import neccessary modules: import pandas as pd import numpy as np -import matplotlib -import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D -import os, sys +from pathlib import Path +import matplotlib.pyplot as plt +import os import obspy -from scipy.signal import find_peaks -from numba import jit, objmode, prange, set_num_threads -import gc -# import multiprocessing as mp -import time -import glob +import datetime +from scipy.signal import find_peaks, hilbert +from numba import jit, objmode, prange +import gc +import logging + +import time +import glob import pickle -from SeisSeeker.processing import lookup_table_manager, location +from SeisSeeker.processing import lookup_table_manager, location +logger = logging.getLogger(__name__) +# GLobal constants: +MAD_SCALE = 1.4826 # Scale factor to convert MAD to std. dev. -#----------------------------------------------- Define main functions ----------------------------------------------- + +# ---------- Define main functions ------------ class CustomError(Exception): pass -def flatten_list(l): - return [item for sublist in l for item in sublist] +def flatten_list(list_to_flatten): + return [item for sublist in list_to_flatten for item in sublist] -def xy_to_rtheta(x,y): +def xy_to_rtheta(x, y): """x,y to r,theta, where x,y in East and North directions. Theta is in degrees from N.""" r = np.sqrt(x**2 + y**2) - theta = np.rad2deg(np.arctan2(x,y)) + theta = np.rad2deg(np.arctan2(x, y)) try: - theta[theta<0] = theta[theta<0] + 360 + theta[theta < 0] = theta[theta < 0] + 360 except TypeError: theta = theta + 360 - return r, theta + return r, theta + + +@jit(nopython=True, parallel=True) +def _fast_freq_domain_array_proc( + data, + min_sl, + max_sl, + n_sl, + min_baz, + max_baz, + n_baz, + fs, + target_freqs, + xx, + yy, + n_stations, + n_t_samp, + remove_autocorr, +): + """ + Performs array processing using method inspired by Bowden et al. (2021). + Performs array processing in polar coordinates. + Designed to be wrapped using Numba to improve performance. + Parameters: + ---------- + data : np.ndarray + 3D numpy array of data to process. Shape must be + (n_windows, n_stations, n_t_samp). + min_sl : float + Minimum slowness to analyse for, in s/km. + max_sl : float + Maximum slowness to analyse for, in s/km. + n_sl : int + Number of slowness values to analyse between min_sl and max_sl. + min_baz : float + Minimum back-azimuth, in degrees. + max_baz : float + Maximum back-azimuth, in degrees. + n_baz : int + Number of back-azimuth values to analyse between min_baz and max_baz. + fs : float + Sampling frequency of the data, in Hz. + target_freqs : list + List of target frequencies to analyse, in Hz. + xx : np.ndarray + 1D numpy array of x-coordinates of station locations, in km. + yy : np.ndarray + 1D numpy array of y-coordinates of station locations, in km. + n_stations : int + Number of stations in the array. + n_t_samp : int + Number of time samples in the data. + remove_autocorr : bool + Whether to remove autocorrelations from the data. -@jit(nopython=True, parallel=True)#, nogil=True) -def _fast_freq_domain_array_proc(data, max_sl, fs, target_freqs, xx, yy, n_stations, n_t_samp, remove_autocorr): - """Function to perform array processing fast due to being designed to - be wrapped using Numba. Function inspired by Bowden et al. (2021). - Performs array processing in polar coordinates. Returns: - Pfreq_all + ---------- + Pfreq_all : np.ndarray + 4D numpy array of processed data. Shape will be + (n_windows, len(target_freqs), n_sl, n_baz). """ # Define grid of slownesses: - # number of pixes in x and y - # (Determines number of phase shifts to perform) - n_ur = 26 #51 #101 - n_utheta = 120 #51 #101 - ur = np.linspace(0,max_sl,n_ur) - utheta = np.linspace(0,360-(360/n_utheta),n_utheta) + ur = np.linspace(min_sl, max_sl, n_sl) + utheta = np.linspace(min_baz, max_baz, n_baz) utheta_rad = np.deg2rad(utheta) - dur=ur[1]-ur[0] - dutheta=utheta[1]-utheta[0] - - # Compute time-shifts once: - # (so that don't have to do it for every frequency) - tlib = np.zeros((n_stations,n_ur,n_utheta), dtype=np.complex128) - for ir in range(0,n_ur): - for itheta in range(0,n_utheta): - # tlib[:,ix,iy] = xx*ux[ix] + yy*uy[iy] # (distance x slowness = distance / velocity = time) - tlib[:,ir,itheta] = xx*ur[ir]*np.sin((utheta_rad[itheta])) + yy*ur[ir]*np.cos((utheta_rad[itheta])) # (distance x slowness = distance / velocity = time) - # Since receivers are relative to the array centre, can shift all receivers back to that centre. - - # Create data stores: - Pfreq_all = np.zeros((data.shape[0],len(target_freqs),n_ur,n_utheta), dtype=np.complex128) # Explicitly create Pxx_all, as otherwise prange won't work correctly. - - # Then loop over windows: + + # # Compute time-shifts once: + tlib = np.zeros((n_stations, n_sl, n_baz), dtype=np.complex128) + # r, theta as this is polar coord system: + for ir in range(n_sl): + for itheta in range(n_baz): + tlib[:, ir, itheta] = xx * ur[ir] * np.sin(utheta_rad[itheta]) + yy * ur[ + ir + ] * np.cos(utheta_rad[itheta]) + + # Vectorized computation of time-shifts for all stations, slowness, and back-azimuth + # xx and yy are (n_stations,) arrays, ur is (n_sl,), utheta_rad is (n_baz,) + # We want tlib shape (n_stations, n_sl, n_baz) + + # xx_ = xx[:, np.newaxis, np.newaxis] # shape (n_stations, 1, 1) + # yy_ = yy[:, np.newaxis, np.newaxis] # shape (n_stations, 1, 1) + # ur_ = ur[np.newaxis, :, np.newaxis] # shape (1, n_sl, 1) + # utheta_rad_ = utheta_rad[np.newaxis, np.newaxis, :] # shape (1, 1, n_baz) + # tlib_tmp = xx_ * ur_ * np.sin(utheta_rad_) + yy_ * ur_ * np.cos(utheta_rad_) + # recast array as complex128 + # tlib = tlib_tmp.astype(np.complex128) + Pfreq_all = np.zeros( + (data.shape[0], len(target_freqs), n_sl, n_baz), dtype=np.complex128 + ) for win_idx in prange(data.shape[0]): - # Calculate spectra: - # Construct data structure: - nfft = (2.0**np.ceil(np.log2(n_t_samp))) - nfft = np.array(nfft, dtype=np.int64) - Pxx_all = np.zeros((np.int((nfft/2)+1), n_stations), dtype=np.complex128) # Power spectra - dt = 1. / fs - df = 1.0/(2.0*nfft*dt) - xf = np.linspace(0.0, 1.0/(2.0*dt), np.int((nfft/2)+1)) + nfft = int(2 ** np.ceil(np.log2(n_t_samp))) + Pxx_all = np.zeros((int(nfft / 2) + 1, n_stations), dtype=np.complex128) + dt = 1.0 / fs + xf = np.linspace(0.0, 1.0 / (2.0 * dt), int(nfft / 2) + 1) # Calculate power spectra for all stations: for sta_idx in range(n_stations): - # Calculate spectra for current station: - ###Pxx_all[:,sta_idx] = np.fft.rfft(data[win_idx,sta_idx,:], n=nfft) # (Use real fft, as input data is real) # DOESN'T WORK WITH NUMBA! - with objmode(Pxx_curr='complex128[:]'): - Pxx_curr = np.fft.rfft(data[win_idx,sta_idx,:], n=nfft) - Pxx_all[:,sta_idx] = Pxx_curr - - # Loop over all freqs, performing phase shifts: - Pfreq=np.zeros((len(target_freqs),n_ur,n_utheta),dtype=np.complex128) - counter_grid = 0 + with objmode(Pxx_curr="complex128[:]"): + Pxx_curr = np.fft.rfft(data[win_idx, sta_idx, :], n=nfft) + Pxx_all[:, sta_idx] = Pxx_curr + + Pfreq = np.zeros((len(target_freqs), n_sl, n_baz), dtype=np.complex128) for ii in range(len(target_freqs)): - # Find closest current freq.: target_f = target_freqs[ii] curr_f_idx = (np.abs(xf - target_f)).argmin() - # Construct a matrix of each station-station correlation before any phase shifts - Rxx=np.zeros((n_stations,n_stations),dtype=np.complex128) - for i1 in range(0,n_stations): - for i2 in range(0,n_stations): - # Remove autocorrelations: + Rxx = np.zeros((n_stations, n_stations), dtype=np.complex128) + for i1 in range(n_stations): + for i2 in range(n_stations): if remove_autocorr: - if not i1 == i2: - Rxx[i1,i2] = np.conj(Pxx_all[curr_f_idx,i1]) * Pxx_all[curr_f_idx,i2] + if i1 != i2: + Rxx[i1, i2] = ( + np.conj(Pxx_all[curr_f_idx, i1]) + * Pxx_all[curr_f_idx, i2] + ) else: - Rxx[i1,i2] = 0 - - # And loop over phase shifts, calculating cross-correlation power: - for ir in range(0,n_ur): - for itheta in range(0,n_utheta): - timeshifts = tlib[:,ir,itheta] # Calculate the "steering vector" (a vector in frequency space, based on phase-shift) - a = np.exp(-1j*2*np.pi*target_f*timeshifts) # (a is a steering vector, to allign all traces with array centre) + Rxx[i1, i2] = 0 + else: + Rxx[i1, i2] = ( + np.conj(Pxx_all[curr_f_idx, i1]) * Pxx_all[curr_f_idx, i2] + ) + + for ir in range(n_sl): + for itheta in range(n_baz): + timeshifts = tlib[:, ir, itheta] + a = np.exp(-1j * 2 * np.pi * target_f * timeshifts) aconj = np.conj(a) - Pfreq[ii,ir,itheta]=np.dot(np.dot(aconj,Rxx),a) # Cross-correlation, with two timeshifts applied to push the two stations to the centre point. - # (This can also be seen as projecting Rxx onto a new basis.) + Pfreq[ii, ir, itheta] = np.dot(np.dot(aconj, Rxx), a) - # And remove any data where stations don't exist: - ###np.nan_to_num(Pfreq, copy=False, nan=0.0) # NOT SUPPORTED BY NUMBA SO DO OUTSIDE NUMBA - - # And append output to datastore: - Pfreq_all[win_idx,:,:,:] = Pfreq + Pfreq_all[win_idx, :, :, :] = Pfreq - return Pfreq_all + return Pfreq_all # @jit(nopython=True, parallel=True)#, nogil=True) -def _phase_associator_core_worker(peaks_Z, peaks_hor, bazis_Z, bazis_hor, bazi_tol, t_Z_secs_after_start, t_hor_secs_after_start, max_phase_sep_s): +def _phase_associator_core_worker( + peaks_Z, + peaks_hor, + bazis_Z, + bazis_hor, + bazi_tol, + t_Z_secs_after_start, + t_hor_secs_after_start, + max_phase_sep_s, + min_phase_sep_s, +): """Function to do the heavy lifting of the phase association.""" # Specify data stores: Z_hor_phase_pair_idxs = [] @@ -143,73 +204,109 @@ def _phase_associator_core_worker(peaks_Z, peaks_hor, bazis_Z, bazis_hor, bazi_t # Loop over phases, seeing if they meet phase association criteria: for i in range(len(peaks_Z)): if i % 1000 == 0: - print(i,"/",len(peaks_Z)) + print(i, "/", len(peaks_Z)) curr_peak_Z_idx = peaks_Z[i] for j in range(len(peaks_hor)): curr_peak_hor_idx = peaks_hor[j] # i. Check if phase arrivals are within specified time limits: - curr_t_phase_diff = t_hor_secs_after_start[curr_peak_hor_idx] - t_Z_secs_after_start[curr_peak_Z_idx] - if curr_t_phase_diff > 0: - if curr_t_phase_diff <= max_phase_sep_s: - # ii. Check if bazis for Z and horizontals current pick match: - if np.abs( bazis_Z[i] - bazis_hor[j] ) < bazi_tol: - match = True - # And deal with if they are close to North: - elif np.abs( bazis_Z[i] - bazis_hor[j] ) > (360. - bazi_tol): - match = True - else: - match = False - + curr_t_phase_diff = ( + t_hor_secs_after_start[curr_peak_hor_idx] + - t_Z_secs_after_start[curr_peak_Z_idx] + ) + if ( + (curr_t_phase_diff > 0) + and (curr_t_phase_diff <= max_phase_sep_s) + and (curr_t_phase_diff >= min_phase_sep_s) + ): + # calc bazi diff between Z and H + bazi_diff = min( + np.abs(bazis_Z[i] - bazis_hor[j]), + 360.0 - np.abs(bazis_Z[i] - bazis_hor[j]), + ) + # ii. Check if bazis for Z and horizontals current pick match: + if bazi_diff < bazi_tol: + match = [curr_peak_Z_idx, curr_peak_hor_idx] # And associate phases and create event data if a match is found: - if match: - # Append pair idxs to data store: - Z_hor_phase_pair_idxs.append([curr_peak_Z_idx, curr_peak_hor_idx]) - + # Append pair idxs to data store: + Z_hor_phase_pair_idxs.append(match) + return Z_hor_phase_pair_idxs -def _phase_associator(t_series_df_Z, t_series_df_hor, peaks_Z, peaks_hor, bazi_tol, filt_phase_assoc_by_max_power, max_phase_sep_s, min_event_sep_s, verbosity=0): +def _phase_associator( + t_series_df_Z, + t_series_df_hor, + peaks_Z, + peaks_hor, + bazi_tol, + filt_phase_assoc_by_max_power, + max_phase_sep_s, + min_phase_sep_s, + min_event_sep_s, + verbosity=0, +): """ Function to perform phase association for numba implementation. """ - # Setup events datastores: - events_df = pd.DataFrame() # Find back-azimuths associated with phase picks: - bazis_Z = t_series_df_Z['back_azi'].values[peaks_Z] - bazis_hor = t_series_df_hor['back_azi'].values[peaks_hor] + bazis_Z = t_series_df_Z["back_azi"].values[peaks_Z] + bazis_hor = t_series_df_hor["back_azi"].values[peaks_hor] - #------------------------------------------------------------------- - # Perform core phae association: + # ------------------------------------------------------------------- + # Perform core phase association: # Prep. data for numba format: if verbosity > 1: - print("Pre-processing time-series") - t_Z_secs_after_start = [] - for index, row in t_series_df_Z.iterrows(): - t_Z_secs_after_start.append(obspy.UTCDateTime(row['t']) - obspy.UTCDateTime(t_series_df_Z['t'][0])) - t_hor_secs_after_start = [] - for index, row in t_series_df_hor.iterrows(): - t_hor_secs_after_start.append(obspy.UTCDateTime(row['t']) - obspy.UTCDateTime(t_series_df_hor['t'][0])) - # Run function: + logger.info("Pre-processing time-series") + t_Z_secs_after_start = np.array(t_series_df_Z["t"] - t_series_df_Z["t"][0]) + t_hor_secs_after_start = np.array(t_series_df_hor["t"] - t_series_df_hor["t"][0]) + if verbosity > 1: - print("Performing phase association") - Z_hor_phase_pair_idxs = _phase_associator_core_worker(peaks_Z, peaks_hor, bazis_Z, bazis_hor, bazi_tol, t_Z_secs_after_start, t_hor_secs_after_start, max_phase_sep_s) + logger.info("Performing phase association") + Z_hor_phase_pair_idxs = _phase_associator_core_worker( + peaks_Z, + peaks_hor, + bazis_Z, + bazis_hor, + bazi_tol, + t_Z_secs_after_start, + t_hor_secs_after_start, + max_phase_sep_s, + min_phase_sep_s, + ) # Organise outputs into useful form: if verbosity > 1: - print("Writing events") + logger.info("Writing events") + + curr_events = { + "t1": [], + "t2": [], + "pow1": [], + "pow2": [], + "slow1": [], + "slow2": [], + "bazi1": [], + "bazi2": [], + } + for event_idx in range(len(Z_hor_phase_pair_idxs)): curr_peak_Z_idx = Z_hor_phase_pair_idxs[event_idx][0] curr_peak_hor_idx = Z_hor_phase_pair_idxs[event_idx][1] - curr_event_df = pd.DataFrame({'t1': [t_series_df_Z['t'][curr_peak_Z_idx]], 't2': [t_series_df_hor['t'][curr_peak_hor_idx]], - 'pow1': [t_series_df_Z['power'][curr_peak_Z_idx]], 'pow2': [t_series_df_hor['power'][curr_peak_hor_idx]], - 'slow1': [t_series_df_Z['slowness'][curr_peak_Z_idx]], 'slow2': [t_series_df_hor['slowness'][curr_peak_hor_idx]], - 'bazi1': [t_series_df_Z['back_azi'][curr_peak_Z_idx]], 'bazi2': [t_series_df_hor['back_azi'][curr_peak_hor_idx]]}) - events_df = events_df.append(curr_event_df) + curr_events["t1"].append(t_series_df_Z["t"][curr_peak_Z_idx]) + curr_events["t2"].append(t_series_df_hor["t"][curr_peak_hor_idx]) + curr_events["pow1"].append(t_series_df_Z["power"][curr_peak_Z_idx]) + curr_events["pow2"].append(t_series_df_hor["power"][curr_peak_hor_idx]) + curr_events["slow1"].append(t_series_df_Z["slowness"][curr_peak_Z_idx]) + curr_events["slow2"].append(t_series_df_hor["slowness"][curr_peak_hor_idx]) + curr_events["bazi1"].append(t_series_df_Z["back_azi"][curr_peak_Z_idx]) + curr_events["bazi2"].append(t_series_df_hor["back_azi"][curr_peak_hor_idx]) + + events_df = pd.DataFrame(curr_events) # And tidy: del t_Z_secs_after_start, t_hor_secs_after_start, Z_hor_phase_pair_idxs gc.collect() - #------------------------------------------------------------------- + # ------------------------------------------------------------------- - # And filter events to only output events with max. power within the max. phase window, + # And filter events to only output events with max. power within the max. phase window, # if speficifed by user: if filt_phase_assoc_by_max_power: # Only process if found some events: @@ -217,46 +314,49 @@ def _phase_associator(t_series_df_Z, t_series_df_hor, peaks_Z, peaks_hor, bazi_t # Calculate max power of P and S for each potential event: # events_overall_powers = events_df['pow1'].values + events_df['pow2'].values # Define datastores: - filt_events_df = pd.DataFrame() + filt_events_lst = [] # And loop over events, selecting only max. power events: tmp_count = 0 - for index, row in events_df.iterrows(): - tmp_count+=1 + for _, row in events_df.iterrows(): + tmp_count += 1 if tmp_count == 1: - tmp_df = pd.DataFrame() - tmp_df = tmp_df.append(row) + tmp_lst = [] + tmp_lst.append(row) else: + # Append event if phase within minimum event separation: - if obspy.UTCDateTime(row['t1']) - obspy.UTCDateTime(tmp_df['t1'].iloc[0]) < min_event_sep_s: + if ( + obspy.UTCDateTime(row["t1"]) - obspy.UTCDateTime(tmp_lst[0].t1) + < min_event_sep_s + ): # Append event to compare: - tmp_df = tmp_df.append(row) + tmp_lst.append(row) else: # Find best event from previous events: - combined_pows_tmp = tmp_df['pow1'].values + tmp_df['pow2'].values - max_power_idx = np.argmax(combined_pows_tmp) - filt_events_df = filt_events_df.append(tmp_df.iloc[max_power_idx]) + max_power_event = _find_max_power_event(tmp_lst) + filt_events_lst.append(max_power_event) # And start acrewing new events: - tmp_df = pd.DataFrame() - tmp_df = tmp_df.append(row) + tmp_lst = [] + tmp_lst.append(row) # And calculate highest power event for final window: - combined_pows_tmp = tmp_df['pow1'].values + tmp_df['pow2'].values - max_power_idx = np.argmax(combined_pows_tmp) - filt_events_df = filt_events_df.append(tmp_df.iloc[max_power_idx]) - + max_power_event = _find_max_power_event(tmp_lst) + filt_events_lst.append(max_power_event) + # Now make new DataFrame + filt_events_df = pd.DataFrame(filt_events_lst) # And sort indices: filt_events_df.reset_index(drop=True, inplace=True) # And remove duplicate S pick associations: # (using same max. power method) # Append summed powers, for sorting: - sum_pows = filt_events_df['pow1'].values + filt_events_df['pow2'].values - sum_pows_df = pd.DataFrame({'sum_pows': sum_pows}) - filt_events_df = filt_events_df.join(sum_pows_df) + sum_pows = filt_events_df["pow1"].values + filt_events_df["pow2"].values + filt_events_df["sum_pows"] = sum_pows # Remove t2 duplicates, keep highest summed power: - filt_events_df = filt_events_df.sort_values('sum_pows').drop_duplicates(subset='t2', keep='last') + filt_events_df.sort_values("sum_pows", inplace=True) + filt_events_df.drop_duplicates(subset="t2", keep="last", inplace=True) # And remove sum_pows column: - filt_events_df = filt_events_df.drop(columns=['sum_pows']) + filt_events_df.drop(columns=["sum_pows"], inplace=True) # And output df: events_df = filt_events_df.copy() @@ -266,10 +366,50 @@ def _phase_associator(t_series_df_Z, t_series_df_hor, peaks_Z, peaks_hor, bazi_t return events_df -def _submit_parallel_fast_freq_domain_array_proc(procnum, return_dict_Pfreq_all, data_curr_run, max_sl, fs, target_freqs, xx, yy, n_stations, n_t_samp, remove_autocorr): +def _find_max_power_event(events): + """ + Find the maximum power event from a list of events + + Parameters: + ---------- + events : list + list (of Dataframe Rows) of event + + """ + pow1_tmp = np.array([event.pow1 for event in events]) + pow2_tmp = np.array([event.pow2 for event in events]) + combined_pows_tmp = pow1_tmp + pow2_tmp + max_power_idx = np.argmax(combined_pows_tmp) + max_power_event = events[max_power_idx] + return max_power_event + + +def _submit_parallel_fast_freq_domain_array_proc( + procnum, + return_dict_Pfreq_all, + data_curr_run, + max_sl, + fs, + target_freqs, + xx, + yy, + n_stations, + n_t_samp, + remove_autocorr, +): """Function to submit parallel runs of _fast_freq_domain_array_proc() function.""" # Run function: - Pfreq_all_curr_run = _fast_freq_domain_array_proc(data_curr_run, max_sl, fs, target_freqs, xx, yy, n_stations, n_t_samp, remove_autocorr) + Pfreq_all_curr_run = _fast_freq_domain_array_proc( + data_curr_run, + max_sl, + fs, + target_freqs, + xx, + yy, + n_stations, + n_t_samp, + remove_autocorr, + ) # And return data return_dict_Pfreq_all[procnum] = Pfreq_all_curr_run @@ -277,11 +417,11 @@ def _submit_parallel_fast_freq_domain_array_proc(procnum, return_dict_Pfreq_all, def _calc_time_shift_from_array_cent(slow, bazi, x_rec, y_rec): """Calculates time shift of signal at receiver from array centre for stacking data. Note: All distances and velocities use km unless otherwise specified. - + Parameters ---------- slow : float - Slowness of arrival, in s/km. + Slowness of arrival, in s/km. bazi : float Back azimuth of arrival in degrees from . x_rec : @@ -292,7 +432,9 @@ def _calc_time_shift_from_array_cent(slow, bazi, x_rec, y_rec): # Calculate time-shift for receiver: # (in polar coord system, for consistency) bazi_rad = np.deg2rad(bazi) - time_shift_curr = x_rec*slow*np.sin(bazi_rad) + y_rec*slow*np.cos(bazi_rad) # (distance x slowness = distance / velocity = time) + time_shift_curr = x_rec * slow * np.sin(bazi_rad) + y_rec * slow * np.cos( + bazi_rad + ) # (distance x slowness = distance / velocity = time) return time_shift_curr @@ -303,55 +445,55 @@ def _create_stacked_data_st(st, Z_all, N_all, E_all): # For Z stacked: tr = st[0].copy() tr.stats.station = "STACK" - tr.stats.channel = st[0].stats.channel[0:2]+"Z" + tr.stats.channel = st[0].stats.channel[0:2] + "Z" tr.data = np.sum(Z_all, axis=1) composite_st.append(tr) # For Z mean: tr = st[0].copy() tr.stats.station = "MEAN" - tr.stats.channel = st[0].stats.channel[0:2]+"Z" + tr.stats.channel = st[0].stats.channel[0:2] + "Z" tr.data = np.mean(Z_all, axis=1) composite_st.append(tr) # For Z stdev: tr = st[0].copy() tr.stats.station = "STDEV" - tr.stats.channel = st[0].stats.channel[0:2]+"Z" + tr.stats.channel = st[0].stats.channel[0:2] + "Z" tr.data = np.std(Z_all, axis=1) composite_st.append(tr) # For N stacked: tr = st[0].copy() tr.stats.station = "STACK" - tr.stats.channel = st[0].stats.channel[0:2]+"N" + tr.stats.channel = st[0].stats.channel[0:2] + "N" tr.data = np.sum(N_all, axis=1) composite_st.append(tr) # For N mean: tr = st[0].copy() tr.stats.station = "MEAN" - tr.stats.channel = st[0].stats.channel[0:2]+"N" + tr.stats.channel = st[0].stats.channel[0:2] + "N" tr.data = np.mean(N_all, axis=1) composite_st.append(tr) # For N stdev: tr = st[0].copy() tr.stats.station = "STDEV" - tr.stats.channel = st[0].stats.channel[0:2]+"N" + tr.stats.channel = st[0].stats.channel[0:2] + "N" tr.data = np.std(N_all, axis=1) composite_st.append(tr) # For E stacked: tr = st[0].copy() tr.stats.station = "STACK" - tr.stats.channel = st[0].stats.channel[0:2]+"E" + tr.stats.channel = st[0].stats.channel[0:2] + "E" tr.data = np.sum(E_all, axis=1) composite_st.append(tr) # For E mean: tr = st[0].copy() tr.stats.station = "MEAN" - tr.stats.channel = st[0].stats.channel[0:2]+"E" + tr.stats.channel = st[0].stats.channel[0:2] + "E" tr.data = np.mean(E_all, axis=1) composite_st.append(tr) # For E stdev: tr = st[0].copy() tr.stats.station = "STDEV" - tr.stats.channel = st[0].stats.channel[0:2]+"E" + tr.stats.channel = st[0].stats.channel[0:2] + "E" tr.data = np.std(E_all, axis=1) composite_st.append(tr) del tr @@ -360,6 +502,52 @@ def _create_stacked_data_st(st, Z_all, N_all, E_all): return composite_st +def _create_phase_weighted_stack_st(st, Z_all, N_all, E_all, degree=1): + """Function to create stacked data st.""" + + Z_analytical_signal = hilbert(Z_all, axis=1) + N_analytical_signal = hilbert(N_all, axis=1) + E_analytical_signal = hilbert(E_all, axis=1) + # . Scipy hilbert transform returns the analytical signal with takes the form: + # . s(t) = x(t) + i * y(t) = A(t) * exp(i * phi(t)) + # . where x(t) is the original signal, y(t) is the Hilbert transform of the signal, + # . A(t) is the instantaneous amplitude (envelope) and phi(t) is the instantaneous phase. + + # . We want to calculate the phase stack, which is given by: + # . c(t) = |(1/N) * sum(exp(i * phi_k(t) | ^v + Z_inst_phase = Z_analytical_signal / np.abs(Z_analytical_signal) + N_inst_phase = N_analytical_signal / np.abs(N_analytical_signal) + E_inst_phase = E_analytical_signal / np.abs(E_analytical_signal) + + Z_coherence = np.absolute(np.mean(Z_inst_phase, axis=1)) + N_coherence = np.absolute(np.mean(N_inst_phase, axis=1)) + E_coherence = np.absolute(np.mean(E_inst_phase, axis=1)) + composite_st = obspy.Stream() + # For Z stacked: + tr = st[0].copy() + tr.stats.station = "PW-STACK" + tr.stats.channel = st[0].stats.channel[0:2] + "Z" + tr.data = np.mean(Z_all, axis=1) * (Z_coherence**degree) + print(tr.data.shape) + composite_st.append(tr) + # For N stacked: + tr = st[0].copy() + tr.stats.station = "PW-STACK" + tr.stats.channel = st[0].stats.channel[0:2] + "N" + tr.data = np.mean(N_all, axis=1) * (N_coherence**degree) + composite_st.append(tr) + # For E stacked: + tr = st[0].copy() + tr.stats.station = "PW-STACK" + tr.stats.channel = st[0].stats.channel[0:2] + "E" + tr.data = np.mean(E_all, axis=1) * (E_coherence**degree) + composite_st.append(tr) + del tr + gc.collect() + + return composite_st + + class setup_detection: """ Class to create detection object, for running array detection algorithm. @@ -368,13 +556,13 @@ class setup_detection: ---------- archivedir : str Path to data archive. Data archive must be of specific format: - /YEAR/JULDAY/YEARJULDAY_*STATION_COMP.* + /YEAR/MONTH/DAY/YYYYMMDDTHHMMSS_*STATION_COMP.* outdir : str Path to directory to save outputs to. stations_fname : str - Path to csv file containing station/receiver locations. Headers need + Path to csv file containing station/receiver locations. Headers need to be of format: Latitude Longitude Elevation Name. starttime : obspy UTCDateTime object @@ -385,12 +573,15 @@ class setup_detection: channels_to_use : list of strs (optional, default = ["??Z"]) List of channels to use for the processing (e.g. HHZ, ??Z or similar). - Note: Currently must be in form: ["??Z"] or ["??Z", "??N", "??E"] - or ["??Z", "??1", "??2"]. + Note: Currently must be in form: ["??Z"] or ["??Z", "??N", "??E"] + or ["??Z", "??1", "??2"]. Attributes ---------- + skip_existing : float + If True, then skip exisitng detection time series in outdir. Default is True. + freqmin : float If specified, lower frequency of bandpass filter, in Hz. Default is None. @@ -398,9 +589,9 @@ class setup_detection: If specified, upper frequency of bandpass filter, in Hz. Default is None. num_freqs : int - Number of discrete frequencies to use between and - in analysis. Default is 100. Note: Reducing this value increases - efficiency linearly, but costs in terms of absolute power. However, + Number of discrete frequencies to use between and + in analysis. Default is 100. Note: Reducing this value increases + efficiency linearly, but costs in terms of absolute power. However, SNR of power time-series remains approximately constant (to a point). Can affect slowness and bazi results, especially within noise, though. @@ -412,21 +603,21 @@ class setup_detection: Default value is 0.1 s. win_step_inc_s : float - The step increment for each window step. Units are seconds. Note: set this value - equal to for no overlap of windows. overlap between windows is given - by win_len_s - win_step_inc_s. For example, if win_step_inc_s = 3 x - then overlap will be 2 x , etc. Greater overlap gives higher frequency + The step increment for each window step. Units are seconds. Note: set this value + equal to for no overlap of windows. overlap between windows is given + by win_len_s - win_step_inc_s. For example, if win_step_inc_s = 3 x + then overlap will be 2 x , etc. Greater overlap gives higher frequency resolution, but at computational cost. Default value is 0.1 s. - + remove_autocorr : bool If True, then will remove autocorrelations. Default is True. norm_pre_stacking : bool - If True, normallises data before stacking. Similar to performing spectral + If True, normallises data before stacking. Similar to performing spectral whitening. Default is False mad_window_length_s : float - Length of time-window, in seconds, to calculate background Median Absolute + Length of time-window, in seconds, to calculate background Median Absolute Deviation (MAD) for triggering events. Default is 3600 s. mad_multiplier : int @@ -436,39 +627,48 @@ class setup_detection: min_event_sep_s : float Minimum separation between event detections, in seconds. Default = 1 s. - bazi_tol : float + bazi_tol : float The back-azimuth tolerance to associate incoming phases with the same event. - Various phases have to fulfil the criteria of having the same back-azimuth + Various phases have to fulfil the criteria of having the same back-azimuth for all phase arrivals, +/- . Units are degrees. Defaul is 20 degrees. max_phase_sep_s : float - The maximum time separation between individual phases. Units are seconds. Default + The maximum time separation between individual phases. Units are seconds. Default is 2.5 s. filt_phase_assoc_by_max_power : bool - If True, filters event phase association within a given window by max. power. I.e. - if True, will only pick highest amplitude phase arrivals on both vertical and + If True, filters event phase association within a given window by max. power. I.e. + if True, will only pick highest amplitude phase arrivals on both vertical and horizontal components, within the time window . Methods ------- """ - - def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, preload_fname=None, channels_to_use=["??Z"]): + + def __init__( + self, + archivedir, + outdir, + stations_fname, + starttime, + endtime, + preload_fname=None, + channels_to_use=["??Z"], + ): """Initiate the class object. Parameters ---------- archivedir : str Path to data archive. Data archive must be of specific format: - /YEAR/JULDAY/YEARJULDAY_*STATION_COMP.* + /YEAR/MONTH/DAY/YYYYMMDDTHHMMSS_*STATION_COMP.* outdir : str Path to directory to save outputs to. stations_fname : str - Path to csv file containing station/receiver locations. Headers need + Path to csv file containing station/receiver locations. Headers need to be of format: Latitude Longitude Elevation Name. starttime : obspy UTCDateTime object @@ -478,7 +678,7 @@ def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, prelo End time of data window to process for. preload_fname : str - Path to previously created detection class object. Optional. Default is + Path to previously created detection class object. Optional. Default is None, which means it doesn't load an existing file. channels_to_use : list of strs (optional, default = ["??Z"]) @@ -486,8 +686,8 @@ def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, prelo """ # Initialise input params: - self.archivedir = archivedir - self.outdir = outdir + self.archivedir = Path(archivedir) + self.outdir = Path(outdir) self.stations_fname = stations_fname self.starttime = starttime self.endtime = endtime @@ -502,12 +702,18 @@ def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, prelo # And define attributes: # For array processing: + self.skip_existing = True self.freqmin = None self.freqmax = None self.num_freqs = 100 + self.min_sl = 0 self.max_sl = 1.0 + self.n_sl = 51 + self.min_baz = 0 + self.max_baz = 360 + self.n_baz = 181 self.win_len_s = 0.1 - self.win_step_inc_s = 0.1 # (Note: Default is to step with no overlap) + self.win_step_inc_s = 0.1 # (Note: Default is to step with no overlap) self.remove_autocorr = True self.norm_pre_stacking = False # self.nproc = 1 @@ -515,11 +721,12 @@ def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, prelo # For detection: self.mad_window_length_s = 3600 self.mad_multiplier = 8 - self.min_event_sep_s = 1. - self.bazi_tol = 20. + self.min_event_sep_s = 1.0 + self.bazi_tol = 20.0 self.max_phase_sep_s = 2.5 + self.min_phase_sep_s = 0 self.filt_phase_assoc_by_max_power = True - self.calc_uncertainties = False + self.calc_uncertainties = False # For location: self.receiver_vp = None self.receiver_vs = None @@ -529,38 +736,210 @@ def __init__(self, archivedir, outdir, stations_fname, starttime, endtime, prelo if preload_fname: self.load(preload_fname) + def run_array_proc(self): + """Function to run core array processing. + Performed in frequency domain. Involves applying phase (equiv. to time) shift + for each frequency, over a range of specified slownesses. + Function inspured by work of D. Bowden (see Bowden et al. (2020)).""" + + # Find number of days to run array processing over + dt_start = self.starttime.date + dt_end = self.endtime.date + if (dt_end - dt_start).days == 0: + # round up to one day at minimum + ndays = 1 + else: + ndays = (dt_end - dt_start).days + + query_dates = [dt_start + datetime.timedelta(days=d) for d in range(0, ndays)] + for date in query_dates: + # Loop over dates within start/end range: + # Loop over channels: + for self.channel_curr in self.channels_to_use: + logger.info("=" * 60) + logger.info( + f"Processing data for day {date}, channel {self.channel_curr}" + ) + # And process for individual hours: + # (to reduce memory usage) + for hour in range(24): + # Loop over every hour in every day.. + # Make outfile + outfile = f"detection_t_series_{date.year:02d}{date.month:02d}{date.day:02d}_{hour:02d}00_ch{self.channel_curr[-1]}.csv" + if ((self.outdir / outfile).is_file()) & (self.skip_existing): + logger.warning(f"{outfile} exists in {self.outdir}") + logger.warning("Move to next hour") + self.out_fnames_array_proc.append(f"{self.outdir}/{outfile}") + continue + if ( + self.starttime + >= obspy.UTCDateTime( + year=date.year, month=date.month, day=date.day, hour=hour + ) + + 3600 + ): + continue + if self.endtime <= obspy.UTCDateTime( + year=date.year, month=date.month, day=date.day, hour=hour + ): + continue + + # Create datastores: + data_store = {"t": [], "power": [], "slowness": [], "back_azi": []} + + # Load data: + try: + st = self._load_data( + year=date.year, month=date.month, day=date.day, hour=hour + ) + except IndexError: + # And skip if no data: + logger.exception("Skipping hour as no data") + del st + gc.collect() + continue + + # And loop over minutes (to save on memory issues): + for minute in range(60): + # Check whether specified window is greater than a minute in duration: + if self.endtime - self.starttime > 60: + # Check time within specified run window: + if ( + self.starttime + >= obspy.UTCDateTime( + year=date.year, + month=date.month, + day=date.day, + hour=hour, + minute=minute, + ) + + 60 + ): + continue + if self.endtime <= obspy.UTCDateTime( + year=date.year, + month=date.month, + day=date.day, + hour=hour, + minute=minute, + ): + continue + elif self.starttime.minute != minute: + continue + + # Trim data: + st_trimmed = st.copy() + if self.win_len_s > self.win_step_inc_s: + self.win_pad_s = self.win_len_s + else: + self.win_pad_s = 0.0 + if self.endtime - self.starttime > 60: + st_trimmed.trim( + starttime=obspy.UTCDateTime( + year=date.year, + month=date.month, + day=date.day, + hour=hour, + minute=minute, + ), + endtime=obspy.UTCDateTime( + year=date.year, + month=date.month, + day=date.day, + hour=hour, + minute=minute, + ) + + 60 + + self.win_pad_s, + ) + else: + st_trimmed.trim( + starttime=self.starttime, + endtime=self.endtime + self.win_pad_s, + ) + time_this_minute_st = st_trimmed[0].stats.starttime + # Run array processing: + # (to get power in slowness space) + Psum_all = self._beamforming(st_trimmed) + del st_trimmed + gc.collect() + + # Calculate time-series outputs (for detection) from data: + t_series, powers, slownesses, back_azis = ( + self._find_time_series(Psum_all) + ) + # And append to data out: + t_series_out = [] + for t_serie in t_series: + t_series_out.append(str(time_this_minute_st + t_serie)) + + data_store["t"].extend(t_series_out) + data_store["power"].extend(powers) + data_store["slowness"].extend(slownesses) + data_store["back_azi"].extend(back_azis) + + # And clear memory: + del Psum_all, t_series, powers, slownesses, back_azis + gc.collect() + + # And save data out: + out_fname = os.path.join(self.outdir, outfile) + # make DataFrame "just-in-time" as it is more efficient this way + store_df = pd.DataFrame(data_store) + store_df.to_csv(out_fname, index=False) + + # And append fname to history: + self.out_fnames_array_proc.append(out_fname) + + return None def _setup_array_receiver_coords(self): - """Function to setup station receiver coords in correct format for + """Function to setup station receiver coords in correct format for array processing.""" - ref_lat = np.mean(self.stations_df['Latitude']) - ref_lon = np.mean(self.stations_df['Longitude']) + ref_lat = np.mean(self.stations_df["Latitude"]) + ref_lon = np.mean(self.stations_df["Longitude"]) # Calculate receiver locations in km grid format: xs = [] ys = [] for index, row in self.stations_df.iterrows(): - lon2 = row['Longitude'] - lat2 = row['Latitude'] + lon2 = row["Longitude"] + lat2 = row["Latitude"] # Calc. inter-station distance: - r, a, b = obspy.geodetics.base.gps2dist_azimuth(ref_lat, ref_lon, lat2, lon2) + r, a, b = obspy.geodetics.base.gps2dist_azimuth( + ref_lat, ref_lon, lat2, lon2 + ) xs.append(r * np.sin(np.deg2rad(a)) / 1000) ys.append(r * np.cos(np.deg2rad(a)) / 1000) - self.stations_df['x (km)'] = xs - self.stations_df['y (km)'] = ys - self.stations_df['z (km)'] = self.stations_df['Elevation'].values / 1000. - self.stations_df['ref_lat'] = np.ones(len(xs)) * ref_lat - self.stations_df['ref_lon'] = np.ones(len(xs)) * ref_lon + self.stations_df["x (km)"] = xs + self.stations_df["y (km)"] = ys + self.stations_df["z (km)"] = self.stations_df["Elevation"].values / 1000.0 + self.stations_df["ref_lat"] = np.ones(len(xs)) * ref_lat + self.stations_df["ref_lon"] = np.ones(len(xs)) * ref_lon # And calculate receiver locations in terms of array centre: - array_centre = np.array([np.mean(self.stations_df['x (km)']), np.mean(self.stations_df['y (km)']), np.mean(self.stations_df['z (km)'])]) - self.stations_df['x_array_coords_km'] = self.stations_df['x (km)'].values - array_centre[0] - self.stations_df['y_array_coords_km'] = self.stations_df['y (km)'].values - array_centre[1] - self.stations_df['z_array_coords_km'] = self.stations_df['z (km)'].values - array_centre[2] + array_centre = np.array( + [ + np.mean(self.stations_df["x (km)"]), + np.mean(self.stations_df["y (km)"]), + np.mean(self.stations_df["z (km)"]), + ] + ) + self.stations_df["x_array_coords_km"] = ( + self.stations_df["x (km)"].values - array_centre[0] + ) + self.stations_df["y_array_coords_km"] = ( + self.stations_df["y (km)"].values - array_centre[1] + ) + self.stations_df["z_array_coords_km"] = ( + self.stations_df["z (km)"].values - array_centre[2] + ) # And in polar coords from N: - self.stations_df['r_array_coords_km'], self.stations_df['theta_array_coords_deg'] = xy_to_rtheta(self.stations_df['x_array_coords_km'], - self.stations_df['y_array_coords_km']) - - + ( + self.stations_df["r_array_coords_km"], + self.stations_df["theta_array_coords_deg"], + ) = xy_to_rtheta( + self.stations_df["x_array_coords_km"], self.stations_df["y_array_coords_km"] + ) def find_min_max_array_sensitivity(self, vel_assumed=3.0): """Function to find array min and max sensitivities. @@ -573,65 +952,105 @@ def find_min_max_array_sensitivity(self, vel_assumed=3.0): inter_station_dists = [] # Loop over first set of stations: for index, row in self.stations_df.iterrows(): - lon1 = row['Longitude'] - lat1 = row['Latitude'] + lon1 = row["Longitude"] + lat1 = row["Latitude"] # Loop over stations again: for index, row in self.stations_df.iterrows(): - lon2 = row['Longitude'] - lat2 = row['Latitude'] + lon2 = row["Longitude"] + lat2 = row["Latitude"] # Calc. inter-station distance: r, a, b = obspy.geodetics.base.gps2dist_azimuth(lat1, lon1, lat2, lon2) inter_station_dists.append(np.abs(r) / 1000) inter_station_dists = np.array(inter_station_dists) inter_station_dists = inter_station_dists[inter_station_dists != 0] - print("="*60) + print("=" * 60) print("Min. inter-station distance:", np.min(inter_station_dists), "km") - print("Therefore, optimal sensitive higher freq.:", vel_assumed / np.min(inter_station_dists), "Hz") + print( + "Therefore, optimal sensitive higher freq.:", + vel_assumed / np.min(inter_station_dists), + "Hz", + ) print("Max. inter-station distance:", np.max(inter_station_dists), "km") - print("Therefore, optimal sensitive lower freq.:", vel_assumed / np.max(inter_station_dists), "Hz") - print("="*60) + print( + "Therefore, optimal sensitive lower freq.:", + vel_assumed / np.max(inter_station_dists), + "Hz", + ) + print("=" * 60) + + def _load_data(self, year, month, day, hour=None, norm=True): + """ + Function to load data. If no hour is specified the whole day will be read in. + Otherwise the hour of data will be loaded. + + Parameters + ---------- + year : int + year to load data for (yyyy) + month : int + month to load data for. Leading 0's will be added. + day : int + day to load data for. Leading 0's will be added. + hour : int, Optional + hour to load data for. Leading 0's will be added. - def _load_day_of_data(self, year, julday, hour=None): - """Function to load a day of data.""" + Returns: + ---------- + data : obspy.Stream + Waveform data for requested date and time. + """ # Load in data: - mseed_dir = os.path.join(self.archivedir, str(year), str(julday).zfill(3)) + mseed_dir = Path( + self.archivedir, str(year), str(month).zfill(2), str(day).zfill(2) + ) + # print(mseed_dir) st = obspy.Stream() for index, row in self.stations_df.iterrows(): - station = row['Name'] + # [J Asplet - think about replacing station DataFrame with StatonXML object] + station = row["Name"] for channel in self.channels_to_use: + if hour is None: + timestamp = f"{year:02d}{month:02d}{day:02d}T*" + else: + timestamp = f"{year:02d}{month:02d}{day:02d}T{hour:02d}0000" try: - if hour: - st_tmp = obspy.read(os.path.join(mseed_dir, ''.join((str(year), str(julday).zfill(3), "_", str(hour).zfill(2), "*", station, "*", channel, "*")))) - else: - st_tmp = obspy.read(os.path.join(mseed_dir, ''.join((str(year), str(julday).zfill(3), "_*", station, "*", channel, "*")))) + st_tmp = obspy.read( + f"{mseed_dir}/{timestamp}_{station}_{channel}.mseed" + ) for tr in st_tmp: st.append(tr) except: - print("No data for "+station+", channel = "+channel+". Skipping this data.") + logger.exception( + f"No data for {station}, channel = {channel}, timestamp {timestamp}. Skipping this data." + ) continue # Merge data: - st.detrend('demean') - st.merge(method=1, fill_value=0.) + st.detrend("demean") + st.detrend("linear") + st.merge(method=1, fill_value=0.0) # And apply filter: if self.freqmin: if self.freqmax: - st.filter('bandpass', freqmin=self.freqmin, freqmax=self.freqmax) - # And trim data, if some lies outside start and end times: + st.filter("bandpass", freqmin=self.freqmin, freqmax=self.freqmax) + # And trim data, if some lies outside start and end time of beamforming period: if self.starttime > st[0].stats.starttime: st.trim(starttime=self.starttime) if self.endtime < st[0].stats.endtime: st.trim(endtime=self.endtime) - return st - + return st + def _convert_st_to_np_data(self, st): """Function to convert data to numpy format for processing.""" - self.n_win = int(((st[0].stats.endtime - self.win_pad_s) - st[0].stats.starttime) / self.win_step_inc_s) # (Note: endtime - self.win_pad_s as pass extra padding via trimmed st) + self.n_win = int( + ((st[0].stats.endtime - self.win_pad_s) - st[0].stats.starttime) + / self.win_step_inc_s + ) # (Note: endtime - self.win_pad_s as pass extra padding via trimmed st) self.fs = st[0].stats.sampling_rate - self.n_t_samp = int(self.win_len_s * self.fs) # num samples in time - station_labels = self.stations_df['Name'].values - self.n_stations = len(station_labels) # num stations + self.n_t_samp = int(self.win_len_s * self.fs) # num samples in time + station_labels = self.stations_df["Name"].values + self.n_stations = len(station_labels) # num stations data = np.zeros((self.n_win, self.n_stations, self.n_t_samp)) for i in range(self.n_win): for j in range(self.n_stations): @@ -639,33 +1058,41 @@ def _convert_st_to_np_data(self, st): win_start_idx = i * int(self.win_step_inc_s * self.fs) win_end_idx = (i * int(self.win_step_inc_s * self.fs)) + self.n_t_samp try: - if win_end_idx < len(st.select(station=station, channel=self.channel_curr)[0].data): - data[i,j,:] = st.select(station=station, channel=self.channel_curr)[0].data[win_start_idx:win_end_idx] + if win_end_idx < len( + st.select(station=station, channel=self.channel_curr)[0].data + ): + data[i, j, :] = st.select( + station=station, channel=self.channel_curr + )[0].data[win_start_idx:win_end_idx] else: # Zero pad data (as insufficient data passed for final window) and print warning: - data[i,j,:] = 0. - print("Warning: Zero-padding as not enough data to fill window overlap ( for win_len_s =", self.win_len_s, "and win_step_inc_s =", self.win_step_inc_s, ")") + data[i, j, :] = 0.0 + logger.warning( + f"Warning: Zero-padding as not enough data to fill window overlap ( for win_len_s = {self.win_len_s}, and win_step_inc_s = {self.win_step_inc_s})" + ) except IndexError: # Deal with if a particular station has no data for given window: - data[i,j,:] = 0. - return data + data[i, j, :] = 0.0 + return data - def _stack_results(self, Pfreq_all): """Function to perform stacking of the results.""" - Psum_all = np.zeros((Pfreq_all.shape[0], Pfreq_all.shape[2], Pfreq_all.shape[3]), dtype=complex) + Psum_all = np.zeros( + (Pfreq_all.shape[0], Pfreq_all.shape[2], Pfreq_all.shape[3]), dtype=complex + ) # Loop over time windows: for i in range(Pfreq_all.shape[0]): if self.norm_pre_stacking: - Pfreq_norm_curr = Pfreq_all[i,:,:,:] / np.sum(np.abs(Pfreq_all[i,:,:,:]), axis=0) - Psum_all[i,:,:] = np.sum(Pfreq_norm_curr,axis=0) + Pfreq_norm_curr = Pfreq_all[i, :, :, :] / np.sum( + np.abs(Pfreq_all[i, :, :, :]), axis=0 + ) + Psum_all[i, :, :] = np.sum(Pfreq_norm_curr, axis=0) else: - Psum_all[i,:,:] = np.sum(Pfreq_all[i,:,:,:],axis=0) + Psum_all[i, :, :] = np.sum(Pfreq_all[i, :, :, :], axis=0) return Psum_all - def _find_time_series(self, Psum_all): - """Function to calculate beamforming time-series outputs, given + """Function to calculate beamforming time-series outputs, given a raw beamforming result. Note that the time-series timestamps are in the middle of the time- wimdows. @@ -674,13 +1101,19 @@ def _find_time_series(self, Psum_all): Returns time-series of coherency (power), slowness and back-azimuth. """ # Calcualte ux, uy: - ur = np.linspace(0, self.max_sl,Psum_all.shape[1]) - utheta = utheta = np.linspace(0,360-(360/Psum_all.shape[2]),Psum_all.shape[2]) - dur=ur[1]-ur[0] - dutheta=utheta[1]-utheta[0] + ur = np.linspace(self.min_sl, self.max_sl, Psum_all.shape[1]) + utheta = np.linspace( + self.min_baz, + self.max_baz - (self.max_baz / Psum_all.shape[2]), + Psum_all.shape[2], + ) # Create time-series: n_win_curr = Psum_all.shape[0] - t_series = np.arange(self.win_step_inc_s/2,(n_win_curr*self.win_step_inc_s) + (self.win_step_inc_s/2), self.win_step_inc_s) + t_series = np.arange( + self.win_step_inc_s / 2, + (n_win_curr * self.win_step_inc_s) + (self.win_step_inc_s / 2), + self.win_step_inc_s, + ) if len(t_series) > n_win_curr: t_series = t_series[0:n_win_curr] # And find power, slowness and back-azimuth: @@ -690,41 +1123,56 @@ def _find_time_series(self, Psum_all): # Loop over windows in time: for i in range(n_win_curr): # Calculate max. power: - powers[i] = np.max(np.abs(Psum_all[i,:,:])) + powers[i] = np.max(np.abs(Psum_all[i, :, :])) # Calculate slowness: - r_idx = np.where(Psum_all[i,:,:] == Psum_all[i,:,:].max())[0][0] - theta_idx = np.where(Psum_all[i,:,:] == Psum_all[i,:,:].max())[1][0] + r_idx = np.where(Psum_all[i, :, :] == Psum_all[i, :, :].max())[0][0] + theta_idx = np.where(Psum_all[i, :, :] == Psum_all[i, :, :].max())[1][0] slownesses[i] = ur[r_idx] # And calculate back-azimuth: back_azis[i] = utheta[theta_idx] - - return t_series, powers, slownesses, back_azis + return t_series, powers, slownesses, back_azis def _beamforming(self, st_trimmed, verbosity=0): - """Function to perform beamforming, given a stream of data for a specific + """Function to perform beamforming, given a stream of data for a specific time-window. Function is primarily called by run_array_proc(). Returns (stacked 2D power-slowness space data).""" # Run heavy array processing algorithm: # Specify various variables needed: # Make a linear spacing of frequencies. One might use periods, logspacing, etc.: # (Note: linspace much less noisy than logspace) - target_freqs = np.linspace(self.freqmin,self.freqmax,self.num_freqs) #np.logspace(self.freqmin,self.freqmax,self.num_freqs) + target_freqs = np.linspace( + self.freqmin, self.freqmax, self.num_freqs + ) # np.logspace(self.freqmin,self.freqmax,self.num_freqs) data = self._convert_st_to_np_data(st_trimmed) # Station locations: - xx = self.stations_df['x_array_coords_km'].values - yy = self.stations_df['y_array_coords_km'].values + xx = self.stations_df["x_array_coords_km"].values + yy = self.stations_df["y_array_coords_km"].values # And run: - if verbosity>1: - print("Performing run for",data.shape[0],"windows") + if verbosity > 1: + logger.info("Performing run for", data.shape[0], "windows") tic = time.time() - Pfreq_all = _fast_freq_domain_array_proc(data, self.max_sl, self.fs, target_freqs, xx, yy, - self.n_stations, self.n_t_samp, self.remove_autocorr) - if verbosity>1: + Pfreq_all = _fast_freq_domain_array_proc( + data, + self.min_sl, + self.max_sl, + self.n_sl, + self.min_baz, + self.max_baz, + self.n_baz, + self.fs, + target_freqs, + xx, + yy, + self.n_stations, + self.n_t_samp, + self.remove_autocorr, + ) + if verbosity > 1: toc = time.time() - print(toc-tic) + logger.info(f"runtime for _beamforming is {toc-tic}") # And tidy: - del data + del data gc.collect() # And remove any data where stations don't exist: @@ -737,331 +1185,301 @@ def _beamforming(self, st_trimmed, verbosity=0): return Psum_all - - def run_array_proc(self): - """Function to run core array processing. - Performed in frequency domain. Involves applying phase (equiv. to time) shift - for each frequency, over a range of specified slownesses. - Function inspured by work of D. Bowden (see Bowden et al. (2020)).""" - # Prep. stations df: - #self._setup_array_receiver_coords() - - # Loop over years: - for year in range(self.starttime.year, self.endtime.year+1): - # Loop over days: - for julday in range(1,367): - # Do some filtering for first and last years: - if year == self.starttime.year: - if julday < self.starttime.julday: - continue # Ignore day, as out of range - if year == self.endtime.year: - if julday > self.endtime.julday: - continue # Ignore day, as out of range - - # And process data: - - # Loop over channels: - for self.channel_curr in self.channels_to_use: - print("="*60) - print("Processing data for year "+str(year)+", day "+str(julday).zfill(3)+", channel "+self.channel_curr) - - # And process for individual hours: - # (to reduce memory usage) - for hour in range(24): - print("Processing for hour", str(hour).zfill(2)) - if self.starttime > obspy.UTCDateTime(year=year, julday=julday, hour=hour) + 3600: - continue - if self.endtime < obspy.UTCDateTime(year=year, julday=julday, hour=hour): - continue - - # Create datastore: - out_df = pd.DataFrame({'t': [], 'power': [], 'slowness': [], 'back_azi': []}) - - # Load data: - st = self._load_day_of_data(year, julday, hour=hour) - # starttime_this_day = obspy.UTCDateTime(year=year, julday=julday) - try: - starttime_this_st = st[0].stats.starttime - except IndexError: - # And skip if no data: - print("Skipping hour as no data") - del st - gc.collect() - continue - - # And loop over minutes (to save on memory issues): - for minute in range(60): - # Check whether specified window is greater than a minute in duration: - if self.endtime - self.starttime > 60: - # Check time within specified run window: - if self.starttime > obspy.UTCDateTime(year=year, julday=julday, hour=hour, minute=minute) + 60: - continue - if self.endtime < obspy.UTCDateTime(year=year, julday=julday, hour=hour, minute=minute): - continue - elif self.starttime.minute != minute: - continue - - # Trim data: - st_trimmed = st.copy() - if self.win_len_s > self.win_step_inc_s: - self.win_pad_s = self.win_len_s - else: - self.win_pad_s = 0. - if self.endtime - self.starttime > 60: - st_trimmed.trim(starttime=obspy.UTCDateTime(year=year, julday=julday, hour=hour, minute=minute), - endtime=obspy.UTCDateTime(year=year, julday=julday, hour=hour, minute=minute)+60+self.win_pad_s) - else: - st_trimmed.trim(starttime=self.starttime, endtime=self.endtime+self.win_pad_s) - - # Run array processing: - # (to get power in slowness space) - Psum_all = self._beamforming(st_trimmed) - del st_trimmed - gc.collect() - - # Calculate time-series outputs (for detection) from data: - t_series, powers, slownesses, back_azis = self._find_time_series(Psum_all) - - # And append to data out: - t_series_out = [] - if self.endtime - self.starttime > 60: - for t_serie in t_series: - t_series_out.append( str(starttime_this_st + (minute*60) + t_serie) ) - else: - for t_serie in t_series: - t_series_out.append( str(starttime_this_st + t_serie) ) - tmp_df = pd.DataFrame({'t': t_series_out, 'power': powers, 'slowness': slownesses, 'back_azi': back_azis}) - out_df = out_df.append(tmp_df) - out_df.reset_index(drop=True, inplace=True) - - # And save data out: - out_fname = os.path.join(self.outdir, ''.join(("detection_t_series_", str(year).zfill(4), str(julday).zfill(3), "_", - str(starttime_this_st.hour).zfill(2), "00", "_ch", self.channel_curr[-1], ".csv"))) - out_df.to_csv(out_fname, index=False) - - # And append fname to history: - self.out_fnames_array_proc.append(out_fname) - - # And clear memory: - del Psum_all, t_series, powers, slownesses, back_azis - gc.collect() - - return None - def _calculate_mad(self, x, scale=1.4826): """ Calculates the Median Absolute Deviation (MAD) of the input array x. - Outputs an array of scaled mean absolute deviation values for the input array, x, - scaled to provide an estimation of the standard deviation of the distribution. + Outputs an array of scaled mean absolute deviation values for the + input array, x, scaled to provide an estimation of the standard + deviation of the distribution. """ # Calculate median and mad values: mad = np.median(np.abs(x - np.median(x))) return scale * mad - - def _calc_uncertainties(self, events_df, t_series_df_Z, t_series_df_hor, verbosity=0): + def plot_polar_slowness_space( + self, beam_power, event_phase_arr_time, component, log=False + ): + + fig = plt.figure() + ax = fig.add_subplot(111, projection="polar") + rad = np.linspace(self.min_sl, self.max_sl, beam_power.shape[0]) + azm = np.linspace( + np.radians(self.min_baz), np.radians(self.max_baz), beam_power.shape[1] + ) + th, r = np.meshgrid(azm, rad) + ax.set_theta_offset(np.pi / 2) + ax.set_theta_direction(-1) + if log: + im = ax.pcolormesh(th, r, np.log(beam_power), cmap="magma") + else: + im = ax.pcolormesh(th, r, beam_power, cmap="magma") + + plt.colorbar(im) + plt.grid() + event_date_stamp = f"{event_phase_arr_time.year:04d}{event_phase_arr_time.month:02d}{event_phase_arr_time.day:02d}" + event_time_stamp = f"{event_phase_arr_time.hour:02d}{event_phase_arr_time.minute:02d}{event_phase_arr_time.second:02d}" + vesp_figpath = Path(self.outdir, "plots", "vespagrams") + vesp_figpath.mkdir( + parents=True, exist_ok=True + ) # makes plots/vespagrams if it doesnt exist + fig.savefig( + f"{vesp_figpath}/Detected_event_{event_date_stamp}_{event_time_stamp}_slow_spac_{component}.png", + dpi=600, + ) + plt.close() + + def _calc_uncertainties( + self, events_df, t_series_df_Z, t_series_df_hor, verbosity=0 + ): """Function to calculate uncertainties for phase-associated event detections. - Calculates uncertainty in t1, t2, slow1, slow2 and bazi1 and bazi2, - assuming Gaussian uncertainties. Uncertainties are defined as the full-width half - maximum (used due to more efficient optimisation than Gaussian fitting, but potentially + Calculates uncertainty in t1, t2, slow1, slow2 and bazi1 and bazi2, + assuming Gaussian uncertainties. Uncertainties are defined as the full-width half + maximum (used due to more efficient optimisation than Gaussian fitting, but potentially over-estimates error).""" # Do preliminary prep. once: # Define temparory datastore: - uncertainties_df = pd.DataFrame({'t1_err': [], 't2_err': [], 'slow1_err': [], 'slow2_err': [], - 'bazi1_err': [], 'bazi2_err': []}) + uncertainties_df = pd.DataFrame( + { + "t1_err": [], + "t2_err": [], + "slow1_err": [], + "slow2_err": [], + "bazi1_err": [], + "bazi2_err": [], + } + ) # Find max. timeshift (for determining beamforming window): - max_t_shift = self.max_sl * ( np.max(np.abs((self.stations_df['x_array_coords_km'].values))) - + np.max(np.abs((self.stations_df['x_array_coords_km'].values))) ) # (effectively d/v) - n_wins_for_max_t_shift = int(np.ceil(max_t_shift / self.win_step_inc_s)) #+ 1 #(+1 just to ensure that window is definitely wide enough) + max_t_shift = self.max_sl * ( + np.max(np.abs((self.stations_df["x_array_coords_km"].values))) + + np.max(np.abs((self.stations_df["x_array_coords_km"].values))) + ) # (effectively d/v) + n_wins_for_max_t_shift = int( + np.ceil(max_t_shift / self.win_step_inc_s) + ) # + 1 #(+1 just to ensure that window is definitely wide enough) # And loop over detected events, calculating uncertainty: count = 0 for index, row in events_df.iterrows(): if count % 10 == 0: if verbosity > 0: - print("Calculating uncertainty for event", count+1, "/", len(events_df)) + logger.info( + "Calculating uncertainty for event", + count + 1, + "/", + len(events_df), + ) # Load in data (if needed): # (done like this to avoid unnneccessary read ins, improving eff.) - event_phase_arr_time = obspy.UTCDateTime(row['t1']) + event_phase_arr_time = obspy.UTCDateTime(row["t1"]) if count == 0: - st = self._load_day_of_data(event_phase_arr_time.year, event_phase_arr_time.julday, hour=event_phase_arr_time.hour) + st = self._load_data( + event_phase_arr_time.year, + event_phase_arr_time.month, + event_phase_arr_time.day, + hour=event_phase_arr_time.hour, + ) # Find uncertainties: # ------- For vertical -------: # Time uncertainty: # Find FWHM for t1 pick: # (only use ascending currently (assume symetric pdf)) - t1_pick_idx = t_series_df_Z.index[t_series_df_Z['t'] == row['t1']][0] - Pxx_curr = t_series_df_Z.iloc[t1_pick_idx]['power'] + t1_pick_idx = t_series_df_Z.index[t_series_df_Z["t"] == row["t1"]][0] + Pxx_curr = t_series_df_Z.iloc[t1_pick_idx]["power"] idx_diff = 0 - while Pxx_curr > t_series_df_Z.iloc[t1_pick_idx]['power'] / 2.: - idx_diff+=1 - Pxx_curr = t_series_df_Z.iloc[t1_pick_idx+idx_diff]['power'] - t1_err = obspy.UTCDateTime(t_series_df_Z.iloc[t1_pick_idx+idx_diff]['t']) - obspy.UTCDateTime(t_series_df_Z.iloc[t1_pick_idx]['t']) - + while Pxx_curr > t_series_df_Z.iloc[t1_pick_idx]["power"] / 2.0: + idx_diff += 1 + Pxx_curr = t_series_df_Z.iloc[t1_pick_idx + idx_diff]["power"] + t1_err = obspy.UTCDateTime( + t_series_df_Z.iloc[t1_pick_idx + idx_diff]["t"] + ) - obspy.UTCDateTime(t_series_df_Z.iloc[t1_pick_idx]["t"]) + # Spatial uncertainty: # (slowness, bazi) # Perform beamforming again around event, to estimate bazi and slowness errs: # Get data: - event_phase_arr_time = obspy.UTCDateTime(row['t1']) + event_phase_arr_time = obspy.UTCDateTime(row["t1"]) # Reload data if needed: - if st[0].stats.starttime > event_phase_arr_time or st[0].stats.endtime < event_phase_arr_time: - st = self._load_day_of_data(event_phase_arr_time.year, event_phase_arr_time.julday, hour=event_phase_arr_time.hour) + if ( + st[0].stats.starttime > event_phase_arr_time + or st[0].stats.endtime < event_phase_arr_time + ): + st = self._load_data( + event_phase_arr_time.year, + event_phase_arr_time.month, + event_phase_arr_time.day, + hour=event_phase_arr_time.hour, + ) st_trimmed = st.copy() - st_trimmed.trim(starttime=event_phase_arr_time-((n_wins_for_max_t_shift+0.5)*self.win_len_s), - endtime=event_phase_arr_time+((n_wins_for_max_t_shift+0.5)*self.win_len_s)) # (Note: 0.5 as windows centred) + st_trimmed.trim( + starttime=event_phase_arr_time + - ((n_wins_for_max_t_shift + 0.5) * self.win_len_s), + endtime=event_phase_arr_time + + ((n_wins_for_max_t_shift + 0.5) * self.win_len_s), + ) # (Note: 0.5 as windows centred) # Run array processing: # (to get power in polar slowness space) # (Note that need to run for a number of windows, to allow for adequate shifting of data) - self.channel_curr = self.channels_to_use[0] # Do for vertical first + self.channel_curr = self.channels_to_use[0] # Do for vertical first Psum_all = self._beamforming(st_trimmed, verbosity=verbosity) - del st_trimmed + del st_trimmed gc.collect() # Find highest power slowness space for event: t_series, powers, slownesses, back_azis = self._find_time_series(Psum_all) max_idx = np.argmax(powers) - Psum_opt = np.abs(Psum_all[max_idx,:,:]) + Psum_opt = np.abs(Psum_all[max_idx, :, :]) # Find FWHM for slowness and bazi: - slow_idx_peak = np.where(Psum_opt==np.max(Psum_opt))[0][0] - bazi_idx_peak = np.where(Psum_opt==np.max(Psum_opt))[1][0] + slow_idx_peak = np.where(Psum_opt == np.max(Psum_opt))[0][0] + bazi_idx_peak = np.where(Psum_opt == np.max(Psum_opt))[1][0] # Slowness: # (go radially outwards for slowness, assumes symetric or sharper gradient inwards) Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak] idx_diff = 0 - while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.: - idx_diff+=1 - if slow_idx_peak+idx_diff < Psum_opt.shape[0]: - Pxx_curr = Psum_opt[slow_idx_peak+idx_diff, bazi_idx_peak] + while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.0: + idx_diff += 1 + if slow_idx_peak + idx_diff < Psum_opt.shape[0]: + Pxx_curr = Psum_opt[slow_idx_peak + idx_diff, bazi_idx_peak] else: - Pxx_curr = 0 # Force exit if reach slowness limits - dslow = self.max_sl / Psum_opt.shape[0] # Assumes linear slowness space + Pxx_curr = 0 # Force exit if reach slowness limits + dslow = self.max_sl / Psum_opt.shape[0] # Assumes linear slowness space slow1_err = idx_diff * dslow # Back-azimuth: # (go clockwise, assuming symetric) Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak] idx_diff = 0 - while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.: - idx_diff+=1 - if bazi_idx_peak+idx_diff < Psum_opt.shape[1]: - Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak+idx_diff] + while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.0: + idx_diff += 1 + if bazi_idx_peak + idx_diff < Psum_opt.shape[1]: + Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak + idx_diff] else: try: - Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak+idx_diff-Psum_opt.shape[1]] # (loop beyond 360 degrees) + Pxx_curr = Psum_opt[ + slow_idx_peak, bazi_idx_peak + idx_diff - Psum_opt.shape[1] + ] # (loop beyond 360 degrees) except IndexError: # Deal with 360 degree error: - Pxx_curr = 0 - idx_diff = Psum_opt.shape[1] + Pxx_curr = 0 + idx_diff = Psum_opt.shape[1] dbazi = 360 / Psum_opt.shape[1] + bazi1_err = idx_diff * dbazi # ------- End vertical ------- # Plot slowness space that used for uncertainty, if specified: - if verbosity > 1: - fig = plt.figure() - Axes3D(fig) - rad = np.linspace(0, self.max_sl, Psum_opt.shape[0]) - azm = np.linspace(0, 2 * np.pi, Psum_opt.shape[1]) - th, r = np.meshgrid(azm, rad) - ax = plt.subplot(projection="polar") - ax.set_theta_offset(np.pi/2) - ax.set_theta_direction(-1) - im = ax.pcolormesh(th, r, Psum_opt, cmap='inferno') - plt.colorbar(im) - plt.grid() - plt.show() - + if verbosity >= 1: + self.plot_polar_slowness_space( + Psum_opt, event_phase_arr_time, component="vert" + ) # ------- For horizontal -------: # And find FWHM for t2 pick: # (only use ascending currently (assume symetric pdf)) - t2_pick_idx = t_series_df_hor.index[t_series_df_hor['t'] == row['t2']][0] - Pxx_curr = t_series_df_hor.iloc[t2_pick_idx]['power'] + t2_pick_idx = t_series_df_hor.index[t_series_df_hor["t"] == row["t2"]][0] + Pxx_curr = t_series_df_hor.iloc[t2_pick_idx]["power"] idx_diff = 0 - while Pxx_curr > t_series_df_hor.iloc[t2_pick_idx]['power'] / 2.: - idx_diff+=1 - Pxx_curr = t_series_df_hor.iloc[t2_pick_idx+idx_diff]['power'] - t2_err = obspy.UTCDateTime(t_series_df_hor.iloc[t2_pick_idx+idx_diff]['t']) - obspy.UTCDateTime(t_series_df_hor.iloc[t2_pick_idx]['t']) + while Pxx_curr > t_series_df_hor.iloc[t2_pick_idx]["power"] / 2.0: + idx_diff += 1 + Pxx_curr = t_series_df_hor.iloc[t2_pick_idx + idx_diff]["power"] + t2_err = obspy.UTCDateTime( + t_series_df_hor.iloc[t2_pick_idx + idx_diff]["t"] + ) - obspy.UTCDateTime(t_series_df_hor.iloc[t2_pick_idx]["t"]) # Spatial uncertainty: # (slowness, bazi) # Perform beamforming again around event, to estimate bazi and slowness errs: # Get data: - event_phase_arr_time = obspy.UTCDateTime(row['t2']) + event_phase_arr_time = obspy.UTCDateTime(row["t2"]) # Reload data if needed: - if st[0].stats.starttime > event_phase_arr_time or st[0].stats.endtime < event_phase_arr_time: - st = self._load_day_of_data(event_phase_arr_time.year, event_phase_arr_time.julday, hour=event_phase_arr_time.hour) + if ( + st[0].stats.starttime > event_phase_arr_time + or st[0].stats.endtime < event_phase_arr_time + ): + st = self._load_data( + event_phase_arr_time.year, + event_phase_arr_time.month, + event_phase_arr_time.day, + hour=event_phase_arr_time.hour, + ) st_trimmed = st.copy() - st_trimmed.trim(starttime=event_phase_arr_time-((n_wins_for_max_t_shift+0.5)*self.win_len_s), - endtime=event_phase_arr_time+((n_wins_for_max_t_shift+0.5)*self.win_len_s)) # (Note: 0.5 as windows centred) + st_trimmed.trim( + starttime=event_phase_arr_time + - ((n_wins_for_max_t_shift + 0.5) * self.win_len_s), + endtime=event_phase_arr_time + + ((n_wins_for_max_t_shift + 0.5) * self.win_len_s), + ) # (Note: 0.5 as windows centred) # Run array processing: # (to get power in polar slowness space) # (Note that need to run for a number of windows, to allow for adequate shifting of data) - self.channel_curr = self.channels_to_use[1] # Do for vertical first + self.channel_curr = self.channels_to_use[1] # Do for vertical first Psum_all_N = self._beamforming(st_trimmed, verbosity=verbosity) - self.channel_curr = self.channels_to_use[2] # Do for vertical first + self.channel_curr = self.channels_to_use[2] # Do for vertical first Psum_all_E = self._beamforming(st_trimmed, verbosity=verbosity) - del st_trimmed + del st_trimmed gc.collect() Psum_all_NE = Psum_all_N + Psum_all_E # Find highest power slowness space for event: - t_series, powers, slownesses, back_azis = self._find_time_series(Psum_all_NE) + t_series, powers, slownesses, back_azis = self._find_time_series( + Psum_all_NE + ) max_idx = np.argmax(powers) - Psum_opt = np.abs(Psum_all_NE[max_idx,:,:]) + Psum_opt = np.abs(Psum_all_NE[max_idx, :, :]) # Find FWHM for slowness and bazi: - slow_idx_peak = np.where(Psum_opt==np.max(Psum_opt))[0][0] - bazi_idx_peak = np.where(Psum_opt==np.max(Psum_opt))[1][0] + slow_idx_peak = np.where(Psum_opt == np.max(Psum_opt))[0][0] + bazi_idx_peak = np.where(Psum_opt == np.max(Psum_opt))[1][0] # Slowness: # (go radially outwards for slowness, assumes symetric or sharper gradient inwards) Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak] idx_diff = 0 - while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.: - idx_diff+=1 - if slow_idx_peak+idx_diff < Psum_opt.shape[0]: - Pxx_curr = Psum_opt[slow_idx_peak+idx_diff, bazi_idx_peak] + while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.0: + idx_diff += 1 + if slow_idx_peak + idx_diff < Psum_opt.shape[0]: + Pxx_curr = Psum_opt[slow_idx_peak + idx_diff, bazi_idx_peak] else: - Pxx_curr = 0 # Force exit if reach slowness limits - dslow = self.max_sl / Psum_opt.shape[0] # Assumes linear slowness space + Pxx_curr = 0 # Force exit if reach slowness limits + dslow = self.max_sl / Psum_opt.shape[0] # Assumes linear slowness space slow2_err = idx_diff * dslow # Back-azimuth: # (go clockwise, assuming symetric) Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak] idx_diff = 0 - while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.: - idx_diff+=1 - if bazi_idx_peak+idx_diff < Psum_opt.shape[1]: - Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak+idx_diff] + while Pxx_curr > Psum_opt[slow_idx_peak, bazi_idx_peak] / 2.0: + idx_diff += 1 + if bazi_idx_peak + idx_diff < Psum_opt.shape[1]: + Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak + idx_diff] else: try: - Pxx_curr = Psum_opt[slow_idx_peak, bazi_idx_peak+idx_diff-Psum_opt.shape[1]] # (loop beyond 360 degrees) + Pxx_curr = Psum_opt[ + slow_idx_peak, bazi_idx_peak + idx_diff - Psum_opt.shape[1] + ] # (loop beyond 360 degrees) except IndexError: # Deal with 360 degree error: - Pxx_curr = 0 + Pxx_curr = 0 idx_diff = Psum_opt.shape[1] - dbazi = 360. / Psum_opt.shape[1] + dbazi = 360.0 / Psum_opt.shape[1] bazi2_err = idx_diff * dbazi # ------- End horizontal ------- - - # Plot slowness space that used for uncertainty, if specified: - if verbosity > 1: - fig = plt.figure() - Axes3D(fig) - rad = np.linspace(0, self.max_sl, Psum_opt.shape[0]) - azm = np.linspace(0, 2 * np.pi, Psum_opt.shape[1]) - th, r = np.meshgrid(azm, rad) - ax = plt.subplot(projection="polar") - ax.set_theta_offset(np.pi/2) - ax.set_theta_direction(-1) - im = ax.pcolormesh(th, r, Psum_opt, cmap='inferno') - plt.colorbar(im) - plt.grid() - plt.show() + # Plot slowness space that used for uncertainty, if specified: + if verbosity >= 1: + self.plot_polar_slowness_space( + Psum_opt, event_phase_arr_time, component="horz" + ) # And append data to overall uncertainties df: - uncertainties_df_curr = pd.DataFrame({'t1_err': [t1_err], 't2_err': [t2_err], 'slow1_err': [slow1_err], - 'slow2_err': [slow2_err], 'bazi1_err': [bazi1_err], 'bazi2_err': [bazi2_err]}) - uncertainties_df = pd.concat([uncertainties_df, uncertainties_df_curr], ignore_index=True) + uncertainties_df_curr = pd.DataFrame( + { + "t1_err": [t1_err], + "t2_err": [t2_err], + "slow1_err": [slow1_err], + "slow2_err": [slow2_err], + "bazi1_err": [bazi1_err], + "bazi2_err": [bazi2_err], + } + ) + uncertainties_df = pd.concat( + [uncertainties_df, uncertainties_df_curr], ignore_index=True + ) # And update count: - count+=1 + count += 1 # And add uncertainties to events_df: events_df = events_df.reset_index(drop=True) @@ -1069,57 +1487,95 @@ def _calc_uncertainties(self, events_df, t_series_df_Z, t_series_df_hor, verbosi events_df = pd.concat([events_df, uncertainties_df], axis=1) return events_df - - - def detect_events(self, verbosity=0): - """Function to detect events, based on the power time-series generated - by run_array_proc(). Note: Currently, only Median Absolute Deviation + + def detect_events(self, verbosity=0, fnames=None): + """Function to detect events, based on the power time-series generated + by run_array_proc(). Note: Currently, only Median Absolute Deviation triggering is implemented. Key attributes used are: - mad_window_length_s - mad_multiplier - min_event_sep_s """ - print("Note: not yet implemented.") + # print("Note: not yet implemented.") # Create datastore: events_df_all = pd.DataFrame() # Loop over array proc outdir data: - for fname in glob.glob(os.path.join(self.outdir, "detection_t_series_*_chZ.csv")): - f_uid = fname[-20:-8] + if fnames is None: + fnames = glob.glob( + os.path.join(self.outdir, "detection_t_series_*_chZ.csv") + ) + for fname in sorted(fnames): + f_uid = fname[-21:-8] # Check if in list to process: if fname in self.out_fnames_array_proc: # And load in data: # Read in vertical data: - t_series_df_Z = pd.read_csv(fname) + t_series_df_Z = pd.read_csv( + fname, converters={"t": lambda x: obspy.UTCDateTime(x)} + ) # And read in horizontals: try: - fname_N = os.path.join(self.outdir, ''.join(( "detection_t_series_", f_uid, "_chN.csv" ))) - t_series_df_N = pd.read_csv(fname_N) + fname_N = os.path.join( + self.outdir, "".join(("detection_t_series_", f_uid, "_chN.csv")) + ) + t_series_df_N = pd.read_csv( + fname_N, converters={"t": lambda x: obspy.UTCDateTime(x)} + ) except FileNotFoundError: - fname_N = os.path.join(self.outdir, ''.join(( "detection_t_series_", f_uid, "_ch1.csv" ))) - t_series_df_N = pd.read_csv(fname_N) + fname_N = os.path.join( + self.outdir, "".join(("detection_t_series_", f_uid, "_ch1.csv")) + ) + t_series_df_N = pd.read_csv( + fname_N, converters={"t": lambda x: obspy.UTCDateTime(x)} + ) try: - fname_E = os.path.join(self.outdir, ''.join(( "detection_t_series_", f_uid, "_chE.csv" ))) - t_series_df_E = pd.read_csv(fname_E) + fname_E = os.path.join( + self.outdir, "".join(("detection_t_series_", f_uid, "_chE.csv")) + ) + t_series_df_E = pd.read_csv( + fname_E, converters={"t": lambda x: obspy.UTCDateTime(x)} + ) except FileNotFoundError: - fname_E = os.path.join(self.outdir, ''.join(( "detection_t_series_", f_uid, "_ch2.csv" ))) - t_series_df_E = pd.read_csv(fname_E) + fname_E = os.path.join( + self.outdir, "".join(("detection_t_series_", f_uid, "_ch2.csv")) + ) + t_series_df_E = pd.read_csv( + fname_E, converters={"t": lambda x: obspy.UTCDateTime(x)} + ) else: - continue # Skip file, as not previously been processed. + logger.warning(f"fname {fname} not in fname_array_proc list") + logger.warning( + f"fname_array_proc first entry looks like this {self.out_fnames_array_proc[0]}" + ) + continue # Skip file, as not previously been processed. # And check to see that t-series exists within file: if len(t_series_df_Z) == 0: - continue + continue if len(t_series_df_N) == 0: - continue + continue if len(t_series_df_E) == 0: - continue + continue # Check if all inputs are same length, and if not, skip file: if not len(t_series_df_Z) == len(t_series_df_N) == len(t_series_df_E): - print("Warning: Files with f uid", f_uid, - "are not of equal length. Therefore using shortest length (will miss some data).") - print("( Lengths are", len(t_series_df_Z) , len(t_series_df_N) , len(t_series_df_E), ")") - min_len = np.min(np.array([len(t_series_df_Z), len(t_series_df_N), len(t_series_df_E)])) + logger.warning( + "Warning: Files with f uid", + f_uid, + "are not of equal length. Therefore using shortest length (will miss some data).", + ) + logger.warning( + "( Lengths are", + len(t_series_df_Z), + len(t_series_df_N), + len(t_series_df_E), + ")", + ) + min_len = np.min( + np.array( + [len(t_series_df_Z), len(t_series_df_N), len(t_series_df_E)] + ) + ) t_series_df_Z = t_series_df_Z.iloc[:min_len] t_series_df_N = t_series_df_N.iloc[:min_len] t_series_df_E = t_series_df_E.iloc[:min_len] @@ -1127,81 +1583,187 @@ def detect_events(self, verbosity=0): # Combine horizontals: # (Using RMS of N and E signals for slowness and average for BAZI) t_series_df_hor = t_series_df_N.copy() - t_series_df_hor["power"] = np.sqrt(t_series_df_N["power"].values**2 + t_series_df_E["power"].values**2) - NE_Pxx_max = np.max(np.concatenate((t_series_df_N["power"].values, t_series_df_E["power"].values))) + t_series_df_hor["power"] = np.sqrt( + t_series_df_N["power"].values ** 2 + t_series_df_E["power"].values ** 2 + ) + NE_Pxx_max = np.max( + np.concatenate( + (t_series_df_N["power"].values, t_series_df_E["power"].values) + ) + ) N_weighting = t_series_df_N["power"].values / NE_Pxx_max E_weighting = t_series_df_E["power"].values / NE_Pxx_max - t_series_df_hor["slowness"] = np.sqrt(np.average(np.vstack((t_series_df_N['slowness']**2, t_series_df_E['slowness']**2)), - axis=0, weights=np.vstack((N_weighting, - E_weighting)))) # Weighted mean (weighted by power) - t_series_df_hor["back_azi"] = np.average(np.vstack((t_series_df_N['back_azi'], t_series_df_E['back_azi'])), - axis=0, weights=np.vstack((N_weighting, - E_weighting))) # Weighted mean (weighted by power) - print("(Weighted horizontal slowness and back-azi using power)") + t_series_df_hor["slowness"] = np.sqrt( + np.average( + np.vstack( + (t_series_df_N["slowness"] ** 2, t_series_df_E["slowness"] ** 2) + ), + axis=0, + weights=np.vstack((N_weighting, E_weighting)), + ) + ) # Weighted mean (weighted by power) + t_series_df_hor["back_azi"] = np.average( + np.vstack((t_series_df_N["back_azi"], t_series_df_E["back_azi"])), + axis=0, + weights=np.vstack((N_weighting, E_weighting)), + ) # Weighted mean (weighted by power) + # print("(Weighted horizontal slowness and back-azi using power)") del N_weighting, E_weighting, t_series_df_N, t_series_df_E gc.collect() # Calculate pick thresholds: - mad_pick_threshold_Z = np.median(t_series_df_Z['power'].values) + (self.mad_multiplier * self._calculate_mad(t_series_df_Z['power'])) - mad_pick_threshold_hor = np.median(t_series_df_hor['power'].values) + (self.mad_multiplier * self._calculate_mad(t_series_df_hor['power'])) - - # Get phase picks: - min_pick_dist = int(self.min_event_sep_s / (obspy.UTCDateTime(t_series_df_Z['t'][1]) - obspy.UTCDateTime(t_series_df_Z['t'][0]))) - peaks_Z, _ = find_peaks(t_series_df_Z['power'].values, height=mad_pick_threshold_Z, distance=min_pick_dist) - peaks_hor, _ = find_peaks(t_series_df_hor['power'].values, height=mad_pick_threshold_hor, distance=min_pick_dist) + mad_pick_threshold_Z = moving_window_mad( + t_series_df_Z["power"].values, + self.mad_window_length_s, + self.mad_multiplier, + ) + + mad_pick_threshold_hor = moving_window_mad( + t_series_df_hor["power"].values, + self.mad_window_length_s, + self.mad_multiplier, + ) - # Phase assoicate by BAZI threshold and max. power: - events_df = _phase_associator(t_series_df_Z, t_series_df_hor, peaks_Z, peaks_hor, - self.bazi_tol, self.filt_phase_assoc_by_max_power, self.max_phase_sep_s, self.min_event_sep_s, verbosity=verbosity) + # Get phase picks: + min_pick_dist = int( + self.min_event_sep_s + / ( + obspy.UTCDateTime(t_series_df_Z["t"][1]) + - obspy.UTCDateTime(t_series_df_Z["t"][0]) + ) + ) + peaks_Z, _ = find_peaks( + t_series_df_Z["power"].values, + height=mad_pick_threshold_Z, + distance=min_pick_dist, + prominence=mad_pick_threshold_Z, + ) + peaks_hor, _ = find_peaks( + t_series_df_hor["power"].values, + height=mad_pick_threshold_hor, + distance=min_pick_dist, + prominence=mad_pick_threshold_hor, + ) + print( + f"Found {len(peaks_Z)} P-phase picks and {len(peaks_hor)} S-phase picks for file with uid {f_uid}" + ) + # Phase associate by BAZI threshold and max. power: + events_df = _phase_associator( + t_series_df_Z, + t_series_df_hor, + peaks_Z, + peaks_hor, + self.bazi_tol, + self.filt_phase_assoc_by_max_power, + self.max_phase_sep_s, + self.min_phase_sep_s, + self.min_event_sep_s, + verbosity=verbosity, + ) # Find uncertainties (in time, bazi, slowness): if self.calc_uncertainties: - events_df = self._calc_uncertainties(events_df, t_series_df_Z, t_series_df_hor, verbosity=verbosity) - - # Append to datastore: - events_df_all = events_df_all.append(events_df) - + events_df = self._calc_uncertainties( + events_df, t_series_df_Z, t_series_df_hor, verbosity=verbosity + ) # Plot detected, phase-associated picks: if verbosity > 1: - print("="*40) - print("Event phase associations:") - print(events_df) - print("="*40) - fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(6,4)) + # print("="*40) + logger.info("Event phase associations:") + # print(events_df) + # print("="*40) + time_after_startZ = np.array(t_series_df_Z["t"]) - t_series_df_Z["t"][0] + time_after_startH = ( + np.array(t_series_df_hor["t"]) - t_series_df_hor["t"][0] + ) + + fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(9, 6)) # Plot power: - ax[0].plot(t_series_df_Z['t'], t_series_df_Z['power'], label="Vertical power") - ax[0].plot(t_series_df_hor['t'], t_series_df_hor['power'], label="Horizontal power") + ax[0].plot( + time_after_startZ, t_series_df_Z["power"], label="Vertical power" + ) + ax[0].plot( + time_after_startH, + t_series_df_hor["power"], + label="Horizontal power", + ) # Plot slowness: - ax[1].plot(t_series_df_Z['t'], t_series_df_Z['slowness'], label="Vertical slowness") - ax[1].plot(t_series_df_hor['t'], t_series_df_hor['slowness'], label="Horizontal slowness") + ax[1].plot( + time_after_startZ, + t_series_df_Z["slowness"], + label="Vertical slowness", + ) + ax[1].plot( + time_after_startH, + t_series_df_hor["slowness"], + label="Horizontal slowness", + ) # Plot back-azimuth: - ax[2].plot(t_series_df_Z['t'], t_series_df_Z['back_azi'], label="Vertical back-azimuth") - ax[2].plot(t_series_df_hor['t'], t_series_df_hor['back_azi'], label="Horizontal back-azimuth") - if len(events_df_all) > 0: - ax[0].scatter(events_df_all['t1'], np.ones(len(events_df_all))*np.max(t_series_df_Z['power']), c='r', label="P phase picks") - ax[0].scatter(events_df_all['t2'], np.ones(len(events_df_all))*np.max(t_series_df_Z['power']), c='b', label="S phase picks") + ax[2].plot( + time_after_startZ, + t_series_df_Z["back_azi"], + label="Vertical back-azimuth", + ) + ax[2].plot( + time_after_startH, + t_series_df_hor["back_azi"], + label="Horizontal back-azimuth", + ) + if len(events_df) > 0: + events_t1_after_start = ( + np.array(events_df["t1"]) - t_series_df_Z["t"][0] + ) + events_t2_after_start = ( + np.array(events_df["t2"]) - t_series_df_Z["t"][0] + ) + ax[0].scatter( + events_t1_after_start, + np.ones(len(events_df)) * np.max(t_series_df_Z["power"]), + c="r", + label="P phase picks", + ) + ax[0].scatter( + events_t2_after_start, + np.ones(len(events_df)) * np.max(t_series_df_Z["power"]), + c="b", + label="S phase picks", + ) else: - print("No events to plot.") + logger.info("No events to plot.") ax[0].legend() ax[2].set_xlabel("Time") ax[0].set_ylabel("Power (arb. units)") ax[1].set_ylabel("Slowness ($km$ $s^{-1}$)") ax[2].set_ylabel("Back-azimuth ($^o$)") - # plt.gca().yaxis.set_major_locator(MaxNLocator(5)) + fig.suptitle( + f"Beamforming time-series for {t_series_df_Z['t'][0]} to {t_series_df_Z['t'].values[-1]} UTC" + ) + # plt.gca().yaxis.set_major_locator(MaxNLocator(5)) for i in range(3): ax[i].xaxis.set_major_locator(plt.MaxNLocator(3)) + figpath = Path(self.outdir, "plots", "detection_t_series") + figpath.mkdir(parents=True, exist_ok=True) + fig.savefig(f"{figpath}/Phase_association_{f_uid}.png", dpi=600) plt.show() - + + # Append to datastore: + events_df_all = pd.concat([events_df_all, events_df]) + events_df_all.reset_index(drop=True, inplace=True) return events_df_all - - - def create_location_LUTs(self, oneD_vel_model_z_df, extent_x_m=4000, dxz=[100,100], array_centre_xz=[0, 0]): + + def create_location_LUTs( + self, + oneD_vel_model_z_df, + extent_x_m=4000, + dxz=[100, 100], + array_centre_xz=[0, 0], + ): """Function to create lookup tables used for location. Lookup tables created are: - P travel-times - S travel-times. - P inclination angles. - S inclination angles. - Note: Array centre is typically defined as the + Note: Array centre is typically defined as the Parameters ---------- @@ -1214,13 +1776,13 @@ def create_location_LUTs(self, oneD_vel_model_z_df, extent_x_m=4000, dxz=[100,10 Extent of lookup table horizontal. Units are metres. Default is 4000 m. array_centre_xz : list - Path to csv file containing station/receiver locations. Headers need - to be of format: Latitude Longitude Elevation Name. No need to change unless - the user has a particularly good reason. Takes a list of length two, containing + Path to csv file containing station/receiver locations. Headers need + to be of format: Latitude Longitude Elevation Name. No need to change unless + the user has a particularly good reason. Takes a list of length two, containing the array centre in x and z, in metres from the LUT grid origin. dxz : list - List of two floats, defining the spatial spacing of the nodes in the LUTs in x + List of two floats, defining the spatial spacing of the nodes in the LUTs in x and z. Units are metres. Default is [100, 100]. Returns @@ -1236,26 +1798,34 @@ def create_location_LUTs(self, oneD_vel_model_z_df, extent_x_m=4000, dxz=[100,10 # Create 2D LUTs: # (for P and S travel-times, and inclination angle) # And get travel times: - trav_times_grid_P, trav_times_grid_S, theta_grid_P, theta_grid_S, vel_model_x_labels, vel_model_z_labels = lookup_table_manager.create_2D_LUT(oneD_vel_model_z_df, - array_centre_xz, extent_x_m=extent_x_m, dxz=dxz) + ( + trav_times_grid_P, + trav_times_grid_S, + theta_grid_P, + theta_grid_S, + vel_model_x_labels, + vel_model_z_labels, + ) = lookup_table_manager.create_2D_LUT( + oneD_vel_model_z_df, array_centre_xz, extent_x_m=extent_x_m, dxz=dxz + ) # And save outputs: LUT_outdir = os.path.join(self.outdir, "LUT") os.makedirs(LUT_outdir, exist_ok=True) LUTs_dict = {} - LUTs_dict['trav_times_grid_P'] = trav_times_grid_P - LUTs_dict['trav_times_grid_S'] = trav_times_grid_S - LUTs_dict['theta_grid_P'] = theta_grid_P - LUTs_dict['theta_grid_S'] = theta_grid_S - LUTs_dict['vel_model_x_labels'] = vel_model_x_labels - LUTs_dict['vel_model_z_labels'] = vel_model_z_labels - LUTs_dict['oneD_vel_model_z_df'] = self.oneD_vel_model_z_df - LUTs_dict['extent_x_m'] = self.extent_x_m - LUTs_dict['dxz'] = self.dxz - LUTs_dict['array_centre_xz'] = self.array_centre_xz - self.LUTs_fname = os.path.join(LUT_outdir, 'LUTs.pkl') - pickle.dump( LUTs_dict, open( self.LUTs_fname, "wb" ) ) + LUTs_dict["trav_times_grid_P"] = trav_times_grid_P + LUTs_dict["trav_times_grid_S"] = trav_times_grid_S + LUTs_dict["theta_grid_P"] = theta_grid_P + LUTs_dict["theta_grid_S"] = theta_grid_S + LUTs_dict["vel_model_x_labels"] = vel_model_x_labels + LUTs_dict["vel_model_z_labels"] = vel_model_z_labels + LUTs_dict["oneD_vel_model_z_df"] = self.oneD_vel_model_z_df + LUTs_dict["extent_x_m"] = self.extent_x_m + LUTs_dict["dxz"] = self.dxz + LUTs_dict["array_centre_xz"] = self.array_centre_xz + self.LUTs_fname = os.path.join(LUT_outdir, "LUTs.pkl") + pickle.dump(LUTs_dict, open(self.LUTs_fname, "wb")) print("Saved LUT to:", self.LUTs_fname) - self.LUTs_dict = LUTs_dict + self.LUTs_dict = LUTs_dict return LUTs_dict def load_location_LUTs(self, LUTs_fname=None): @@ -1263,23 +1833,22 @@ def load_location_LUTs(self, LUTs_fname=None): Parameters ---------- LUTs_fname : str - Path to LUT file to load. Optional. Default is to use the attribute - . + Path to LUT file to load. Optional. Default is to use the attribute + . """ # Assign attribute, if not already specified correctly: if LUTs_fname: self.LUTs_fname = LUTs_fname # And load LUTs: - LUTs_dict = pickle.load( open( self.LUTs_fname, "rb" ) ) + LUTs_dict = pickle.load(open(self.LUTs_fname, "rb")) # And assign relevent attributes: - self.oneD_vel_model_z_df = LUTs_dict['oneD_vel_model_z_df'] - self.extent_x_m = LUTs_dict['extent_x_m'] - self.dxz = LUTs_dict['dxz'] - self.array_centre_xz = LUTs_dict['array_centre_xz'] - self.LUTs_dict = LUTs_dict + self.oneD_vel_model_z_df = LUTs_dict["oneD_vel_model_z_df"] + self.extent_x_m = LUTs_dict["extent_x_m"] + self.dxz = LUTs_dict["dxz"] + self.array_centre_xz = LUTs_dict["array_centre_xz"] + self.LUTs_dict = LUTs_dict return LUTs_dict - def locate_events(self, events_df, verbosity=0): """Function to locate events using LUT.""" # Perform tests to check that various required attributes are specified: @@ -1299,44 +1868,58 @@ def locate_events(self, events_df, verbosity=0): # Locate events: if not exit_bool: - events_df = location.locate_events_from_P_and_S_array_arrivals(events_df, self.LUTs_dict, self.array_latlon, - self.receiver_vp, self.receiver_vs, - verbosity=verbosity) - + events_df = location.locate_events_from_P_and_S_array_arrivals( + events_df, + self.LUTs_dict, + self.array_latlon, + self.receiver_vp, + self.receiver_vs, + verbosity=verbosity, + ) + return events_df - def save(self, out_fname=None): """Function to save class object to file. Parameters ---------- out_fname : str - Path to save class object to. Optional. If not specified, then save to + Path to save class object to. Optional. If not specified, then save to /detect_obj.pkl. """ # Save class to file: if not out_fname: out_fname = os.path.join(self.outdir, "detect_obj.pkl") - f = open(out_fname, 'wb') + f = open(out_fname, "wb") pickle.dump(self.__dict__, f) f.close() print("Saved detection instance to:", out_fname) - + def load(self, preload_fname): """try load self.name.txt""" - f = open(preload_fname, 'rb') - self.__dict__ = pickle.load(f) + f = open(preload_fname, "rb") + self.__dict__ = pickle.load(f) f.close() print("Loaded detection instance from:", preload_fname) - - def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slows_1_2, t_before_s=10, t_after_s=10, st_out_fname='out.m', return_streams=False): - """Function to find array stacked stream from back-azimuth and slowness. Returns average amplitude + def get_composite_array_st_from_bazi_slowness( + self, + arrival_time, + bazis_1_2, + slows_1_2, + t_before_s=10, + t_after_s=10, + st_out_fname="out.m", + return_streams=False, + method="linear", + degree=1, + ): + """Function to find array stacked stream from back-azimuth and slowness. Returns average amplitude time-series seismogram of stacked array data, for all three componets. Parameters ---------- arrival_time : obspy UTCDateTime object - Arrival time at array. Will output window around this time, as defined by + Arrival time at array. Will output window around this time, as defined by and . bazis_1_2 : list of 2 floats @@ -1346,26 +1929,32 @@ def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slo Slowness of event phase arrivals in vertical and horizontal, in seconds/km. t_before_s : float - Time, in seconds, before arrival time to include in window. Optional. Default + Time, in seconds, before arrival time to include in window. Optional. Default is 10 s. t_after_s : float - Time, in seconds, after arrival time to include in window. Optional. Default - is 10 s. - + Time, in seconds, after arrival time to include in window. Optional. Default + is 10 s. + st_out_fname : str Filename of to save mseed data stream to. Optional. Default is out.m. return_streams : bool - If True, returns st and composite_st. Optional. Default = False. + If True, returns st and composite_st. Optional. Default = False. """ # Load in raw mseed data: - st = self._load_day_of_data(arrival_time.year, arrival_time.julday, hour=arrival_time.hour) + st = self._load_data( + arrival_time.year, + arrival_time.month, + arrival_time.day, + hour=arrival_time.hour, + norm=False, + ) # And trim data: - st.trim(starttime=arrival_time-t_before_s, endtime=arrival_time+t_after_s) + st.trim(starttime=arrival_time - t_before_s, endtime=arrival_time + t_after_s) # Upsample data soas to provide best time shift later: - st.interpolate(sampling_rate=10*st[0].stats.sampling_rate) + st.interpolate(sampling_rate=10 * st[0].stats.sampling_rate) # Find and perform time shifts for all receivers: # Loop over stations: @@ -1373,31 +1962,35 @@ def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slo # Find time shift: # Get current station location (relative to array centre): curr_station = st[i].stats.station - x_rec = self.stations_df.loc[self.stations_df['Name'] == curr_station]['x_array_coords_km'].values[0] - y_rec = self.stations_df.loc[self.stations_df['Name'] == curr_station]['y_array_coords_km'].values[0] + x_rec = self.stations_df.loc[self.stations_df["Name"] == curr_station][ + "x_array_coords_km" + ].values[0] + y_rec = self.stations_df.loc[self.stations_df["Name"] == curr_station][ + "y_array_coords_km" + ].values[0] # Select either vertical or horizontal slowness and back-azimuth, depending on component: comp = st[i].stats.channel[-1] - if comp == 'Z': + if comp == "Z": bazi = bazis_1_2[0] slow = slows_1_2[0] - elif comp == 'N' or comp == 'E' or comp == '1' or comp == '2': + elif comp == "N" or comp == "E" or comp == "1" or comp == "2": bazi = bazis_1_2[1] slow = slows_1_2[1] # Calculate arrival time shift relative to centre of array: - time_shift_curr_s = _calc_time_shift_from_array_cent(slow, bazi, x_rec, y_rec) - + time_shift_curr_s = _calc_time_shift_from_array_cent( + slow, bazi, x_rec, y_rec + ) + # And perform time shift on data: n_samp_to_shift = round(time_shift_curr_s * st[i].stats.sampling_rate) st[i].data = np.roll(st[i].data, n_samp_to_shift) - + # And find stacked, mean and stdev of data: n_stat = len(st.select(channel="??Z")) # Get unique channels: chan_labels_tmp = [] - chan_labels_unique = [] for tr in st: chan_labels_tmp.append(tr.stats.channel) - chan_labels_unique = list(set(chan_labels_tmp)) # Create datastores to save to: max_st_len = 0 for i in range(len(st)): @@ -1411,36 +2004,55 @@ def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slo for i in range(n_stat): # Z: if len(st.select(channel="??Z")[i].data) == max_st_len: - Z_all[:,i] = st.select(channel="??Z")[i].data + Z_all[:, i] = st.select(channel="??Z")[i].data else: - Z_all[:len(st.select(channel="??Z")[i].data),i] = st.select(channel="??Z")[i].data - try: + Z_all[: len(st.select(channel="??Z")[i].data), i] = st.select( + channel="??Z" + )[i].data + try: # N: if len(st.select(channel="??N")[i].data) == max_st_len: - N_all[:,i] = st.select(channel="??N")[i].data + N_all[:, i] = st.select(channel="??N")[i].data else: - N_all[:len(st.select(channel="??N")[i].data),i] = st.select(channel="??N")[i].data + N_all[: len(st.select(channel="??N")[i].data), i] = st.select( + channel="??N" + )[i].data # E: if len(st.select(channel="??E")[i].data) == max_st_len: - E_all[:,i] = st.select(channel="??E")[i].data + E_all[:, i] = st.select(channel="??E")[i].data else: - E_all[:len(st.select(channel="??E")[i].data),i] = st.select(channel="??E")[i].data + E_all[: len(st.select(channel="??E")[i].data), i] = st.select( + channel="??E" + )[i].data except IndexError: # And write if uses 1 and 2 labels rather than N and E: # N: if len(st.select(channel="??1")[i].data) == max_st_len: - N_all[:,i] = st.select(channel="??1")[i].data + N_all[:, i] = st.select(channel="??1")[i].data else: - N_all[:len(st.select(channel="??1")[i].data),i] = st.select(channel="??1")[i].data + N_all[: len(st.select(channel="??1")[i].data), i] = st.select( + channel="??1" + )[i].data # E: if len(st.select(channel="??2")[i].data) == max_st_len: - E_all[:,i] = st.select(channel="??2")[i].data + E_all[:, i] = st.select(channel="??2")[i].data else: - E_all[:len(st.select(channel="??2")[i].data),i] = st.select(channel="??2")[i].data - + E_all[: len(st.select(channel="??2")[i].data), i] = st.select( + channel="??2" + )[i].data + except: + print("Date gap, continue... for now") + continue # And create stacked data stream: - composite_st = _create_stacked_data_st(st, Z_all, N_all, E_all) - + if method == "linear": + composite_st = _create_stacked_data_st(st, Z_all, N_all, E_all) + elif method == "pws": + composite_st = _create_phase_weighted_stack_st( + st, Z_all, N_all, E_all, degree + ) + elif method == "nth_root": + pass + # composite_st = _create_nth_root_stack_st(st, Z_all, N_all, E_all, degree) # And decimate data back down to original sampling rate: st.decimate(10, no_filter=True) composite_st.decimate(10, no_filter=True) @@ -1448,7 +2060,7 @@ def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slo # And save data: if st_out_fname: st.write(st_out_fname, format="MSEED") - composite_st_out_fname = st_out_fname.split('.')[0] + "_composite.m" + composite_st_out_fname = st_out_fname.split(".")[0] + "_composite.m" composite_st.write(composite_st_out_fname, format="MSEED") if return_streams: @@ -1458,18 +2070,45 @@ def get_composite_array_st_from_bazi_slowness(self, arrival_time, bazis_1_2, slo gc.collect() +def moving_window_mad(trace, window_len, mad_multiplier): + """ + Function to calculate the median absolute deviation (MAD) using a + moving window over a beam power time-series. + For the start and end of the time-series, where a full window is not + available, the window is truncated to fit within the time-series. + Note that the MAD is scaled to be equivalent to the standard deviation + for a Gaussian distribution, by multiplying by the constant 1.4826. + See: https://en.wikipedia.org/wiki/Median_absolute_deviation + Parameters + ---------- + trace : np.ndarray + 1D numpy array containing beam power time-series. + window_len : int + Length of moving window, in number of samples. - - - - - - -#----------------------------------------------- End: Define main functions ----------------------------------------------- - - - + mad_multiplier : float + Multiplier for MAD to set detection threshold. For example, a value of + 2 would set the threshold to be 2 times the MAD above the median. This is scaled + by 1.4826 within the function to be equivalent to standard deviations for a + Gaussian distribution. + """ + # Create datastore: + mad_thresholds = np.zeros(trace.shape) + half_win = int(window_len / 2) + for i in range(len(trace)): + start = max(0, i - half_win) + end = min(len(trace), i + half_win) + window = trace[start:end] + mad = np.median(np.abs(window - np.median(window))) + # add one to mad_multiplier as the multiplier represents the + # number of standard deviations above the median. + mad_thresholds[i] = MAD_SCALE * (1 + mad_multiplier) * mad + + return mad_thresholds + + +# ------------- End: Define main functions -------------------- diff --git a/SeisSeeker/processing/location.py b/SeisSeeker/processing/location.py index 1fe7894..5126083 100755 --- a/SeisSeeker/processing/location.py +++ b/SeisSeeker/processing/location.py @@ -1,35 +1,36 @@ #!/usr/bin/python -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Script Description: # Script to perform earthquake location using array processing methods. # Created by Tom Hudson, 10th August 2022 -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Import neccessary modules: import pandas as pd import numpy as np -import matplotlib.pyplot as plt +import matplotlib.pyplot as plt import matplotlib import os, sys import obspy from scipy.signal import find_peaks from numba import jit, objmode, prange, set_num_threads -import gc -# import multiprocessing as mp -import time +import gc +# import multiprocessing as mp +import time -#----------------------------------------------- Define main functions ----------------------------------------------- +# ----------------------------------------------- Define main functions ----------------------------------------------- class CustomError(Exception): pass - -def locate_events_from_P_and_S_array_arrivals(events_df, LUTs_dict, array_latlon, receiver_vp, receiver_vs, verbosity=0): +def locate_events_from_P_and_S_array_arrivals( + events_df, LUTs_dict, array_latlon, receiver_vp, receiver_vs, verbosity=0 +): """Function to locate events from P and S phase arrivals.""" # Append rows to events_df to save locations: events_df["x_km"] = "" @@ -39,62 +40,66 @@ def locate_events_from_P_and_S_array_arrivals(events_df, LUTs_dict, array_latlon events_df["lon"] = "" # Calculate various objects only once, for efficiency: - delta_tp_ts_grid = LUTs_dict['trav_times_grid_S'] - LUTs_dict['trav_times_grid_P'] + delta_tp_ts_grid = LUTs_dict["trav_times_grid_S"] - LUTs_dict["trav_times_grid_P"] # Loop over events, processing to find locations: - count=0 + count = 0 for index, row in events_df.iterrows(): - print("Locating event", count+1, '/', len(events_df['bazi1'])) + print("Locating event", count + 1, "/", len(events_df["bazi1"])) # 1. Find effective radius within LUT grid by minimising: # delta_tp_ts - (ts_tt - tp_tt) # To create PDF of uncertainty in grids of LUT for location # Calc. delta_tp_ts: - delta_tp_ts = obspy.UTCDateTime(row['t2']) - obspy.UTCDateTime(row['t1']) + delta_tp_ts = obspy.UTCDateTime(row["t2"]) - obspy.UTCDateTime(row["t1"]) # Find minimum: abs_min_arr = np.abs(delta_tp_ts - delta_tp_ts_grid) tp_ts_res_pdf = 1 - (abs_min_arr / np.max(abs_min_arr)) # 2. Use inclination angle from slowness (of P and/or S) to search from LUT for possible cells: - if row['slow1'] > 0 and row['slow2'] > 0: + if row["slow1"] > 0 and row["slow2"] > 0: # To create PDF for location of cells vertically in LUT # 2.i. For P: # Calculate inclination angle: # v = v_app * sin(inc_angle_from_vert) ?! # Therefore, theta = arcsin( v / v_app ) ?! - v_app_P = 1. / row['slow1'] - inc_angle_P = np.rad2deg(np.arcsin( v_app_P / receiver_vp )) + v_app_P = 1.0 / row["slow1"] + inc_angle_P = np.rad2deg(np.arcsin(v_app_P / receiver_vp)) # And calculate inc. angle pdf: - abs_min_arr = np.abs(LUTs_dict['theta_grid_P'] - inc_angle_P) + abs_min_arr = np.abs(LUTs_dict["theta_grid_P"] - inc_angle_P) P_inc_angle_res_pdf = 1 - (abs_min_arr / np.max(abs_min_arr)) # 2.ii. For S: # Calculate inclination angle: # v = v_app * sin(inc_angle_from_vert) ?! # Therefore, theta = arcsin( v / v_app ) ?! - v_app_S = 1. / row['slow2'] - inc_angle_S = np.rad2deg(np.arcsin( v_app_S / receiver_vs )) + v_app_S = 1.0 / row["slow2"] + inc_angle_S = np.rad2deg(np.arcsin(v_app_S / receiver_vs)) # And calculate inc. angle pdf: - abs_min_arr = np.abs(LUTs_dict['theta_grid_S'] - inc_angle_S) + abs_min_arr = np.abs(LUTs_dict["theta_grid_S"] - inc_angle_S) S_inc_angle_res_pdf = 1 - (abs_min_arr / np.max(abs_min_arr)) # 2.iii. Stack P and S inc. angle PDFs: - PS_stack_inc_angle_res_pdf = (P_inc_angle_res_pdf + S_inc_angle_res_pdf) / 2. + PS_stack_inc_angle_res_pdf = ( + P_inc_angle_res_pdf + S_inc_angle_res_pdf + ) / 2.0 # 3. Combine radius PDF and inc-angle PDF to get best result location within 2D LUT: # Stack pdfs: - stacked_pdf = (tp_ts_res_pdf + PS_stack_inc_angle_res_pdf) / 2. + stacked_pdf = (tp_ts_res_pdf + PS_stack_inc_angle_res_pdf) / 2.0 # And get best result: max_xz_idxs = np.argwhere(stacked_pdf == np.max(stacked_pdf)) x_idx_curr = max_xz_idxs[0][0] z_idx_curr = max_xz_idxs[0][1] - event_x_coord_km = LUTs_dict['vel_model_x_labels'][x_idx_curr] - event_z_coord_km = LUTs_dict['vel_model_z_labels'][z_idx_curr] + event_x_coord_km = LUTs_dict["vel_model_x_labels"][x_idx_curr] + event_z_coord_km = LUTs_dict["vel_model_z_labels"][z_idx_curr] else: event_x_coord_km = np.nan event_z_coord_km = np.nan # Plot workings, if specified: if verbosity > 1: - fig, axes = plt.subplots(nrows=4, figsize=(4,16)) - Z, X = np.meshgrid(LUTs_dict['vel_model_z_labels'], LUTs_dict['vel_model_x_labels']) + fig, axes = plt.subplots(nrows=4, figsize=(4, 16)) + Z, X = np.meshgrid( + LUTs_dict["vel_model_z_labels"], LUTs_dict["vel_model_x_labels"] + ) # Plot P-S pdf: im = axes[0].pcolormesh(X, Z, tp_ts_res_pdf, cmap="Greys_r", vmin=0, vmax=1) plt.colorbar(im, ax=axes[0], label="PDF, $t_{P-S}$") @@ -103,16 +108,20 @@ def locate_events_from_P_and_S_array_arrivals(events_df, LUTs_dict, array_latlon # z_idx_curr = min_xz_idxs[i][1] # axes[0].scatter(LUTs_dict['vel_model_x_labels'][x_idx_curr], LUTs_dict['vel_model_z_labels'][z_idx_curr], c='k') # And plot P inc. angle pdf: - im2 = axes[1].pcolormesh(X, Z, P_inc_angle_res_pdf, cmap="Greys_r", vmin=0, vmax=1) + im2 = axes[1].pcolormesh( + X, Z, P_inc_angle_res_pdf, cmap="Greys_r", vmin=0, vmax=1 + ) plt.colorbar(im2, ax=axes[1], label="PDF, $\\theta_P$") # And plot S inc. angle pdf: - im3 = axes[2].pcolormesh(X, Z, S_inc_angle_res_pdf, cmap="Greys_r", vmin=0, vmax=1) + im3 = axes[2].pcolormesh( + X, Z, S_inc_angle_res_pdf, cmap="Greys_r", vmin=0, vmax=1 + ) plt.colorbar(im3, ax=axes[2], label="PDF, $\\theta_S$") # And plot stacked pdf: im4 = axes[3].pcolormesh(X, Z, stacked_pdf, cmap="Greys_r", vmin=0, vmax=1) plt.colorbar(im4, ax=axes[3], label="stacked PDF") # And plot minimum idx on stacked pdf: - axes[3].scatter(event_x_coord_km, event_z_coord_km, c='g') + axes[3].scatter(event_x_coord_km, event_z_coord_km, c="g") for i in range(len(axes)): axes[i].invert_yaxis() axes[i].set_xlabel("X (m)") @@ -122,28 +131,32 @@ def locate_events_from_P_and_S_array_arrivals(events_df, LUTs_dict, array_latlon # 4. Convert 2D LUT result into 3D cartesian location by combining with bazi.: # (and append to events_df) r_hor_km = event_x_coord_km / 1000 - mean_bazi = (row['bazi1'] + row['bazi2']) / 2 - events_df.iloc[count, events_df.columns.get_loc('x_km')] = r_hor_km * np.sin( np.deg2rad( mean_bazi ) ) - events_df.iloc[count, events_df.columns.get_loc('y_km')] = r_hor_km * np.cos( np.deg2rad( mean_bazi ) ) - events_df.iloc[count, events_df.columns.get_loc('z_km')] = event_z_coord_km / 1000 + mean_bazi = (row["bazi1"] + row["bazi2"]) / 2 + events_df.iloc[count, events_df.columns.get_loc("x_km")] = r_hor_km * np.sin( + np.deg2rad(mean_bazi) + ) + events_df.iloc[count, events_df.columns.get_loc("y_km")] = r_hor_km * np.cos( + np.deg2rad(mean_bazi) + ) + events_df.iloc[count, events_df.columns.get_loc("z_km")] = ( + event_z_coord_km / 1000 + ) # Calculate event lat and lons relative to array centre: - events_df.iloc[count, events_df.columns.get_loc('lat')] = array_latlon[0] + obspy.geodetics.base.kilometers2degrees(events_df.iloc[count, events_df.columns.get_loc('y_km')]) - events_df.iloc[count, events_df.columns.get_loc('lon')] = array_latlon[1] + obspy.geodetics.base.kilometers2degrees(events_df.iloc[count, events_df.columns.get_loc('x_km')]) + events_df.iloc[count, events_df.columns.get_loc("lat")] = array_latlon[ + 0 + ] + obspy.geodetics.base.kilometers2degrees( + events_df.iloc[count, events_df.columns.get_loc("y_km")] + ) + events_df.iloc[count, events_df.columns.get_loc("lon")] = array_latlon[ + 1 + ] + obspy.geodetics.base.kilometers2degrees( + events_df.iloc[count, events_df.columns.get_loc("x_km")] + ) # Update event count: - count+=1 + count += 1 return events_df - - - - - - - -#----------------------------------------------- End: Define main functions ----------------------------------------------- - - - +# ----------------------------------------------- End: Define main functions ----------------------------------------------- diff --git a/SeisSeeker/processing/lookup_table_manager.py b/SeisSeeker/processing/lookup_table_manager.py index f283f11..ce967ef 100755 --- a/SeisSeeker/processing/lookup_table_manager.py +++ b/SeisSeeker/processing/lookup_table_manager.py @@ -1,5 +1,5 @@ #!/Users/eart0504/opt/anaconda3/bin/python -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Script Description: # Script to calculate travel times lookup tables for various seismic phases. @@ -10,24 +10,26 @@ # Created by Tom Hudson, 17th August 2022 -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Import neccessary modules: import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -import skfmm # For fast-marching travel-time lookup tables +import skfmm # For fast-marching travel-time lookup tables + # import ttcrpy.rgrid as ttcrpy_rgrid # For ray-tracing based incidence angles -import gc -import pykonal # For ray-tracing based incidence angles +import gc + +# import pykonal # For ray-tracing based incidence angles -#----------------------------------------------- Define constants and parameters ----------------------------------------------- +# ----------------------------------------------- Define constants and parameters ----------------------------------------------- -#----------------------------------------------- End: Define constants and parameters ----------------------------------------------- +# ----------------------------------------------- End: Define constants and parameters ----------------------------------------------- -#----------------------------------------------- Define main functions ----------------------------------------------- -def read_1D_vel_model_to_3D_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000,3000], dxyz=[1.,1.,1.]): +# ----------------------------------------------- Define main functions ----------------------------------------------- +def read_1D_vel_model_to_3D_model_for_fmm( + oneD_vel_model_z_df, extent_xy_m=[3000, 3000], dxyz=[1.0, 1.0, 1.0] +): """Function to create specific velocity model. Inputs: oneD_vel_model_z_df - Pandas DataFrame of shape (z_extent, 3) with the columns corresponding to: depth, vp, vs @@ -37,12 +39,12 @@ def read_1D_vel_model_to_3D_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000 vel_model_arr_S - 3D Velocity grid for S wave ix, iy, iz - Grids containing index labels """ - dx = dxyz[0] # Grid spacing in x dir - dy = dxyz[1] # Grid spacing in x dir - dz = dxyz[2] # Grid spacing in z dir - vel_model_x_labels = np.arange(0., extent_xy_m[0]+dx, dx) - vel_model_y_labels = np.arange(0., extent_xy_m[1]+dy, dy) - vel_model_z_labels = oneD_vel_model_z_df['depth'].values + dx = dxyz[0] # Grid spacing in x dir + dy = dxyz[1] # Grid spacing in x dir + dz = dxyz[2] # Grid spacing in z dir + vel_model_x_labels = np.arange(0.0, extent_xy_m[0] + dx, dx) + vel_model_y_labels = np.arange(0.0, extent_xy_m[1] + dy, dy) + vel_model_z_labels = oneD_vel_model_z_df["depth"].values # Create the index grids: iy, ix, iz = np.meshgrid(vel_model_y_labels, vel_model_x_labels, vel_model_z_labels) # Create the P and S wave model: @@ -50,9 +52,15 @@ def read_1D_vel_model_to_3D_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000 vel_model_arr_S = np.zeros(np.shape(ix), dtype=float) for i in range(vel_model_arr_P.shape[0]): for j in range(vel_model_arr_P.shape[1]): - vel_model_arr_P[i,j,:] = oneD_vel_model_z_df['vp'] - vel_model_arr_S[i,j,:] = oneD_vel_model_z_df['vs'] - return vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_y_labels, vel_model_z_labels + vel_model_arr_P[i, j, :] = oneD_vel_model_z_df["vp"] + vel_model_arr_S[i, j, :] = oneD_vel_model_z_df["vs"] + return ( + vel_model_arr_P, + vel_model_arr_S, + vel_model_x_labels, + vel_model_y_labels, + vel_model_z_labels, + ) def vec_norm(x): @@ -60,7 +68,9 @@ def vec_norm(x): return np.sqrt(x.dot(x)) -def read_1D_vel_model_to_2D_model_for_fmm(oneD_vel_model_z_df, extent_x_m=3000, dxz=[1.,1.]): +def read_1D_vel_model_to_2D_model_for_fmm( + oneD_vel_model_z_df, extent_x_m=3000, dxz=[1.0, 1.0] +): """Function to create specific velocity model in 2D (x, z). Inputs: oneD_vel_model_z_df - Pandas DataFrame of shape (z_extent, 3) with the columns corresponding to: depth, vp, vs @@ -69,26 +79,26 @@ def read_1D_vel_model_to_2D_model_for_fmm(oneD_vel_model_z_df, extent_x_m=3000, vel_model_arr_P - 3D Velocity grid for P wave vel_model_arr_S - 3D Velocity grid for S wave """ - dx = dxz[0] # Grid spacing in x dir - dz = dxz[1] # Grid spacing in z dir - vel_model_x_labels = np.arange(0., extent_x_m+dx, dx) - vel_model_z_labels = oneD_vel_model_z_df['depth'].values + dx = dxz[0] # Grid spacing in x dir + dz = dxz[1] # Grid spacing in z dir + vel_model_x_labels = np.arange(0.0, extent_x_m + dx, dx) + vel_model_z_labels = oneD_vel_model_z_df["depth"].values # Create the index grids: iz, ix = np.meshgrid(vel_model_z_labels, vel_model_x_labels) # Create the P and S wave model: vel_model_arr_P = np.zeros(np.shape(ix), dtype=float) vel_model_arr_S = np.zeros(np.shape(iz), dtype=float) for i in range(vel_model_arr_P.shape[0]): - vel_model_arr_P[i,:] = oneD_vel_model_z_df['vp'] - vel_model_arr_S[i,:] = oneD_vel_model_z_df['vs'] + vel_model_arr_P[i, :] = oneD_vel_model_z_df["vp"] + vel_model_arr_S[i, :] = oneD_vel_model_z_df["vs"] return vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_z_labels -def ray_tracer(v_model, dxyz, src_loc=[0,0,0], rx_loc=[0,0,0]): +def ray_tracer(v_model, dxyz, src_loc=[0, 0, 0], rx_loc=[0, 0, 0]): """Function to perform ray-tracing, using pykonal package. All units are in km. src_loc, rx_loc are in km x,y,z.""" # Define the solver: - solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian") + solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian") # Define the computational domain: solver.vv.min_coords = 0, 0, 0 solver.vv.node_intervals = dxyz[0], dxyz[1], dxyz[2] @@ -111,19 +121,24 @@ def find_nearest(array, value): idx = (np.abs(array - value)).argmin() return array[idx], idx -def create_2D_LUT(oneD_vel_model_z_df, array_centre_xz, extent_x_m=3000., dxz=[1.,1.], n_threads=1): - """Get theoretical P and S travel times for all nodes in a lookup table, relative to the centre of the array + +def create_2D_LUT( + oneD_vel_model_z_df, array_centre_xz, extent_x_m=3000.0, dxz=[1.0, 1.0], n_threads=1 +): + """Get theoretical P and S travel times for all nodes in a lookup table, relative to the centre of the array (using Eikonal method). Inputs: oneD_vel_model_z_df - Pandas DataFrame of shape (z_extent, 3) with the columns corresponding to: depth, vp, vs. array_centre_xz - The location of the array centre within the proposed lookup table grid, in metres. """ # Get 2D velocity model: - dx = dxz[0] # Grid spacing in x dir - dz = dxz[1] # Grid spacing in z dir - vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_z_labels = read_1D_vel_model_to_2D_model_for_fmm(oneD_vel_model_z_df, - extent_x_m=extent_x_m, - dxz=dxz) + dx = dxz[0] # Grid spacing in x dir + dz = dxz[1] # Grid spacing in z dir + vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_z_labels = ( + read_1D_vel_model_to_2D_model_for_fmm( + oneD_vel_model_z_df, extent_x_m=extent_x_m, dxz=dxz + ) + ) # Create travel-time lookup table grids: # (from fast marching method) @@ -133,10 +148,10 @@ def create_2D_LUT(oneD_vel_model_z_df, array_centre_xz, extent_x_m=3000., dxz=[1 val, x_idx = find_nearest(vel_model_x_labels, array_centre_xz[0]) val, z_idx = find_nearest(vel_model_z_labels, array_centre_xz[1]) # Set array_centre_xyz location: - phi[x_idx, z_idx] = 1. + phi[x_idx, z_idx] = 1.0 # Calculate travel times array to all the points in the grid: - trav_times_grid_P = skfmm.travel_time(phi,vel_model_arr_P,dx=dxz) - trav_times_grid_S = skfmm.travel_time(phi,vel_model_arr_S,dx=dxz) + trav_times_grid_P = skfmm.travel_time(phi, vel_model_arr_P, dx=dxz) + trav_times_grid_S = skfmm.travel_time(phi, vel_model_arr_S, dx=dxz) # And create incidence angle lookup table grids: # (from ray tracing) @@ -145,8 +160,12 @@ def create_2D_LUT(oneD_vel_model_z_df, array_centre_xz, extent_x_m=3000., dxz=[1 # z_node_labels = vel_model_z_labels # rgrid = ttcrpy_rgrid.Grid2d(x_node_labels, z_node_labels, cell_slowness=False, n_threads=n_threads) # Specify velocity model for ray tracing: - vel_model_arr_P_3D = vel_model_arr_P.reshape(1, vel_model_arr_P.shape[0], vel_model_arr_P.shape[1]) - vel_model_arr_S_3D = vel_model_arr_S.reshape(1, vel_model_arr_S.shape[0], vel_model_arr_S.shape[1]) + vel_model_arr_P_3D = vel_model_arr_P.reshape( + 1, vel_model_arr_P.shape[0], vel_model_arr_P.shape[1] + ) + vel_model_arr_S_3D = vel_model_arr_S.reshape( + 1, vel_model_arr_S.shape[0], vel_model_arr_S.shape[1] + ) # Calculate rays for each point in the grid: theta_grid_P = np.zeros(trav_times_grid_P.shape) theta_grid_S = np.zeros(trav_times_grid_P.shape) @@ -157,96 +176,99 @@ def create_2D_LUT(oneD_vel_model_z_df, array_centre_xz, extent_x_m=3000., dxz=[1 for j in range(theta_grid_P.shape[1]): # Print progress: if count % 100 == 0: - print("Processing for ray", count, "/", theta_grid_P.shape[0]*theta_grid_P.shape[1]) + print( + "Processing for ray", + count, + "/", + theta_grid_P.shape[0] * theta_grid_P.shape[1], + ) count += 1 # Calculate ray-tracing and incidence angle for P wave: # Calculate rays for current event: # (Note that coords are converted into 3D and km) - node_coords = np.array([0, vel_model_x_labels[i], vel_model_z_labels[j]], dtype=float) / 1000 - node_coords = node_coords + 0.001 # Add 1 metre, so that node coords are non-zero - receiver_coords = np.array([0, array_centre_xz[0], array_centre_xz[1]], dtype=float) / 1000 + node_coords = ( + np.array([0, vel_model_x_labels[i], vel_model_z_labels[j]], dtype=float) + / 1000 + ) + node_coords = ( + node_coords + 0.001 + ) # Add 1 metre, so that node coords are non-zero + receiver_coords = ( + np.array([0, array_centre_xz[0], array_centre_xz[1]], dtype=float) + / 1000 + ) # tt, rays = rgrid.raytrace(event_coords, receiver_coords, 1./vel_model_arr_P, return_rays=True) # Perform ray tracing: dxyz = np.array([dx, dx, dz], dtype=float) / 1000 - ray = ray_tracer(vel_model_arr_P_3D / 1000, dxyz, src_loc=node_coords, rx_loc=receiver_coords) + ray = ray_tracer( + vel_model_arr_P_3D / 1000, + dxyz, + src_loc=node_coords, + rx_loc=receiver_coords, + ) # Calculate incidence angle at array: try: - ray_vec_at_receiver = - np.array([ ray[-1,1] - ray[-2,1], ray[-1,2] - ray[-2,2] ]) # Note minus sign, as defining as vector out from array - vert_vec = np.array([ 0, 1]) - theta_curr = np.arccos( ( ray_vec_at_receiver.dot(vert_vec) ) / ( vec_norm(ray_vec_at_receiver) * vec_norm(vert_vec) ) ) # cos(theta) = a.b / |a| |b| + ray_vec_at_receiver = -np.array( + [ray[-1, 1] - ray[-2, 1], ray[-1, 2] - ray[-2, 2]] + ) # Note minus sign, as defining as vector out from array + vert_vec = np.array([0, 1]) + theta_curr = np.arccos( + (ray_vec_at_receiver.dot(vert_vec)) + / (vec_norm(ray_vec_at_receiver) * vec_norm(vert_vec)) + ) # cos(theta) = a.b / |a| |b| theta_curr = np.rad2deg(theta_curr) - theta_grid_P[i,j] = theta_curr + theta_grid_P[i, j] = theta_curr except IndexError: # Or assign previous value, if failed to trace rays for some reason: - if j>0: - theta_grid_P[i,j] = theta_grid_P[i,j-1] + if j > 0: + theta_grid_P[i, j] = theta_grid_P[i, j - 1] else: - theta_grid_P[i,j] = theta_grid_P[i-1,j] + theta_grid_P[i, j] = theta_grid_P[i - 1, j] # And clear up: - del ray + del ray gc.collect() # Calculate ray-tracing and incidence angle for S wave: # Calculate rays for current event: # (Uses some parameters as for P wave) # Perform ray tracing: - ray = ray_tracer(vel_model_arr_S_3D / 1000, dxyz, src_loc=node_coords, rx_loc=receiver_coords) + ray = ray_tracer( + vel_model_arr_S_3D / 1000, + dxyz, + src_loc=node_coords, + rx_loc=receiver_coords, + ) # Calculate incidence angle at array: try: - ray_vec_at_receiver = - np.array([ ray[-1,1] - ray[-2,1], ray[-1,2] - ray[-2,2] ]) # Note minus sign, as defining as vector out from array - vert_vec = np.array([ 0, 1]) - theta_curr = np.arccos( ( ray_vec_at_receiver.dot(vert_vec) ) / ( vec_norm(ray_vec_at_receiver) * vec_norm(vert_vec) ) ) # cos(theta) = a.b / |a| |b| + ray_vec_at_receiver = -np.array( + [ray[-1, 1] - ray[-2, 1], ray[-1, 2] - ray[-2, 2]] + ) # Note minus sign, as defining as vector out from array + vert_vec = np.array([0, 1]) + theta_curr = np.arccos( + (ray_vec_at_receiver.dot(vert_vec)) + / (vec_norm(ray_vec_at_receiver) * vec_norm(vert_vec)) + ) # cos(theta) = a.b / |a| |b| theta_curr = np.rad2deg(theta_curr) - theta_grid_S[i,j] = theta_curr + theta_grid_S[i, j] = theta_curr except IndexError: # Or assign previous value, if failed to trace rays for some reason: - if j>0: - theta_grid_S[i,j] = theta_grid_S[i,j-1] + if j > 0: + theta_grid_S[i, j] = theta_grid_S[i, j - 1] else: - theta_grid_S[i,j] = theta_grid_S[i-1,j] + theta_grid_S[i, j] = theta_grid_S[i - 1, j] # And clear up: - del ray + del ray gc.collect() - - - - return trav_times_grid_P, trav_times_grid_S, theta_grid_P, theta_grid_S, vel_model_x_labels, vel_model_z_labels - - - - - - -#----------------------------------------------- End: Define main functions ----------------------------------------------- - - - - - - - - - - - - - - - - - - - - - - - - - - - - + return ( + trav_times_grid_P, + trav_times_grid_S, + theta_grid_P, + theta_grid_S, + vel_model_x_labels, + vel_model_z_labels, + ) +# ----------------------------------------------- End: Define main functions ----------------------------------------------- diff --git a/SeisSeeker/processing/lookup_table_manager_3D_backup.py b/SeisSeeker/processing/lookup_table_manager_3D_backup.py index d199e18..c7d0abc 100755 --- a/SeisSeeker/processing/lookup_table_manager_3D_backup.py +++ b/SeisSeeker/processing/lookup_table_manager_3D_backup.py @@ -1,5 +1,5 @@ #!/Users/eart0504/opt/anaconda3/bin/python -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Script Description: # Script to calculate travel times lookup tables for various seismic phases. @@ -10,21 +10,23 @@ # Created by Tom Hudson, 17th August 2022 -#----------------------------------------------------------------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------------------------------------------------------------------- # Import neccessary modules: import numpy as np -import pandas as pd +import pandas as pd import matplotlib.pyplot as plt -import skfmm +import skfmm -#----------------------------------------------- Define constants and parameters ----------------------------------------------- +# ----------------------------------------------- Define constants and parameters ----------------------------------------------- -#----------------------------------------------- End: Define constants and parameters ----------------------------------------------- +# ----------------------------------------------- End: Define constants and parameters ----------------------------------------------- -#----------------------------------------------- Define main functions ----------------------------------------------- -def read_1D_vel_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000,3000], dxyz=[1.,1.,1.]): +# ----------------------------------------------- Define main functions ----------------------------------------------- +def read_1D_vel_model_for_fmm( + oneD_vel_model_z_df, extent_xy_m=[3000, 3000], dxyz=[1.0, 1.0, 1.0] +): """Function to create specific ice velocity model. Inputs: oneD_vel_model_z_df - Pandas DataFrame of shape (z_extent, 3) with the columns corresponding to: depth, vp, vs @@ -34,12 +36,12 @@ def read_1D_vel_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000,3000], dxyz vel_model_arr_S - 3D Velocity grid for S wave ix, iy, iz - Grids containing index labels """ - dx = dxyz[0] # Grid spacing in x dir - dy = dxyz[1] # Grid spacing in x dir - dz = dxyz[2] # Grid spacing in z dir - vel_model_x_labels = np.arange(0., extent_xy_m[0]+dx, dx) - vel_model_y_labels = np.arange(0., extent_xy_m[1]+dy, dy) - vel_model_z_labels = oneD_vel_model_z_df['depth'].values + dx = dxyz[0] # Grid spacing in x dir + dy = dxyz[1] # Grid spacing in x dir + dz = dxyz[2] # Grid spacing in z dir + vel_model_x_labels = np.arange(0.0, extent_xy_m[0] + dx, dx) + vel_model_y_labels = np.arange(0.0, extent_xy_m[1] + dy, dy) + vel_model_z_labels = oneD_vel_model_z_df["depth"].values # Create the index grids: iy, ix, iz = np.meshgrid(vel_model_y_labels, vel_model_x_labels, vel_model_z_labels) # Create the P and S wave model: @@ -47,17 +49,30 @@ def read_1D_vel_model_for_fmm(oneD_vel_model_z_df, extent_xy_m=[3000,3000], dxyz vel_model_arr_S = np.zeros(np.shape(ix), dtype=float) for i in range(vel_model_arr_P.shape[0]): for j in range(vel_model_arr_P.shape[1]): - vel_model_arr_P[i,j,:] = oneD_vel_model_z_df['vp'] - vel_model_arr_S[i,j,:] = oneD_vel_model_z_df['vs'] - return vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_y_labels, vel_model_z_labels + vel_model_arr_P[i, j, :] = oneD_vel_model_z_df["vp"] + vel_model_arr_S[i, j, :] = oneD_vel_model_z_df["vs"] + return ( + vel_model_arr_P, + vel_model_arr_S, + vel_model_x_labels, + vel_model_y_labels, + vel_model_z_labels, + ) + def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return array[idx], idx -def create_LUT(oneD_vel_model_z_df, array_centre_xyz, extent_xy_m=[3000.,3000.], dxyz=[1.,1.,1.]): - """Get theoretical P and S travel times for all nodes in a lookup table, relative to the centre of the array + +def create_LUT( + oneD_vel_model_z_df, + array_centre_xyz, + extent_xy_m=[3000.0, 3000.0], + dxyz=[1.0, 1.0, 1.0], +): + """Get theoretical P and S travel times for all nodes in a lookup table, relative to the centre of the array (using Eikonal method). Inputs: oneD_vel_model_z_df - Pandas DataFrame of shape (z_extent, 3) with the columns corresponding to: depth, vp, vs. @@ -65,13 +80,19 @@ def create_LUT(oneD_vel_model_z_df, array_centre_xyz, extent_xy_m=[3000.,3000.], src_depth_xyz - The source depth in x,y,z in metres. """ # Get velocity model: - dx = dxyz[0] # Grid spacing in x dir - dy = dxyz[1] # Grid spacing in x dir - dz = dxyz[2] # Grid spacing in z dir - vel_model_arr_P, vel_model_arr_S, vel_model_x_labels, vel_model_y_labels, vel_model_z_labels = read_1D_vel_model_for_fmm(oneD_vel_model_z_df, - extent_xy_m=extent_xy_m, - dxyz=dxyz) - + dx = dxyz[0] # Grid spacing in x dir + dy = dxyz[1] # Grid spacing in x dir + dz = dxyz[2] # Grid spacing in z dir + ( + vel_model_arr_P, + vel_model_arr_S, + vel_model_x_labels, + vel_model_y_labels, + vel_model_z_labels, + ) = read_1D_vel_model_for_fmm( + oneD_vel_model_z_df, extent_xy_m=extent_xy_m, dxyz=dxyz + ) + # Calculate travel times for all points in the grid: # Setup array for masking/showing where array centre is is: phi = -np.ones(vel_model_arr_P.shape) @@ -79,35 +100,43 @@ def create_LUT(oneD_vel_model_z_df, array_centre_xyz, extent_xy_m=[3000.,3000.], val, y_idx = find_nearest(vel_model_y_labels, array_centre_xyz[1]) val, z_idx = find_nearest(vel_model_z_labels, array_centre_xyz[2]) # Set array_centre_xyz location: - phi[x_idx, y_idx, z_idx] = 1. + phi[x_idx, y_idx, z_idx] = 1.0 # Calculate travel times array to all the points in the grid: - trav_times_grid_P = skfmm.travel_time(phi,vel_model_arr_P,dx=dxyz) - trav_times_grid_S = skfmm.travel_time(phi,vel_model_arr_S,dx=dxyz) - - return trav_times_grid_P, trav_times_grid_S - - - + trav_times_grid_P = skfmm.travel_time(phi, vel_model_arr_P, dx=dxyz) + trav_times_grid_S = skfmm.travel_time(phi, vel_model_arr_S, dx=dxyz) + return trav_times_grid_P, trav_times_grid_S -#----------------------------------------------- End: Define main functions ----------------------------------------------- +# ----------------------------------------------- End: Define main functions ----------------------------------------------- -#----------------------------------------------- Run script ----------------------------------------------- +# ----------------------------------------------- Run script ----------------------------------------------- if __name__ == "__main__": # Example of usage: # (For a 2200 x 1 x 2100 m grid with receivers in a line from the surface to 2 km then along horizontally for 1 km) # Create velocity model: - oneD_vel_model_z_df = pd.DataFrame({'depth': np.arange(0,3000,1.), 'vp': 3500*np.ones(2100), 'vs': 2000*np.ones(2100)}) - station_xs = np.concatenate((1100*np.ones(len(np.arange(1,2000))), np.arange(100,1100)[::-1])) + oneD_vel_model_z_df = pd.DataFrame( + { + "depth": np.arange(0, 3000, 1.0), + "vp": 3500 * np.ones(2100), + "vs": 2000 * np.ones(2100), + } + ) + station_xs = np.concatenate( + (1100 * np.ones(len(np.arange(1, 2000))), np.arange(100, 1100)[::-1]) + ) station_ys = np.zeros(len(station_xs)) - station_zs = np.concatenate((np.arange(1,2000), 2000*np.ones(1000))) - station_xyz_coords_df = pd.DataFrame({'x_m': station_xs, 'y_m': station_ys, 'z_m': station_zs}) - src_xyz = [600., 0., 1500.] + station_zs = np.concatenate((np.arange(1, 2000), 2000 * np.ones(1000))) + station_xyz_coords_df = pd.DataFrame( + {"x_m": station_xs, "y_m": station_ys, "z_m": station_zs} + ) + src_xyz = [600.0, 0.0, 1500.0] # And get travel times: - trav_times_P, trav_times_S = calc_travel_times_fmm(oneD_vel_model_z_df, station_xyz_coords_df, src_xyz, extent_xy_m=[2200.,1.]) + trav_times_P, trav_times_S = calc_travel_times_fmm( + oneD_vel_model_z_df, station_xyz_coords_df, src_xyz, extent_xy_m=[2200.0, 1.0] + ) # And plot result: plt.figure() @@ -116,32 +145,3 @@ def create_LUT(oneD_vel_model_z_df, array_centre_xyz, extent_xy_m=[3000.,3000.], plt.show() print("Finished") - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/rutford_icequake_example/rutford_icequake_example.ipynb b/examples/rutford_icequake_example/rutford_icequake_example.ipynb index de923ab..0954ac8 100644 --- a/examples/rutford_icequake_example/rutford_icequake_example.ipynb +++ b/examples/rutford_icequake_example/rutford_icequake_example.ipynb @@ -40,7 +40,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -59,7 +59,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -77,17 +77,17 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Min. inter-station distance: 0.0193294780285 km\n", - "Therefore, optimal sensitive higher freq.: 98.2954633955 Hz\n", - "Max. inter-station distance: 0.0923463197348 km\n", - "Therefore, optimal sensitive lower freq.: 20.5747235565 Hz\n" + "Min. inter-station distance: 0.019329478028457248 km\n", + "Therefore, optimal sensitive higher freq.: 98.29546339548236 Hz\n", + "Max. inter-station distance: 0.09234631973482979 km\n", + "Therefore, optimal sensitive lower freq.: 20.574723556453616 Hz\n" ] } ], @@ -119,39 +119,6 @@ " " ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Define a suitable event for using in this example:" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "# # Specify data to load:\n", - "# archive = \"inputs/mseed_data_archive\"\n", - "# year = 2020\n", - "# julday = 1\n", - "# hour = 1\n", - "# station_codes_to_use = \"A\"\n", - "# freqmin = 10\n", - "# freqmax = 150\n", - "# mseed_dir = os.path.join(archive, str(year), str(julday).zfill(3))\n", - "\n", - "# # Load data:\n", - "# st = obspy.read(os.path.join(mseed_dir, ''.join((str(year), str(julday).zfill(3), \"_\", str(hour).zfill(2), \"*\", station_codes_to_use, \"*\"))))\n", - "\n", - "# # And filter:\n", - "# st.filter('bandpass', freqmin=freqmin, freqmax=freqmax)\n", - "\n", - "# # And trim:\n", - "# st.trim(starttime=starttime, endtime=endtime)\n" - ] - }, { "cell_type": "code", "execution_count": null, @@ -168,7 +135,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -1296,9 +1263,7 @@ { "cell_type": "code", "execution_count": 18, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -2342,9 +2307,7 @@ { "cell_type": "code", "execution_count": 21, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -4360,9 +4323,7 @@ { "cell_type": "code", "execution_count": 24, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -5375,9 +5336,7 @@ { "cell_type": "code", "execution_count": 25, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -10349,7 +10308,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -10363,7 +10322,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.6" + "version": "3.9.16" } }, "nbformat": 4,