Skip to content

Noise Toolkit PDF PSD package V.2

Manoch Bahavar edited this page Jan 28, 2021 · 5 revisions

DESCRIPTION (V.2.0):

The PDF/PSD package provides three highly configurable Python 3 scripts to calculate waveform spectra. This package takes advantage of FDSN Web service client for ObsPy to retrieve necessary waveform data and it also allows users to process waveform data from their local files. This package provides PSD file collections similar to the PQLX package (McNamara and Boaz, 2005, https://www.usgs.gov/software/pqlx-a-software-tool-evaluate-seismic-station-performance).

The scripts included in this package are:

  • ntk_computePSD.py – requests waveforms and response data for given station(s) using the ObsPy FDSN client OR
    reads user’s waveform data files in SAC, MSEED, CSS, etc. formats and computes PSDs and populates a file-based
    PSD database

    NOTE: ntk_computePSD.py first identifies the appropriate FDSN data provider for the requested station using the Fedcatalog service from IRIS (service.iris.edu/irisws/fedcatalog/1/ ) and then requests waveform/response data for that station using ObsPy’s FDSN client

  • ntk_extractPsdHour.py – extracts PSDs for a given channel and bounding parameters from the PSD database
    (the output is similar to PQLX’s exPSDhour script, https://pubs.usgs.gov/of/2010/1292/ )
  • ntk_binPsdDay.py – bins PSD’s to daily files for a given channel and bounding parameter

QUICK START

Installation:

See also the package’s INSTALL.txt file
  • install Python 3.8 or higher on your computer
  • additional required Python modules are:
    . obspy 1.2.2
    . matplotlib 3.3.3
    . numpy 6.2
    . scipy 1.5.3
  • download the package under the installation directory.
Note: if you are installing this package on a PC, edit the param/shared.py file and set the namingConvention parameter to ‘WINDOWS

Running the scripts:

  • under the root directory of the Noise Toolkit execute the following command to plot PSD of TA.O18A.—.BHZ for 12 pm to 1 pm on 2008-08-14

python bin/ntk_computePSD.py net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T13:00:00 xtype=period plot=1

  • review all parameter files to get familiar with the script capabilities
  • use examples given below to run/test the scripts and become familiar with them
  • copy and edit the parameter files under the param directory to customize parameters

NOTES

  • for more information visit the IRIS DMC Noise Toolkit Data Product web page at:
    ds.iris.edu/ds/products/noise-toolkit
  • examples below may turn on the verbose mode. Once you have configured the script, you can turn it off

COMMENTS/QUESTIONS/SUGGESTIONS:

  • We welcome patches and enhancements to this software. When developing patches, please pay particular attention to ease of use and maintenance and also keep dependencies to a minimum.
  • for issues, file a ticket under Issues

HISTORY

  • 2020-11-16 V.2.0.0: Python 3, use of Fedcatalog, adoption of CSD changes and adoption of PEP 8 style guide.
  • 2019-09-09 Robert Anthony (USGS, Albuquerque Seismological Laboratory):
    Using CSD to compute the cross spectral density of two signals
  • 2017-01-18 V.0.9.5: support for reading data and metadata from files only with no Internet requirement
  • 2016-11-01 V.0.9.0: added support for obtaining channel responses from local station XML response files by introducing the following two functions in tsLIB.py:
    getResponseInventory – to build a list of response inventories under a given met directory
    getResponseFromFile – to find response for a given Network, Station, Location and channel
    added respDirectory to common.py parameter file to disable looking for response files on local drives, set this parameter to None. Otherwise, set it to the response directory path
  • 2016-01-25 V.0.8.2: added polarization parameters to common.par and also made changes to some libraries in support of the polarization package
    added user and password parameters to common.par and ntk_computePSD.py in support of restricted data access
  • 2015-04-30 V.0.8.1: added powerDirectory and imageDirectory parameters to common.par in support of ME package.
  • 2015-04-07 V.0.8.0: now produces clear messages for missing parameters in the parameter file. Addresses reported maximum period values. Minor bugs and enhancements
  • 2015-02-24 V.0.7.0: introduced two new parameters (performInstrumentCorrection, applyScale) to allow user avoid instrument correction also now user can turn od decon. filter
  • 2014-11-24 V.0.6.2: documentation update and expansion of the common.py file to support Microseism Energy package
  • 2014-10-16 V.0.6.0: ntk_extractPsdHour.py output file name now includes the x-axis type
  • 2014-10-01 V.0.5.0: initial Beta release. Made the smoothing configurable, reorganized parameters for easier maintenance
  • 2014-05-20 modified the output format to:
    1. compute hourly PSDs and store them under the psdDb directory (this will be similar to PQLX’s database) as basis for script outputs
    2. provide the extractPsdHourly.py script to extract PSD’s for a given channel and bounding parameters. The output is similar to PQLX’s exPSDhour script
    3. provide the extractPdf.py that will extract PSDs to create (hourly and daily) PDFs for a given channel and bounding parameters.
    The output is similar to the current PDF output from the PDF-PSD data product (ds.iris.edu/ds/products/pdf-psd/)
Users should be able to use the above outputs with their existing scripts that work with the PQLX output
  • 2014-03-19: Added data file input option
  • 2014-03-07: Initial Alpha release

USER’S GUIDE

Use below examples to run/test the scripts and become familiar with them and edit the parameter files under the param directory to change parameters based on your needs.

Parameters

All parameters are under the param directory. The “shared.py” parameter file contains basic parameters shared by all parameter files. Please review ALL parameter to familiarize yourself with the script capabilities. For more information about this package, visit the IRIS DMC Noise Toolkit Data Product web page at:

http://ds.iris.edu/ds/products/noise-toolkit/

SCRIPT: ntk_computePSD.py

A Python 3 script that calculates average power spectral densities for a given station. The script:
* identifies the FDSN data provider for the requested station using the Fedcatalog service
from IRIS (https://service.iris.edu/irisws/fedcatalog/1/)
* requests waveform and response data for the given station(s)/channels(s) using ObsPy’s FDSN client
OR
* reads user-supplied waveform data files in SAC, MSEED, CSS, etc. format from a local disk

Then
* computes PSDs for the waveform data and populates a file-based PSD database

USAGE:

	python bin/ntk_computePSD.py to display the usage message
	  OR
	python bin/ntk_computePSD.py param=FileName client=[FDSN|FILES] net=network sta=station loc=location chan=channel(s) start=YYYY-MM-DDTHH:MM:SS end=YYYY-MM-DDTHH:MM:SS xtype=[period|frequency] plot=[0|1] verbose=[0|1] timing=[0|1]


	to perform computations where:
	  param		[default: computePSD] the configuration file name 
	  client	[default: FDSN] client to use to make data/metadata requests (FDSN or FILES) 
	  net		[required] network code
	  sta		[required] station code
	  loc		[required] location ID
	  chan		[default: BH?] channel ID(s); separate multiple channel codes by comma (no space)
	  xtype		[period or frequency, default: period] X-axis type (period or frequency for outputs and plots)
	  start		[required] start date-time (UTC) of the interval for which PSDs to be computed (format YYYY-MM-DDTHH:MM:SS)
	  end		[required] end date-time (UTC) of the interval for which PSDs to be computed (format YYYY-MM-DDTHH:MM:SS)
	  sw_width	[default: 0.25] Smoothing window width in octave
	  sw_shift	[default: 0.125] Smoothing window shift in fraction of octave
	  plotnm	[0 or 1, default 1] plot the New High/Low Noise Models [0|1]
	  plot		[0 or 1, default: 0] to run in plot mode set to 1
	  timing	[0 or 1, default: 0] to run in timing mode (set to 1 to output run times for different segments of the script)
	  verbose	[0 or 1, default: 0] to run in verbose mode set to 1

NOTES & EXAMPLES:

Use the run in plot (plot=1) or verbose mode (verbose=1_ to tune the parameters before a production run (verbose=0 plot=0):

python bin/ntk_computePSD.py net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T13:00:00 xtype=period plot=1

Use the x-axis type to not only define the axis type but also ask the script to create regular sampling in the corresponding log period or log frequency The first sample is placed at xStart second or xStart Hz as defined in the parameter file and log-regular samples are created on both sides. User may set xStart=‘Nyquist’ in the parameter file to set the Nyquist frequency as the reference. The former option will result the same period or frequency samples regardless of the station sampling interval

The following example limits channel to BHZ, processes only one hour (default) and plot the PSD with more smoothing (sw_width=0.5 and sw_shift=0.25):

python bin/ntk_computePSD.py param=computePSD net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T13:00:00 xtype=period sw_width=0.5 sw_shift=0.25 plot=1


With even more smoothing (sw_width=1):

python bin/ntk_computePSD.py param=computePSD net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T13:30:00 xtype=period sw_width=1 plot=1


Your request could extend over multiple days:

python bin/ntk_computePSD.py param=computePSD net=NM sta=SLM loc=DASH chan=BHZ start=2009-03-01T00:00:00 end=2009-03-31T00:00:00 xtype=period plot=0 verbose=0

If you already have your waveform data saved on a local storage under file formats such as SAC, MSEED, CSS, etc. (see the supported-formats at docs.obspy.org/packages/autogen/obspy.core.stream.read.html#obspy.core.stream.read) then you can instruct the script to get the waveform data from these files and connect to Web Services for the response information only. To do this you need to update two parameters in your parameter file (such pas param/computePSD.py):
set requestClient = "FILES" to flag the script that the waveform data are coming from file and fileTag = "{IRIS_NTK_PSD}/data/TEST/SAC/W*.SAC" to tell it which files to look at

python bin/ntk_computePSD.py net=TA sta=W5 loc=DASH start=2014-03-17T04:30:00 end=2014-03-17T05:30:00 type=period client=FILES plot=1

Requesting data from a FDSN data center, other than IRIS, is automatic and depends on the station you are requesting. For example requesting BHZ channel for GR.BFO:
python bin/ntk_computePSD.py param=computePSD net=GR sta=BFO loc=DASH chan=BHZ start=2020-10-01T00:00:00 end=2020-10-01T01:00:00 xtype=period plot=1

SCRIPT: ntk_extractPsdHour.py

a Python script to extract PSDs for a given channel and bounding parameters.

USAGE:

	ntk_extractPsdHour.py to display the usage message 

	  OR

	python bin/ntk_extractPsdHour.py param=FileName net=network sta=station loc=location chan=channel(s) start=YYYY-MM-DDTHH:MM:SS end=YYYY-MM-DDTHH:MM:SS xtype=[period|frequency] verbose=[0|1]

	to perform binning:
	  param		[default: extractPsdHour] the configuration file name 
	  net		[required] network code
	  sta		[required] station code
	  loc		[required] location ID
	  chan		[required] channel ID 
	  xtype		[required] X-axis type (period or frequency for output)
	  start		[required] start date-time (UTC) of the interval for which PSDs to be computed (format YYYY-MM-DDTHH:MM:SS)
	  end		[required] end date-time (UTC) of the interval for which PSDs to be computed (format YYYY-MM-DDTHH:MM:SS)
	  verbose	[0 or 1, default: 0] to run in verbose mode set to 1

INPUT:

The ntk_computePSD.py output files

OUTPUT FORMAT:

Date Hour X-value (period/frequency) Power with values separated using the “separator” character specified in the parameter file.

OUTPUT FILE:

	Full path to the  output data file is provided at the end of the run. 

	The output file name has the form:
		net.sta.loc.chan.start.end.xtype.txt
	for example:
		TA.O18A.--.BHZ.2008-08-14.2008-08-14.period.txt

EXAMPLES:

  • usage:
    python bin/ntk_extractPsdHour.py
  • Assuming that you already have tried the following ntk_compute_PSD.py example:
    python bin/ntk_computePSD.py net=TA sta=O18A loc=DASH start=2008-08-14T12:00:00 end=2008-08-14T13:30:00
  • Then you can perform PSD extraction via:
    python bin/ntk_extractPsdHour.py net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T12:30:00 xtype=period

SCRIPT: ntk_binPsdDay.py

A Python 3 script to bin PSDs to daily files for a given channel and bounding parameters.

USAGE:

	ntk_binPsdDay.py to display the usage message 

	  OR

	ntk_binPsdDay.py param=FileName net=network sta=station loc=location chan=channel(s) start=YYYY-MM-DDTHH:MM:SS end=YYYY-MM-DDTHH:MM:SS xtype=[period|frequency] verbose=[0|1]

	where:
	  param		[default: binPsdDay] the configuration file name 
	  net		[required] network code
	  sta		[required] station code
	  loc		[required] location ID
	  chan		[required] channel ID 
	  xtype		[required] X-axis type for output (period or frequency)
	  start		[required] start date-time (UTC) of the binning interval (format YYYY-MM-DDTHH:MM:SS)
	  end		[required] end date-time (UTC) of the binning interval (format YYYY-MM-DDTHH:MM:SS)
	  verbose	[0 or 1, default: 0] to run in verbose mode set to 1

INPUT:

The ntk_computePSD.py output files

OUTPUT FORMAT:

hour X-value (period/frequency) Power with values separated using the “separator” character specified in the parameter file.
h3. OUTPUT FILE:

Full paths to the daily and hourly output data files are provided at the end of the run. You may turn off hourly file output via the configuration file.

The daily and hourly output file names have the following form respectively:
D?.bin and H?.bin (??? represents 3-digit day of the year)
for example:
D227.bin and H227.bin

NOTES & EXAMPLES:

* usage:
python bin/ntk_binPsdDay.py

* Assuming that you already have tried the following ntk_compute_PSD.py example:
python bin/ntk_computePSD.py net=TA sta=O18A loc=DASH start=2008-08-14T12:00:00 end=2008-08-14T13:30:00

* Then you can perform binning via:
python bin/ntk_binPsdDay.py net=TA sta=O18A loc=DASH chan=BHZ start=2008-08-14T12:00:00 end=2008-08-14T13:30:00 xtype=period

COPYRIGHT © 2020 Product Team, IRIS Data Management Center

This package is provided by the IRIS DMC Data Products Team WITHOUT ANY WARRANTY AND/OR SUPPORT This is a free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. This script is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License (GNU-LGPL) for more details. The GNU-LGPL and further information can be found here: www.gnu.org You should have received a copy of the GNU Affero General Public License along with this program. If not, see http://www.gnu.org/licenses.