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

1745 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Basic signal processing for seismic waveforms. 

8''' 

9 

10import time 

11import math 

12import copy 

13import logging 

14import hashlib 

15from collections import defaultdict 

16 

17import numpy as num 

18from scipy import signal 

19 

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

21from pyrocko.util import reuse 

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

23 StringChoice, Timestamp 

24from pyrocko.guts_array import Array 

25from pyrocko.model import squirrel_content 

26 

27# backward compatibility 

28from pyrocko.util import UnavailableDecimation # noqa 

29from pyrocko.response import ( # noqa 

30 FrequencyResponse, Evalresp, InverseEvalresp, PoleZeroResponse, 

31 ButterworthResponse, SampledResponse, IntegrationResponse, 

32 DifferentiationResponse, MultiplyResponse) 

33 

34 

35guts_prefix = 'pf' 

36 

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

38 

39 

40g_tapered_coeffs_cache = {} 

41g_one_response = FrequencyResponse() 

42 

43 

44@squirrel_content 

45class Trace(Object): 

46 

47 ''' 

48 Create new trace object. 

49 

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

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

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

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

54 and channel code). 

55 

56 :param network: network code 

57 :param station: station code 

58 :param location: location code 

59 :param channel: channel code 

60 :param extra: extra code 

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

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

63 computed from length of ``ydata``) 

64 :param deltat: sampling interval in [s] 

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

66 ``tmax`` is not ``None``) 

67 :param mtime: optional modification time 

68 

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

70 library) 

71 

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

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

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

75 silently truncated when the trace is stored 

76 ''' 

77 

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

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

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

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

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

83 

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

85 tmax = Timestamp.T() 

86 

87 deltat = Float.T(default=1.0) 

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

89 

90 mtime = Timestamp.T(optional=True) 

91 

92 cached_frequencies = {} 

93 

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

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

96 meta=None, extra=''): 

97 

98 Object.__init__(self, init_props=False) 

99 

100 time_float = util.get_time_float() 

101 

102 if not isinstance(tmin, time_float): 

103 tmin = Trace.tmin.regularize_extra(tmin) 

104 

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

106 tmax = Trace.tmax.regularize_extra(tmax) 

107 

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

109 mtime = Trace.mtime.regularize_extra(mtime) 

110 

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

112 ydata = Trace.ydata.regularize_extra(ydata) 

113 

114 self._growbuffer = None 

115 

116 tmin = time_float(tmin) 

117 if tmax is not None: 

118 tmax = time_float(tmax) 

119 

120 if mtime is None: 

121 mtime = time_float(time.time()) 

122 

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

124 self.extra = [ 

125 reuse(x) for x in ( 

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

127 

128 self.tmin = tmin 

129 self.deltat = deltat 

130 

131 if tmax is None: 

132 if ydata is not None: 

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

134 else: 

135 raise Exception( 

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

137 else: 

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

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

140 

141 self.meta = meta 

142 self.ydata = ydata 

143 self.mtime = mtime 

144 self._update_ids() 

145 self.file = None 

146 self._pchain = None 

147 

148 def __str__(self): 

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

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

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

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

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

154 

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

156 if self.meta: 

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

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

159 return s 

160 

161 @property 

162 def codes(self): 

163 from pyrocko.squirrel import model 

164 return model.CodesNSLCE( 

165 self.network, self.station, self.location, self.channel, 

166 self.extra) 

167 

168 @property 

169 def time_span(self): 

170 return self.tmin, self.tmax 

171 

172 @property 

173 def summary_entries(self): 

174 return ( 

175 self.__class__.__name__, 

176 str(self.codes), 

177 self.str_time_span, 

178 '%g' % (1.0/self.deltat)) 

179 

180 @property 

181 def summary_stats_entries(self): 

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

183 self.ydata.min(), 

184 self.ydata.max(), 

185 self.ydata.mean(), 

186 self.ydata.std())) 

187 

188 @property 

189 def summary(self): 

190 return util.fmt_summary(self.summary_entries, (10, 20, 55, 0)) 

191 

192 @property 

193 def summary_stats(self): 

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

195 self.summary_stats_entries, (12, 12, 12, 12)) 

196 

197 def __getstate__(self): 

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

199 self.tmin, self.tmax, self.deltat, self.mtime, 

200 self.ydata, self.meta, self.extra) 

201 

202 def __setstate__(self, state): 

203 if len(state) == 11: 

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

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

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

207 

208 elif len(state) == 12: 

209 # backward compatibility 0 

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

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

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

213 

214 elif len(state) == 10: 

215 # backward compatibility 1 

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

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

218 self.ydata, self.meta = state 

219 

220 self.extra = '' 

221 

222 else: 

223 # backward compatibility 2 

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

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

226 self.ydata = None 

227 self.meta = None 

228 self.extra = '' 

229 

230 self._growbuffer = None 

231 self._update_ids() 

232 

233 def hash(self, unsafe=False): 

234 sha1 = hashlib.sha1() 

235 for code in self.nslc_id: 

236 sha1.update(code.encode()) 

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

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

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

240 

241 if self.ydata is not None: 

242 if not unsafe: 

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

244 else: 

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

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

247 

248 return sha1.hexdigest() 

249 

250 def name(self): 

251 ''' 

252 Get a short string description. 

253 ''' 

254 

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

256 util.time_to_str(self.tmin), 

257 util.time_to_str(self.tmax))) 

258 

259 return s 

260 

261 def __eq__(self, other): 

262 return ( 

263 isinstance(other, Trace) 

264 and self.network == other.network 

265 and self.station == other.station 

266 and self.location == other.location 

267 and self.channel == other.channel 

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

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

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

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

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

273 

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

275 return ( 

276 self.network == other.network 

277 and self.station == other.station 

278 and self.location == other.location 

279 and self.channel == other.channel 

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

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

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

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

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

285 

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

287 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

304 

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

306 'trace values differ' 

307 

308 def __hash__(self): 

309 return id(self) 

310 

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

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

313 if clip: 

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

315 else: 

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

317 raise IndexError() 

318 

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

320 

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

322 ''' 

323 Value of trace between supporting points through linear interpolation. 

324 

325 :param t: time instant 

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

327 ''' 

328 

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

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

331 if t0 == t1: 

332 return y0 

333 else: 

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

335 

336 def index_clip(self, i): 

337 ''' 

338 Clip index to valid range. 

339 ''' 

340 

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

342 

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

344 ''' 

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

346 

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

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

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

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

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

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

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

354 match. 

355 ''' 

356 

357 if interpolate: 

358 assert self.deltat <= other.deltat \ 

359 or same_sampling_rate(self, other) 

360 

361 other_xdata = other.get_xdata() 

362 xdata = self.get_xdata() 

363 self.ydata += num.interp( 

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

365 else: 

366 assert self.deltat == other.deltat 

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

368 ibeg = max(0, ioff) 

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

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

371 

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

373 ''' 

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

375 

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

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

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

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

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

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

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

383 match. 

384 ''' 

385 

386 if interpolate: 

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

388 same_sampling_rate(self, other) 

389 

390 other_xdata = other.get_xdata() 

391 xdata = self.get_xdata() 

392 self.ydata *= num.interp( 

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

394 else: 

395 assert self.deltat == other.deltat 

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

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

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

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

400 

401 ibeg1 = self.index_clip(ibeg1) 

402 iend1 = self.index_clip(iend1) 

403 ibeg2 = self.index_clip(ibeg2) 

404 iend2 = self.index_clip(iend2) 

405 

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

407 

408 def max(self): 

409 ''' 

410 Get time and value of data maximum. 

411 ''' 

412 

413 i = num.argmax(self.ydata) 

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

415 

416 def min(self): 

417 ''' 

418 Get time and value of data minimum. 

419 ''' 

420 

421 i = num.argmin(self.ydata) 

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

423 

424 def absmax(self): 

425 ''' 

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

427 ''' 

428 

429 tmi, mi = self.min() 

430 tma, ma = self.max() 

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

432 return tmi, abs(mi) 

433 else: 

434 return tma, abs(ma) 

435 

436 def set_codes( 

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

438 extra=None): 

439 

440 ''' 

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

442 ''' 

443 

444 if network is not None: 

445 self.network = network 

446 if station is not None: 

447 self.station = station 

448 if location is not None: 

449 self.location = location 

450 if channel is not None: 

451 self.channel = channel 

452 if extra is not None: 

453 self.extra = extra 

454 

455 self._update_ids() 

456 

457 def set_network(self, network): 

458 self.network = network 

459 self._update_ids() 

460 

461 def set_station(self, station): 

462 self.station = station 

463 self._update_ids() 

464 

465 def set_location(self, location): 

466 self.location = location 

467 self._update_ids() 

468 

469 def set_channel(self, channel): 

470 self.channel = channel 

471 self._update_ids() 

472 

473 def overlaps(self, tmin, tmax): 

474 ''' 

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

476 ''' 

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

478 

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

480 ''' 

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

482 condition callback. (internal use) 

483 ''' 

484 

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

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

487 

488 def _update_ids(self): 

489 ''' 

490 Update dependent ids. 

491 ''' 

492 

493 self.full_id = ( 

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

495 self.nslc_id = reuse( 

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

497 

498 def prune_from_reuse_cache(self): 

499 util.deuse(self.nslc_id) 

500 util.deuse(self.network) 

501 util.deuse(self.station) 

502 util.deuse(self.location) 

503 util.deuse(self.channel) 

504 

505 def set_mtime(self, mtime): 

506 ''' 

507 Set modification time of the trace. 

508 ''' 

509 

510 self.mtime = mtime 

511 

512 def get_xdata(self): 

513 ''' 

514 Create array for time axis. 

515 ''' 

516 

517 if self.ydata is None: 

518 raise NoData() 

519 

520 return self.tmin \ 

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

522 

523 def get_ydata(self): 

524 ''' 

525 Get data array. 

526 ''' 

527 

528 if self.ydata is None: 

529 raise NoData() 

530 

531 return self.ydata 

532 

533 def set_ydata(self, new_ydata): 

534 ''' 

535 Replace data array. 

536 ''' 

537 

538 self.drop_growbuffer() 

539 self.ydata = new_ydata 

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

541 

542 def data_len(self): 

543 if self.ydata is not None: 

544 return self.ydata.size 

545 else: 

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

547 

548 def drop_data(self): 

549 ''' 

550 Forget data, make dataless trace. 

551 ''' 

552 

553 self.drop_growbuffer() 

554 self.ydata = None 

555 

556 def drop_growbuffer(self): 

557 ''' 

558 Detach the traces grow buffer. 

559 ''' 

560 

561 self._growbuffer = None 

562 self._pchain = None 

563 

564 def copy(self, data=True): 

565 ''' 

566 Make a deep copy of the trace. 

567 ''' 

568 

569 tracecopy = copy.copy(self) 

570 tracecopy.drop_growbuffer() 

571 if data: 

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

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

574 return tracecopy 

575 

576 def crop_zeros(self): 

577 ''' 

578 Remove any zeros at beginning and end. 

579 ''' 

580 

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

582 if indices.size == 0: 

583 raise NoData() 

584 

585 ibeg = indices[0] 

586 iend = indices[-1]+1 

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

588 return 

589 

590 self.drop_growbuffer() 

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

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

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

594 self._update_ids() 

595 

596 def append(self, data): 

597 ''' 

598 Append data to the end of the trace. 

599 

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

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

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

603 currently filled portion of the grow buffer array. 

604 ''' 

605 

606 assert self.ydata.dtype == data.dtype 

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

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

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

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

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

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

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

614 

615 def chop( 

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

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

618 

619 ''' 

620 Cut the trace to given time span. 

621 

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

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

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

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

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

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

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

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

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

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

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

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

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

635 span. 

636 ''' 

637 

638 if want_incomplete: 

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

640 raise NoData() 

641 else: 

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

643 raise NoData() 

644 

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

646 iplus = 0 

647 if include_last: 

648 iplus = 1 

649 

650 iend = min( 

651 self.data_len(), 

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

653 

654 if ibeg >= iend: 

655 raise NoData() 

656 

657 obj = self 

658 if not inplace: 

659 obj = self.copy(data=False) 

660 

661 self.drop_growbuffer() 

662 if self.ydata is not None: 

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

664 else: 

665 obj.ydata = None 

666 

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

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

669 

670 obj._update_ids() 

671 

672 return obj 

673 

674 def downsample( 

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

676 cut=False): 

677 

678 ''' 

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

680 

681 Antialiasing filter details: 

682 

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

684 ripple of 0.05 dB in the passband. 

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

686 well-behaved phase response. 

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

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

689 

