1#!/usr/bin/env python 2# -*- coding: utf-8 -*- 3# ------------------------------------------------------------------- 4# Filename: cross_correlation.py 5# Author: Moritz Beyreuther, Tobias Megies, Tom Eulenfeld 6# Email: megies@geophysik.uni-muenchen.de 7# 8# Copyright (C) 2008-2019 Moritz Beyreuther, Tobias Megies, Tom Eulenfeld 9# ------------------------------------------------------------------ 10""" 11Signal processing routines based on cross correlation techniques. 12 13:copyright: 14 The ObsPy Development Team (devs@obspy.org) 15:license: 16 GNU Lesser General Public License, Version 3 17 (https://www.gnu.org/copyleft/lesser.html) 18""" 19from __future__ import (absolute_import, division, print_function, 20 unicode_literals) 21from future.builtins import * # NOQA 22from future.utils import native_str 23 24from bisect import bisect_left 25from copy import copy 26import ctypes as C # NOQA 27from distutils.version import LooseVersion 28import warnings 29 30import numpy as np 31import scipy 32 33from obspy import Stream, Trace 34from obspy.core.util.misc import MatplotlibBackend 35from obspy.signal.headers import clibsignal 36from obspy.signal.invsim import cosine_taper 37 38 39def _pad_zeros(a, num, num2=None): 40 """Pad num zeros at both sides of array a""" 41 if num2 is None: 42 num2 = num 43 hstack = [np.zeros(num, dtype=a.dtype), a, np.zeros(num2, dtype=a.dtype)] 44 return np.hstack(hstack) 45 46 47def _call_scipy_correlate(a, b, mode, method): 48 """ 49 Call the correct correlate function depending on Scipy version and method. 50 """ 51 if LooseVersion(scipy.__version__) >= LooseVersion('0.19'): 52 cc = scipy.signal.correlate(a, b, mode=mode, method=method) 53 elif method in ('fft', 'auto'): 54 cc = scipy.signal.fftconvolve(a, b[::-1], mode=mode) 55 elif method == 'direct': 56 cc = scipy.signal.correlate(a, b, mode=mode) 57 else: 58 msg = "method keyword has to be one of ('auto', 'fft', 'direct')" 59 raise ValueError(msg) 60 return cc 61 62 63def _xcorr_padzeros(a, b, shift, method): 64 """ 65 Cross-correlation using SciPy with mode='valid' and precedent zero padding. 66 """ 67 if shift is None: 68 shift = (len(a) + len(b) - 1) // 2 69 dif = len(a) - len(b) - 2 * shift 70 if dif > 0: 71 b = _pad_zeros(b, dif // 2) 72 else: 73 a = _pad_zeros(a, -dif // 2) 74 return _call_scipy_correlate(a, b, 'valid', method) 75 76 77def _xcorr_slice(a, b, shift, method): 78 """ 79 Cross-correlation using SciPy with mode='full' and subsequent slicing. 80 """ 81 mid = (len(a) + len(b) - 1) // 2 82 if shift is None: 83 shift = mid 84 if shift > mid: 85 # Such a large shift is not possible without zero padding 86 return _xcorr_padzeros(a, b, shift, method) 87 cc = _call_scipy_correlate(a, b, 'full', method) 88 return cc[mid - shift:mid + shift + len(cc) % 2] 89 90 91def correlate(a, b, shift, demean=True, normalize='naive', method='auto', 92 domain=None): 93 """ 94 Cross-correlation of two signals up to a specified maximal shift. 95 96 This function only allows 'naive' normalization with the overall 97 standard deviations. This is a reasonable approximation for signals of 98 similar length and a relatively small shift parameter 99 (e.g. noise cross-correlation). 100 If you are interested in the full cross-correlation function better use 101 :func:`~obspy.signal.cross_correlation.correlate_template` which also 102 provides correct normalization. 103 104 :type a: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 105 :param a: first signal 106 :type b: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 107 :param b: second signal to correlate with first signal 108 :param int shift: Number of samples to shift for cross correlation. 109 The cross-correlation will consist of ``2*shift+1`` or 110 ``2*shift`` samples. The sample with zero shift will be in the middle. 111 :param bool demean: Demean data beforehand. 112 :param normalize: Method for normalization of cross-correlation. 113 One of ``'naive'`` or ``None`` 114 (``True`` and ``False`` are supported for backwards compatibility). 115 ``'naive'`` normalizes by the overall standard deviation. 116 ``None`` does not normalize. 117 :param str method: Method to use to calculate the correlation. 118 ``'direct'``: The correlation is determined directly from sums, 119 the definition of correlation. 120 ``'fft'`` The Fast Fourier Transform is used to perform the 121 correlation more quickly. 122 ``'auto'`` Automatically chooses direct or Fourier method based on an 123 estimate of which is faster. (Only availlable for SciPy versions >= 124 0.19. For older Scipy version method defaults to ``'fft'``.) 125 :param str domain: Deprecated. Please use the method argument. 126 127 :return: cross-correlation function. 128 129 To calculate shift and value of the maximum of the returned 130 cross-correlation function use 131 :func:`~obspy.signal.cross_correlation.xcorr_max`. 132 133 .. note:: 134 135 For most input parameters cross-correlation using the FFT is much 136 faster. 137 Only for small values of ``shift`` (approximately less than 100) 138 direct time domain cross-correlation migth save some time. 139 140 .. note:: 141 142 If the signals have different length, they will be aligned around 143 their middle. The sample with zero shift in the cross-correlation 144 function corresponds to this correlation: 145 146 :: 147 148 --aaaa-- 149 bbbbbbbb 150 151 For odd ``len(a)-len(b)`` the cross-correlation function will 152 consist of only ``2*shift`` samples because a shift of 0 153 corresponds to the middle between two samples. 154 155 .. rubric:: Example 156 157 >>> from obspy import read 158 >>> a = read()[0][450:550] 159 >>> b = a[:-2] 160 >>> cc = correlate(a, b, 2) 161 >>> cc 162 array([ 0.62390515, 0.99630851, 0.62187106, -0.05864797, -0.41496995]) 163 >>> shift, value = xcorr_max(cc) 164 >>> shift 165 -1 166 >>> round(value, 3) 167 0.996 168 """ 169 if normalize is False: 170 normalize = None 171 if normalize is True: 172 normalize = 'naive' 173 if domain is not None: 174 if domain == 'freq': 175 method = 'fft' 176 elif domain == 'time': 177 method = 'direct' 178 from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning 179 msg = ("'domain' keyword of correlate function is deprecated and will " 180 "be removed in a subsequent ObsPy release. " 181 "Please use the 'method' keyword.") 182 warnings.warn(msg, ObsPyDeprecationWarning) 183 # if we get Trace objects, use their data arrays 184 if isinstance(a, Trace): 185 a = a.data 186 if isinstance(b, Trace): 187 b = b.data 188 a = np.asarray(a) 189 b = np.asarray(b) 190 if demean: 191 a = a - np.mean(a) 192 b = b - np.mean(b) 193 # choose the usually faster xcorr function for each method 194 _xcorr = _xcorr_padzeros if method == 'direct' else _xcorr_slice 195 cc = _xcorr(a, b, shift, method) 196 if normalize == 'naive': 197 norm = (np.sum(a ** 2) * np.sum(b ** 2)) ** 0.5 198 if norm <= np.finfo(float).eps: 199 # norm is zero 200 # => cross-correlation function will have only zeros 201 cc[:] = 0 202 elif cc.dtype == float: 203 cc /= norm 204 else: 205 cc = cc / norm 206 elif normalize is not None: 207 raise ValueError("normalize has to be one of (None, 'naive'))") 208 return cc 209 210 211def _window_sum(data, window_len): 212 """Rolling sum of data.""" 213 window_sum = np.cumsum(data) 214 # in-place equivalent of 215 # window_sum = window_sum[window_len:] - window_sum[:-window_len] 216 # return window_sum 217 np.subtract(window_sum[window_len:], window_sum[:-window_len], 218 out=window_sum[:-window_len]) 219 return window_sum[:-window_len] 220 221 222def correlate_template(data, template, mode='valid', normalize='full', 223 demean=True, method='auto'): 224 """ 225 Normalized cross-correlation of two signals with specified mode. 226 227 If you are interested only in a part of the cross-correlation function 228 around zero shift consider using function 229 :func:`~obspy.signal.cross_correlation.correlate` which allows to 230 explicetly specify the maximum shift. 231 232 :type data: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 233 :param data: first signal 234 :type template: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 235 :param template: second signal to correlate with first signal. 236 Its length must be smaller or equal to the length of ``data``. 237 :param str mode: correlation mode to use. 238 It is passed to the used correlation function. 239 See :func:`scipy.signal.correlate` for possible options. 240 The parameter determines the length of the correlation function. 241 :param normalize: 242 One of ``'naive'``, ``'full'`` or ``None``. 243 ``'full'`` normalizes every correlation properly, 244 whereas ``'naive'`` normalizes by the overall standard deviations. 245 ``None`` does not normalize. 246 :param demean: Demean data beforehand. For ``normalize='full'`` data is 247 demeaned in different windows for each correlation value. 248 :param str method: Method to use to calculate the correlation. 249 ``'direct'``: The correlation is determined directly from sums, 250 the definition of correlation. 251 ``'fft'`` The Fast Fourier Transform is used to perform the 252 correlation more quickly. 253 ``'auto'`` Automatically chooses direct or Fourier method based on an 254 estimate of which is faster. (Only availlable for SciPy versions >= 255 0.19. For older Scipy version method defaults to ``'fft'``.) 256 257 :return: cross-correlation function. 258 259 .. note:: 260 Calling the function with ``demean=True, normalize='full'`` (default) 261 returns the zero-normalized cross-correlation function. 262 Calling the function with ``demean=False, normalize='full'`` 263 returns the normalized cross-correlation function. 264 265 .. rubric:: Example 266 267 >>> from obspy import read 268 >>> data = read()[0] 269 >>> template = data[450:550] 270 >>> cc = correlate_template(data, template) 271 >>> index = np.argmax(cc) 272 >>> index 273 450 274 >>> round(cc[index], 9) 275 1.0 276 """ 277 # if we get Trace objects, use their data arrays 278 if isinstance(data, Trace): 279 data = data.data 280 if isinstance(template, Trace): 281 template = template.data 282 data = np.asarray(data) 283 template = np.asarray(template) 284 lent = len(template) 285 if len(data) < lent: 286 raise ValueError('Data must not be shorter than template.') 287 if demean: 288 template = template - np.mean(template) 289 if normalize != 'full': 290 data = data - np.mean(data) 291 cc = _call_scipy_correlate(data, template, mode, method) 292 if normalize is not None: 293 tnorm = np.sum(template ** 2) 294 if normalize == 'naive': 295 norm = (tnorm * np.sum(data ** 2)) ** 0.5 296 if norm <= np.finfo(float).eps: 297 cc[:] = 0 298 elif cc.dtype == float: 299 cc /= norm 300 else: 301 cc = cc / norm 302 elif normalize == 'full': 303 pad = len(cc) - len(data) + lent 304 if mode == 'same': 305 pad1, pad2 = (pad + 2) // 2, (pad - 1) // 2 306 else: 307 pad1, pad2 = (pad + 1) // 2, pad // 2 308 data = _pad_zeros(data, pad1, pad2) 309 # in-place equivalent of 310 # if demean: 311 # norm = ((_window_sum(data ** 2, lent) - 312 # _window_sum(data, lent) ** 2 / lent) * tnorm) ** 0.5 313 # else: 314 # norm = (_window_sum(data ** 2, lent) * tnorm) ** 0.5 315 # cc = cc / norm 316 if demean: 317 norm = _window_sum(data, lent) ** 2 318 if norm.dtype == float: 319 norm /= lent 320 else: 321 norm = norm / lent 322 np.subtract(_window_sum(data ** 2, lent), norm, out=norm) 323 else: 324 norm = _window_sum(data ** 2, lent) 325 norm *= tnorm 326 if norm.dtype == float: 327 np.sqrt(norm, out=norm) 328 else: 329 norm = np.sqrt(norm) 330 mask = norm <= np.finfo(float).eps 331 if cc.dtype == float: 332 cc[~mask] /= norm[~mask] 333 else: 334 cc = cc / norm 335 cc[mask] = 0 336 else: 337 msg = "normalize has to be one of (None, 'naive', 'full')" 338 raise ValueError(msg) 339 return cc 340 341 342def xcorr(tr1, tr2, shift_len, full_xcorr=False): 343 """ 344 Cross correlation of tr1 and tr2 in the time domain using window_len. 345 346 .. note:: 347 Please use the :func:`~obspy.signal.cross_correlation.correlate` 348 function for new code. 349 350 :: 351 352 Mid Sample 353 | 354 |AAAAAAAAAAAAAAA|AAAAAAAAAAAAAAA|AAAAAAAAAAAAAAA|AAAAAAAAAAAAAAA| 355 |BBBBBBBBBBBBBBB|BBBBBBBBBBBBBBB|BBBBBBBBBBBBBBB|BBBBBBBBBBBBBBB| 356 |<-shift_len/2->| <- region of support -> |<-shift_len/2->| 357 358 359 :type tr1: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 360 :param tr1: Trace 1 361 :type tr2: :class:`~numpy.ndarray`, :class:`~obspy.core.trace.Trace` 362 :param tr2: Trace 2 to correlate with trace 1 363 :type shift_len: int 364 :param shift_len: Total length of samples to shift for cross correlation. 365 :type full_xcorr: bool 366 :param full_xcorr: If ``True``, the complete xcorr function will be 367 returned as :class:`~numpy.ndarray` 368 :return: **index, value[, fct]** - Index of maximum xcorr value and the 369 value itself. The complete xcorr function is returned only if 370 ``full_xcorr=True``. 371 372 .. note:: 373 As shift_len gets higher the window supporting the cross correlation 374 actually gets smaller. So with shift_len=0 you get the correlation 375 coefficient of both traces as a whole without any shift applied. As the 376 xcorr function works in time domain and does not zero pad at all, with 377 higher shifts allowed the window of support gets smaller so that the 378 moving windows shifted against each other do not run out of the 379 timeseries bounds at high time shifts. Of course there are other 380 possibilities to do cross correlations e.g. in frequency domain. 381 382 .. seealso:: 383 `ObsPy-users mailing list 384 <http://lists.obspy.org/pipermail/obspy-users/2011-March/000056.html>`_ 385 and `issue #249 <https://github.com/obspy/obspy/issues/249>`_. 386 """ 387 from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning 388 msg = ('Call to deprecated function xcorr(). Please use the correlate and ' 389 'xcorr_max functions.') 390 warnings.warn(msg, ObsPyDeprecationWarning) 391 392 # if we get Trace objects, use their data arrays 393 for tr in [tr1, tr2]: 394 if isinstance(tr, Trace): 395 tr = tr.data 396 397 # check if shift_len parameter is in an acceptable range. 398 # if not the underlying c code tampers with shift_len and uses shift_len/2 399 # instead. we want to avoid this silent automagic and raise an error in the 400 # python layer right here. 401 # see ticket #249 and src/xcorr.c lines 43-57 402 if min(len(tr1), len(tr2)) - 2 * shift_len <= 0: 403 msg = "shift_len too large. The underlying C code would silently " + \ 404 "use shift_len/2 which we want to avoid." 405 raise ValueError(msg) 406 # be nice and adapt type if necessary 407 tr1 = np.ascontiguousarray(tr1, np.float32) 408 tr2 = np.ascontiguousarray(tr2, np.float32) 409 corp = np.empty(2 * shift_len + 1, dtype=np.float64, order='C') 410 411 shift = C.c_int() 412 coe_p = C.c_double() 413 414 res = clibsignal.X_corr(tr1, tr2, corp, shift_len, len(tr1), len(tr2), 415 C.byref(shift), C.byref(coe_p)) 416 if res: 417 raise MemoryError 418 419 if full_xcorr: 420 return shift.value, coe_p.value, corp 421 else: 422 return shift.value, coe_p.value 423 424 425def xcorr_3c(st1, st2, shift_len, components=["Z", "N", "E"], 426 full_xcorr=False, abs_max=True): 427 """ 428 Calculates the cross correlation on each of the specified components 429 separately, stacks them together and estimates the maximum and shift of 430 maximum on the stack. 431 432 Basically the same as :func:`~obspy.signal.cross_correlation.xcorr` but 433 for (normally) three components, please also take a look at the 434 documentation of that function. Useful e.g. for estimation of waveform 435 similarity on a three component seismogram. 436 437 :type st1: :class:`~obspy.core.stream.Stream` 438 :param st1: Stream 1, containing one trace for Z, N, E component (other 439 component_id codes are ignored) 440 :type st2: :class:`~obspy.core.stream.Stream` 441 :param st2: Stream 2, containing one trace for Z, N, E component (other 442 component_id codes are ignored) 443 :type shift_len: int 444 :param shift_len: Total length of samples to shift for cross correlation. 445 :type components: list of str 446 :param components: List of components to use in cross-correlation, defaults 447 to ``['Z', 'N', 'E']``. 448 :type full_xcorr: bool 449 :param full_xcorr: If ``True``, the complete xcorr function will be 450 returned as :class:`~numpy.ndarray`. 451 :param bool abs_max: *shift* will be calculated for maximum or 452 absolute maximum. 453 :return: **index, value[, fct]** - index of maximum xcorr value and the 454 value itself. The complete xcorr function is returned only if 455 ``full_xcorr=True``. 456 """ 457 streams = [st1, st2] 458 # check if we can actually use the provided streams safely 459 for st in streams: 460 if not isinstance(st, Stream): 461 raise TypeError("Expected Stream object but got %s." % type(st)) 462 for component in components: 463 if not len(st.select(component=component)) == 1: 464 msg = "Expected exactly one %s trace in stream" % component + \ 465 " but got %s." % len(st.select(component=component)) 466 raise ValueError(msg) 467 ndat = len(streams[0].select(component=components[0])[0]) 468 if False in [len(st.select(component=component)[0]) == ndat 469 for st in streams for component in components]: 470 raise ValueError("All traces have to be the same length.") 471 # everything should be ok with the input data... 472 corp = np.zeros(2 * shift_len + 1, dtype=np.float64, order='C') 473 for component in components: 474 xx = correlate(streams[0].select(component=component)[0], 475 streams[1].select(component=component)[0], 476 shift_len) 477 corp += xx 478 corp /= len(components) 479 shift, value = xcorr_max(corp, abs_max=abs_max) 480 if full_xcorr: 481 return shift, value, corp 482 else: 483 return shift, value 484 485 486def xcorr_max(fct, abs_max=True): 487 """ 488 Return shift and value of the maximum of the cross-correlation function. 489 490 :type fct: :class:`~numpy.ndarray` 491 :param fct: Cross-correlation function e.g. returned by correlate. 492 :param bool abs_max: Determines if the absolute maximum should be used. 493 :return: **shift, value** - Shift and value of maximum of 494 cross-correlation. 495 496 .. rubric:: Example 497 498 >>> fct = np.zeros(101) 499 >>> fct[50] = -1.0 500 >>> xcorr_max(fct) 501 (0, -1.0) 502 >>> fct[50], fct[60] = 0.0, 1.0 503 >>> xcorr_max(fct) 504 (10, 1.0) 505 >>> fct[60], fct[40] = 0.0, -1.0 506 >>> xcorr_max(fct) 507 (-10, -1.0) 508 >>> fct[60], fct[40] = 0.5, -1.0 509 >>> xcorr_max(fct, abs_max=True) 510 (-10, -1.0) 511 >>> xcorr_max(fct, abs_max=False) 512 (10, 0.5) 513 >>> xcorr_max(fct[:-1], abs_max=False) 514 (10.5, 0.5) 515 """ 516 mid = (len(fct) - 1) / 2 517 if len(fct) % 2 == 1: 518 mid = int(mid) 519 index = np.argmax(np.abs(fct) if abs_max else fct) 520 # float() call is workaround for future package 521 # see https://travis-ci.org/obspy/obspy/jobs/174284750 522 return index - mid, float(fct[index]) 523 524 525def xcorr_pick_correction(pick1, trace1, pick2, trace2, t_before, t_after, 526 cc_maxlag, filter=None, filter_options={}, 527 plot=False, filename=None): 528 """ 529 Calculate the correction for the differential pick time determined by cross 530 correlation of the waveforms in narrow windows around the pick times. 531 For details on the fitting procedure refer to [Deichmann1992]_. 532 533 The parameters depend on the epicentral distance and magnitude range. For 534 small local earthquakes (Ml ~0-2, distance ~3-10 km) with consistent manual 535 picks the following can be tried:: 536 537 t_before=0.05, t_after=0.2, cc_maxlag=0.10, 538 filter="bandpass", filter_options={'freqmin': 1, 'freqmax': 20} 539 540 The appropriate parameter sets can and should be determined/verified 541 visually using the option `plot=True` on a representative set of picks. 542 543 To get the corrected differential pick time calculate: ``((pick2 + 544 pick2_corr) - pick1)``. To get a corrected differential travel time using 545 origin times for both events calculate: ``((pick2 + pick2_corr - ot2) - 546 (pick1 - ot1))`` 547 548 :type pick1: :class:`~obspy.core.utcdatetime.UTCDateTime` 549 :param pick1: Time of pick for `trace1`. 550 :type trace1: :class:`~obspy.core.trace.Trace` 551 :param trace1: Waveform data for `pick1`. Add some time at front/back. 552 The appropriate part of the trace is used automatically. 553 :type pick2: :class:`~obspy.core.utcdatetime.UTCDateTime` 554 :param pick2: Time of pick for `trace2`. 555 :type trace2: :class:`~obspy.core.trace.Trace` 556 :param trace2: Waveform data for `pick2`. Add some time at front/back. 557 The appropriate part of the trace is used automatically. 558 :type t_before: float 559 :param t_before: Time to start cross correlation window before pick times 560 in seconds. 561 :type t_after: float 562 :param t_after: Time to end cross correlation window after pick times in 563 seconds. 564 :type cc_maxlag: float 565 :param cc_maxlag: Maximum lag/shift time tested during cross correlation in 566 seconds. 567 :type filter: str 568 :param filter: `None` for no filtering or name of filter type 569 as passed on to :meth:`~obspy.core.Trace.trace.filter` if filter 570 should be used. To avoid artifacts in filtering provide 571 sufficiently long time series for `trace1` and `trace2`. 572 :type filter_options: dict 573 :param filter_options: Filter options that get passed on to 574 :meth:`~obspy.core.Trace.trace.filter` if filtering is used. 575 :type plot: bool 576 :param plot: If `True`, a plot window illustrating the alignment of the two 577 traces at best cross correlation will be shown. This can and should be 578 used to verify the used parameters before running automatedly on large 579 data sets. 580 :type filename: str 581 :param filename: If plot option is selected, specifying a filename here 582 (e.g. 'myplot.pdf' or 'myplot.png') will output the plot to a file 583 instead of opening a plot window. 584 :rtype: (float, float) 585 :returns: Correction time `pick2_corr` for `pick2` pick time as a float and 586 corresponding correlation coefficient. 587 """ 588 # perform some checks on the traces 589 if trace1.stats.sampling_rate != trace2.stats.sampling_rate: 590 msg = "Sampling rates do not match: %s != %s" % \ 591 (trace1.stats.sampling_rate, trace2.stats.sampling_rate) 592 raise Exception(msg) 593 if trace1.id != trace2.id: 594 msg = "Trace ids do not match: %s != %s" % (trace1.id, trace2.id) 595 warnings.warn(msg) 596 samp_rate = trace1.stats.sampling_rate 597 # don't modify existing traces with filters 598 if filter: 599 trace1 = trace1.copy() 600 trace2 = trace2.copy() 601 # check data, apply filter and take correct slice of traces 602 slices = [] 603 for _i, (t, tr) in enumerate(((pick1, trace1), (pick2, trace2))): 604 start = t - t_before - (cc_maxlag / 2.0) 605 end = t + t_after + (cc_maxlag / 2.0) 606 duration = end - start 607 # check if necessary time spans are present in data 608 if tr.stats.starttime > start: 609 msg = "Trace %s starts too late." % _i 610 raise Exception(msg) 611 if tr.stats.endtime < end: 612 msg = "Trace %s ends too early." % _i 613 raise Exception(msg) 614 if filter and start - tr.stats.starttime < duration: 615 msg = "Artifacts from signal processing possible. Trace " + \ 616 "%s should have more additional data at the start." % _i 617 warnings.warn(msg) 618 if filter and tr.stats.endtime - end < duration: 619 msg = "Artifacts from signal processing possible. Trace " + \ 620 "%s should have more additional data at the end." % _i 621 warnings.warn(msg) 622 # apply signal processing and take correct slice of data 623 if filter: 624 tr.data = tr.data.astype(np.float64) 625 tr.detrend(type='demean') 626 tr.data *= cosine_taper(len(tr), 0.1) 627 tr.filter(type=filter, **filter_options) 628 slices.append(tr.slice(start, end)) 629 # cross correlate 630 shift_len = int(cc_maxlag * samp_rate) 631 cc = correlate(slices[0].data, slices[1].data, shift_len, method='direct') 632 _cc_shift, cc_max = xcorr_max(cc) 633 cc_curvature = np.concatenate((np.zeros(1), np.diff(cc, 2), np.zeros(1))) 634 cc_convex = np.ma.masked_where(np.sign(cc_curvature) >= 0, cc) 635 cc_concave = np.ma.masked_where(np.sign(cc_curvature) < 0, cc) 636 # check results of cross correlation 637 if cc_max < 0: 638 msg = "Absolute maximum is negative: %.3f. " % cc_max + \ 639 "Using positive maximum: %.3f" % max(cc) 640 warnings.warn(msg) 641 cc_max = max(cc) 642 if cc_max < 0.8: 643 msg = "Maximum of cross correlation lower than 0.8: %s" % cc_max 644 warnings.warn(msg) 645 # make array with time shifts in seconds corresponding to cc function 646 cc_t = np.linspace(-cc_maxlag, cc_maxlag, shift_len * 2 + 1) 647 # take the subportion of the cross correlation around the maximum that is 648 # convex and fit a parabola. 649 # use vertex as subsample resolution best cc fit. 650 peak_index = cc.argmax() 651 first_sample = peak_index 652 # XXX this could be improved.. 653 while first_sample > 0 and cc_curvature[first_sample - 1] <= 0: 654 first_sample -= 1 655 last_sample = peak_index 656 while last_sample < len(cc) - 1 and cc_curvature[last_sample + 1] <= 0: 657 last_sample += 1 658 if first_sample == 0 or last_sample == len(cc) - 1: 659 msg = "Fitting at maximum lag. Maximum lag time should be increased." 660 warnings.warn(msg) 661 # work on subarrays 662 num_samples = last_sample - first_sample + 1 663 if num_samples < 3: 664 msg = "Less than 3 samples selected for fit to cross " + \ 665 "correlation: %s" % num_samples 666 raise Exception(msg) 667 if num_samples < 5: 668 msg = "Less than 5 samples selected for fit to cross " + \ 669 "correlation: %s" % num_samples 670 warnings.warn(msg) 671 # quadratic fit for small subwindow 672 coeffs, residual = scipy.polyfit( 673 cc_t[first_sample:last_sample + 1], 674 cc[first_sample:last_sample + 1], deg=2, full=True)[:2] 675 # check results of fit 676 if coeffs[0] >= 0: 677 msg = "Fitted parabola opens upwards!" 678 warnings.warn(msg) 679 if residual > 0.1: 680 msg = "Residual in quadratic fit to cross correlation maximum " + \ 681 "larger than 0.1: %s" % residual 682 warnings.warn(msg) 683 # X coordinate of vertex of parabola gives time shift to correct 684 # differential pick time. Y coordinate gives maximum correlation 685 # coefficient. 686 dt = -coeffs[1] / 2.0 / coeffs[0] 687 coeff = (4 * coeffs[0] * coeffs[2] - coeffs[1] ** 2) / (4 * coeffs[0]) 688 # this is the shift to apply on the time axis of `trace2` to align the 689 # traces. Actually we do not want to shift the trace to align it but we 690 # want to correct the time of `pick2` so that the traces align without 691 # shifting. This is the negative of the cross correlation shift. 692 dt = -dt 693 pick2_corr = dt 694 # plot the results if selected 695 if plot is True: 696 with MatplotlibBackend(filename and "AGG" or None, sloppy=True): 697 import matplotlib.pyplot as plt 698 fig = plt.figure() 699 ax1 = fig.add_subplot(211) 700 tmp_t = np.linspace(0, len(slices[0]) / samp_rate, len(slices[0])) 701 ax1.plot(tmp_t, slices[0].data / float(slices[0].data.max()), "k", 702 label="Trace 1") 703 ax1.plot(tmp_t, slices[1].data / float(slices[1].data.max()), "r", 704 label="Trace 2") 705 ax1.plot(tmp_t - dt, slices[1].data / float(slices[1].data.max()), 706 "g", label="Trace 2 (shifted)") 707 ax1.legend(loc="lower right", prop={'size': "small"}) 708 ax1.set_title("%s" % slices[0].id) 709 ax1.set_xlabel("time [s]") 710 ax1.set_ylabel("norm. amplitude") 711 ax2 = fig.add_subplot(212) 712 ax2.plot(cc_t, cc_convex, ls="", marker=".", color="k", 713 label="xcorr (convex)") 714 ax2.plot(cc_t, cc_concave, ls="", marker=".", color="0.7", 715 label="xcorr (concave)") 716 ax2.plot(cc_t[first_sample:last_sample + 1], 717 cc[first_sample:last_sample + 1], "b.", 718 label="used for fitting") 719 tmp_t = np.linspace(cc_t[first_sample], cc_t[last_sample], 720 num_samples * 10) 721 ax2.plot(tmp_t, scipy.polyval(coeffs, tmp_t), "b", label="fit") 722 ax2.axvline(-dt, color="g", label="vertex") 723 ax2.axhline(coeff, color="g") 724 ax2.set_xlabel("%.2f at %.3f seconds correction" % (coeff, -dt)) 725 ax2.set_ylabel("correlation coefficient") 726 ax2.set_ylim(-1, 1) 727 ax2.set_xlim(cc_t[0], cc_t[-1]) 728 ax2.legend(loc="lower right", prop={'size': "x-small"}) 729 # plt.legend(loc="lower left") 730 if filename: 731 fig.savefig(filename) 732 else: 733 plt.show() 734 735 return (pick2_corr, coeff) 736 737 738def templates_max_similarity(st, time, streams_templates): 739 """ 740 Compares all event templates in the streams_templates list of streams 741 against the given stream around the time of the suspected event. The stream 742 that is being checked has to include all trace ids that are included in 743 template events. One component streams can be checked as well as multiple 744 components simultaneously. In case of multiple components it is made sure, 745 that all three components are shifted together. The traces in any stream 746 need to have a reasonable common starting time. The stream to check should 747 have some additional data to left/right of suspected event, the event 748 template streams should be cut to the portion of the event that should be 749 compared. Also see :func:`obspy.signal.trigger.coincidence_trigger` and the 750 corresponding example in the 751 `Trigger/Picker Tutorial 752 <https://tutorial.obspy.org/code_snippets/trigger_tutorial.html>`_. 753 754 - computes cross correlation on each component (one stream serves as 755 template, one as a longer search stream) 756 - stack all three and determine best shift in stack 757 - normalization is a bit problematic so compute the correlation coefficient 758 afterwards for the best shift to make sure the result is between 0 and 1. 759 760 >>> from obspy import read, UTCDateTime 761 >>> import numpy as np 762 >>> np.random.seed(123) # make test reproducible 763 >>> st = read() 764 >>> t = UTCDateTime(2009, 8, 24, 0, 20, 7, 700000) 765 >>> templ = st.copy().slice(t, t+5) 766 >>> for tr in templ: 767 ... tr.data += np.random.random(len(tr)) * tr.data.max() * 0.5 768 >>> print(templates_max_similarity(st, t, [templ])) 769 0.922536411468 770 771 :param time: Time around which is checked for a similarity. Cross 772 correlation shifts of around template event length are checked. 773 :type time: :class:`~obspy.core.utcdatetime.UTCDateTime` 774 :param st: One or multi-component Stream to check against event templates. 775 Should have additional data left/right of suspected event (around half 776 the length of template events). 777 :type st: :class:`~obspy.core.stream.Stream` 778 :param streams_templates: List of streams with template events to check for 779 waveform similarity. Each template has to include data for all 780 channels present in stream to check. 781 :type streams_templates: list of :class:`~obspy.core.stream.Stream` 782 :returns: Best correlation coefficient obtained by the comparison against 783 all template events (0 to 1). 784 """ 785 values = [] 786 for st_tmpl in streams_templates: 787 ids = [tr.id for tr in st_tmpl] 788 duration = st_tmpl[0].stats.endtime - st_tmpl[0].stats.starttime 789 st_ = st.slice(time - (duration * 0.5), 790 time + (duration * 1.5)) 791 cc = None 792 for id_ in reversed(ids): 793 if not st_.select(id=id_): 794 msg = "Skipping trace %s in template correlation " + \ 795 "(not present in stream to check)." 796 warnings.warn(msg % id_) 797 ids.remove(id_) 798 if not ids: 799 msg = ("Skipping template(s) for station '{}': No common SEED IDs " 800 "when comparing template ({}) and data streams ({}).") 801 warnings.warn(msg.format( 802 st_tmpl[0].stats.station, 803 ', '.join(sorted(set(tr.id for tr in st_tmpl))), 804 ', '.join(sorted(set(tr.id for tr in st_))))) 805 continue 806 # determine best (combined) shift of multi-component data 807 for id_ in ids: 808 tr1 = st_.select(id=id_)[0] 809 tr2 = st_tmpl.select(id=id_)[0] 810 if len(tr1) > len(tr2): 811 data_short = tr2.data 812 data_long = tr1.data 813 else: 814 data_short = tr1.data 815 data_long = tr2.data 816 data_short = (data_short - data_short.mean()) / data_short.std() 817 data_long = (data_long - data_long.mean()) / data_long.std() 818 tmp = np.correlate(data_long, data_short, native_str("valid")) 819 try: 820 cc += tmp 821 except TypeError: 822 cc = tmp 823 except ValueError: 824 cc = None 825 break 826 if cc is None: 827 msg = "Skipping template(s) for station %s due to problems in " + \ 828 "three component correlation (gappy traces?)" 829 warnings.warn(msg % st_tmpl[0].stats.station) 830 continue 831 ind = cc.argmax() 832 ind2 = ind + len(data_short) 833 coef = 0.0 834 # determine correlation coefficient of best shift as the mean of all 835 # components 836 for id_ in ids: 837 tr1 = st_.select(id=id_)[0] 838 tr2 = st_tmpl.select(id=id_)[0] 839 if len(tr1) > len(tr2): 840 data_short = tr2.data 841 data_long = tr1.data 842 else: 843 data_short = tr1.data 844 data_long = tr2.data 845 coef += np.corrcoef(data_short, data_long[ind:ind2])[0, 1] 846 coef /= len(ids) 847 values.append(coef) 848 if values: 849 return max(values) 850 else: 851 return 0 852 853 854def _prep_streams_correlate(stream, template, template_time=None): 855 """ 856 Prepare stream and template for cross-correlation. 857 858 Select traces in stream and template with the same seed id and trim 859 stream to correct start and end times. 860 """ 861 if len({tr.stats.sampling_rate for tr in stream + template}) > 1: 862 raise ValueError('Traces have different sampling rate') 863 ids = {tr.id for tr in stream} & {tr.id for tr in template} 864 if len(ids) == 0: 865 raise ValueError('No traces with matching ids in template and stream') 866 stream = copy(stream) 867 template = copy(template) 868 stream.traces = [tr for tr in stream if tr.id in ids] 869 template.traces = [tr for tr in template if tr.id in ids] 870 template.sort() 871 stream.sort() 872 if len(stream) != len(template): 873 msg = ('Length of prepared template stream and data stream are ' 874 'different. Make sure the data does not contain gaps.') 875 raise ValueError(msg) 876 starttime = max(tr.stats.starttime for tr in stream) 877 endtime = min(tr.stats.endtime for tr in stream) 878 starttime_template = min(tr.stats.starttime for tr in template) 879 len_templ = max(tr.stats.endtime - tr.stats.starttime for tr in template) 880 if template_time is None: 881 template_offset = 0 882 else: 883 template_offset = template_time - starttime_template 884 # trim traces 885 trim1 = [trt.stats.starttime - starttime_template for trt in template] 886 trim2 = [trt.stats.endtime - starttime_template - len_templ 887 for trt in template] 888 trim1 = [t - min(trim1) for t in trim1] 889 trim2 = [t - max(trim2) for t in trim2] 890 for i, tr in enumerate(stream): 891 tr = tr.slice(starttime + trim1[i], endtime + trim2[i]) 892 tr.stats.starttime = starttime + template_offset 893 stream.traces[i] = tr 894 return stream, template 895 896 897def _correlate_prepared_stream_template(stream, template, **kwargs): 898 """ 899 Calculate cross-correlation of traces in stream with traces in template. 900 901 Operates on prepared streams. 902 """ 903 for tr, trt in zip(stream, template): 904 tr.data = correlate_template(tr, trt, mode='valid', **kwargs) 905 # make sure xcorrs have the same length, can differ by one sample 906 lens = {len(tr) for tr in stream} 907 if len(lens) > 1: 908 warnings.warn('Samples of traces are slightly misaligned. ' 909 'Use Stream.interpolate if this is not intended.') 910 if max(lens) - min(lens) > 1: 911 msg = 'This should not happen. Please contact the developers.' 912 raise RuntimeError(msg) 913 for tr in stream: 914 tr.data = tr.data[:min(lens)] 915 return stream 916 917 918def correlate_stream_template(stream, template, template_time=None, **kwargs): 919 """ 920 Calculate cross-correlation of traces in stream with traces in template. 921 922 Only matching seed ids are correlated, other traces are silently discarded. 923 The template stream and data stream might have traces of different 924 length and different start times. 925 The data stream must not have gaps and will be sliced as necessary. 926 927 :param stream: Stream with data traces. 928 :param template: Stream with template traces (should be shorter than data). 929 :param template_time: UTCDateTime associated with template event 930 (e.g. origin time, default is the start time of the template stream). 931 The start times of the returned Stream will be shifted by the given 932 template time minus the template start time. 933 :param kwargs: kwargs are passed to 934 :func:`~obspy.signal.cross_correlation.correlate_template` function. 935 936 :return: Stream with cross-correlations. 937 938 .. note:: 939 940 Use :func:`~obspy.signal.cross_correlation.correlation_detector` 941 for detecting events based on their similarity. 942 The returned stream of cross-correlations is suitable for 943 use with :func:`~obspy.signal.trigger.coincidence_trigger`, though. 944 945 .. rubric:: Example 946 947 >>> from obspy import read, UTCDateTime 948 >>> data = read().filter('highpass', freq=5) 949 >>> pick = UTCDateTime('2009-08-24T00:20:07.73') 950 >>> template = data.slice(pick, pick + 10) 951 >>> ccs = correlate_stream_template(data, template) 952 >>> print(ccs) # doctest: +ELLIPSIS 953 3 Trace(s) in Stream: 954 BW.RJOB..EHE | 2009-08-24T00:20:03.000000Z - ... | 100.0 Hz, 2000 samples 955 BW.RJOB..EHN | 2009-08-24T00:20:03.000000Z - ... | 100.0 Hz, 2000 samples 956 BW.RJOB..EHZ | ... - 2009-08-24T00:20:22.990000Z | 100.0 Hz, 2000 samples 957 """ 958 stream, template = _prep_streams_correlate(stream, template, 959 template_time=template_time) 960 return _correlate_prepared_stream_template(stream, template, **kwargs) 961 962 963def _calc_mean(stream): 964 """ 965 Return trace with mean of traces in stream. 966 """ 967 if len(stream) == 0: 968 return stream 969 matrix = np.array([tr.data for tr in stream]) 970 header = dict(sampling_rate=stream[0].stats.sampling_rate, 971 starttime=stream[0].stats.starttime) 972 return Trace(data=np.mean(matrix, axis=0), header=header) 973 974 975def _find_peaks(data, height, holdon_samples, holdoff_samples): 976 """ 977 Peak finding function used for Scipy versions smaller than 1.1. 978 """ 979 cond = data >= height 980 # loop through True values in cond array and guarantee hold time 981 similarity_cond = data[cond] 982 cindices = np.nonzero(cond)[0] 983 detections_index = [] 984 i = 0 985 while True: 986 try: 987 cindex = cindices[i] 988 except IndexError: 989 break 990 # look for maximum inside holdon time 991 j = bisect_left(cindices, cindex + holdon_samples, lo=i) 992 k = i + np.argmax(similarity_cond[i:j]) 993 cindex = cindices[k] 994 detections_index.append(cindex) 995 # wait holdoff time after detection 996 i = bisect_left(cindices, cindex + holdoff_samples, lo=j) 997 return detections_index 998 999 1000def _similarity_detector(similarity, height, distance, 1001 details=False, cross_correlations=None, **kwargs): 1002 """ 1003 Detector based on the similarity of waveforms. 1004 """ 1005 starttime = similarity.stats.starttime 1006 dt = similarity.stats.delta 1007 if distance is not None: 1008 distance = int(round(distance / dt)) 1009 try: 1010 from scipy.signal import find_peaks 1011 except ImportError: 1012 indices = _find_peaks(similarity.data, height, distance, distance) 1013 properties = {} 1014 else: 1015 indices, properties = find_peaks(similarity.data, height, 1016 distance=distance, **kwargs) 1017 detections = [] 1018 for i, index in enumerate(indices): 1019 detection = {'time': starttime + index * dt, 1020 'similarity': similarity.data[index]} 1021 if details and cross_correlations is not None: 1022 detection['cc_values'] = {tr.id: tr.data[index] 1023 for tr in cross_correlations} 1024 if details: 1025 for k, v in properties.items(): 1026 if k != 'peak_heights': 1027 detection[k[:-1] if k.endswith('s') else k] = v[i] 1028 detections.append(detection) 1029 return detections 1030 1031 1032def _insert_amplitude_ratio(detections, stream, template, template_time=None, 1033 template_magnitude=None): 1034 """ 1035 Insert amplitude ratio and magnitude into detections. 1036 """ 1037 stream, template = _prep_streams_correlate(stream, template, 1038 template_time=template_time) 1039 ref_amp = np.mean([np.mean(np.abs(tr.data)) for tr in template]) 1040 for detection in detections: 1041 t = detection['time'] 1042 ratio = np.mean([np.mean(np.abs(tr.slice(t).data[:len(trt)])) 1043 for tr, trt in zip(stream, template)]) / ref_amp 1044 detection['amplitude_ratio'] = ratio 1045 if template_magnitude is not None: 1046 magdiff = 4 / 3 * np.log10(ratio) 1047 detection['magnitude'] = template_magnitude + magdiff 1048 return detections 1049 1050 1051def _get_item(list_, index): 1052 if isinstance(list_, str): 1053 return list_ 1054 try: 1055 return list_[index] 1056 except TypeError: 1057 return list_ 1058 1059 1060def _plot_detections(detections, similarities, stream=None, heights=None, 1061 template_names=None): 1062 """ 1063 Plot detections together with similarity traces and data stream. 1064 """ 1065 import matplotlib.pyplot as plt 1066 from obspy.imaging.util import _set_xaxis_obspy_dates 1067 if stream in (True, None): 1068 stream = [] 1069 akw = dict(xy=(0.02, 0.95), xycoords='axes fraction', va='top') 1070 num1 = len(stream) 1071 num2 = len(similarities) 1072 fig, ax = plt.subplots(num1 + num2, 1, sharex=True) 1073 if num1 + num2 == 1: 1074 ax = [ax] 1075 for detection in detections: 1076 tid = detection.get('template_id', 0) 1077 color = 'C{}'.format((tid + 1) % 10) 1078 for i in list(range(num1)) + [num1 + tid]: 1079 ax[i].axvline(detection['time'].matplotlib_date, color=color) 1080 for i, tr in enumerate(stream): 1081 ax[i].plot(tr.times('matplotlib'), tr.data, 'k') 1082 ax[i].annotate(tr.id, **akw) 1083 for i, tr in enumerate(similarities): 1084 if tr is not None: 1085 ax[num1 + i].plot(tr.times('matplotlib'), tr.data, 'k') 1086 height = _get_item(heights, i) 1087 if isinstance(height, (float, int)): 1088 ax[num1 + i].axhline(height) 1089 template_name = _get_item(template_names, i) 1090 text = ('similarity {}'.format(template_name) if template_name else 1091 'similarity' if num2 == 1 else 1092 'similarity template {}'.format(i)) 1093 ax[num1 + i].annotate(text, **akw) 1094 try: 1095 _set_xaxis_obspy_dates(ax[-1]) 1096 except ValueError: 1097 # work-around for python 2.7, minimum dependencies, see 1098 # https://travis-ci.org/obspy/obspy/jobs/508313177 1099 # can be safely removed later 1100 pass 1101 plt.show() 1102 1103 1104def correlation_detector(stream, templates, heights, distance, 1105 template_times=None, template_magnitudes=None, 1106 template_names=None, 1107 similarity_func=_calc_mean, details=None, 1108 plot=None, **kwargs): 1109 """ 1110 Detector based on the cross-correlation of waveforms. 1111 1112 This detector cross-correlates the stream with each of the template 1113 streams (compare with 1114 :func:`~obspy.signal.cross_correlation.correlate_stream_template`). 1115 A similarity is defined, by default it is the mean of all 1116 cross-correlation functions for each template. 1117 If the similarity exceeds the `height` threshold a detection is triggered. 1118 This peak finding utilizes the SciPy function 1119 :func:`~scipy.signal.find_peaks` with parameters `height` and `distance`. 1120 For a SciPy version smaller than 1.1 it uses a custom function 1121 for peak finding. 1122 1123 :param stream: Stream with data traces. 1124 :param templates: List of streams with template traces. 1125 Each template stream should be shorter than the data stream. 1126 This argument can also be a single template stream. 1127 :param heights: Similarity values to trigger a detection, 1128 one for each template. This argument can also be a single value. 1129 :param distance: The distance in seconds between two detections. 1130 :param template_times: UTCDateTimes associated with template event 1131 (e.g. origin times, 1132 default are the start times of the template streams). 1133 This argument can also be a single value. 1134 :param template_magnitudes: Magnitudes of the template events. 1135 If provided, amplitude ratios between templates and detections will 1136 be calculated and the magnitude of detections will be estimated. 1137 This argument can also be a single value. 1138 This argument can be set to `True`, 1139 then only amplitude ratios will be calculated. 1140 :param template_names: List of template names, the corresponding 1141 template name will be inserted into the detection. 1142 :param similarity_func: By default, the similarity will be calculated by 1143 the mean of cross-correlations. If provided, `similarity_func` will be 1144 called with the stream of cross correlations and the returned trace 1145 will be used as similarity. See the tutorial for an example. 1146 :param details: If set to True detections include detailed information. 1147 :param plot: Plot detections together with the data of the 1148 supplied stream. The default `plot=None` does not plot 1149 anything. `plot=True` plots the similarity traces together 1150 with the detections. If a stream is passed as argument, the traces 1151 in the stream will be plotted together with the similarity traces and 1152 detections. 1153 :param kwargs: Suitable kwargs are passed to 1154 :func:`~obspy.signal.cross_correlation.correlate_template` function. 1155 All other kwargs are passed to :func:`~scipy.signal.find_peaks`. 1156 1157 :return: List of event detections sorted chronologically and 1158 list of similarity traces - one for each template. 1159 Each detection is a dictionary with the following keys: 1160 time, similarity, template_id, 1161 amplitude_ratio, magnitude (if template_magnitudes is provided), 1162 template_name (if template_names is provided), 1163 cross-correlation values, properties returned by find_peaks 1164 (if details are requested) 1165 1166 .. rubric:: Example 1167 1168 >>> from obspy import read, UTCDateTime 1169 >>> data = read().filter('highpass', freq=5) 1170 >>> pick = UTCDateTime('2009-08-24T00:20:07.73') 1171 >>> template = data.slice(pick, pick + 10) 1172 >>> detections, sims = correlation_detector(data, template, 0.5, 10) 1173 >>> print(detections) # doctest: +SKIP 1174 [{'time': UTCDateTime(2009, 8, 24, 0, 20, 7, 730000), 1175 'similarity': 0.99999999999999944, 1176 'template_id': 0}] 1177 1178 A more advanced :ref:`tutorial <correlation-detector-tutorial>` 1179 is available. 1180 """ 1181 if isinstance(templates, Stream): 1182 templates = [templates] 1183 cckeys = ('normalize', 'demean', 'method') 1184 cckwargs = {k: v for k, v in kwargs.items() if k in cckeys} 1185 pfkwargs = {k: v for k, v in kwargs.items() if k not in cckeys} 1186 possible_detections = [] 1187 similarities = [] 1188 for template_id, template in enumerate(templates): 1189 template_time = _get_item(template_times, template_id) 1190 try: 1191 ccs = correlate_stream_template(stream, template, 1192 template_time=template_time, 1193 **cckwargs) 1194 except ValueError as ex: 1195 msg = '{} -> do not use template {}'.format(ex, template_id) 1196 warnings.warn(msg) 1197 similarities.append(None) 1198 continue 1199 similarity = similarity_func(ccs) 1200 height = _get_item(heights, template_id) 1201 detections_template = _similarity_detector( 1202 similarity, height, distance, details=details, 1203 cross_correlations=ccs, **pfkwargs) 1204 for d in detections_template: 1205 template_name = _get_item(template_names, template_id) 1206 if template_name is not None: 1207 d['template_name'] = template_name 1208 d['template_id'] = template_id 1209 if template_magnitudes is True: 1210 template_magnitude = None 1211 else: 1212 template_magnitude = _get_item(template_magnitudes, template_id) 1213 if template_magnitudes is not None: 1214 _insert_amplitude_ratio(detections_template, stream, template, 1215 template_time=template_time, 1216 template_magnitude=template_magnitude) 1217 possible_detections.extend(detections_template) 1218 similarities.append(similarity) 1219 # discard detections with small distance, prefer those with high 1220 # similarity 1221 if len(templates) == 1: 1222 detections = possible_detections 1223 else: 1224 detections = [] 1225 times = [] 1226 for pd in sorted(possible_detections, key=lambda d: -d['similarity']): 1227 if all(abs(pd['time'] - t) > distance for t in times): 1228 times.append(pd['time']) 1229 detections.append(pd) 1230 detections = sorted(detections, key=lambda d: d['time']) 1231 if plot is not None: 1232 _plot_detections(detections, similarities, stream=plot, 1233 heights=heights, template_names=template_names) 1234 return detections, similarities 1235 1236 1237if __name__ == '__main__': 1238 import doctest 1239 doctest.testmod(exclude_empty=True) 1240