# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
'''
Multi-component waveform data model.
A :py:class:`Carpet` holds the samples of several components on a common
time axis in a single 2D array. Typical uses are the waveforms of a
multi-component sensor or of a sensor array, and spectrograms, where each
component is a frequency band (see :py:attr:`Carpet.component_axes`).
Carpets can be stored in files of the ``.rug`` format (see
:py:mod:`pyrocko.io.rug`) and are handled by Squirrel like other content
(see :py:meth:`pyrocko.squirrel.base.Squirrel.get_carpets`).
.. note::
This module is experimental. Its interface may still change and a warning
is shown when it is used.
'''
import logging
import warnings
import re
from collections import defaultdict
import numpy as num
import numpy.ma as ma
from scipy import signal
from . import trace, util
from .trace import Trace, AboveNyquist, _get_cached_filter_coeffs, costaper
from .guts import Object, Float, Timestamp, List, Int, Dict, String, Any
from .guts_array import Array
from .model.codes import CodesNSLCE
guts_prefix = 'pf'
logger = logging.getLogger('pyrocko.carpet')
EMPTY_CODES = CodesNSLCE()
def _normalize_mask(mask, shape):
'''
Check a mask and convert it to the canonical form.
Returns ``None`` for no mask (also for a mask without any masked samples)
or a boolean array of the given shape.
'''
if mask is None or mask is ma.nomask:
return None
mask = num.asarray(mask, dtype=bool)
if mask.shape != tuple(shape):
raise ValueError(
'Carpet: shape of mask %s does not match shape of data %s.' % (
mask.shape, tuple(shape)))
if not num.any(mask):
return None
return mask
def _split_masked(data, mask=None):
'''
Split a possibly masked array into plain data and mask.
The mask of a :py:class:`numpy.ma.MaskedArray` is combined with ``mask``.
Returns ``(data, mask)`` with ``data`` a plain :py:class:`numpy.ndarray`
(sharing the samples) or ``None``, and ``mask`` as returned by
:py:func:`_normalize_mask`.
'''
if data is None:
return None, mask
if isinstance(data, ma.MaskedArray):
data_mask = ma.getmask(data)
data = data.data
if data_mask is not ma.nomask:
mask = data_mask if mask is None \
else num.logical_or(data_mask, mask)
return data, _normalize_mask(mask, data.shape)
[docs]
class CarpetError(Exception):
'''
Base class for errors raised by :py:class:`Carpet` operations.
'''
pass
[docs]
class CarpetOverlapError(CarpetError):
'''
Raised by some :py:class:`Carpet` operations when carpets overlap but
should not (see :py:func:`check_overlaps`).
'''
pass
[docs]
class CarpetResampleError(CarpetError):
'''
Raised when :py:meth:`Carpet.resample_band` cannot be applied.
'''
pass
[docs]
class CarpetStringFiller:
'''
Placeholder values for :py:meth:`Carpet.fill_template`.
Gives access to the following placeholders by name:
``network``, ``station``, ``location``, ``channel``, ``extra``
Codes of the carpet (:py:attr:`Carpet.codes`).
``network_safe``, ..., ``extra_safe``, ``nslce_safe``, ``nslc_safe``
The same, with characters other than ``0-9``, ``A-Z``, ``a-z``, ``_``
and ``-`` replaced by ``_``, for use in file names. ``nslce_safe`` and
``nslc_safe`` join the codes with ``.``.
``network_dsafe``, ..., ``extra_dsafe``
Like the ``_safe`` variants, but ``_`` instead of an empty string,
for use as directory names.
``tmin``, ``tmax``
Start and end time as ``YYYY-MM-DD_HH-MM-SS``.
``tmin_ms``, ``tmax_ms``, ``tmin_us``, ``tmax_us``
The same, with milli- or microseconds.
``tmin_year``, ``tmin_month``, ``tmin_day``, ``tmin_jday``, ``tmax_...``
Parts of the start and end time (``jday``: day of the year).
``julianday``
Day of the year of the start time, without zero padding.
Any other name is looked up in ``additional``.
'''
def __init__(self, carpet, additional={}):
self.carpet = carpet
self.codes = carpet.codes
self.additional = additional
def __getitem__(self, k):
if k in ('network', 'station', 'location', 'channel', 'extra'):
return getattr(self.codes, k)
method = getattr(self, 'get_' + k, None)
if method:
return method()
return self.additional[k]
def _filename_safe(self, s):
return re.sub(r'[^0-9A-Za-z_-]', '_', s)
def get_network_safe(self):
return self._filename_safe(self.codes.network)
def get_station_safe(self):
return self._filename_safe(self.codes.station)
def get_location_safe(self):
return self._filename_safe(self.codes.location)
def get_channel_safe(self):
return self._filename_safe(self.codes.channel)
def get_extra_safe(self):
return self._filename_safe(self.codes.extra)
def get_nslce_safe(self):
return self._filename_safe('.'.join(self.codes))
def get_nslc_safe(self):
return self._filename_safe('.'.join(self.codes[:4]))
def get_network_dsafe(self):
return self._filename_safe(self.codes.network) or '_'
def get_station_dsafe(self):
return self._filename_safe(self.codes.station) or '_'
def get_location_dsafe(self):
return self._filename_safe(self.codes.location) or '_'
def get_channel_dsafe(self):
return self._filename_safe(self.codes.channel) or '_'
def get_extra_dsafe(self):
return self._filename_safe(self.codes.extra) or '_'
def get_tmin(self):
return util.time_to_str(
self.carpet.tmin, format='%Y-%m-%d_%H-%M-%S')
def get_tmax(self):
return util.time_to_str(
self.carpet.tmax, format='%Y-%m-%d_%H-%M-%S')
def get_tmin_ms(self):
return util.time_to_str(
self.carpet.tmin, format='%Y-%m-%d_%H-%M-%S.3FRAC')
def get_tmax_ms(self):
return util.time_to_str(
self.carpet.tmax, format='%Y-%m-%d_%H-%M-%S.3FRAC')
def get_tmin_us(self):
return util.time_to_str(
self.carpet.tmin, format='%Y-%m-%d_%H-%M-%S.6FRAC')
def get_tmax_us(self):
return util.time_to_str(
self.carpet.tmax, format='%Y-%m-%d_%H-%M-%S.6FRAC')
def get_tmin_year(self):
return util.time_to_str(self.carpet.tmin, format='%Y')
def get_tmin_month(self):
return util.time_to_str(self.carpet.tmin, format='%m')
def get_tmin_day(self):
return util.time_to_str(self.carpet.tmin, format='%d')
def get_tmax_year(self):
return util.time_to_str(self.carpet.tmax, format='%Y')
def get_tmax_month(self):
return util.time_to_str(self.carpet.tmax, format='%m')
def get_tmax_day(self):
return util.time_to_str(self.carpet.tmax, format='%d')
def get_julianday(self):
return str(util.julian_day_of_year(self.carpet.tmin))
def get_tmin_jday(self):
return util.time_to_str(self.carpet.tmin, format='%j')
def get_tmax_jday(self):
return util.time_to_str(self.carpet.tmax, format='%j')
[docs]
class CarpetStats(Object):
'''
Summary statistics of the data samples of a carpet.
Masked samples are ignored. Non-finite values are ignored in the
percentiles and counted in :py:gattr:`nan_count`.
'''
min = Float.T(help='Minimum value.')
p10 = Float.T(help='10th percentile.')
median = Float.T(help='Median.')
p90 = Float.T(help='90th percentile.')
max = Float.T(help='Maximum value.')
nan_count = Int.T(help='Number of non-finite values (NaN, inf).')
[docs]
@classmethod
def make(cls, data, mask=None):
'''
Compute statistics of a data array.
:param data:
Data samples.
:type data:
:py:class:`numpy.ndarray`
:param mask:
Mask of samples to ignore (``True`` = masked).
:type mask:
:py:class:`numpy.ndarray` of :py:class:`bool`
:returns:
New statistics object.
:rtype:
:py:class:`CarpetStats`
'''
data = num.asarray(data, dtype=float)
finite = num.isfinite(data)
if mask is not None:
nan_count = num.sum(~finite & ~mask)
finite &= ~mask
else:
nan_count = num.sum(~finite)
values = data[finite]
if values.size == 0:
min = p10 = median = p90 = max = num.nan
else:
min, p10, median, p90, max = num.percentile(
values, [0., 10., 50., 90., 100.])
return cls(
min=min,
p10=p10,
median=median,
p90=p90,
max=max,
nan_count=nan_count)
def _smooth3(data, mask):
'''
Weighted average over three samples, weights ``(0.25, 0.5, 0.25)``,
leaving out masked samples and samples beyond the ends.
'''
if mask is None:
weight = num.ones(data.shape)
else:
weight = (~mask).astype(float)
def conv(x):
y = 0.5 * x
y[:, 1:] += 0.25 * x[:, :-1]
y[:, :-1] += 0.25 * x[:, 1:]
return y
sum_weight = conv(weight)
out = num.ones_like(data, dtype=float)
num.divide(conv(data * weight), sum_weight, out=out, where=sum_weight > 0)
return out
[docs]
class Carpet(Object):
'''
Container for multi-component waveforms with common time span and sampling.
Instances of this class can be used to efficiently represent
multi-component waveforms of a single sensor or of a sensor array, or
spectrograms. The data samples are stored in a single 2D array where the
first index runs over components and the second index over time. Metadata
contains sampling interval, start time and
:py:class:`~pyrocko.model.codes.CodesNSLCE` identifiers for the carpet as
a whole (:py:gattr:`codes`) and for each component
(:py:gattr:`component_codes`). Further coordinates of the components,
e.g. the frequencies of a spectrogram, can be attached as
:py:gattr:`component_axes`.
The samples are held in :py:gattr:`data`, a plain
:py:class:`numpy.ndarray`. Missing or invalid samples, e.g. gaps or
samples polluted by artifacts, can be marked in the optional boolean
:py:gattr:`mask` (``True`` = masked). Without a mask (``None``), all
samples are valid. The values of masked samples are undefined. Methods
carry the mask along: filters treat masked samples as zeros and keep them
masked, and :py:meth:`bleed` extends the masked regions to hide edge
effects. :py:attr:`masked_data` gives a :py:class:`numpy.ma.MaskedArray`
view, for use with :py:mod:`numpy.ma` functions.
Some carpets, e.g. spectrograms, mark missing values with NaN instead.
:py:meth:`mask_invalid` converts NaN (and inf) into masked samples.
A carpet can also be created without data ("dataless"), e.g. to hold
metadata only. In this case, ``component_codes`` and ``nsamples`` must be
given.
Most methods returning a carpet create a new one, which may share the
data with the original (:py:meth:`chop`, :py:meth:`chopper`,
:py:meth:`crop`). The filter methods (:py:meth:`lowpass`,
:py:meth:`highpass`, :py:meth:`smooth`, :py:meth:`whiten`, ...) modify
the carpet in place.
:param traces:
If given, construct the carpet from these single-component waveforms,
see ``assemble``. ``data``, ``component_codes``, ``tmin`` and
``deltat`` are then taken from the traces.
:type traces:
:py:class:`list` of :py:class:`~pyrocko.trace.Trace`
:param assemble:
How to combine ``traces``: ``'concatenate'`` (default) stacks the
traces in the given order, one component per trace. They must have
the same time span, sampling interval and data type (see
:py:func:`~pyrocko.trace.get_traces_data_as_array`). ``'merge'``
combines traces with the same codes into one component each (sorted
by codes), covering the joint time span of all traces. Samples not
covered by any trace are masked (see
:py:func:`~pyrocko.trace.merge_traces_data_as_array`). The traces must
have the same sampling interval and sampling instants.
:type assemble:
str
:param data:
Data samples, indexed as ``(icomponent, isample)``. If a
:py:class:`numpy.ma.MaskedArray` is given, its mask is combined with
``mask``.
:type data:
:py:class:`numpy.ndarray` or :py:class:`numpy.ma.MaskedArray`
:param mask:
Mask of invalid samples, same shape as ``data`` (``True`` = masked).
:type mask:
:py:class:`numpy.ndarray` of :py:class:`bool`
:param nsamples:
Number of samples. Required for a dataless carpet, otherwise taken
from ``data`` (and checked, if given).
:type nsamples:
int
:param codes:
Codes identifying the carpet as a whole.
:type codes:
:py:class:`~pyrocko.model.codes.CodesNSLCE`
:param component_codes:
Codes of the components, one per row of ``data``. Default: empty
codes for all components.
:type component_codes:
:py:class:`list` of :py:class:`~pyrocko.model.codes.CodesNSLCE`
:param component_axes:
Auxiliary coordinate axes of the components, by name, e.g.
``{'frequency': frequencies}``. Each array must have one value per
component.
:type component_axes:
:py:class:`dict` of :py:class:`str` to :py:class:`numpy.ndarray`
:param tmin:
Time of the first sample [s].
:type tmin:
:py:func:`~pyrocko.util.get_time_float`
:param deltat:
Sampling interval [s].
:type deltat:
float
:param meta:
Additional free-form metadata.
:type meta:
:py:class:`dict`
:raises:
:py:exc:`ValueError` if the arguments are inconsistent, e.g. if the
number of component codes does not match the number of rows of
``data``.
'''
codes = CodesNSLCE.T(
default=CodesNSLCE.D(),
help='Codes identifying the carpet as a whole.')
component_codes = List.T(
CodesNSLCE.T(),
help='List of codes identifying the components.')
component_axes = Dict.T(
String.T(),
Array.T(shape=(None,), serialize_as='base64+meta'),
help='Auxiliary component coordinate axes, by name, e.g. the '
'frequencies of a spectrogram. Each axis has one value per '
'component.')
nsamples = Int.T(
help='Number of samples.')
data = Array.T(
optional=True,
shape=(None, None),
serialize_as='base64+meta',
help='Array containing the data samples indexed as '
'``(icomponent, isample)``.')
mask = Array.T(
optional=True,
shape=(None, None),
dtype=bool,
serialize_as='base64+meta',
help='Mask of invalid samples, same shape as :py:gattr:`data` '
'(``True`` = masked). ``None`` if all samples are valid.')
tmin = Timestamp.T(
default=Timestamp.D('1970-01-01 00:00:00'),
help='Time of first sample [s].')
deltat = Float.T(
default=1.0,
help='Sampling interval [s].')
meta = Dict.T(
String.T(), Any.T(),
help='Additional free-form metadata.')
stats__ = CarpetStats.T(
optional=True,
help='Summary statistics of the data samples. Computed on first '
'access and cached; the cache is not updated when the data '
'changes.')
def __init__(
self,
traces=None,
assemble='concatenate',
data=None,
nsamples=None,
codes=None,
component_codes=None,
component_axes=None,
tmin=None,
deltat=None,
mask=None,
meta=None,
stats=None):
util.experimental_feature_used('pyrocko.carpet')
if data is not None and not isinstance(data, num.ndarray):
data = self.T.get_property('data').regularize_extra(data)
if mask is not None and not isinstance(mask, num.ndarray):
mask = self.T.get_property('mask').regularize_extra(mask)
if traces is not None:
if len(traces) == 0:
data = num.zeros((0, 0))
else:
if assemble == 'merge':
data, component_codes, tmin, deltat \
= trace.merge_traces_data_as_array(traces)
elif assemble == 'concatenate':
data = trace.get_traces_data_as_array(traces)
component_codes = [tr.codes for tr in traces]
tmin = traces[0].tmin
deltat = traces[0].deltat
data, mask = _split_masked(data, mask)
if nsamples is not None and data is not None \
and data.shape[1] != nsamples:
raise ValueError(
'Carpet construction: mismatch between expected number of '
'samples and number of samples in data array.')
if data is None:
if component_codes is None or nsamples is None:
raise ValueError(
'Dataless Carpet must be constructed with '
'`component_codes` ' 'and `nsamples` set.')
if mask is not None:
raise ValueError(
'Carpet construction: mask given without data.')
self.ncomponents = len(component_codes)
else:
self.ncomponents, nsamples = data.shape
if component_codes is None:
component_codes = [CodesNSLCE()] * self.ncomponents
if len(component_codes) != self.ncomponents:
raise ValueError(
'Carpet construction: mismatch between number of traces '
'and number of component codes given.')
if deltat is None:
deltat = self.T.get_property('deltat').default()
if tmin is None:
tmin = self.T.get_property('tmin').default()
if codes is None:
codes = CodesNSLCE()
if component_axes is None:
component_axes = {}
if meta is None:
meta = {}
Object.__init__(
self,
codes=codes,
component_codes=component_codes,
component_axes=component_axes,
data=data,
mask=mask,
tmin=tmin,
nsamples=nsamples,
deltat=deltat,
meta=meta,
stats=stats)
@property
def stats(self):
if self._stats is None:
self._stats = CarpetStats.make(self.data, self.mask)
return self._stats
@stats.setter
def stats(self, stats):
self._stats = stats
@property
def masked_data(self):
'''
The data as :py:class:`numpy.ma.MaskedArray` (read-only).
Shares the data samples with the carpet. Changing the mask of the
returned array does not change the carpet's :py:gattr:`mask`, use
:py:meth:`set_mask` for that.
'''
return ma.MaskedArray(
self.data,
mask=ma.nomask if self.mask is None else self.mask.copy())
[docs]
def get_mask(self):
'''
Get the mask, also if all samples are valid.
:returns:
Boolean array of the shape of :py:gattr:`data`, ``True`` for
masked samples. A new array of ``False`` values if the carpet has
no mask.
:rtype:
:py:class:`numpy.ndarray`
'''
if self.mask is None:
return num.zeros(self.data.shape, dtype=bool)
else:
return self.mask
[docs]
def set_mask(self, mask):
'''
Set the mask.
:param mask:
Boolean array of the shape of :py:gattr:`data` (``True`` =
masked), or ``None`` to mark all samples as valid. A mask without
any masked samples is stored as ``None``.
:type mask:
:py:class:`numpy.ndarray`
'''
self.mask = _normalize_mask(mask, self.data.shape)
[docs]
def mask_invalid(self):
'''
Mask non-finite samples (NaN, inf), in place.
The samples already masked stay masked.
'''
if self.data.dtype.kind in 'fc':
invalid = ~num.isfinite(self.data)
if num.any(invalid):
self.set_mask(self.get_mask() | invalid)
@property
def summary_component_codes(self):
'''
Short textual representation of the component codes, e.g. first and
last, for :py:attr:`summary`.
'''
if all(codes == EMPTY_CODES for codes in self.component_codes):
return ''
elif len(self.component_codes) == 1:
return str(self.component_codes[0])
elif len(self.component_codes) == 2:
return '%s, %s' % (
self.component_codes[0],
self.component_codes[-1])
else:
return '%s, ..., %s' % (
self.component_codes[0],
self.component_codes[-1])
@property
def summary_entries(self):
'''
Entries of :py:attr:`summary`, as a tuple of strings: class name,
codes, number of components, number of samples, data type, sampling
interval, start and end time, and component codes.
'''
return (
self.__class__.__name__,
str(self.codes),
str(self.ncomponents),
str(self.nsamples),
str(self.data.dtype) if self.data is not None else '<unknown>',
str(self.deltat),
util.time_to_str(self.tmin),
util.time_to_str(self.tmax),
self.summary_component_codes)
@property
def summary(self):
'''
One-line textual summary of the carpet's metadata attributes.
'''
return util.fmt_summary(
self.summary_entries, (10, 20, 5, 7, 10, 10, 25, 25, 50))
def __len__(self):
'''
Get number of components.
'''
return self.ncomponents
def __getitem__(self, i):
'''
Get single component waveform (shared data).
See :py:meth:`get_trace`. Iterating over a carpet yields its
components as traces.
:param i:
Component index.
:type i:
int
:rtype:
:py:class:`~pyrocko.trace.Trace`
'''
return self.get_trace(i)
[docs]
def itime(self, time, out_of_bounds='raise'):
'''
Get index of the sample nearest to a given time.
:param time:
Time [s].
:type time:
:py:func:`~pyrocko.util.get_time_float`
:param out_of_bounds:
What to do if ``time`` is outside the carpet's time span:
``'raise'`` raises :py:exc:`ValueError`, ``'clip'`` returns the
index of the first or last sample, ``'ignore'`` returns the
index anyway (it may then be negative or too large).
:type out_of_bounds:
str
:returns:
Sample index.
:rtype:
int
'''
itime = int(round((time - self.tmin) / self.deltat))
if itime < 0 or itime > self.nsamples-1:
if out_of_bounds == 'raise':
raise ValueError('Out of bounds.')
elif out_of_bounds == 'clip':
itime = min(max(0, itime), self.nsamples-1)
elif out_of_bounds == 'ignore':
pass
else:
raise ValueError(
'Invalid choice for `out_of_bounds`: %s' % out_of_bounds)
return itime
[docs]
def fill_template(self, template, **additional):
'''
Fill placeholders in a template string with the carpet's attributes.
Used e.g. to create file names. Placeholders can be given in
``str.format`` style, e.g. ``'{station}_{tmin}.rug'``, or in
``%`` style, e.g. ``'%(station)s_%(tmin)s.rug'``. See
:py:class:`CarpetStringFiller` for the available names.
:param template:
Template string.
:type template:
str
:param additional:
Values for additional placeholders.
:returns:
Filled string. A template without placeholders is returned
unchanged.
:rtype:
str
'''
if '{' in template:
return template.format_map(
CarpetStringFiller(self, additional=additional))
elif '%' in template:
return template % CarpetStringFiller(self, additional=additional)
else:
return template
[docs]
def copy(self, data='copy'):
'''
Create a copy.
Codes, component codes, component axes, start time, sampling
interval and :py:gattr:`meta` are copied (shallow).
:param data:
``'copy'`` to deeply copy the data and mask, or ``'reference'`` to
create a shallow copy, referencing the original data and mask, or
``'drop'`` to create a dataless copy. An array given here is used
as the data of the copy instead, with the mask of a
:py:class:`numpy.ma.MaskedArray` or else no mask.
:type data:
str or :py:class:`numpy.ndarray`
:rtype:
:py:class:`Carpet`
'''
mask = None
if isinstance(data, str):
if data == 'drop':
data = None
elif data == 'copy':
data = self.data.copy()
mask = None if self.mask is None else self.mask.copy()
elif data == 'reference':
data = self.data
mask = self.mask
else:
raise ValueError('Carpet.copy(): invalid argument for `data`.')
elif not isinstance(data, num.ndarray):
raise ValueError('Carpet.copy(): invalid argument for `data`.')
return Carpet(
data=data,
mask=mask,
codes=self.codes,
component_codes=list(self.component_codes),
component_axes=self.component_axes,
nsamples=self.nsamples if data is None else None,
tmin=self.tmin,
deltat=self.deltat,
meta=dict(self.meta))
[docs]
def chop(self, tmin, tmax, snap=(round, round)):
'''
Cut out a time span.
Unlike :py:meth:`pyrocko.trace.Trace.chop`, this does not modify the
carpet but returns a new one, sharing the data with the original.
The new carpet contains the samples with indices from
``snap[0]((tmin - self.tmin) / deltat)`` up to, but excluding,
``snap[1]((tmax - self.tmin) / deltat)``, limited to the available
samples. With the default rounding, ``tmax`` is therefore exclusive:
``chop(t, t + n*deltat)`` gives ``n`` samples.
:param tmin:
Start time [s].
:type tmin:
:py:func:`~pyrocko.util.get_time_float`
:param tmax:
End time [s].
:type tmax:
:py:func:`~pyrocko.util.get_time_float`
:param snap:
Rounding functions to convert ``tmin`` and ``tmax`` to sample
indices, e.g. ``(math.floor, math.ceil)`` to include partially
covered samples.
:type snap:
:py:class:`tuple` of two functions
:raises:
:py:exc:`~pyrocko.trace.NoData` if no samples remain.
:rtype:
:py:class:`Carpet`
'''
istart = int(snap[0]((tmin - self.tmin) / self.deltat))
iend = int(snap[1]((tmax - self.tmin) / self.deltat))
istart = max(0, istart)
iend = min(self.nsamples, iend)
if istart >= iend:
raise trace.NoData()
return self._new_time_window(istart, iend)
def _new_time_window(self, istart, iend):
return Carpet(
data=self.data[:, istart:iend],
mask=None if self.mask is None else self.mask[:, istart:iend],
codes=self.codes,
component_codes=self.component_codes,
tmin=self.tmin + istart * self.deltat,
deltat=self.deltat,
component_axes=self.component_axes)
[docs]
def chopper(self, tinc):
'''
Iterate over consecutive time windows.
The windows start at :py:gattr:`tmin` and have a length of
``floor(tinc / deltat)`` samples each. The last window may be
shorter. The yielded carpets share the data with the original.
:param tinc:
Window length [s].
:type tinc:
float
:yields:
:py:class:`Carpet` objects.
'''
nwindows = int(num.floor((self.tmax - self.tmin) / tinc)) + 1
nsamples = int(num.floor(tinc / self.deltat))
for iwindow in range(nwindows):
istart = int(num.floor((iwindow * tinc) / self.deltat))
iend = istart + nsamples
yield self._new_time_window(istart, iend)
[docs]
def overlaps(self, *args):
'''
Check if this carpet overlaps with a given time span or other carpet.
Call as ``overlaps(tmin, tmax)`` or ``overlaps(other_carpet)``. Time
spans touching at a single instant count as overlapping. Only times
are compared, not codes.
:returns:
``True`` if the time spans overlap.
:rtype:
bool
'''
if len(args) == 1 and isinstance(args[0], Carpet):
tmin, tmax = args[0].tmin, args[0].tmax
elif len(args) == 2:
tmin, tmax = args
return not (tmax < self.tmin or self.tmax < tmin)
@property
def tmax(self):
'''
End time (time of last sample, read-only).
'''
return self.tmin + (self.nsamples - 1) * self.deltat
@property
def times(self):
'''
Times of the samples [s] (read-only).
:rtype:
:py:class:`numpy.ndarray`
'''
return self.tmin + num.arange(self.nsamples) * self.deltat
[docs]
def get_trace(self, i, span=slice(None)):
'''
Get single component waveform (shared data).
The trace refers to the raw sample values: masked samples are
included with whatever value they hold. Use
:py:meth:`iter_valid_traces` to get the unmasked parts only.
:param i:
Component index.
:type i:
int
:param span:
Range of sample indices to include, default: all.
:type span:
:py:class:`slice`
:rtype:
:py:class:`~pyrocko.trace.Trace`
'''
network, station, location, channel, extra = self.component_codes[i]
return Trace(
network=network,
station=station,
location=location,
channel=channel,
extra=extra,
tmin=self.tmin + (span.start or 0) * self.deltat,
deltat=self.deltat,
ydata=self.data[i, span])
[docs]
def iter_valid_traces(self):
'''
Iterate over the unmasked parts of the components, as traces.
Each contiguous range of unmasked samples of a component gives one
trace (shared data).
:yields:
:py:class:`~pyrocko.trace.Trace` objects.
'''
if self.mask is None:
yield from self
else:
for irow in range(self.ncomponents):
for span in self._valid_spans(irow):
yield self.get_trace(irow, span)
def _valid_spans(self, irow):
if self.mask is None:
return [slice(0, self.nsamples)]
return ma.clump_unmasked(
ma.MaskedArray(self.data[irow], mask=self.mask[irow]))
[docs]
def get_traces(self):
'''
Get all components as traces, including masked samples.
See :py:meth:`get_trace`.
:rtype:
:py:class:`list` of :py:class:`~pyrocko.trace.Trace`
'''
return list(self)
[docs]
def get_valid_traces(self):
'''
Get the unmasked parts of the components as traces.
See :py:meth:`iter_valid_traces`.
:rtype:
:py:class:`list` of :py:class:`~pyrocko.trace.Trace`
'''
return list(self.iter_valid_traces())
[docs]
def get_component_axis(self, name=None):
'''
Get coordinates of the components.
:param name:
Name of an entry in :py:gattr:`component_axes`, or ``None`` to get
the component indices.
:type name:
str
:rtype:
:py:class:`numpy.ndarray`
'''
if name is None:
return num.arange(self.ncomponents)
else:
return self.component_axes[name]
[docs]
def mpl_draw(
self,
axes,
component_axis=None,
fslice=slice(None, None),
**kwargs):
'''
Draw the carpet into Matplotlib axes, with
:py:meth:`matplotlib.axes.Axes.pcolormesh`.
Time runs along the x-axis, the components along the y-axis. Masked
samples are left blank.
:param axes:
Axes to draw into.
:type axes:
:py:class:`matplotlib.axes.Axes`
:param component_axis:
Coordinates for the y-axis: ``None`` for the component indices,
``'default'`` for the only entry of :py:gattr:`component_axes` if
there is exactly one (indices otherwise), or the name of an entry
of :py:gattr:`component_axes`.
:type component_axis:
str
:param fslice:
Components to draw.
:type fslice:
:py:class:`slice`
:param kwargs:
Passed to :py:meth:`~matplotlib.axes.Axes.pcolormesh`.
:returns:
Object returned by :py:meth:`~matplotlib.axes.Axes.pcolormesh`.
'''
if component_axis is None:
ys = num.arange(self.ncomponents)
elif component_axis == 'default':
if len(self.component_axes) == 1:
ys = list(self.component_axes.values())[0]
else:
ys = num.arange(self.ncomponents)
else:
ys = self.component_axes[component_axis]
return axes.pcolormesh(
self.times,
ys[fslice],
self.masked_data[fslice, :],
**kwargs)
[docs]
def plot(
self,
component_axis='default',
component_axis_scale='linear',
fslice=slice(None, None),
path=None, **kwargs):
'''
Plot the carpet in a new Matplotlib figure.
See :py:meth:`mpl_draw` for ``component_axis``, ``fslice`` and
``kwargs``.
:param component_axis_scale:
Scale of the y-axis, e.g. ``'linear'`` or ``'log'``.
:type component_axis_scale:
str
:param path:
If given, save the figure to this file instead of showing it.
:type path:
str
'''
from pyrocko import plot
import matplotlib.pyplot as plt
fig = plt.figure(figsize=plot.mpl_papersize('a4', 'landscape'))
axes = fig.add_subplot(1, 1, 1)
self.mpl_draw(
axes,
component_axis=component_axis,
fslice=fslice, **kwargs)
axes.set_yscale(component_axis_scale)
plot.mpl_time_axis(axes)
if path is None:
plt.show()
else:
fig.savefig(path)
[docs]
def snuffle(self, what='valid'):
'''
Show the components in Snuffler.
:param what:
``'valid'`` to show only the unmasked parts (see
:py:meth:`get_valid_traces`), ``'raw'`` to show all samples.
:type what:
str
'''
assert what in ('valid', 'raw')
if what == 'valid':
trace.snuffle(self.get_valid_traces())
else:
trace.snuffle(list(self))
[docs]
def bleed(self, t):
'''
Extend masked regions and mask the edges, in place.
Masks ``t`` seconds at both ends of every contiguous range of
unmasked samples, and at the start and end of the carpet. Useful to
hide edge effects, e.g. of filters, near gaps.
:param t:
Duration to mask [s]. Nothing is done if it is shorter than one
sample.
:type t:
float
'''
nt = int(num.round(abs(t)/self.deltat))
if nt < 1:
return
mask = self.get_mask().copy()
for irow in range(self.ncomponents):
for span in self._valid_spans(irow):
mask[irow, span.start:span.start+nt] = True
mask[irow, max(0, span.stop-nt):span.stop] = True
mask[:, :nt] = True
mask[:, -nt:] = True
self.set_mask(mask)
[docs]
def set_data(self, data):
'''
Replace the data samples, keeping the shape.
:param data:
New data samples, with the same shape as the current ones. If a
plain :py:class:`numpy.ndarray` is given, the current mask is
kept. The mask of a :py:class:`numpy.ma.MaskedArray` replaces the
current mask.
:type data:
:py:class:`numpy.ndarray` or :py:class:`numpy.ma.MaskedArray`
'''
if data is self.data:
return
if data.shape != self.data.shape:
raise ValueError(
'Carpet.set_data(): shape of new data %s differs from shape '
'of current data %s.' % (data.shape, self.data.shape))
if isinstance(data, ma.MaskedArray):
self.data, self.mask = _split_masked(data)
else:
self.data = data
[docs]
def apply(self, f):
'''
Apply a function to the data samples, in place.
:param f:
Function taking and returning a data array of the same shape,
see :py:meth:`set_data`. It gets the plain data array; masked
samples hold undefined values.
:type f:
callable
'''
self.set_data(f(self.data))
[docs]
def reduce(self, f, component_codes):
'''
Combine the components with a function, in place.
Example: ``carpet.reduce(lambda data: data.sum(axis=0), codes)``
replaces all components by their sum.
:param f:
Function taking the data array and returning a new one with the
same number of samples, either 2D ``(icomponent, isample)`` or
1D (a single component). If the carpet has a mask, ``f`` gets
:py:attr:`masked_data`, so that :py:mod:`numpy.ma` semantics
apply (e.g. ``sum`` skips masked samples), and the mask of the
result becomes the new mask.
:type f:
callable
:param component_codes:
Codes of the new components.
:type component_codes:
:py:class:`~pyrocko.model.codes.CodesNSLCE` or :py:class:`list`
of these
:py:gattr:`component_axes` is kept if the number of components
stays the same, and cleared otherwise.
'''
data = f(self.data if self.mask is None else self.masked_data)
if data.ndim == 1:
data = data[num.newaxis, :]
if isinstance(component_codes, CodesNSLCE):
component_codes = [component_codes]
if data.ndim != 2 or data.shape[1] != self.nsamples \
or len(component_codes) != data.shape[0]:
raise ValueError(
'Carpet.reduce(): result of shape %s does not match number '
'of samples (%i) and number of component codes (%i).' % (
data.shape, self.nsamples, len(component_codes)))
if data.shape[0] != self.ncomponents:
self.component_axes = {}
self.data, self.mask = _split_masked(data)
self.component_codes = component_codes
self.ncomponents = data.shape[0]
[docs]
def nyquist_check(
self,
frequency,
intro='Corner frequency',
warn=True,
raise_exception=False):
'''
Check if a given frequency is above the Nyquist frequency of the
carpet.
:param frequency:
Frequency to check [Hz].
:type frequency:
float
:param intro:
String used to introduce the warning/error message.
:type intro:
str
:param warn:
Whether to emit a warning message.
:type warn:
bool
:param raise_exception:
Whether to raise :py:exc:`~pyrocko.trace.AboveNyquist`.
:type raise_exception:
bool
'''
if frequency >= 0.5/self.deltat:
message = '%s (%g Hz) is equal to or higher than nyquist ' \
'frequency (%g Hz). (%s)' \
% (intro, frequency, 0.5/self.deltat, self.summary)
if warn:
logger.warning(message)
if raise_exception:
raise AboveNyquist(message)
[docs]
def lfilter(self, b, a, demean=True):
'''
Filter all components with :py:func:`scipy.signal.lfilter`, in place.
Sample data is converted to type :py:class:`float`, possibly demeaned
and filtered along the time axis using :py:func:`scipy.signal.lfilter`.
Masked samples are set to zero before filtering and stay masked. Use
:py:meth:`bleed` afterwards to mask the filter's transients at gaps.
:param b:
Numerator coefficients.
:type b:
:py:class:`numpy.ndarray`
:param a:
Denominator coefficients.
:type a:
:py:class:`numpy.ndarray`
:param demean:
Subtract the mean of each component (of the unmasked samples)
before filtering.
:type demean:
bool
'''
def filt(data):
data = self._get_filled_float_data(demean)
return signal.lfilter(b, a, data)
self.apply(filt)
def _get_filled_float_data(self, demean=False):
'''
Data as float array, masked samples set to zero, optionally demeaned
(mean of the unmasked samples).
'''
data = self.data.astype(num.float64)
if demean:
if self.mask is None:
data -= num.mean(data, axis=1)[:, num.newaxis]
else:
mean = ma.MaskedArray(data, mask=self.mask).mean(axis=1)
data -= ma.filled(mean, 0.0)[:, num.newaxis]
if self.mask is not None:
data[self.mask] = 0.0
return data
[docs]
def lowpass(self, order, corner, nyquist_warn=True,
nyquist_exception=False, demean=True):
'''
Filter all components using a Butterworth lowpass, in place.
Sample data is converted to type :py:class:`float`, possibly demeaned
and filtered using :py:func:`scipy.signal.lfilter`. Filter coefficients
are generated with :py:func:`scipy.signal.butter`.
:param order:
Order of the filter.
:type order:
int
:param corner:
Corner frequency of the filter [Hz].
:type corner:
float
:param demean:
Subtract mean before filtering.
:type demean:
bool
:param nyquist_warn:
Warn if corner frequency is greater than Nyquist frequency.
:type nyquist_warn:
bool
:param nyquist_exception:
Raise :py:exc:`pyrocko.trace.AboveNyquist` if corner frequency is
greater than Nyquist frequency.
:type nyquist_exception:
bool
'''
self.nyquist_check(
corner, 'Corner frequency of lowpass', nyquist_warn,
nyquist_exception)
(b, a) = _get_cached_filter_coeffs(
order, [corner*2.0*self.deltat], btype='low')
if len(a) != order+1 or len(b) != order+1:
logger.warning(
'Erroneous filter coefficients returned by '
'scipy.signal.butter(). Should downsample before filtering.')
self.lfilter(b, a, demean=demean)
[docs]
def highpass(self, order, corner, nyquist_warn=True,
nyquist_exception=False, demean=True):
'''
Filter all components using a Butterworth highpass, in place.
Sample data is converted to type :py:class:`float`, possibly demeaned
and filtered using :py:func:`scipy.signal.lfilter`. Filter coefficients
are generated with :py:func:`scipy.signal.butter`.
:param order:
Order of the filter.
:type order:
int
:param corner:
Corner frequency of the filter [Hz].
:type corner:
float
:param demean:
Subtract mean before filtering.
:type demean:
bool
:param nyquist_warn:
Warn if corner frequency is greater than Nyquist frequency.
:type nyquist_warn:
bool
:param nyquist_exception:
Raise :py:exc:`~pyrocko.trace.AboveNyquist` if corner frequency is
greater than Nyquist frequency.
:type nyquist_exception:
bool
'''
self.nyquist_check(
corner, 'Corner frequency of highpass', nyquist_warn,
nyquist_exception)
(b, a) = _get_cached_filter_coeffs(
order, [corner*2.0*self.deltat], btype='high')
if len(a) != order+1 or len(b) != order+1:
logger.warning(
'Erroneous filter coefficients returned by '
'scipy.signal.butter(). Should downsample before filtering.')
self.lfilter(b, a, demean=demean)
[docs]
def fold(self, t, method=num.mean):
'''
Downsample by combining blocks of samples.
Blocks of ``nfold = round(t / deltat)`` consecutive samples are
combined into one with ``method``. The blocks are aligned to
multiples of the new sampling interval ``nfold * deltat`` in absolute
time, so that folded carpets of different time spans fit together.
Incomplete blocks at the start and end are dropped. If the carpet has
a mask, ``method`` gets a :py:class:`numpy.ma.MaskedArray`: with
:py:func:`numpy.mean` (or other functions supporting masked arrays),
masked samples are skipped and a block is only masked if all its
samples are masked.
:param t:
New sampling interval [s], rounded to a multiple of
:py:gattr:`deltat`.
:type t:
float
:param method:
Function to combine the samples of a block, called as
``method(data, axis=2)`` on an array of shape
``(ncomponents, nblocks, nfold)``, e.g. :py:func:`numpy.mean`
(default) or :py:func:`numpy.max`.
:type method:
callable
:raises:
:py:exc:`~pyrocko.trace.NoData` if no complete block remains.
:returns:
New carpet.
:rtype:
:py:class:`Carpet`
'''
nfold = max(1, int(num.round(t / self.deltat)))
itmin = int(round(self.tmin / self.deltat))
itmax = itmin + self.nsamples
ibmin = (((itmin - 1) // nfold) + 1) * nfold
ibmax = (itmax // nfold) * nfold
ibmin -= itmin
ibmax -= itmin
if ibmin >= ibmax:
raise trace.NoData()
data = self.data if self.mask is None else self.masked_data
data = data[:, ibmin:ibmax].reshape((
self.ncomponents,
(ibmax - ibmin) // nfold,
nfold))
with warnings.catch_warnings():
warnings.simplefilter('ignore')
data = method(data, axis=2)
return Carpet(
data=data,
codes=self.codes,
component_axes=self.component_axes,
component_codes=self.component_codes,
tmin=self.tmin + self.deltat * ibmin,
deltat=nfold*self.deltat)
[docs]
def smooth(self, t, window=num.hanning):
'''
Smooth along the time axis with a moving weighted average, in place.
Each sample is replaced by the average of the valid samples around
it, weighted with the window. Masked samples and samples beyond the
ends of the carpet are left out of the average, and the weights of
the remaining samples are renormalized. So, a constant signal is left
unchanged, also next to gaps and at the ends. The mask is not changed.
:param t:
Window length [s], rounded to an odd number of samples.
:type t:
float
:param window:
Function returning the window for a given number of samples, e.g.
:py:func:`numpy.hanning` (default).
:type window:
callable
'''
self.set_data(self._smoothed(
self._get_filled_float_data(), t, window))
def _smoothed(self, data, t, window):
'''
Mask-aware moving weighted average of ``data`` (masked samples must
be set to zero), see :py:meth:`smooth`.
'''
n = (int(num.round(t / self.deltat)) // 2) * 2 + 1
taper = num.asarray(window(n), dtype=float)
nh = n // 2
# long enough that the convolution does not wrap around
ntrans = trace.nextpow2(self.nsamples + n)
taper_pad = num.zeros(ntrans)
taper_pad[:nh+1] = taper[nh:]
if nh:
taper_pad[ntrans-nh:] = taper[:nh]
taper_fd = num.fft.rfft(taper_pad)
def convolve(x):
return num.fft.irfft(
num.fft.rfft(x, ntrans, axis=1) * taper_fd[num.newaxis, :],
ntrans, axis=1)[:, :self.nsamples]
if self.mask is None:
valid = num.ones((1, self.nsamples))
else:
valid = (~self.mask).astype(float)
weight = convolve(valid)
smoothed = convolve(data)
# The weight is at least taper[nh] at valid samples. Where it is
# (numerically) zero, all samples in the window are masked.
ok = weight > 1e-6 * num.max(taper)
out = num.zeros_like(smoothed)
num.divide(smoothed, weight, out=out, where=ok)
return out
[docs]
def filter_frequency_cos_taper(self, f1, f2, f3, f4):
'''
Bandpass filter with a cosine taper in the frequency domain, in place.
The spectrum is multiplied by a taper which rises from 0 at ``f1`` to
1 at ``f2``, and falls from 1 at ``f3`` to 0 at ``f4`` (see
:py:func:`pyrocko.trace.costaper`). See :py:meth:`apply_via_fft`.
:param f1:
Frequency [Hz] where the taper starts to rise from 0.
:type f1:
float
:param f2:
Frequency [Hz] where the taper reaches 1.
:type f2:
float
:param f3:
Frequency [Hz] where the taper starts to fall from 1.
:type f3:
float
:param f4:
Frequency [Hz] where the taper reaches 0.
:type f4:
float
'''
def multiply_taper(frequency_delta, ntrans, spectrum):
taper = costaper(
f1, f2, f3, f4, spectrum.shape[1], frequency_delta)
spectrum *= taper[num.newaxis, :]
return spectrum
self.apply_via_fft(
multiply_taper)
[docs]
def whiten(self, deltaf, window=num.hanning):
'''
Spectral whitening, in place.
The spectrum of each component is divided by its amplitude spectrum,
smoothed over a frequency band of width ``deltaf``. See
:py:meth:`apply_via_fft`. The smoothing window is normalized to unit
sum, so that the resulting amplitude spectrum is close to one.
:param deltaf:
Width of the smoothing window [Hz].
:type deltaf:
float
:param window:
Function returning the smoothing window for a given number of
samples.
:type window:
callable
'''
def smooth(frequency_delta, ntrans, spectrum):
n = (int(num.round(deltaf / frequency_delta)) // 2) * 2 + 1
taper = num.asarray(window(n), dtype=float)
taper /= num.sum(taper)
amp_spec = num.abs(spectrum)
amp_spec_smooth = signal.fftconvolve(
amp_spec, taper[num.newaxis, :], mode='same', axes=1)
spectrum /= amp_spec_smooth
return spectrum
self.apply_via_fft(smooth)
[docs]
def normalize(self, deltat, window=num.hanning):
'''
Time-domain normalization, in place.
Each component is divided by its running mean absolute amplitude,
i.e. its sample-wise absolute value, smoothed over a time window of
length ``deltat`` (see :py:meth:`smooth`, masked samples are left
out). This is the running absolute mean normalization commonly used
in ambient noise processing. Samples where this amplitude is zero are
set to zero.
:param deltat:
Length of the smoothing window [s].
:type deltat:
float
:param window:
Function returning the smoothing window for a given number of
samples.
:type window:
callable
'''
amp = self._smoothed(
num.abs(self._get_filled_float_data()), deltat, window)
data = num.zeros_like(amp)
num.divide(self.data, amp, out=data, where=amp > 0.0)
self.set_data(data)
[docs]
def apply_via_fft(self, f, ntrans_min=0):
'''
Apply a function to the spectra of the components, in place.
The spectra are computed with :py:meth:`get_spectrum`, passed to
``f``, and the result is transformed back to the time domain and cut
to the original number of samples.
:param f:
Function called as ``f(frequency_delta, ntrans, spectrum)``,
returning the modified spectrum (see :py:meth:`get_spectrum` for
the arguments).
:type f:
callable
:param ntrans_min:
Minimum FFT length.
:type ntrans_min:
int
'''
frequency_delta, ntrans, spectrum = self.get_spectrum(ntrans_min)
spectrum = f(frequency_delta, ntrans, spectrum)
data_new = num.fft.irfft(spectrum)[:, :self.nsamples]
self.set_data(data_new)
[docs]
def get_spectrum(self, ntrans_min=0):
'''
Get the spectra of the components.
Masked samples are replaced by zeros. The data is zero-padded to the
next power of two of at least ``ntrans_min`` and the number of
samples, and transformed with :py:func:`numpy.fft.rfft`.
:param ntrans_min:
Minimum FFT length.
:type ntrans_min:
int
:returns:
``(frequency_delta, ntrans, spectrum)``: frequency spacing [Hz],
FFT length, and complex spectra indexed as ``(icomponent,
ifrequency)``.
:rtype:
:py:class:`tuple`
'''
ntrans = trace.nextpow2(max(ntrans_min, self.nsamples))
data = self._get_filled_float_data()
spectrum = num.fft.rfft(data, ntrans)
frequency_delta = 1.0 / (self.deltat * ntrans)
return frequency_delta, ntrans, spectrum
[docs]
def get_cross_spectrum(self, ntrans_min=0):
'''
Get the cross spectra of all pairs of components.
:param ntrans_min:
Minimum FFT length, see :py:meth:`get_spectrum`.
:type ntrans_min:
int
:returns:
``(frequency_delta, ntrans, cross_spectrum)``, where
``cross_spectrum[i, j, k]`` is ``S_i(f_k) * conj(S_j(f_k))`` for
the spectra ``S`` of components ``i`` and ``j``.
:rtype:
:py:class:`tuple`
'''
frequency_delta, ntrans, spectrum = self.get_spectrum(ntrans_min)
return (
frequency_delta,
ntrans,
num.einsum('ik,jk->ijk', spectrum, num.conj(spectrum)))
[docs]
def get_component_codes_grouped(self, codes_projection):
'''
Group the components by their codes.
:param codes_projection:
Determines the grouping, see
:py:meth:`~pyrocko.squirrel.operators.base.CodesProjection.group_key`.
:type codes_projection:
:py:class:`~pyrocko.squirrel.operators.base.CodesProjection`
:returns:
Component indices by group key.
:rtype:
:py:class:`dict` of :py:class:`list` of :py:class:`int`
'''
key = codes_projection.group_key
groups = defaultdict(list)
for irow, component_codes in enumerate(self.component_codes):
groups[key(component_codes)].append(irow)
return groups
[docs]
def get_energy(
self,
codes_projection=None,
postprocessing=None):
'''
Get the energy, i.e. the sum of squared samples, of groups of
components.
The components are grouped with ``codes_projection``, e.g. the three
components of a sensor. Masked input samples do not contribute.
Samples with zero energy, e.g. where all input samples are masked,
are masked (and set to 1).
:param codes_projection:
Grouping of the components and codes of the resulting components.
Default: one group per sensor, i.e. by network, station,
location and channel code without the component character. The
resulting components get these codes and the extra code ``P``.
:type codes_projection:
:py:class:`~pyrocko.squirrel.operators.base.CodesProjection`
:param postprocessing:
Function applied to the result, see :py:meth:`apply`.
:type postprocessing:
callable
:returns:
New carpet, one component per group.
:rtype:
:py:class:`Carpet`
'''
from pyrocko.squirrel.operators.base import CodesProjection, Outlet
if codes_projection is None:
codes_projection = CodesProjection(
template='{i.network}.{i.station}.{i.location}'
'.{i.channel_no_component}.P')
groups = self.get_component_codes_grouped(codes_projection)
data = self._get_filled_float_data()
data **= 2
data3 = num.empty((len(groups), self.nsamples))
component_codes = []
for irow_out, irows_in in enumerate(groups.values()):
codes_group_in = [self.component_codes[i] for i in irows_in]
data3[irow_out, :] = data[irows_in, :].sum(axis=0)
codes_out, = codes_projection.project(
None, codes_group_in, [Outlet()])
component_codes.append(codes_out)
mask3 = data3 == 0.0
data3[mask3] = 1.0
energy = Carpet(
data=data3,
mask=mask3,
codes=self.codes,
component_codes=component_codes,
tmin=self.tmin,
deltat=self.deltat)
if postprocessing is not None:
energy.apply(postprocessing)
return energy
[docs]
def get_rms(self, codes_projection=None):
'''
Get the RMS amplitude, i.e. the square root of the energy, of groups
of components.
Like :py:meth:`get_energy`, but with the square root applied.
:param codes_projection:
Grouping of the components, see :py:meth:`get_energy`.
:type codes_projection:
:py:class:`~pyrocko.squirrel.operators.base.CodesProjection`
:returns:
New carpet, one component per group.
:rtype:
:py:class:`Carpet`
'''
return self.get_energy(
codes_projection,
postprocessing=lambda data: num.sqrt(data, out=data))
[docs]
def get_log_rms(self, codes_projection=None):
'''
Get the natural logarithm of the RMS amplitude of groups of
components.
Like :py:meth:`get_rms`, but logarithmic. Before taking the
logarithm, the energy is smoothed over three samples with weights
``(0.25, 0.5, 0.25)``. Masked samples and samples beyond the ends are
left out of the average.
:param codes_projection:
Grouping of the components, see :py:meth:`get_energy`.
:type codes_projection:
:py:class:`~pyrocko.squirrel.operators.base.CodesProjection`
:returns:
New carpet, one component per group. Masked samples are set to
zero.
:rtype:
:py:class:`Carpet`
'''
return self._get_log_rms(codes_projection, 0.5)
[docs]
def get_log10_rms(self, codes_projection=None):
'''
Get the decadic logarithm of the RMS amplitude of groups of
components.
Like :py:meth:`get_log_rms`, but with the decadic logarithm.
:param codes_projection:
Grouping of the components, see :py:meth:`get_energy`.
:type codes_projection:
:py:class:`~pyrocko.squirrel.operators.base.CodesProjection`
:returns:
New carpet, one component per group. Masked samples are set to
zero.
:rtype:
:py:class:`Carpet`
'''
return self._get_log_rms(codes_projection, 0.5 / num.log(10.0))
def _get_log_rms(self, codes_projection, factor):
energy = self.get_energy(codes_projection)
data = _smooth3(energy.data, energy.mask)
if energy.mask is not None:
data[energy.mask] = 1.0
num.log(data, out=data)
data *= factor
energy.set_data(data)
return energy
[docs]
def crop(self, fslice=slice(None, None)):
'''
Select a range of components.
:param fslice:
Components to keep. The component axes are cut accordingly.
:type fslice:
:py:class:`slice`
:returns:
New carpet, sharing the data with the original.
:rtype:
:py:class:`Carpet`
'''
return Carpet(
tmin=self.tmin,
deltat=self.deltat,
codes=self.codes,
component_codes=self.component_codes[fslice],
component_axes=dict(
(k, v[fslice]) for (k, v) in self.component_axes.items()),
data=self.data[fslice, :],
mask=None if self.mask is None else self.mask[fslice, :])
[docs]
def resample_band(
self, fmin, fmax, nf,
registration='cell',
scale='log',
component_axis='frequency'):
'''
Resample the components to regularly spaced bins of a component axis.
Typically used to bring a spectrogram to ``nf`` frequency bands,
regularly spaced on a logarithmic or linear frequency scale. The
value of each new band is the average of the data over the band,
computed by integrating (trapezoidal rule, on the chosen scale) along
the component axis, which must be increasing.
Components whose axis value does not meet the scale's condition, e.g.
frequency 0 for logarithmic scale, are dropped first. Masked and NaN
input samples count as zero in the integration. Bands outside the
range of the component axis, or with a result of exactly zero, are
set to NaN. The result has no mask.
:param fmin:
Lower end of the range of the new bands. Default: lower end of
the component axis.
:type fmin:
float
:param fmax:
Upper end of the range of the new bands. Default: upper end of
the component axis.
:type fmax:
float
:param nf:
Number of new bands.
:type nf:
int
:param registration:
``'cell'``: ``fmin`` and ``fmax`` are the outer edges of the first
and last band. ``'node'``: they are (approximately) the centres of
the first and last band.
:type registration:
str
:param scale:
``'log'``, ``'lin'``, or a tuple ``(f, f_inverse, condition)`` of
functions for a custom scale, e.g. ``(numpy.log, numpy.exp,
lambda y: y > 0)``.
:type scale:
str or :py:class:`tuple`
:param component_axis:
Name of the component axis in :py:gattr:`component_axes`. If it
is not present, the component indices are used.
:type component_axis:
str
:raises:
:py:exc:`CarpetResampleError` if the scale is invalid, the
component axis is not increasing, or no data is in the requested
range.
:returns:
New carpet with ``nf`` components, empty component codes, and the
band centres as component axis ``component_axis``.
:rtype:
:py:class:`Carpet`
'''
if scale == 'log':
log, exp, condition = num.log, num.exp, lambda y: y > 0
elif scale == 'lin':
log, exp, condition = (
lambda y: y,
lambda y: y,
lambda y: num.ones_like(y, dtype=bool))
elif isinstance(scale, tuple):
log, exp, condition = scale
else:
raise CarpetResampleError(
'Carpet.resample_band: Scale must be "lin", "log" or a tuple '
'with scaling, inverse scaling, and condition functions, e.g. '
'`(log, exp, lambda y: y > 0)`.')
def get_frequencies(cp):
frequencies = cp.component_axes.get(component_axis)
if frequencies is None:
frequencies = num.arange(cp.ncomponents, dtype=float)
return frequencies
frequencies = get_frequencies(self)
if not num.all(num.diff(frequencies) > 0.0):
raise CarpetResampleError(
'Carpet.resample_band: Component axis must be monotonically '
'increasing.')
iok = num.where(condition(frequencies))[0]
if iok.size < self.ncomponents:
if iok.size == 0:
raise CarpetResampleError(
'Carpet.resample_band: No elements of component axis meet '
'scaling condition (e.g. y > 0 for log scaling)')
self = self.crop(fslice=slice(iok[0], None))
frequencies = get_frequencies(self)
if fmin is None:
fmin = frequencies[0]
if fmax is None:
fmax = frequencies[-1]
if registration == 'cell':
log_frequencies_out = num.linspace(
log(fmin), log(fmax), nf+1)
elif registration == 'node':
d_log_f = log(fmin + (fmax - fmin) / nf) - log(fmin)
log_frequencies_out = num.linspace(
log(fmin)-0.5*d_log_f,
log(fmax)+0.5*d_log_f,
nf+1)
log_frequencies = log(frequencies)
iok = num.where(num.logical_and(
log_frequencies[0] <= log_frequencies_out,
log_frequencies_out <= log_frequencies[-1]))[0]
if iok.size == 0:
raise CarpetResampleError(
'No frequency values in requested range.')
fslice = slice(iok[0], iok[-1]+1)
fslice_out = slice(iok[0], iok[-1])
d_log_frequencies_half = num.diff(log_frequencies) * 0.5
values = self.data.astype(num.float64)
if self.mask is not None:
values[self.mask] = num.nan
int_values = num.zeros_like(values)
num.nancumsum(
(values[1:, :] + values[:-1, :])
* d_log_frequencies_half[:, num.newaxis],
axis=0,
out=int_values[1:, :])
i = num.searchsorted(log_frequencies, log_frequencies_out[fslice])
w1 = (log_frequencies_out[fslice] - log_frequencies[i-1]) \
/ (log_frequencies[i] - log_frequencies[i-1])
w0 = 1.0 - w1
int_values_out = int_values[i-1, :] * w0[:, num.newaxis] \
+ int_values[i, :] * w1[:, num.newaxis]
values_out = num.full((nf, self.data.shape[1]), num.nan)
values_out[fslice_out, :] = num.diff(int_values_out, axis=0) \
/ num.diff(log_frequencies_out[fslice])[:, num.newaxis]
values_out[values_out == 0.0] = num.nan
frequencies_out = exp(
0.5 * (log_frequencies_out[1:] + log_frequencies_out[:-1]))
return Carpet(
codes=self.codes,
tmin=self.tmin,
deltat=self.deltat,
component_axes={component_axis: frequencies_out},
data=values_out)
[docs]
def correlate(a, b, mode='valid', normalization=None, use_fft=False):
'''
Cross-correlate traces and/or carpets.
Each of ``a`` and ``b`` can be a :py:class:`~pyrocko.trace.Trace` or a
:py:class:`Carpet`. Components are correlated pairwise with
:py:func:`pyrocko.trace.correlate`, see there for the other arguments.
:returns:
For two traces, the correlation trace. Otherwise a carpet with one
component per pair of components (for two carpets, ordered by the
components of ``a``, then those of ``b``).
:rtype:
:py:class:`~pyrocko.trace.Trace` or :py:class:`Carpet`
'''
if isinstance(a, Trace) and isinstance(b, Trace):
return trace.correlate(
a, b, mode=mode, normalization=normalization, use_fft=use_fft)
elif isinstance(a, Trace) and isinstance(b, Carpet):
return Carpet([
trace.correlate(
a, b_,
mode=mode, normalization=normalization, use_fft=use_fft)
for b_ in b])
elif isinstance(a, Carpet) and isinstance(b, Trace):
return Carpet([
trace.correlate(
a_, b,
mode=mode, normalization=normalization, use_fft=use_fft)
for a_ in a])
elif isinstance(a, Carpet) and isinstance(b, Carpet):
return Carpet([
trace.correlate(
a_, b_,
mode=mode, normalization=normalization, use_fft=use_fft)
for a_ in a for b_ in b])
[docs]
def join(carpets):
'''
Join carpets which are contiguous in time.
Carpets with the same codes, number of components, sampling interval and
data type are concatenated if one starts right after the other ends.
Carpets without samples are dropped. Component codes and component axes
are taken from the first carpet of each joined group. Masks are joined
as well.
:param carpets:
Carpets to join.
:type carpets:
:py:class:`list` of :py:class:`Carpet`
:returns:
Joined carpets, sorted by codes and start time.
:rtype:
:py:class:`list` of :py:class:`Carpet`
'''
if not carpets:
return []
eps = 1e-4
carpets = sorted(carpets, key=lambda cp: (cp.codes, cp.tmin))
groups = []
for current in carpets:
if current.nsamples == 0:
continue
lastgroup = groups[-1] if groups else None
last = lastgroup[-1] if lastgroup else None
deltat = current.deltat
if (last
and last.deltat == deltat
and abs(last.tmax + deltat - current.tmin) < eps * deltat
and last.ncomponents == current.ncomponents
and all(
a == b
for (a, b)
in zip(last.codes, current.codes))
and last.data.dtype == current.data.dtype):
lastgroup.append(current)
else:
groups.append([current])
carpets_out = []
for group in groups:
if len(group) > 1:
data = num.hstack([cp.data for cp in group])
if all(cp.mask is None for cp in group):
mask = None
else:
mask = num.hstack([cp.get_mask() for cp in group])
carpet = group[0].copy(data=data)
carpet.set_mask(mask)
carpets_out.append(carpet)
else:
carpets_out.append(group[0])
return carpets_out
[docs]
def check_same_components(a, b):
'''
Check that two carpets have the same components.
Compares the number of components, the component codes and the component
axes.
:raises:
:py:exc:`CarpetError` if the components differ.
'''
if a.ncomponents != b.ncomponents:
raise CarpetError('Different number of components.')
for i, (a_codes, b_codes) in enumerate(
zip(a.component_codes, b.component_codes)):
if a_codes != b_codes:
raise CarpetError(
'Component codes differ at component %i: %s != %s.' % (
i, a_codes, b_codes))
ks = set(a.component_axes.keys())
ks.update(b.component_axes.keys())
for k in ks:
if k not in a.component_axes or k not in b.component_axes:
raise CarpetError('Component axis not present: %s' % k)
av = a.component_axes[k]
bv = b.component_axes[k]
if not num.all(av == bv):
raise CarpetError('Component axis values differ.')
[docs]
def check_compatible(carpets):
'''
Check that carpets have the same sampling and components.
The carpets must have the same sampling interval, start time and number
of samples, and the same components (see :py:func:`check_same_components`).
:raises:
:py:exc:`CarpetError` if the carpets are not compatible.
'''
vecs = [
(carpet.deltat, carpet.tmin, carpet.nsamples) for carpet in carpets]
if not all(vec == vecs[0] for vec in vecs[1:]):
raise CarpetError(
'Carpets differ in sampling interval, onset time'
' or number of samples.')
for carpet in carpets[1:]:
check_same_components(carpets[0], carpet)
[docs]
def sum(carpets, nan_aware=False, normalize=False, dtype=float):
'''
Sum carpets sample by sample.
The carpets must be compatible (see :py:func:`check_compatible`).
:param carpets:
Carpets to sum.
:type carpets:
:py:class:`list` of :py:class:`Carpet`
:param nan_aware:
If ``True``, non-finite values (NaN, inf) are skipped. Otherwise they
propagate into the sum. Masked samples are always skipped.
:type nan_aware:
bool
:param normalize:
If ``True``, divide by the number of finite, unmasked values summed,
i.e. compute the average.
:type normalize:
bool
:param dtype:
Data type of the sum.
:type dtype:
:py:class:`numpy.dtype`
:raises:
:py:exc:`CarpetError` if no carpets are given, a carpet has no data,
or the carpets are not compatible.
:returns:
New carpet, with codes ``.SUM...`` and the components of the first
carpet. Samples with no finite, unmasked values are NaN. The result
has no mask.
:rtype:
:py:class:`Carpet`
'''
if not carpets:
raise CarpetError('Need at least one carpet for summation.')
for carpet in carpets:
if carpet.data is None:
raise CarpetError(
'Carpets supplied for summation must contain data.')
try:
check_compatible(carpets)
except CarpetError as e:
raise CarpetError(
'Incompatible carpets supplied for summation (%s):\n%s' % (
str(e),
'\n'.join(' %s' % carpet.summary for carpet in carpets)))
data = num.zeros(carpets[0].data.shape, dtype=dtype)
count = num.zeros(carpets[0].data.shape, dtype=int)
for carpet in carpets:
ok = num.isfinite(carpet.data)
if carpet.mask is not None:
ok &= ~carpet.mask
if nan_aware:
data += num.where(ok, carpet.data, 0.0)
elif carpet.mask is not None:
data += num.where(carpet.mask, 0.0, carpet.data)
else:
data += carpet.data
count += ok
if normalize:
with num.errstate(divide='ignore', invalid='ignore'):
data /= count
data[count == 0] = num.nan
carpet_sum = carpets[0].copy(data='drop')
carpet_sum.codes = CodesNSLCE('.SUM...')
carpet_sum.data = data
return carpet_sum
[docs]
def check_overlaps(carpets_a, carpets_b=None, message='Carpets overlap.'):
'''
Check that carpets with the same codes do not overlap in time.
:param carpets_a:
Carpets to check.
:type carpets_a:
:py:class:`list` of :py:class:`Carpet`
:param carpets_b:
If given, check each carpet of ``carpets_a`` against these instead of
against the others in ``carpets_a``.
:type carpets_b:
:py:class:`list` of :py:class:`Carpet`
:param message:
Beginning of the error message.
:type message:
str
:raises:
:py:exc:`CarpetOverlapError` on the first overlap found.
'''
for ia, cp_a in enumerate(carpets_a):
for cp_b in carpets_a[ia+1:] if carpets_b is None else carpets_b:
if cp_a.codes == cp_b.codes and cp_a.overlaps(cp_b):
raise CarpetOverlapError(
message + '\n Carpet 1: %s\n Carpet 2: %s' % (
cp_a.summary, cp_b.summary))
#: Remove overlaps between carpets with the same codes, see
#: :py:func:`pyrocko.trace.deoverlap`, which works for carpets as well.
deoverlap = trace.deoverlap