690 Comparison of the digital filters: 

691 

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

693 :width: 60% 

694 :alt: Comparison of the downsampling filters. 

695 

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

697 

698 :param ndecimate: 

699 Decimation factor, avoid values larger than 8. 

700 :type ndecimate: 

701 int 

702 

703 :param snap: 

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

705 the sampling rate (according to absolute time). 

706 :type snap: 

707 bool 

708 

709 :param demean: 

710 Whether to demean the signal before filtering. 

711 :type demean: 

712 bool 

713 

714 :param ftype: 

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

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

717 

718 :param cut: 

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

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

721 :type cut: 

722 bool 

723 ''' 

724 newdeltat = self.deltat*ndecimate 

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

726 if snap: 

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

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

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

730 else: 

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

732 

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

734 if data.size != 0: 

735 if demean: 

736 data -= num.mean(data) 

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

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

739 else: 

740 self.ydata = data 

741 

742 self.tmin += ilag * self.deltat 

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

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

745 self._update_ids() 

746 

747 def downsample_to( 

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

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

750 

751 ''' 

752 Downsample to given sampling rate. 

753 

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

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

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

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

758 possible downsampling ratios. 

759 

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

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

762 

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

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

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

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

767 given downsampling configuration can be estimated with 

768 :py:func:`downsample_tpad`. 

769 

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

771 

772 :param deltat: 

773 Desired sampling interval in [s]. 

774 :type deltat: 

775 float 

776 

777 :param allow_upsample_max: 

778 Maximum allowed upsampling factor of the signal to achieve the 

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

780 :type allow_upsample_max: 

781 int 

782 

783 :param snap: 

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

785 the sampling rate (according to absolute time). 

786 :type snap: 

787 bool 

788 

789 :param demean: 

790 Whether to demean the signal before filtering. 

791 :type demean: 

792 bool 

793 

794 :param ftype: 

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

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

797 

798 :param cut: 

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

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

801 :type demean: 

802 bool 

803 ''' 

804 

805 upsratio, deci_seq = _configure_downsampling( 

806 self.deltat, deltat, allow_upsample_max) 

807 

808 if demean: 

809 self.drop_growbuffer() 

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

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

812 

813 if upsratio > 1: 

814 self.drop_growbuffer() 

815 ydata = self.ydata 

816 self.ydata = num.zeros( 

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

818 self.ydata[::upsratio] = ydata 

819 for i in range(1, upsratio): 

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

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

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

823 self.deltat = self.deltat/upsratio 

824 

825 for i, ndecimate in enumerate(deci_seq): 

826 self.downsample( 

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

828 

829 def resample(self, deltat): 

830 ''' 

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

832 

833 Resampling is performed in the frequency domain. 

834 ''' 

835 

836 ndata = self.ydata.size 

837 ntrans = nextpow2(ndata) 

838 fntrans2 = ntrans * self.deltat/deltat 

839 ntrans2 = int(round(fntrans2)) 

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

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

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

843 logger.warning( 

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

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

846 

847 data = self.ydata 

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

849 data_pad[:ndata] = data 

850 fdata = num.fft.rfft(data_pad) 

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

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

853 fdata2[:n] = fdata[:n] 

854 data2 = num.fft.irfft(fdata2) 

855 data2 = data2[:ndata2] 

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

857 self.deltat = deltat2 

858 self.set_ydata(data2) 

859 

860 def resample_simple(self, deltat): 

861 tyear = 3600*24*365. 

862 

863 if deltat == self.deltat: 

864 return 

865 

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

867 logger.warning( 

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

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

870 return 

871 

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

873 if abs(ninterval) < 20: 

874 logger.error( 

875 'resample_simple: sample insertion/deletion interval less ' 

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

877 raise ResamplingFailed() 

878 

879 delete = False 

880 if ninterval < 0: 

881 ninterval = - ninterval 

882 delete = True 

883 

884 tyearbegin = util.year_start(self.tmin) 

885 

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

887 

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

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

890 if nindices > 0: 

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

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

893 data = [] 

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

895 if delete: 

896 ln = ln[:-1] 

897 

898 data.append(ln) 

899 if not delete: 

900 if ln.size == 0: 

901 v = h[0] 

902 else: 

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

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

905 

906 data.append(data_split[-1]) 

907 

908 ydata_new = num.concatenate(data) 

909 

910 self.tmin = tyearbegin + nmin * deltat 

911 self.deltat = deltat 

912 self.set_ydata(ydata_new) 

913 

914 def stretch(self, tmin_new, tmax_new): 

915 ''' 

916 Stretch signal while preserving sample rate using sinc interpolation. 

917 

918 :param tmin_new: new time of first sample 

919 :param tmax_new: new time of last sample 

920 

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

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

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

924 that by the approximations used. 

925 ''' 

926 

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

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

929 

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

931 n_new = int(round(r)) 

932 if abs(n_new - r) > 0.001: 

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

934 

935 assert n_new >= 2 

936 

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

938 

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

940 signal_ext.antidrift(i_control, t_control, 

941 self.ydata.astype(float), 

942 tmin_new, self.deltat, ydata_new) 

943 

944 self.tmin = tmin_new 

945 self.set_ydata(ydata_new) 

946 self._update_ids() 

947 

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

949 raise_exception=False): 

950 

951 ''' 

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

953 

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

955 :param warn: whether to emit a warning 

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

957 exception. 

958 ''' 

959 

960 if frequency >= 0.5/self.deltat: 

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

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

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

964 if warn: 

965 logger.warning(message) 

966 if raise_exception: 

967 raise AboveNyquist(message) 

968 

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

970 nyquist_exception=False, demean=True): 

971 

972 ''' 

973 Apply Butterworth lowpass to the trace. 

974 

975 :param order: order of the filter 

976 :param corner: corner frequency of the filter 

977 

978 Mean is removed before filtering. 

979 ''' 

980 

981 self.nyquist_check( 

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

983 nyquist_exception) 

984 

985 (b, a) = _get_cached_filter_coeffs( 

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

987 

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

989 logger.warning( 

990 'Erroneous filter coefficients returned by ' 

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

992 'signal before filtering.') 

993 

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

995 if demean: 

996 data -= num.mean(data) 

997 self.drop_growbuffer() 

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

999 

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

1001 nyquist_exception=False, demean=True): 

1002 

1003 ''' 

1004 Apply butterworth highpass to the trace. 

1005 

1006 :param order: order of the filter 

1007 :param corner: corner frequency of the filter 

1008 

1009 Mean is removed before filtering. 

1010 ''' 

1011 

1012 self.nyquist_check( 

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

1014 nyquist_exception) 

1015 

1016 (b, a) = _get_cached_filter_coeffs( 

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

1018 

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

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

1021 logger.warning( 

1022 'Erroneous filter coefficients returned by ' 

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

1024 'signal before filtering.') 

1025 if demean: 

1026 data -= num.mean(data) 

1027 self.drop_growbuffer() 

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

1029 

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

1031 ''' 

1032 Apply butterworth bandpass to the trace. 

1033 

1034 :param order: order of the filter 

1035 :param corner_hp: lower corner frequency of the filter 

1036 :param corner_lp: upper corner frequency of the filter 

1037 

1038 Mean is removed before filtering. 

1039 ''' 

1040 

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

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

1043 (b, a) = _get_cached_filter_coeffs( 

1044 order, 

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

1046 btype='band') 

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

1048 if demean: 

1049 data -= num.mean(data) 

1050 self.drop_growbuffer() 

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

1052 

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

1054 ''' 

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

1056 

1057 :param order: order of the filter 

1058 :param corner_hp: lower corner frequency of the filter 

1059 :param corner_lp: upper corner frequency of the filter 

1060 

1061 Mean is removed before filtering. 

1062 ''' 

1063 

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

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

1066 (b, a) = _get_cached_filter_coeffs( 

1067 order, 

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

1069 btype='bandstop') 

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

1071 if demean: 

1072 data -= num.mean(data) 

1073 self.drop_growbuffer() 

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

1075 

1076 def envelope(self, inplace=True): 

1077 ''' 

1078 Calculate the envelope of the trace. 

1079 

1080 :param inplace: calculate envelope in place 

1081 

1082 The calculation follows: 

1083 

1084 .. math:: 

1085 

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

1087 

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

1089 ''' 

1090 

1091 y = self.ydata.astype(float) 

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

1093 if inplace: 

1094 self.drop_growbuffer() 

1095 self.ydata = env 

1096 else: 

1097 tr = self.copy(data=False) 

1098 tr.ydata = env 

1099 return tr 

1100 

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

1102 ''' 

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

1104 

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

1106 :param inplace: apply taper inplace 

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

1108 trace 

1109 ''' 

1110 

1111 if not inplace: 

1112 tr = self.copy() 

1113 else: 

1114 tr = self 

1115 

1116 if chop: 

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

1118 tr.shift(i*tr.deltat) 

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

1120 

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

1122 

1123 if not inplace: 

1124 return tr 

1125 

1126 def whiten(self, order=6): 

1127 ''' 

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

1129 

1130 :param order: order of the autoregression process 

1131 ''' 

1132 

1133 b, a = self.whitening_coefficients(order) 

1134 self.drop_growbuffer() 

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

1136 

1137 def whitening_coefficients(self, order=6): 

1138 ar = yulewalker(self.ydata, order) 

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

1140 return b, a 

1141 

1142 def ampspec_whiten( 

1143 self, 

1144 width, 

1145 td_taper='auto', 

1146 fd_taper='auto', 

1147 pad_to_pow2=True, 

1148 demean=True): 

1149 

1150 ''' 

1151 Whiten signal via frequency domain using moving average on amplitude 

1152 spectra. 

1153 

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

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

1156 ``None`` or ``'auto'``. 

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

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

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

1160 of 2^n 

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

1162 

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

1164 the spectrum is calculated and inversely weighted with a smoothed 

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

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

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

1168 time domain. 

1169 

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

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

1172 ''' 

1173 

1174 ndata = self.data_len() 

1175 

1176 if pad_to_pow2: 

1177 ntrans = nextpow2(ndata) 

1178 else: 

1179 ntrans = ndata 

1180 

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

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

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

1184 raise TraceTooShort( 

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

1186 

1187 if td_taper == 'auto': 

1188 td_taper = CosFader(1./width) 

1189 

1190 if fd_taper == 'auto': 

1191 fd_taper = CosFader(width) 

1192 

1193 if td_taper: 

1194 self.taper(td_taper) 

1195 

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

1197 if demean: 

1198 ydata -= ydata.mean() 

1199 

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

1201 

1202 amp = num.abs(spec) 

1203 nspec = amp.size 

1204 csamp = num.cumsum(amp) 

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

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

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

1208 amp_smoothed[:n1] = amp_smoothed[n1] 

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

1210 

1211 denom = amp_smoothed * amp 

1212 numer = amp 

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

1214 if eps == 0.0: 

1215 eps = 1e-9 

1216 

1217 numer += eps 

1218 denom += eps 

1219 spec *= numer/denom 

1220 

1221 if fd_taper: 

1222 fd_taper(spec, 0., df) 

1223 

1224 ydata = num.fft.irfft(spec) 

1225 self.set_ydata(ydata[:ndata]) 

1226 

1227 def _get_cached_freqs(self, nf, deltaf): 

1228 ck = (nf, deltaf) 

1229 if ck not in Trace.cached_frequencies: 

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

1231 nf, dtype=float) 

1232 

1233 return Trace.cached_frequencies[ck] 

1234 

1235 def bandpass_fft(self, corner_hp, corner_lp): 

1236 ''' 

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

1238 ''' 

1239 

1240 n = len(self.ydata) 

1241 n2 = nextpow2(n) 

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

1243 data[:n] = self.ydata 

1244 fdata = num.fft.rfft(data) 

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

1246 fdata[0] = 0.0 

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

1248 data = num.fft.irfft(fdata) 

1249 self.drop_growbuffer() 

1250 self.ydata = data[:n] 

1251 

1252 def shift(self, tshift): 

1253 ''' 

1254 Time shift the trace. 

1255 ''' 

1256 

1257 self.tmin += tshift 

1258 self.tmax += tshift 

1259 self._update_ids() 

1260 

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

1262 ''' 

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

1264 

1265 :param inplace: (boolean) snap traces inplace 

1266 

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

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

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

