Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/trace.py: 77%

1# https://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

4# ---|P------/S----------~Lg----------

6'''

7Basic signal processing for seismic waveforms.

8'''

10import time

11import math

12import copy

13import logging

14import hashlib

15import re

16from collections import defaultdict

18import numpy as num

19from scipy import signal

21from pyrocko import util, orthodrome, pchain, model, signal_ext

22from pyrocko.util import reuse

23from pyrocko.guts import Object, Float, Int, String, List, \

24 StringChoice, Timestamp

25from pyrocko.guts_array import Array

26from pyrocko.model import squirrel_content

28# backward compatibility

29from pyrocko.util import UnavailableDecimation # noqa

30from pyrocko.response import ( # noqa

31 FrequencyResponse, Evalresp, InverseEvalresp, PoleZeroResponse,

32 ButterworthResponse, SampledResponse, IntegrationResponse,

33 DifferentiationResponse, MultiplyResponse)

36guts_prefix = 'pf'

38logger = logging.getLogger('pyrocko.trace')

41g_tapered_coeffs_cache = {}

42g_one_response = FrequencyResponse()

45@squirrel_content

46class Trace(Object):

48 '''

49 Create new trace object.

51 A ``Trace`` object represents a single continuous strip of evenly sampled

52 time series data. It is built from a 1D NumPy array containing the data

53 samples and some attributes describing its beginning and ending time, its

54 sampling rate and four string identifiers (its network, station, location

55 and channel code).

57 :param network: network code

58 :param station: station code

59 :param location: location code

60 :param channel: channel code

61 :param extra: extra code

62 :param tmin: system time of first sample in [s]

63 :param tmax: system time of last sample in [s] (if set to ``None`` it is

64 computed from length of ``ydata``)

65 :param deltat: sampling interval in [s]

66 :param ydata: 1D numpy array with data samples (can be ``None`` when

67 ``tmax`` is not ``None``)

68 :param mtime: optional modification time

70 :param meta: additional meta information (not used, but maintained by the

71 library)

73 The length of the network, station, location and channel codes is not

74 resricted by this software, but data formats like SAC, Mini-SEED or GSE

75 have different limits on the lengths of these codes. The codes set here are

76 silently truncated when the trace is stored

77 '''

79 network = String.T(default='', help='Deployment/network code (1-8)')

80 station = String.T(default='STA', help='Station code (1-5)')

81 location = String.T(default='', help='Location code (0-2)')

82 channel = String.T(default='', help='Channel code (3)')

83 extra = String.T(default='', help='Extra/custom code')

85 tmin = Timestamp.T(default=Timestamp.D('1970-01-01 00:00:00'))

86 tmax = Timestamp.T()

88 deltat = Float.T(default=1.0)

89 ydata = Array.T(optional=True, shape=(None,), serialize_as='base64+meta')

91 mtime = Timestamp.T(optional=True)

93 cached_frequencies = {}

95 def __init__(self, network='', station='STA', location='', channel='',

96 tmin=0., tmax=None, deltat=1., ydata=None, mtime=None,

97 meta=None, extra=''):

99 Object.__init__(self, init_props=False)

100

101 time_float = util.get_time_float()

102

103 if not isinstance(tmin, time_float):

104 tmin = Trace.tmin.regularize_extra(tmin)

105

106 if tmax is not None and not isinstance(tmax, time_float):

107 tmax = Trace.tmax.regularize_extra(tmax)

108

109 if mtime is not None and not isinstance(mtime, time_float):

110 mtime = Trace.mtime.regularize_extra(mtime)

111

112 if ydata is not None and not isinstance(ydata, num.ndarray):

113 ydata = Trace.ydata.regularize_extra(ydata)

114

115 self._growbuffer = None

116

117 tmin = time_float(tmin)

118 if tmax is not None:

119 tmax = time_float(tmax)

120

121 if mtime is None:

122 mtime = time_float(time.time())

123

124 self.network, self.station, self.location, self.channel, \

125 self.extra = [

126 reuse(x) for x in (

127 network, station, location, channel, extra)]

128

129 self.tmin = tmin

130 self.deltat = deltat

131

132 if tmax is None:

133 if ydata is not None:

134 self.tmax = self.tmin + (ydata.size-1)*self.deltat

135 else:

136 raise Exception(

137 'fixme: trace must be created with tmax or ydata')

138 else:

139 if deltat != 0:

140 n = int(round((tmax - self.tmin) / self.deltat)) + 1

141 self.tmax = self.tmin + (n - 1) * self.deltat

142 else:

143 self.tmax = self.tmin

144

145 self.meta = meta

146 self.ydata = ydata

147 self.mtime = mtime

148 self._update_ids()

149 self.file = None

150 self._pchain = None

151

152 def __str__(self):

153 fmt = min(9, max(0, -int(math.floor(math.log10(self.deltat)))))

154 s = 'Trace (%s, %s, %s, %s)\n' % self.nslc_id

155 s += ' timerange: %s - %s\n' % (

156 util.time_to_str(self.tmin, format=fmt),

157 util.time_to_str(self.tmax, format=fmt))

158

159 s += ' delta t: %g\n' % self.deltat

160 if self.meta:

161 for k in sorted(self.meta.keys()):

162 s += ' %s: %s\n' % (k, self.meta[k])

163 return s

164

165 @property

166 def codes(self):

167 from pyrocko.squirrel import model

168 return model.CodesNSLCE(

169 self.network, self.station, self.location, self.channel,

170 self.extra)

171

172 @property

173 def time_span(self):

174 return self.tmin, self.tmax

175

176 @property

177 def summary_entries(self):

178 return (

179 self.__class__.__name__,

180 str(self.codes),

181 self.str_time_span,

182 '%8g' % (1.0/self.deltat if self.deltat != 0.0 else 0.0),

183 '%9i' % self.data_len())

184

185 @property

186 def summary_stats_entries(self):

187 return tuple('%13.7g' % v for v in (

188 self.ydata.min(),

189 self.ydata.max(),

190 self.ydata.mean(),

191 self.ydata.std()))

192

193 @property

194 def summary(self):

195 return util.fmt_summary(self.summary_entries, (10, 20, 55, 8, 9))

196

197 @property

198 def summary_stats(self):

199 return self.summary + ' | ' + util.fmt_summary(

200 self.summary_stats_entries, (12, 12, 12, 12))

201

202 def __getstate__(self):

203 return (self.network, self.station, self.location, self.channel,

204 self.tmin, self.tmax, self.deltat, self.mtime,

205 self.ydata, self.meta, self.extra)

206

207 def __setstate__(self, state):

208 if len(state) == 11:

209 self.network, self.station, self.location, self.channel, \

210 self.tmin, self.tmax, self.deltat, self.mtime, \

211 self.ydata, self.meta, self.extra = state

212

213 elif len(state) == 12:

214 # backward compatibility 0

215 self.network, self.station, self.location, self.channel, \

216 self.tmin, self.tmax, self.deltat, self.mtime, \

217 self.ydata, self.meta, _, self.extra = state

218

219 elif len(state) == 10:

220 # backward compatibility 1

221 self.network, self.station, self.location, self.channel, \

222 self.tmin, self.tmax, self.deltat, self.mtime, \

223 self.ydata, self.meta = state

224

225 self.extra = ''

226

227 else:

228 # backward compatibility 2

229 self.network, self.station, self.location, self.channel, \

230 self.tmin, self.tmax, self.deltat, self.mtime = state

231 self.ydata = None

232 self.meta = None

233 self.extra = ''

234

235 self._growbuffer = None

236 self._update_ids()

237

238 def hash(self, unsafe=False):

239 sha1 = hashlib.sha1()

240 for code in self.nslc_id:

241 sha1.update(code.encode())

242 sha1.update(self.tmin.hex().encode())

243 sha1.update(self.tmax.hex().encode())

244 sha1.update(self.deltat.hex().encode())

245

246 if self.ydata is not None:

247 if not unsafe:

248 sha1.update(memoryview(self.ydata))

249 else:

250 sha1.update(memoryview(self.ydata[:32]))

251 sha1.update(memoryview(self.ydata[-32:]))

252

253 return sha1.hexdigest()

254

255 def name(self):

256 '''

257 Get a short string description.

258 '''

259

260 s = '%s.%s.%s.%s, %s, %s' % (self.nslc_id + (

261 util.time_to_str(self.tmin),

262 util.time_to_str(self.tmax)))

263

264 return s

265

266 def __eq__(self, other):

267 return (

268 isinstance(other, Trace)

269 and self.network == other.network

270 and self.station == other.station

271 and self.location == other.location

272 and self.channel == other.channel

273 and (abs(self.deltat - other.deltat)

274 < (self.deltat + other.deltat)*1e-6)

275 and abs(self.tmin-other.tmin) < self.deltat*0.01

276 and abs(self.tmax-other.tmax) < self.deltat*0.01

277 and num.all(self.ydata == other.ydata))

278

279 def almost_equal(self, other, rtol=1e-5, atol=0.0):

280 return (

281 self.network == other.network

282 and self.station == other.station

283 and self.location == other.location

284 and self.channel == other.channel

285 and (abs(self.deltat - other.deltat)

286 < (self.deltat + other.deltat)*1e-6)

287 and abs(self.tmin-other.tmin) < self.deltat*0.01

288 and abs(self.tmax-other.tmax) < self.deltat*0.01

289 and num.allclose(self.ydata, other.ydata, rtol=rtol, atol=atol))

290

291 def assert_almost_equal(self, other, rtol=1e-5, atol=0.0):

292

293 assert self.network == other.network, \

294 'network codes differ: %s, %s' % (self.network, other.network)

295 assert self.station == other.station, \

296 'station codes differ: %s, %s' % (self.station, other.station)

297 assert self.location == other.location, \

298 'location codes differ: %s, %s' % (self.location, other.location)

299 assert self.channel == other.channel, 'channel codes differ'

300 assert (abs(self.deltat - other.deltat)

301 < (self.deltat + other.deltat)*1e-6), \

302 'sampling intervals differ %g, %g' % (self.deltat, other.delta)

303 assert abs(self.tmin-other.tmin) < self.deltat*0.01, \

304 'start times differ: %s, %s' % (

305 util.time_to_str(self.tmin), util.time_to_str(other.tmin))

306 assert abs(self.tmax-other.tmax) < self.deltat*0.01, \

307 'end times differ: %s, %s' % (

308 util.time_to_str(self.tmax), util.time_to_str(other.tmax))

309

310 assert num.allclose(self.ydata, other.ydata, rtol=rtol, atol=atol), \

311 'trace values differ'

312

313 def __hash__(self):

314 return id(self)

315

316 def __call__(self, t, clip=False, snap=round):

317 it = int(snap((t-self.tmin)/self.deltat))

318 if clip:

319 it = max(0, min(it, self.ydata.size-1))

320 else:

321 if it < 0 or self.ydata.size <= it:

322 raise IndexError()

323

324 return self.tmin+it*self.deltat, self.ydata[it]

325

326 def interpolate(self, t, clip=False):

327 '''

328 Value of trace between supporting points through linear interpolation.

329

330 :param t: time instant

331 :param clip: whether to clip indices to trace ends

332 '''

333

334 t0, y0 = self(t, clip=clip, snap=math.floor)

335 t1, y1 = self(t, clip=clip, snap=math.ceil)

336 if t0 == t1:

337 return y0

338 else:

339 return y0+(t-t0)/(t1-t0)*(y1-y0)

340

341 def index_clip(self, i):

342 '''

343 Clip index to valid range.

344 '''

345

346 return min(max(0, i), self.ydata.size)

347

348 def add(self, other, interpolate=True, left=0., right=0.):

349 '''

350 Add values of other trace (self += other).

351

352 Add values of ``other`` trace to the values of ``self``, where it

353 intersects with ``other``. This method does not change the extent of

354 ``self``. If ``interpolate`` is ``True`` (the default), the values of

355 ``other`` to be added are interpolated at sampling instants of

356 ``self``. Linear interpolation is performed. In this case the sampling

357 rate of ``other`` must be equal to or lower than that of ``self``. If

358 ``interpolate`` is ``False``, the sampling rates of the two traces must

359 match.

360 '''

361

362 if interpolate:

363 assert self.deltat <= other.deltat \

364 or same_sampling_rate(self, other)

365

366 other_xdata = other.get_xdata()

367 xdata = self.get_xdata()

368 self.ydata += num.interp(

369 xdata, other_xdata, other.ydata, left=left, right=left)

370 else:

371 assert self.deltat == other.deltat

372 ioff = int(round((other.tmin-self.tmin)/self.deltat))

373 ibeg = max(0, ioff)

374 iend = min(self.data_len(), ioff+other.data_len())

375 self.ydata[ibeg:iend] += other.ydata[ibeg-ioff:iend-ioff]

376

377 def mult(self, other, interpolate=True):

378 '''

379 Muliply with values of other trace ``(self *= other)``.

380

381 Multiply values of ``other`` trace to the values of ``self``, where it

382 intersects with ``other``. This method does not change the extent of

383 ``self``. If ``interpolate`` is ``True`` (the default), the values of

384 ``other`` to be multiplied are interpolated at sampling instants of

385 ``self``. Linear interpolation is performed. In this case the sampling

386 rate of ``other`` must be equal to or lower than that of ``self``. If

387 ``interpolate`` is ``False``, the sampling rates of the two traces must

388 match.

389 '''

390

391 if interpolate:

392 assert self.deltat <= other.deltat or \

393 same_sampling_rate(self, other)

394

395 other_xdata = other.get_xdata()

396 xdata = self.get_xdata()

397 self.ydata *= num.interp(

398 xdata, other_xdata, other.ydata, left=0., right=0.)

399 else:

400 assert self.deltat == other.deltat

401 ibeg1 = int(round((other.tmin-self.tmin)/self.deltat))

402 ibeg2 = int(round((self.tmin-other.tmin)/self.deltat))

403 iend1 = int(round((other.tmax-self.tmin)/self.deltat))+1

404 iend2 = int(round((self.tmax-other.tmin)/self.deltat))+1

405

406 ibeg1 = self.index_clip(ibeg1)

407 iend1 = self.index_clip(iend1)

408 ibeg2 = self.index_clip(ibeg2)

409 iend2 = self.index_clip(iend2)

410

411 self.ydata[ibeg1:iend1] *= other.ydata[ibeg2:iend2]

412

413 def max(self):

414 '''

415 Get time and value of data maximum.

416 '''

417

418 i = num.argmax(self.ydata)

419 return self.tmin + i*self.deltat, self.ydata[i]

420

421 def min(self):

422 '''

423 Get time and value of data minimum.

424 '''

425

426 i = num.argmin(self.ydata)

427 return self.tmin + i*self.deltat, self.ydata[i]

428

429 def absmax(self):

430 '''

431 Get time and value of maximum of the absolute of data.

432 '''

433

434 tmi, mi = self.min()

435 tma, ma = self.max()

436 if abs(mi) > abs(ma):

437 return tmi, abs(mi)

438 else:

439 return tma, abs(ma)

440

441 def set_codes(

442 self, network=None, station=None, location=None, channel=None,

443 extra=None):

444

445 '''

446 Set network, station, location, and channel codes.

447 '''

448

449 if network is not None:

450 self.network = network

451 if station is not None:

452 self.station = station

453 if location is not None:

454 self.location = location

455 if channel is not None:

456 self.channel = channel

457 if extra is not None:

458 self.extra = extra

459

460 self._update_ids()

461

462 def set_network(self, network):

463 self.network = network

464 self._update_ids()

465

466 def set_station(self, station):

467 self.station = station

468 self._update_ids()

469

470 def set_location(self, location):

471 self.location = location

472 self._update_ids()

473

474 def set_channel(self, channel):

475 self.channel = channel

476 self._update_ids()

477

478 def overlaps(self, tmin, tmax):

479 '''

480 Check if trace has overlap with a given time span.

481 '''

482 return not (tmax < self.tmin or self.tmax < tmin)

483

484 def is_relevant(self, tmin, tmax, selector=None):

485 '''

486 Check if trace has overlap with a given time span and matches a

487 condition callback. (internal use)

488 '''

489

490 return not (tmax <= self.tmin or self.tmax < tmin) \

491 and (selector is None or selector(self))

492

493 def _update_ids(self):

494 '''

495 Update dependent ids.

496 '''

497

498 self.full_id = (

499 self.network, self.station, self.location, self.channel, self.tmin)

500 self.nslc_id = reuse(

501 (self.network, self.station, self.location, self.channel))

502

503 def prune_from_reuse_cache(self):

504 util.deuse(self.nslc_id)

505 util.deuse(self.network)

506 util.deuse(self.station)

507 util.deuse(self.location)

508 util.deuse(self.channel)

509

510 def set_mtime(self, mtime):

511 '''

512 Set modification time of the trace.

513 '''

514

515 self.mtime = mtime

516

517 def get_xdata(self):

518 '''

519 Create array for time axis.

520 '''

521

522 if self.ydata is None:

523 raise NoData()

524

525 return self.tmin \

526 + num.arange(len(self.ydata), dtype=num.float64) * self.deltat

527

528 def get_ydata(self):

529 '''

530 Get data array.

531 '''

532

533 if self.ydata is None:

534 raise NoData()

535

536 return self.ydata

537

538 def set_ydata(self, new_ydata):

539 '''

540 Replace data array.

541 '''

542

543 self.drop_growbuffer()

544 self.ydata = new_ydata

545 self.tmax = self.tmin+(len(self.ydata)-1)*self.deltat

546

547 def data_len(self):

548 if self.ydata is not None:

549 return self.ydata.size

550 else:

551 return int(round((self.tmax-self.tmin)/self.deltat)) + 1

552

553 def drop_data(self):

554 '''

555 Forget data, make dataless trace.

556 '''

557

558 self.drop_growbuffer()

559 self.ydata = None

