From 7235a780e2ea9a3456877c1d877708c269feb476 Mon Sep 17 00:00:00 2001 From: Bryant Chow Date: Thu, 1 Dec 2022 16:04:21 -0900 Subject: [PATCH] Inherit pyatoa gathering and I/O routines (#58) * moved all pyatoa moment tensor related functionality into a new utility 'mt' * changed name of util function 'read_synthetics' to 'read_sem' to match Pyatoa original naming (and keep things short for quicker imports) * added in or fixed up I/O routines from Pyatoa which have been fully migrated here * added tests and other routines from pyatoa that need to be refactored and organized * Update README.md * Update README.md * fixed failing tests and added a few tests for new mt capabilities * small formatting changes * fixed typo in util tests, all tests passing --- README.md | 479 +------------------------------ pysep/pysep.py | 7 +- pysep/recsec.py | 14 +- pysep/tests/test_process.py | 1 + pysep/tests/test_utils.py | 46 ++- pysep/utils/fmt.py | 29 ++ pysep/utils/io.py | 253 ++++++++++++----- pysep/utils/mt.py | 543 ++++++++++++++++++++++++++++++++++++ 8 files changed, 821 insertions(+), 551 deletions(-) create mode 100644 pysep/utils/mt.py diff --git a/README.md b/README.md index fb04d4b..784739e 100644 --- a/README.md +++ b/README.md @@ -3,479 +3,12 @@ Python Seismogram Extraction and Processing [![SCOPED](https://img.shields.io/endpoint?url=https://runkit.io/wangyinz/scoped/branches/master/pysep)](https://github.com/SeisSCOPED/container/pkgs/container/pysep) -`PySEP` is a waveform and metadata download tool. The package also contains -`RecSec`, a record section plotting tool. As part of the -[SCOPED toolkit](https://github.com/SeisSCOPED), `PySEP` has been +`PySEP` uses ObsPy routines to request and standardize seismic data and metadata, outputting formatted SAC files that are ready for moment tensor and waveform inversion techniques. The package also contains `RecSec` a **rec**ord **sec**tion plotter for rapid visualization of +waveform data. + +- Documentation can found on the GitHub Wiki Page: https://github.com/adjtomo/pysep/wiki +- Found a bug, requesting a new feature, or wanting to know how to use`PySEP`? [Feel free to open a GitHub issue](https://github.com/adjtomo/pysep/issues). +- As part of the [SCOPED toolkit](https://github.com/SeisSCOPED), `PySEP` has been [containerized](https://github.com/SeisSCOPED/pysep/pkgs/container/pysep) using Docker. -## Introduction - -`PySEP` uses ObsPy tools to request seismic data and metadata, process, -standardize and format waveform data, and write out SAC files with the underlying -motivation of preparing data for moment tensor inversions, although the -written files can be used for other purposes. - -The main processing steps taken by `PySEP` are: - -1. Create or gather event metadata (QuakeML) with user-defined event parameters -2. Gather station metdata (StationXML) in a bounding box surrounding an event, -3. Curtail station list due to distance or azimuth constraints -4. Gather three-component waveform data for chosen station catalog -5. Quality check waveforms: remove gaps, null out missing components, address - clipped amplitudes -6. Preprocess waveforms: detrend, remove response, amplitude scaling, - standardizing time series. -7. Rotate streams to desired components: ZNE, RTZ, UVW (triaxial orthogonal) -8. Append SAC headers with event and station metadata, and TauP arrivals -9. Write per-component SAC files, and StationXML, QuakeML and MSEED files -10. Write CAP (Cut-and-Paste) weight files for moment tensor inversions -11. Write config YAML files which can be used to re-run data gathering/processing -12. Plot a record section and source-receiver map - -`PySEP` also has the capacity to: - -* Interface with a Lawrence Livermore National Laboratory database of nuclear - explosion and earthquake waveforms (see note) -* Allow access to embargoed waveform data and PASSCAL HDF5 (PH5) datasets -* Access pre-defined configuration files for data used in previous studies -* Input custom TauP models for arrival time estimation - -## Installation - -We recommend installing PySEP into a Conda environment. Dependencies are -installed via Conda where possible, with Pip used to install PySEP itself. -This will set two command line tools `pysep` and `recsec` -```bash -conda create -n pysep python=3.10 -conda activate pysep -git clone https://github.com/uafgeotools/pysep.git -cd pysep -conda install --file requirements.txt -pip install -e . -``` - -## Running Tests - -PySEP comes with unit testing which should be run before and after making any -changes to see if your changes have altered the expected code behavior. -```bash -cd pysep/tests -pytest -``` - -## Gallery - -An example record section produced by the `recsec` tool inside PySEP - -![](docs/images/record_section.png) - -An example station map generated from collected metadata during a PySEP run - -![](docs/images/station_map.png) - --------------------------------------------------------------------------------- - -## Command Line Usage - -Normal users should use PySEP as a command line tool. - -To bring up the command line tool help message: - -```bash -pysep -h -``` - -To list out available pre-defined configuration files - -```bash -pysep -l # or pysep --list - --p/--preset -e/--event --p MTUQ2022_workshop -e 2009-04-07T201255_ANCHORAGE.yaml --p MTUQ2022_workshop -e 2014-08-25T161903_ICELAND.yaml --p MTUQ2022_workshop -e 2017-09-03T033001_NORTH_KOREA.yaml --p MTUQ2022_workshop -e 2020-04-04T015318_SOUTHERN_CALIFORNIA.yaml -``` - -To run one of the pre-defined configuration files - -``` bash -pysep -p MTUQ2022_workshop -e 2017-09-03T033001_NORTH_KOREA.yaml -``` - -To create a template configuration file that you must fill out with your own -parameters - -```bash -pysep -W # or pysep --write -``` - -To run this newly created configuration file - -```bash -pysep -c pysep_config.yaml -``` - --------------------------------------------------------------------------------- -### Multiple Event Input - -To use the same configuration file with multiple events, you can use an event -file passed to PySEP through the command line. - -When using this option, the event parameters inside the config file will be -ignored, but all the other parameters will be used to gather data and metadata. - -Event input files should be text files where each row describes one event with -the following parameters as columns: - -> ORIGIN_TIME LONGITUDE LATITUDE DEPTH[KM] MAGNITUDE - -For an example event input file called 'event_input.txt', call structure is: - -```bash -pysep -c pysep_config.yaml -E event_input.txt -``` - --------------------------------------------------------------------------------- -### Legacy Filenaming Schema - -The new version of PySEP uses a file naming schema that is incompatible with -previous versions, which may lead to problems in established workflows. - -To honor the legacy naming schema of PySEP, simply use the ``--legacy_naming`` -argument in the command line. This will change how the event tag is formatted, -how the output directory is structured, and how the output SAC files are named. - -```bash -pysep -c pysep_config.yaml --legacy_naming -``` - --------------------------------------------------------------------------------- -### Output Filename Control - -The event tag used to name the output directory and written SAC files can be set -manually by the user using the ``--event_tag`` argument. If not given, the tag -will default to a string consisting of event origin time and Flinn-Engdahl -region (or just origin time if ``--legacy_naming`` is used). Other output files -such as the config file and ObsPy objects can be set as in the following: - -```bash -pysep -c pysep_config.yaml \ - --overwrite_event_tag event_abc \ - --config_fid event_abc.yaml \ - --stations_fid event_abc_stations.txt \ - --inv_fid event_abc_inv.xml \ - --event_fid event_abc_event.xml \ - --stream_fid event_abc_st.ms -``` - -Please note: the output SAC file names are hardcoded and cannot be changed -by the user. If this is a required feature, please open up a GitHub issue, and -the developers will address this need. - --------------------------------------------------------------------------------- - -## Record Section plotter - -PySEP also comes with a pretty sophisticated record section tool, which plots -seismic data acquired by PySEP. When you have successfully collected your data, -it will reside in the /SAC folder of the PySEP output directory. - - -To see available record section plotting commands - -```bash -recsec -h # RECordSECtion -``` - -To plot the waveform data in a record section with default parameters - -```bash -recsec --pysep_path ./SAC -``` - -To plot a record section with a 7km/s move out, high-pass filtered at 1s -s -```bash -recsec --pysep_path ./SAC --move_out 7 --min_period_s 1 -``` - -To sort your record section by azimuth and not distance (default sorting) - -```bash -recsec --pysep_path ./SAC --sort_by azimuth -``` - -Have a look at the -h/--help message and the docstring at the top of `recsec.py` -for more options. - -### Customizing RecSec figures - -Much of the aesthetic look of RecSec figures is hardcoded, however there are -some keyword arguments that you can provide as flags which may help. The -following parameters correspond to Matplotlib plot adjustments. - -- ytick_fontsize: Fontsize of the Y-axis tick labels -- xtick_fontsize: Fontsize of the X-axis tick labels -- tick_linewidth: Linewidth of axes tick marks -- tick_length: Length of axes tick marks -- label_fontsize: Fontsize of X and Y axis labels -- axis_linewidth: Linewidth of border around figure -- title_fontsize: Fontsize of the title -- xtick_minor: Frequency of minor ticks on the X axis -- xtick_major: Frequency of major ticks on the X axis - -To set one of these parameters, just set as a flag, e.g., - -```bash -recsec --pysep_path ./SAC --xtick_minor 100 --xtick_major 500 -``` - -### Plotting SPECFEM synthetics - -RecSec can plot SPECFEM-generated synthetic seismograms in ASCII format. Here -the domain should be defined by geographic coordinates (latitude/longitude). If -your domain is defined in Cartesian, see below. - -Record sections can be plotted standalone, or alongside observed seismograms -to look at data-synthetic misfit. - -To access metadata, RecSec requires the CMTSOLUTION and STATIONS file that were -used by SPECFEM to generate the synthetics. Based on a standard -SPECFEM3D_Cartesian working directory, plotting synthetics only would have -the following call structure: - -```bash -recsec --syn_path OUTPUT_FILES/ --cmtsolution DATA/CMTSOLUTION --stations DATA/STATIONS -``` - -To compare observed and synthetic data, you would have name the --pysep_path -and tell RecSec to preprocess both data streams identically - -```bash -recsec --pysep_path ./SAC --syn_path OUTPUT_FILES/ --cmtsolution DATA/CMTSOLUTION --stations DATA/STATIONS --preprocess both -``` - -Preprocessing flags such as filtering and move out will be applied to both -observed and synthetic data. - -### Plotting SPECFEM synthetics in Cartesian - -Under the hood, RecSec uses some of ObsPy's geodetic -functions to calculate distances and azimuths. Because of this, RecSec will -fail if coordinates are defined in a Cartesian coordinate system (XYZ), which -may often be the case when working in SPECFEM2D or in a local domain of -SPECFEM3D_Cartesian. - -To circumvent this, RecSec has a flag `--cartesian`, which will swap out the -read functions to work with a Cartesian coordinate system. The call is very -similar to the above: - -For SPECFEM3D_Cartesian this would look like - -```bash -recsec --syn_path OUTPUT_FILES --cmtsolution DATA/CMTSOLUTION --stations DATA/STATIONS --cartesian -``` - -For SPECFEM2D, the source file may not be a CMTSOLUTION. Additionally, the -default seismogram components may not be defined in ZNE - -```bash -recsec --syn_path OUTPUT_FILES --cmtsolution DATA/SOURCE --stations DATA/STATIONS --components Y --cartesian -``` - -### Plotting two sets of synthetics (synsyn) - -It is often useful to compare two sets of synthetic seismograms, where one set -represents 'data', while the other represents synthetics. For example, during -a tomographic inversion, a Target model may be used to generate data. - -RecSec can plot two sets of synthetics in a similar vein as plotting -data and synthetics together. The User only needs to add the `--synsyn` flag -and provide paths to both `--pysep_path` and `--syn_path`. - ->__NOTE__: RecSec makes the assumption that both sets of synthetics share the -> same metadata provided in the `--cmtsolution` and `--stations` flags. - -Let's say you've stored your 'data' in a directory called 'observed/' and your -synthetics in a directory called 'synthetic/' - -```bash -recsec --pysep_path observed/ --syn_path synthetic/ --cmtsolution DATA/CMTSOLUTION --stations DATA/STATIONS --synsyn -``` - -### Scale by: Geometric Spreading - -Rather than normalizing traces, it is possible to scale peak amplitudes using -a simple geometric spreading equation, which takes into account amplitude -loss due to propagation distance. See the docstring of -`recsec.RecordSection.calculate\_geometric\_spreading()` for more details on -how this is implemented. - -To apply geometric spreading to data collected by PySEP for your record section: - -```bash -recsec --pysep_path ./SAC --scale_by geometric_spreading - --geometric_spreading_exclude STA1 STA2 STA3 - --geometric_spreading_save ./geometric_spread.png -``` - -Where `--geometric_spreading_exclude` is a list of stations that should *not* be -included in the spreading equation (here, e.g., STA1 through STA3). - -Note that this option will generate a scatterplot showing the line fit to -the geometric spreading function saved to with figure name -'./geometric_spreading.png' - -Have a look at the RecSec init docstring for 'scale\_by' to see related -parameters (which start with *geometric_spreading_*) and how they should be used. - --------------------------------------------------------------------------------- - -## Scripting PySEP - -More advanced users can use PySEP as a scripting tool rather than a command -line tool. - -All of the gathered data/metadata are saved as attributes of the Pysep class -with typical ObsPy naming schema - -```python ->>> from pysep import Pysep ->>> sep = Pysep(config_file='pysep_config.yaml') ->>> sep.run() ->>> print(sep.st[0].stats.sac) ->>> sep.inv ->>> sep.event -``` - -Although not the preferred method of interacting with PySEP, you can forgo the -config file and pass parameters directly to the instantiation of the PySEP -class, making PySEP a bit more flexible. - -```python ->>> from pysep import Pysep ->>> sep = Pysep(origin_time="2000-01-01T00:00:00", event_latitude=64.8596, - event_longitude=-147.8498, event_depth_km=15., .... - ) -``` - -Check out the Pysep.run() function for other API options for using PySEP. - -### Scripting RecSec - -The RECord SECtion tool can also be scripted. It simply requires an ObsPy Stream -object as input. Tunable parameters can be fed in as input variables. - -```python ->>> from obspy import read ->>> from pysep.recsec import plotw_rs ->>> st = read() ->>> plotw_rs(st=st, sort_by="distance") -``` - -## Miscellaneous Functionality - -### Append SAC headers to existing Streams - -To append SAC headers to your own seismic data, you can directly use the -`PySEP` utility functions `append_sac_headers` and -`format_sac_header_w_taup_traveltimes` - -```python ->>> from pysep.utils.cap_sac import append_sac_headers, format_sac_header_w_taup_traveltimes ->>> from obspy import read, read_events, read_inventory ->>> st = read() ->>> inv = read_inventory() ->>> event = read_events()[0] ->>> st = append_sac_headers(st=st, inv=inv, event=event) ->>> st = format_sac_header_w_taup_traveltimes(st=st, model="ak135") -``` - -### Converting SPECFEM-generated synthetics - -PySEP contains a utility function `read_synthetics` to read in -SPECFEM-generated synthetics with appropriately crafted SAC headers. -Given a standard SPECFEM3D working directory, reading in SPECFEM synthetics -and saving them as SAC files might look something like: - -```python ->>> from pysep.utils.io import read_synthetics ->>> st = read_synthetics(fid="OUTPUT_FILES/NZ.BFZ.BXE.semd", - cmtsolution="DATA/CMTSOLUTION", - stations="DATA/STATIONS") ->>> for tr in st: ->>> st.write(f"{tr.get_id()}.sac", format="SAC") -``` - -### Pointing PySEP to custom, local databases -Data are often stored in custom databases that we cannot predict the -structure of. To point PySEP at your local databases, the simplest method would -be to find a way to read your data and metadata into ObsPy objects, which -you can then feed into the PySEP machinery. - -```python ->>> from pysep import Pysep ->>> from obspy import read, read_events, read_inventory ->>> st = read() ->>> inv = read_inventory() ->>> event = read_events()[0] ->>> sep = Pysep(st=st, inv=inv, event=event, config_file="config.yaml") -``` - - - --------------------------------------------------------------------------------- - -## Running PySEP on UAF Chinook - -Chinook is University of Alaska Fairbanks' (UAF) high performance computer. -We can run PySEP on Chinook using Docker containers through -Singularity/Apptainer. - -PySEP has been containerized directly, and any changes pushed to the repo will -trigger the container to rebuild, keeping everything up-to-date. -The following code snippet downloads the correct -[SCOPED container](https://github.com/SeisSCOPED/pysep/pkgs/container/pysep). -and runs the PySEP help message on Chinook. - -```bash -module load singularity -singularity pull ghcr.io/seisscoped/pysep:centos7 -singularity exec -c pysep_centos7.sif pysep -h -``` - -To run a data download we will need to mount the local filesystem into the -container using the ``--bind`` command. Using the Anchorage example event: - -```bash -singularity exec -c --bind $(pwd):/home1 pysep_centos7.sif \ - bash -c "cd /home1/; pysep -p mtuq_workshop_2022 -e 2009-04-07T201255_ANCHORAGE.yaml" -``` - -In the above example, the `-c/--contain` flag preserves the internal container -filesystem, the `-B/--bind` flag binds the current working directory on Chinook -(i.e., pwd)to a directory called */home1* within the container, and then the -`bash -c` command changes to this */home1* directory and runs PySEP. Files are -subsequently saved to the local filesystem in our current working directory. - -`RecSec`, the record section plotting tool, can be run from the command line -using a similar format. With the Anchorange example files we just generated: - -```bash -cd 2009-04-07T201255_SOUTHERN_ALASKA/ -singularity exec -c --bind $(pwd):/home1 ../pysep_centos7.sif \ - bash -c "cd /home1; recsec --pysep_path SAC/ --min_period 10 --save record_section_tmin10.png" -``` - -## LLNL Note - -`PySEP` interfaces with the databases of: - -* W. Walter et al. (2006) - An assembled western United States dataset for regional seismic analysis - ISSO 9660 CD, LLNL release UCRL-MI-222502 - - Download link: https://ds.iris.edu/mda/18-001 - diff --git a/pysep/pysep.py b/pysep/pysep.py index d7429c9..e8d1dcc 100644 --- a/pysep/pysep.py +++ b/pysep/pysep.py @@ -30,7 +30,7 @@ quality_check_waveforms_before_processing, quality_check_waveforms_after_processing) from pysep.utils.fmt import format_event_tag, format_event_tag_legacy, get_codes -from pysep.utils.io import read_yaml, read_event_file, write_stations_file +from pysep.utils.io import read_yaml, read_event_file, write_pysep_stations_file from pysep.utils.llnl import scale_llnl_waveform_amplitudes from pysep.utils.process import (merge_and_trim_start_end_times, resample_data, format_streams_for_rotation, rotate_to_uvw, @@ -906,7 +906,8 @@ def rotate_streams(self): for sta in stations: _st = st_rtz.select(station=sta) _st.rotate(method="NE->RT") # in place rot. - logger.debug(f"{sta}: BAz={_st[0].stats.back_azimuth}") + if hasattr(_st[0].stats, "back_azimuth"): + logger.debug(f"{sta}: BAz={_st[0].stats.back_azimuth}") st_out += st_rtz st_out = format_sac_headers_post_rotation(st_out) @@ -994,7 +995,7 @@ def write(self, write_files=None, _return_filenames=False, stations_fid or "stations_list.txt") logger.info("writing stations file") logger.debug(fid) - write_stations_file(self.inv, self.event, fid) + write_pysep_stations_file(self.inv, self.event, fid) if "inv" in write_files or "all" in write_files: fid = os.path.join(self.output_dir, inv_fid or f"inv.xml") diff --git a/pysep/recsec.py b/pysep/recsec.py index f073b22..a5ff5b5 100755 --- a/pysep/recsec.py +++ b/pysep/recsec.py @@ -97,7 +97,7 @@ from pysep import logger from pysep.utils.cap_sac import origin_time_from_sac_header -from pysep.utils.io import read_synthetics, read_synthetics_cartesian +from pysep.utils.io import read_sem, read_sem_cartesian from pysep.utils.curtail import subset_streams from pysep.utils.plot import plot_geometric_spreading, set_plot_aesthetic @@ -491,15 +491,13 @@ def _generate_synthetic_stream(self, syn_path, source, stations, for fid in fids: logger.debug(fid) if not cartesian: - st_syn += read_synthetics(fid=fid, - cmtsolution=source, - stations=stations) + st_syn += read_sem(fid=fid, source=source, + stations=stations) else: # If we are using SPECFEM2D synthetics, trying to read # the SOURCE file will - st_syn += read_synthetics_cartesian(fid=fid, - source=source, - stations=stations) + st_syn += read_sem_cartesian(fid=fid, source=source, + stations=stations) return st_syn def check_parameters(self): @@ -1481,7 +1479,7 @@ def _plot_trace(self, idx, y_index, choice="st"): x -= self.zero_pad_s[0] # index 0 is start, index 1 is end # Synthetics will already have a time offset stat from the - # 'read_synthetics' function which grabs T0 value from ASCII + # 'read_sem' function which grabs T0 value from ASCII # Data will have a time offset relative to event origin time if hasattr(tr.stats, "time_offset"): x += tr.stats.time_offset diff --git a/pysep/tests/test_process.py b/pysep/tests/test_process.py index e53c00a..d6e9315 100644 --- a/pysep/tests/test_process.py +++ b/pysep/tests/test_process.py @@ -143,6 +143,7 @@ def test_rotate_streams_fail(test_st, test_inv, test_event): sep.event = test_event sep.st = sep.preprocess() # make sure that streams are formatted correctly + # No back aizmuth attribute found with pytest.raises(TypeError): sep.rotate_streams() diff --git a/pysep/tests/test_utils.py b/pysep/tests/test_utils.py index 64dade0..54fc1df 100644 --- a/pysep/tests/test_utils.py +++ b/pysep/tests/test_utils.py @@ -18,7 +18,9 @@ from pysep.utils.fmt import format_event_tag, format_event_tag_legacy from pysep.utils.process import estimate_prefilter_corners from pysep.utils.plot import plot_source_receiver_map -from pysep.utils.io import read_synthetics +from pysep.utils.io import read_sem +from pysep.utils.mt import get_gcmt_moment_tensor, get_usgs_moment_tensor + @pytest.fixture @@ -137,7 +139,7 @@ def test_plot_map(test_event, test_inv): show=False) -def test_read_synthetics(): +def test_read_sem(): """ Test reading SPECFEM-generated synthetics in as SAC formatted files """ @@ -147,8 +149,8 @@ def test_read_synthetics(): st = Stream() for test_synthetic in test_synthetics: - st += read_synthetics(fid=test_synthetic, cmtsolution=test_cmtsolution, - stations=test_stations) + st += read_sem(fid=test_synthetic, source=test_cmtsolution, + stations=test_stations) assert(st[0].stats.sac.evla == -40.5405) @@ -195,3 +197,39 @@ def test_subset_streams(test_st): # all station ids should be the same for tr_out, tr_smaller_out in zip(test_st_out, test_st_smaller_out): assert(tr_out.get_id() == tr_smaller_out.get_id()) + + +def test_get_usgs_moment_tensor(): + """ + Just ensure that getting via USGS works as intended using an example event + """ + event = test_get_gcmt_moment_tensor() + del event.focal_mechanisms + + cat = get_usgs_moment_tensor(event=event) + assert(len(cat) == 1) + event = cat[0] + assert hasattr(event, "focal_mechanisms") + + m_rr = event.preferred_focal_mechanism().moment_tensor.tensor.m_rr + assert(pytest.approx(m_rr, 1E-20) == 4.81E20) + + +@pytest.mark.skip(reason="test already called by 'test_get_usgs_moment_tensor'") +def test_get_gcmt_moment_tensor(): + """ + Just ensure that getting via GCMT works as intended using an example event + """ + # Kaikoura Earthquake + origintime = "2016-11-13T11:02:00" + magnitude = 7.8 + + cat = get_gcmt_moment_tensor(event=None, origintime=origintime, + magnitude=magnitude) + assert(len(cat) == 1) + event = cat[0] + assert hasattr(event, "focal_mechanisms") + m_rr = event.preferred_focal_mechanism().moment_tensor.tensor.m_rr + assert(pytest.approx(m_rr, 1E-20) == 3.56E20) + + return event diff --git a/pysep/utils/fmt.py b/pysep/utils/fmt.py index f2c1abd..62d7380 100644 --- a/pysep/utils/fmt.py +++ b/pysep/utils/fmt.py @@ -4,6 +4,35 @@ from obspy.clients.iris import Client +def channel_code(dt): + """ + Specfem outputs seismograms with channel band codes set by IRIS. Instrument + codes are always X for synthetics, but band code will vary with the sampling + rate of the data, return the correct code given a sampling rate. + Taken from Appenix B of the Specfem3D cartesian manual (June 15, 2018) + + :type dt: float + :param dt: sampling rate of the data in seconds + :rtype: str + :return: band code as specified by Iris + :raises KeyError: when dt is specified incorrectly + """ + if dt >= 1: + return "L" # long period + elif 0.1 < dt < 1: + return "M" # mid period + elif 0.0125 < dt <= 0.1: + return "B" # broad band + elif 0.001 <= dt <= 0.0125: + return "H" # high broad band + elif 0.004 <= dt < 0.001: + return "C" + elif 0.001 <= dt < 0.0002: + return "F" + else: + raise KeyError("Channel code does not exist for this value of 'dt'") + + def get_codes(st=None, choice=None, suffix=None, up_to=True): """ Get station codes from the internal stream attribute, where station diff --git a/pysep/utils/io.py b/pysep/utils/io.py index 870ba67..00ff4c5 100644 --- a/pysep/utils/io.py +++ b/pysep/utils/io.py @@ -12,14 +12,51 @@ from obspy.geodetics import gps2dist_azimuth from obspy.core.event import Event, Origin, Magnitude -from pysep.utils.fmt import format_event_tag_legacy +from pysep import logger +from pysep.utils.mt import moment_magnitude, seismic_moment, Source +from pysep.utils.fmt import format_event_tag_legacy, channel_code from pysep.utils.cap_sac import append_sac_headers -def read_synthetics(fid, cmtsolution, stations, location="", precision=4): +def read_yaml(fid): + """ + Read a YAML file and return a dictionary + + :type fid: str + :param fid: YAML file to read from + :rtype: dict + :return: YAML keys and variables in a dictionary """ - Specfem3D outputs seismograms to ASCII (.sem?) files. Converts SPECFEM - .sem? files into Stream objects with the correct header + # work around PyYAML bugs + yaml.SafeLoader.add_implicit_resolver( + u'tag:yaml.org,2002:float', + re.compile(u'''^(?: + [-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)? + |[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+) + |\\.[0-9_]+(?:[eE][-+][0-9]+)? + |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\\.[0-9_]* + |[-+]?\\.(?:inf|Inf|INF) + |\\.(?:nan|NaN|NAN))$''', re.X), + list(u'-+0123456789.')) + + with open(fid, 'r') as f: + config = yaml.safe_load(f) + + # Replace 'None' and 'inf' values to match expectations + for key, val in config.items(): + if val == "None": + config[key] = None + if val == "inf": + config[key] = np.inf + + return config + + +def read_sem(fid, origintime=None, source=None, stations=None, location="", + precision=4): + """ + Specfem3D outputs seismograms to ASCII (.sem? or .sem.ascii) files. + Converts SPECFEM synthetics into ObsPy Stream objects with the correct header information. .. note :: @@ -27,14 +64,20 @@ def read_synthetics(fid, cmtsolution, stations, location="", precision=4): :type fid: str :param fid: path of the given ascii file - :type cmtsolution: str - :param cmtsolution: CMTSOLUTION file defining the event which generated - the synthetics. Used to grab event information. + :type origintime: obspy.UTCDateTime + :param origintime: UTCDatetime object for the origintime of the event + :type source: str + :param source: optional SPECFEM source file (e.g., CMTSOLUTION, SOURCE) + defining the event which generated the synthetics. Used to grab event + information and append as SAC headers to the ObsPy Stream :type stations: str - :param stations: STATIONS file defining the station locations for the - SPECFEM generated synthetics + :param stations: optional STATIONS file defining the station locations for + the SPECFEM generated synthetics, used to generate SAC headers :type location: str :param location: location value for a given station/component + :type precision: int + :param precision: dt precision determined by differencing two + adjancent time steps in the underlying ascii text file. :rtype st: obspy.Stream.stream :return st: stream containing header and data info taken from ascii file """ @@ -69,13 +112,18 @@ def read_synthetics(fid, cmtsolution, stations, location="", precision=4): delta = round(times[1] - times[0], precision) # Get metadata information from CMTSOLUTION and STATIONS files - event = read_events(cmtsolution, format="CMTSOLUTION")[0] - inv = read_stations(stations) - - starttime = event.preferred_origin().time + if origintime is None and source is None: + logger.warning("no `origintime` or `event` given, setting dummy " + "starttime: 1970-01-01T00:00:00") + origintime = UTCDateTime("1970-01-01T00:00:00") + elif source: + # !!! Add logic here to try different read options + event = read_events(source, format="CMTSOLUTION")[0] + origintime = event.preferred_origin().time + logger.info(f"reading origintime from event: {origintime}") # Honor that Specfem doesn't start exactly on 0 - starttime += times[0] + origintime += times[0] # Write out the header information try: @@ -86,19 +134,27 @@ def read_synthetics(fid, cmtsolution, stations, location="", precision=4): net, sta, cha, fmt, suffix = os.path.basename(fid).split(".") stats = {"network": net, "station": sta, "location": location, - "channel": cha, "starttime": starttime, "npts": len(data), + "channel": cha, "starttime": origintime, "npts": len(data), "delta": delta, "mseed": {"dataquality": 'D'}, "format": fmt } st = Stream([Trace(data=data, header=stats)]) - st = append_sac_headers(st, event, inv) + + if source and stations: + try: + inv = read_stations(stations) + st = append_sac_headers(st, event, inv) + # Broad catch here as this is an optional step that might not always + # work or be possible + except Exception as e: + logger.warning(f"could not append SAC header to trace because {e}") return st -def read_synthetics_cartesian(fid, source, stations, location="", precision=4): +def read_sem_cartesian(fid, source, stations, location="", precision=4): """ Specfem2D and Specfem3D may have domains defined in a Cartesian coordinate - system. Because of this, the read_synthetics() function will fail because + system. Because of this, the read_sem() function will fail because the intermediate ObsPy objects and functions expect geographic coordinates. This function bypasses these checks with some barebones objects which mimic their behavior. Only used for RecSec to plot record sections. @@ -289,27 +345,33 @@ def _get_resource_id(name, res_type, tag=None): def read_specfem2d_source(path_to_source, origin_time=None): """ Create a barebones ObsPy Event object from a SPECFEM2D Source file, which - contains information on location. The remainder are left as dummy values - because this will only be required by RecSec to access source location. + only contains information required by Pyatoa. - .. note:: - This is modified from Pyatoa.utils.read.read_specfem2d_source + Only requires access to: event.preferred_origin(), + event.preferred_magnitude() and event.preferred_moment_tensor(). + Moment tensor is wrapped in try-except so we only need origin and magnitude. + + Modified from: + https://docs.obspy.org/master/_modules/obspy/io/cmtsolution/ + core.html#_internal_read_cmtsolution .. note:: Source files do not provide origin times so we just provide an arbitrary value but allow user to set time """ + def _get_resource_id(name, res_type, tag=None): - """ - Helper function to create consistent resource ids. From ObsPy - """ + """Helper function to create consistent resource ids. From ObsPy""" res_id = f"smi:local/source/{name:s}/{res_type:s}" if tag is not None: res_id += "#" + tag return res_id + # First set dummy origin time if origin_time is None: - origin_time = "2000-01-01T00:00:00" + origin_time = "1970-01-01T00:00:00" + logger.warning("no origin time set for SPECFEM2D source, setting " + f"dummy value: {origin_time}") with open(path_to_source, "r") as f: lines = f.readlines() @@ -332,9 +394,21 @@ def _get_resource_id(name, res_type, tag=None): latitude=source_dict["xs"], depth=source_dict["zs"] ) + # Attempt to calculate the moment magnitude from the moment tensor + # components. This might fail because the values don't make sense or are + # not provided + try: + moment = seismic_moment(mt=[float(source_dict["Mxx"]), + float(source_dict["Mzz"]), + float(source_dict["Mxz"]) + ]) + moment_mag = moment_magnitude(moment=moment) + except Exception as e: + moment_mag = None + magnitude = Magnitude( resource_id=_get_resource_id(source_name, "magnitude"), - mag=None, magnitude_type="Mw", origin_id=origin.resource_id.id + mag=moment_mag, magnitude_type="Mw", origin_id=origin.resource_id.id ) event = Event(resource_id=_get_resource_id(name=source_name, @@ -348,16 +422,41 @@ def _get_resource_id(name, res_type, tag=None): return event +def read_forcesolution(path_to_source, default_time="2000-01-01T00:00:00"): + """ + Create a barebones Pyatoa Source object from a FORCESOLUTION Source file, + which mimics the behavior of the more complex ObsPy Event object + """ + with open(path_to_source, "r") as f: + lines = f.readlines() + + # Place values from file into dict + source_dict = {} + for line in lines[:]: + if ":" not in line: + continue + key, val = line.strip().split(":") + source_dict[key] = val + + # First line has the source name + source_dict["source_id"] = lines[0].strip().split()[-1] + + event = Source( + resource_id=f"pysep:source/{source_dict['source_id']}", + origin_time=default_time, latitude=source_dict["latorUTM"], + longitude=source_dict["longorUTM"], depth=source_dict["depth"] + ) + + return event + + def read_stations(path_to_stations): """ - Convert a Specfem3D STATIONS file into an ObsPy Inventory object. + Convert a SPECFEM STATIONS file into an ObsPy Inventory object. Specfem3D STATION files contain no channel or location information, so the inventory can only go down to the station level. - .. note:: - This is copied verbatim from Pyatoa.utils.read.read_stations() - .. note:: This assumes a row structure for the station file is STA, NET, LAT [deg], LON [deg], ELEVATION [m], BURIAL [m] @@ -436,43 +535,34 @@ def read_event_file(fid): return list_out -def read_yaml(fid): +def write_sem(st, unit, path="./", time_offset=0): """ - Read a YAML file and return a dictionary - - :type fid: str - :param fid: YAML file to read from - :rtype: dict - :return: YAML keys and variables in a dictionary + Write an ObsPy Stream in the two-column ASCII format that Specfem uses + + :type st: obspy.core.stream.Stream + :param st: stream containing synthetic data to be written + :type unit: str + :param unit: the units of the synthetic data, used for file extension, must + be 'd', 'v', 'a' for displacement, velocity, acceleration, resp. + :type path: str + :param path: path to save data to, defaults to cwd + :type time_offset: float + :param time_offset: optional argument to offset the time array. Sign matters + e.g. time_offset=-20 means t0=-20 """ - # work around PyYAML bugs - yaml.SafeLoader.add_implicit_resolver( - u'tag:yaml.org,2002:float', - re.compile(u'''^(?: - [-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)? - |[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+) - |\\.[0-9_]+(?:[eE][-+][0-9]+)? - |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\\.[0-9_]* - |[-+]?\\.(?:inf|Inf|INF) - |\\.(?:nan|NaN|NAN))$''', re.X), - list(u'-+0123456789.')) - - with open(fid, 'r') as f: - config = yaml.safe_load(f) - - # Replace 'None' and 'inf' values to match expectations - for key, val in config.items(): - if val == "None": - config[key] = None - if val == "inf": - config[key] = np.inf - - return config + assert(unit.lower() in ["d", "v", "a"]), "'unit' must be 'd', 'v' or 'a'" + for tr in st: + s = tr.stats + fid = f"{s.network}.{s.station}.{channel_code(s.delta)}X{s.channel[-1]}" + fid = os.path.join(path, f"{fid}.sem{unit.lower()}") + data = np.vstack((tr.times() + time_offset, tr.data)).T + np.savetxt(fid, data, ["%13.7f", "%17.7f"]) -def write_stations_file(inv, event, fid="./stations_list.txt"): +def write_pysep_stations_file(inv, event, fid="./stations_list.txt"): """ - Write a list of station codes, distances, etc. + Write a list of station codes, distances, etc. useful for understanding + characterization of all collected stations :type event: obspy.core.event.Event :param event: optional user-provided event object which will force a @@ -498,3 +588,40 @@ def write_stations_file(inv, event, fid="./stations_list.txt"): f.write(f"{sta.code:<6} {net.code:<2} " f"{sta.latitude:9.4f} {sta.longitude:9.4f} " f"{dist_km:8.3f} {az:6.2f}\n") + + +def write_specfem_stations_file(inv, fid="./STATIONS", elevation=False, + burial=0.): + """ + Write a SPECFEM3D STATIONS file given an ObsPy inventory object + + .. note:: + If topography is implemented in your mesh, elevation values should be + set to 0 which means 'surface' in SPECFEM. + + .. note:: + burial information is not contained in the ObsPy inventory so it is + always set to a constant value, which can be given as input. default 0 + + :type inv: obspy.core.inventory.Inventory + :param inv: Inventory object with station locations to write + :type fid: str + :param fid: path and file id to save the STATIONS file to. + :type elevation: bool + :param elevation: if True, sets the actual elevation value from the + inventory. if False, sets elevation to 0. + :type burial: float + :param burial: the constant value to set burial values to. defaults to 0. + """ + with open(fid, "w") as f: + for net in inv: + for sta in net: + lat = sta.latitude + lon = sta.longitude + if elevation: + elv = sta.elevation + else: + elv = 0. + + f.write(f"{sta.code:>6}{net.code:>6}{lat:12.4f}{lon:12.4f}" + f"{elv:7.1f}{burial:7.1f}\n") \ No newline at end of file diff --git a/pysep/utils/mt.py b/pysep/utils/mt.py new file mode 100644 index 0000000..fa20365 --- /dev/null +++ b/pysep/utils/mt.py @@ -0,0 +1,543 @@ +""" +Moment Tensor related functions for grabbing moment tensors, appending them +to existing event objects, and writing them out in specific formats such as the +CMTSOLUTION format required by SPECFEM +""" +import csv +import requests +import numpy as np +from obspy.core.event import source, Event +from obspy.core.event.base import Comment +from obspy.core.event.source import Tensor +from urllib.error import HTTPError +from obspy import Catalog, UTCDateTime, read_events +from obspy.clients.fdsn import Client + +from pysep import logger + + +def seismic_moment(mt): + """ + Return the seismic moment based on a moment tensor. + Can take a list of tensor components, or a Tensor object from ObsPy. + + :type mt: list of floats or obspy.core.event.source.Tensor + :param mt: the components of the moment tensor M_ij + :rtype: float + :return: the seismic moment, in units of N*m + """ + if isinstance(mt, Tensor): + # Little one liner to spit out moment tensor components into a list + mt_temp = [getattr(mt, key) for key in mt.keys() + if not key.endswith("errors")] + assert (len(mt_temp) == 6), "Moment tensor should have 6 components" + mt = mt_temp + return 1 / np.sqrt(2) * np.sqrt(sum([_ ** 2 for _ in mt])) + + +def moment_magnitude(moment, c=10.7): + """ + Return the moment magitude, M_w, based on a seismic moment. Equation from + Hanks & Kanamori (1979) + + :type c: float + :param c: correction factor for conversion, 10.7 for units of N*m, + 16.1 for units of dyne*cm + :type moment: float + :param moment: the seismic moment, in units of N*m + :rtype: float + :return: moment magnitude, M_w + """ + return 2 / 3 * np.log10(moment) - c + + +def half_duration_from_m0(moment): + """ + Empirical formula for half duration used by Harvard CMT, stated in + Dahlen and Tromp (1998, p.178). + + :type moment: float + :param moment: seismic moment in N*m + :rtype: float + :return: empirically scaled half duration in unit seconds + """ + return 2.4E-6 * moment**(1/3) + + +def mt_transform(mt, method): + """ + Transform moment tensor between XYZ and RTP coordinates. Used primarily + to transform GeoNet (John Ristau) moment tensors into the correct coordinate + system for use in SPECFEM + + Based on Equation from 'Aki and Richards: Quantitative Seismology' book + + .. note:: + Acceptable formats for the parameter `mt`: + + 1) [m11, m22, m33, m12, m13, m23] + 2) [mxx, myy, mzz, mxy, mxz, myz] + 3) [mrr, mtt, mpp, mrt, mrp, mtp] + + :type mt: dict + :param mt: moment tensor in format above + :type method: str + :param method: type of conversion, "rtp2xyz" or "xyz2rtp" + :rtype: dict + :return: converted moment tensor dictionary + """ + assert(method in ["xyz2rtp", "rtp2xyz"]), \ + "method must be 'xyz2rtp' or 'rtp2xyz'" + + if method == "xyz2rtp": + if "m_xx" not in mt.keys(): + print("for xyz2rtp, dict must have keys in xyz") + m_rr = mt["m_zz"] + m_tt = mt["m_xx"] + m_pp = mt["m_yy"] + m_rt = mt["m_xz"] + m_rp = -1 * mt["m_yz"] + m_tp = -1 * mt["m_xy"] + return {"m_rr": m_rr, "m_tt": m_tt, "m_pp": m_pp, "m_rt": m_rt, + "m_rp": m_rp, "m_tp": m_tp} + elif method == "rtp2xyz": + if "m_tt" not in mt.keys(): + print("for rtp2xyz, dict must have keys in rtp") + m_xx = mt["m_tt"] + m_yy = mt["m_pp"] + m_zz = mt["m_rr"] + m_xy = -1 * mt["m_tp"] + m_xz = mt["m_rt"] + m_yz = -1 * mt["m_rp"] + return {"m_xx": m_xx, "m_yy": m_yy, "m_zz": m_zz, "m_xy": m_xy, + "m_xz": m_xz, "m_yz": m_yz} + + +def get_gcmt_moment_tensor(event=None, origintime=None, magnitude=None, + time_wiggle_sec=120, magnitude_wiggle=0.5): + """ + Query online GCMT moment tensor catalog via URL access for moment tensor + components of a given event. Searches based on origin time and magnitude + of an event with a given amount of wiggle room for catalog mismatch of + origin time and magnitude. + + .. note:: + input is either `event` OR `origintime` AND `magnitude` + + :type event: obspy.core.event.Event + :param event: Event to use to query for moment tensor + :type origintime: UTCDateTime or str + :param origintime: event origin time + :type magnitude: float + :param magnitude: centroid moment magnitude for event lookup + :type time_wiggle_sec: int + :param time_wiggle_sec: padding on catalog filtering criteria realted to + event origin time + :type magnitude_wiggle: float + :param magnitude_wiggle: padding on catalog filter for magnitude + :rtype: obspy.core.event.Event + :return: event object for given earthquake + """ + if event is None: + assert(origintime is not None and magnitude is not None), ( + "GCMT moment tensor query requires `event` or `origintime` " + "and `magnitude" + ) + origintime = UTCDateTime(origintime) + else: + origintime = event.preferred_origin().time + magnitude = event.preferred_magnitude().mag + + # Determine filename using datetime properties + month = origintime.strftime('%b').lower() # e.g. 'jul' + year_short = origintime.strftime('%y') # e.g. '19' + year_long = origintime.strftime('%Y') # e.g. '2019' + + fid = f"{month}{year_short}.ndk" + try: + cat = read_events( + "https://www.ldeo.columbia.edu/~gcmt/projects/CMT/" + f"catalog/NEW_MONTHLY/{year_long}/{fid}" + ) + except HTTPError: + cat = read_events( + "http://www.ldeo.columbia.edu/~gcmt/projects/CMT/" + "catalog/NEW_QUICK/qcmt.ndk" + ) + # GCMT catalogs contain all events for a span of time + # filter catalogs using ObsPy to find events with our specifications. + # Magnitudes and origintimes are not always in agreement between agencies + # So allow for some wiggle room + cat_filt = cat.filter(f"time > {str(origintime - time_wiggle_sec)}", + f"time < {str(origintime + time_wiggle_sec)}", + f"magnitude >= {magnitude - magnitude_wiggle}", + f"magnitude <= {magnitude + magnitude_wiggle}", + ) + + return cat_filt + + +def get_usgs_moment_tensor(event, time_wiggle_sec=120., magnitude_wiggle=.5, + latitude_wiggle_deg=1., longitude_wiggle_deg=1., + depth_wiggle_km=5., **kwargs): + """ + Query FDSN webservices USGS client for moment tensors using the current + event definition, which may or may not have been collected via USGS. + + Kwargs passed to Client.get_events() for additional event constraint pars + + :type event: obspy.core.event.Event + :param event: Event to use to query for moment tensor + :type time_wiggle_sec: float + :param time_wiggle_sec: padding on catalog filtering criteria realted to + event origin time + :type magnitude_wiggle: float + :param magnitude_wiggle: +/- padding on magnitude search + :type latitude_wiggle_deg: float + :param latitude_wiggle_deg: +/- padding on latitude search + :type longitude_wiggle_deg: float + :param longitude_wiggle_deg: +/- padding on longitude search + :type depth_wiggle_km: float + :param depth_wiggle_km: +/- padding on depth search + :rtype: obspy.core.event.Event + :return: event object for given earthquake + """ + c = Client("USGS") + origintime = event.preferred_origin().time + magnitude = event.preferred_magnitude().mag + latitude = event.preferred_origin().latitude + longitude = event.preferred_origin().longitude + depth = event.preferred_origin().depth * 1E-3 + + # Assuming that time, magnitude, and hypocenter are enough to uniquely + # identify a given earthquake + try: + cat = c.get_events(starttime=origintime - time_wiggle_sec, + endtime=origintime + time_wiggle_sec, + minmagnitude=magnitude - magnitude_wiggle, + maxmagnitude=magnitude + magnitude_wiggle, + mindepth=depth - depth_wiggle_km, + maxdepth=depth + depth_wiggle_km, + minlatitude=latitude - latitude_wiggle_deg, + maxlatitude=latitude + latitude_wiggle_deg, + minlongitude=longitude - longitude_wiggle_deg, + maxlongitude=longitude + longitude_wiggle_deg, + includeallorigins=True, **kwargs + ) + # Broad failure criteria but these are usually FDSNExceptions from ObsPy + except Exception as e: + logger.warning(e) + cat = None + return cat + + +def get_geonet_mt(event_id, units, csv_fid=None): + """ + Focal mechanisms created by John Ristau are written to a .csv file + located on Github. This function will append information from the .csv file + onto the Obspy event object so that all the information can be located in a + single object + + :type event_id: str + :param event_id: unique event identifier + :type units: str + :param units: output units of the focal mechanism, either: + 'dynecm': for dyne*cm or + 'nm': for Newton*meter + :type csv_fid: str + :param csv_fid: optional local path to .csv file containing Ristau catalog + of moment tensors. If not given, will search a default online URL + where the catalog is assumed to be stored + :rtype focal_mechanism: obspy.core.event.FocalMechanism + :return focal_mechanism: generated focal mechanism + """ + assert(units in ["dynecm", "nm"]), "units must be 'dynecm' or 'nm'" + + mtlist = query_geonet_mt_catalog(event_id, csv_fid=csv_fid) + + # Match the identifier with Goenet + id_template = f"smi:local/geonetcsv/{mtlist['PublicID']}/{{}}" + + # Generate the Nodal Plane objects containing strike-dip-rake + nodal_plane_1 = source.NodalPlane( + strike=mtlist['strike1'], dip=mtlist['dip1'], rake=mtlist['rake1'] + ) + nodal_plane_2 = source.NodalPlane( + strike=mtlist['strike2'], dip=mtlist['dip2'], rake=mtlist['rake2'] + ) + nodal_planes = source.NodalPlanes( + nodal_plane_1, nodal_plane_2, preferred_plane=1 + ) + + # Create the Principal Axes as Axis objects + tension_axis = source.Axis( + azimuth=mtlist['Taz'], plunge=mtlist['Tpl'], length=mtlist['Tva'] + ) + null_axis = source.Axis( + azimuth=mtlist['Naz'], plunge=mtlist['Npl'], length=mtlist['Nva'] + ) + pressure_axis = source.Axis( + azimuth=mtlist['Paz'], plunge=mtlist['Ppl'], length=mtlist['Pva'] + ) + principal_axes = source.PrincipalAxes( + t_axis=tension_axis, p_axis=pressure_axis, n_axis=null_axis + ) + + # Create the Moment Tensor object with correct units and scaling + if units == "nm": + c = 1E-7 # conversion from dyne*cm to N*m + logger.debug(f"GeoNet moment tensor is in units of Newton*meters") + elif units == "dynecm": + c = 1 + logger.debug(f"GeoNet moment tensor is in units of dyne*cm") + + # CV is the conversion from non-units to the desired output units + cv = 1E20 * c + seismic_moment_in_nm = mtlist['Mo'] * c + + # Convert the XYZ coordinate system of GeoNet to an RTP coordinate system + # expected in the CMTSOLUTION file of Specfem + rtp = mt_transform(mt={"m_xx": mtlist['Mxx']*cv, "m_yy": mtlist['Myy']*cv, + "m_zz": mtlist['Mzz']*cv, "m_xy": mtlist['Mxy']*cv, + "m_xz": mtlist['Mxz']*cv, "m_yz": mtlist['Myz']*cv + }, + method="xyz2rtp" + ) + tensor = source.Tensor(m_rr=rtp['m_rr'], m_tt=rtp['m_tt'], + m_pp=rtp['m_pp'], m_rt=rtp['m_rt'], + m_rp=rtp['m_rp'], m_tp=rtp['m_tp'] + ) + # Create the source time function + source_time_function = source.SourceTimeFunction( + duration=2 * half_duration_from_m0(seismic_moment_in_nm) + ) + + # Generate a comment for provenance + comment = Comment(force_resource_id=True, + text="Automatically generated by Pyatoa via GeoNet MT CSV" + ) + + # Fill the moment tensor object + moment_tensor = source.MomentTensor( + force_resource_id=True, tensor=tensor, + source_time_function=source_time_function, + # !!! + # This doesn't play nice with obspy.Catalog.write(format='CMTSOLUTION') + # so ignore the origin id + # derived_origin_id=id_template.format('origin#ristau'), + scalar_moment=seismic_moment_in_nm, double_couple=mtlist['DC']/100, + variance_reduction=mtlist['VR'], comment=comment + ) + + # Finally, assemble the Focal Mechanism. Force a resource id so that + # the event can identify its preferred focal mechanism + focal_mechanism = source.FocalMechanism( + force_resource_id=True, nodal_planes=nodal_planes, + moment_tensor=moment_tensor, principal_axes=principal_axes, + comments=[comment] + ) + + return focal_mechanism + + +def query_geonet_mt_catalog(event_id, csv_fid=None): + """ + Get moment tensor information from a internal csv file, + or from an external github repository query. + Only relevant to the new zealand tomography problem. + Geonet moment tensors stored with a specific column format. + + :type event_id: str + :param event_id: unique event identifier + :type csv_fid: str + :param csv_fid: optional path to GeoNet CMT solution file that is stored + locally on disk, will be accessed before querying web service + :rtype moment_tensor: dict + :return moment_tensor: dictionary created from rows of csv file + """ + reader = None + if csv_fid is not None: + try: + reader = csv.reader(open(csv_fid, 'r'), delimiter=',') + except FileNotFoundError: + pass + + if reader is None: + # Request and open the CSV file. Assumed that GeoNet will keep their + # moment-tensor information in their GitHub repository + # Last accessed 23.6.19 + geonet_mt_csv = ( + "https://raw.githubusercontent.com/GeoNet/data/master/" + "moment-tensor/GeoNet_CMT_solutions.csv" + ) + response = requests.get(geonet_mt_csv) + if not response.ok: + raise FileNotFoundError(f"Response from {geonet_mt_csv} not ok") + + reader = csv.reader(response.text.splitlines(), delimiter=',') + + # Parse the CSV file + for i, row in enumerate(reader): + # First row contains header information + if i == 0: + tags = row + # First column gives event ids + if row[0] == event_id: + values = [] + # Grab the relevant information from the file + for t, v in zip(tags, row): + if t == "Date": + values.append(UTCDateTime(v)) + elif t == "PublicID": + values.append(v) + else: + values.append(float(v)) + + moment_tensor = dict(zip(tags, values)) + logger.info(f"geonet moment tensor found for: {event_id}") + return moment_tensor + else: + raise AttributeError(f"no geonet moment tensor found for: {event_id}") + + +def append_focal_mechanism_to_event(event, method="all", overwrite_focmec=False, + overwrite_event=False, client=None): + """ + Attempt to find focal mechanism information with a given ObsPy Event object. + + .. note:: + FDSN fetched events are devoid of a few bits of information that are + useful for our applications, e.g. moment tensor, focal mechanisms. + This function will perform the conversions and append the necessary + information to the event located in the dataset. + + :type event: obspy.core.event.Event + :param event: Event object to append a focal mechanism to. + :type method: bool + :param method: try to find correspondig focal mechanism + using various public catalogs. Currently available: + 'all': Try all available options in order until MT is found + 'USGS': Search the USGS moment tensor catalog + 'GCMT': Search the GCMT moment tensor catalog + False: Don't attempt to search for moment tensors + :type client: str + :param client: Specific `client`s come built-in with specific MT catalogs + If matching client, will ignore other MT choices: + 'GEONET': will search John Ristau catalog for moment tensors, + :type overwrite_focmec: bool + :param overwrite_focmec: If the event already has a focal mechanism, + overwrite the existing focal mechanism. + :type overwrite_event: bool + :param overwrite_event: A new event object is usually retrieved when + gathering MT from USGS or GCMT. Often the locations/timing of this event + are less accurate than the input event (which is usually sourced from + a regional catalog). This parameter controls which event object is + taken. If `True`, takes the USGS or GCMT catalog information, if `False` + only takes the focal mechanism attribute. + :rtype event: obspy.core.event.Event + :return event: Event with a new focal mechanism if one was found + :raises TypeError: if event is not provided as an obspy.core.event.Event + """ + if not isinstance(event, Event): + raise TypeError(f"`event` must be an ObsPy Event object, " + f"not: {type(event)}") + # If the event already has a focal mechanism attribute, don't gather + elif hasattr(event, "focal_mechanisms") and \ + event.focal_mechanisms and not overwrite_focmec: + logger.debug("event already has focal mechanism, will not attempt to" + "append new focal mechanism") + return event + # Only gather moment tensors if we're already trying to do FDSN stuff + elif client is None: + logger.debug("client not specified, will not attempt gathering " + "moment tensor") + return event + + method = method.upper() + event_id = event.resource_id.id # assuming datacenter tags ID with event id + cat = Catalog() + focal_mechanism = None + if client.upper() == "GEONET": + logger.info("querying GeoNet moment tensor catalog") + focal_mechanism = get_geonet_mt(event_id=event_id, units="nm") + else: + # Try 1: Look at USGS catalog + if method in ["ALL", "USGS"]: + logger.debug("querying USGS database for moment tensor") + cat = get_usgs_moment_tensor(event=event) + # Try 2: Look at GCMT catalog if USGS catalog did not return + elif (method in ["ALL", "GCMT"]) and len(cat) == 0: + logger.debug("querying GCMT database for moment tensor") + cat = get_gcmt_moment_tensor(event=event) + # Try ?: Add options below for more catalog selection + # +++++++++++++++++++++++++++++++++++++++++++++++++++ + # If multiple events found for a given set of event criteria, pick first + if cat is not None: + if len(cat) > 1: + logger.warning(f"multiple ({len(cat)}) events found, " + f"picking zeroth index") + # Distinguish `event_new` from `event`, sometimes you still want the + # catalog location, not the one from USGS or GCMT. Or if nothing was + # found, then we will return the same event + event_new = cat[0] + focal_mechanism = event_new.preferred_focal_mechanism() + # Append or overwrite focal mechanism or event + if focal_mechanism is None: + event_out = event + else: + if overwrite_event: + logger.debug("overwriting input event object with newly gathered " + "event containing focal mechanism") + event_out = event_new + else: + logger.debug("appending gathered focal mechanism to current event") + event_out = event.copy() + event_out.focal_mechanisms = [focal_mechanism] + event_out.preferred_focal_mechanism_id = focal_mechanism.resource_id + + return event_out + + +class Source: + """ + A generic Source object to characterize FORCESOLUTION files and SPECFEM2D + SOURCE files without breaking the architechture of Pyatoa which was built + around CMTSOLUTIONs and ObsPy Event objects + + Essentially this class tries to mimic the ObsPy Event object and return + required information that is queried throughout a Pyatoa workflow + """ + def __init__(self, resource_id, origin_time, longitude, latitude, depth): + """ + Only define the essential values required of a source + + :type resource_id: str + :param resource_id: unique label for the event + :type origin_time: str or UTCDateTime + :param origin_time: origin time for the event + :type longitude: float + :param longitude: longitude or X value of the event in the domain + :type latitude: float + :param latitude: latitude or Y value of the event in the domain + :type depth: float + :param depth: depth in km, inverted Z axis, positive values means deeper + """ + self.id = resource_id + self.time = UTCDateTime(origin_time) + self.longitude = float(longitude) + self.latitude = float(latitude) + self.depth = float(depth) + + def preferred_origin(self): + """ + Convenience function to mimic behavior of ObsPy Event object + """ + return self + + @property + def resource_id(self): + """ + Convenenience function to mimic behavior of ObsPy Event object + """ + return self +