obspy.signal.spectral_estimation.PPSD

class PPSD(stats, metadata, skip_on_gaps=False, db_bins=(-200, -50, 1.0), ppsd_length=3600.0, overlap=0.5, special_handling=None, period_smoothing_width_octaves=1.0, period_step_octaves=0.125, period_limits=None, **kwargs)[source]

Bases: object

Class to compile probabilistic power spectral densities for one combination of network/station/location/channel/sampling_rate.

Calculations are based on the routine used by [McNamara2004]. For information on New High/Low Noise Model see [Peterson1993].

Basic Usage

>>> from obspy import read
>>> from obspy.signal import PPSD
>>> st = read()
>>> tr = st.select(channel="EHZ")[0]
>>> paz = {'gain': 60077000.0,
...        'poles': [-0.037004+0.037016j, -0.037004-0.037016j,
...                  -251.33+0j, -131.04-467.29j, -131.04+467.29j],
...        'sensitivity': 2516778400.0,
...        'zeros': [0j, 0j]}
>>> ppsd = PPSD(tr.stats, paz)
>>> print(ppsd.id)
BW.RJOB..EHZ
>>> print(ppsd.times_processed)
[]

Now we could add data to the probabilistic psd (all processing like demeaning, tapering and so on is done internally) and plot it like …

>>> ppsd.add(st) 
>>> print(ppsd.times_processed) 
[UTCDateTime(...), UTCDateTime(...), ..., UTCDateTime(...)]
>>> ppsd.plot() 

… but the example stream is too short and does not contain enough data.

Note

For a real world example see the ObsPy Tutorial.

Saving and Loading

The PPSD object supports saving to a numpy npz compressed binary file:

>>> ppsd.save_npz("myfile.npz") 

The saved PPSD can then be loaded again using the static method load_npz(), e.g. to plot the results again. If additional data is to be processed (note that another option is to combine multiple npz files using add_npz()), metadata must be provided again, since they are not stored in the numpy npz file:

>>> ppsd = PPSD.load_npz("myfile.npz")  

Note

When using metadata from an Inventory, a Parser instance or from a RESP file, information on metadata will be correctly picked for the respective starttime of the data trace. This means that instrument changes are correctly taken into account during response removal. This is obviously not the case for a static PAZ dictionary!

Attributes

NPZ_SIMPLE_TYPE_MAP

NPZ_SIMPLE_TYPE_MAP_R

NPZ_STORE_KEYS

NPZ_STORE_KEYS_ARRAY_TYPES

NPZ_STORE_KEYS_LIST_TYPES

NPZ_STORE_KEYS_SIMPLE_TYPES

NPZ_STORE_KEYS_VERSION_NUMBERS

channel

current_histogram

current_histogram_count

current_histogram_cumulative

current_times_used

db_bin_centers

db_bin_edges

delta

len

Trace length for one psd segment.

location

merge_method

network

nfft

nlap

period_bin_centers

Return centers of period bins (geometric mean of left and right edge of period smoothing ranges).

period_bin_left_edges

Returns left edges of period bins (same length as number of bins).

period_bin_right_edges

Returns right edges of period bins (same length as number of bins).

period_xedges

Returns edges of period histogram bins (one element longer than number of bins).

psd_frequencies

psd_periods

psd_values

Returns all individual smoothed psd arrays as a list.

station

step

Time step between start times of adjacent psd segments in seconds (assuming gap-less data).

times_data

times_gaps

times_processed

Public Methods

add

Process all traces with compatible information and add their spectral estimates to the histogram containing the probabilistic psd.

add_npz

Add previously computed PPSD results to current PPSD instance.

calculate_histogram

Calculate and set current 2D histogram stack, optionally with start- and endtime and time of day restrictions.

extract_psd_values

Extract PSD values for given period in seconds.

get_mean