560

561 def drop_growbuffer(self):

562 '''

563 Detach the traces grow buffer.

564 '''

565

566 self._growbuffer = None

567 self._pchain = None

568

569 def copy(self, data=True):

570 '''

571 Make a deep copy of the trace.

572 '''

573

574 tracecopy = copy.copy(self)

575 tracecopy.drop_growbuffer()

576 if data:

577 tracecopy.ydata = self.ydata.copy()

578 tracecopy.meta = copy.deepcopy(self.meta)

579 return tracecopy

580

581 def crop_zeros(self):

582 '''

583 Remove any zeros at beginning and end.

584 '''

585

586 indices = num.where(self.ydata != 0.0)[0]

587 if indices.size == 0:

588 raise NoData()

589

590 ibeg = indices[0]

591 iend = indices[-1]+1

592 if ibeg == 0 and iend == self.ydata.size-1:

593 return

594

595 self.drop_growbuffer()

596 self.ydata = self.ydata[ibeg:iend].copy()

597 self.tmin = self.tmin+ibeg*self.deltat

598 self.tmax = self.tmin+(len(self.ydata)-1)*self.deltat

599 self._update_ids()

600

601 def append(self, data):

602 '''

603 Append data to the end of the trace.

604

605 To make this method efficient when successively very few or even single

606 samples are appended, a larger grow buffer is allocated upon first

607 invocation. The traces data is then changed to be a view into the

608 currently filled portion of the grow buffer array.

609 '''

610

611 assert self.ydata.dtype == data.dtype

612 newlen = data.size + self.ydata.size

613 if self._growbuffer is None or self._growbuffer.size < newlen:

614 self._growbuffer = num.empty(newlen*2, dtype=self.ydata.dtype)

615 self._growbuffer[:self.ydata.size] = self.ydata

616 self._growbuffer[self.ydata.size:newlen] = data

617 self.ydata = self._growbuffer[:newlen]

618 self.tmax = self.tmin + (newlen-1)*self.deltat

619

620 def chop(

621 self, tmin, tmax, inplace=True, include_last=False,

622 snap=(round, round), want_incomplete=True):

623

624 '''

625 Cut the trace to given time span.

626

627 If the ``inplace`` argument is True (the default) the trace is cut in

628 place, otherwise a new trace with the cut part is returned. By

629 default, the indices where to start and end the trace data array are

630 determined by rounding of ``tmin`` and ``tmax`` to sampling instances

631 using Python's :py:func:`round` function. This behaviour can be changed

632 with the ``snap`` argument, which takes a tuple of two functions (one

633 for the lower and one for the upper end) to be used instead of

634 :py:func:`round`. The last sample is by default not included unless

635 ``include_last`` is set to True. If the given time span exceeds the

636 available time span of the trace, the available part is returned,

637 unless ``want_incomplete`` is set to False - in that case, a

638 :py:exc:`NoData` exception is raised. This exception is always raised,

639 when the requested time span does dot overlap with the trace's time

640 span.

641 '''

642

643 if self.deltat == 0:

644 assert self.tmin == self.tmax

645 if self.tmin < tmin or self.tmin > tmax:

646 raise NoData()

647 else:

648 if not inplace:

649 return self.copy()

650 else:

651 return

652

653 if want_incomplete:

654 if tmax <= self.tmin-self.deltat or self.tmax+self.deltat < tmin:

655 raise NoData()

656 else:

657 if tmin < self.tmin or self.tmax < tmax:

658 raise NoData()

659

660 ibeg = max(0, t2ind(tmin-self.tmin, self.deltat, snap[0]))

661 iplus = 0

662 if include_last:

663 iplus = 1

664

665 iend = min(

666 self.data_len(),

667 t2ind(tmax-self.tmin, self.deltat, snap[1])+iplus)

668

669 if ibeg >= iend:

670 raise NoData()

671

672 obj = self

673 if not inplace:

674 obj = self.copy(data=False)

675

676 self.drop_growbuffer()

677 if self.ydata is not None:

678 obj.ydata = self.ydata[ibeg:iend].copy()

679 else:

680 obj.ydata = None

681

682 obj.tmin = obj.tmin+ibeg*obj.deltat

683 obj.tmax = obj.tmin+((iend-ibeg)-1)*obj.deltat

684

685 obj._update_ids()

686

687 return obj

688

689 def downsample(

690 self, ndecimate, snap=False, demean=False, ftype='fir-remez',

691 cut=False):

692

693 '''

694 Downsample (decimate) trace by a given integer factor.

695

696 Antialiasing filter details:

697

698 * ``iir``: A Chebyshev type I digital filter of order 8 with maximum

699 ripple of 0.05 dB in the passband.

700 * ``fir``: A FIR filter using a Hamming window and 31 taps and a

701 well-behaved phase response.

702 * ``fir-remez``: A FIR filter based on the Remez exchange algorithm

703 with ``45*ndecimate`` taps and a cutoff at 75% Nyquist frequency.

704

705 Comparison of the digital filters:

706

707 .. figure :: ../../static/downsampling-filter-comparison.png

708 :width: 60%

709 :alt: Comparison of the downsampling filters.

710

711 See also :py:meth:`Trace.downsample_to`.

712

713 :param ndecimate:

714 Decimation factor, avoid values larger than 8.

715 :type ndecimate:

716 int

717

718 :param snap:

719 Whether to put the new sampling instants closest to multiples of

720 the sampling rate (according to absolute time).

721 :type snap:

722 bool

723

724 :param demean:

725 Whether to demean the signal before filtering.

726 :type demean:

727 bool

728

729 :param ftype:

730 Which FIR filter to use, choose from ``'iir'``, ``'fir'``,

731 ``'fir-remez'``. Default is ``'fir-remez'``.

732

733 :param cut:

734 Whether to cut off samples in the beginning of the trace which

735 are polluted by artifacts of the anti-aliasing filter.

736 :type cut:

737 bool

738 '''

739 newdeltat = self.deltat*ndecimate

740 b, a, n = util.decimate_coeffs(ndecimate, None, ftype)

741 if snap:

