diff --git a/data/examples/20210809074550.tgz b/data/examples/20210809074550.tgz new file mode 100644 index 000000000..47830d421 Binary files /dev/null and b/data/examples/20210809074550.tgz differ diff --git a/data/examples/unpack.bash b/data/examples/unpack.bash index fa640c856..da3897878 100755 --- a/data/examples/unpack.bash +++ b/data/examples/unpack.bash @@ -9,7 +9,7 @@ cd $(dirname ${BASH_SOURCE[0]}) wd=$PWD for filename in \ - 20090407201255351.tgz SPECFEM3D_SGT.tgz SPECFEM3D_SAC.tgz; + 20090407201255351.tgz 20210809074550.tgz 20SPECFEM3D_SGT.tgz SPECFEM3D_SAC.tgz; do cd $wd cd $(dirname $filename) diff --git a/examples/GridSearch.Force.py b/examples/GridSearch.Force.py new file mode 100755 index 000000000..15cc1f649 --- /dev/null +++ b/examples/GridSearch.Force.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python + +import os +import numpy as np + +from mtuq import read, open_db, download_greens_tensors +from mtuq.event import Origin +from mtuq.graphics import plot_data_greens1, plot_misfit_force +from mtuq.grid import ForceGridRegular +from mtuq.grid_search import grid_search +from mtuq.misfit import Misfit +from mtuq.process_data import ProcessData +from mtuq.util import fullpath, merge_dicts, save_json +from mtuq.util.cap import parse_station_codes, Trapezoid +from mtuq.misfit.waveform import calculate_norm_data + + + +if __name__=='__main__': + # + # Carries out grid search over 64,000 double couple moment tensors + # + # USAGE + # mpirun -n python GridSearch.DoubleCouple.py + # + # For a simpler example, see SerialGridSearch.DoubleCouple.py, + # which runs the same inversion in serial + # + + + # + # We will investigate the source process of an Mw~4 earthquake using data + # from a regional seismic array + # + + path_data= fullpath('data/examples/20210809074550/*[ZRT].sac') + path_weights= fullpath('data/examples/20210809074550/weights.dat') + event_id= '20210809074550' + model= 'ak135' + + + # + # We are only using surface waves in this example. Check out the DC or FMT examples for multi-mode inversions. + # + + process_sw = ProcessData( + filter_type='Bandpass', + freq_min=0.025, + freq_max=0.0625, + pick_type='taup', + taup_model=model, + window_type='surface_wave', + window_length=150., + capuaf_file=path_weights, + ) + + + # + # For our objective function, we will use the L2 norm of the misfit between + # observed and synthetic waveforms. + # + + misfit_sw = Misfit( + norm='L2', + time_shift_min=-35., + time_shift_max=+35., + time_shift_groups=['ZR','T'], + ) + + + # + # User-supplied weights control how much each station contributes to the + # objective function + # + + station_id_list = parse_station_codes(path_weights) + + + # + # Next, we specify the moment tensor grid and source-time function + # + + grid = ForceGridRegular(magnitudes_in_N=10**np.arange(9,12.1,0.1), npts_per_axis=90) + + + # In this example, we use a simple trapezoidal source-time function with a + # rise time of 4.5 second and a half-duration of 6.75 seconds, obtained from + # earthquake scaling law. + wavelet = Trapezoid( + magnitude=8) + + + # + # Origin time and location will be fixed. For an example in which they + # vary, see examples/GridSearch.DoubleCouple+Magnitude+Depth.py + # + # See also Dataset.get_origins(), which attempts to create Origin objects + # from waveform metadata + # + + origin = Origin({ + 'time': '2021-08-09T07:45:50.000000Z', + 'latitude': 61.24, + 'longitude': -147.96, + 'depth_in_m': 0, + }) + + + from mpi4py import MPI + comm = MPI.COMM_WORLD + + + # + # The main I/O work starts now + # + + if comm.rank==0: + print('Reading data...\n') + data = read(path_data, format='sac', + event_id=event_id, + station_id_list=station_id_list, + tags=['units:m', 'type:velocity']) + + + data.sort_by_distance() + stations = data.get_stations() + + + print('Processing data...\n') + data_sw = data.map(process_sw) + + + print('Reading Greens functions...\n') + greens = download_greens_tensors(stations, origin, model, include_mt=False, include_force=True) + + print('Processing Greens functions...\n') + greens.convolve(wavelet) + greens_sw = greens.map(process_sw) + + + else: + stations = None + data_sw = None + greens_sw = None + + + stations = comm.bcast(stations, root=0) + data_sw = comm.bcast(data_sw, root=0) + greens_sw = comm.bcast(greens_sw, root=0) + + # + # The main computational work starts now + # + + if comm.rank==0: + print('Evaluating surface wave misfit...\n') + + results_sw = grid_search( + data_sw, greens_sw, misfit_sw, origin, grid) + + if comm.rank==0: + # Computing the norm of the data for the misfit normalization + norm_sw = calculate_norm_data(data_sw, misfit_sw.norm, ['Z','R','T']) + + results = results_sw/norm_sw + + # + # Collect information about best-fitting source + # + + # `grid` index corresponding to minimum misfit force + idx = results.source_idxmin() + + # Force object + best_force = grid.get(idx) + + # dictionary of best force direction (F0, phi, h) + direction_dict = grid.get_dict(idx) + + # dictionary of Fi parameters + force_dict = best_force.as_dict() + + merged_dict = merge_dicts( + direction_dict, force_dict, origin) + + # + # Generate figures and save results + # + + print('Generating figures...\n') + + plot_data_greens1(event_id+'FMT_waveforms1.png', + data_sw, greens_sw, process_sw, misfit_sw, stations, origin, best_force, direction_dict) + + plot_misfit_force(event_id+'DC_misfit_force.png', results) + + print('Saving results...\n') + + # save best-fitting source + save_json(event_id+'Force_solution.json', merged_dict) + + # save misfit surface + results.save(event_id+'DC_misfit.nc') + + print('\nFinished\n') \ No newline at end of file diff --git a/mtuq/io/clients/syngine.py b/mtuq/io/clients/syngine.py index 842f26323..e5a22e71d 100644 --- a/mtuq/io/clients/syngine.py +++ b/mtuq/io/clients/syngine.py @@ -1,6 +1,7 @@ import obspy import numpy as np +import os from obspy.core import Stream from mtuq.greens_tensor.syngine import GreensTensor @@ -9,7 +10,7 @@ from mtuq.util import unzip from mtuq.util.syngine import download_unzip_mt_response, download_force_response,\ resolve_model,\ - GREENS_TENSOR_FILENAMES, SYNTHETICS_FILENAMES + GREENS_TENSOR_FILENAMES, SYNTHETICS_FILENAMES_BX, SYNTHETICS_FILENAMES_MX @@ -95,11 +96,23 @@ def _get_greens_tensor(self, station=None, origin=None): dirnames += [unzip(filename)] for _i, dirname in enumerate(dirnames): - for filename in SYNTHETICS_FILENAMES: - stream += obspy.read(dirname+'/'+filename, format='sac') + # Attempt to read using the BX naming convention + files_read = False + for filename in SYNTHETICS_FILENAMES_BX: + file_path = os.path.join(dirname, filename) + if os.path.isfile(file_path): + stream += obspy.read(file_path, format='sac') + stream[-1].stats.channel = stream[-1].stats.channel[-1] + str(_i) + files_read = True + + # If no files were read with the BX naming convention, try the MX naming convention + if not files_read: + for filename in SYNTHETICS_FILENAMES_MX: + file_path = os.path.join(dirname, filename) + if os.path.isfile(file_path): + stream += obspy.read(file_path, format='sac') + stream[-1].stats.channel = stream[-1].stats.channel[-1] + str(_i) - # overwrite channel attribute - stream[-1].stats.channel = stream[-1].stats.channel[-1]+str(_i) # what are the start and end times of the data? diff --git a/mtuq/util/syngine.py b/mtuq/util/syngine.py index efe4cfb9a..2b8c51f68 100644 --- a/mtuq/util/syngine.py +++ b/mtuq/util/syngine.py @@ -24,12 +24,18 @@ 'greensfunction_XX.GF001..TDS.sac', ] -SYNTHETICS_FILENAMES = [ +SYNTHETICS_FILENAMES_BX = [ 'XX.S0001.SE.BXZ.sac', 'XX.S0001.SE.BXR.sac', 'XX.S0001.SE.BXT.sac', ] +SYNTHETICS_FILENAMES_MX = [ + 'XX.S0001.SE.MXZ.sac', + 'XX.S0001.SE.MXR.sac', + 'XX.S0001.SE.MXT.sac', + ] + SYNGINE_MODELS = [ 'ak135f_2s', 'ak135f_5s',