-
Notifications
You must be signed in to change notification settings - Fork 22
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Structural Changes for SSST-relocation * SSST Relocation * Added a new class for manipulating parametric-data * Updated picking routines that now use the new parametric-data class * Entire file-content is no longer loaded into memory. * Arrivals are now labelled by event-source, e.g GA, ISC etc. * Updating ellip-corr hash * Updating ellip-corr * New Classes for SSST-Relocation * Added a new class for travel-time interpolation * Added a barebone ssst-relocator class * Added a parallel function for tt interpolation * Removing MPI-parallel Catalogue Ingestion MPI-Allgatherv is prone to data corruption when tens of processors are used to gather arrays on all ranks amounting to hundreds of MB of data. Catalogue ingestion now takes place on rank 0 and events and arrivals are saved to disk -- subsequently, all other ranks simply load these events and arrivals from disk. * Minor typos * Disk-based Sync Across MPI Ranks * Added a disk-based sync approach for variables across MPI-ranks. The rationale behind is documented in Jira. * Added methods for computing ellipticity-correction * New class for vectorized computation of ellipticity corrections. * Added a function to coalesce discrepant network codes * Added a function for parallel residual computation * Adding ellipticity correction table * Adding support for masked local computations on each rank * Parallelized phase-redefinition * Added parallel ssst-correction calculator * Moving function for coalescing network codes to base class * Adding some documentation * Working draft of SSST-relocator * Refactored SSST-relocator * refactored and cleaned up SSST-relocator * minor changes to travel_time and ellipticity modules * Config Parser and Miscellaneous Changes * Added config parser for relocation/redefinition * Added support for dropping automatic phases during that do not meet quality criteria during relocation * Adding SST Relocation * Added SST relocation as a potential first step * Optimized parallel sync * Miscellaneous Optimizations * Arrivals whose phases are redefined without their origins being relocated are now handled such that SST corrections (if available) are factored in during phase redefinition in each SSST-relocation iteration. * Adding a cli * Minor cleanup * Adding Plotting Routines * Minor cleanups in ssst_relocator and moved CLI to ssst_relocate * Added routines for generating detailed plots in ssst_relocate * Added Arrivals Exporter * Some minor changes in ssst_relocator * Added export_arrivals for exporting cleansed arrivals that meet user-defined criteria * Miscellaneous Changes * cwb2asdf.py can now convert semi-permanent mseed files into ASDF * extract_event_traces.py now reflect recent changes made to the picking workflow * added more plots to export_arrivals.py
- Loading branch information
Showing
32 changed files
with
5,233 additions
and
149 deletions.
There are no files selected for viewing
Submodule ellip-corr
updated
8 files
| +8 −0 | Makefile | |
| +49 −0 | PyEllipCorr.py | |
| +10 −28 | README.md | |
| +0 −0 | __init__.py | |
| +3 −3 | ellip/Makefile | |
| +185 −145 | ellip/pyellip.f | |
| +0 −86 | setup.py | |
| +49 −36 | tests/test_ellipticity_corr.py |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,162 @@ | ||
| from collections import defaultdict | ||
| import numpy as np | ||
| from scipy.interpolate import RegularGridInterpolator | ||
| import traceback | ||
| import os | ||
|
|
||
| class Phase: | ||
| def __init__(self, ecdist, depth, tau0, tau1, tau2): | ||
| self.ecdist = ecdist | ||
| self.depth = depth | ||
| self.tau0 = tau0 | ||
| self.tau1 = tau1 | ||
| self.tau2 = tau2 | ||
|
|
||
| self.tau0io = RegularGridInterpolator((self.ecdist, self.depth), self.tau0, method='linear', | ||
| bounds_error=False, fill_value=0.) | ||
| self.tau1io = RegularGridInterpolator((self.ecdist, self.depth), self.tau1, method='linear', | ||
| bounds_error=False, fill_value=0.) | ||
| self.tau2io = RegularGridInterpolator((self.ecdist, self.depth), self.tau2, method='linear', | ||
| bounds_error=False, fill_value=0.) | ||
| # end func | ||
| # end class | ||
|
|
||
| class Ellipticity: | ||
| def __init__(self): | ||
| self._elcordir_tbl = os.path.dirname(os.path.abspath(__file__)) + '/ellip/elcordir.tbl' | ||
| self.depth = np.array([0.0, 100.0, 200.0, 300.0, 500.0, 700.0], dtype='f4') | ||
| self.phases = defaultdict(list) | ||
| self.phase_alias = { b'Pg': b'Pup', | ||
| b'Sg': b'Sup', | ||
| b'pPg': b'pP', | ||
| b'sPg': b'sP', | ||
| b'pSg': b'pS', | ||
| b'sSg': b'sS', | ||
| b'Pb': b'P', | ||
| b'Sb': b'S', | ||
| b'pPb': b'pP', | ||
| b'sPb': b'sP', | ||
| b'pSb': b'pS', | ||
| b'sSb': b'sS', | ||
| b'Pn': b'P', | ||
| b'Sn': b'S', | ||
| b'pPn': b'pP', | ||
| b'sPn': b'sP', | ||
| b'pSn': b'pS', | ||
| b'sSn': b'sS', | ||
| b'SPn': b'SP', | ||
| b'SPb': b'SP', | ||
| b'SPg': b'SP', | ||
| b'SnP': b'SP', | ||
| b'PSn': b'PS', | ||
| b'PnPn': b'PP', | ||
| b'SnSn': b'SS', | ||
| b'p': b'Pup', | ||
| b's': b'Sup', | ||
| b'Pdif': b'Pdiff', | ||
| b'Sdif': b'Sdiff' } | ||
| self._parse_table() | ||
| # end func | ||
|
|
||
| def _parse_table(self): | ||
| line = open(self._elcordir_tbl).read() | ||
| ll = len(line) | ||
|
|
||
| indices = np.arange(0, ll+80, 80) | ||
| rows = list() | ||
|
|
||
| for i in range(len(indices)-1): | ||
| row = line[indices[i]:indices[i+1]].split() | ||
| rows.append(row) | ||
| #end for | ||
|
|
||
| phase = None | ||
| ecdist = list() | ||
| tau = 0 | ||
| tau0 = list() | ||
| tau1 = list() | ||
| tau2 = list() | ||
| for irow, row in enumerate(rows): | ||
| if len(row) == 4: | ||
| if irow: | ||
| ecdist = np.array(ecdist).astype('f4') | ||
| tau0 = np.array(tau0).astype('f4') | ||
| tau1 = np.array(tau1).astype('f4') | ||
| tau2 = np.array(tau2).astype('f4') | ||
| self.phases[phase] = Phase(ecdist, self.depth, tau0, tau1, tau2) | ||
| ecdist = list() | ||
| tau0 = list() | ||
| tau1 = list() | ||
| tau2 = list() | ||
| #end if | ||
| phase = (row[0]).encode() # convert to byte-string | ||
| elif len(row) == 1: | ||
| ecdist.append(float(row[0])) | ||
| elif len(row) == 6 and tau == 0: | ||
| tau0.append(row) | ||
| tau = (tau + 1) % 3 | ||
| elif len(row) == 6 and tau == 1: | ||
| tau1.append(row) | ||
| tau = (tau + 1) % 3 | ||
| elif len(row) == 6 and tau == 2: | ||
| tau2.append(row) | ||
| tau = (tau + 1) % 3 | ||
| #end if | ||
| #end for | ||
|
|
||
| for pnew, palias in self.phase_alias.items(): | ||
| if(palias in self.phases.keys()): self.phases[pnew] = self.phases[palias] | ||
| # end for | ||
| # end func | ||
|
|
||
| def get_correction(self, phase, ecdist, depth_km, elat, azim): | ||
| azim = np.radians(azim) | ||
| ecolat = np.radians(90 - elat) | ||
| s3 = np.sqrt(3.0) / 2.0 | ||
|
|
||
| try: | ||
| if(type(phase) == np.ndarray): | ||
| result = np.zeros(len(phase), dtype='f4') | ||
|
|
||
| corr_count = 0 | ||
| for cphase in self.phases.keys(): | ||
| indices = np.argwhere(cphase == phase).flatten() | ||
|
|
||
| if(len(indices)): | ||
| sc0 = 0.25 * (1.0 + 3.0 * np.cos(2.0 * ecolat[indices])) | ||
| sc1 = s3 * np.sin(2.0 * ecolat[indices]) | ||
| sc2 = s3 * np.sin(ecolat[indices]) * np.sin(ecolat[indices]) | ||
|
|
||
| tau0 = self.phases[cphase].tau0io((ecdist[indices], depth_km[indices])) | ||
| tau1 = self.phases[cphase].tau1io((ecdist[indices], depth_km[indices])) | ||
| tau2 = self.phases[cphase].tau2io((ecdist[indices], depth_km[indices])) | ||
|
|
||
| result[indices] = sc0*tau0 + sc1*np.cos(azim[indices])*tau1 + \ | ||
| sc2*np.cos(2.0*azim[indices])*tau2 | ||
| corr_count += len(indices) | ||
| # end if | ||
| # end for | ||
| #if(corr_count < len(phase)): print('Warning: some phases {} not found..'.format(set(list(phase)) - set(self.phases.keys()))) | ||
|
|
||
| return result | ||
| else: | ||
| result = 0. | ||
| sc0 = 0.25 * (1.0 + 3.0 * np.cos(2.0 * ecolat)) | ||
| sc1 = s3 * np.sin(2.0 * ecolat) | ||
| sc2 = s3 * np.sin(ecolat) * np.sin(ecolat) | ||
|
|
||
| tau0 = self.phases[phase].tau0io((ecdist, depth_km)) | ||
| tau1 = self.phases[phase].tau1io((ecdist, depth_km)) | ||
| tau2 = self.phases[phase].tau2io((ecdist, depth_km)) | ||
|
|
||
| result = sc0 * tau0 + sc1 * np.cos(azim) * tau1 + \ | ||
| sc2 * np.cos(2.0 * azim) * tau2 | ||
|
|
||
| return result | ||
| # end if | ||
| except Exception as e: | ||
| print(traceback.format_exc()) | ||
| raise ValueError('Phase not found..') | ||
| # end try | ||
| # end func | ||
| # end class |
Oops, something went wrong.