742 ilag = int(round((math.ceil(

743 (self.tmin+(n//2 if cut else 0)*self.deltat) /

744 newdeltat) * newdeltat - self.tmin) / self.deltat))

745 else:

746 ilag = (n//2 if cut else 0)

747

748 data = self.ydata.astype(num.float64)

749 if data.size != 0:

750 if demean:

751 data -= num.mean(data)

752 y = signal.lfilter(b, a, data)

753 self.ydata = y[n//2+ilag::ndecimate].copy()

754 else:

755 self.ydata = data

756

757 self.tmin += ilag * self.deltat

758 self.deltat = reuse(self.deltat*ndecimate)

759 self.tmax = self.tmin+(len(self.ydata)-1)*self.deltat

760 self._update_ids()

761

762 def downsample_to(

763 self, deltat, snap=False, allow_upsample_max=1, demean=False,

764 ftype='fir-remez', cut=False):

765

766 '''

767 Downsample to given sampling rate.

768

769 Tries to downsample the trace to a target sampling interval of

770 ``deltat``. This runs :py:meth:`downsample` one or several times. If

771 ``allow_upsample_max`` is set to a value larger than 1, intermediate

772 upsampling steps are allowed, in order to increase the number of

773 possible downsampling ratios.

774

775 If the requested ratio is not supported, an exception of type

776 :py:exc:`pyrocko.util.UnavailableDecimation` is raised.

777

778 The downsampled trace will be shorter than the input trace because the

779 anti-aliasing filter introduces edge effects. If `cut` is selected,

780 additionally, polluted samples in the beginning of the trace are

781 removed. The approximate amount of cutoff which will be produced by a

782 given downsampling configuration can be estimated with

783 :py:func:`downsample_tpad`.

784

785 See also: :meth:`Trace.downsample`.

786

787 :param deltat:

788 Desired sampling interval in [s].

789 :type deltat:

790 float

791

792 :param allow_upsample_max:

793 Maximum allowed upsampling factor of the signal to achieve the

794 desired sampling interval. Default is ``1``.

795 :type allow_upsample_max:

796 int

797

798 :param snap:

799 Whether to put the new sampling instants closest to multiples of

800 the sampling rate (according to absolute time).

801 :type snap:

802 bool

803

804 :param demean:

805 Whether to demean the signal before filtering.

806 :type demean:

807 bool

808

809 :param ftype:

810 Which FIR filter to use, choose from ``'iir'``, ``'fir'``,

811 ``'fir-remez'``. Default is ``'fir-remez'``.

812

813 :param cut:

814 Whether to cut off samples in the beginning of the trace which

815 are polluted by artifacts of the anti-aliasing filter.

816 :type demean:

817 bool

818 '''

819

820 upsratio, deci_seq = _configure_downsampling(

821 self.deltat, deltat, allow_upsample_max)

822

823 if demean:

824 self.drop_growbuffer()

825 self.ydata = self.ydata.astype(num.float64)

826 self.ydata -= num.mean(self.ydata)

827

828 if upsratio > 1:

829 self.drop_growbuffer()

830 ydata = self.ydata

831 self.ydata = num.zeros(

832 ydata.size*upsratio-(upsratio-1), ydata.dtype)

833 self.ydata[::upsratio] = ydata

834 for i in range(1, upsratio):

835 self.ydata[i::upsratio] = \

836 float(i)/upsratio * ydata[:-1] \

837 + float(upsratio-i)/upsratio * ydata[1:]

838 self.deltat = self.deltat/upsratio

839

840 for i, ndecimate in enumerate(deci_seq):

841 self.downsample(

842 ndecimate, snap=snap, demean=False, ftype=ftype, cut=cut)

843

844 def resample(self, deltat):

845 '''

846 Resample to given sampling rate ``deltat``.

847

848 Resampling is performed in the frequency domain.

849 '''

850

851 ndata = self.ydata.size

852 ntrans = nextpow2(ndata)

853 fntrans2 = ntrans * self.deltat/deltat

854 ntrans2 = int(round(fntrans2))

855 deltat2 = self.deltat * float(ntrans)/float(ntrans2)

856 ndata2 = int(round(ndata*self.deltat/deltat2))

857 if abs(fntrans2 - ntrans2) > 1e-7:

858 logger.warning(

859 'resample: requested deltat %g could not be matched exactly: '

860 '%g' % (deltat, deltat2))

861

862 data = self.ydata

863 data_pad = num.zeros(ntrans, dtype=float)

864 data_pad[:ndata] = data

865 fdata = num.fft.rfft(data_pad)

866 fdata2 = num.zeros((ntrans2+1)//2, dtype=fdata.dtype)

867 n = min(fdata.size, fdata2.size)

868 fdata2[:n] = fdata[:n]

869 data2 = num.fft.irfft(fdata2)

870 data2 = data2[:ndata2]

871 data2 *= float(ntrans2) / float(ntrans)

872 self.deltat = deltat2

873 self.set_ydata(data2)

874

875 def resample_simple(self, deltat):

876 tyear = 3600*24*365.

877

878 if deltat == self.deltat:

879 return

880

881 if abs(self.deltat - deltat) * tyear/deltat < deltat:

882 logger.warning(

883 'resample_simple: less than one sample would have to be '

884 'inserted/deleted per year. Doing nothing.')

885 return

886

887 ninterval = int(round(deltat / (self.deltat - deltat)))

888 if abs(ninterval) < 20:

889 logger.error(

890 'resample_simple: sample insertion/deletion interval less '

891 'than 20. results would be erroneous.')

892 raise ResamplingFailed()

893

894 delete = False

895 if ninterval < 0:

896 ninterval = - ninterval

897 delete = True

898

899 tyearbegin = util.year_start(self.tmin)

900

901 nmin = int(round((self.tmin - tyearbegin)/deltat))

902

903 ibegin = (((nmin-1)//ninterval)+1) * ninterval - nmin

904 nindices = (len(self.ydata) - ibegin - 1) / ninterval + 1

905 if nindices > 0:

906 indices = ibegin + num.arange(nindices) * ninterval

907 data_split = num.split(self.ydata, indices)

908 data = []

909 for ln, h in zip(data_split[:-1], data_split[1:]):

910 if delete:

911 ln = ln[:-1]

912

913 data.append(ln)

914 if not delete:

915 if ln.size == 0:

916 v = h[0]

917 else:

918 v = 0.5*(ln[-1] + h[0])

919 data.append(num.array([v], dtype=ln.dtype))

920

921 data.append(data_split[-1])

922

923 ydata_new = num.concatenate(data)

924

925 self.tmin = tyearbegin + nmin * deltat

926 self.deltat = deltat

927 self.set_ydata(ydata_new)

928

929 def stretch(self, tmin_new, tmax_new):

930 '''

931 Stretch signal while preserving sample rate using sinc interpolation.

932

933 :param tmin_new: new time of first sample

934 :param tmax_new: new time of last sample

935

936 This method can be used to correct for a small linear time drift or to

937 introduce sub-sample time shifts. The amount of stretching is limited

938 to 10% by the implementation and is expected to be much smaller than

939 that by the approximations used.

940 '''

941

942 i_control = num.array([0, self.ydata.size-1], dtype=num.int64)

943 t_control = num.array([tmin_new, tmax_new], dtype=float)

944

945 r = (tmax_new - tmin_new) / self.deltat + 1.0

946 n_new = int(round(r))

947 if abs(n_new - r) > 0.001:

948 n_new = int(math.floor(r))

949

950 assert n_new >= 2

951

952 tmax_new = tmin_new + (n_new-1) * self.deltat

953

954 ydata_new = num.empty(n_new, dtype=float)

955 signal_ext.antidrift(i_control, t_control,

956 self.ydata.astype(float),

957 tmin_new, self.deltat, ydata_new)

958

959 self.tmin = tmin_new

960 self.set_ydata(ydata_new)

961 self._update_ids()

962

963 def nyquist_check(self, frequency, intro='Corner frequency', warn=True,

964 raise_exception=False):

965

966 '''

967 Check if a given frequency is above the Nyquist frequency of the trace.

968

969 :param intro: string used to introduce the warning/error message

970 :param warn: whether to emit a warning

971 :param raise_exception: whether to raise an :py:exc:`AboveNyquist`

972 exception.

973 '''

974

975 if frequency >= 0.5/self.deltat:

976 message = '%s (%g Hz) is equal to or higher than nyquist ' \

977 'frequency (%g Hz). (Trace %s)' \

978 % (intro, frequency, 0.5/self.deltat, self.name())

979 if warn:

980 logger.warning(message)

981 if raise_exception:

982 raise AboveNyquist(message)

983

984 def lowpass(self, order, corner, nyquist_warn=True,

985 nyquist_exception=False, demean=True):

986

987 '''

988 Apply Butterworth lowpass to the trace.

989

990 :param order: order of the filter

991 :param corner: corner frequency of the filter

992

993 Mean is removed before filtering.

994 '''

995

996 self.nyquist_check(

997 corner, 'Corner frequency of lowpass', nyquist_warn,

998 nyquist_exception)

999

1000 (b, a) = _get_cached_filter_coeffs(

1001 order, [corner*2.0*self.deltat], btype='low')

1002

1003 if len(a) != order+1 or len(b) != order+1:

1004 logger.warning(

1005 'Erroneous filter coefficients returned by '

1006 'scipy.signal.butter(). You may need to downsample the '

1007 'signal before filtering.')

1008

1009 data = self.ydata.astype(num.float64)

1010 if demean:

1011 data -= num.mean(data)

1012 self.drop_growbuffer()

1013 self.ydata = signal.lfilter(b, a, data)

1014

1015 def highpass(self, order, corner, nyquist_warn=True,

1016 nyquist_exception=False, demean=True):

1017

1018 '''

1019 Apply butterworth highpass to the trace.

1020

1021 :param order: order of the filter

1022 :param corner: corner frequency of the filter

1023

1024 Mean is removed before filtering.

1025 '''

1026

1027 self.nyquist_check(

1028 corner, 'Corner frequency of highpass', nyquist_warn,

1029 nyquist_exception)

1030

1031 (b, a) = _get_cached_filter_coeffs(

1032 order, [corner*2.0*self.deltat], btype='high')

1033

1034 data = self.ydata.astype(num.float64)

1035 if len(a) != order+1 or len(b) != order+1:

1036 logger.warning(

1037 'Erroneous filter coefficients returned by '

1038 'scipy.signal.butter(). You may need to downsample the '

1039 'signal before filtering.')

1040 if demean:

1041 data -= num.mean(data)

1042 self.drop_growbuffer()

1043 self.ydata = signal.lfilter(b, a, data)

1044

1045 def bandpass(self, order, corner_hp, corner_lp, demean=True):

1046 '''

1047 Apply butterworth bandpass to the trace.

1048

1049 :param order: order of the filter

1050 :param corner_hp: lower corner frequency of the filter

1051 :param corner_lp: upper corner frequency of the filter

1052

1053 Mean is removed before filtering.

1054 '''

1055

1056 self.nyquist_check(corner_hp, 'Lower corner frequency of bandpass')

1057 self.nyquist_check(corner_lp, 'Higher corner frequency of bandpass')

1058 (b, a) = _get_cached_filter_coeffs(

1059 order,

1060 [corner*2.0*self.deltat for corner in (corner_hp, corner_lp)],

1061 btype='band')

1062 data = self.ydata.astype(num.float64)

1063 if demean:

1064 data -= num.mean(data)

1065 self.drop_growbuffer()

1066 self.ydata = signal.lfilter(b, a, data)

1067

1068 def bandstop(self, order, corner_hp, corner_lp, demean=True):

1069 '''

1070 Apply bandstop (attenuates frequencies in band) to the trace.

1071

1072 :param order: order of the filter

1073 :param corner_hp: lower corner frequency of the filter

1074 :param corner_lp: upper corner frequency of the filter

1075

1076 Mean is removed before filtering.

1077 '''

1078

1079 self.nyquist_check(corner_hp, 'Lower corner frequency of bandstop')

1080 self.nyquist_check(corner_lp, 'Higher corner frequency of bandstop')

1081 (b, a) = _get_cached_filter_coeffs(

1082 order,

1083 [corner*2.0*self.deltat for corner in (corner_hp, corner_lp)],

1084 btype='bandstop')

1085 data = self.ydata.astype(num.float64)

1086 if demean:

1087 data -= num.mean(data)

1088 self.drop_growbuffer()

1089 self.ydata = signal.lfilter(b, a, data)

1090

1091 def envelope(self, inplace=True):

1092 '''

1093 Calculate the envelope of the trace.

1094

1095 :param inplace: calculate envelope in place

1096

1097 The calculation follows:

1098

1099 .. math::

1100

1101 Y' = \\sqrt{Y^2+H(Y)^2}

1102

1103 where H is the Hilbert-Transform of the signal Y.

1104 '''

1105

1106 y = self.ydata.astype(float)

1107 env = num.abs(hilbert(y))

1108 if inplace:

1109 self.drop_growbuffer()

1110 self.ydata = env

1111 else:

1112 tr = self.copy(data=False)

1113 tr.ydata = env

1114 return tr

1115

1116 def taper(self, taperer, inplace=True, chop=False):

1117 '''

1118 Apply a :py:class:`Taper` to the trace.

1119

1120 :param taperer: instance of :py:class:`Taper` subclass

1121 :param inplace: apply taper inplace

1122 :param chop: if ``True``: exclude tapered parts from the resulting

1123 trace

1124 '''

1125

1126 if not inplace:

1127 tr = self.copy()

1128 else:

1129 tr = self

1130

1131 if chop:

1132 i, n = taperer.span(tr.ydata, tr.tmin, tr.deltat)

1133 tr.shift(i*tr.deltat)

1134 tr.set_ydata(tr.ydata[i:i+n])

1135

1136 taperer(tr.ydata, tr.tmin, tr.deltat)

1137

1138 if not inplace:

1139 return tr

1140

1141 def whiten(self, order=6):

1142 '''

1143 Whiten signal in time domain using autoregression and recursive filter.

1144

1145 :param order: order of the autoregression process

1146 '''

1147

1148 b, a = self.whitening_coefficients(order)

1149 self.drop_growbuffer()

1150 self.ydata = signal.lfilter(b, a, self.ydata)

1151

1152 def whitening_coefficients(self, order=6):

1153 ar = yulewalker(self.ydata, order)

1154 b, a = [1.] + ar.tolist(), [1.]

1155 return b, a

1156

1157 def ampspec_whiten(

1158 self,

1159 width,

1160 td_taper='auto',

1161 fd_taper='auto',

1162 pad_to_pow2=True,

1163 demean=True):

1164

1165 '''

1166 Whiten signal via frequency domain using moving average on amplitude

1167 spectra.

1168

1169 :param width: width of smoothing kernel [Hz]

1170 :param td_taper: time domain taper, object of type :py:class:`Taper` or

1171 ``None`` or ``'auto'``.

1172 :param fd_taper: frequency domain taper, object of type

1173 :py:class:`Taper` or ``None`` or ``'auto'``.

1174 :param pad_to_pow2: whether to pad the signal with zeros up to a length

1175 of 2^n

1176 :param demean: whether to demean the signal before tapering

1177

1178 The signal is first demeaned and then tapered using ``td_taper``. Then,

1179 the spectrum is calculated and inversely weighted with a smoothed

1180 version of its amplitude spectrum. A moving average is used for the

1181 smoothing. The smoothed spectrum is then tapered using ``fd_taper``.

1182 Finally, the smoothed and tapered spectrum is back-transformed into the

1183 time domain.

1184

1185 If ``td_taper`` is set to ``'auto'``, ``CosFader(1.0/width)`` is used.

1186 If ``fd_taper`` is set to ``'auto'``, ``CosFader(width)`` is used.

1187 '''

1188

1189 ndata = self.data_len()

1190

1191 if pad_to_pow2:

1192 ntrans = nextpow2(ndata)

1193 else:

1194 ntrans = ndata

1195

1196 df = 1./(ntrans*self.deltat)

1197 nw = int(round(width/df))

1198 if ndata//2+1 <= nw:

1199 raise TraceTooShort(

1200 'Samples in trace: %s, samples needed: %s' % (ndata, nw))

1201

1202 if td_taper == 'auto':

1203 td_taper = CosFader(1./width)

1204

1205 if fd_taper == 'auto':

1206 fd_taper = CosFader(width)

1207

1208 if td_taper:

1209 self.taper(td_taper)

1210

1211 ydata = self.get_ydata().astype(float)

1212 if demean:

1213 ydata -= ydata.mean()

1214

1215 spec = num.fft.rfft(ydata, ntrans)

1216

1217 amp = num.abs(spec)

1218 nspec = amp.size

1219 csamp = num.cumsum(amp)

1220 amp_smoothed = num.empty(nspec, dtype=csamp.dtype)

1221 n1, n2 = nw//2, nw//2 + nspec - nw

1222 amp_smoothed[n1:n2] = (csamp[nw:] - csamp[:-nw]) / nw

1223 amp_smoothed[:n1] = amp_smoothed[n1]

1224 amp_smoothed[n2:] = amp_smoothed[n2-1]

1225

1226 denom = amp_smoothed * amp

1227 numer = amp

1228 eps = num.mean(denom) * 1e-9

1229 if eps == 0.0:

1230 eps = 1e-9

1231

1232 numer += eps

1233 denom += eps

1234 spec *= numer/denom

1235

1236 if fd_taper:

1237 fd_taper(spec, 0., df)

1238

1239 ydata = num.fft.irfft(spec)

1240 self.set_ydata(ydata[:ndata])

1241

1242 def _get_cached_freqs(self, nf, deltaf):

1243 ck = (nf, deltaf)

1244 if ck not in Trace.cached_frequencies:

1245 Trace.cached_frequencies[ck] = deltaf * num.arange(

1246 nf, dtype=float)

1247

1248 return Trace.cached_frequencies[ck]

1249

1250 def bandpass_fft(self, corner_hp, corner_lp):

1251 '''

1252 Apply boxcar bandbpass to trace (in spectral domain).

1253 '''

1254

1255 n = len(self.ydata)

1256 n2 = nextpow2(n)

1257 data = num.zeros(n2, dtype=num.float64)

1258 data[:n] = self.ydata

1259 fdata = num.fft.rfft(data)

1260 freqs = self._get_cached_freqs(len(fdata), 1./(self.deltat*n2))

1261 fdata[0] = 0.0

1262 fdata *= num.logical_and(corner_hp < freqs, freqs < corner_lp)

1263 data = num.fft.irfft(fdata)

1264 self.drop_growbuffer()

1265 self.ydata = data[:n]

1266

1267 def shift(self, tshift):

1268 '''

1269 Time shift the trace.

1270 '''

1271

1272 self.tmin += tshift

1273 self.tmax += tshift

1274 self._update_ids()

1275

1276 def snap(self, inplace=True, interpolate=False):

1277 '''

1278 Shift trace samples to nearest even multiples of the sampling rate.

1279

1280 :param inplace: (boolean) snap traces inplace

1281

1282 If ``inplace`` is ``False`` and the difference of tmin and tmax of

1283 both, the snapped and the original trace is smaller than 0.01 x deltat,

1284 :py:func:`snap` returns the unsnapped instance of the original trace.

1285 '''

1286

1287 tmin = round(self.tmin/self.deltat)*self.deltat

1288 tmax = tmin + (self.ydata.size-1)*self.deltat

1289

1290 if inplace:

1291 xself = self

1292 else:

1293 if abs(self.tmin - tmin) < 1e-2*self.deltat and \

1294 abs(self.tmax - tmax) < 1e-2*self.deltat:

1295 return self

1296

1297 xself = self.copy()

1298

1299 n = xself.data_len()

1300 if interpolate and n > 2:

1301 ydata_new = num.empty(n, dtype=float)

1302 i_control = num.array([0, n-1], dtype=num.int64)

1303 tref = tmin

1304 t_control = num.array(

1305 [float(xself.tmin-tref), float(xself.tmax-tref)],

1306 dtype=float)

1307

1308 signal_ext.antidrift(i_control, t_control,

1309 xself.ydata.astype(float),

1310 float(tmin-tref), xself.deltat, ydata_new)

1311

1312 xself.ydata = ydata_new

1313

1314 xself.tmin = tmin

1315 xself.tmax = tmax

1316 xself._update_ids()

1317

1318 return xself

1319

1320 def fix_deltat_rounding_errors(self):

1321 '''

1322 Try to undo sampling rate rounding errors.

1323

1324 See :py:func:`fix_deltat_rounding_errors`.

1325 '''

1326

1327 self.deltat = fix_deltat_rounding_errors(self.deltat)

1328 self.tmax = self.tmin + (self.data_len() - 1) * self.deltat

1329

1330 def sta_lta_centered(self, tshort, tlong, quad=True, scalingmethod=1):

1331 '''

1332 Run special STA/LTA filter where the short time window is centered on

1333 the long time window.

1334

1335 :param tshort: length of short time window in [s]

1336 :param tlong: length of long time window in [s]

1337 :param quad: whether to square the data prior to applying the STA/LTA

1338 filter

1339 :param scalingmethod: integer key to select how output values are

1340 scaled / normalized (``1``, ``2``, or ``3``)

1341

1342 ============= ====================================== ===========

1343 Scalingmethod Implementation Range

1344 ============= ====================================== ===========

1345 ``1`` As/Al* Ts/Tl [0,1]

1346 ``2`` (As/Al - 1) / (Tl/Ts - 1) [-Ts/Tl,1]

1347 ``3`` Like ``2`` but clipping range at zero [0,1]

1348 ============= ====================================== ===========

1349

1350 '''

1351

1352 nshort = int(round(tshort/self.deltat))

1353 nlong = int(round(tlong/self.deltat))

1354

1355 assert nshort < nlong

1356 if nlong > len(self.ydata):

1357 raise TraceTooShort(

1358 'Samples in trace: %s, samples needed: %s'

1359 % (len(self.ydata), nlong))

1360

1361 if quad:

1362 sqrdata = self.ydata**2

1363 else:

1364 sqrdata = self.ydata

1365

1366 mavg_short = moving_avg(sqrdata, nshort)

1367 mavg_long = moving_avg(sqrdata, nlong)

1368

1369 self.drop_growbuffer()

1370

1371 if scalingmethod not in (1, 2, 3):

1372 raise Exception('Invalid argument to scalingrange argument.')

1373

1374 if scalingmethod == 1:

1375 self.ydata = mavg_short/mavg_long * float(nshort)/float(nlong)

1376 elif scalingmethod in (2, 3):

1377 self.ydata = (mavg_short/mavg_long - 1.) \

1378 / ((float(nlong)/float(nshort)) - 1)

1379

1380 if scalingmethod == 3:

1381 self.ydata = num.maximum(self.ydata, 0.)

1382

1383 def sta_lta_right(self, tshort, tlong, quad=True, scalingmethod=1):

1384 '''

1385 Run special STA/LTA filter where the short time window is overlapping

1386 with the last part of the long time window.

1387

1388 :param tshort: length of short time window in [s]

1389 :param tlong: length of long time window in [s]

1390 :param quad: whether to square the data prior to applying the STA/LTA

1391 filter

1392 :param scalingmethod: integer key to select how output values are

1393 scaled / normalized (``1``, ``2``, or ``3``)

1394

1395 ============= ====================================== ===========

1396 Scalingmethod Implementation Range

1397 ============= ====================================== ===========

1398 ``1`` As/Al* Ts/Tl [0,1]

1399 ``2`` (As/Al - 1) / (Tl/Ts - 1) [-Ts/Tl,1]

1400 ``3`` Like ``2`` but clipping range at zero [0,1]

1401 ============= ====================================== ===========

1402

1403 With ``scalingmethod=1``, the values produced by this variant of the

1404 STA/LTA are equivalent to

1405

1406 .. math::

1407 s_i = \\frac{s}{l} \\frac{\\frac{1}{s}\\sum_{j=i}^{i+s-1} f_j}

1408 {\\frac{1}{l}\\sum_{j=i+s-l}^{i+s-1} f_j}

1409

1410 where :math:`f_j` are the input samples, :math:`s` are the number of

1411 samples in the short time window and :math:`l` are the number of

1412 samples in the long time window.

1413 '''

1414

1415 n = self.data_len()

1416 tmin = self.tmin

1417

1418 nshort = max(1, int(round(tshort/self.deltat)))

1419 nlong = max(1, int(round(tlong/self.deltat)))

1420

1421 assert nshort < nlong

1422

1423 if nlong > len(self.ydata):

1424 raise TraceTooShort(

1425 'Samples in trace: %s, samples needed: %s'

1426 % (len(self.ydata), nlong))

1427

1428 if quad:

1429 sqrdata = self.ydata**2

1430 else:

1431 sqrdata = self.ydata

1432

1433 nshift = int(math.floor(0.5 * (nlong - nshort)))

1434 if nlong % 2 != 0 and nshort % 2 == 0:

1435 nshift += 1

1436

1437 mavg_short = moving_avg(sqrdata, nshort)[nshift:]

1438 mavg_long = moving_avg(sqrdata, nlong)[:sqrdata.size-nshift]

1439

1440 self.drop_growbuffer()

1441

1442 if scalingmethod not in (1, 2, 3):

1443 raise Exception('Invalid argument to scalingrange argument.')

1444

1445 if scalingmethod == 1:

1446 ydata = mavg_short/mavg_long * nshort/nlong

1447 elif scalingmethod in (2, 3):

1448 ydata = (mavg_short/mavg_long - 1.) \

1449 / ((float(nlong)/float(nshort)) - 1)

1450

1451 if scalingmethod == 3:

1452 ydata = num.maximum(ydata, 0.)

1453

1454 self.set_ydata(ydata)

1455

1456 self.shift((math.ceil(0.5*nlong) - nshort + 1) * self.deltat)

1457

1458 self.chop(

1459 tmin + (nlong - nshort) * self.deltat,

1460 tmin + (n - nshort) * self.deltat)

1461

1462 def peaks(self, threshold, tsearch,

1463 deadtime=False,

1464 nblock_duration_detection=100):

1465

1466 '''

1467 Detect peaks above a given threshold (method 1).

1468

1469 From every instant, where the signal rises above ``threshold``, a time

1470 length of ``tsearch`` seconds is searched for a maximum. A list with

1471 tuples (time, value) for each detected peak is returned. The

1472 ``deadtime`` argument turns on a special deadtime duration detection

1473 algorithm useful in combination with recursive STA/LTA filters.

1474 '''

1475

1476 y = self.ydata

1477 above = num.where(y > threshold, 1, 0)

1478 deriv = num.zeros(y.size, dtype=num.int8)

1479 deriv[1:] = above[1:]-above[:-1]

1480 itrig_positions = num.nonzero(deriv > 0)[0]

1481 tpeaks = []

1482 apeaks = []

1483 tzeros = []

1484 tzero = self.tmin

1485

1486 for itrig_pos in itrig_positions:

1487 ibeg = itrig_pos

1488 iend = min(

1489 len(self.ydata),

1490 itrig_pos + int(math.ceil(tsearch/self.deltat)))

1491 ipeak = num.argmax(y[ibeg:iend])

1492 tpeak = self.tmin + (ipeak+ibeg)*self.deltat

1493 apeak = y[ibeg+ipeak]

1494

1495 if tpeak < tzero:

1496 continue

1497

1498 if deadtime:

1499 ibeg = itrig_pos

1500 iblock = 0

1501 nblock = nblock_duration_detection

1502 totalsum = 0.

1503 while True:

1504 if ibeg+iblock*nblock >= len(y):

1505 tzero = self.tmin + (len(y)-1) * self.deltat

1506 break

1507

1508 logy = num.log(

1509 y[ibeg+iblock*nblock:ibeg+(iblock+1)*nblock])

1510 logy[0] += totalsum

1511 ysum = num.cumsum(logy)

1512 totalsum = ysum[-1]

1513 below = num.where(ysum <= 0., 1, 0)

1514 deriv = num.zeros(ysum.size, dtype=num.int8)

1515 deriv[1:] = below[1:]-below[:-1]

1516 izero_positions = num.nonzero(deriv > 0)[0] + iblock*nblock

1517 if len(izero_positions) > 0:

1518 tzero = self.tmin + self.deltat * (

1519 ibeg + izero_positions[0])

1520 break

1521 iblock += 1

1522 else:

1523 tzero = ibeg*self.deltat + self.tmin + tsearch

1524

1525 tpeaks.append(tpeak)

1526 apeaks.append(apeak)

1527 tzeros.append(tzero)

1528

1529 if deadtime:

1530 return tpeaks, apeaks, tzeros

1531 else:

1532 return tpeaks, apeaks

1533

1534 def peaks2(self, threshold, tsearch):

1535

1536 '''

1537 Detect peaks above a given threshold (method 2).

1538

1539 This variant of peak detection is a bit more robust (and slower) than

1540 the one implemented in :py:meth:`Trace.peaks`. First all samples with

1541 ``a[i-1] < a[i] > a[i+1]`` are masked as potential peaks. From these,

1542 iteratively the one with the maximum amplitude ``a[j]`` and time

1543 ``t[j]`` is choosen and potential peaks within

1544 ``t[j] - tsearch, t[j] + tsearch``

1545 are discarded. The algorithm stops, when ``a[j] < threshold`` or when

1546 no more potential peaks are left.

1547 '''

1548

1549 a = num.copy(self.ydata)

1550

1551 amin = num.min(a)

1552

1553 a[0] = amin

1554 a[1: -1][num.diff(a, 2) <= 0.] = amin

1555 a[-1] = amin

1556

1557 data = []

1558 while True:

1559 imax = num.argmax(a)

1560 amax = a[imax]

1561

1562 if amax < threshold or amax == amin:

1563 break

1564

1565 data.append((self.tmin + imax * self.deltat, amax))

1566

1567 ntsearch = int(round(tsearch / self.deltat))

1568 a[max(imax-ntsearch//2, 0):min(imax+ntsearch//2, a.size)] = amin

1569

1570 if data:

1571 data.sort()

1572 tpeaks, apeaks = list(zip(*data))

1573 else:

1574 tpeaks, apeaks = [], []

1575

1576 return tpeaks, apeaks

1577

1578 def extend(self, tmin=None, tmax=None, fillmethod='zeros'):

1579 '''

1580 Extend trace to given span.

1581

1582 :param tmin: begin time of new span

1583 :param tmax: end time of new span

1584 :param fillmethod: ``'zeros'``, ``'repeat'``, ``'mean'``, or

1585 ``'median'``

1586 '''

1587

1588 nold = self.ydata.size

1589

1590 if tmin is not None:

1591 nl = min(0, int(round((tmin-self.tmin)/self.deltat)))

1592 else:

1593 nl = 0

1594

1595 if tmax is not None:

1596 nh = max(nold - 1, int(round((tmax-self.tmin)/self.deltat)))

1597 else:

1598 nh = nold - 1

1599

1600 n = nh - nl + 1

1601 data = num.zeros(n, dtype=self.ydata.dtype)

1602 data[-nl:-nl + nold] = self.ydata

1603 if self.ydata.size >= 1:

1604 if fillmethod == 'repeat':

1605 data[:-nl] = self.ydata[0]

1606 data[-nl + nold:] = self.ydata[-1]

1607 elif fillmethod == 'median':

1608 v = num.median(self.ydata)

1609 data[:-nl] = v

1610 data[-nl + nold:] = v

1611 elif fillmethod == 'mean':

1612 v = num.mean(self.ydata)

1613 data[:-nl] = v

1614 data[-nl + nold:] = v

1615

1616 self.drop_growbuffer()

1617 self.ydata = data

1618

1619 self.tmin += nl * self.deltat

1620 self.tmax = self.tmin + (self.ydata.size - 1) * self.deltat

1621

1622 self._update_ids()

1623

1624 def transfer(self,

1625 tfade=0.,

1626 freqlimits=None,

1627 transfer_function=None,

1628 cut_off_fading=True,

1629 demean=True,

1630 invert=False):

1631

1632 '''

1633 Return new trace with transfer function applied (convolution).

1634

1635 :param tfade: rise/fall time in seconds of taper applied in timedomain

1636 at both ends of trace.

1637 :param freqlimits: 4-tuple with corner frequencies in Hz.

1638 :param transfer_function: FrequencyResponse object; must provide a

1639 method 'evaluate(freqs)', which returns the transfer function

1640 coefficients at the frequencies 'freqs'.

1641 :param cut_off_fading: whether to cut off rise/fall interval in output

1642 trace.

1643 :param demean: remove mean before applying transfer function

1644 :param invert: set to True to do a deconvolution

1645 '''

1646

1647 if transfer_function is None:

1648 transfer_function = g_one_response

1649

1650 if self.tmax - self.tmin <= tfade*2.:

1651 raise TraceTooShort(

1652 'Trace %s.%s.%s.%s too short for fading length setting. '

1653 'trace length = %g, fading length = %g'

1654 % (self.nslc_id + (self.tmax-self.tmin, tfade)))

1655

1656 if freqlimits is None and (

1657 transfer_function is None or transfer_function.is_scalar()):

1658

1659 # special case for flat responses

1660

1661 output = self.copy()

1662 data = self.ydata

1663 ndata = data.size

1664

1665 if transfer_function is not None:

1666 c = num.abs(transfer_function.evaluate(num.ones(1))[0])

1667

1668 if invert:

1669 c = 1.0/c

1670

1671 data *= c

1672

1673 if tfade != 0.0:

1674 data *= costaper(

1675 0., tfade, self.deltat*(ndata-1)-tfade, self.deltat*ndata,

1676 ndata, self.deltat)

1677

1678 output.ydata = data

1679

1680 else:

1681 ndata = self.ydata.size

1682 ntrans = nextpow2(ndata*1.2)

1683 coeffs = self._get_tapered_coeffs(

1684 ntrans, freqlimits, transfer_function, invert=invert,

1685 demean=demean)

1686

1687 data = self.ydata

1688

1689 data_pad = num.zeros(ntrans, dtype=float)

1690 data_pad[:ndata] = data

1691 if demean:

1692 data_pad[:ndata] -= data.mean()

1693

1694 if tfade != 0.0:

1695 data_pad[:ndata] *= costaper(

1696 0., tfade, self.deltat*(ndata-1)-tfade, self.deltat*ndata,

1697 ndata, self.deltat)

1698

1699 fdata = num.fft.rfft(data_pad)

1700 fdata *= coeffs

1701 ddata = num.fft.irfft(fdata)

1702 output = self.copy()

1703 output.ydata = ddata[:ndata]

1704

1705 if cut_off_fading and tfade != 0.0:

1706 try:

1707 output.chop(output.tmin+tfade, output.tmax-tfade, inplace=True)

1708 except NoData:

1709 raise TraceTooShort(

1710 'Trace %s.%s.%s.%s too short for fading length setting. '

1711 'trace length = %g, fading length = %g'

1712 % (self.nslc_id + (self.tmax-self.tmin, tfade)))

1713 else:

1714 output.ydata = output.ydata.copy()

1715

1716 return output

1717

1718 def differentiate(self, n=1, order=4, inplace=True):

1719 '''

1720 Approximate first or second derivative of the trace.

1721

1722 :param n: 1 for first derivative, 2 for second

1723 :param order: order of the approximation 2 and 4 are supported

1724 :param inplace: if ``True`` the trace is differentiated in place,

1725 otherwise a new trace object with the derivative is returned.

1726

1727 Raises :py:exc:`ValueError` for unsupported `n` or `order`.

1728

1729 See :py:func:`~pyrocko.util.diff_fd` for implementation details.

1730 '''

1731

1732 ddata = util.diff_fd(n, order, self.deltat, self.ydata)

1733

1734 if inplace:

1735 self.ydata = ddata

1736 else:

1737 output = self.copy(data=False)

1738 output.set_ydata(ddata)

1739 return output

1740

1741 def drop_chain_cache(self):

1742 if self._pchain:

1743 self._pchain.clear()

1744

1745 def init_chain(self):

1746 self._pchain = pchain.Chain(

1747 do_downsample,

1748 do_extend,

1749 do_pre_taper,

1750 do_fft,

1751 do_filter,

1752 do_ifft)

1753

1754 def run_chain(self, tmin, tmax, deltat, setup, nocache):

1755 if setup.domain == 'frequency_domain':

1756 _, _, data = self._pchain(

1757 (self, deltat),

1758 (tmin, tmax),

1759 (setup.taper,),

1760 (setup.filter,),

1761 (setup.filter,),

1762 nocache=nocache)

1763

1764 return num.abs(data), num.abs(data)

1765

1766 else:

1767 processed = self._pchain(

1768 (self, deltat),

1769 (tmin, tmax),

1770 (setup.taper,),

1771 (setup.filter,),

1772 (setup.filter,),

1773 (),

1774 nocache=nocache)

1775

1776 if setup.domain == 'time_domain':

1777 data = processed.get_ydata()

1778

1779 elif setup.domain == 'envelope':

1780 processed = processed.envelope(inplace=False)

1781

1782 elif setup.domain == 'absolute':

1783 processed.set_ydata(num.abs(processed.get_ydata()))

1784

1785 return processed.get_ydata(), processed

1786

1787 def misfit(self, candidate, setup, nocache=False, debug=False):

1788 '''

1789 Calculate misfit and normalization factor against candidate trace.

1790

1791 :param candidate: :py:class:`Trace` object

1792 :param setup: :py:class:`MisfitSetup` object

1793 :returns: tuple ``(m, n)``, where m is the misfit value and n is the

1794 normalization divisor

1795

1796 If the sampling rates of ``self`` and ``candidate`` differ, the trace

1797 with the higher sampling rate will be downsampled.

1798 '''

1799

1800 a = self

1801 b = candidate

1802

1803 for tr in (a, b):

1804 if not tr._pchain:

1805 tr.init_chain()

1806

1807 deltat = max(a.deltat, b.deltat)

1808 tmin = min(a.tmin, b.tmin) - deltat

1809 tmax = max(a.tmax, b.tmax) + deltat

1810

1811 adata, aproc = a.run_chain(tmin, tmax, deltat, setup, nocache)

1812 bdata, bproc = b.run_chain(tmin, tmax, deltat, setup, nocache)

1813

1814 if setup.domain != 'cc_max_norm':

1815 m, n = Lx_norm(bdata, adata, norm=setup.norm)

1816 else:

1817 ctr = correlate(aproc, bproc, mode='full', normalization='normal')

1818 ccmax = ctr.max()[1]

1819 m = 0.5 - 0.5 * ccmax

1820 n = 0.5

1821

1822 if debug:

1823 return m, n, aproc, bproc

1824 else:

1825 return m, n

1826

1827 def spectrum(self, pad_to_pow2=False, tfade=None, ntrans_min=None):

1828 '''

1829 Get FFT spectrum of trace.

1830

1831 :param pad_to_pow2: whether to zero-pad the data to next larger

1832 power-of-two length

1833 :param tfade: ``None`` or a time length in seconds, to apply cosine

1834 shaped tapers to both

1835

1836 :returns: a tuple with (frequencies, values)

1837 '''

1838

1839 if ntrans_min is None:

1840 ndata = self.ydata.size

1841 else:

1842 ndata = ntrans_min

1843

1844 if pad_to_pow2:

1845 ntrans = nextpow2(ndata)

1846 else:

1847 ntrans = ndata

1848

1849 if tfade is None:

1850 ydata = self.ydata

1851 else:

1852 ydata = self.ydata * costaper(

1853 0., tfade, self.deltat*(ndata-1)-tfade, self.deltat*ndata,

1854 ndata, self.deltat)

1855

1856 fydata = num.fft.rfft(ydata, ntrans)

1857 df = 1./(ntrans*self.deltat)

1858 fxdata = num.arange(len(fydata))*df

1859 return fxdata, fydata

1860

1861 def multi_filter(self, filter_freqs, bandwidth):

1862

1863 class Gauss(FrequencyResponse):

1864 f0 = Float.T()

1865 a = Float.T(default=1.0)

1866

1867 def __init__(self, f0, a=1.0, **kwargs):

1868 FrequencyResponse.__init__(self, f0=f0, a=a, **kwargs)

1869

1870 def evaluate(self, freqs):

1871 omega0 = 2.*math.pi*self.f0

1872 omega = 2.*math.pi*freqs

1873 return num.exp(-((omega-omega0)

1874 / (self.a*omega0))**2)

1875

1876 freqs, coeffs = self.spectrum()

1877 n = self.data_len()

1878 nfilt = len(filter_freqs)

1879 signal_tf = num.zeros((nfilt, n))

1880 centroid_freqs = num.zeros(nfilt)

1881 for ifilt, f0 in enumerate(filter_freqs):

1882 taper = Gauss(f0, a=bandwidth)

1883 weights = taper.evaluate(freqs)

1884 nhalf = freqs.size

1885 analytic_spec = num.zeros(n, dtype=complex)

1886 analytic_spec[:nhalf] = coeffs*weights

1887

1888 enorm = num.abs(analytic_spec[:nhalf])**2

1889 enorm /= num.sum(enorm)

1890

1891 if n % 2 == 0:

1892 analytic_spec[1:nhalf-1] *= 2.

1893 else:

1894 analytic_spec[1:nhalf] *= 2.

1895

1896 analytic = num.fft.ifft(analytic_spec)

1897 signal_tf[ifilt, :] = num.abs(analytic)

1898

1899 enorm = num.abs(analytic_spec[:nhalf])**2

1900 enorm /= num.sum(enorm)

1901 centroid_freqs[ifilt] = num.sum(freqs*enorm)

1902

1903 return centroid_freqs, signal_tf

1904

1905 def _get_tapered_coeffs(

1906 self, ntrans, freqlimits, transfer_function, invert=False,

1907 demean=True):

1908

1909 cache_key = (

1910 ntrans, self.deltat, freqlimits, transfer_function.uuid, invert,

1911 demean)

1912

1913 if cache_key in g_tapered_coeffs_cache:

1914 return g_tapered_coeffs_cache[cache_key]

1915

1916 deltaf = 1./(self.deltat*ntrans)

1917 nfreqs = ntrans//2 + 1

1918 transfer = num.ones(nfreqs, dtype=complex)

1919 hi = snapper(nfreqs, deltaf)

1920 if freqlimits is not None:

1921 kmin, kmax = hi(freqlimits[0]), hi(freqlimits[3])

1922 freqs = num.arange(kmin, kmax)*deltaf

1923 coeffs = transfer_function.evaluate(freqs)

1924 if invert:

1925 if num.any(coeffs == 0.0):

1926 raise InfiniteResponse('%s.%s.%s.%s' % self.nslc_id)

1927

1928 transfer[kmin:kmax] = 1.0 / coeffs

1929 else:

1930 transfer[kmin:kmax] = coeffs

1931

1932 tapered_transfer = costaper(*freqlimits, nfreqs, deltaf) * transfer

1933 else:

1934 if invert:

1935 raise Exception(

1936 'transfer: `freqlimits` must be given when `invert` is '

1937 'set to `True`')

1938

1939 freqs = num.arange(nfreqs) * deltaf

1940 tapered_transfer = transfer_function.evaluate(freqs)

1941

1942 g_tapered_coeffs_cache[cache_key] = tapered_transfer

1943

1944 if demean:

1945 tapered_transfer[0] = 0.0 # don't introduce static offsets

1946

1947 return tapered_transfer

1948

1949 def fill_template(self, template, **additional):

1950 '''

1951 Fill string template with trace metadata.

1952

1953 Uses normal python '%(placeholder)s' string templates. The following

1954 placeholders are considered: ``network``, ``station``, ``location``,

1955 ``channel``, ``tmin`` (time of first sample), ``tmax`` (time of last

1956 sample), ``tmin_ms``, ``tmax_ms``, ``tmin_us``, ``tmax_us``,

1957 ``tmin_year``, ``tmax_year``, ``tmin_month``, ``tmax_month``,

1958 ``tmin_day``, ``tmax_day``, ``julianday``. The variants ending with

1959 ``'_ms'`` include milliseconds, those with ``'_us'`` include

1960 microseconds, those with ``'_year'`` contain only the year.

1961 '''

1962

1963 template = template.replace('%n', '%(network)s')\

1964 .replace('%s', '%(station)s')\

1965 .replace('%l', '%(location)s')\

1966 .replace('%c', '%(channel)s')\

1967 .replace('%b', '%(tmin)s')\

1968 .replace('%e', '%(tmax)s')\

1969 .replace('%j', '%(julianday)s')

1970

1971 return template % TraceStringFiller(self, additional=additional)

1972

1973 def plot(self):

1974 '''

1975 Show trace with matplotlib.

1976

1977 See also: :py:meth:`Trace.snuffle`.

1978 '''

1979

1980 import pylab

1981 pylab.plot(self.get_xdata(), self.get_ydata())

1982 name = '%s %s %s - %s' % (

1983 self.channel,

1984 self.station,

1985 time.strftime('%d-%m-%y %H:%M:%S', time.gmtime(self.tmin)),

1986 time.strftime('%d-%m-%y %H:%M:%S', time.gmtime(self.tmax)))

1987

1988 pylab.title(name)

1989 pylab.show()

1990

1991 def snuffle(self, **kwargs):

1992 '''

1993 Show trace in a snuffler window.

1994

1995 :param stations: list of :py:class:`pyrocko.model.station.Station`

1996 objects or ``None``

1997 :param events: list of :py:class:`pyrocko.model.event.Event` objects or

1998 ``None``

1999 :param markers: list of :py:class:`pyrocko.gui.snuffler.marker.Marker`

2000 objects or ``None``

2001 :param ntracks: float, number of tracks to be shown initially (default:

2002 12)

2003 :param follow: time interval (in seconds) for real time follow mode or

2004 ``None``

2005 :param controls: bool, whether to show the main controls (default:

2006 ``True``)

2007 :param opengl: bool, whether to use opengl (default: ``False``)

2008 '''

2009

2010 return snuffle([self], **kwargs)

2011

2012

2013def snuffle(traces, **kwargs):

2014 '''

2015 Show traces in a snuffler window.

2016

2017 :param stations: list of :py:class:`pyrocko.model.station.Station` objects

2018 or ``None``

2019 :param events: list of :py:class:`pyrocko.model.event.Event` objects or

2020 ``None``

2021 :param markers: list of :py:class:`pyrocko.gui.snuffler.marker.Marker`

2022 objects or ``None``

2023 :param ntracks: int, number of tracks to be shown initially (default: 12)

2024 :param follow: time interval (in seconds) for real time follow mode or

2025 ``None``

2026 :param controls: bool, whether to show the main controls (default:

2027 ``True``)

2028 :param opengl: bool, whether to use opengl (default: ``False``)

2029 '''

2030

2031 from pyrocko import pile

2032 from pyrocko.gui.snuffler import snuffler

2033 p = pile.Pile()

2034 if traces:

2035 trf = pile.MemTracesFile(None, traces)

2036 p.add_file(trf)

2037 return snuffler.snuffle(p, **kwargs)

2038

2039

2040def downsample_tpad(

2041 deltat_in, deltat_out, allow_upsample_max=1, ftype='fir-remez'):

2042 '''

2043 Get approximate amount of cutoff which will be produced by downsampling.

2044

2045 The :py:meth:`Trace.downsample_to` method removes some samples at the

2046 beginning and end of the trace which is downsampled. This function

2047 estimates the approximate length [s] which will be cut off for a given set

2048 of parameters supplied to :py:meth:`Trace.downsample_to`.

2049

2050 :param deltat_in:

2051 Input sampling interval [s].

2052 :type deltat_in:

2053 float

2054

2055 :param deltat_out:

2056 Output samling interval [s].

2057 :type deltat_out:

2058 float

2059

2060 :returns:

2061 Approximate length [s] which will be cut off.

2062

2063 See :py:meth:`Trace.downsample_to` for details.

2064 '''

2065

2066 upsratio, deci_seq = _configure_downsampling(

2067 deltat_in, deltat_out, allow_upsample_max)

2068

2069 tpad = 0.0

2070 deltat = deltat_in / upsratio

2071 for deci in deci_seq:

2072 b, a, n = util.decimate_coeffs(deci, None, ftype)

2073 # n//2 for the antialiasing

2074 # +deci for possible snap to multiples

2075 # +1 for rounding errors

2076 tpad += (n//2 + deci + 1) * deltat

2077 deltat = deltat * deci

2078

2079 return tpad

2080

2081

2082def _configure_downsampling(deltat_in, deltat_out, allow_upsample_max):

2083 for upsratio in range(1, allow_upsample_max+1):

2084 dratio = (upsratio/deltat_in) / (1./deltat_out)

2085 deci_seq = util.decitab(int(round(dratio)))

2086 if abs(dratio - round(dratio)) / dratio < 0.0001 and deci_seq:

2087 return upsratio, [deci for deci in deci_seq if deci != 1]

2088

2089 raise util.UnavailableDecimation('ratio = %g' % (deltat_out / deltat_in))

2090

2091

2092def _all_same(xs):

2093 return all(x == xs[0] for x in xs)

2094

2095

2096def _incompatibilities(traces):

2097 if not traces:

2098 return None

2099

2100 params = [

2101 (tr.ydata.size, tr.ydata.dtype, tr.deltat, tr.tmin)

2102 for tr in traces]

2103

2104 if not _all_same(params):

2105 return params

2106 else:

2107 return None

2108

2109

2110def _raise_incompatible_traces(params):

2111 raise IncompatibleTraces(

2112 'Given traces are incompatible. Sampling rate, start time, '

2113 'number of samples and data type must match.\n%s\n%s' % (

2114 ' %10s %-10s %12s %22s' % (

2115 'nsamples', 'dtype', 'deltat', 'tmin'),

2116 '\n'.join(

2117 ' %10i %-10s %12.5e %22s' % (

2118 nsamples, dtype, deltat, util.time_to_str(tmin))

2119 for (nsamples, dtype, deltat, tmin) in params)))

2120

2121

2122def _ensure_compatible(traces):

2123 params = _incompatibilities(traces)

2124 if params:

2125 _raise_incompatible_traces(params)

2126

2127

2128def _almost_equal(a, b, atol):

2129 return abs(a-b) < atol

2130

2131

2132def get_traces_data_as_array(traces):

2133 '''

2134 Merge data samples from multiple traces into a 2D array.

2135

2136 :param traces:

2137 Input waveforms.

2138 :type traces:

2139 list of :py:class:`pyrocko.Trace <pyrocko.trace.Trace>` objects

2140

2141 :raises:

2142 :py:class:`IncompatibleTraces` if traces have different time

2143 span, sample rate or data type, or if traces is an empty list.

2144

2145 :returns:

2146 2D array as ``data[itrace, isample]``.

2147 :rtype:

2148 :py:class:`numpy.ndarray`

2149 '''

2150

2151 if not traces:

2152 raise IncompatibleTraces('Need at least one trace.')

2153

2154 _ensure_compatible(traces)

2155

2156 return num.vstack([tr.ydata for tr in traces])

2157

2158

2159def make_traces_compatible(

2160 traces,

2161 dtype=None,

2162 deltat=None,

2163 enforce_global_snap=True,

2164 warn_snap=False):

2165

2166 '''

2167 Homogenize sampling rate, time span, sampling instants, and data type.

2168

2169 This function takes a group of traces and tries to make them compatible in

2170 terms of data type and sampling rate, time span, and sampling instants of

2171 time.

2172

2173 If necessary, traces are (in order):

2174

2175 - casted to the same data type.

2176 - downsampled to a common sampling rate, using decimation cascades.

2177 - resampled to common sampling instants in time, using Sinc interpolation.

2178 - cut to the same time span. The longest time span covered by all traces is

2179 used.

2180

2181 :param traces:

2182 Input waveforms.

2183 :type traces:

2184 :py:class:`list` of :py:class:`Trace`

2185

2186 :param dtype:

2187 Force traces to be casted to the given data type. If not specified, the

2188 traces are cast to :py:class:`float`.

2189 :type dtype:

2190 :py:class:`numpy.dtype`

2191

2192 :param deltat:

2193 Sampling interval [s]. If not specified, the longest sampling interval

2194 among the input traces is chosen.

2195 :type deltat:

2196 float

2197

2198 :param enforce_global_snap:

2199 If ``True``, choose sampling instants to be even multiples of the

2200 sampling rate in system time. When set to ``False`` traces are still

2201 resampled to common time instants (if necessary), but they may be

2202 offset to the system time sampling rate multiples.

2203 :type enforce_global_snap:

2204 bool

2205

2206 :param warn_snap:

2207 If set to ``True`` warn, when resampling has to be performed.

2208 :type warn_snap:

2209 bool

2210 '''

2211

2212 eps_snap = 1e-3

2213

2214 if not traces:

2215 return []

2216

2217 traces = list(traces)

2218

2219 dtypes = [tr.ydata.dtype for tr in traces]

2220 if not _all_same(dtypes) or dtype is not None:

2221

2222 if dtype is None:

2223 dtype = float

2224 logger.warning(

2225 'make_traces_compatible: Inconsistent data types - converting '

2226 'sample datatype to %s.' % str(dtype))

2227

2228 for itr, tr in enumerate(traces):

2229 tr_copy = tr.copy(data=False)

2230 tr_copy.set_ydata(tr.ydata.astype(dtype))

2231 traces[itr] = tr_copy

2232

2233 deltats = [tr.deltat for tr in traces]

2234 if not _all_same(deltats) or deltat is not None:

2235 if deltat is None:

2236 deltat = max(deltats)

2237 logger.warning(

2238 'make_traces_compatible: Inconsistent sampling rates - '

2239 'downsampling to lowest rate among input traces: %g Hz.'

2240 % (1.0 / deltat))

2241

2242 for itr, tr in enumerate(traces):

2243 if tr.deltat != deltat:

2244 tr_copy = tr.copy()

2245 tr_copy.downsample_to(deltat, snap=True, cut=True)

2246 traces[itr] = tr_copy

2247

2248 tmins = num.array([tr.tmin for tr in traces])

2249 is_aligned = num.abs(num.round(tmins / deltat) * deltat - tmins) \

2250 > deltat * eps_snap

2251

2252 if enforce_global_snap or any(is_aligned):

2253 tref = util.to_time_float(0.0)

2254 else:

2255 # to keep a common subsample shift

2256 tref = num.max(tmins)

2257

2258 tmins_snap = num.round((tmins - tref) / deltat) * deltat + tref

2259 need_snap = num.abs(tmins_snap - tmins) > deltat * eps_snap

2260 if num.any(need_snap):

2261 if warn_snap:

2262 logger.warning(

2263 'make_traces_compatible: Misaligned sampling - introducing '

2264 'subsample shifts for proper alignment.')

2265

2266 for itr, tr in enumerate(traces):

2267 if need_snap[itr]:

2268 tr_copy = tr.copy()

2269 if tref != 0.0:

2270 tr_copy.shift(-tref)

2271

2272 tr_copy.snap(interpolate=True)

2273 if tref != 0.0:

2274 tr_copy.shift(tref)

2275

2276 traces[itr] = tr_copy

2277

2278 tmins = num.array([tr.tmin for tr in traces])

2279 nsamples = num.array([tr.ydata.size for tr in traces])

2280 tmaxs = tmins + (nsamples - 1) * deltat

2281

2282 tmin = num.max(tmins)

2283 tmax = num.min(tmaxs)

2284

2285 if tmin > tmax:

2286 raise IncompatibleTraces('Traces do not overlap.')

2287

2288 nsamples_must = int(round((tmax - tmin) / deltat)) + 1

2289 for itr, tr in enumerate(traces):

2290 if not (_almost_equal(tr.tmin, tmin, deltat*eps_snap)

2291 and _almost_equal(tr.tmax, tmax, deltat*eps_snap)):

2292

2293 traces[itr] = tr.chop(

2294 tmin, tmax,

2295 inplace=False,

2296 want_incomplete=False,

2297 include_last=True)

2298

2299 xtr = traces[itr]

2300 assert _almost_equal(xtr.tmin, tmin, deltat*eps_snap)

2301 assert int(round((xtr.tmax - xtr.tmin) / deltat)) + 1 == nsamples_must

2302 xtr.tmin = tmin

2303 xtr.tmax = tmax

2304 xtr.deltat = deltat

2305 xtr._update_ids()

2306

2307 return traces

2308

2309

2310class IncompatibleTraces(Exception):

2311 '''

2312 Raised when traces have incompatible sampling rate, time span or data type.

2313 '''

2314

2315

2316class InfiniteResponse(Exception):

2317 '''

2318 This exception is raised by :py:class:`Trace` operations when deconvolution

2319 of a frequency response (instrument response transfer function) would

2320 result in a division by zero.

2321 '''

2322

2323

2324class MisalignedTraces(Exception):

2325 '''

2326 This exception is raised by some :py:class:`Trace` operations when tmin,

2327 tmax or number of samples do not match.

2328 '''

2329

2330 pass

2331

2332

2333class OverlappingTraces(Exception):

2334 '''

2335 This exception is raised by some :py:class:`Trace` operations when traces

2336 overlap but should not.

2337 '''

2338

2339 pass

2340

2341

2342class NoData(Exception):

2343 '''

2344 This exception is raised by some :py:class:`Trace` operations when no or

2345 not enough data is available.

2346 '''

2347

2348 pass

2349

2350

2351class AboveNyquist(Exception):

2352 '''

2353 This exception is raised by some :py:class:`Trace` operations when given

2354 frequencies are above the Nyquist frequency.

2355 '''

2356

2357 pass

2358

2359

2360class TraceTooShort(Exception):

2361 '''

2362 This exception is raised by some :py:class:`Trace` operations when the

2363 trace is too short.

2364 '''

2365

2366 pass

2367

2368

2369class ResamplingFailed(Exception):

2370 pass

2371

2372

2373class TraceStringFiller:

2374

2375 def __init__(self, tr, additional={}):

2376 self.tr = tr

2377 self.additional = additional

2378

2379 def __getitem__(self, k):

2380 if k in ('network', 'station', 'location', 'channel', 'extra'):

2381 return getattr(self.tr, k)

2382

2383 method = getattr(self, 'get_' + k, None)

2384 if method:

2385 return method()

2386

2387 return self.additional[k]

2388

2389 def _filename_safe(self, s):

2390 return re.sub(r'[^0-9A-Za-z_-]', '_', s)

2391

2392 def get_network_safe(self):

2393 return self._filename_safe(self.tr.network)

2394

2395 def get_station_safe(self):

2396 return self._filename_safe(self.tr.station)

2397

2398 def get_location_safe(self):

2399 return self._filename_safe(self.tr.location)

2400

2401 def get_channel_safe(self):

2402 return self._filename_safe(self.tr.channel)

2403

2404 def get_extra_safe(self):

2405 return self._filename_safe(self.tr.extra)

2406

2407 def get_network_dsafe(self):

2408 return self._filename_safe(self.tr.network) or '_'

2409

2410 def get_station_dsafe(self):

2411 return self._filename_safe(self.tr.station) or '_'

2412

2413 def get_location_dsafe(self):

2414 return self._filename_safe(self.tr.location) or '_'

2415

2416 def get_channel_dsafe(self):

2417 return self._filename_safe(self.tr.channel) or '_'

2418

2419 def get_extra_dsafe(self):

2420 return self._filename_safe(self.tr.extra) or '_'

2421

2422 def get_tmin(self):

2423 return util.time_to_str(self.tr.tmin, format='%Y-%m-%d_%H-%M-%S')

2424

2425 def get_tmax(self):

2426 return util.time_to_str(self.tr.tmax, format='%Y-%m-%d_%H-%M-%S')

2427

2428 def get_tmin_ms(self):

2429 return util.time_to_str(self.tr.tmin, format='%Y-%m-%d_%H-%M-%S.3FRAC')

2430

2431 def get_tmax_ms(self):

2432 return util.time_to_str(self.tr.tmax, format='%Y-%m-%d_%H-%M-%S.3FRAC')

2433

2434 def get_tmin_us(self):

2435 return util.time_to_str(self.tr.tmin, format='%Y-%m-%d_%H-%M-%S.6FRAC')

2436

2437 def get_tmax_us(self):

2438 return util.time_to_str(self.tr.tmax, format='%Y-%m-%d_%H-%M-%S.6FRAC')

2439

2440 def get_tmin_year(self):

2441 return util.time_to_str(self.tr.tmin, format='%Y')

2442

2443 def get_tmin_month(self):

2444 return util.time_to_str(self.tr.tmin, format='%m')

2445

2446 def get_tmin_day(self):

2447 return util.time_to_str(self.tr.tmin, format='%d')

2448

2449 def get_tmax_year(self):

2450 return util.time_to_str(self.tr.tmax, format='%Y')

2451

2452 def get_tmax_month(self):

2453 return util.time_to_str(self.tr.tmax, format='%m')

2454

2455 def get_tmax_day(self):

2456 return util.time_to_str(self.tr.tmax, format='%d')

2457

2458 def get_julianday(self):

2459 return str(util.julian_day_of_year(self.tr.tmin))

2460

2461 def get_tmin_jday(self):

2462 return util.time_to_str(self.tr.tmin, format='%j')

2463

2464 def get_tmax_jday(self):

2465 return util.time_to_str(self.tr.tmax, format='%j')

2466

2467

2468def minmax(traces, key=None, mode='minmax', outer_mode='minmax'):

2469

2470 '''

2471 Get data range given traces grouped by selected pattern.

2472

2473 :param key: a callable which takes as single argument a trace and returns a

2474 key for the grouping of the results. If this is ``None``, the default,

2475 ``lambda tr: (tr.network, tr.station, tr.location, tr.channel)`` is

2476 used.

2477 :param mode: ``'minmax'`` or floating point number. If this is

2478 ``'minmax'``, minimum and maximum of the traces are used, if it is a

2479 number, mean +- standard deviation times ``mode`` is used.

2480

2481 param outer_mode: ``'minmax'`` to use mininum and maximum of the

2482 single-trace ranges, or ``'robust'`` to use the interval to discard 10%

2483 extreme values on either end.

2484

2485 :returns: a dict with the combined data ranges.

2486

2487 Examples::

2488

2489 ranges = minmax(traces, lambda tr: tr.channel)

2490 print ranges['N'] # print min & max of all traces with channel == 'N'

2491 print ranges['E'] # print min & max of all traces with channel == 'E'

2492

2493 ranges = minmax(traces, lambda tr: (tr.network, tr.station))

2494 print ranges['GR', 'HAM3'] # print min & max of all traces with

2495 # network == 'GR' and station == 'HAM3'

2496

2497 ranges = minmax(traces, lambda tr: None)

2498 print ranges[None] # prints min & max of all traces

2499 '''

2500

2501 if key is None:

2502 key = _default_key

2503

2504 ranges = defaultdict(list)

2505 for trace in traces:

2506 if isinstance(mode, str) and mode == 'minmax':

2507 mi, ma = num.nanmin(trace.ydata), num.nanmax(trace.ydata)

2508 else:

2509 mean = trace.ydata.mean()

2510 std = trace.ydata.std()

2511 mi, ma = mean-std*mode, mean+std*mode

2512

2513 k = key(trace)

2514 ranges[k].append((mi, ma))

2515

2516 for k in ranges:

2517 mins, maxs = num.array(ranges[k]).T

2518 if outer_mode == 'minmax':

2519 ranges[k] = num.nanmin(mins), num.nanmax(maxs)

2520 elif outer_mode == 'robust':

2521 ranges[k] = num.percentile(mins, 10.), num.percentile(maxs, 90.)

2522

2523 return ranges

2524

2525

2526def minmaxtime(traces, key=None):

2527

2528 '''

2529 Get time range given traces grouped by selected pattern.

2530

2531 :param key: a callable which takes as single argument a trace and returns a

2532 key for the grouping of the results. If this is ``None``, the default,

2533 ``lambda tr: (tr.network, tr.station, tr.location, tr.channel)`` is

2534 used.

2535

2536 :returns: a dict with the combined data ranges.

2537 '''

2538

2539 if key is None:

2540 key = _default_key

2541

2542 ranges = {}

2543 for trace in traces:

2544 mi, ma = trace.tmin, trace.tmax

2545 k = key(trace)

2546 if k not in ranges:

2547 ranges[k] = mi, ma

2548 else:

2549 tmi, tma = ranges[k]

2550 ranges[k] = min(tmi, mi), max(tma, ma)

2551

2552 return ranges

2553

2554

2555def degapper(

2556 traces,

2557 maxgap=5,

2558 fillmethod='interpolate',

2559 deoverlap='use_second',

2560 maxlap=None):

2561

2562 '''

2563 Try to connect traces and remove gaps.

2564

2565 This method will combine adjacent traces, which match in their network,

2566 station, location and channel attributes. Overlapping parts are handled

2567 according to the ``deoverlap`` argument.

2568

2569 :param traces: input traces, must be sorted by their full_id attribute.

2570 :param maxgap: maximum number of samples to interpolate.

2571 :param fillmethod: what to put into the gaps: 'interpolate' or 'zeros'.

2572 :param deoverlap: how to handle overlaps: 'use_second' to use data from

2573 second trace (default), 'use_first' to use data from first trace,

2574 'crossfade_cos' to crossfade with cosine taper, 'add' to add amplitude

2575 values.

2576 :param maxlap: maximum number of samples of overlap which are removed

2577

2578 :returns: list of traces

2579 '''

2580

2581 in_traces = []

2582 out_traces = []

2583 for tr in traces:

2584 if tr.deltat == 0:

2585 out_traces.append(tr)

2586 else:

2587 in_traces.append(tr)

2588

2589 if not in_traces:

2590 return out_traces

2591

2592 out_traces.append(in_traces.pop(0))

2593

2594 while in_traces:

2595

2596 a = out_traces[-1]

2597 b = in_traces.pop(0)

2598

2599 avirt, bvirt = a.ydata is None, b.ydata is None

2600 assert avirt == bvirt, \

2601 'traces given to degapper() must either all have data or have ' \

2602 'no data.'

2603

2604 virtual = avirt and bvirt

2605

2606 if (a.nslc_id == b.nslc_id and a.deltat == b.deltat

2607 and a.data_len() >= 1 and b.data_len() >= 1

2608 and (virtual or a.ydata.dtype == b.ydata.dtype)):

2609

2610 dist = (b.tmin-(a.tmin+(a.data_len()-1)*a.deltat))/a.deltat

2611 idist = int(round(dist))

2612 if abs(dist - idist) > 0.05 and idist <= maxgap:

2613 # logger.warning('Cannot degap traces with displaced sampling '

2614 # '(%s, %s, %s, %s)' % a.nslc_id)

2615

2616 if idist <= 0 and (maxlap is None or -maxlap < idist):

2617 # still cut off overlap (cut off on first trace)

2618 try:

2619 a.chop(a.tmin, max(a.tmin, b.tmin-b.deltat))

2620 except NoData:

2621 pass

2622

2623 pass

2624

2625 else:

2626 if 1 < idist <= maxgap:

2627 if not virtual:

2628 if fillmethod == 'interpolate':

2629 filler = a.ydata[-1] + (

2630 ((1.0 + num.arange(idist-1, dtype=float))

2631 / idist) * (b.ydata[0]-a.ydata[-1])

2632 ).astype(a.ydata.dtype)

2633 elif fillmethod == 'zeros':

2634 filler = num.zeros(idist-1, dtype=a.ydata.dtype)

2635 a.ydata = num.concatenate((a.ydata, filler, b.ydata))

2636 a.tmax = b.tmax

2637 if a.mtime and b.mtime:

2638 a.mtime = max(a.mtime, b.mtime)

2639 continue

2640

2641 elif idist == 1:

2642 if not virtual:

2643 a.ydata = num.concatenate((a.ydata, b.ydata))

2644 a.tmax = b.tmax

2645 if a.mtime and b.mtime:

2646 a.mtime = max(a.mtime, b.mtime)

2647 continue

2648

2649 elif idist <= 0 and (maxlap is None or -maxlap < idist):

2650 if b.tmax > a.tmax:

2651 if not virtual:

2652 na = a.ydata.size

2653 n = -idist+1

2654 if deoverlap == 'use_second':

2655 a.ydata = num.concatenate(

2656 (a.ydata[:-n], b.ydata))

2657 elif deoverlap in ('use_first', 'crossfade_cos'):

2658 a.ydata = num.concatenate(

2659 (a.ydata, b.ydata[n:]))

2660 elif deoverlap == 'add':

2661 a.ydata[-n:] += b.ydata[:n]

2662 a.ydata = num.concatenate(

2663 (a.ydata, b.ydata[n:]))

2664 else:

2665 assert False, 'unknown deoverlap method'

2666

2667 if deoverlap == 'crossfade_cos':

2668 n = -idist+1

2669 taper = 0.5-0.5*num.cos(

2670 (1.+num.arange(n))/(1.+n)*num.pi)

2671 a.ydata[na-n:na] *= 1.-taper

2672 a.ydata[na-n:na] += b.ydata[:n] * taper

2673

2674 a.tmax = b.tmax

2675 if a.mtime and b.mtime:

2676 a.mtime = max(a.mtime, b.mtime)

2677 continue

2678 else:

2679 # make short second trace vanish

2680 continue

2681

2682 if b.data_len() >= 1:

2683 out_traces.append(b)

2684

2685 for tr in out_traces:

2686 tr._update_ids()

2687

2688 return out_traces

2689

2690

2691def deoverlap(traces):

2692 groups = util.group_by(lambda tr: tr.codes, traces)

2693 traces_out = []

2694 for codes, group in groups.items():

2695 group_out = []

2696 group.sort(key=lambda tr: -(tr.tmax - tr.tmin))

2697 for tr in group:

2698 keep = tr

2699 for have in group_out:

2700 if tr.tmax < have.tmin or tr.tmin > have.tmax:

2701 pass

2702 elif tr.tmin < have.tmin and tr.tmax >= have.tmin:

2703 if keep is tr:

2704 keep = keep.copy()

2705 try:

2706 keep.chop(tr.tmin, have.tmin)

2707 except NoData:

2708 keep = None

2709 break

2710 elif tr.tmin >= have.tmin and tr.tmax <= have.tmax:

2711 keep = None

2712 break

2713 elif tr.tmin <= have.tmax and tr.tmax > have.tmax:

2714 if keep is tr:

2715 keep = keep.copy()

2716 try:

2717 keep.chop(have.tmax+tr.deltat, tr.tmax+tr.deltat)

2718 except NoData:

2719 keep = None

2720 break

2721 else:

2722 print(tr.tmin, tr.tmax, have.tmin, have.tmax)

2723 assert False

2724

2725 if keep is not None:

2726 group_out.append(keep)

2727

2728 traces_out.extend(sorted(group_out, key=lambda tr: tr.tmin))

2729

2730 return traces_out

2731

2732

2733def rotate(traces, azimuth, in_channels, out_channels):

2734 '''

2735 2D rotation of traces.

2736

2737 :param traces: list of input traces

2738 :param azimuth: difference of the azimuths of the component directions

2739 (azimuth of out_channels[0]) - (azimuth of in_channels[0])

2740 :param in_channels: names of the input channels (e.g. 'N', 'E')

2741 :param out_channels: names of the output channels (e.g. 'R', 'T')

2742 :returns: list of rotated traces

2743 '''

2744

2745 phi = azimuth/180.*math.pi

2746 cphi = math.cos(phi)

2747 sphi = math.sin(phi)

2748 rotated = []

2749 in_channels = tuple(_channels_to_names(in_channels))

2750 out_channels = tuple(_channels_to_names(out_channels))

2751 for a in traces:

2752 for b in traces:

2753 if ((a.channel, b.channel) == in_channels and

2754 a.nslc_id[:3] == b.nslc_id[:3] and

2755 abs(a.deltat-b.deltat) < a.deltat*0.001):

2756 tmin = max(a.tmin, b.tmin)

2757 tmax = min(a.tmax, b.tmax)

2758

2759 if tmin < tmax:

2760 ac = a.chop(tmin, tmax, inplace=False, include_last=True)

2761 bc = b.chop(tmin, tmax, inplace=False, include_last=True)

2762 if abs(ac.tmin - bc.tmin) > ac.deltat*0.01:

2763 logger.warning(

2764 'Cannot rotate traces with displaced sampling '

2765 '(%s, %s, %s, %s)' % a.nslc_id)

2766 continue

2767

2768 acydata = ac.get_ydata()*cphi+bc.get_ydata()*sphi

2769 bcydata = -ac.get_ydata()*sphi+bc.get_ydata()*cphi

2770 ac.set_ydata(acydata)

2771 bc.set_ydata(bcydata)

2772

2773 ac.set_codes(channel=out_channels[0])

2774 bc.set_codes(channel=out_channels[1])

2775 rotated.append(ac)

2776 rotated.append(bc)

2777

2778 return rotated

2779

2780

2781def rotate_to_rt(n, e, source, receiver, out_channels=('R', 'T')):

2782 '''

2783 Rotate traces from NE to RT system.

2784

2785 :param n:

2786 North trace.

2787 :type n:

2788 :py:class:`~pyrocko.trace.Trace`

2789

2790 :param e:

2791 East trace.

2792 :type e:

2793 :py:class:`~pyrocko.trace.Trace`

2794

2795 :param source:

2796 Source of the recorded signal.

2797 :type source:

2798 :py:class:`pyrocko.gf.seismosizer.Source`

2799

2800 :param receiver:

2801 Receiver of the recorded signal.

2802 :type receiver:

2803 :py:class:`pyrocko.model.location.Location`

2804

2805 :param out_channels:

2806 Channel codes of the output channels (radial, transversal).

2807 Default is ('R', 'T').

2808

2809 :type out_channels

2810 optional, tuple[str, str]

2811

2812 :returns:

2813 Rotated traces (radial, transversal).

2814 :rtype:

2815 tuple[

2816 :py:class:`~pyrocko.trace.Trace`,

2817 :py:class:`~pyrocko.trace.Trace`]

2818 '''

2819 azimuth = orthodrome.azimuth(receiver, source) + 180.

2820 in_channels = n.channel, e.channel

2821 out = rotate(

2822 [n, e], azimuth,

2823 in_channels=in_channels,

2824 out_channels=out_channels)

2825

2826 assert len(out) == 2

2827 for tr in out:

2828 if tr.channel == out_channels[0]:

2829 r = tr

2830 elif tr.channel == out_channels[1]:

2831 t = tr

2832 else:

2833 assert False

2834

2835 return r, t

2836

2837

2838def rotate_to_lqt(traces, backazimuth, incidence, in_channels,

2839 out_channels=('L', 'Q', 'T')):

2840 '''

2841 Rotate traces from ZNE to LQT system.

2842

2843 :param traces: list of traces in arbitrary order

2844 :param backazimuth: backazimuth in degrees clockwise from north

2845 :param incidence: incidence angle in degrees from vertical

2846 :param in_channels: input channel names

2847 :param out_channels: output channel names (default: ('L', 'Q', 'T'))

2848 :returns: list of transformed traces

2849 '''

2850 i = incidence/180.*num.pi

2851 b = backazimuth/180.*num.pi

2852

2853 ci = num.cos(i)

2854 cb = num.cos(b)

2855 si = num.sin(i)

2856 sb = num.sin(b)

2857

2858 rotmat = num.array(

2859 [[ci, -cb*si, -sb*si], [si, cb*ci, sb*ci], [0., sb, -cb]])

2860 return project(traces, rotmat, in_channels, out_channels)

2861

2862

2863def _decompose(a):

2864 '''

2865 Decompose matrix into independent submatrices.

2866 '''

2867

2868 def depends(iout, a):

2869 row = a[iout, :]

2870 return set(num.arange(row.size).compress(row != 0.0))

2871

2872 def provides(iin, a):

2873 col = a[:, iin]

2874 return set(num.arange(col.size).compress(col != 0.0))

2875

2876 a = num.asarray(a)

2877 outs = set(range(a.shape[0]))

2878 systems = []

2879 while outs:

2880 iout = outs.pop()

2881

2882 gout = set()

2883 for iin in depends(iout, a):

2884 gout.update(provides(iin, a))

2885

2886 if not gout:

2887 continue

2888

2889 gin = set()

2890 for iout2 in gout:

2891 gin.update(depends(iout2, a))

2892

2893 if not gin:

2894 continue

2895

2896 for iout2 in gout:

2897 if iout2 in outs:

2898 outs.remove(iout2)

2899

2900 gin = list(gin)

2901 gin.sort()

2902 gout = list(gout)

2903 gout.sort()

2904

2905 systems.append((gin, gout, a[gout, :][:, gin]))

2906

2907 return systems

2908

2909

2910def _channels_to_names(channels):

2911 from pyrocko import squirrel

2912 names = []

2913 for ch in channels:

2914 if isinstance(ch, model.Channel):

2915 names.append(ch.name)

2916 elif isinstance(ch, squirrel.Channel):

2917 names.append(ch.codes.channel)

2918 else:

2919 names.append(ch)

2920

2921 return names

2922

2923

2924def project(traces, matrix, in_channels, out_channels):

2925 '''

2926 Affine transform of three-component traces.

2927

2928 Compute matrix-vector product of three-component traces, to e.g. rotate

2929 traces into a different basis. The traces are distinguished and ordered by

2930 their channel attribute. The tranform is applied to overlapping parts of

2931 any appropriate combinations of the input traces. This should allow this

2932 function to be robust with data gaps. It also tries to apply the

2933 tranformation to subsets of the channels, if this is possible, so that, if

2934 for example a vertical compontent is missing, horizontal components can

2935 still be rotated.

2936

2937 :param traces: list of traces in arbitrary order

2938 :param matrix: tranformation matrix

2939 :param in_channels: input channel names

2940 :param out_channels: output channel names

2941 :returns: list of transformed traces

2942 '''

2943

2944 in_channels = tuple(_channels_to_names(in_channels))

2945 out_channels = tuple(_channels_to_names(out_channels))

2946 systems = _decompose(matrix)

2947

2948 # fallback to full matrix if some are not quadratic

2949 for iins, iouts, submatrix in systems:

2950 if submatrix.shape[0] != submatrix.shape[1]:

2951 if len(in_channels) != 3 or len(out_channels) != 3:

2952 return []

2953 else:

2954 return _project3(traces, matrix, in_channels, out_channels)

2955

2956 projected = []

2957 for iins, iouts, submatrix in systems:

2958 in_cha = tuple([in_channels[iin] for iin in iins])

2959 out_cha = tuple([out_channels[iout] for iout in iouts])

2960 if submatrix.shape[0] == 1:

2961 projected.extend(_project1(traces, submatrix, in_cha, out_cha))

2962 elif submatrix.shape[1] == 2:

2963 projected.extend(_project2(traces, submatrix, in_cha, out_cha))

2964 else:

2965 projected.extend(_project3(traces, submatrix, in_cha, out_cha))

2966

2967 return projected

2968

2969

2970def project_dependencies(matrix, in_channels, out_channels):

2971 '''

2972 Figure out what dependencies project() would produce.

2973 '''

2974

2975 in_channels = tuple(_channels_to_names(in_channels))

2976 out_channels = tuple(_channels_to_names(out_channels))

2977 systems = _decompose(matrix)

2978

2979 subpro = []

2980 for iins, iouts, submatrix in systems:

2981 if submatrix.shape[0] != submatrix.shape[1]:

2982 subpro.append((matrix, in_channels, out_channels))

2983

2984 if not subpro:

2985 for iins, iouts, submatrix in systems:

2986 in_cha = tuple([in_channels[iin] for iin in iins])

2987 out_cha = tuple([out_channels[iout] for iout in iouts])

2988 subpro.append((submatrix, in_cha, out_cha))

2989

2990 deps = {}

2991 for mat, in_cha, out_cha in subpro:

2992 for oc in out_cha:

2993 if oc not in deps:

2994 deps[oc] = []

2995

2996 for ic in in_cha:

2997 deps[oc].append(ic)

2998

2999 return deps

3000

3001

3002def _project1(traces, matrix, in_channels, out_channels):

3003 assert len(in_channels) == 1

3004 assert len(out_channels) == 1

3005 assert matrix.shape == (1, 1)

3006

3007 projected = []

3008 for a in traces:

3009 if not (a.channel,) == in_channels:

3010 continue

3011

3012 ac = a.copy()

3013 ac.set_ydata(matrix[0, 0]*a.get_ydata())

3014 ac.set_codes(channel=out_channels[0])

3015 projected.append(ac)

3016

3017 return projected

3018

3019

3020def _project2(traces, matrix, in_channels, out_channels):

3021 assert len(in_channels) == 2

3022 assert len(out_channels) == 2

3023 assert matrix.shape == (2, 2)

3024 projected = []

3025 for a in traces:

3026 for b in traces:

3027 if not ((a.channel, b.channel) == in_channels and

3028 a.nslc_id[:3] == b.nslc_id[:3] and

3029 abs(a.deltat-b.deltat) < a.deltat*0.001):

3030 continue

3031

3032 tmin = max(a.tmin, b.tmin)

3033 tmax = min(a.tmax, b.tmax)

3034

3035 if tmin > tmax:

3036 continue

3037

3038 ac = a.chop(tmin, tmax, inplace=False, include_last=True)

3039 bc = b.chop(tmin, tmax, inplace=False, include_last=True)

3040 if abs(ac.tmin - bc.tmin) > ac.deltat*0.01:

3041 logger.warning(

3042 'Cannot project traces with displaced sampling '

3043 '(%s, %s, %s, %s)' % a.nslc_id)

3044 continue

3045

3046 acydata = num.dot(matrix[0], (ac.get_ydata(), bc.get_ydata()))

3047 bcydata = num.dot(matrix[1], (ac.get_ydata(), bc.get_ydata()))

3048

3049 ac.set_ydata(acydata)

3050 bc.set_ydata(bcydata)

3051

3052 ac.set_codes(channel=out_channels[0])

3053 bc.set_codes(channel=out_channels[1])

3054

3055 projected.append(ac)

3056 projected.append(bc)

3057

3058 return projected

3059

3060

3061def _project3(traces, matrix, in_channels, out_channels):

3062 assert len(in_channels) == 3

3063 assert len(out_channels) == 3

3064 assert matrix.shape == (3, 3)

3065 projected = []

3066 for a in traces:

3067 for b in traces:

3068 for c in traces:

3069 if not ((a.channel, b.channel, c.channel) == in_channels

3070 and a.nslc_id[:3] == b.nslc_id[:3]

3071 and b.nslc_id[:3] == c.nslc_id[:3]

3072 and abs(a.deltat-b.deltat) < a.deltat*0.001

3073 and abs(b.deltat-c.deltat) < b.deltat*0.001):

3074

3075 continue

3076

3077 tmin = max(a.tmin, b.tmin, c.tmin)

3078 tmax = min(a.tmax, b.tmax, c.tmax)

3079

3080 if tmin >= tmax:

3081 continue

3082

3083 ac = a.chop(tmin, tmax, inplace=False, include_last=True)

3084 bc = b.chop(tmin, tmax, inplace=False, include_last=True)

3085 cc = c.chop(tmin, tmax, inplace=False, include_last=True)

3086 if (abs(ac.tmin - bc.tmin) > ac.deltat*0.01

3087 or abs(bc.tmin - cc.tmin) > bc.deltat*0.01):

3088

3089 logger.warning(

3090 'Cannot project traces with displaced sampling '

3091 '(%s, %s, %s, %s)' % a.nslc_id)

3092 continue

3093

3094 acydata = num.dot(

3095 matrix[0],

3096 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata()))

3097 bcydata = num.dot(

3098 matrix[1],

3099 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata()))

3100 ccydata = num.dot(

3101 matrix[2],

3102 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata()))

3103

3104 ac.set_ydata(acydata)

3105 bc.set_ydata(bcydata)

3106 cc.set_ydata(ccydata)

3107

3108 ac.set_codes(channel=out_channels[0])

3109 bc.set_codes(channel=out_channels[1])

3110 cc.set_codes(channel=out_channels[2])

3111

3112 projected.append(ac)

3113 projected.append(bc)

3114 projected.append(cc)

3115

3116 return projected

3117

3118

3119def correlate(a, b, mode='valid', normalization=None, use_fft=False):

3120 '''

3121 Cross correlation of two traces.

3122

3123 :param a,b: input traces

3124 :param mode: ``'valid'``, ``'full'``, or ``'same'``

3125 :param normalization: ``'normal'``, ``'gliding'``, or ``None``

3126 :param use_fft: bool, whether to do cross correlation in spectral domain

3127

3128 :returns: trace containing cross correlation coefficients

3129

3130 This function computes the cross correlation between two traces. It

3131 evaluates the discrete equivalent of

3132

3133 .. math::

3134

3135 c(t) = \\int_{-\\infty}^{\\infty} a^{\\ast}(\\tau) b(t+\\tau) d\\tau

3136

3137 where the star denotes complex conjugate. Note, that the arguments here are

3138 swapped when compared with the :py:func:`numpy.correlate` function,

3139 which is internally called. This function should be safe even with older

3140 versions of NumPy, where the correlate function has some problems.

3141

3142 A trace containing the cross correlation coefficients is returned. The time

3143 information of the output trace is set so that the returned cross

3144 correlation can be viewed directly as a function of time lag.

3145

3146 Example::

3147

3148 # align two traces a and b containing a time shifted similar signal:

3149 c = pyrocko.trace.correlate(a,b)

3150 t, coef = c.max() # get time and value of maximum

3151 b.shift(-t) # align b with a

3152

3153 '''

3154

3155 assert_same_sampling_rate(a, b)

3156

3157 ya, yb = a.ydata, b.ydata

3158

3159 # need reversed order here:

3160 yc = numpy_correlate_fixed(yb, ya, mode=mode, use_fft=use_fft)

3161 kmin, kmax = numpy_correlate_lag_range(yb, ya, mode=mode, use_fft=use_fft)

3162

3163 if normalization == 'normal':

3164 normfac = num.sqrt(num.sum(ya**2))*num.sqrt(num.sum(yb**2))

3165 yc = yc/normfac

3166

3167 elif normalization == 'gliding':

3168 if mode != 'valid':

3169 assert False, 'gliding normalization currently only available ' \

3170 'with "valid" mode.'

3171

3172 if ya.size < yb.size:

3173 yshort, ylong = ya, yb

3174 else:

3175 yshort, ylong = yb, ya

3176

3177 epsilon = 0.00001

3178 normfac_short = num.sqrt(num.sum(yshort**2))

3179 normfac = normfac_short * num.sqrt(

3180 moving_sum(ylong**2, yshort.size, mode='valid')) \

3181 + normfac_short*epsilon

3182

3183 if yb.size <= ya.size:

3184 normfac = normfac[::-1]

3185

3186 yc /= normfac

3187

3188 c = a.copy()

3189 c.set_ydata(yc)

3190 c.set_codes(*merge_codes(a, b, '~'))

3191 c.shift(-c.tmin + b.tmin-a.tmin + kmin * c.deltat)

3192

3193 return c

3194

3195

3196def deconvolve(

3197 a, b, waterlevel,

3198 tshift=0.,

3199 pad=0.5,

3200 fd_taper=None,

3201 pad_to_pow2=True):

3202

3203 same_sampling_rate(a, b)

3204 assert abs(a.tmin - b.tmin) < a.deltat * 0.001

3205 deltat = a.deltat

3206 npad = int(round(a.data_len()*pad + tshift / deltat))

3207

3208 ndata = max(a.data_len(), b.data_len())

3209 ndata_pad = ndata + npad

3210

3211 if pad_to_pow2:

3212 ntrans = nextpow2(ndata_pad)

3213 else:

3214 ntrans = ndata

3215

3216 aspec = num.fft.rfft(a.ydata, ntrans)

3217 bspec = num.fft.rfft(b.ydata, ntrans)

3218

3219 out = aspec * num.conj(bspec)

3220

3221 bautocorr = bspec*num.conj(bspec)

3222 denom = num.maximum(bautocorr, waterlevel * bautocorr.max())

3223

3224 out /= denom

3225 df = 1/(ntrans*deltat)

3226

3227 if fd_taper is not None:

3228 fd_taper(out, 0.0, df)

3229

3230 ydata = num.roll(num.fft.irfft(out), int(round(tshift/deltat)))

3231 c = a.copy(data=False)

3232 c.set_ydata(ydata[:ndata])

3233 c.set_codes(*merge_codes(a, b, '/'))

3234 return c

3235

3236

3237def assert_same_sampling_rate(a, b, eps=1.0e-6):

3238 assert same_sampling_rate(a, b, eps), \

3239 'Sampling rates differ: %g != %g' % (a.deltat, b.deltat)

3240

3241

3242def same_sampling_rate(a, b, eps=1.0e-6):

3243 '''

3244 Check if two traces have the same sampling rate.

3245

3246 :param a,b: input traces

3247 :param eps: relative tolerance

3248 '''

3249

3250 return abs(a.deltat - b.deltat) < (a.deltat + b.deltat)*eps

3251

3252

3253def fix_deltat_rounding_errors(deltat):

3254 '''

3255 Try to undo sampling rate rounding errors.

3256

3257 Fix rounding errors of sampling intervals when these are read from single

3258 precision floating point values.

3259

3260 Assumes that the true sampling rate or sampling interval was an integer

3261 value. No correction will be applied if this would change the sampling

3262 rate by more than 0.001%.

3263 '''

3264

3265 if deltat <= 1.0:

3266 deltat_new = 1.0 / round(1.0 / deltat)

3267 else:

3268 deltat_new = round(deltat)

3269

3270 if abs(deltat_new - deltat) / deltat > 1e-5:

3271 deltat_new = deltat

3272

3273 return deltat_new

3274

3275

3276def merge_codes(a, b, sep='-'):

3277 '''

3278 Merge network-station-location-channel codes of a pair of traces.

3279 '''

3280

3281 o = []

3282 for xa, xb in zip(a.nslc_id, b.nslc_id):

3283 if xa == xb:

3284 o.append(xa)

3285 else:

3286 o.append(sep.join((xa, xb)))

3287 return o

3288

3289

3290class Taper(Object):

3291 '''

3292 Base class for tapers.

3293

3294 Does nothing by default.

3295 '''

3296

3297 def __call__(self, y, x0, dx):

3298 pass

3299

3300

3301class CosTaper(Taper):

3302 '''

3303 Cosine Taper.

3304

3305 :param a: start of fading in

3306 :param b: end of fading in

3307 :param c: start of fading out

3308 :param d: end of fading out

3309 '''

3310

3311 a = Float.T()

3312 b = Float.T()

3313 c = Float.T()

3314 d = Float.T()

3315

3316 def __init__(self, a, b, c, d):

3317 Taper.__init__(self, a=a, b=b, c=c, d=d)

3318

3319 def __call__(self, y, x0, dx):

3320

3321 if y.dtype == num.dtype(float):

3322 _apply_costaper = signal_ext.apply_costaper

3323 else:

3324 _apply_costaper = apply_costaper

3325

3326 _apply_costaper(self.a, self.b, self.c, self.d, y, x0, dx)

3327

3328 def span(self, y, x0, dx):

3329 return span_costaper(self.a, self.b, self.c, self.d, y, x0, dx)

3330

3331 def time_span(self):

3332 return self.a, self.d

3333

3334

3335class CosFader(Taper):

3336 '''

3337 Cosine Fader.

3338

3339 :param xfade: fade in and fade out time in seconds (optional)

3340 :param xfrac: fade in and fade out as fraction between 0. and 1. (optional)

3341

3342 Only one argument can be set. The other should to be ``None``.

3343 '''

3344

3345 xfade = Float.T(optional=True)

3346 xfrac = Float.T(optional=True)

3347

3348 def __init__(self, xfade=None, xfrac=None):

3349 Taper.__init__(self, xfade=xfade, xfrac=xfrac)

3350 assert (xfade is None) != (xfrac is None)

3351 self._xfade = xfade

3352 self._xfrac = xfrac

3353

3354 def __call__(self, y, x0, dx):

3355

3356 xfade = self._xfade

3357

3358 xlen = (y.size - 1)*dx

3359 if xfade is None:

3360 xfade = xlen * self._xfrac

3361

3362 a = x0

3363 b = x0 + xfade

3364 c = x0 + xlen - xfade

3365 d = x0 + xlen

3366

3367 apply_costaper(a, b, c, d, y, x0, dx)

3368

3369 def span(self, y, x0, dx):

3370 return 0, y.size

3371

3372 def time_span(self):

3373 return None, None

3374

3375

3376def none_min(li):

3377 if None in li:

3378 return None

3379 else:

3380 return min(x for x in li if x is not None)

3381

3382

3383def none_max(li):

3384 if None in li:

3385 return None

3386 else:

3387 return max(x for x in li if x is not None)

3388

3389

3390class MultiplyTaper(Taper):

3391 '''

3392 Multiplication of several tapers.

3393 '''

3394

3395 tapers = List.T(Taper.T())

3396

3397 def __init__(self, tapers=None):

3398 if tapers is None:

3399 tapers = []

3400

3401 Taper.__init__(self, tapers=tapers)

3402

3403 def __call__(self, y, x0, dx):

3404 for taper in self.tapers:

3405 taper(y, x0, dx)

3406

3407 def span(self, y, x0, dx):

3408 spans = []

3409 for taper in self.tapers:

3410 spans.append(taper.span(y, x0, dx))

3411

3412 mins, maxs = list(zip(*spans))

3413 return min(mins), max(maxs)

3414

3415 def time_span(self):

3416 spans = []

3417 for taper in self.tapers:

3418 spans.append(taper.time_span())

3419

3420 mins, maxs = list(zip(*spans))

3421 return none_min(mins), none_max(maxs)

3422

3423

3424class GaussTaper(Taper):

3425 '''

3426 Frequency domain Gaussian filter.

3427 '''

3428

3429 alpha = Float.T()

3430

3431 def __init__(self, alpha):

3432 Taper.__init__(self, alpha=alpha)

3433 self._alpha = alpha

3434

3435 def __call__(self, y, x0, dx):

3436 f = x0 + num.arange(y.size)*dx

3437 y *= num.exp(-num.pi**2 / (self._alpha**2) * f**2)

3438

3439

3440cached_coefficients = {}

3441

3442

3443def _get_cached_filter_coeffs(order, corners, btype):

3444 ck = (order, tuple(corners), btype)

3445 if ck not in cached_coefficients:

3446 if len(corners) == 1:

3447 corners = corners[0]

3448

3449 cached_coefficients[ck] = signal.butter(

3450 order, corners, btype=btype)

3451

3452 return cached_coefficients[ck]

3453

3454

3455class _globals(object):

3456 _numpy_has_correlate_flip_bug = None

3457

3458

3459def _default_key(tr):

3460 return (tr.network, tr.station, tr.location, tr.channel)

3461

3462

3463def numpy_has_correlate_flip_bug():

3464 '''

3465 Check if NumPy's correlate function reveals old behaviour.

3466 '''

3467

3468 if _globals._numpy_has_correlate_flip_bug is None:

3469 a = num.array([0, 0, 1, 0, 0, 0, 0])

3470 b = num.array([0, 0, 0, 0, 1, 0, 0, 0])

3471 ab = num.correlate(a, b, mode='same')

3472 ba = num.correlate(b, a, mode='same')

3473 _globals._numpy_has_correlate_flip_bug = num.all(ab == ba)

3474

3475 return _globals._numpy_has_correlate_flip_bug

3476

3477

3478def numpy_correlate_fixed(a, b, mode='valid', use_fft=False):

3479 '''

3480 Call :py:func:`numpy.correlate` with fixes.

3481

3482 c[k] = sum_i a[i+k] * conj(b[i])

3483

3484 Note that the result produced by newer numpy.correlate is always flipped

3485 with respect to the formula given in its documentation (if ascending k

3486 assumed for the output).

3487 '''

3488

3489 if use_fft:

3490 if a.size < b.size:

3491 c = signal.fftconvolve(b[::-1], a, mode=mode)

3492 else:

3493 c = signal.fftconvolve(a, b[::-1], mode=mode)

3494 return c

3495

3496 else:

3497 buggy = numpy_has_correlate_flip_bug()

3498

3499 a = num.asarray(a)

3500 b = num.asarray(b)

3501

3502 if buggy:

3503 b = num.conj(b)

3504

3505 c = num.correlate(a, b, mode=mode)

3506

3507 if buggy and a.size < b.size:

3508 return c[::-1]

3509 else:

3510 return c

3511

3512

3513def numpy_correlate_emulate(a, b, mode='valid'):

3514 '''

3515 Slow version of :py:func:`numpy.correlate` for comparison.

3516 '''

3517

3518 a = num.asarray(a)

3519 b = num.asarray(b)

3520 kmin = -(b.size-1)

3521 klen = a.size-kmin

3522 kmin, kmax = numpy_correlate_lag_range(a, b, mode=mode)

3523 kmin = int(kmin)

3524 kmax = int(kmax)

3525 klen = kmax - kmin + 1

3526 c = num.zeros(klen, dtype=num.promote_types(b.dtype, a.dtype))

3527 for k in range(kmin, kmin+klen):

3528 imin = max(0, -k)

3529 ilen = min(b.size, a.size-k) - imin

3530 c[k-kmin] = num.sum(

3531 a[imin+k:imin+ilen+k] * num.conj(b[imin:imin+ilen]))

3532

3533 return c

3534

3535

3536def numpy_correlate_lag_range(a, b, mode='valid', use_fft=False):

3537 '''

3538 Get range of lags for which :py:func:`numpy.correlate` produces values.

3539 '''

3540

3541 a = num.asarray(a)

3542 b = num.asarray(b)

3543

3544 kmin = -(b.size-1)

3545 if mode == 'full':

3546 klen = a.size-kmin

3547 elif mode == 'same':

3548 klen = max(a.size, b.size)

3549 kmin += (a.size+b.size-1 - max(a.size, b.size)) // 2 + \

3550 int(not use_fft and a.size % 2 == 0 and b.size > a.size)

3551 elif mode == 'valid':

3552 klen = abs(a.size - b.size) + 1

3553 kmin += min(a.size, b.size) - 1

3554

3555 return kmin, kmin + klen - 1

3556

3557

3558def autocorr(x, nshifts):

3559 '''

3560 Compute biased estimate of the first autocorrelation coefficients.

3561

3562 :param x: input array

3563 :param nshifts: number of coefficients to calculate

3564 '''

3565

3566 mean = num.mean(x)

3567 std = num.std(x)

3568 n = x.size

3569 xdm = x - mean

3570 r = num.zeros(nshifts)

3571 for k in range(nshifts):

3572 r[k] = 1./((n-num.abs(k))*std) * num.sum(xdm[:n-k] * xdm[k:])

3573

3574 return r

3575

3576

3577def yulewalker(x, order):

3578 '''

3579 Compute autoregression coefficients using Yule-Walker method.

3580

3581 :param x: input array

3582 :param order: number of coefficients to produce

3583

3584 A biased estimate of the autocorrelation is used. The Yule-Walker equations

3585 are solved by :py:func:`numpy.linalg.inv` instead of Levinson-Durbin

3586 recursion which is normally used.

3587 '''

3588

3589 gamma = autocorr(x, order+1)

3590 d = gamma[1:1+order]

3591 a = num.zeros((order, order))

3592 gamma2 = num.concatenate((gamma[::-1], gamma[1:order]))

3593 for i in range(order):

3594 ioff = order-i

3595 a[i, :] = gamma2[ioff:ioff+order]

3596

3597 return num.dot(num.linalg.inv(a), -d)

3598

3599

3600def moving_avg(x, n):

3601 n = int(n)

3602 cx = x.cumsum()

3603 nn = len(x)

3604 y = num.zeros(nn, dtype=cx.dtype)

3605 y[n//2:n//2+(nn-n)] = (cx[n:]-cx[:-n])/n

3606 y[:n//2] = y[n//2]

3607 y[n//2+(nn-n):] = y[n//2+(nn-n)-1]

3608 return y

3609

3610

3611def moving_sum(x, n, mode='valid'):

3612 n = int(n)

3613 cx = x.cumsum()

3614 nn = len(x)

3615

3616 if mode == 'valid':

3617 if nn-n+1 <= 0:

3618 return num.zeros(0, dtype=cx.dtype)

3619 y = num.zeros(nn-n+1, dtype=cx.dtype)

3620 y[0] = cx[n-1]

3621 y[1:nn-n+1] = cx[n:nn]-cx[0:nn-n]

3622

3623 if mode == 'full':

3624 y = num.zeros(nn+n-1, dtype=cx.dtype)

3625 if n <= nn:

3626 y[0:n] = cx[0:n]

3627 y[n:nn] = cx[n:nn]-cx[0:nn-n]

3628 y[nn:nn+n-1] = cx[-1]-cx[nn-n:nn-1]

3629 else:

3630 y[0:nn] = cx[0:nn]

3631 y[nn:n] = cx[nn-1]

3632 y[n:nn+n-1] = cx[nn-1] - cx[0:nn-1]

3633

3634 if mode == 'same':

3635 n1 = (n-1)//2

3636 y = num.zeros(nn, dtype=cx.dtype)

3637 if n <= nn:

3638 y[0:n-n1] = cx[n1:n]

3639 y[n-n1:nn-n1] = cx[n:nn]-cx[0:nn-n]

3640 y[nn-n1:nn] = cx[nn-1] - cx[nn-n:nn-n+n1]

3641 else:

3642 y[0:max(0, nn-n1)] = cx[min(n1, nn):nn]

3643 y[max(nn-n1, 0):min(n-n1, nn)] = cx[nn-1]

3644 y[min(n-n1, nn):nn] = cx[nn-1] - cx[0:max(0, nn-(n-n1))]

3645

3646 return y

3647

3648

3649def nextpow2(i):

3650 return 2**int(math.ceil(math.log(i)/math.log(2.)))

3651

3652

3653def snapper_w_offset(nmax, offset, delta, snapfun=math.ceil):

3654 def snap(x):

3655 return max(0, min(int(snapfun((x-offset)/delta)), nmax))

3656 return snap

3657

3658

3659def snapper(nmax, delta, snapfun=math.ceil):

3660 def snap(x):

3661 return max(0, min(int(snapfun(x/delta)), nmax))

3662 return snap

3663

3664

3665def apply_costaper(a, b, c, d, y, x0, dx):

3666 abcd = num.array((a, b, c, d), dtype=float)

3667 ja, jb, jc, jd = num.clip(num.ceil((abcd-x0)/dx).astype(int), 0, y.size)

3668 y[:ja] = 0.

3669 y[ja:jb] *= 0.5 \

3670 - 0.5*num.cos((dx*num.arange(ja, jb)-(a-x0))/(b-a)*num.pi)

3671 y[jc:jd] *= 0.5 \

3672 + 0.5*num.cos((dx*num.arange(jc, jd)-(c-x0))/(d-c)*num.pi)

3673 y[jd:] = 0.

3674

3675

3676def span_costaper(a, b, c, d, y, x0, dx):

3677 hi = snapper_w_offset(y.size, x0, dx)

3678 return hi(a), hi(d) - hi(a)

3679

3680

3681def costaper(a, b, c, d, nfreqs, deltaf):

3682 hi = snapper(nfreqs, deltaf)

3683 tap = num.zeros(nfreqs)

3684 tap[hi(a):hi(b)] = 0.5 \

3685 - 0.5*num.cos((deltaf*num.arange(hi(a), hi(b))-a)/(b-a)*num.pi)

3686 tap[hi(b):hi(c)] = 1.

3687 tap[hi(c):hi(d)] = 0.5 \

3688 + 0.5*num.cos((deltaf*num.arange(hi(c), hi(d))-c)/(d-c)*num.pi)

3689

3690 return tap

3691

3692

3693def t2ind(t, tdelta, snap=round):

3694 return int(snap(t/tdelta))

3695

3696

3697def hilbert(x, N=None):

3698 '''

3699 Return the hilbert transform of x of length N.

3700

3701 (from scipy.signal, but changed to use fft and ifft from numpy.fft)

3702 '''

3703

3704 x = num.asarray(x)

3705 if N is None:

3706 N = len(x)

3707 if N <= 0:

3708 raise ValueError('N must be positive.')

3709 if num.iscomplexobj(x):

3710 logger.warning('imaginary part of x ignored.')

3711 x = num.real(x)

3712

3713 Xf = num.fft.fft(x, N, axis=0)

3714 h = num.zeros(N)

3715 if N % 2 == 0:

3716 h[0] = h[N//2] = 1

3717 h[1:N//2] = 2

3718 else:

3719 h[0] = 1

3720 h[1:(N+1)//2] = 2

3721

3722 if len(x.shape) > 1:

3723 h = h[:, num.newaxis]

3724 x = num.fft.ifft(Xf*h)

3725 return x

3726

3727

3728def near(a, b, eps):

3729 return abs(a-b) < eps

3730

3731

3732def coroutine(func):

3733 def wrapper(*args, **kwargs):

3734 gen = func(*args, **kwargs)

3735 next(gen)

3736 return gen

3737

3738 wrapper.__name__ = func.__name__

3739 wrapper.__dict__ = func.__dict__

3740 wrapper.__doc__ = func.__doc__

3741 return wrapper

3742

3743

3744class States(object):

3745 '''

3746 Utility to store channel-specific state in coroutines.

3747 '''

3748

3749 def __init__(self):

3750 self._states = {}

3751

3752 def get(self, tr):

3753 k = tr.nslc_id

3754 if k in self._states:

3755 tmin, deltat, dtype, value = self._states[k]

3756 if (near(tmin, tr.tmin, deltat/100.)

3757 and near(deltat, tr.deltat, deltat/10000.)

3758 and dtype == tr.ydata.dtype):

3759

3760 return value

3761

3762 return None

3763

3764 def set(self, tr, value):

3765 k = tr.nslc_id

3766 if k in self._states and self._states[k][-1] is not value:

3767 self.free(self._states[k][-1])

3768

3769 self._states[k] = (tr.tmax+tr.deltat, tr.deltat, tr.ydata.dtype, value)

3770

3771 def free(self, value):

3772 pass

3773

3774

3775@coroutine

3776def co_list_append(list):

3777 while True:

3778 list.append((yield))

3779

3780

3781class ScipyBug(Exception):

3782 pass

3783

3784

3785@coroutine

3786def co_lfilter(target, b, a):

3787 '''

3788 Successively filter broken continuous trace data (coroutine).

3789

3790 Create coroutine which takes :py:class:`Trace` objects, filters their data

3791 through :py:func:`scipy.signal.lfilter` and sends new :py:class:`Trace`

3792 objects containing the filtered data to target. This is useful, if one

3793 wants to filter a long continuous time series, which is split into many

3794 successive traces without producing filter artifacts at trace boundaries.

3795

3796 Filter states are kept *per channel*, specifically, for each (network,

3797 station, location, channel) combination occuring in the input traces, a

3798 separate state is created and maintained. This makes it possible to filter

3799 multichannel or multistation data with only one :py:func:`co_lfilter`

3800 instance.

3801

3802 Filter state is reset, when gaps occur.

3803

3804 Use it like this::

3805

3806 from pyrocko.trace import co_lfilter, co_list_append

3807

3808 filtered_traces = []

3809 pipe = co_lfilter(co_list_append(filtered_traces), a, b)

3810 for trace in traces:

3811 pipe.send(trace)

3812

3813 pipe.close()

3814

3815 '''

3816

3817 try:

3818 states = States()

3819 output = None

3820 while True:

3821 input = (yield)

3822

3823 zi = states.get(input)

3824 if zi is None:

3825 zi = num.zeros(max(len(a), len(b))-1, dtype=float)

3826

3827 output = input.copy(data=False)

3828 try:

3829 ydata, zf = signal.lfilter(b, a, input.get_ydata(), zi=zi)

3830 except ValueError:

3831 raise ScipyBug(

3832 'signal.lfilter failed: could be related to a bug '

3833 'in some older scipy versions, e.g. on opensuse42.1')

3834

3835 output.set_ydata(ydata)

3836 states.set(input, zf)

3837 target.send(output)

3838

3839 except GeneratorExit:

3840 target.close()

3841

3842

3843def co_antialias(target, q, n=None, ftype='fir'):

3844 b, a, n = util.decimate_coeffs(q, n, ftype)

3845 anti = co_lfilter(target, b, a)

3846 return anti

3847

3848

3849@coroutine

3850def co_dropsamples(target, q, nfir):

3851 try:

3852 states = States()

3853 while True:

3854 tr = (yield)

3855 newdeltat = q * tr.deltat

3856 ioffset = states.get(tr)

3857 if ioffset is None:

3858 # for fir filter, the first nfir samples are pulluted by

3859 # boundary effects; cut it off.

3860 # for iir this may be (much) more, we do not correct for that.

3861 # put sample instances to a time which is a multiple of the

3862 # new sampling interval.

3863 newtmin_want = math.ceil(

3864 (tr.tmin+(nfir+1)*tr.deltat) / newdeltat) * newdeltat \

3865 - (nfir/2*tr.deltat)

3866 ioffset = int(round((newtmin_want - tr.tmin)/tr.deltat))

3867 if ioffset < 0:

3868 ioffset = ioffset % q

3869

3870 newtmin_have = tr.tmin + ioffset * tr.deltat

3871 newtr = tr.copy(data=False)

3872 newtr.deltat = newdeltat

3873 # because the fir kernel shifts data by nfir/2 samples:

3874 newtr.tmin = newtmin_have - (nfir/2*tr.deltat)

3875 newtr.set_ydata(tr.get_ydata()[ioffset::q].copy())

3876 states.set(tr, (ioffset % q - tr.data_len() % q) % q)

3877 target.send(newtr)

3878

3879 except GeneratorExit:

3880 target.close()

3881

3882

3883def co_downsample(target, q, n=None, ftype='fir'):

3884 '''

3885 Successively downsample broken continuous trace data (coroutine).

3886

3887 Create coroutine which takes :py:class:`Trace` objects, downsamples their

3888 data and sends new :py:class:`Trace` objects containing the downsampled

3889 data to target. This is useful, if one wants to downsample a long

3890 continuous time series, which is split into many successive traces without

3891 producing filter artifacts and gaps at trace boundaries.

3892

3893 Filter states are kept *per channel*, specifically, for each (network,

3894 station, location, channel) combination occuring in the input traces, a

3895 separate state is created and maintained. This makes it possible to filter

3896 multichannel or multistation data with only one :py:func:`co_lfilter`

3897 instance.

3898

3899 Filter state is reset, when gaps occur. The sampling instances are choosen

3900 so that they occur at (or as close as possible) to even multiples of the

3901 sampling interval of the downsampled trace (based on system time).

3902 '''

3903

3904 b, a, n = util.decimate_coeffs(q, n, ftype)

3905 return co_antialias(co_dropsamples(target, q, n), q, n, ftype)

3906

3907

3908@coroutine

3909def co_downsample_to(target, deltat):

3910

3911 decimators = {}

3912 try:

3913 while True:

3914 tr = (yield)

3915 ratio = deltat / tr.deltat

3916 rratio = round(ratio)

3917 if abs(rratio - ratio)/ratio > 0.0001:

3918 raise util.UnavailableDecimation('ratio = %g' % ratio)

3919

3920 deci_seq = tuple(x for x in util.decitab(int(rratio)) if x != 1)

3921 if deci_seq not in decimators:

3922 pipe = target

3923 for q in deci_seq[::-1]:

3924 pipe = co_downsample(pipe, q)

3925

3926 decimators[deci_seq] = pipe

3927

3928 decimators[deci_seq].send(tr)

3929

3930 except GeneratorExit:

3931 for g in decimators.values():

3932 g.close()

3933

3934

3935class DomainChoice(StringChoice):

3936 choices = [

3937 'time_domain',

3938 'frequency_domain',

3939 'envelope',

3940 'absolute',

3941 'cc_max_norm']

3942

3943

3944class MisfitSetup(Object):

3945 '''

3946 Contains misfit setup to be used in :py:meth:`Trace.misfit`

3947

3948 :param description: Description of the setup

3949 :param norm: L-norm classifier

3950 :param taper: Object of :py:class:`Taper`

3951 :param filter: Object of :py:class:`~pyrocko.response.FrequencyResponse`

3952 :param domain: ['time_domain', 'frequency_domain', 'envelope', 'absolute',

3953 'cc_max_norm']

3954

3955 Can be dumped to a yaml file.

3956 '''

3957

3958 xmltagname = 'misfitsetup'

3959 description = String.T(optional=True)

3960 norm = Int.T(optional=False)

3961 taper = Taper.T(optional=False)

3962 filter = FrequencyResponse.T(optional=True)

3963 domain = DomainChoice.T(default='time_domain')

3964

3965

3966def equalize_sampling_rates(trace_1, trace_2):

3967 '''

3968 Equalize sampling rates of two traces (reduce higher sampling rate to

3969 lower).

3970

3971 :param trace_1: :py:class:`Trace` object

3972 :param trace_2: :py:class:`Trace` object

3973

3974 Returns a copy of the resampled trace if resampling is needed.

3975 '''

3976

3977 if same_sampling_rate(trace_1, trace_2):

3978 return trace_1, trace_2

3979

3980 if trace_1.deltat < trace_2.deltat:

3981 t1_out = trace_1.copy()

3982 t1_out.downsample_to(deltat=trace_2.deltat, snap=True)

3983 logger.debug('Trace downsampled (return copy of trace): %s'

3984 % '.'.join(t1_out.nslc_id))

3985 return t1_out, trace_2

3986

3987 elif trace_1.deltat > trace_2.deltat:

3988 t2_out = trace_2.copy()

3989 t2_out.downsample_to(deltat=trace_1.deltat, snap=True)

3990 logger.debug('Trace downsampled (return copy of trace): %s'

3991 % '.'.join(t2_out.nslc_id))

3992 return trace_1, t2_out

3993

3994

3995def Lx_norm(u, v, norm=2):

3996 '''

3997 Calculate the misfit denominator *m* and the normalization divisor *n*

3998 according to norm.

3999

4000 The normalization divisor *n* is calculated from ``v``.

4001

4002 :param u: :py:class:`numpy.ndarray`

4003 :param v: :py:class:`numpy.ndarray`

4004 :param norm: (default = 2)

4005

4006 ``u`` and ``v`` must be of same size.

4007 '''

4008

4009 if norm == 1:

4010 return (

4011 num.sum(num.abs(v-u)),

4012 num.sum(num.abs(v)))

4013

4014 elif norm == 2:

4015 return (

4016 num.sqrt(num.sum((v-u)**2)),

4017 num.sqrt(num.sum(v**2)))

4018

4019 else:

4020 return (

4021 num.power(num.sum(num.abs(num.power(v - u, norm))), 1./norm),

4022 num.power(num.sum(num.abs(num.power(v, norm))), 1./norm))

4023

4024

4025def do_downsample(tr, deltat):

4026 if abs(tr.deltat - deltat) / tr.deltat > 1e-6:

4027 tr = tr.copy()

4028 tr.downsample_to(deltat, snap=True, demean=False)

4029 else:

4030 if tr.tmin/tr.deltat > 1e-6 or tr.tmax/tr.deltat > 1e-6:

4031 tr = tr.copy()

4032 tr.snap()

4033 return tr

4034

4035

4036def do_extend(tr, tmin, tmax):

4037 if tmin < tr.tmin or tmax > tr.tmax:

4038 tr = tr.copy()

4039 tr.extend(tmin=tmin, tmax=tmax, fillmethod='repeat')

4040

4041 return tr

4042

4043

4044def do_pre_taper(tr, taper):

4045 return tr.taper(taper, inplace=False, chop=True)

4046

4047

4048def do_fft(tr, filter):

4049 if filter is None:

4050 return tr

4051 else:

4052 ndata = tr.ydata.size

4053 nfft = nextpow2(ndata)

4054 padded = num.zeros(nfft, dtype=float)

4055 padded[:ndata] = tr.ydata

4056 spectrum = num.fft.rfft(padded)

4057 df = 1.0 / (tr.deltat * nfft)

4058 frequencies = num.arange(spectrum.size)*df

4059 return [tr, frequencies, spectrum]

4060

4061

4062def do_filter(inp, filter):

4063 if filter is None:

4064 return inp

4065 else:

4066 tr, frequencies, spectrum = inp

4067 spectrum *= filter.evaluate(frequencies)

4068 return [tr, frequencies, spectrum]

4069

4070

4071def do_ifft(inp):

4072 if isinstance(inp, Trace):

4073 return inp

4074 else:

4075 tr, _, spectrum = inp

4076 ndata = tr.ydata.size

4077 tr = tr.copy(data=False)

4078 tr.set_ydata(num.fft.irfft(spectrum)[:ndata])

4079 return tr

4080

4081

4082def check_alignment(t1, t2):

4083 if abs(t1.tmin-t2.tmin) > t1.deltat * 1e-4 or \

4084 abs(t1.tmax - t2.tmax) > t1.deltat * 1e-4 or \

4085 t1.ydata.shape != t2.ydata.shape:

4086 raise MisalignedTraces(

4087 'Cannot calculate misfit of %s and %s due to misaligned '

4088 'traces.' % ('.'.join(t1.nslc_id), '.'.join(t2.nslc_id)))

4089

4090

4091def check_overlaps(traces_a, traces_b=None, message='Traces overlap.'):

4092 for ia, tr_a in enumerate(traces_a):

4093 for tr_b in traces_a[ia+1:] if traces_b is None else traces_b:

4094 if tr_a.nslc_id == tr_b.nslc_id and tr_a.overlaps(

4095 tr_b.tmin, tr_b.tmax):

4096 raise OverlappingTraces(

4097 message + '\n Trace 1: %s\n Trace 2: %s' % (

4098 tr_a.summary, tr_b.summary))