Returns periods and mean psd values (i.e.

get_mode

Returns periods and mode psd values (i.e.

get_percentile

Returns periods and approximate psd values for given percentile value.

load_npz

Load previously computed PPSD results.

plot

Plot the 2D histogram of the current PPSD.

plot_coverage

Plot the data coverage of the histogram of the current PPSD.

plot_spectrogram

Plot the temporal evolution of the PSD in a spectrogram-like plot.

plot_temporal

Plot the evolution of PSD value of one (or more) period bins over time.

save_npz

Saves the PPSD as a compressed numpy binary (npz format).

Private Methods

Warning

Private methods are mainly for internal/developer use and their API might change without notice.

PPSD._add_npz(filename, allow_pickle=False)[source]

See PPSD.add_npz().

PPSD._get_gapless_psd()[source]

Helper routine to get a list of 2-tuples with gapless portions of processed PPSD time ranges. This means that PSD time history is split whenever to adjacent PSD timestamps are not separated by exactly self.ppsd_length * (1 - self.overlap).

PPSD._get_plot_title()[source]
PPSD._get_response(tr)[source]
PPSD._get_response_from_inventory(tr)[source]
PPSD._get_response_from_parser(tr)[source]
PPSD._get_response_from_paz_dict(tr)[source]
PPSD._get_response_from_resp(tr)[source]
PPSD._get_times_all_details()[source]
PPSD._plot_histogram(fig, draw=False, filename=None)[source]

Reuse a previously created figure returned by plot(show=False) and plot the current histogram stack (pre-computed using calculate_histogram()) into the figure. If a filename is provided, the figure will be saved to a local file. Note that many aspects of the plot are statically set during the first plot() call, so this routine can only be used to update with data from a new stack.

PPSD._setup_period_binning(period_smoothing_width_octaves, period_step_octaves, period_limits)[source]

Set up period binning.

PPSD._split_lists(times, psds)[source]
PPSD._stack_selection(starttime=None, endtime=None, time_of_weekday=None, year=None, month=None, isoweek=None, callback=None)[source]

For details on restrictions see calculate_histogram().

Return type:

numpy.ndarray of bool

Returns:

Boolean array of which psd pieces should be included in the stack.

Special Methods

PPSD.__check_histogram()
PPSD.__check_time_present(utcdatetime)

Checks if the given UTCDateTime is already part of the current PPSD instance. That is, checks if from utcdatetime to utcdatetime plus ppsd_length there is already data in the PPSD. Returns True if adding ppsd_length starting at the given time would result in an overlap of the ppsd data base, False if it is OK to insert this piece of data.

PPSD.__insert_data_times(stream)

Gets gap information of stream and adds the encountered gaps to the gap list of the PPSD instance.

PPSD.__insert_gap_times(stream)

Gets gap information of stream and adds the encountered gaps to the gap list of the PPSD instance.

PPSD.__insert_processed_data(utcdatetime, spectrum)

Inserts the given UTCDateTime and processed/octave-binned spectrum at the right position in the lists, keeping the order intact.

Replaces old obspy.signal.spectral_estimation.PPSD.__insert_used_time() private method and the addition ot the histogram stack that was performed directly in obspy.signal.spectral_estimation.PPSD.__process().

PPSD.__invalidate_histogram()
PPSD.__plot_coverage(ax)

Helper function to plot coverage into given axes.

PPSD.__process(tr)

Processes a segment of data and save the psd information. Whether Trace is compatible (station, channel, …) has to checked beforehand.

Parameters:

tr (Trace) – Compatible Trace with data of one PPSD segment

Returns:

True if segment was successfully processed, False otherwise.

PPSD.__sanity_check(trace)

Checks if trace is compatible for use in the current PPSD instance. Returns True if trace can be used or False if not.

PPSD.__delattr__(name, /)

Implement delattr(self, name).

PPSD.__dir__()

Default dir() implementation.

PPSD.__eq__(value, /)

Return self==value.

PPSD.__format__(format_spec, /)

Default object formatter.

PPSD.__ge__(value, /)

Return self>=value.

PPSD.__getattribute__(name, /)

Return getattr(self, name).

PPSD.__gt__(value, /)

Return self>value.

PPSD.__hash__()

Return hash(self).

PPSD.__init__(stats, metadata, skip_on_gaps=False, db_bins=(-200, -50, 1.0), ppsd_length=3600.0, overlap=0.5, special_handling=None, period_smoothing_width_octaves=1.0, period_step_octaves=0.125, period_limits=None, **kwargs)[source]

Initialize the PPSD object setting all fixed information on the station that should not change afterwards to guarantee consistent spectral estimates. The necessary instrument response information can be provided in several ways using the metadata keyword argument:

  • Providing an Inventory object (e.g. read from a StationXML file using read_inventory() or fetched from a FDSN webservice).

  • Providing an obspy.io.xseed.parser.Parser, (e.g. containing metadata from a Dataless SEED file).

  • Providing the filename/path to a local RESP file.

  • Providing a dictionary containing poles and zeros information. Be aware that this leads to wrong results if the instrument’s response is changing over the timespans that are added to the PPSD. Use with caution!

Note:

When using special_handling=”ringlaser” the applied processing steps are changed. Differentiation of data (converting velocity to acceleration data) will be omitted and a flat instrument response is assumed, leaving away response removal and only dividing by metadata[‘sensitivity’] specified in the provided metadata dictionary (other keys do not have to be present then). For scaling factors that are usually multiplied to the data remember to use the inverse as metadata[‘sensitivity’].

Parameters:
  • stats (Stats) – Stats of the station/instrument to process

  • metadata (Inventory or Parser or str or dict) – Response information of instrument. See above notes for details.

  • skip_on_gaps (bool, optional) – Determines whether time segments with gaps should be skipped entirely. [McNamara2004] merge gappy traces by filling with zeros. This results in a clearly identifiable outlier psd line in the PPSD visualization. Select skip_on_gaps=True for not filling gaps with zeros which might result in some data segments shorter than ppsd_length not used in the PPSD.

  • db_bins (tuple(int, int, int) or tuple(float, float, float)) – Specify the lower and upper boundary and the width of the db bins. The bin width might get adjusted to fit a number of equally spaced bins in between the given boundaries.

  • ppsd_length (float, optional) – Length of data segments passed to psd in seconds. In the paper by [McNamara2004] a value of 3600 (1 hour) was chosen. Longer segments increase the upper limit of analyzed periods but decrease the number of analyzed segments.

  • overlap (float, optional) – Overlap of segments passed to psd. Overlap may take values between 0 and 1 and is given as fraction of the length of one segment, e.g. ppsd_length=3600 and overlap=0.5 result in an overlap of 1800s of the segments.

  • special_handling (str, optional) – Switches on customized handling for data other than seismometer recordings. Can be one of: ‘ringlaser’ (no instrument correction, just division by metadata[“sensitivity”] of provided metadata dictionary), ‘hydrophone’ (no differentiation after instrument correction).

  • period_smoothing_width_octaves (float) – Determines over what period/frequency range the psd is smoothed around every central period/frequency. Given in fractions of octaves (default of 1 means the psd is averaged over a full octave at each central frequency).

  • period_step_octaves (float) – Step length on frequency axis in fraction of octaves (default of 0.125 means one smoothed psd value on the frequency axis is measured every 1/8 of an octave).

  • period_limits (tuple or list[float, float]) – Set custom lower and upper end of period range (e.g. (0.01, 100)). The specified lower end of period range will be set as the central period of the first bin (geometric mean of left/right edges of smoothing interval). At the upper end of the specified period range, no more additional bins will be added after the bin whose center frequency exceeds the given upper end for the first time.

PPSD.__init_subclass__()

This method is called when a class is subclassed.

The default implementation does nothing. It may be overridden to extend subclasses.

PPSD.__le__(value, /)

Return self<=value.

PPSD.__lt__(value, /)

Return self<value.

PPSD.__ne__(value, /)

Return self!=value.

PPSD.__new__(**kwargs)
PPSD.__reduce__()

Helper for pickle.

PPSD.__reduce_ex__(protocol, /)

Helper for pickle.

PPSD.__repr__()

Return repr(self).

PPSD.__setattr__(name, value, /)

Implement setattr(self, name, value).

PPSD.__sizeof__()

Size of object in memory, in bytes.

PPSD.__str__()

Return str(self).

PPSD.__subclasshook__()

Abstract classes can override this to customize issubclass().

This is invoked early on by abc.ABCMeta.__subclasscheck__(). It should return True, False or NotImplemented. If it returns NotImplemented, the normal algorithm is used. Otherwise, it overrides the normal algorithm (and the outcome is cached).