"""
Functions for network matched-filter detection of seismic data.
Designed to cross-correlate templates generated by template_gen function
with data and output the detections.
:copyright:
EQcorrscan developers.
:license:
GNU Lesser General Public License, Version 3
(https://www.gnu.org/copyleft/lesser.html)
"""
import copy
import glob
import os
import shutil
import pickle
import tarfile
import tempfile
import logging
from os.path import join
import numpy as np
from obspy import Catalog, read_events, Stream
from eqcorrscan.core.match_filter.family import _write_family, _read_family
from eqcorrscan.core.match_filter.matched_filter import MatchFilterError
from eqcorrscan.core.match_filter.template import Template, group_templates
from eqcorrscan.core.match_filter.family import Family
from eqcorrscan.core.match_filter.detection import write_detections
from eqcorrscan.core.match_filter.helpers import (
_total_microsec, temporary_directory, _safemembers, _templates_match)
from eqcorrscan.utils.catalog_utils import _get_origin
from eqcorrscan.utils.findpeaks import decluster, decluster_distance_time
from eqcorrscan.utils.plotting import cumulative_detections
Logger = logging.getLogger(__name__)
[docs]
class Party(object):
"""
Container for multiple Family objects.
:type families: list of Family
:type families: The family objects that make up the Party.
"""
[docs]
def __init__(self, families=None):
"""Instantiate the Party object."""
self.families = []
if isinstance(families, Family):
families = [families]
if families:
self.families.extend(families)
def __repr__(self):
"""
Print short info about the Party.
:return: str
.. rubric:: Example
>>> print(Party())
Party of 0 Families.
"""
print_str = ('Party of %s Families.' % len(self.families))
return print_str
def __iadd__(self, other):
"""
Method for in-place addition '+='. Uses the Party.__add__() method.
:type other: `Party`
:param other: Another party to merge with the current family.
:return: Works in place on self.
.. rubric:: Example
>>> from eqcorrscan import Family, Template
>>> party_a = Party(families=[Family(template=Template(name='a'),
... detections=[])])
>>> party_b = Party(families=[Family(template=Template(name='b'),
... detections=[])])
>>> party_a += party_b
>>> print(party_a)
Party of 2 Families.
"""
if isinstance(other, Family):
families = [other]
elif isinstance(other, Party):
families = other.families
else:
raise NotImplementedError(
'Ambiguous add, only allowed Party or Family additions.')
for oth_fam in families:
fam = self.select(oth_fam.template.name)
# This check is taken care of by Family.__iadd__
# assert fam.template == oth_fam.template, (
# "Matching template names, but different templates")
if fam is not None:
fam += oth_fam
else:
self.families.append(oth_fam)
return self
def __add__(self, other):
"""
Method for addition '+'.
:type other: `Party` or `Family`
:param other: Another party to merge with the current family.
:return: Works in place on self.
.. note:: Works in place on party, will alter this original party.
.. rubric:: Example
Addition of two parties together:
>>> from eqcorrscan import Family, Template
>>> party_a = Party(families=[Family(template=Template(name='a'),
... detections=[])])
>>> party_b = Party(families=[Family(template=Template(name='b'),
... detections=[])])
>>> party_c = party_a + party_b
>>> print(party_c)
Party of 2 Families.
Addition of a family to a party:
>>> party_a = Party(families=[Family(template=Template(name='a'),
... detections=[])])
>>> family_b = Family(template=Template(name='b'), detections=[])
>>> party_c = party_a + family_b
>>> print(party_c)
Party of 2 Families.
Addition of a party with some families using the same templates:
>>> party_a = Party(families=[Family(template=Template(name='a'),
... detections=[])])
>>> party_b = Party(families=[Family(template=Template(name='a'),
... detections=[])])
>>> party_c = party_a + party_b
>>> print(party_c)
Party of 1 Families.
Addition of non Family or Party objects is not allows:
>>> party_a = Party(families=[Family(template=Template(name='a'))])
>>> misc = 1.0
>>> party_c = party_a + misc # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
NotImplementedError: Ambiguous add, only allowed Party or Family \
additions
"""
return self.copy().__iadd__(other)
def __eq__(self, other, verbose=False):
"""
Equality testing, rich comparison '=='.
:param other: Another object.
:return: bool
.. rubric:: Example
Compare equal Parties:
>>> from eqcorrscan import Family, Template
>>> party_a = Party(families=[Family(template=Template(name='a'))])
>>> party_b = Party(families=[Family(template=Template(name='a'))])
>>> party_a == party_b
True
Compare Parties with different templates:
>>> party_a = Party(families=[Family(template=Template(name='a'))])
>>> party_b = Party(families=[Family(template=Template(name='b'))])
>>> party_a == party_b
False
Compare a Party with a Family:
>>> party = Party(families=[Family(template=Template(name='a'))])
>>> family = Family(template=Template(name='a'))
>>> party == family
False
"""
if not isinstance(other, Party):
return False
for family, oth_fam in zip(self.sort().families,
other.sort().families):
if not family.__eq__(oth_fam, verbose=verbose):
return False
else:
return True
def __ne__(self, other):
"""
Rich comparison operator '!='.
:param other: other object
:return: bool
.. rubric:: Example
Compare two equal Parties:
>>> from eqcorrscan import Family, Template
>>> party_a = Party(families=[Family(template=Template(name='a'))])
>>> party_b = Party(families=[Family(template=Template(name='a'))])
>>> party_a != party_b
False
"""
return not self.__eq__(other)
def __getitem__(self, index):
"""
Get families from the Party. Can accept either an index or slice.
:param index: Family number or name, or slice.
:return: Party (if a slice is given) or a single Family
.. rubric:: Example
Extract a single family:
>>> from eqcorrscan import Family, Template
>>> party = Party(families=[Family(template=Template(name='a')),
... Family(template=Template(name='b')),
... Family(template=Template(name='c'))])
>>> party[1]
Family of 0 detections from template b
Extract a list of families by giving a slice:
>>> party = Party(families=[Family(template=Template(name='a')),
... Family(template=Template(name='b')),
... Family(template=Template(name='c'))])
>>> party[1:]
Party of 2 Families.
Extract a single family by template name:
>>> party = Party(families=[Family(template=Template(name='a')),
... Family(template=Template(name='b')),
... Family(template=Template(name='c'))])
>>> party['b']
Family of 0 detections from template b
"""
if len(self.families) == 0:
return None
if isinstance(index, slice):
return self.__class__(families=self.families.__getitem__(index))
elif isinstance(index, int):
return self.families.__getitem__(index)
else:
_index = [i for i, family in enumerate(self.families)
if family.template.name == index]
try:
return self.families.__getitem__(_index[0])
except IndexError:
Logger.warning('Family: %s not in party' % index)
return []
def __iter__(self):
return self.families.__iter__()
def __next__(self):
return self.families.__next__()
def __len__(self):
"""
Get total number of detections in Party.
:return: length, int
.. rubric:: Example
>>> from eqcorrscan import Family, Template
>>> party = Party(families=[Family(template=Template(name='a')),
... Family(template=Template(name='b'))])
>>> len(party)
0
"""
length = 0
for family in self.families:
length += len(family)
return length
@property
def _template_dict(self):
return {family.template.name: family for family in self}
[docs]
def select(self, template_name):
"""
Select a specific family from the party.
:type template_name: str
:param template_name: Template name of Family to select from a party.
:returns: Family
"""
return self._template_dict.get(template_name)
[docs]
def sort(self):
"""
Sort the families by template name.
.. rubric:: Example
>>> from eqcorrscan import Family, Template
>>> party = Party(families=[Family(template=Template(name='b')),
... Family(template=Template(name='a'))])
>>> party[0]
Family of 0 detections from template b
>>> party.sort()[0]
Family of 0 detections from template a
"""
self.families.sort(key=lambda x: x.template.name)
return self
[docs]
def filter(self, dates=None, min_dets=1):
"""
Return a new Party filtered according to conditions.
Return a new Party with only detections within a date range and
only families with a minimum number of detections.
:type dates: list of obspy.core.UTCDateTime objects
:param dates: A start and end date for the new Party
:type min_dets: int
:param min_dets: Minimum number of detections per family
.. rubric:: Example
>>> from obspy import UTCDateTime
>>> Party().read().filter(dates=[UTCDateTime(2016, 1, 1),
... UTCDateTime(2017, 1, 1)],
... min_dets=30) # doctest: +SKIP
"""
if dates is None:
raise MatchFilterError('Need a list defining a date range')
new_party = Party()
for fam in self.families:
new_fam = Family(
template=fam.template,
detections=[det for det in fam if
dates[0] < det.detect_time < dates[1]])
if len(new_fam) >= min_dets:
new_party.families.append(new_fam)
return new_party
[docs]
def plot(self, plot_grouped=False, dates=None, min_dets=1, rate=False,
**kwargs):
"""
Plot the cumulative detections in time.
:type plot_grouped: bool
:param plot_grouped:
Whether to plot all families together (plot_grouped=True), or each
as a separate line.
:type dates: list
:param dates: list of obspy.core.UTCDateTime objects bounding the
plot. The first should be the start date, the last the end date.
:type min_dets: int
:param min_dets: Plot only families with this number of detections
or more.
:type rate: bool
:param rate: Whether or not to plot the daily rate of detection as
opposed to cumulative number. Only works with plot_grouped=True.
:param \**kwargs: Any other arguments accepted by
:func:`eqcorrscan.utils.plotting.cumulative_detections`
.. rubric:: Examples
Plot cumulative detections for all templates individually:
>>> Party().read().plot() # doctest: +SKIP
Plot cumulative detections for all templates grouped together:
>>> Party().read().plot(plot_grouped=True) # doctest: +SKIP
Plot the rate of detection for all templates grouped together:
>>> Party().read().plot(plot_grouped=True, rate=True) # doctest: +SKIP
Plot cumulative detections for all templates with more than five
detections between June 1st, 2012 and July 31st, 2012:
>>> from obspy import UTCDateTime
>>> Party().read().plot(dates=[UTCDateTime(2012, 6, 1),
... UTCDateTime(2012, 7, 31)],
... min_dets=5) # doctest: +SKIP
"""
all_dets = []
if dates:
new_party = self.filter(dates=dates, min_dets=min_dets)
for fam in new_party.families:
all_dets.extend(fam.detections)
else:
for fam in self.families:
all_dets.extend(fam.detections)
fig = cumulative_detections(detections=all_dets,
plot_grouped=plot_grouped,
rate=rate, **kwargs)
return fig
[docs]
def rethreshold(self, new_threshold, new_threshold_type='MAD',
abs_values=False):
"""
Remove detections from the Party that are below a new threshold.
.. Note:: threshold can only be set higher.
.. Warning::
Works in place on Party.
:type new_threshold: float
:param new_threshold: New threshold level
:type new_threshold_type: str
:param new_threshold_type: Either 'MAD', 'absolute' or 'av_chan_corr'
:type abs_values: bool
:param abs_values:
Whether to compare the absolute value of the detection-value.
.. rubric:: Examples
Using the MAD threshold on detections made using the MAD threshold:
>>> party = Party().read()
>>> len(party)
4
>>> party = party.rethreshold(10.0)
>>> len(party)
4
>>> # Note that all detections are self detections
Using the absolute thresholding method on the same Party:
>>> party = Party().read().rethreshold(5.9, 'absolute')
>>> len(party)
1
Using the av_chan_corr method on the same Party:
>>> party = Party().read().rethreshold(0.9, 'av_chan_corr')
>>> len(party)
4
"""
for family in self.families:
rethresh_detections = []
for d in family.detections:
if new_threshold_type == 'MAD' and d.threshold_type == 'MAD':
new_thresh = (d.threshold /
d.threshold_input) * new_threshold
elif new_threshold_type == 'MAD' and d.threshold_type != 'MAD':
raise MatchFilterError(
'Cannot recalculate MAD level, '
'use another threshold type')
elif new_threshold_type == 'absolute':
new_thresh = new_threshold
elif new_threshold_type == 'av_chan_corr':
new_thresh = new_threshold * d.no_chans
else:
raise MatchFilterError(
'new_threshold_type %s is not recognised' %
str(new_threshold_type))
rethresh = False
if abs_values:
if abs(float(d.detect_val)) >= new_thresh:
rethresh = True
else:
if float(d.detect_val) >= new_thresh:
rethresh = True
if rethresh:
d.threshold = new_thresh
d.threshold_input = new_threshold
d.threshold_type = new_threshold_type
rethresh_detections.append(d)
family.detections = rethresh_detections
return self
[docs]
def decluster(self, trig_int, timing='detect', metric='avg_cor',
hypocentral_separation=None, min_chans=0,
absolute_values=False, num_threads=None):
"""
De-cluster a Party of detections by enforcing a detection separation.
De-clustering occurs between events detected by different (or the same)
templates. If multiple detections occur within trig_int then the
preferred detection will be determined by the metric argument. This
can be either the average single-station correlation coefficient which
is calculated as Detection.detect_val / Detection.no_chans, or the
raw cross channel correlation sum which is simply Detection.detect_val.
:type trig_int: float
:param trig_int: Minimum detection separation in seconds.
:type metric: str
:param metric: What metric to sort peaks by. Either 'avg_cor' which
takes the single station average correlation, 'cor_sum' which
takes the total correlation sum across all channels, or
'thresh_exc' which takes the factor by how much the detection
exceeded the input threshold.
:type timing: str
:param timing:
Either 'detect' or 'origin' to decluster based on either the
detection time or the origin time.
:type hypocentral_separation: float
:param hypocentral_separation:
Maximum inter-event separation in km to decluster events within.
If an event happens within this distance of another event, and
within the trig_int defined time, then they will be declustered.
:type min_chans: int
:param min_chans:
Minimum number of channels for a detection to be retained. If a
detection was made on very few channels, then the MAD detection
statistics become much less robust.
:type absolute_values: bool
:param absolute_values:
Use the absolute value of the metric to choose the preferred
detection.
:type num_threads: int
:param num_threads:
Number of threads to use for internal c-funcs if available.
Only valid if hypocentral_separation used.
.. Warning::
Works in place on object, if you need to keep the original safe
then run this on a copy of the object!
.. rubric:: Example
>>> party = Party().read()
>>> len(party)
4
>>> declustered = party.decluster(20)
>>> len(party)
3
.. rubric:: Example using epicentral-distance
>>> party = Party().read()
>>> len(party)
4
>>> # Calculate the expected origins based on the template origin.
>>> for family in party:
... for detection in family:
... _ = detection._calculate_event(template=family.template)
>>> declustered = party.decluster(20, hypocentral_separation=1)
>>> len(party)
4
>>> declustered = party.decluster(20, hypocentral_separation=100)
>>> len(party)
3
"""
if self.__len__() == 0:
return self
all_detections = [d for fam in self.families for d in fam.detections]
catalog = None
if hypocentral_separation:
catalog = Catalog([d.event for fam in self.families
for d in fam.detections if d.event])
if len(catalog) == 0:
Logger.warning("Detections do not have events, cannot use "
"hypocentral separation")
catalog = None
if min_chans > 0:
for d in all_detections.copy():
if d.no_chans < min_chans:
all_detections.remove(d)
if len(all_detections) == 0:
return Party()
assert metric in ('avg_cor', 'cor_sum', 'thresh_exc'), \
'metric is not cor_sum, avg_cor or thresh_exc'
assert timing in ('detect', 'origin'), 'timing is not detect or origin'
if timing == 'detect':
if metric == 'avg_cor':
detect_info = [(d.detect_time, d.detect_val / d.no_chans)
for d in all_detections]
elif metric == 'cor_sum':
detect_info = [(d.detect_time, d.detect_val)
for d in all_detections]
elif metric == 'thresh_exc':
detect_info = [(d.detect_time, d.detect_val / d.threshold)
for d in all_detections]
else:
if metric == 'avg_cor':
detect_info = [(_get_origin(d.event).time,
d.detect_val / d.no_chans)
for d in all_detections]
elif metric == 'cor_sum':
detect_info = [(_get_origin(d.event).time, d.detect_val)
for d in all_detections]
elif metric == 'thresh_exc':
detect_info = [(_get_origin(d.event).time,
d.detect_val / d.threshold)
for d in all_detections]
if absolute_values:
for j, tup in enumerate(detect_info):
detect_info[j] = (tup[0], abs(tup[1]))
min_det = sorted([d[0] for d in detect_info])[0]
detect_vals = np.array([d[1] for d in detect_info], dtype=np.float32)
detect_times = np.array([
_total_microsec(d[0].datetime, min_det.datetime)
for d in detect_info])
# Trig_int must be converted from seconds to micro-seconds
if hypocentral_separation and catalog:
peaks_out = decluster_distance_time(
peaks=detect_vals, index=detect_times,
trig_int=trig_int * 10 ** 6, catalog=catalog,
hypocentral_separation=hypocentral_separation,
num_threads=num_threads)
else:
peaks_out = decluster(
peaks=detect_vals, index=detect_times,
trig_int=trig_int * 10 ** 6)
# Need to match both the time and the detection value
declustered_detections = []
for ind in peaks_out:
matching_time_indices = np.where(detect_times == ind[-1])[0]
matches = matching_time_indices[
np.where(detect_vals[matching_time_indices] == ind[0])[0][0]]
declustered_detections.append(all_detections[matches])
# Convert this list into families
template_names = list(set([d.template_name
for d in declustered_detections]))
new_families = []
for template_name in template_names:
template = [fam.template for fam in self.families
if fam.template.name == template_name][0]
new_families.append(Family(
template=template,
detections=[d for d in declustered_detections
if d.template_name == template_name]))
# TODO: this might be better changing the list of families in place.
self.families = new_families
return self
[docs]
def copy(self):
"""
Returns a copy of the Party.
:return: Copy of party
.. rubric:: Example
>>> from eqcorrscan import Family, Template
>>> party = Party(families=[Family(template=Template(name='a'))])
>>> party_b = party.copy()
>>> party == party_b
True
"""
return copy.deepcopy(self)
[docs]
def write(self, filename, format='tar', write_detection_catalog=True,
catalog_format="QUAKEML", overwrite=False):
"""
Write Family out, select output format.
:type format: str
:param format:
One of either 'tar', 'csv', or any obspy supported
catalog output. See note below on formats
:type filename: str
:param filename: Path to write file to.
:type write_detection_catalog: bool
:param write_detection_catalog:
Whether to write the detection catalog object or not - writing
large catalog files can be slow, and catalogs can be reconstructed
from the Tribe.
:type catalog_format: str
:param catalog_format:
What format to write the detection-catalog with. Only Nordic,
SC3ML, QUAKEML are supported. Note that not all information is
written for all formats (QUAKEML is the most complete, but is
slow for IO).
:type overwrite: bool
:param overwrite:
Specifies whether detection-files are overwritten if they exist
already. By default, no files are overwritten.
.. NOTE::
csv format will write out detection objects, all other
outputs will write the catalog. These cannot be rebuilt into
a Family object. The only format that can be read back into
Family objects is the 'tar' type.
.. NOTE::
We recommend writing to the 'tar' format, which will write out
all the template information (wavefiles as miniseed and metadata)
alongside the detections and store these in a tar archive. This
is readable by other programs and maintains all information
required for further study.
.. rubric:: Example
>>> party = Party().read()
>>> party.write('test_tar_write', format='tar')
Party of 4 Families.
>>> party.write('test_csv_write.csv', format='csv')
Party of 4 Families.
>>> party.write('test_quakeml.xml', format='quakeml')
Party of 4 Families.
"""
from eqcorrscan.core.match_filter.tribe import Tribe
from eqcorrscan.core.match_filter import CAT_EXT_MAP
if catalog_format not in CAT_EXT_MAP.keys():
raise TypeError("{0} is not supported".format(catalog_format))
if format.lower() == 'csv':
if os.path.isfile(filename) and not overwrite:
raise MatchFilterError(
'Will not overwrite existing file: %s' % filename)
if os.path.isfile(filename) and overwrite:
os.remove(filename)
for family in self.families:
write_detections(fname=filename, detections=family.detections,
mode="a")
elif format.lower() == 'tar':
if not filename.endswith('.tgz'):
filename = filename + ".tgz"
if os.path.exists(filename) and not overwrite:
raise IOError('Will not overwrite existing file: %s'
% filename)
# os.makedirs(filename)
with temporary_directory() as temp_dir:
Tribe([f.template for f in self.families]).write(
filename=temp_dir, compress=False,
catalog_format=catalog_format)
if write_detection_catalog:
all_cat = Catalog()
for family in self.families:
all_cat += family.catalog
if not len(all_cat) == 0:
all_cat.write(
join(temp_dir, 'catalog.{0}'.format(
CAT_EXT_MAP[catalog_format])),
format=catalog_format)
for i, family in enumerate(self.families):
Logger.debug('Writing family %i' % i)
name = family.template.name + '_detections.csv'
name_to_write = join(temp_dir, name)
_write_family(family=family, filename=name_to_write)
with tarfile.open(filename, "w:gz") as tar:
tar.add(temp_dir, arcname=os.path.basename(filename))
else:
Logger.warning('Writing only the catalog component, metadata '
'will not be preserved')
self.get_catalog().write(filename=filename, format=format)
return self
[docs]
def read(self, filename=None, read_detection_catalog=True,
estimate_origin=True):
"""
Read a Party from a file.
:type filename: str
:param filename:
File to read from - can be a list of files, and can contain
wildcards.
:type read_detection_catalog: bool
:param read_detection_catalog:
Whether to read the detection catalog or not, if False, catalog
will be regenerated - for large catalogs this can be faster.
:type estimate_origin: bool
:param estimate_origin:
If True and no catalog is found, or read_detection_catalog is False
then new events with origins estimated from the template origin
time will be created.
.. rubric:: Example
>>> Party().read()
Party of 4 Families.
"""
from eqcorrscan.core.match_filter.tribe import Tribe
tribe = Tribe()
families = []
if filename is None:
# If there is no filename given, then read the example.
filename = os.path.join(
os.path.dirname(__file__), '..', '..', 'tests', 'test_data',
'test_party.tgz')
if isinstance(filename, list):
filenames = []
for _filename in filename:
# Expand wildcards
filenames.extend(glob.glob(_filename))
else:
# Expand wildcards
filenames = glob.glob(filename)
for _filename in filenames:
Logger.info(f"Reading from {_filename}")
# Cope with pickled files
if _filename.endswith('.pkl'):
with open(_filename, "rb") as _f:
chunk_party = pickle.load(_f)
self.__iadd__(chunk_party)
continue
with tarfile.open(_filename, "r:*") as arc:
temp_dir = tempfile.mkdtemp()
arc.extractall(path=temp_dir, members=_safemembers(arc))
# Read in the detections first, this way, if we read from multiple
# files then we can just read in extra templates as needed.
# Read in families here!
party_dir = glob.glob(temp_dir + os.sep + '*')[0]
tribe._read_from_folder(dirname=party_dir)
det_cat_file = glob.glob(os.path.join(party_dir, "catalog.*"))
if len(det_cat_file) != 0 and read_detection_catalog:
try:
all_cat = read_events(det_cat_file[0])
except TypeError as e:
Logger.error(e)
pass
else:
all_cat = Catalog()
for family_file in glob.glob(join(party_dir, '*_detections.csv')):
template = [
t for t in tribe if _templates_match(t, family_file)]
family = Family(template=template[0] or Template())
new_family = True
if family.template.name in [f.template.name for f in families]:
family = [
f for f in families if
f.template.name == family.template.name][0]
new_family = False
family.detections += _read_family(
fname=family_file, all_cat=all_cat, template=template[0],
estimate_origin=estimate_origin)
if new_family:
families.append(family)
shutil.rmtree(temp_dir)
self.families = families
return self
[docs]
def lag_calc(self, stream, pre_processed, shift_len=0.2, min_cc=0.4,
min_cc_from_mean_cc_factor=None, vertical_chans=['Z'],
horizontal_chans=['E', 'N', '1', '2'],
cores=1, interpolate=False, plot=False, plotdir=None,
parallel=True, process_cores=None, ignore_length=False,
ignore_bad_data=False, export_cc=False, cc_dir=None,
**kwargs):
"""
Compute picks based on cross-correlation alignment.
Works in-place on events in Party.
:type stream: obspy.core.stream.Stream
:param stream:
All the data needed to cut from - can be a gappy Stream.
:type pre_processed: bool
:param pre_processed:
Whether the stream has been pre-processed or not to match the
templates. See note below.
:type shift_len: float
:param shift_len:
Shift length allowed for the pick in seconds, will be plus/minus
this amount - default=0.2
:type min_cc: float
:param min_cc:
Minimum cross-correlation value to be considered a pick,
default=0.4.
:type min_cc_from_mean_cc_factor: float
:param min_cc_from_mean_cc_factor:
If set to a value other than None, then the minimum cross-
correlation value for a trace is set individually for each
detection based on:
min(detect_val / n_chans * min_cc_from_mean_cc_factor, min_cc).
:type horizontal_chans: list
:param horizontal_chans:
List of channel endings for horizontal-channels, on which S-picks
will be made.
:type vertical_chans: list
:param vertical_chans:
List of channel endings for vertical-channels, on which P-picks
will be made.
:type cores: int
:param cores:
Number of cores to use in parallel processing, defaults to one.
:type interpolate: bool
:param interpolate:
Interpolate the correlation function to achieve sub-sample
precision.
:type plot: bool
:param plot:
To generate a plot for every detection or not, defaults to False
:type plotdir: str
:param plotdir:
The path to save plots to. If `plotdir=None` (default) then the
figure will be shown on screen.
:type export_cc: bool
:param export_cc:
To generate a binary file in NumPy for every detection or not,
defaults to False
:type cc_dir: str
:param cc_dir:
Path to saving folder, NumPy files will be output here.
:type parallel: bool
:param parallel: Turn parallel processing on or off.
:type process_cores: int
:param process_cores:
Number of processes to use for pre-processing (if different to
`cores`).
:type ignore_length: bool
:param ignore_length:
If using daylong=True, then dayproc will try check that the data
are there for at least 80% of the day, if you don't want this check
(which will raise an error if too much data are missing) then set
ignore_length=True. This is not recommended!
:type ignore_bad_data: bool
:param ignore_bad_data:
If False (default), errors will be raised if data are excessively
gappy or are mostly zeros. If True then no error will be raised,
but an empty trace will be returned (and not used in detection).
:returns:
Catalog of events with picks. No origin information is included.
These events can then be written out via
:func:`obspy.core.event.Catalog.write`, or to Nordic Sfiles using
:func:`eqcorrscan.utils.sfile_util.eventtosfile` and located
externally.
:rtype: obspy.core.event.Catalog
.. Note::
Note on pre-processing: You can provide a pre-processed stream,
which may be beneficial for detections over large time periods
(the stream can have gaps, which reduces memory usage). However,
in this case the processing steps are not checked, so you must
ensure that all the template in the Party have the same sampling
rate and filtering as the stream.
If pre-processing has not be done then the data will be processed
according to the parameters in the templates, in this case
templates will be grouped by processing parameters and run with
similarly processed data. In this case, all templates do not have
to have the same processing parameters.
.. Note::
Picks are corrected for the template pre-pick time.
"""
process_cores = process_cores or cores
template_groups = group_templates(
[_f.template for _f in self.families
if len(_f) > 0]) # Fix for #341
catalog = Catalog()
for template_group in template_groups:
family = [_f for _f in self.families
if _f.template == template_group[0]][0]
group_seed_ids = {tr.id for template in template_group
for tr in template.st}
template_stream = Stream()
for seed_id in group_seed_ids:
net, sta, loc, chan = seed_id.split('.')
template_stream += stream.select(
network=net, station=sta, location=loc, channel=chan)
# Process once and only once for each group.
processed_stream = family._process_streams(
stream=template_stream, pre_processed=pre_processed,
process_cores=process_cores, parallel=parallel,
ignore_bad_data=ignore_bad_data, ignore_length=ignore_length,
select_used_chans=False)
for template in template_group:
family = [_f for _f in self.families
if _f.template == template][0]
catalog += family.lag_calc(
stream=processed_stream, pre_processed=True,
shift_len=shift_len, min_cc=min_cc,
min_cc_from_mean_cc_factor=min_cc_from_mean_cc_factor,
horizontal_chans=horizontal_chans,
vertical_chans=vertical_chans, cores=cores,
interpolate=interpolate, plot=plot, plotdir=plotdir,
export_cc=export_cc, cc_dir=cc_dir,
parallel=parallel, process_cores=process_cores,
ignore_bad_data=ignore_bad_data,
ignore_length=ignore_length, **kwargs)
return catalog
@staticmethod
def relative_magnitudes(*args, **kwargs):
print("This function is not functional, please keep an eye out for "
"this in future releases.")
# def relative_magnitudes(self, stream, pre_processed, process_cores=1,
# ignore_bad_data=False, parallel=False,
# min_cc=0.4, **kwargs):
# """
# Compute relative magnitudes for the detections.
#
# Works in place on events in the Family
#
# :type stream: obspy.core.stream.Stream
# :param stream:
# All the data needed to cut from - can be a gappy Stream.
# :type pre_processed: bool
# :param pre_processed:
# Whether the stream has been pre-processed or not to match the
# templates. See note below.
# :param parallel: Turn parallel processing on or off.
# :type process_cores: int
# :param process_cores:
# Number of processes to use for pre-processing (if different to
# `cores`).
# :type ignore_bad_data: bool
# :param ignore_bad_data:
# If False (default), errors will be raised if data are excessively
# gappy or are mostly zeros. If True then no error will be raised,
# but an empty trace will be returned (and not used in detection).
# :type min_cc: float
# :param min_cc: Minimum correlation for magnitude to be computed.
# :param kwargs:
# Keyword arguments passed to `utils.mag_calc.relative_mags`
#
# .. Note::
# Note on pre-processing: You can provide a pre-processed stream,
# which may be beneficial for detections over large time periods
# (the stream can have gaps, which reduces memory usage). However,
# in this case the processing steps are not checked, so you must
# ensure that the template in the Family has the same sampling
# rate and filtering as the stream.
# If pre-processing has not be done then the data will be processed
# according to the parameters in the template.
# """
# template_groups = group_templates(
# [_f.template for _f in self.families])
# for template_group in template_groups:
# family = [_f for _f in self.families
# if _f.template == template_group[0]][0]
# processed_stream = family._process_streams(
# stream=stream, pre_processed=pre_processed,
# process_cores=process_cores, parallel=parallel,
# ignore_bad_data=ignore_bad_data)
# for template in template_group:
# family = [_f for _f in self.families
# if _f.template == template][0]
# family.relative_magnitudes(
# stream=processed_stream, pre_processed=True,
# min_cc=min_cc, parallel=parallel,
# process_cores=process_cores,
# ignore_bad_data=ignore_bad_data,
# **kwargs)
# return self.get_catalog()
[docs]
def get_catalog(self):
"""
Get an obspy catalog object from the party.
:returns: :class:`obspy.core.event.Catalog`
.. rubric:: Example
>>> party = Party().read()
>>> cat = party.get_catalog()
>>> print(len(cat))
4
"""
catalog = Catalog()
for fam in self.families:
if len(fam.catalog) != 0:
catalog.events.extend(fam.catalog.events)
return catalog
[docs]
def min_chans(self, min_chans):
"""
Remove detections using min_chans or fewer channels.
:type min_chans: int
:param min_chans:
Detections using more than this number of channels are maintained.
Note that this is a strict `if detection.no_chans > min_chans:`
rather than >=. Maintained for backwards compatability.
:return: Party
.. Note:: Works in place on Party.
.. rubric:: Example
>>> party = Party().read()
>>> print(len(party))
4
>>> party = party.min_chans(5)
>>> print(len(party))
1
"""
declustered = Party()
for family in self.families:
fam = Family(family.template)
for d in family.detections:
if d.no_chans > min_chans:
fam.detections.append(d)
declustered.families.append(fam)
self.families = declustered.families
return self
[docs]
def read_party(fname=None, read_detection_catalog=True, *args, **kwargs):
"""
Read detections and metadata from a tar archive.
:type fname: str
:param fname:
Filename to read from, if this contains a single Family, then will
return a party of length = 1
:type read_detection_catalog: bool
:param read_detection_catalog:
Whether to read the detection catalog or not, if False, catalog
will be regenerated - for large catalogs this can be faster.
:return: :class:`eqcorrscan.core.match_filter.Party`
"""
party = Party()
party.read(filename=fname, read_detection_catalog=read_detection_catalog,
*args, **kwargs)
return party
if __name__ == "__main__":
import doctest
doctest.testmod()
# List files to be removed after doctest
cleanup = ['test_tar_write.tgz', 'test_csv_write.csv', 'test_quakeml.xml']
for f in cleanup:
if os.path.isfile(f):
os.remove(f)
elif os.path.isdir(f):
shutil.rmtree(f)