PySEP v0.3.0 (#78)

* Fix various bugs and adds features to code (#60) * renamed specfem station writing function, and added an ordering component * bumping version in setup and changing url and author * small mt typo change * added an example config that gathers a small amount of data for testing/dev purposes * added a new kwarg 'order_stations_list_by' which sets the order of the output stations list, related to #36 * bugfix typo parameter call * bugfix rotating was not actually rotating streams in place for arbitrary components ->ENZ * bugfix fixed rotation testing, which was not actually evaluating as expected. some confusion about editing streams and inventories in place which was not actually happening. Fixed and now rotate test works as expected * added 12Z test data for rotating 12Z -> ENZ which works now, added test to cover this case * moved 2012 central alaska IU gather script (which gathers 12Z component) into examples * bugfix pysep was not setting log level properly from the config file. now in the load() function, config log level is allowed to override logger level * changed config examples dir name to test, to be used as a test bed for checking specific data gathering features such as data gap removal * bugfix #57, reading event origin time from sac header was shifting times unexpectedly * added utility function to convert obspy catalog to event input list * Bugfix missing rotation metadata throwing TypeError (#63) * bugfix clipped amplitude check was not properly checking the data array * added error catch and release during rotate code to avoid processing breakdown. also rotate now goes by channel rather than the entire stream at once * rearranged rotation algorithm to be a bit cleaner, previously too many different streams floating around and being fed to one another. * pysep rotation now has additional checks for station metadata and error catching to get around any rotation bugs that used to cause the rotation to fail added new tests to cover rotation and fixed a rotation test that was failing due to change added a config for testing the rotation bug, providing station info * updating bug config file * feature mass download (#64) * fixed exclude string making and added a bool check for deleting tmpdir in mass download function * removed accidental debug statement * Remove llnl_db_client hard dependency (#65) * moved llnl_db_client as a hard dep of pysep and moved it into an optional dependency hid llnl_db_client import statement behind a try-except block and threw a check in the init of pysep only requiring this dep if the client choice is LLNL * Added setup.py legacy file * removed llnl dep from conda env yml file * tested install with and without llnl dep and tested check statement for missing llnl import * Feature: event declustering and source receiver weighting (#66) * started adding declustering rewrite with test and test data * finished cartesian gridding decluster script, added radial gridding decluster script, added tests to cover both and plotting * decluster catalog added feature to threshold events by min magnitude and data availability separate from declustering added plot feture to connect source receiver pairs based on data availability added feature to allow sorting by data availability * adjusting default figure names for declust plotting scripts * cleaned up plotting functions into functions to avoid repeat plotting calls, added removed events to plotting routines * small typo updte declust * moved some declustering functions into util * declust started srcrcv weighting function * individual source or receiver list weighting working with the smart scan feature next up is to implement the entire weighting scheme with normalization and based on data availability * declust slow progression * finished srcrcv weighting scheme * finished declustering and weighting algos with additional plotting and text file writing for weights * added basic test for srcrcv weight calc * clean up docstrings of pysep/recsec, bugfix pysep restrictions (#68) * reformatted PySEP init docstring, added missing parameters and categorized parameters to make things easier to find * removed hard requirement that user provides event depth and magnitude if event selection is default * added boolean flag to toggle insufficient length checker bugfix: added remove_clipped boolean toggle on the actual processing step, which was not there before * bugfix phase list passed from pysep into util function for phase getting * last minute touch ups on docstring * recsec removed unneeded myround function from top of script * categorized recsec init docstring for easier readability * bugfix: rotate parameter check function was setting rotate parameter as an empty list, but it was expected as NoneType within the main processing function. removed this type conversion * bugfix: log message failing when magnitude or depth set to NoneType * bugfix allow NoneType event magnitude and depth, ass ign dummy values to cap header because these are required by record section * Update docs (#69) * moved docstring of pysep into class description for autoapi * instantiated sphinx for doc creation and autoapi building, reorganized pysep and recsec docstrings to be well formatted for autoapi * migrated wiki docs into docs directory * fixed docstrings and corrected links for API references in docs * fixed up typos and cleaned up docs * Delete jekyll-gh-pages.yml * Bugfix: unable to set event_depth_km or event_magnitude as NoneType (#72) * allow taup arrival time get to be skipped if event depth not provided * sac header functions that required 'evdp' are now skipped over if 'evdp' is not present in the sac header * Update pysep.py bugfix: removing debug statement left in pysep * Feature recsec tmarks (#73) * replace plotw_rs with record section run command in pysep plotting, sets default sorting to 'distance' and not 'distance_r' for pysep-generated record section * working on ordering of multiple pages * added feature tmark to add static lines at given time values, does not address time shifts or move out * Bugfix: plotw rs sort (#76) * reformatting plotw_rs page separation logic * sort order of multi-page record sections is now reversed to be more natural * Update CHANGELOG.md * version bump v0.3.0
adjtomo · Feb 18, 2023 · 939758f · 939758f
1 parent e7e1643
commit 939758f
Show file tree

Hide file tree

Showing 40 changed files with 39,666 additions and 349 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -0,0 +1,19 @@
+# PySEP Changelog
+
+## Version 0.2.0 (Master)
+
+- RecSec plot aesthetics altered to provide more information efficiently  
+- Shifted documentation to GitHub Wiki page
+- Added moment tensor gathering routines from Pyatoa
+
+## Version 0.3.0 (Devel)
+
+- Removed llnl_db_client dependency from package and made it optional 
+- Re-instated 'mass_downloader' option 
+- Added `Declust` class for declustering and weighting
+- Completed and re-organized PySEP and RecSec docstrings 
+- Bugfix: rotation was not able to be set as null 
+- Bugfix: some flag check parameters were not being used during processing
+- Removed hard restriction on requiring event depth and magnitude for default event selection
+- Add feature 'tmarks' to record section to plot vertical lines at reference times
+- Fix ordering of multi-page record sections when using plotw_rs 
diff --git a/docs/Makefile b/docs/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/docs/conf.py b/docs/conf.py
@@ -0,0 +1,60 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# For the full list of built-in configuration values, see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+def skip(app, what, name, obj, would_skip, options):
+    if name == "__init__":
+        return False
+    return would_skip
+
+def setup(app):
+    app.connect("autodoc-skip-member", skip)
+
+# -- Project information -----------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
+
+project = 'PySEP'
+copyright = '2023, adjTomo Dev Team'
+author = 'adjTomo Dev Team'
+release = ''
+version = '0.2.0'
+
+# -- General configuration ---------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
+
+extensions = [
+        'sphinx.ext.doctest',
+        'sphinx.ext.intersphinx',
+        'sphinx.ext.todo',
+        'sphinx.ext.coverage',
+        'sphinx.ext.viewcode',
+        'sphinx.ext.napoleon',
+        'sphinx.ext.autosummary',
+        "autoapi.extension",
+        "myst_parser"
+        ]
+
+templates_path = ['_templates']
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# Need to tell the autoapi that our source code is one level up
+autoapi_type = "python"
+autoapi_dirs = [
+        "../pysep",
+        "../pysep/utils",
+        "../pysep/tests"
+        ]
+autoapi_add_toctree_entry = True
+autoapi_python_class_content = 'both'
+# autoclass_content = 'both'  # show init
+
+source_suffix = {'.rst': 'restructuredtext',
+                 '.md': 'markdown'}
+
+
+# -- Options for HTML output -------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
+
+html_theme = 'furo'
+html_static_path = ['_static']
diff --git a/docs/cookbook.md b/docs/cookbook.md
@@ -0,0 +1,250 @@
+# Cookbook
+
+Here we provide some code snippets using PySEP routines or utilities to perform
+tasks useful for manipulating or preparing data for moment tensor and waveform 
+inversion software.
+
+## Create STATIONS file for SPECFEM
+
+SPECFEM requires a STATIONS file that defines geographical coordinates of stations that
+will be used in a simulation. Using PySEP and ObsPy, you can quickly generate a STATIONS
+file for a given region
+
+```python
+from obspy import UTCDateTime
+from obspy.clients.fdsn import Client
+from pysep.utils.io import write_stations_file
+
+c = Client("IRIS")
+# This is for the region of northern Alaska, broadband seismometers only
+# collecting online stations from Jan. 1, 2000 until present time
+inv = c.get_stations(network="*", station="*", 
+                     channel="BH?,BL?,HH?,HL?",
+                     minlatitude=64.5, maxlatitude=72., 
+                     minlongitude=-168., maxlongitude=-140.,
+                     level="station", 
+                     starttime=UTCDateTime("2000-01-01T00:00:00"), 
+                     endtime=UTCDateTime()
+                     )
+
+# Write 'STATIONS' file, ordered alphabetically by station name
+write_specfem_stations_file(inv, fid="STATIONS", order_by="station")
+```
+
+```
+$ head STATIONS
+  A19K    AK     70.2043   -161.0713    0.0    0.0
+  A21K    AK     71.3221   -156.6175    0.0    0.0
+  A22K    AK     71.0033   -154.9742    0.0    0.0
+   ANM    AK     64.5646   -165.3732    0.0    0.0
+  B18K    AK     69.3641   -161.8016    0.0    0.0
+  B20K    AK     70.0079   -157.1599    0.0    0.0
+  B22K    AK     70.3400   -153.4196    0.0    0.0
+  C16K    AK     68.2746   -165.3436    0.0    0.0
+  C18K    AK     68.6483   -161.1943    0.0    0.0
+  C19K    AK     69.1049   -159.5874    0.0    0.0
+
+```
+
+## Fetch moment tensors as CMTSOLUTION files
+
+It's often useful to gather a catalog of events and their moment tensors as a starting point 
+for waveform inversion. PySEP has some built in moment tensor routines that can help us 
+with that.
+
+### From the GCMT catalog
+
+This routine queries the GCMT catalog for moment tensors matching a given origin time and
+magnitude. Given an example event, we can return a Catalog object containing the GCMT 
+moment tensor. Taking this [M6.4 in northern Alaska](https://earthquake.alaska.edu/event/018aap2cqu) as an example.
+
+```python
+from obspy import UTCDateTime
+from pysep.utils.mt import get_gcmt_moment_tensor
+
+origintime = UTCDateTime("2018-08-12T14:58:53")
+magnitude=6.4
+cat = get_gcmt_moment_tensor(origintime=origintime, magnitude=magnitude)
+cat.write("CMTSOLUTION_GCMT", format="CMTSOLUTION")
+```
+
+```
+$ cat CMTSOLUTION_GCMT
+ PDE 2018 08 12 14 58 54.30   69.5600 -145.3000   2.2 6.4 6.4 NORTHERN ALASKA
+event name:   C201808121458A
+time shift:           7.9000
+half duration:        4.0000
+latitude:            69.7400
+longitude:         -144.7800
+depth:               12.0000
+Mrr:           -7.690000E+24
+Mtt:           -1.990000E+25
+Mpp:            2.760000E+25
+Mrt:            4.890000E+24
+Mrp:           -1.510000E+25
+Mtp:           -4.730000E+25
+```
+
+### From USGS Catalog
+
+Using the event we grabbed in the previous section (from GCMT), we can query 
+the USGS catalog for their moment tensor solution.
+
+>__NOTE__: USGS moment tensor searches are more strict and require an existing 
+  Event object from which it takes the hypocentral location, origin time and magnitude
+
+```python
+# Continuing from codeblock above
+from pysep.utils.mt import get_usgs_moment_tensor
+
+cat = get_usgs_moment_tensor(event=cat[0])
+cat.write("CMTSOLUTION_USGS", format="CMTSOLUTION")
+UserWarning: Moment tensor has no source time function. Half duration will be set to 1.0.
+```
+
+```
+$ cat CMTSOLUTION_USGS
+ PDE 2018 08 12 14 58 53.50   69.5762 -145.2910  15.8 6.4 6.4 NORTHERN ALASKA    
+event name:89 km SW of Kaktovik, Alaska                                          
+time shift:          -0.2030                                                     
+half duration:        1.0000                                                     
+latitude:            69.6864                                                     
+longitude:         -144.8187                                                     
+depth:               11.5000                                                     
+Mrr:           -1.015100E+25                                                     
+Mtt:           -1.644300E+25                                                     
+Mpp:            2.659300E+25                                                     
+Mrt:           -8.487000E+24                                                     
+Mrp:           -1.468700E+25                                                     
+Mtp:           -4.792400E+25  
+```
+
+### From GeoNet Catalog (New Zealand)
+
+For those working with New Zealand data, it is possible to grab 
+moment tensor data from [John Ristau's catalog](https://github.com/GeoNet/data/tree/main/moment-tensor). 
+For these events you will need to provide a corresponding GeoNet event id. 
+We'll take [event 2018p130600](https://www.geonet.org.nz/earthquake/2018p130600) as an example.
+
+```python
+from obspy.clients.fdsn import Client
+from pysep.utils.mt import get_geonet_mt
+
+c = Client("GEONET")
+cat = c.get_events(eventid="2018p130600")
+focal_mechanism = get_geonet_mt(event_id="2018p130600", units="nm")  # units can also be dynecm
+cat[0].focal_mechanisms = [focal_mechanism]
+cat.write("CMTSOLUTION_GEONET", format="CMTSOLUTION")
+```
+```
+$ cat CMTSOLUTION_GEONET
+ PDE 2018 02 18 07 43 48.13  -39.9490  176.2995  20.6 5.2 5.2 NORTH ISLAND, NEW ZEALAND
+event name:           8E23C1
+time shift:           0.0000
+half duration:        0.6989
+latitude:           -39.9490
+longitude:          176.2995
+depth:               20.5946
+Mrr:           -2.479380E+23
+Mtt:            1.314880E+23
+Mpp:            1.164500E+23
+Mrt:            5.032500E+22
+Mrp:            6.607700E+22
+Mtp:            9.359300E+22
+```
+-------------------------
+
+## Convert SPECFEM-generated synthetics to SAC files
+
+Some may find it useful to convert synthetic seismograms into SAC files to simplify later processing. 
+The following code snippet will write out SAC files on a per-component basis for synthetics generated 
+in a geographic coordinate system (lat/lon) from SPECFEM2D, 3D or 3D_GLOBE.
+
+>__NOTE:__ The `read_sem` function automatically appends SAC headers using the provided metadata
+
+```python
+# cd path/to/specfem_workdir
+from glob import glob
+from pysep.utils.io import read_sem
+
+for fid in glob("./OUTPUT_FILES/*sem*):
+    st = read_sem(fid=fid, source="./DATA/CMTSOLUTION", stations="./DATA/STATIONS")
+    st.write(os.path.basename(fid), format="SAC")
+```
+
+If your synthetics were generated in a Cartesian coordinate system (XYZ) you will need to use a 
+separate function, as ObsPy does not like coordinate systems that do not adhere to WGS84
+
+```python
+# cd path/to/specfem_workdir
+from glob import glob
+from pysep.utils.io import read_sem_cartesian
+
+for fid in glob("./OUTPUT_FILES/*sem*):
+    st = read_sem_cartesian(fid=fid, source="./DATA/CMTSOLUTION", stations="./DATA/STATIONS")
+    st.write(f"{st[0].get_id()}.sac", format="SAC")
+```
+
+-------------------------
+
+## 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 obspy import read, read_events, read_inventory
+from pysep.utils.cap_sac import append_sac_headers, format_sac_header_w_taup_traveltimes
+
+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")
+```
+
+-------------------------
+
+##  Set custom order on output station file 
+
+PySEP creates a text file for stations gathered. This list contains station IDs, coordinates,
+epicentral distance and azimuth from a given event. By default, this list is ordered by the 
+internal order of the station metadata (usually alphabetical). 
+
+Users can set the order of this station list by specifying the `order_stations_by` parameter.
+This is provided as a keyword argument to the initiation of PySEP and can be set with:
+
+```bash
+pysep -c pysep_config.yaml --order_stations_list_by distance
+```
+
+Acceptable arguments for this are: 'network', 'station', 'latitude', 'longitude', 'distance', 'azimuth'
+
+You can also call this from a Python environment
+
+```python
+from pysep import Pysep
+
+sep = Pysep(..., order_stations_list_by="distance")
+```
+
+
+-------------------------
+
+##  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")
+```
+