From e51f65b985f84c7ad3734cfd582a8dce95d6e7c7 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 00:01:11 -0700 Subject: [PATCH 1/8] Updated URLs --- docs/install/issues.rst | 2 +- tests/check_urls.bash | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/install/issues.rst b/docs/install/issues.rst index 14fc4664e..8ee8368a3 100644 --- a/docs/install/issues.rst +++ b/docs/install/issues.rst @@ -37,7 +37,7 @@ Older versions of the conda package manager can be very slow. For a significant conda update -n base conda -For reference, the largest potential speed up comes from the new `mamba `_ dependency solver, which was `adopted `_ in the 23.10 release. +For reference, the largest potential speed up comes from the new `mamba `_ dependency solver, which was `adopted `_ in the 23.10 release. diff --git a/tests/check_urls.bash b/tests/check_urls.bash index 9053377f3..9bef6b62d 100755 --- a/tests/check_urls.bash +++ b/tests/check_urls.bash @@ -20,6 +20,8 @@ URLS="\ https://docs.xarray.dev/en/stable/generated/xarray.DataArray.html\ https://github.com/uafgeotools/mtuq/blob/master/docs/user_guide/05/code/gallery_mt.py https://github.com/uafgeotools/mtuq/blob/master/docs/user_guide/05/code/gallery_force.py + https://conda.org/blog/2023-11-06-conda-23-10-0-release\ + https://www.anaconda.com/blog/a-faster-conda-for-a-growing-community\ " From ef33c26737cb5943d89126e7b628e157993e2fdb Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 13:19:22 -0700 Subject: [PATCH 2/8] New Greens function docs pages --- docs/user_guide/03.rst | 39 ++++++++++++++++------------ docs/user_guide/03/receiver_side.rst | 15 +++++++++++ docs/user_guide/03/source_side.rst | 15 +++++++++++ 3 files changed, 53 insertions(+), 16 deletions(-) create mode 100644 docs/user_guide/03/receiver_side.rst create mode 100644 docs/user_guide/03/source_side.rst diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index d66655088..4918a826b 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -11,6 +11,24 @@ Role of Green's functions By combining Green's functions, it is possible to obtain synthetics from any moment tensor or force source. Generating synthetics in this way is useful because it provides cost savings compared with on-the-fly wavefield simulations. Synthetics can then be compared with data to determine a best-fitting source. +Green's function conventions +---------------------------- + +MTUQ consistently uses an `up-south-east` `convention `_ for internally storing Green's functions (as well as moment tensors and forces). + +A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. + +Under the hood, MTUQ deals with a variety of externel conventions related to + +- the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media) + +- the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`; see `ObsPy documentation `_ for more information) + +- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) + +A better sense of what's involved can be obtained by browsing the `source code `_ and references therein. + + `GreensTensor` objects ---------------------- @@ -81,7 +99,7 @@ Computing Green's functions from 3D Earth models MTUQ currently supports 3D Green's functions from SPECFEM3D/3D_GLOBE. -To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: +To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: .. code :: @@ -114,27 +132,16 @@ To open an SPECFEM3D database client in MTUQ: db = open_db(path_SPECFEM3D_output_directory, format="SPECFEM3D") -Once opened, a SPECFEM3D database client can be used to generate `GreensTensor `_ objects as follows: +Once opened, a SPECFEM3D/3D_GLOBE database client can be used to generate `GreensTensor `_ objects as follows: .. code:: greens_tensors = db.get_greens_tensors(stations, origin) +For more information, please see: -Green's function conventions ----------------------------- - -A variety of Green's function conventions exist. Figuring out which are used in a particular application can be challenging because it depends on - -- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) - -- the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media) - -- the choice of local Cartesian basis conventions (for example, some authors employ `up-south-east`, others `north-east-down`; see `ObsPy documentation `_ for more information) - -A major goal is to avoid exposing users to unnecessary basis complexity. MTUQ accomplishes this by understanding external formats and converting to a common internal format that works for both 1D and 3D media. - -For internally storing moment tensors, forces, and Green's functions, MTUQ consistently uses an `up-south-east` Cartesian convention. +`Detailed guide to source-side 3D Green's functions (under construction) `_ +`Detailed guide to receiver-side 3D Green's functions (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst new file mode 100644 index 000000000..f1a4be8a7 --- /dev/null +++ b/docs/user_guide/03/receiver_side.rst @@ -0,0 +1,15 @@ + +.. warning:: + + This page is currently just as stub. + + To fill in missing documentation, feel free to submit a pull request. + + +Working example +=============== + +A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: + +`test_greens_SPECFEM3D_SGT.py `_ + diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst new file mode 100644 index 000000000..56fdfc872 --- /dev/null +++ b/docs/user_guide/03/source_side.rst @@ -0,0 +1,15 @@ + +.. warning:: + + This page is currently just as stub. + + To fill in missing documentation, feel free to submit a pull request. + + +Working example +=============== + +A working of example using source-side 3D Green's functions from SPECFEM3D_GLOBE in a moment tensor inversion: + +`test_greens_SPECFEM3D_SAC.py `_ + From dc5116b9c92f36303a2f7890bce93825c95540f2 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 17:33:58 -0700 Subject: [PATCH 3/8] Filled in source-side Greens function notes --- docs/user_guide/03.rst | 33 ++++------- docs/user_guide/03/receiver_side.rst | 4 +- docs/user_guide/03/source_side.rst | 85 +++++++++++++++++++++++++++- 3 files changed, 98 insertions(+), 24 deletions(-) diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index 4918a826b..ce045e12f 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -18,15 +18,7 @@ MTUQ consistently uses an `up-south-east` `convention `_ for more information) - -- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) - -A better sense of what's involved can be obtained by browsing the `source code `_ and references therein. +Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. @@ -35,8 +27,7 @@ A better sense of what's involved can be obtained by browsing the `source code < `GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given station and origin. -Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics for the given station. - +Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. @@ -52,7 +43,7 @@ To download AK135 Green's functions and generate MTUQ `GreensTensor` objects: from mtuq import download_greens_functions greens_tensors = download_greens_functions(stations, origins, model="ak135f_2s") -A limitation of syngine is that Green's functions can be downloaded only on a station-by-station basis, not over an entire subset of Earth. An alternative that avoids this limitation is to compute your own Green's functions. +A limitation of syngine is that Green's functions can be downloaded only on a station-by-station basis, not over an entire area or volume. An alternative that avoids this limitation is to compute your own Green's functions. @@ -62,7 +53,7 @@ Computing Green's functions from 1D Earth models MTUQ supports the following 1D Green's functions formats: AxiSEM (preferred) and FK. -`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF formated needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. +`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. To open an AxiSEM database client in MTUQ: @@ -99,11 +90,11 @@ Computing Green's functions from 3D Earth models MTUQ currently supports 3D Green's functions from SPECFEM3D/3D_GLOBE. -To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: +To generate a full set of Green's functions for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: .. code :: - {basedir}/{event_id}/ + {event_id}/ {net}.{sta}.{loc}.Z.Mrr.sac {net}.{sta}.{loc}.Z.Mtt.sac {net}.{sta}.{loc}.Z.Mpp.sac @@ -124,24 +115,24 @@ To generate a complete Green's function database for a given hypocenter and dept {net}.{sta}.{loc}.T.Mtp.sac -To open an SPECFEM3D database client in MTUQ: +To open a SPECFEM3D/3D_GLOBE database client: -.. code :: +.. code:: from mtuq import open_db db = open_db(path_SPECFEM3D_output_directory, format="SPECFEM3D") -Once opened, a SPECFEM3D/3D_GLOBE database client can be used to generate `GreensTensor `_ objects as follows: +Once opened, the database client can be used to generate `GreensTensor `_ objects as follows: .. code:: greens_tensors = db.get_greens_tensors(stations, origin) -For more information, please see: +For more information, also see: -`Detailed guide to source-side 3D Green's functions (under construction) `_ +`Source-side SPECFEM3D/3D_GLOBE Green's functions `_ -`Detailed guide to receiver-side 3D Green's functions (under construction) `_ +`Receiver-side SPECFEM3D/3D_GLOBE Green's functions (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index f1a4be8a7..deaa470d3 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -6,8 +6,8 @@ To fill in missing documentation, feel free to submit a pull request. -Working example -=============== +Notes on receiver-side 3D Green's functions +=========================================== A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index 56fdfc872..605db9db1 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -6,8 +6,91 @@ To fill in missing documentation, feel free to submit a pull request. +Notes on source-side 3D Green's functions +========================================= + +**File format** + +- Currently, MTUQ reads source-side 3D Green's functions from SAC files + +- SAC output format is natively supported by SPECFEM3D_GLOBE, but not SPECFEM3D Cartesian (output from the latter can be manually converted) + + +**Units convention** + +- MTUQ treats each individual SAC file as the response in meters to a unit (1 Newton-meter) force couple + +- Users must ensure that source-side Green's functions are normalized according to the above units + + +**Basis convention** + +- For moment tensors, MTUQ uses an `up-south-east` convention, identical to one the used by SPECFEM3D_GLOBE + + +**Origin time convention** + +- For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time + +- MTUQ uses the `REF_TIME` header from the SPECFEM3D_GLOBE SAC output files, which gives the peak excitation of the source relative to the simulation start time + +- MTUQ ignores the origin time given in the CMTSOLUTION file and `ORIGIN_TIME` header + + +**Depth searches (experimental)** + +- Only depth searches are possible with source-side 3D Green's functions (no other hypocenter parameters) + +- The current `depth search `_ implementation is somewhat crude and experimental; consider local modifcations to suit your needs + +- To allow depth searches, create subdirectories for each centroid depth, as below + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + ... + + Working example -=============== +--------------- A working of example using source-side 3D Green's functions from SPECFEM3D_GLOBE in a moment tensor inversion: From bcebc6b5a8fbb1bd3de77c952467fa942b14e17d Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 9 Feb 2024 20:24:01 -0700 Subject: [PATCH 4/8] Revised comments/docstrings --- docs/install/issues.rst | 2 +- docs/user_guide/03.rst | 6 ++++-- docs/user_guide/03/receiver_side.rst | 2 +- docs/user_guide/03/source_side.rst | 2 +- mtuq/greens_tensor/base.py | 2 +- tests/test_time_shifts.py | 9 ++++----- 6 files changed, 12 insertions(+), 11 deletions(-) diff --git a/docs/install/issues.rst b/docs/install/issues.rst index 8ee8368a3..5ac380216 100644 --- a/docs/install/issues.rst +++ b/docs/install/issues.rst @@ -31,7 +31,7 @@ We note that some versions of GMT and ObsPy do not plot `full moment tensors `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. +`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. + +To generate AxiSEM synthetics in a format MTUQ supports, follow the instructions in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." To open an AxiSEM database client in MTUQ: diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index deaa470d3..10994583b 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -1,7 +1,7 @@ .. warning:: - This page is currently just as stub. + This page is currently just a stub. To fill in missing documentation, feel free to submit a pull request. diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index 605db9db1..df27bdaf9 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -1,7 +1,7 @@ .. warning:: - This page is currently just as stub. + This page is currently just a stub. To fill in missing documentation, feel free to submit a pull request. diff --git a/mtuq/greens_tensor/base.py b/mtuq/greens_tensor/base.py index 34c5cfea1..3ec868d4a 100644 --- a/mtuq/greens_tensor/base.py +++ b/mtuq/greens_tensor/base.py @@ -471,7 +471,7 @@ def sort_by_azimuth(self, reverse=False): def sort_by_function(self, function, reverse=False): - """ Sorts in-place using the python built-in `sort` + """ Sorts in-place by user-supplied function """ self.sort(key=function, reverse=reverse) diff --git a/tests/test_time_shifts.py b/tests/test_time_shifts.py index 0a6b31374..2c5db6a25 100644 --- a/tests/test_time_shifts.py +++ b/tests/test_time_shifts.py @@ -161,13 +161,16 @@ def __call__(self, traces): return super(ProcessData, self).__call__( traces, station=station, origin=origin) + + # + # plotting functions + # def _plot_dat(axis, t, data, attrs, pathspec='-k'): stream = data[0] trace = data[0][0] stats = trace.stats axis.plot(t, trace.data, pathspec) - def _plot_syn(axis, t, data, attrs, pathspec='-r'): stream = data[0] trace = data[0][0] @@ -176,7 +179,6 @@ def _plot_syn(axis, t, data, attrs, pathspec='-r'): idx2 = attrs.idx_stop axis.plot(t, trace.data[idx1:idx2], pathspec) - def _annotate(axis, attrs): text = 'static_shift: %.1f' % attrs.static_shift pyplot.text(0.02, 0.85, text, transform=axis.transAxes) @@ -187,15 +189,12 @@ def _annotate(axis, attrs): text = 'total_shift: %.1f' % attrs.total_shift pyplot.text(0.02, 0.55, text, transform=axis.transAxes) - def _get_synthetics(greens, mt): return greens.get_synthetics(mt, components=['Z']) - def _get_attrs(data, greens, misfit): return misfit.collect_attributes(data, greens, mt)[0]['Z'] - def _get_time_sampling(data): stream = data[0] t1 = float(stream[0].stats.starttime) From a138cf7c2db3613af51de59501460181b22bd945 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 11 Mar 2024 16:26:40 -0600 Subject: [PATCH 5/8] Added CPS GreensTensor subclass --- mtuq/greens_tensor/CPS.py | 109 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 mtuq/greens_tensor/CPS.py diff --git a/mtuq/greens_tensor/CPS.py b/mtuq/greens_tensor/CPS.py new file mode 100644 index 000000000..ba229c6e0 --- /dev/null +++ b/mtuq/greens_tensor/CPS.py @@ -0,0 +1,109 @@ + +import obspy +import numpy as np + +from mtuq.greens_tensor.base import GreensTensor as GreensTensorBase + + + +print('WARNING: CPS Greens functions are not fully tested yet') + + + +class GreensTensor(GreensTensorBase): + """ + FK Green's tensor object + + Overloads base class with machinery for working with CPS-style + Green's functions + + """ + def __init__(self, *args, **kwargs): + super(GreensTensor, self).__init__(*args, **kwargs) + + if 'type:greens' not in self.tags: + self.tags += ['type:greens'] + + if 'type:velocity' not in self.tags: + self.tags += ['type:velocity'] + + if 'units:m' not in self.tags: + self.tags += ['units:m'] + + def _precompute(self): + """ Computes NumPy arrays used by get_synthetics + """ + if self.include_mt: + self._precompute_mt() + + if self.include_force: + self._precompute_force() + + + def _precompute_mt(self): + """ Recombines CPS time series so they can be used in straightforward + linear combination with Mrr,Mtt,Mpp,Mrt,Mrp,Mtp + """ + + array = self._array + phi = np.deg2rad(self.azimuth) + _j = 0 + + # The formulas below were obtained by reverse engineering FK + + for _i, component in enumerate(self.components): + if component=='Z': + ZSS = self.select(channel="ZSS")[0].data + ZDS = self.select(channel="ZDS")[0].data + ZDD = self.select(channel="ZDD")[0].data + ZEP = self.select(channel="ZEP")[0].data + array[_i, _j+0, :] = ZSS/2. * np.cos(2*phi) - ZDD/6. + ZEP/3. + array[_i, _j+1, :] = -ZSS/2. * np.cos(2*phi) - ZDD/6. + ZEP/3. + array[_i, _j+2, :] = ZDD/3. + ZEP/3. + array[_i, _j+3, :] = ZSS * np.sin(2*phi) + array[_i, _j+4, :] = ZDS * np.cos(phi) + array[_i, _j+5, :] = ZDS * np.sin(phi) + + elif component=='R': + RSS = self.select(channel="RSS")[0].data + RDS = self.select(channel="RDS")[0].data + RDD = self.select(channel="RDD")[0].data + REP = self.select(channel="REP")[0].data + array[_i, _j+0, :] = RSS/2. * np.cos(2*phi) - RDD/6. + REP/3. + array[_i, _j+1, :] = -RSS/2. * np.cos(2*phi) - RDD/6. + REP/3. + array[_i, _j+2, :] = RDD/3. + REP/3. + array[_i, _j+3, :] = RSS * np.sin(2*phi) + array[_i, _j+4, :] = RDS * np.cos(phi) + array[_i, _j+5, :] = RDS * np.sin(phi) + + elif component=='T': + TSS = self.select(channel="TSS")[0].data + TDS = self.select(channel="TDS")[0].data + array[_i, _j+0, :] = TSS/2. * np.sin(2*phi) + array[_i, _j+1, :] = -TSS/2. * np.sin(2*phi) + array[_i, _j+2, :] = 0. + array[_i, _j+3, :] = -TSS * np.cos(2*phi) + array[_i, _j+4, :] = TDS * np.sin(phi) + array[_i, _j+5, :] = -TDS * np.cos(phi) + + else: + raise ValueError + + # + # CPS uses a north-east-down basis convention, while mtuq uses an + # up-south-east basis convention, so a permutation is necessary + # + array_copy = array.copy() + array[:, 0, :] = array_copy[:, 2, :] + array[:, 1, :] = array_copy[:, 0, :] + array[:, 2, :] = array_copy[:, 1, :] + array[:, 3, :] = array_copy[:, 4, :] + array[:, 4, :] = -array_copy[:, 5, :] + array[:, 5, :] = -array_copy[:, 3, :] + + + def _precompute_force(self): + raise NotImplementedError() + + + From 9f162f3dba3d3b8c9a52e214d9e8ff496a7af7e1 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Wed, 13 Mar 2024 14:35:21 -0600 Subject: [PATCH 6/8] Typo fixes --- docs/README | 2 +- mtuq/misfit/polarity.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/README b/docs/README index 732d3375c..5e1130515 100644 --- a/docs/README +++ b/docs/README @@ -14,7 +14,7 @@ HOW TO UPDATE DOCS >> ./build.bash -3. If build fails with "Encountered unknown tag 'do', then try the following: +3. If build fails with "Encountered unknown tag 'do'", try the following: Modify the sphinx autosummary source code by adding 'jinja2.ext.do' to the jinja Environment() extensions list diff --git a/mtuq/misfit/polarity.py b/mtuq/misfit/polarity.py index 401284de3..08021f4ed 100644 --- a/mtuq/misfit/polarity.py +++ b/mtuq/misfit/polarity.py @@ -58,9 +58,9 @@ class PolarityMisfit(object): .. rubric:: Other input arguments that may be required, depending on the above ``taup_model`` (`str`): Name of built-in ObsPy TauP model or path to custom - ObsPy TauP model, required for `type=taup` + ObsPy TauP model, required for `method=taup` - ``FK_database`` (`str`): Path to FK database, required for for `type=FK_metadata`. + ``FK_database`` (`str`): Path to FK database, required for for `method=FK_metadata`. .. note:: From 0af5838416c2794c79af3ac75e89d1ca0cfb38fa Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Mar 2024 15:49:49 -0600 Subject: [PATCH 7/8] Revised documentation --- docs/user_guide/03.rst | 62 ++++++------ docs/user_guide/03/receiver_side.rst | 9 +- docs/user_guide/03/source_side.rst | 120 ++++++++++++++++++++---- docs/user_guide/05/gallery_force.rst | 4 + docs/user_guide/05/gallery_mt.rst | 6 ++ docs/user_guide/06/trace_attributes.rst | 2 +- 6 files changed, 149 insertions(+), 54 deletions(-) diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index 405acc3d1..90b119556 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -10,24 +10,21 @@ Role of Green's functions By combining Green's functions, it is possible to obtain synthetics from any moment tensor or force source. Generating synthetics in this way is useful because it provides cost savings compared with on-the-fly wavefield simulations. Synthetics can then be compared with data to determine a best-fitting source. +`GreensTensor` objects +---------------------- -Green's function conventions ----------------------------- - -MTUQ consistently uses an `up-south-east` `convention `_ for internally storing Green's functions (as well as moment tensors and forces). - -A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. +`GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given hypocenter and station. -Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. +Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. -`GreensTensor` objects ----------------------- +Green's function conventions +---------------------------- -`GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given station and origin. +A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. Specifically, MTUQ uses an `up-south-east` `convention `_ for internally storing impulse responses corresponding to individual force couples. (Moment tensors and forces are internally represented using the same `up-south-east` basis.) -Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. +Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. @@ -97,24 +94,25 @@ To generate a full set of Green's functions for a given hypocenter and depth, si .. code :: {event_id}/ - {net}.{sta}.{loc}.Z.Mrr.sac - {net}.{sta}.{loc}.Z.Mtt.sac - {net}.{sta}.{loc}.Z.Mpp.sac - {net}.{sta}.{loc}.Z.Mrt.sac - {net}.{sta}.{loc}.Z.Mrp.sac - {net}.{sta}.{loc}.Z.Mtp.sac - {net}.{sta}.{loc}.R.Mrr.sac - {net}.{sta}.{loc}.R.Mtt.sac - {net}.{sta}.{loc}.R.Mpp.sac - {net}.{sta}.{loc}.R.Mrt.sac - {net}.{sta}.{loc}.R.Mrp.sac - {net}.{sta}.{loc}.R.Mtp.sac - {net}.{sta}.{loc}.T.Mrr.sac - {net}.{sta}.{loc}.T.Mtt.sac - {net}.{sta}.{loc}.T.Mpp.sac - {net}.{sta}.{loc}.T.Mrt.sac - {net}.{sta}.{loc}.T.Mrp.sac - {net}.{sta}.{loc}.T.Mtp.sac + {depth_in_m}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac To open a SPECFEM3D/3D_GLOBE database client: @@ -132,9 +130,9 @@ Once opened, the database client can be used to generate `GreensTensor `_ +`Source-side Green's function details (under construction) `_ -`Receiver-side SPECFEM3D/3D_GLOBE Green's functions (under construction) `_ +`Receiver-side Green's function details (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index 10994583b..c13fdc2c6 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -1,13 +1,12 @@ .. warning:: - This page is currently just a stub. + This page is still under construction. To help improve the + documentation, feel free to submit a pull request. - To fill in missing documentation, feel free to submit a pull request. - -Notes on receiver-side 3D Green's functions -=========================================== +Receiver-side 3D Green's functions +================================== A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index df27bdaf9..b43c17f12 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -1,49 +1,137 @@ .. warning:: - This page is currently just a stub. + This page is still under construction. To help improve the + documentation, feel free to submit a pull request. - To fill in missing documentation, feel free to submit a pull request. +Source-side 3D Green's functions +================================ +Generating source-side 3D Green's functions using SPECFEM3D/3D_GLOBE +-------------------------------------------------------------------- -Notes on source-side 3D Green's functions -========================================= +In principle, any 3D solver can be used to generate Green's functions, as long as the `requirements `_ below are satisfied. + +So far, however, the source-side Green's function machinery has been tested using only SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary. + + +**Generate SAC binary files** + +SAC binary output format is natively supported by SPECFEM3D_GLOBE (in the parameter file, set `OUTPUT_SEISMOS_SAC_BINARY = .true.`). + +Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it becomes necessary to manually convert SPECFEM3D output to SAC binary format. + + +**Convert to SI units** + +MTUQ uses the fully SI convention described below. In contrast, SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions. + + +**Additional amplitude scaling** + +In addition to converting to SI units, it is also necessary to account for any scaling factors in the SPECFEM3D/3D_GLOBE input files. Such scaling factors can enter, for example, through the `M_rr,M_tt,M_pp,M_rt,M_rp,M_tp` values in the moment tensor input file or through the `scaling factor `_ in the force input file. + + +**Rotate to vertical, radial, transverse components** + +Conveniently, SPECFEM3D_GLOBE can be made to automatically rotate output seismograms into vertical, radial and transverse components (set `ROTATE_SEISMOGRAMS_RT = .true.` in the parameter file). + +No modifications are necessary on account of moment tensor basis convention, since MTUQ's `up-south-east` convention matches SPECFEM3D/3D_GLOBE's. + + + + +Requirements for MTUQ source-side 3D Green's functions +------------------------------------------------------ **File format** -- Currently, MTUQ reads source-side 3D Green's functions from SAC files +Individual Green's functions must be written to SAC binary files. + +A total 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion). -- SAC output format is natively supported by SPECFEM3D_GLOBE, but not SPECFEM3D Cartesian (output from the latter can be manually converted) **Units convention** -- MTUQ treats each individual SAC file as the response in meters to a unit (1 Newton-meter) force couple +For a moment tensor inversion, each SAC binary file must give the response in meters to a 1 Newton-meter force couple. + +For a force inversion, each SAC binary file must give the response in meters to a 1 Newton force. + +In both cases, MTUQ uses a fully SI units convention (compare with SPECFEM3D/3D_GLOBE notes, above). -- Users must ensure that source-side Green's functions are normalized according to the above units **Basis convention** -- For moment tensors, MTUQ uses an `up-south-east` convention, identical to one the used by SPECFEM3D_GLOBE +MTUQ uses an `up-south-east` basis convention in which `r` denotes up, `t` denotes south, and `p` denotes east. + +Green's functions must be rotated into into vertical `Z`, radial `R` and transverse `T` components relative to the source-receiver backazimuth. + +Place all seismograms for the same hypocenter in a single directory as follows: + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + ... + +The corresponding convention for force responses is: + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Fr.sac + {net}.{sta}.{loc}.Z.Ft.sac + {net}.{sta}.{loc}.Z.Fp.sac + {net}.{sta}.{loc}.R.Fr.sac + {net}.{sta}.{loc}.R.Ft.sac + {net}.{sta}.{loc}.R.Fp.sac + {net}.{sta}.{loc}.T.Fr.sac + {net}.{sta}.{loc}.T.Fp.sac + {net}.{sta}.{loc}.T.Fp.sac + ... **Origin time convention** -- For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time +For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time. + +MTUQ uses the begin time (`B`) and end time (`E`) headers from the SAC binary files to align the Green's functions relative to centroid origin time. + +Currently, these are the only SAC headers used in reading Green's functions. -- MTUQ uses the `REF_TIME` header from the SPECFEM3D_GLOBE SAC output files, which gives the peak excitation of the source relative to the simulation start time +(Note that different `SAC headers `_ are required in reading `observed data `_.) -- MTUQ ignores the origin time given in the CMTSOLUTION file and `ORIGIN_TIME` header -**Depth searches (experimental)** +Hypocenter searches (experimental) +---------------------------------- -- Only depth searches are possible with source-side 3D Green's functions (no other hypocenter parameters) +Currently, only searches over source depth are possible with source-side 3D Green's functions (no other hypocenter parameters). -- The current `depth search `_ implementation is somewhat crude and experimental; consider local modifcations to suit your needs +The current `depth search `_ implementation is especially crude and experimental (consider local modifications to suit your needs). -- To allow depth searches, create subdirectories for each centroid depth, as below +To allow depth searches, create subdirectories for each centroid depth as follows: .. code :: diff --git a/docs/user_guide/05/gallery_force.rst b/docs/user_guide/05/gallery_force.rst index deaa33576..1383e0757 100644 --- a/docs/user_guide/05/gallery_force.rst +++ b/docs/user_guide/05/gallery_force.rst @@ -33,6 +33,10 @@ For a grid of regulary-spaced forces, `ds` may look something like: * origin_idx (origin_idx) int64 0 +.. note:: + + Force grids are implemented using parameters `F0, phi, h`, which are related to `r, phi, theta` spherical coordinates (`physics convention `_) by `F0 = r`, `phi = phi`, `h = cos(theta)`. In addition, `F0, phi, h` are related to geographic directions by these `formulas `_. + Misfit values """"""""""""" diff --git a/docs/user_guide/05/gallery_mt.rst b/docs/user_guide/05/gallery_mt.rst index 8b57c8d84..ff9880e58 100644 --- a/docs/user_guide/05/gallery_mt.rst +++ b/docs/user_guide/05/gallery_mt.rst @@ -37,6 +37,12 @@ For a grid of regulary-spaced moment tensors, `ds` may look something like: +.. note:: + + Moment tensor grids are implemented using the `rho, v, w, kappa, sigma, h` parameterization of `Tape2012 `_ and `Tape2015 `_, from which `formulas `_ converting to parameterizations can be derived. + + + Misfit values """"""""""""" diff --git a/docs/user_guide/06/trace_attributes.rst b/docs/user_guide/06/trace_attributes.rst index 41014330e..5e267d2f9 100644 --- a/docs/user_guide/06/trace_attributes.rst +++ b/docs/user_guide/06/trace_attributes.rst @@ -2,7 +2,7 @@ Trace attributes ================ -At various points during an inversion, waveform differences, phase shifts, and other values are calculated from observed and synthetic seismic traces. Such `trace attribute` quantities provide important information about how data misfit varies by geographic location and seismic component. +Waveform differences, time-shift corrections, and other values are calculated during an inversion on a trace-by-trace basis. Such `trace attributes` provide important information about how data misfit varies by geographic location and seismic component. Collecting trace attributes From 3d33d54be4732a4b731767719187dd38756a4741 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Mar 2024 15:53:48 -0600 Subject: [PATCH 8/8] Updated check_urls.bash --- tests/check_urls.bash | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/check_urls.bash b/tests/check_urls.bash index 9bef6b62d..f712d96b5 100755 --- a/tests/check_urls.bash +++ b/tests/check_urls.bash @@ -10,8 +10,9 @@ URLS="\ https://www.eas.slu.edu/People/LZhu/home.html\ https://github.com/geodynamics/axisem\ https://github.com/Liang-Ding/seisgen\ - http://ds.iris.edu/ds/products/syngine\ - http://ds.iris.edu/ds/products/syngine/#models\ + https://ds.iris.edu/ds/products/syngine\ + https://ds.iris.edu/ds/products/syngine/#models\ + https://ds.iris.edu/files/sac-manual/manual/file_format.html\ https://instaseis.net\ https://docs.obspy.org/tutorial/index.html\ https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.html\