1270 ''' 

1271 

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

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

1274 

1275 if inplace: 

1276 xself = self 

1277 else: 

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

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

1280 return self 

1281 

1282 xself = self.copy() 

1283 

1284 if interpolate: 

1285 n = xself.data_len() 

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

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

1288 tref = tmin 

1289 t_control = num.array( 

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

1291 dtype=float) 

1292 

1293 signal_ext.antidrift(i_control, t_control, 

1294 xself.ydata.astype(float), 

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

1296 

1297 xself.ydata = ydata_new 

1298 

1299 xself.tmin = tmin 

1300 xself.tmax = tmax 

1301 xself._update_ids() 

1302 

1303 return xself 

1304 

1305 def fix_deltat_rounding_errors(self): 

1306 ''' 

1307 Try to undo sampling rate rounding errors. 

1308 

1309 See :py:func:`fix_deltat_rounding_errors`. 

1310 ''' 

1311 

1312 self.deltat = fix_deltat_rounding_errors(self.deltat) 

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

1314 

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

1316 ''' 

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

1318 the long time window. 

1319 

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

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

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

1323 filter 

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

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

1326 

1327 ============= ====================================== =========== 

1328 Scalingmethod Implementation Range 

1329 ============= ====================================== =========== 

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

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

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

1333 ============= ====================================== =========== 

1334 

1335 ''' 

1336 

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

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

1339 

1340 assert nshort < nlong 

1341 if nlong > len(self.ydata): 

1342 raise TraceTooShort( 

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

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

1345 

1346 if quad: 

1347 sqrdata = self.ydata**2 

1348 else: 

1349 sqrdata = self.ydata 

1350 

1351 mavg_short = moving_avg(sqrdata, nshort) 

1352 mavg_long = moving_avg(sqrdata, nlong) 

1353 

1354 self.drop_growbuffer() 

1355 

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

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

1358 

1359 if scalingmethod == 1: 

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

1361 elif scalingmethod in (2, 3): 

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

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

1364 

1365 if scalingmethod == 3: 

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

1367 

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

1369 ''' 

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

1371 with the last part of the long time window. 

1372 

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

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

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

1376 filter 

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

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

1379 

1380 ============= ====================================== =========== 

1381 Scalingmethod Implementation Range 

1382 ============= ====================================== =========== 

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

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

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

1386 ============= ====================================== =========== 

1387 

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

1389 STA/LTA are equivalent to 

1390 

1391 .. math:: 

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

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

1394 

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

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

1397 samples in the long time window. 

1398 ''' 

1399 

1400 n = self.data_len() 

1401 tmin = self.tmin 

1402 

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

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

1405 

1406 assert nshort < nlong 

1407 

1408 if nlong > len(self.ydata): 

1409 raise TraceTooShort( 

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

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

1412 

1413 if quad: 

1414 sqrdata = self.ydata**2 

1415 else: 

1416 sqrdata = self.ydata 

1417 

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

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

1420 nshift += 1 

1421 

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

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

1424 

1425 self.drop_growbuffer() 

1426 

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

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

1429 

1430 if scalingmethod == 1: 

1431 ydata = mavg_short/mavg_long * nshort/nlong 

1432 elif scalingmethod in (2, 3): 

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

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

1435 

1436 if scalingmethod == 3: 

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

1438 

1439 self.set_ydata(ydata) 

1440 

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

1442 

1443 self.chop( 

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

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

1446 

1447 def peaks(self, threshold, tsearch, 

1448 deadtime=False, 

1449 nblock_duration_detection=100): 

1450 

1451 ''' 

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

1453 

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

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

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

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

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

1459 ''' 

1460 

1461 y = self.ydata 

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

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

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

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

1466 tpeaks = [] 

1467 apeaks = [] 

1468 tzeros = [] 

1469 tzero = self.tmin 

1470 

1471 for itrig_pos in itrig_positions: 

1472 ibeg = itrig_pos 

1473 iend = min( 

1474 len(self.ydata), 

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

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

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

1478 apeak = y[ibeg+ipeak] 

1479 

1480 if tpeak < tzero: 

1481 continue 

1482 

1483 if deadtime: 

1484 ibeg = itrig_pos 

1485 iblock = 0 

1486 nblock = nblock_duration_detection 

1487 totalsum = 0. 

1488 while True: 

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

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

1491 break 

1492 

1493 logy = num.log( 

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

1495 logy[0] += totalsum 

1496 ysum = num.cumsum(logy) 

1497 totalsum = ysum[-1] 

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

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

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

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

1502 if len(izero_positions) > 0: 

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

1504 ibeg + izero_positions[0]) 

1505 break 

1506 iblock += 1 

1507 else: 

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

1509 

1510 tpeaks.append(tpeak) 

1511 apeaks.append(apeak) 

1512 tzeros.append(tzero) 

1513 

1514 if deadtime: 

1515 return tpeaks, apeaks, tzeros 

1516 else: 

1517 return tpeaks, apeaks 

1518 

1519 def peaks2(self, threshold, tsearch): 

1520 

1521 ''' 

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

1523 

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

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

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

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

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

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

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

1531 no more potential peaks are left. 

1532 ''' 

1533 

1534 a = num.copy(self.ydata) 

1535 

1536 amin = num.min(a) 

1537 

1538 a[0] = amin 

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

1540 a[-1] = amin 

1541 

1542 data = [] 

1543 while True: 

1544 imax = num.argmax(a) 

1545 amax = a[imax] 

1546 

1547 if amax < threshold or amax == amin: 

1548 break 

1549 

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

1551 

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

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

1554 

1555 if data: 

1556 data.sort() 

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

1558 else: 

1559 tpeaks, apeaks = [], [] 

1560 

1561 return tpeaks, apeaks 

1562 

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

1564 ''' 

1565 Extend trace to given span. 

1566 

1567 :param tmin: begin time of new span 

1568 :param tmax: end time of new span 

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

1570 ``'median'`` 

1571 ''' 

1572 

1573 nold = self.ydata.size 

1574 

1575 if tmin is not None: 

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

1577 else: 

1578 nl = 0 

1579 

1580 if tmax is not None: 

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

1582 else: 

1583 nh = nold - 1 

1584 

1585 n = nh - nl + 1 

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

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

1588 if self.ydata.size >= 1: 

1589 if fillmethod == 'repeat': 

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

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

1592 elif fillmethod == 'median': 

1593 v = num.median(self.ydata) 

1594 data[:-nl] = v 

1595 data[-nl + nold:] = v 

1596 elif fillmethod == 'mean': 

1597 v = num.mean(self.ydata) 

1598 data[:-nl] = v 

1599 data[-nl + nold:] = v 

1600 

1601 self.drop_growbuffer() 

1602 self.ydata = data 

1603 

1604 self.tmin += nl * self.deltat 

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

1606 

1607 self._update_ids() 

1608 

1609 def transfer(self, 

1610 tfade=0., 

1611 freqlimits=None, 

1612 transfer_function=None, 

1613 cut_off_fading=True, 

1614 demean=True, 

1615 invert=False): 

1616 

1617 ''' 

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

1619 

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

1621 at both ends of trace. 

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

1623 :param transfer_function: FrequencyResponse object; must provide a 

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

1625 coefficients at the frequencies 'freqs'. 

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

1627 trace. 

1628 :param demean: remove mean before applying transfer function 

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

1630 ''' 

1631 

1632 if transfer_function is None: 

1633 transfer_function = g_one_response 

1634 

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

1636 raise TraceTooShort( 

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

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

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

1640 

1641 if freqlimits is None and ( 

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

1643 

1644 # special case for flat responses 

1645 

1646 output = self.copy() 

1647 data = self.ydata 

1648 ndata = data.size 

1649 

1650 if transfer_function is not None: 

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

1652 

1653 if invert: 

1654 c = 1.0/c 

1655 

1656 data *= c 

1657 

1658 if tfade != 0.0: 

1659 data *= costaper( 

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

1661 ndata, self.deltat) 

1662 

1663 output.ydata = data 

1664 

1665 else: 

1666 ndata = self.ydata.size 

1667 ntrans = nextpow2(ndata*1.2) 

1668 coeffs = self._get_tapered_coeffs( 

1669 ntrans, freqlimits, transfer_function, invert=invert, 

1670 demean=demean) 

1671 

1672 data = self.ydata 

1673 

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

1675 data_pad[:ndata] = data 

1676 if demean: 

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

1678 

1679 if tfade != 0.0: 

1680 data_pad[:ndata] *= costaper( 

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

1682 ndata, self.deltat) 

1683 

1684 fdata = num.fft.rfft(data_pad) 

1685 fdata *= coeffs 

1686 ddata = num.fft.irfft(fdata) 

1687 output = self.copy() 

1688 output.ydata = ddata[:ndata] 

1689 

1690 if cut_off_fading and tfade != 0.0: 

1691 try: 

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

1693 except NoData: 

1694 raise TraceTooShort( 

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

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

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

1698 else: 

1699 output.ydata = output.ydata.copy() 

1700 

1701 return output 

1702 

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

1704 ''' 

1705 Approximate first or second derivative of the trace. 

1706 

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

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

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

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

1711 

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

1713 

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

1715 ''' 

1716 

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

1718 

1719 if inplace: 

1720 self.ydata = ddata 

1721 else: 

1722 output = self.copy(data=False) 

1723 output.set_ydata(ddata) 

1724 return output 

1725 

1726 def drop_chain_cache(self): 

1727 if self._pchain: 

1728 self._pchain.clear() 

1729 

1730 def init_chain(self): 

1731 self._pchain = pchain.Chain( 

1732 do_downsample, 

1733 do_extend, 

1734 do_pre_taper, 

1735 do_fft, 

1736 do_filter, 

1737 do_ifft) 

1738 

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

1740 if setup.domain == 'frequency_domain': 

1741 _, _, data = self._pchain( 

1742 (self, deltat), 

1743 (tmin, tmax), 

1744 (setup.taper,), 

1745 (setup.filter,), 

1746 (setup.filter,), 

1747 nocache=nocache) 

1748 

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

1750 

1751 else: 

1752 processed = self._pchain( 

1753 (self, deltat), 

1754 (tmin, tmax), 

1755 (setup.taper,), 

1756 (setup.filter,), 

1757 (setup.filter,), 

1758 (), 

1759 nocache=nocache) 

1760 

1761 if setup.domain == 'time_domain': 

1762 data = processed.get_ydata() 

1763 

1764 elif setup.domain == 'envelope': 

1765 processed = processed.envelope(inplace=False) 

1766 

1767 elif setup.domain == 'absolute': 

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

1769 

1770 return processed.get_ydata(), processed 

1771 

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

1773 ''' 

1774 Calculate misfit and normalization factor against candidate trace. 

1775 

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

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

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

1779 normalization divisor 

1780 

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

1782 with the higher sampling rate will be downsampled. 

1783 ''' 

1784 

1785 a = self 

1786 b = candidate 

1787 

1788 for tr in (a, b): 

1789 if not tr._pchain: 

1790 tr.init_chain() 

1791 

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

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

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

1795 

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

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

1798 

1799 if setup.domain != 'cc_max_norm': 

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

1801 else: 

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

1803 ccmax = ctr.max()[1] 

1804 m = 0.5 - 0.5 * ccmax 

1805 n = 0.5 

1806 

1807 if debug: 

1808 return m, n, aproc, bproc 

1809 else: 

1810 return m, n 

1811 

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

1813 ''' 

1814 Get FFT spectrum of trace. 

1815 

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

1817 power-of-two length 

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

1819 shaped tapers to both 

1820 

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

1822 ''' 

1823 

1824 ndata = self.ydata.size 

1825 

1826 if pad_to_pow2: 

1827 ntrans = nextpow2(ndata) 

1828 else: 

1829 ntrans = ndata 

1830 

1831 if tfade is None: 

1832 ydata = self.ydata 

1833 else: 

1834 ydata = self.ydata * costaper( 

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

1836 ndata, self.deltat) 

1837 

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

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

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

1841 return fxdata, fydata 

1842 

1843 def multi_filter(self, filter_freqs, bandwidth): 

1844 

1845 class Gauss(FrequencyResponse): 

1846 f0 = Float.T() 

1847 a = Float.T(default=1.0) 

1848 

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

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

1851 

1852 def evaluate(self, freqs): 

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

1854 omega = 2.*math.pi*freqs 

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

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

1857 

1858 freqs, coeffs = self.spectrum() 

1859 n = self.data_len() 

1860 nfilt = len(filter_freqs) 

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

1862 centroid_freqs = num.zeros(nfilt) 

1863 for ifilt, f0 in enumerate(filter_freqs): 

1864 taper = Gauss(f0, a=bandwidth) 

1865 weights = taper.evaluate(freqs) 

1866 nhalf = freqs.size 

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

1868 analytic_spec[:nhalf] = coeffs*weights 

1869 

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

1871 enorm /= num.sum(enorm) 

1872 

1873 if n % 2 == 0: 

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

1875 else: 

1876 analytic_spec[1:nhalf] *= 2. 

1877 

1878 analytic = num.fft.ifft(analytic_spec) 

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

1880 

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

1882 enorm /= num.sum(enorm) 

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

1884 

1885 return centroid_freqs, signal_tf 

1886 

1887 def _get_tapered_coeffs( 

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

1889 demean=True): 

1890 

1891 cache_key = ( 

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

1893 demean) 

1894 

1895 if cache_key in g_tapered_coeffs_cache: 

1896 return g_tapered_coeffs_cache[cache_key] 

1897 

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

1899 nfreqs = ntrans//2 + 1 

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

1901 hi = snapper(nfreqs, deltaf) 

1902 if freqlimits is not None: 

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

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

1905 coeffs = transfer_function.evaluate(freqs) 

1906 if invert: 

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

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

1909 

1910 transfer[kmin:kmax] = 1.0 / coeffs 

1911 else: 

1912 transfer[kmin:kmax] = coeffs 

1913 

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

1915 else: 

1916 if invert: 

1917 raise Exception( 

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

1919 'set to `True`') 

1920 

1921 freqs = num.arange(nfreqs) * deltaf 

1922 tapered_transfer = transfer_function.evaluate(freqs) 

1923 

1924 g_tapered_coeffs_cache[cache_key] = tapered_transfer 

1925 

1926 if demean: 

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

1928 

1929 return tapered_transfer 

1930 

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

1932 ''' 

1933 Fill string template with trace metadata. 

1934 

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

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

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

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

1939 ``tmin_year``, ``tmax_year``, ``tmin_month``, ``tmax_month``, 

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

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

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

1943 ''' 

1944 

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

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

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

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

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

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

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

1952 

1953 params = dict( 

1954 zip(('network', 'station', 'location', 'channel'), self.nslc_id)) 

1955 params['tmin'] = util.time_to_str( 

1956 self.tmin, format='%Y-%m-%d_%H-%M-%S') 

1957 params['tmax'] = util.time_to_str( 

1958 self.tmax, format='%Y-%m-%d_%H-%M-%S') 

1959 params['tmin_ms'] = util.time_to_str( 

1960 self.tmin, format='%Y-%m-%d_%H-%M-%S.3FRAC') 

1961 params['tmax_ms'] = util.time_to_str( 

1962 self.tmax, format='%Y-%m-%d_%H-%M-%S.3FRAC') 

1963 params['tmin_us'] = util.time_to_str( 

1964 self.tmin, format='%Y-%m-%d_%H-%M-%S.6FRAC') 

1965 params['tmax_us'] = util.time_to_str( 

1966 self.tmax, format='%Y-%m-%d_%H-%M-%S.6FRAC') 

1967 params['tmin_year'], params['tmin_month'], params['tmin_day'] \ 

1968 = util.time_to_str(self.tmin, format='%Y-%m-%d').split('-') 

1969 params['tmax_year'], params['tmax_month'], params['tmax_day'] \ 

1970 = util.time_to_str(self.tmax, format='%Y-%m-%d').split('-') 

1971 params['julianday'] = util.julian_day_of_year(self.tmin) 

1972 params.update(additional) 

1973 return template % params 

1974 

1975 def plot(self): 

1976 ''' 

1977 Show trace with matplotlib. 

1978 

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

1980 ''' 

1981 

1982 import pylab 

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

1984 name = '%s %s %s - %s' % ( 

1985 self.channel, 

1986 self.station, 

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

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

1989 

1990 pylab.title(name) 

1991 pylab.show() 

1992 

1993 def snuffle(self, **kwargs): 

1994 ''' 

1995 Show trace in a snuffler window. 

1996 

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

1998 objects or ``None`` 

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

2000 ``None`` 

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

2002 objects or ``None`` 

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

2004 12) 

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

2006 ``None`` 

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

2008 ``True``) 

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

2010 ''' 

2011 

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

2013 

2014 

2015def snuffle(traces, **kwargs): 

2016 ''' 

2017 Show traces in a snuffler window. 

2018 

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

2020 or ``None`` 

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

2022 ``None`` 

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

2024 objects or ``None`` 

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

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

2027 ``None`` 

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

2029 ``True``) 

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

2031 ''' 

2032 

2033 from pyrocko import pile 

2034 from pyrocko.gui.snuffler import snuffler 

2035 p = pile.Pile() 

2036 if traces: 

2037 trf = pile.MemTracesFile(None, traces) 

2038 p.add_file(trf) 

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

2040 

2041 

2042def downsample_tpad( 

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

2044 ''' 

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

2046 

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

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

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

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

2051 

2052 :param deltat_in: 

2053 Input sampling interval [s]. 

2054 :type deltat_in: 

2055 float 

2056 

2057 :param deltat_out: 

2058 Output samling interval [s]. 

2059 :type deltat_out: 

2060 float 

2061 

2062 :returns: 

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

2064 

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

2066 ''' 

2067 

2068 upsratio, deci_seq = _configure_downsampling( 

2069 deltat_in, deltat_out, allow_upsample_max) 

2070 

2071 tpad = 0.0 

2072 deltat = deltat_in / upsratio 

2073 for deci in deci_seq: 

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

2075 # n//2 for the antialiasing 

2076 # +deci for possible snap to multiples 

2077 # +1 for rounding errors 

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

2079 deltat = deltat * deci 

2080 

2081 return tpad 

2082 

2083 

2084def _configure_downsampling(deltat_in, deltat_out, allow_upsample_max): 

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

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

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

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

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

2090 

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

2092 

2093 

2094def _all_same(xs): 

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

2096 

2097 

2098def _incompatibilities(traces): 

2099 if not traces: 

2100 return None 

2101 

2102 params = [ 

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

2104 for tr in traces] 

2105 

2106 if not _all_same(params): 

2107 return params 

2108 else: 

2109 return None 

2110 

2111 

2112def _raise_incompatible_traces(params): 

2113 raise IncompatibleTraces( 

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

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

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

2117 'nsamples', 'dtype', 'deltat', 'tmin'), 

2118 '\n'.join( 

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

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

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

2122 

2123 

2124def _ensure_compatible(traces): 

2125 params = _incompatibilities(traces) 

2126 if params: 

2127 _raise_incompatible_traces(params) 

2128 

2129 

2130def _almost_equal(a, b, atol): 

2131 return abs(a-b) < atol 

2132 

2133 

2134def get_traces_data_as_array(traces): 

2135 ''' 

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

2137 

2138 :param traces: 

2139 Input waveforms. 

2140 :type traces: 

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

2142 

2143 :raises: 

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

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

2146 

2147 :returns: 

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

2149 :rtype: 

2150 :py:class:`numpy.ndarray` 

2151 ''' 

2152 

2153 if not traces: 

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

2155 

2156 _ensure_compatible(traces) 

2157 

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

2159 

2160 

2161def make_traces_compatible( 

2162 traces, 

2163 dtype=None, 

2164 deltat=None, 

2165 enforce_global_snap=True, 

2166 warn_snap=False): 

2167 

2168 ''' 

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

2170 

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

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

2173 time. 

2174 

2175 If necessary, traces are (in order): 

2176 

2177 - casted to the same data type. 

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

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

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

2181 used. 

2182 

2183 :param traces: 

2184 Input waveforms. 

2185 :type traces: 

2186 :py:class:`list` of :py:class:`Trace` 

2187 

2188 :param dtype: 

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

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

2191 :type dtype: 

2192 :py:class:`numpy.dtype` 

2193 

2194 :param deltat: 

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

2196 among the input traces is chosen. 

2197 :type deltat: 

2198 float 

2199 

2200 :param enforce_global_snap: 

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

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

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

2204 offset to the system time sampling rate multiples. 

2205 :type enforce_global_snap: 

2206 bool 

2207 

2208 :param warn_snap: 

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

2210 :type warn_snap: 

2211 bool 

2212 ''' 

2213 

2214 eps_snap = 1e-3 

2215 

2216 if not traces: 

2217 return [] 

2218 

2219 traces = list(traces) 

2220 

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

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

2223 

2224 if dtype is None: 

2225 dtype = float 

2226 logger.warning( 

2227 'make_traces_compatible: Inconsistent data types - converting ' 

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

2229 

2230 for itr, tr in enumerate(traces): 

2231 tr_copy = tr.copy(data=False) 

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

2233 traces[itr] = tr_copy 

2234 

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

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

2237 if deltat is None: 

2238 deltat = max(deltats) 

2239 logger.warning( 

2240 'make_traces_compatible: Inconsistent sampling rates - ' 

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

2242 % (1.0 / deltat)) 

2243 

2244 for itr, tr in enumerate(traces): 

2245 if tr.deltat != deltat: 

2246 tr_copy = tr.copy() 

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

2248 traces[itr] = tr_copy 

2249 

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

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

2252 > deltat * eps_snap 

2253 

2254 if enforce_global_snap or any(is_aligned): 

2255 tref = util.to_time_float(0.0) 

2256 else: 

2257 # to keep a common subsample shift 

2258 tref = num.max(tmins) 

2259 

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

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

2262 if num.any(need_snap): 

2263 if warn_snap: 

2264 logger.warning( 

2265 'make_traces_compatible: Misaligned sampling - introducing ' 

2266 'subsample shifts for proper alignment.') 

2267 

2268 for itr, tr in enumerate(traces): 

2269 if need_snap[itr]: 

2270 tr_copy = tr.copy() 

2271 if tref != 0.0: 

2272 tr_copy.shift(-tref) 

2273 

2274 tr_copy.snap(interpolate=True) 

2275 if tref != 0.0: 

2276 tr_copy.shift(tref) 

2277 

2278 traces[itr] = tr_copy 

2279 

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

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

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

2283 

2284 tmin = num.max(tmins) 

2285 tmax = num.min(tmaxs) 

2286 

2287 if tmin > tmax: 

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

2289 

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

2291 for itr, tr in enumerate(traces): 

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

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

2294 

2295 traces[itr] = tr.chop( 

2296 tmin, tmax, 

2297 inplace=False, 

2298 want_incomplete=False, 

2299 include_last=True) 

2300 

2301 xtr = traces[itr] 

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

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

2304 xtr.tmin = tmin 

2305 xtr.tmax = tmax 

2306 xtr.deltat = deltat 

2307 xtr._update_ids() 

2308 

2309 return traces 

2310 

2311 

2312class IncompatibleTraces(Exception): 

2313 ''' 

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

2315 ''' 

2316 

2317 

2318class InfiniteResponse(Exception): 

2319 ''' 

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

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

2322 result in a division by zero. 

2323 ''' 

2324 

2325 

2326class MisalignedTraces(Exception): 

2327 ''' 

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

2329 tmax or number of samples do not match. 

2330 ''' 

2331 

2332 pass 

2333 

2334 

2335class NoData(Exception): 

2336 ''' 

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

2338 not enough data is available. 

2339 ''' 

2340 

2341 pass 

2342 

2343 

2344class AboveNyquist(Exception): 

2345 ''' 

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

2347 frequencies are above the Nyquist frequency. 

2348 ''' 

2349 

2350 pass 

2351 

2352 

2353class TraceTooShort(Exception): 

2354 ''' 

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

2356 trace is too short. 

2357 ''' 

2358 

2359 pass 

2360 

2361 

2362class ResamplingFailed(Exception): 

2363 pass 

2364 

2365 

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

2367 

2368 ''' 

2369 Get data range given traces grouped by selected pattern. 

2370 

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

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

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

2374 used. 

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

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

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

2378 

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

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

2381 extreme values on either end. 

2382 

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

2384 

2385 Examples:: 

2386 

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

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

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

2390 

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

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

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

2394 

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

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

2397 ''' 

2398 

2399 if key is None: 

2400 key = _default_key 

2401 

2402 ranges = defaultdict(list) 

2403 for trace in traces: 

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

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

2406 else: 

2407 mean = trace.ydata.mean() 

2408 std = trace.ydata.std() 

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

2410 

2411 k = key(trace) 

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

2413 

2414 for k in ranges: 

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

2416 if outer_mode == 'minmax': 

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

2418 elif outer_mode == 'robust': 

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

2420 

2421 return ranges 

2422 

2423 

2424def minmaxtime(traces, key=None): 

2425 

2426 ''' 

2427 Get time range given traces grouped by selected pattern. 

2428 

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

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

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

2432 used. 

2433 

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

2435 ''' 

2436 

2437 if key is None: 

2438 key = _default_key 

2439 

2440 ranges = {} 

2441 for trace in traces: 

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

2443 k = key(trace) 

2444 if k not in ranges: 

2445 ranges[k] = mi, ma 

2446 else: 

2447 tmi, tma = ranges[k] 

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

2449 

2450 return ranges 

2451 

2452 

2453def degapper( 

2454 traces, 

2455 maxgap=5, 

2456 fillmethod='interpolate', 

2457 deoverlap='use_second', 

2458 maxlap=None): 

2459 

2460 ''' 

2461 Try to connect traces and remove gaps. 

2462 

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

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

2465 according to the ``deoverlap`` argument. 

2466 

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

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

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

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

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

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

2473 values. 

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

2475 

2476 :returns: list of traces 

2477 ''' 

2478 

2479 in_traces = traces 

2480 out_traces = [] 

2481 if not in_traces: 

2482 return out_traces 

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

2484 while in_traces: 

2485 

2486 a = out_traces[-1] 

2487 b = in_traces.pop(0) 

2488 

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

2490 assert avirt == bvirt, \ 

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

2492 'no data.' 

2493 

2494 virtual = avirt and bvirt 

2495 

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

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

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

2499 

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

2501 idist = int(round(dist)) 

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

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

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

2505 pass 

2506 else: 

2507 if 1 < idist <= maxgap: 

2508 if not virtual: 

2509 if fillmethod == 'interpolate': 

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

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

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

2513 ).astype(a.ydata.dtype) 

2514 elif fillmethod == 'zeros': 

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

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

2517 a.tmax = b.tmax 

2518 if a.mtime and b.mtime: 

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

2520 continue 

2521 

2522 elif idist == 1: 

2523 if not virtual: 

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

2525 a.tmax = b.tmax 

2526 if a.mtime and b.mtime: 

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

2528 continue 

2529 

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

2531 if b.tmax > a.tmax: 

2532 if not virtual: 

2533 na = a.ydata.size 

2534 n = -idist+1 

2535 if deoverlap == 'use_second': 

2536 a.ydata = num.concatenate( 

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

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

2539 a.ydata = num.concatenate( 

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

2541 elif deoverlap == 'add': 

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

2543 a.ydata = num.concatenate( 

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

2545 else: 

2546 assert False, 'unknown deoverlap method' 

2547 

2548 if deoverlap == 'crossfade_cos': 

2549 n = -idist+1 

2550 taper = 0.5-0.5*num.cos( 

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

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

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

2554 

2555 a.tmax = b.tmax 

2556 if a.mtime and b.mtime: 

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

2558 continue 

2559 else: 

2560 # make short second trace vanish 

2561 continue 

2562 

2563 if b.data_len() >= 1: 

2564 out_traces.append(b) 

2565 

2566 for tr in out_traces: 

2567 tr._update_ids() 

2568 

2569 return out_traces 

2570 

2571 

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

2573 ''' 

2574 2D rotation of traces. 

2575 

2576 :param traces: list of input traces 

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

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

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

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

2581 :returns: list of rotated traces 

2582 ''' 

2583 

2584 phi = azimuth/180.*math.pi 

2585 cphi = math.cos(phi) 

2586 sphi = math.sin(phi) 

2587 rotated = [] 

2588 in_channels = tuple(_channels_to_names(in_channels)) 

2589 out_channels = tuple(_channels_to_names(out_channels)) 

2590 for a in traces: 

2591 for b in traces: 

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

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

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

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

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

2597 

2598 if tmin < tmax: 

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

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

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

2602 logger.warning( 

2603 'Cannot rotate traces with displaced sampling ' 

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

2605 continue 

2606 

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

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

2609 ac.set_ydata(acydata) 

2610 bc.set_ydata(bcydata) 

2611 

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

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

2614 rotated.append(ac) 

2615 rotated.append(bc) 

2616 

2617 return rotated 

2618 

2619 

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

2621 ''' 

2622 Rotate traces from NE to RT system. 

2623 

2624 :param n: 

2625 North trace. 

2626 :type n: 

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

2628 

2629 :param e: 

2630 East trace. 

2631 :type e: 

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

2633 

2634 :param source: 

2635 Source of the recorded signal. 

2636 :type source: 

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

2638 

2639 :param receiver: 

2640 Receiver of the recorded signal. 

2641 :type receiver: 

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

2643 

2644 :param out_channels: 

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

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

2647 

2648 :type out_channels 

2649 optional, tuple[str, str] 

2650 

2651 :returns: 

2652 Rotated traces (radial, transversal). 

2653 :rtype: 

2654 tuple[ 

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

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

2657 ''' 

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

2659 in_channels = n.channel, e.channel 

2660 out = rotate( 

2661 [n, e], azimuth, 

2662 in_channels=in_channels, 

2663 out_channels=out_channels) 

2664 

2665 assert len(out) == 2 

2666 for tr in out: 

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

2668 r = tr 

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

2670 t = tr 

2671 else: 

2672 assert False 

2673 

2674 return r, t 

2675 

2676 

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

2678 out_channels=('L', 'Q', 'T')): 

2679 ''' 

2680 Rotate traces from ZNE to LQT system. 

2681 

2682 :param traces: list of traces in arbitrary order 

2683 :param backazimuth: backazimuth in degrees clockwise from north 

2684 :param incidence: incidence angle in degrees from vertical 

2685 :param in_channels: input channel names 

2686 :param out_channels: output channel names (default: ('L', 'Q', 'T')) 

2687 :returns: list of transformed traces 

2688 ''' 

2689 i = incidence/180.*num.pi 

2690 b = backazimuth/180.*num.pi 

2691 

2692 ci = num.cos(i) 

2693 cb = num.cos(b) 

2694 si = num.sin(i) 

2695 sb = num.sin(b) 

2696 

2697 rotmat = num.array( 

2698 [[ci, -cb*si, -sb*si], [si, cb*ci, sb*ci], [0., sb, -cb]]) 

2699 return project(traces, rotmat, in_channels, out_channels) 

2700 

2701 

2702def _decompose(a): 

2703 ''' 

2704 Decompose matrix into independent submatrices. 

2705 ''' 

2706 

2707 def depends(iout, a): 

2708 row = a[iout, :] 

2709 return set(num.arange(row.size).compress(row != 0.0)) 

2710 

2711 def provides(iin, a): 

2712 col = a[:, iin] 

2713 return set(num.arange(col.size).compress(col != 0.0)) 

2714 

2715 a = num.asarray(a) 

2716 outs = set(range(a.shape[0])) 

2717 systems = [] 

2718 while outs: 

2719 iout = outs.pop() 

2720 

2721 gout = set() 

2722 for iin in depends(iout, a): 

2723 gout.update(provides(iin, a)) 

2724 

2725 if not gout: 

2726 continue 

2727 

2728 gin = set() 

2729 for iout2 in gout: 

2730 gin.update(depends(iout2, a)) 

2731 

2732 if not gin: 

2733 continue 

2734 

2735 for iout2 in gout: 

2736 if iout2 in outs: 

2737 outs.remove(iout2) 

2738 

2739 gin = list(gin) 

2740 gin.sort() 

2741 gout = list(gout) 

2742 gout.sort() 

2743 

2744 systems.append((gin, gout, a[gout, :][:, gin])) 

2745 

2746 return systems 

2747 

2748 

2749def _channels_to_names(channels): 

2750 from pyrocko import squirrel 

2751 names = [] 

2752 for ch in channels: 

2753 if isinstance(ch, model.Channel): 

2754 names.append(ch.name) 

2755 elif isinstance(ch, squirrel.Channel): 

2756 names.append(ch.codes.channel) 

2757 else: 

2758 names.append(ch) 

2759 

2760 return names 

2761 

2762 

2763def project(traces, matrix, in_channels, out_channels): 

2764 ''' 

2765 Affine transform of three-component traces. 

2766 

2767 Compute matrix-vector product of three-component traces, to e.g. rotate 

2768 traces into a different basis. The traces are distinguished and ordered by 

2769 their channel attribute. The tranform is applied to overlapping parts of 

2770 any appropriate combinations of the input traces. This should allow this 

2771 function to be robust with data gaps. It also tries to apply the 

2772 tranformation to subsets of the channels, if this is possible, so that, if 

2773 for example a vertical compontent is missing, horizontal components can 

2774 still be rotated. 

2775 

2776 :param traces: list of traces in arbitrary order 

2777 :param matrix: tranformation matrix 

2778 :param in_channels: input channel names 

2779 :param out_channels: output channel names 

2780 :returns: list of transformed traces 

2781 ''' 

2782 

2783 in_channels = tuple(_channels_to_names(in_channels)) 

2784 out_channels = tuple(_channels_to_names(out_channels)) 

2785 systems = _decompose(matrix) 

2786 

2787 # fallback to full matrix if some are not quadratic 

2788 for iins, iouts, submatrix in systems: 

2789 if submatrix.shape[0] != submatrix.shape[1]: 

2790 if len(in_channels) != 3 or len(out_channels) != 3: 

2791 return [] 

2792 else: 

2793 return _project3(traces, matrix, in_channels, out_channels) 

2794 

2795 projected = [] 

2796 for iins, iouts, submatrix in systems: 

2797 in_cha = tuple([in_channels[iin] for iin in iins]) 

2798 out_cha = tuple([out_channels[iout] for iout in iouts]) 

2799 if submatrix.shape[0] == 1: 

2800 projected.extend(_project1(traces, submatrix, in_cha, out_cha)) 

2801 elif submatrix.shape[1] == 2: 

2802 projected.extend(_project2(traces, submatrix, in_cha, out_cha)) 

2803 else: 

2804 projected.extend(_project3(traces, submatrix, in_cha, out_cha)) 

2805 

2806 return projected 

2807 

2808 

2809def project_dependencies(matrix, in_channels, out_channels): 

2810 ''' 

2811 Figure out what dependencies project() would produce. 

2812 ''' 

2813 

2814 in_channels = tuple(_channels_to_names(in_channels)) 

2815 out_channels = tuple(_channels_to_names(out_channels)) 

2816 systems = _decompose(matrix) 

2817 

2818 subpro = [] 

2819 for iins, iouts, submatrix in systems: 

2820 if submatrix.shape[0] != submatrix.shape[1]: 

2821 subpro.append((matrix, in_channels, out_channels)) 

2822 

2823 if not subpro: 

2824 for iins, iouts, submatrix in systems: 

2825 in_cha = tuple([in_channels[iin] for iin in iins]) 

2826 out_cha = tuple([out_channels[iout] for iout in iouts]) 

2827 subpro.append((submatrix, in_cha, out_cha)) 

2828 

2829 deps = {} 

2830 for mat, in_cha, out_cha in subpro: 

2831 for oc in out_cha: 

2832 if oc not in deps: 

2833 deps[oc] = [] 

2834 

2835 for ic in in_cha: 

2836 deps[oc].append(ic) 

2837 

2838 return deps 

2839 

2840 

2841def _project1(traces, matrix, in_channels, out_channels): 

2842 assert len(in_channels) == 1 

2843 assert len(out_channels) == 1 

2844 assert matrix.shape == (1, 1) 

2845 

2846 projected = [] 

2847 for a in traces: 

2848 if not (a.channel,) == in_channels: 

2849 continue 

2850 

2851 ac = a.copy() 

2852 ac.set_ydata(matrix[0, 0]*a.get_ydata()) 

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

2854 projected.append(ac) 

2855 

2856 return projected 

2857 

2858 

2859def _project2(traces, matrix, in_channels, out_channels): 

2860 assert len(in_channels) == 2 

2861 assert len(out_channels) == 2 

2862 assert matrix.shape == (2, 2) 

2863 projected = [] 

2864 for a in traces: 

2865 for b in traces: 

2866 if not ((a.channel, b.channel) == in_channels and 

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

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

2869 continue 

2870 

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

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

2873 

2874 if tmin > tmax: 

2875 continue 

2876 

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

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

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

2880 logger.warning( 

2881 'Cannot project traces with displaced sampling ' 

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

2883 continue 

2884 

2885 acydata = num.dot(matrix[0], (ac.get_ydata(), bc.get_ydata())) 

2886 bcydata = num.dot(matrix[1], (ac.get_ydata(), bc.get_ydata())) 

2887 

2888 ac.set_ydata(acydata) 

2889 bc.set_ydata(bcydata) 

2890 

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

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

2893 

2894 projected.append(ac) 

2895 projected.append(bc) 

2896 

2897 return projected 

2898 

2899 

2900def _project3(traces, matrix, in_channels, out_channels): 

2901 assert len(in_channels) == 3 

2902 assert len(out_channels) == 3 

2903 assert matrix.shape == (3, 3) 

2904 projected = [] 

2905 for a in traces: 

2906 for b in traces: 

2907 for c in traces: 

2908 if not ((a.channel, b.channel, c.channel) == in_channels 

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

2910 and b.nslc_id[:3] == c.nslc_id[:3] 

2911 and abs(a.deltat-b.deltat) < a.deltat*0.001 

2912 and abs(b.deltat-c.deltat) < b.deltat*0.001): 

2913 

2914 continue 

2915 

2916 tmin = max(a.tmin, b.tmin, c.tmin) 

2917 tmax = min(a.tmax, b.tmax, c.tmax) 

2918 

2919 if tmin >= tmax: 

2920 continue 

2921 

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

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

2924 cc = c.chop(tmin, tmax, inplace=False, include_last=True) 

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

2926 or abs(bc.tmin - cc.tmin) > bc.deltat*0.01): 

2927 

2928 logger.warning( 

2929 'Cannot project traces with displaced sampling ' 

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

2931 continue 

2932 

2933 acydata = num.dot( 

2934 matrix[0], 

2935 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata())) 

2936 bcydata = num.dot( 

2937 matrix[1], 

2938 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata())) 

2939 ccydata = num.dot( 

2940 matrix[2], 

2941 (ac.get_ydata(), bc.get_ydata(), cc.get_ydata())) 

2942 

2943 ac.set_ydata(acydata) 

2944 bc.set_ydata(bcydata) 

2945 cc.set_ydata(ccydata) 

2946 

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

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

2949 cc.set_codes(channel=out_channels[2]) 

2950 

2951 projected.append(ac) 

2952 projected.append(bc) 

2953 projected.append(cc) 

2954 

2955 return projected 

2956 

2957 

2958def correlate(a, b, mode='valid', normalization=None, use_fft=False): 

2959 ''' 

2960 Cross correlation of two traces. 

2961 

2962 :param a,b: input traces 

2963 :param mode: ``'valid'``, ``'full'``, or ``'same'`` 

2964 :param normalization: ``'normal'``, ``'gliding'``, or ``None`` 

2965 :param use_fft: bool, whether to do cross correlation in spectral domain 

2966 

2967 :returns: trace containing cross correlation coefficients 

2968 

2969 This function computes the cross correlation between two traces. It 

2970 evaluates the discrete equivalent of 

2971 

2972 .. math:: 

2973 

2974 c(t) = \\int_{-\\infty}^{\\infty} a^{\\ast}(\\tau) b(t+\\tau) d\\tau 

2975 

2976 where the star denotes complex conjugate. Note, that the arguments here are 

2977 swapped when compared with the :py:func:`numpy.correlate` function, 

2978 which is internally called. This function should be safe even with older 

2979 versions of NumPy, where the correlate function has some problems. 

2980 

2981 A trace containing the cross correlation coefficients is returned. The time 

2982 information of the output trace is set so that the returned cross 

2983 correlation can be viewed directly as a function of time lag. 

2984 

2985 Example:: 

2986 

2987 # align two traces a and b containing a time shifted similar signal: 

2988 c = pyrocko.trace.correlate(a,b) 

2989 t, coef = c.max() # get time and value of maximum 

2990 b.shift(-t) # align b with a 

2991 

2992 ''' 

2993 

2994 assert_same_sampling_rate(a, b) 

2995 

2996 ya, yb = a.ydata, b.ydata 

2997 

2998 # need reversed order here: 

2999 yc = numpy_correlate_fixed(yb, ya, mode=mode, use_fft=use_fft) 

3000 kmin, kmax = numpy_correlate_lag_range(yb, ya, mode=mode, use_fft=use_fft) 

3001 

3002 if normalization == 'normal': 

3003 normfac = num.sqrt(num.sum(ya**2))*num.sqrt(num.sum(yb**2)) 

3004 yc = yc/normfac 

3005 

3006 elif normalization == 'gliding': 

3007 if mode != 'valid': 

3008 assert False, 'gliding normalization currently only available ' \ 

3009 'with "valid" mode.' 

3010 

3011 if ya.size < yb.size: 

3012 yshort, ylong = ya, yb 

3013 else: 

3014 yshort, ylong = yb, ya 

3015 

3016 epsilon = 0.00001 

3017 normfac_short = num.sqrt(num.sum(yshort**2)) 

3018 normfac = normfac_short * num.sqrt( 

3019 moving_sum(ylong**2, yshort.size, mode='valid')) \ 

3020 + normfac_short*epsilon 

3021 

3022 if yb.size <= ya.size: 

3023 normfac = normfac[::-1] 

3024 

3025 yc /= normfac 

3026 

3027 c = a.copy() 

3028 c.set_ydata(yc) 

3029 c.set_codes(*merge_codes(a, b, '~')) 

3030 c.shift(-c.tmin + b.tmin-a.tmin + kmin * c.deltat) 

3031 

3032 return c 

3033 

3034 

3035def deconvolve( 

3036 a, b, waterlevel, 

3037 tshift=0., 

3038 pad=0.5, 

3039 fd_taper=None, 

3040 pad_to_pow2=True): 

3041 

3042 same_sampling_rate(a, b) 

3043 assert abs(a.tmin - b.tmin) < a.deltat * 0.001 

3044 deltat = a.deltat 

3045 npad = int(round(a.data_len()*pad + tshift / deltat)) 

3046 

3047 ndata = max(a.data_len(), b.data_len()) 

3048 ndata_pad = ndata + npad 

3049 

3050 if pad_to_pow2: 

3051 ntrans = nextpow2(ndata_pad) 

3052 else: 

3053 ntrans = ndata 

3054 

3055 aspec = num.fft.rfft(a.ydata, ntrans) 

3056 bspec = num.fft.rfft(b.ydata, ntrans) 

3057 

3058 out = aspec * num.conj(bspec) 

3059 

3060 bautocorr = bspec*num.conj(bspec) 

3061 denom = num.maximum(bautocorr, waterlevel * bautocorr.max()) 

3062 

3063 out /= denom 

3064 df = 1/(ntrans*deltat) 

3065 

3066 if fd_taper is not None: 

3067 fd_taper(out, 0.0, df) 

3068 

3069 ydata = num.roll(num.fft.irfft(out), int(round(tshift/deltat))) 

3070 c = a.copy(data=False) 

3071 c.set_ydata(ydata[:ndata]) 

3072 c.set_codes(*merge_codes(a, b, '/')) 

3073 return c 

3074 

3075 

3076def assert_same_sampling_rate(a, b, eps=1.0e-6): 

3077 assert same_sampling_rate(a, b, eps), \ 

3078 'Sampling rates differ: %g != %g' % (a.deltat, b.deltat) 

3079 

3080 

3081def same_sampling_rate(a, b, eps=1.0e-6): 

3082 ''' 

3083 Check if two traces have the same sampling rate. 

3084 

3085 :param a,b: input traces 

3086 :param eps: relative tolerance 

3087 ''' 

3088 

3089 return abs(a.deltat - b.deltat) < (a.deltat + b.deltat)*eps 

3090 

3091 

3092def fix_deltat_rounding_errors(deltat): 

3093 ''' 

3094 Try to undo sampling rate rounding errors. 

3095 

3096 Fix rounding errors of sampling intervals when these are read from single 

3097 precision floating point values. 

3098 

3099 Assumes that the true sampling rate or sampling interval was an integer 

3100 value. No correction will be applied if this would change the sampling 

3101 rate by more than 0.001%. 

3102 ''' 

3103 

3104 if deltat <= 1.0: 

3105 deltat_new = 1.0 / round(1.0 / deltat) 

3106 else: 

3107 deltat_new = round(deltat) 

3108 

3109 if abs(deltat_new - deltat) / deltat > 1e-5: 

3110 deltat_new = deltat 

3111 

3112 return deltat_new 

3113 

3114 

3115def merge_codes(a, b, sep='-'): 

3116 ''' 

3117 Merge network-station-location-channel codes of a pair of traces. 

3118 ''' 

3119 

3120 o = [] 

3121 for xa, xb in zip(a.nslc_id, b.nslc_id): 

3122 if xa == xb: 

3123 o.append(xa) 

3124 else: 

3125 o.append(sep.join((xa, xb))) 

3126 return o 

3127 

3128 

3129class Taper(Object): 

3130 ''' 

3131 Base class for tapers. 

3132 

3133 Does nothing by default. 

3134 ''' 

3135 

3136 def __call__(self, y, x0, dx): 

3137 pass 

3138 

3139 

3140class CosTaper(Taper): 

3141 ''' 

3142 Cosine Taper. 

3143 

3144 :param a: start of fading in 

3145 :param b: end of fading in 

3146 :param c: start of fading out 

3147 :param d: end of fading out 

3148 ''' 

3149 

3150 a = Float.T() 

3151 b = Float.T() 

3152 c = Float.T() 

3153 d = Float.T() 

3154 

3155 def __init__(self, a, b, c, d): 

3156 Taper.__init__(self, a=a, b=b, c=c, d=d) 

3157 

3158 def __call__(self, y, x0, dx): 

3159 

3160 if y.dtype == num.dtype(float): 

3161 _apply_costaper = signal_ext.apply_costaper 

3162 else: 

3163 _apply_costaper = apply_costaper 

3164 

3165 _apply_costaper(self.a, self.b, self.c, self.d, y, x0, dx) 

3166 

3167 def span(self, y, x0, dx): 

3168 return span_costaper(self.a, self.b, self.c, self.d, y, x0, dx) 

3169 

3170 def time_span(self): 

3171 return self.a, self.d 

3172 

3173 

3174class CosFader(Taper): 

3175 ''' 

3176 Cosine Fader. 

3177 

3178 :param xfade: fade in and fade out time in seconds (optional) 

3179 :param xfrac: fade in and fade out as fraction between 0. and 1. (optional) 

3180 

3181 Only one argument can be set. The other should to be ``None``. 

3182 ''' 

3183 

3184 xfade = Float.T(optional=True) 

3185 xfrac = Float.T(optional=True) 

3186 

3187 def __init__(self, xfade=None, xfrac=None): 

3188 Taper.__init__(self, xfade=xfade, xfrac=xfrac) 

3189 assert (xfade is None) != (xfrac is None) 

3190 self._xfade = xfade 

3191 self._xfrac = xfrac 

3192 

3193 def __call__(self, y, x0, dx): 

3194 

3195 xfade = self._xfade 

3196 

3197 xlen = (y.size - 1)*dx 

3198 if xfade is None: 

3199 xfade = xlen * self._xfrac 

3200 

3201 a = x0 

3202 b = x0 + xfade 

3203 c = x0 + xlen - xfade 

3204 d = x0 + xlen 

3205 

3206 apply_costaper(a, b, c, d, y, x0, dx) 

3207 

3208 def span(self, y, x0, dx): 

3209 return 0, y.size 

3210 

3211 def time_span(self): 

3212 return None, None 

3213 

3214 

3215def none_min(li): 

3216 if None in li: 

3217 return None 

3218 else: 

3219 return min(x for x in li if x is not None) 

3220 

3221 

3222def none_max(li): 

3223 if None in li: 

3224 return None 

3225 else: 

3226 return max(x for x in li if x is not None) 

3227 

3228 

3229class MultiplyTaper(Taper): 

3230 ''' 

3231 Multiplication of several tapers. 

3232 ''' 

3233 

3234 tapers = List.T(Taper.T()) 

3235 

3236 def __init__(self, tapers=None): 

3237 if tapers is None: 

3238 tapers = [] 

3239 

3240 Taper.__init__(self, tapers=tapers) 

3241 

3242 def __call__(self, y, x0, dx): 

3243 for taper in self.tapers: 

3244 taper(y, x0, dx) 

3245 

3246 def span(self, y, x0, dx): 

3247 spans = [] 

3248 for taper in self.tapers: 

3249 spans.append(taper.span(y, x0, dx)) 

3250 

3251 mins, maxs = list(zip(*spans)) 

3252 return min(mins), max(maxs) 

3253 

3254 def time_span(self): 

3255 spans = [] 

3256 for taper in self.tapers: 

3257 spans.append(taper.time_span()) 

3258 

3259 mins, maxs = list(zip(*spans)) 

3260 return none_min(mins), none_max(maxs) 

3261 

3262 

3263class GaussTaper(Taper): 

3264 ''' 

3265 Frequency domain Gaussian filter. 

3266 ''' 

3267 

3268 alpha = Float.T() 

3269 

3270 def __init__(self, alpha): 

3271 Taper.__init__(self, alpha=alpha) 

3272 self._alpha = alpha 

3273 

3274 def __call__(self, y, x0, dx): 

3275 f = x0 + num.arange(y.size)*dx 

3276 y *= num.exp(-num.pi**2 / (self._alpha**2) * f**2) 

3277 

3278 

3279cached_coefficients = {} 

3280 

3281 

3282def _get_cached_filter_coeffs(order, corners, btype): 

3283 ck = (order, tuple(corners), btype) 

3284 if ck not in cached_coefficients: 

3285 if len(corners) == 1: 

3286 corners = corners[0] 

3287 

3288 cached_coefficients[ck] = signal.butter( 

3289 order, corners, btype=btype) 

3290 

3291 return cached_coefficients[ck] 

3292 

3293 

3294class _globals(object): 

3295 _numpy_has_correlate_flip_bug = None 

3296 

3297 

3298def _default_key(tr): 

3299 return (tr.network, tr.station, tr.location, tr.channel) 

3300 

3301 

3302def numpy_has_correlate_flip_bug(): 

3303 ''' 

3304 Check if NumPy's correlate function reveals old behaviour. 

3305 ''' 

3306 

3307 if _globals._numpy_has_correlate_flip_bug is None: 

3308 a = num.array([0, 0, 1, 0, 0, 0, 0]) 

3309 b = num.array([0, 0, 0, 0, 1, 0, 0, 0]) 

3310 ab = num.correlate(a, b, mode='same') 

3311 ba = num.correlate(b, a, mode='same') 

3312 _globals._numpy_has_correlate_flip_bug = num.all(ab == ba) 

3313 

3314 return _globals._numpy_has_correlate_flip_bug 

3315 

3316 

3317def numpy_correlate_fixed(a, b, mode='valid', use_fft=False): 

3318 ''' 

3319 Call :py:func:`numpy.correlate` with fixes. 

3320 

3321 c[k] = sum_i a[i+k] * conj(b[i]) 

3322 

3323 Note that the result produced by newer numpy.correlate is always flipped 

3324 with respect to the formula given in its documentation (if ascending k 

3325 assumed for the output). 

3326 ''' 

3327 

3328 if use_fft: 

3329 if a.size < b.size: 

3330 c = signal.fftconvolve(b[::-1], a, mode=mode) 

3331 else: 

3332 c = signal.fftconvolve(a, b[::-1], mode=mode) 

3333 return c 

3334 

3335 else: 

3336 buggy = numpy_has_correlate_flip_bug() 

3337 

3338 a = num.asarray(a) 

3339 b = num.asarray(b) 

3340 

3341 if buggy: 

3342 b = num.conj(b) 

3343 

3344 c = num.correlate(a, b, mode=mode) 

3345 

3346 if buggy and a.size < b.size: 

3347 return c[::-1] 

3348 else: 

3349 return c 

3350 

3351 

3352def numpy_correlate_emulate(a, b, mode='valid'): 

3353 ''' 

3354 Slow version of :py:func:`numpy.correlate` for comparison. 

3355 ''' 

3356 

3357 a = num.asarray(a) 

3358 b = num.asarray(b) 

3359 kmin = -(b.size-1) 

3360 klen = a.size-kmin 

3361 kmin, kmax = numpy_correlate_lag_range(a, b, mode=mode) 

3362 kmin = int(kmin) 

3363 kmax = int(kmax) 

3364 klen = kmax - kmin + 1 

3365 c = num.zeros(klen, dtype=num.promote_types(b.dtype, a.dtype)) 

3366 for k in range(kmin, kmin+klen): 

3367 imin = max(0, -k) 

3368 ilen = min(b.size, a.size-k) - imin 

3369 c[k-kmin] = num.sum( 

3370 a[imin+k:imin+ilen+k] * num.conj(b[imin:imin+ilen])) 

3371 

3372 return c 

3373 

3374 

3375def numpy_correlate_lag_range(a, b, mode='valid', use_fft=False): 

3376 ''' 

3377 Get range of lags for which :py:func:`numpy.correlate` produces values. 

3378 ''' 

3379 

3380 a = num.asarray(a) 

3381 b = num.asarray(b) 

3382 

3383 kmin = -(b.size-1) 

3384 if mode == 'full': 

3385 klen = a.size-kmin 

3386 elif mode == 'same': 

3387 klen = max(a.size, b.size) 

3388 kmin += (a.size+b.size-1 - max(a.size, b.size)) // 2 + \ 

3389 int(not use_fft and a.size % 2 == 0 and b.size > a.size) 

3390 elif mode == 'valid': 

3391 klen = abs(a.size - b.size) + 1 

3392 kmin += min(a.size, b.size) - 1 

3393 

3394 return kmin, kmin + klen - 1 

3395 

3396 

3397def autocorr(x, nshifts): 

3398 ''' 

3399 Compute biased estimate of the first autocorrelation coefficients. 

3400 

3401 :param x: input array 

3402 :param nshifts: number of coefficients to calculate 

3403 ''' 

3404 

3405 mean = num.mean(x) 

3406 std = num.std(x) 

3407 n = x.size 

3408 xdm = x - mean 

3409 r = num.zeros(nshifts) 

3410 for k in range(nshifts): 

3411 r[k] = 1./((n-num.abs(k))*std) * num.sum(xdm[:n-k] * xdm[k:]) 

3412 

3413 return r 

3414 

3415 

3416def yulewalker(x, order): 

3417 ''' 

3418 Compute autoregression coefficients using Yule-Walker method. 

3419 

3420 :param x: input array 

3421 :param order: number of coefficients to produce 

3422 

3423 A biased estimate of the autocorrelation is used. The Yule-Walker equations 

3424 are solved by :py:func:`numpy.linalg.inv` instead of Levinson-Durbin 

3425 recursion which is normally used. 

3426 ''' 

3427 

3428 gamma = autocorr(x, order+1) 

3429 d = gamma[1:1+order] 

3430 a = num.zeros((order, order)) 

3431 gamma2 = num.concatenate((gamma[::-1], gamma[1:order])) 

3432 for i in range(order): 

3433 ioff = order-i 

3434 a[i, :] = gamma2[ioff:ioff+order] 

3435 

3436 return num.dot(num.linalg.inv(a), -d) 

3437 

3438 

3439def moving_avg(x, n): 

3440 n = int(n) 

3441 cx = x.cumsum() 

3442 nn = len(x) 

3443 y = num.zeros(nn, dtype=cx.dtype) 

3444 y[n//2:n//2+(nn-n)] = (cx[n:]-cx[:-n])/n 

3445 y[:n//2] = y[n//2] 

3446 y[n//2+(nn-n):] = y[n//2+(nn-n)-1] 

3447 return y 

3448 

3449 

3450def moving_sum(x, n, mode='valid'): 

3451 n = int(n) 

3452 cx = x.cumsum() 

3453 nn = len(x) 

3454 

3455 if mode == 'valid': 

3456 if nn-n+1 <= 0: 

3457 return num.zeros(0, dtype=cx.dtype) 

3458 y = num.zeros(nn-n+1, dtype=cx.dtype) 

3459 y[0] = cx[n-1] 

3460 y[1:nn-n+1] = cx[n:nn]-cx[0:nn-n] 

3461 

3462 if mode == 'full': 

3463 y = num.zeros(nn+n-1, dtype=cx.dtype) 

3464 if n <= nn: 

3465 y[0:n] = cx[0:n] 

3466 y[n:nn] = cx[n:nn]-cx[0:nn-n] 

3467 y[nn:nn+n-1] = cx[-1]-cx[nn-n:nn-1] 

3468 else: 

3469 y[0:nn] = cx[0:nn] 

3470 y[nn:n] = cx[nn-1] 

3471 y[n:nn+n-1] = cx[nn-1] - cx[0:nn-1] 

3472 

3473 if mode == 'same': 

3474 n1 = (n-1)//2 

3475 y = num.zeros(nn, dtype=cx.dtype) 

3476 if n <= nn: 

3477 y[0:n-n1] = cx[n1:n] 

3478 y[n-n1:nn-n1] = cx[n:nn]-cx[0:nn-n] 

3479 y[nn-n1:nn] = cx[nn-1] - cx[nn-n:nn-n+n1] 

3480 else: 

3481 y[0:max(0, nn-n1)] = cx[min(n1, nn):nn] 

3482 y[max(nn-n1, 0):min(n-n1, nn)] = cx[nn-1] 

3483 y[min(n-n1, nn):nn] = cx[nn-1] - cx[0:max(0, nn-(n-n1))] 

3484 

3485 return y 

3486 

3487 

3488def nextpow2(i): 

3489 return 2**int(math.ceil(math.log(i)/math.log(2.))) 

3490 

3491 

3492def snapper_w_offset(nmax, offset, delta, snapfun=math.ceil): 

3493 def snap(x): 

3494 return max(0, min(int(snapfun((x-offset)/delta)), nmax)) 

3495 return snap 

3496 

3497 

3498def snapper(nmax, delta, snapfun=math.ceil): 

3499 def snap(x): 

3500 return max(0, min(int(snapfun(x/delta)), nmax)) 

3501 return snap 

3502 

3503 

3504def apply_costaper(a, b, c, d, y, x0, dx): 

3505 abcd = num.array((a, b, c, d), dtype=float) 

3506 ja, jb, jc, jd = num.clip(num.ceil((abcd-x0)/dx).astype(int), 0, y.size) 

3507 y[:ja] = 0. 

3508 y[ja:jb] *= 0.5 \ 

3509 - 0.5*num.cos((dx*num.arange(ja, jb)-(a-x0))/(b-a)*num.pi) 

3510 y[jc:jd] *= 0.5 \ 

3511 + 0.5*num.cos((dx*num.arange(jc, jd)-(c-x0))/(d-c)*num.pi) 

3512 y[jd:] = 0. 

3513 

3514 

3515def span_costaper(a, b, c, d, y, x0, dx): 

3516 hi = snapper_w_offset(y.size, x0, dx) 

3517 return hi(a), hi(d) - hi(a) 

3518 

3519 

3520def costaper(a, b, c, d, nfreqs, deltaf): 

3521 hi = snapper(nfreqs, deltaf) 

3522 tap = num.zeros(nfreqs) 

3523 tap[hi(a):hi(b)] = 0.5 \ 

3524 - 0.5*num.cos((deltaf*num.arange(hi(a), hi(b))-a)/(b-a)*num.pi) 

3525 tap[hi(b):hi(c)] = 1. 

3526 tap[hi(c):hi(d)] = 0.5 \ 

3527 + 0.5*num.cos((deltaf*num.arange(hi(c), hi(d))-c)/(d-c)*num.pi) 

3528 

3529 return tap 

3530 

3531 

3532def t2ind(t, tdelta, snap=round): 

3533 return int(snap(t/tdelta)) 

3534 

3535 

3536def hilbert(x, N=None): 

3537 ''' 

3538 Return the hilbert transform of x of length N. 

3539 

3540 (from scipy.signal, but changed to use fft and ifft from numpy.fft) 

3541 ''' 

3542 

3543 x = num.asarray(x) 

3544 if N is None: 

3545 N = len(x) 

3546 if N <= 0: 

3547 raise ValueError('N must be positive.') 

3548 if num.iscomplexobj(x): 

3549 logger.warning('imaginary part of x ignored.') 

3550 x = num.real(x) 

3551 

3552 Xf = num.fft.fft(x, N, axis=0) 

3553 h = num.zeros(N) 

3554 if N % 2 == 0: 

3555 h[0] = h[N//2] = 1 

3556 h[1:N//2] = 2 

3557 else: 

3558 h[0] = 1 

3559 h[1:(N+1)//2] = 2 

3560 

3561 if len(x.shape) > 1: 

3562 h = h[:, num.newaxis] 

3563 x = num.fft.ifft(Xf*h) 

3564 return x 

3565 

3566 

3567def near(a, b, eps): 

3568 return abs(a-b) < eps 

3569 

3570 

3571def coroutine(func): 

3572 def wrapper(*args, **kwargs): 

3573 gen = func(*args, **kwargs) 

3574 next(gen) 

3575 return gen 

3576 

3577 wrapper.__name__ = func.__name__ 

3578 wrapper.__dict__ = func.__dict__ 

3579 wrapper.__doc__ = func.__doc__ 

3580 return wrapper 

3581 

3582 

3583class States(object): 

3584 ''' 

3585 Utility to store channel-specific state in coroutines. 

3586 ''' 

3587 

3588 def __init__(self): 

3589 self._states = {} 

3590 

3591 def get(self, tr): 

3592 k = tr.nslc_id 

3593 if k in self._states: 

3594 tmin, deltat, dtype, value = self._states[k] 

3595 if (near(tmin, tr.tmin, deltat/100.) 

3596 and near(deltat, tr.deltat, deltat/10000.) 

3597 and dtype == tr.ydata.dtype): 

3598 

3599 return value 

3600 

3601 return None 

3602 

3603 def set(self, tr, value): 

3604 k = tr.nslc_id 

3605 if k in self._states and self._states[k][-1] is not value: 

3606 self.free(self._states[k][-1]) 

3607 

3608 self._states[k] = (tr.tmax+tr.deltat, tr.deltat, tr.ydata.dtype, value) 

3609 

3610 def free(self, value): 

3611 pass 

3612 

3613 

3614@coroutine 

3615def co_list_append(list): 

3616 while True: 

3617 list.append((yield)) 

3618 

3619 

3620class ScipyBug(Exception): 

3621 pass 

3622 

3623 

3624@coroutine 

3625def co_lfilter(target, b, a): 

3626 ''' 

3627 Successively filter broken continuous trace data (coroutine). 

3628 

3629 Create coroutine which takes :py:class:`Trace` objects, filters their data 

3630 through :py:func:`scipy.signal.lfilter` and sends new :py:class:`Trace` 

3631 objects containing the filtered data to target. This is useful, if one 

3632 wants to filter a long continuous time series, which is split into many 

3633 successive traces without producing filter artifacts at trace boundaries. 

3634 

3635 Filter states are kept *per channel*, specifically, for each (network, 

3636 station, location, channel) combination occuring in the input traces, a 

3637 separate state is created and maintained. This makes it possible to filter 

3638 multichannel or multistation data with only one :py:func:`co_lfilter` 

3639 instance. 

3640 

3641 Filter state is reset, when gaps occur. 

3642 

3643 Use it like this:: 

3644 

3645 from pyrocko.trace import co_lfilter, co_list_append 

3646 

3647 filtered_traces = [] 

3648 pipe = co_lfilter(co_list_append(filtered_traces), a, b) 

3649 for trace in traces: 

3650 pipe.send(trace) 

3651 

3652 pipe.close() 

3653 

3654 ''' 

3655 

3656 try: 

3657 states = States() 

3658 output = None 

3659 while True: 

3660 input = (yield) 

3661 

3662 zi = states.get(input) 

3663 if zi is None: 

3664 zi = num.zeros(max(len(a), len(b))-1, dtype=float) 

3665 

3666 output = input.copy(data=False) 

3667 try: 

3668 ydata, zf = signal.lfilter(b, a, input.get_ydata(), zi=zi) 

3669 except ValueError: 

3670 raise ScipyBug( 

3671 'signal.lfilter failed: could be related to a bug ' 

3672 'in some older scipy versions, e.g. on opensuse42.1') 

3673 

3674 output.set_ydata(ydata) 

3675 states.set(input, zf) 

3676 target.send(output) 

3677 

3678 except GeneratorExit: 

3679 target.close() 

3680 

3681 

3682def co_antialias(target, q, n=None, ftype='fir'): 

3683 b, a, n = util.decimate_coeffs(q, n, ftype) 

3684 anti = co_lfilter(target, b, a) 

3685 return anti 

3686 

3687 

3688@coroutine 

3689def co_dropsamples(target, q, nfir): 

3690 try: 

3691 states = States() 

3692 while True: 

3693 tr = (yield) 

3694 newdeltat = q * tr.deltat 

3695 ioffset = states.get(tr) 

3696 if ioffset is None: 

3697 # for fir filter, the first nfir samples are pulluted by 

3698 # boundary effects; cut it off. 

3699 # for iir this may be (much) more, we do not correct for that. 

3700 # put sample instances to a time which is a multiple of the 

3701 # new sampling interval. 

3702 newtmin_want = math.ceil( 

3703 (tr.tmin+(nfir+1)*tr.deltat) / newdeltat) * newdeltat \ 

3704 - (nfir/2*tr.deltat) 

3705 ioffset = int(round((newtmin_want - tr.tmin)/tr.deltat)) 

3706 if ioffset < 0: 

3707 ioffset = ioffset % q 

3708 

3709 newtmin_have = tr.tmin + ioffset * tr.deltat 

3710 newtr = tr.copy(data=False) 

3711 newtr.deltat = newdeltat 

3712 # because the fir kernel shifts data by nfir/2 samples: 

3713 newtr.tmin = newtmin_have - (nfir/2*tr.deltat) 

3714 newtr.set_ydata(tr.get_ydata()[ioffset::q].copy()) 

3715 states.set(tr, (ioffset % q - tr.data_len() % q) % q) 

3716 target.send(newtr) 

3717 

3718 except GeneratorExit: 

3719 target.close() 

3720 

3721 

3722def co_downsample(target, q, n=None, ftype='fir'): 

3723 ''' 

3724 Successively downsample broken continuous trace data (coroutine). 

3725 

3726 Create coroutine which takes :py:class:`Trace` objects, downsamples their 

3727 data and sends new :py:class:`Trace` objects containing the downsampled 

3728 data to target. This is useful, if one wants to downsample a long 

3729 continuous time series, which is split into many successive traces without 

3730 producing filter artifacts and gaps at trace boundaries. 

3731 

3732 Filter states are kept *per channel*, specifically, for each (network, 

3733 station, location, channel) combination occuring in the input traces, a 

3734 separate state is created and maintained. This makes it possible to filter 

3735 multichannel or multistation data with only one :py:func:`co_lfilter` 

3736 instance. 

3737 

3738 Filter state is reset, when gaps occur. The sampling instances are choosen 

3739 so that they occur at (or as close as possible) to even multiples of the 

3740 sampling interval of the downsampled trace (based on system time). 

3741 ''' 

3742 

3743 b, a, n = util.decimate_coeffs(q, n, ftype) 

3744 return co_antialias(co_dropsamples(target, q, n), q, n, ftype) 

3745 

3746 

3747@coroutine 

3748def co_downsample_to(target, deltat): 

3749 

3750 decimators = {} 

3751 try: 

3752 while True: 

3753 tr = (yield) 

3754 ratio = deltat / tr.deltat 

3755 rratio = round(ratio) 

3756 if abs(rratio - ratio)/ratio > 0.0001: 

3757 raise util.UnavailableDecimation('ratio = %g' % ratio) 

3758 

3759 deci_seq = tuple(x for x in util.decitab(int(rratio)) if x != 1) 

3760 if deci_seq not in decimators: 

3761 pipe = target 

3762 for q in deci_seq[::-1]: 

3763 pipe = co_downsample(pipe, q) 

3764 

3765 decimators[deci_seq] = pipe 

3766 

3767 decimators[deci_seq].send(tr) 

3768 

3769 except GeneratorExit: 

3770 for g in decimators.values(): 

3771 g.close() 

3772 

3773 

3774class DomainChoice(StringChoice): 

3775 choices = [ 

3776 'time_domain', 

3777 'frequency_domain', 

3778 'envelope', 

3779 'absolute', 

3780 'cc_max_norm'] 

3781 

3782 

3783class MisfitSetup(Object): 

3784 ''' 

3785 Contains misfit setup to be used in :py:meth:`Trace.misfit` 

3786 

3787 :param description: Description of the setup 

3788 :param norm: L-norm classifier 

3789 :param taper: Object of :py:class:`Taper` 

3790 :param filter: Object of :py:class:`~pyrocko.response.FrequencyResponse` 

3791 :param domain: ['time_domain', 'frequency_domain', 'envelope', 'absolute', 

3792 'cc_max_norm'] 

3793 

3794 Can be dumped to a yaml file. 

3795 ''' 

3796 

3797 xmltagname = 'misfitsetup' 

3798 description = String.T(optional=True) 

3799 norm = Int.T(optional=False) 

3800 taper = Taper.T(optional=False) 

3801 filter = FrequencyResponse.T(optional=True) 

3802 domain = DomainChoice.T(default='time_domain') 

3803 

3804 

3805def equalize_sampling_rates(trace_1, trace_2): 

3806 ''' 

3807 Equalize sampling rates of two traces (reduce higher sampling rate to 

3808 lower). 

3809 

3810 :param trace_1: :py:class:`Trace` object 

3811 :param trace_2: :py:class:`Trace` object 

3812 

3813 Returns a copy of the resampled trace if resampling is needed. 

3814 ''' 

3815 

3816 if same_sampling_rate(trace_1, trace_2): 

3817 return trace_1, trace_2 

3818 

3819 if trace_1.deltat < trace_2.deltat: 

3820 t1_out = trace_1.copy() 

3821 t1_out.downsample_to(deltat=trace_2.deltat, snap=True) 

3822 logger.debug('Trace downsampled (return copy of trace): %s' 

3823 % '.'.join(t1_out.nslc_id)) 

3824 return t1_out, trace_2 

3825 

3826 elif trace_1.deltat > trace_2.deltat: 

3827 t2_out = trace_2.copy() 

3828 t2_out.downsample_to(deltat=trace_1.deltat, snap=True) 

3829 logger.debug('Trace downsampled (return copy of trace): %s' 

3830 % '.'.join(t2_out.nslc_id)) 

3831 return trace_1, t2_out 

3832 

3833 

3834def Lx_norm(u, v, norm=2): 

3835 ''' 

3836 Calculate the misfit denominator *m* and the normalization divisor *n* 

3837 according to norm. 

3838 

3839 The normalization divisor *n* is calculated from ``v``. 

3840 

3841 :param u: :py:class:`numpy.ndarray` 

3842 :param v: :py:class:`numpy.ndarray` 

3843 :param norm: (default = 2) 

3844 

3845 ``u`` and ``v`` must be of same size. 

3846 ''' 

3847 

3848 if norm == 1: 

3849 return ( 

3850 num.sum(num.abs(v-u)), 

3851 num.sum(num.abs(v))) 

3852 

3853 elif norm == 2: 

3854 return ( 

3855 num.sqrt(num.sum((v-u)**2)), 

3856 num.sqrt(num.sum(v**2))) 

3857 

3858 else: 

3859 return ( 

3860 num.power(num.sum(num.abs(num.power(v - u, norm))), 1./norm), 

3861 num.power(num.sum(num.abs(num.power(v, norm))), 1./norm)) 

3862 

3863 

3864def do_downsample(tr, deltat): 

3865 if abs(tr.deltat - deltat) / tr.deltat > 1e-6: 

3866 tr = tr.copy() 

3867 tr.downsample_to(deltat, snap=True, demean=False) 

3868 else: 

3869 if tr.tmin/tr.deltat > 1e-6 or tr.tmax/tr.deltat > 1e-6: 

3870 tr = tr.copy() 

3871 tr.snap() 

3872 return tr 

3873 

3874 

3875def do_extend(tr, tmin, tmax): 

3876 if tmin < tr.tmin or tmax > tr.tmax: 

3877 tr = tr.copy() 

3878 tr.extend(tmin=tmin, tmax=tmax, fillmethod='repeat') 

3879 

3880 return tr 

3881 

3882 

3883def do_pre_taper(tr, taper): 

3884 return tr.taper(taper, inplace=False, chop=True) 

3885 

3886 

3887def do_fft(tr, filter): 

3888 if filter is None: 

3889 return tr 

3890 else: 

3891 ndata = tr.ydata.size 

3892 nfft = nextpow2(ndata) 

3893 padded = num.zeros(nfft, dtype=float) 

3894 padded[:ndata] = tr.ydata 

3895 spectrum = num.fft.rfft(padded) 

3896 df = 1.0 / (tr.deltat * nfft) 

3897 frequencies = num.arange(spectrum.size)*df 

3898 return [tr, frequencies, spectrum] 

3899 

3900 

3901def do_filter(inp, filter): 

3902 if filter is None: 

3903 return inp 

3904 else: 

3905 tr, frequencies, spectrum = inp 

3906 spectrum *= filter.evaluate(frequencies) 

3907 return [tr, frequencies, spectrum] 

3908 

3909 

3910def do_ifft(inp): 

3911 if isinstance(inp, Trace): 

3912 return inp 

3913 else: 

3914 tr, _, spectrum = inp 

3915 ndata = tr.ydata.size 

3916 tr = tr.copy(data=False) 

3917 tr.set_ydata(num.fft.irfft(spectrum)[:ndata]) 

3918 return tr 

3919 

3920 

3921def check_alignment(t1, t2): 

3922 if abs(t1.tmin-t2.tmin) > t1.deltat * 1e-4 or \ 

3923 abs(t1.tmax - t2.tmax) > t1.deltat * 1e-4 or \ 

3924 t1.ydata.shape != t2.ydata.shape: 

3925 raise MisalignedTraces( 

3926 'Cannot calculate misfit of %s and %s due to misaligned ' 

3927 'traces.' % ('.'.join(t1.nslc_id), '.'.join(t2.nslc_id)))