1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Data model and content types handled by the Squirrel framework. 

8 

9Squirrel uses flat content types to represent waveform, station, channel, 

10response, event, and a few other objects. A common subset of the information in 

11these objects is indexed in the database, currently: kind, codes, time interval 

12and sampling rate. The :py:class:`Nut` objects encapsulate this information 

13together with information about where and how to get the associated bulk data. 

14 

15Further content types are defined here to handle waveform orders, waveform 

16promises, data coverage and sensor information. 

17''' 

18 

19from __future__ import absolute_import, print_function 

20 

21import re 

22import fnmatch 

23import hashlib 

24import numpy as num 

25from os import urandom 

26from base64 import urlsafe_b64encode 

27from collections import defaultdict, namedtuple 

28 

29from pyrocko import util 

30from pyrocko.guts import Object, SObject, String, Timestamp, Float, Int, \ 

31 Unicode, Tuple, List, StringChoice, Any, Dict 

32from pyrocko.model import squirrel_content, Location 

33from pyrocko.response import FrequencyResponse, MultiplyResponse, \ 

34 IntegrationResponse, DifferentiationResponse, simplify_responses, \ 

35 FrequencyResponseCheckpoint 

36 

37from .error import ConversionError, SquirrelError 

38 

39 

40guts_prefix = 'squirrel' 

41 

42 

43g_codes_pool = {} 

44 

45 

46class CodesError(SquirrelError): 

47 pass 

48 

49 

50class Codes(SObject): 

51 pass 

52 

53 

54def normalize_nslce(*args, **kwargs): 

55 if args and kwargs: 

56 raise ValueError('Either *args or **kwargs accepted, not both.') 

57 

58 if len(args) == 1: 

59 if isinstance(args[0], str): 

60 args = tuple(args[0].split('.')) 

61 elif isinstance(args[0], tuple): 

62 args = args[0] 

63 else: 

64 raise ValueError('Invalid argument type: %s' % type(args[0])) 

65 

66 nargs = len(args) 

67 if nargs == 5: 

68 t = args 

69 

70 elif nargs == 4: 

71 t = args + ('',) 

72 

73 elif nargs == 0: 

74 d = dict( 

75 network='', 

76 station='', 

77 location='', 

78 channel='', 

79 extra='') 

80 

81 d.update(kwargs) 

82 t = tuple(kwargs.get(k, '') for k in ( 

83 'network', 'station', 'location', 'channel', 'extra')) 

84 

85 else: 

86 raise CodesError( 

87 'Does not match NSLC or NSLCE codes pattern: %s' % '.'.join(args)) 

88 

89 if '.'.join(t).count('.') != 4: 

90 raise CodesError( 

91 'Codes may not contain a ".": "%s", "%s", "%s", "%s", "%s"' % t) 

92 

93 return t 

94 

95 

96CodesNSLCEBase = namedtuple( 

97 'CodesNSLCEBase', [ 

98 'network', 'station', 'location', 'channel', 'extra']) 

99 

100 

101class CodesNSLCE(CodesNSLCEBase, Codes): 

102 ''' 

103 Codes denominating a seismic channel (NSLC or NSLCE). 

104 

105 FDSN/SEED style NET.STA.LOC.CHA is accepted or NET.STA.LOC.CHA.EXTRA, where 

106 the EXTRA part in the latter form can be used to identify a custom 

107 processing applied to a channel. 

108 ''' 

109 

110 __slots__ = () 

111 __hash__ = CodesNSLCEBase.__hash__ 

112 

113 as_dict = CodesNSLCEBase._asdict 

114 

115 def __new__(cls, *args, safe_str=None, **kwargs): 

116 nargs = len(args) 

117 if nargs == 1 and isinstance(args[0], CodesNSLCE): 

118 return args[0] 

119 elif nargs == 1 and isinstance(args[0], CodesNSL): 

120 t = (args[0].nsl) + ('*', '*') 

121 elif nargs == 1 and isinstance(args[0], CodesX): 

122 t = ('*', '*', '*', '*', '*') 

123 elif safe_str is not None: 

124 t = safe_str.split('.') 

125 else: 

126 t = normalize_nslce(*args, **kwargs) 

127 

128 x = CodesNSLCEBase.__new__(cls, *t) 

129 return g_codes_pool.setdefault(x, x) 

130 

131 def __init__(self, *args, **kwargs): 

132 Codes.__init__(self) 

133 

134 def __str__(self): 

135 if self.extra == '': 

136 return '.'.join(self[:-1]) 

137 else: 

138 return '.'.join(self) 

139 

140 def __eq__(self, other): 

141 if not isinstance(other, CodesNSLCE): 

142 other = CodesNSLCE(other) 

143 

144 return CodesNSLCEBase.__eq__(self, other) 

145 

146 @property 

147 def safe_str(self): 

148 return '.'.join(self) 

149 

150 @property 

151 def nslce(self): 

152 return self[:4] 

153 

154 @property 

155 def nslc(self): 

156 return self[:4] 

157 

158 @property 

159 def nsl(self): 

160 return self[:3] 

161 

162 @property 

163 def ns(self): 

164 return self[:2] 

165 

166 def as_tuple(self): 

167 return tuple(self) 

168 

169 def replace(self, **kwargs): 

170 x = CodesNSLCEBase._replace(self, **kwargs) 

171 return g_codes_pool.setdefault(x, x) 

172 

173 

174def normalize_nsl(*args, **kwargs): 

175 if args and kwargs: 

176 raise ValueError('Either *args or **kwargs accepted, not both.') 

177 

178 if len(args) == 1: 

179 if isinstance(args[0], str): 

180 args = tuple(args[0].split('.')) 

181 elif isinstance(args[0], tuple): 

182 args = args[0] 

183 else: 

184 raise ValueError('Invalid argument type: %s' % type(args[0])) 

185 

186 nargs = len(args) 

187 if nargs == 3: 

188 t = args 

189 

190 elif nargs == 0: 

191 d = dict( 

192 network='', 

193 station='', 

194 location='') 

195 

196 d.update(kwargs) 

197 t = tuple(kwargs.get(k, '') for k in ( 

198 'network', 'station', 'location')) 

199 

200 else: 

201 raise CodesError( 

202 'Does not match NSL codes pattern: %s' % '.'.join(args)) 

203 

204 if '.'.join(t).count('.') != 2: 

205 raise CodesError( 

206 'Codes may not contain a ".": "%s", "%s", "%s"' % t) 

207 

208 return t 

209 

210 

211CodesNSLBase = namedtuple( 

212 'CodesNSLBase', [ 

213 'network', 'station', 'location']) 

214 

215 

216class CodesNSL(CodesNSLBase, Codes): 

217 ''' 

218 Codes denominating a seismic station (NSL). 

219 

220 NET.STA.LOC is accepted, slightly different from SEED/StationXML, where 

221 LOC is part of the channel. By setting location='*' is possible to get 

222 compatible behaviour in most cases. 

223 ''' 

224 

225 __slots__ = () 

226 __hash__ = CodesNSLBase.__hash__ 

227 

228 as_dict = CodesNSLBase._asdict 

229 

230 def __new__(cls, *args, safe_str=None, **kwargs): 

231 nargs = len(args) 

232 if nargs == 1 and isinstance(args[0], CodesNSL): 

233 return args[0] 

234 elif nargs == 1 and isinstance(args[0], CodesNSLCE): 

235 t = args[0].nsl 

236 elif nargs == 1 and isinstance(args[0], CodesX): 

237 t = ('*', '*', '*') 

238 elif safe_str is not None: 

239 t = safe_str.split('.') 

240 else: 

241 t = normalize_nsl(*args, **kwargs) 

242 

243 x = CodesNSLBase.__new__(cls, *t) 

244 return g_codes_pool.setdefault(x, x) 

245 

246 def __init__(self, *args, **kwargs): 

247 Codes.__init__(self) 

248 

249 def __str__(self): 

250 return '.'.join(self) 

251 

252 def __eq__(self, other): 

253 if not isinstance(other, CodesNSL): 

254 other = CodesNSL(other) 

255 

256 return CodesNSLBase.__eq__(self, other) 

257 

258 @property 

259 def safe_str(self): 

260 return '.'.join(self) 

261 

262 @property 

263 def ns(self): 

264 return self[:2] 

265 

266 @property 

267 def nsl(self): 

268 return self[:3] 

269 

270 def as_tuple(self): 

271 return tuple(self) 

272 

273 def replace(self, **kwargs): 

274 x = CodesNSLBase._replace(self, **kwargs) 

275 return g_codes_pool.setdefault(x, x) 

276 

277 

278CodesXBase = namedtuple( 

279 'CodesXBase', [ 

280 'name']) 

281 

282 

283class CodesX(CodesXBase, Codes): 

284 ''' 

285 General purpose codes for anything other than channels or stations. 

286 ''' 

287 

288 __slots__ = () 

289 __hash__ = CodesXBase.__hash__ 

290 __eq__ = CodesXBase.__eq__ 

291 

292 as_dict = CodesXBase._asdict 

293 

294 def __new__(cls, name='', safe_str=None): 

295 if isinstance(name, CodesX): 

296 return name 

297 elif isinstance(name, (CodesNSLCE, CodesNSL)): 

298 name = '*' 

299 elif safe_str is not None: 

300 name = safe_str 

301 else: 

302 if '.' in name: 

303 raise CodesError('Code may not contain a ".": %s' % name) 

304 

305 x = CodesXBase.__new__(cls, name) 

306 return g_codes_pool.setdefault(x, x) 

307 

308 def __init__(self, *args, **kwargs): 

309 Codes.__init__(self) 

310 

311 def __str__(self): 

312 return '.'.join(self) 

313 

314 @property 

315 def safe_str(self): 

316 return '.'.join(self) 

317 

318 @property 

319 def ns(self): 

320 return self[:2] 

321 

322 def as_tuple(self): 

323 return tuple(self) 

324 

325 def replace(self, **kwargs): 

326 x = CodesXBase._replace(self, **kwargs) 

327 return g_codes_pool.setdefault(x, x) 

328 

329 

330g_codes_patterns = {} 

331 

332 

333def match_codes(pattern, codes): 

334 spattern = pattern.safe_str 

335 scodes = codes.safe_str 

336 if spattern not in g_codes_patterns: 

337 rpattern = re.compile(fnmatch.translate(spattern), re.I) 

338 g_codes_patterns[spattern] = rpattern 

339 

340 rpattern = g_codes_patterns[spattern] 

341 return bool(rpattern.match(scodes)) 

342 

343 

344g_content_kinds = [ 

345 'undefined', 

346 'waveform', 

347 'station', 

348 'channel', 

349 'response', 

350 'event', 

351 'waveform_promise'] 

352 

353 

354g_codes_classes = [ 

355 CodesX, 

356 CodesNSLCE, 

357 CodesNSL, 

358 CodesNSLCE, 

359 CodesNSLCE, 

360 CodesX, 

361 CodesNSLCE] 

362 

363g_codes_classes_ndot = { 

364 0: CodesX, 

365 2: CodesNSL, 

366 3: CodesNSLCE, 

367 4: CodesNSLCE} 

368 

369 

370def to_codes_simple(kind_id, codes_safe_str): 

371 return g_codes_classes[kind_id](safe_str=codes_safe_str) 

372 

373 

374def to_codes(kind_id, obj): 

375 return g_codes_classes[kind_id](obj) 

376 

377 

378def to_codes_guess(s): 

379 try: 

380 return g_codes_classes_ndot[s.count('.')](s) 

381 except KeyError: 

382 raise CodesError('Cannot guess codes type: %s' % s) 

383 

384 

385g_content_kind_ids = ( 

386 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT, 

387 WAVEFORM_PROMISE) = range(len(g_content_kinds)) 

388 

389 

390g_tmin, g_tmax = util.get_working_system_time_range()[:2] 

391 

392 

393try: 

394 g_tmin_queries = max(g_tmin, util.str_to_time_fillup('1900-01-01')) 

395except Exception: 

396 g_tmin_queries = g_tmin 

397 

398 

399def to_kind(kind_id): 

400 return g_content_kinds[kind_id] 

401 

402 

403def to_kinds(kind_ids): 

404 return [g_content_kinds[kind_id] for kind_id in kind_ids] 

405 

406 

407def to_kind_id(kind): 

408 return g_content_kinds.index(kind) 

409 

410 

411def to_kind_ids(kinds): 

412 return [g_content_kinds.index(kind) for kind in kinds] 

413 

414 

415g_kind_mask_all = 0xff 

416 

417 

418def to_kind_mask(kinds): 

419 if kinds: 

420 return sum(1 << kind_id for kind_id in to_kind_ids(kinds)) 

421 else: 

422 return g_kind_mask_all 

423 

424 

425def str_or_none(x): 

426 if x is None: 

427 return None 

428 else: 

429 return str(x) 

430 

431 

432def float_or_none(x): 

433 if x is None: 

434 return None 

435 else: 

436 return float(x) 

437 

438 

439def int_or_none(x): 

440 if x is None: 

441 return None 

442 else: 

443 return int(x) 

444 

445 

446def int_or_g_tmin(x): 

447 if x is None: 

448 return g_tmin 

449 else: 

450 return int(x) 

451 

452 

453def int_or_g_tmax(x): 

454 if x is None: 

455 return g_tmax 

456 else: 

457 return int(x) 

458 

459 

460def tmin_or_none(x): 

461 if x == g_tmin: 

462 return None 

463 else: 

464 return x 

465 

466 

467def tmax_or_none(x): 

468 if x == g_tmax: 

469 return None 

470 else: 

471 return x 

472 

473 

474def time_or_none_to_str(x): 

475 if x is None: 

476 return '...'.ljust(17) 

477 else: 

478 return util.time_to_str(x) 

479 

480 

481def codes_to_str_abbreviated(codes, indent=' '): 

482 codes = [str(x) for x in codes] 

483 

484 if len(codes) > 20: 

485 scodes = '\n' + util.ewrap(codes[:10], indent=indent) \ 

486 + '\n%s[%i more]\n' % (indent, len(codes) - 20) \ 

487 + util.ewrap(codes[-10:], indent=' ') 

488 else: 

489 scodes = '\n' + util.ewrap(codes, indent=indent) \ 

490 if codes else '<none>' 

491 

492 return scodes 

493 

494 

495g_offset_time_unit_inv = 1000000000 

496g_offset_time_unit = 1.0 / g_offset_time_unit_inv 

497 

498 

499def tsplit(t): 

500 if t is None: 

501 return None, 0.0 

502 

503 t = util.to_time_float(t) 

504 if type(t) is float: 

505 t = round(t, 5) 

506 else: 

507 t = round(t, 9) 

508 

509 seconds = num.floor(t) 

510 offset = t - seconds 

511 return int(seconds), int(round(offset * g_offset_time_unit_inv)) 

512 

513 

514def tjoin(seconds, offset): 

515 if seconds is None: 

516 return None 

517 

518 return util.to_time_float(seconds) \ 

519 + util.to_time_float(offset*g_offset_time_unit) 

520 

521 

522tscale_min = 1 

523tscale_max = 365 * 24 * 3600 # last edge is one above 

524tscale_logbase = 20 

525 

526tscale_edges = [tscale_min] 

527while True: 

528 tscale_edges.append(tscale_edges[-1]*tscale_logbase) 

529 if tscale_edges[-1] >= tscale_max: 

530 break 

531 

532 

533tscale_edges = num.array(tscale_edges) 

534 

535 

536def tscale_to_kscale(tscale): 

537 

538 # 0 <= x < tscale_edges[1]: 0 

539 # tscale_edges[1] <= x < tscale_edges[2]: 1 

540 # ... 

541 # tscale_edges[len(tscale_edges)-1] <= x: len(tscale_edges) 

542 

543 return int(num.searchsorted(tscale_edges, tscale)) 

544 

545 

546@squirrel_content 

547class Station(Location): 

548 ''' 

549 A seismic station. 

550 ''' 

551 

552 codes = CodesNSL.T() 

553 

554 tmin = Timestamp.T(optional=True) 

555 tmax = Timestamp.T(optional=True) 

556 

557 description = Unicode.T(optional=True) 

558 

559 def __init__(self, **kwargs): 

560 kwargs['codes'] = CodesNSL(kwargs['codes']) 

561 Location.__init__(self, **kwargs) 

562 

563 @property 

564 def time_span(self): 

565 return (self.tmin, self.tmax) 

566 

567 def get_pyrocko_station(self): 

568 from pyrocko import model 

569 return model.Station(*self._get_pyrocko_station_args()) 

570 

571 def _get_pyrocko_station_args(self): 

572 return ( 

573 self.codes.network, 

574 self.codes.station, 

575 self.codes.location, 

576 self.lat, 

577 self.lon, 

578 self.elevation, 

579 self.depth, 

580 self.north_shift, 

581 self.east_shift) 

582 

583 

584class Sensor(Location): 

585 ''' 

586 Representation of a channel group. 

587 ''' 

588 

589 codes = CodesNSLCE.T() 

590 

591 tmin = Timestamp.T(optional=True) 

592 tmax = Timestamp.T(optional=True) 

593 

594 deltat = Float.T(optional=True) 

595 

596 @property 

597 def time_span(self): 

598 return (self.tmin, self.tmax) 

599 

600 def __init__(self, **kwargs): 

601 kwargs['codes'] = CodesNSLCE(kwargs['codes']) 

602 Location.__init__(self, **kwargs) 

603 

604 def _get_sensor_args(self): 

605 def getattr_rep(k): 

606 if k == 'codes': 

607 return self.codes.replace( 

608 channel=self.codes.channel[:-1] + '?') 

609 else: 

610 return getattr(self, k) 

611 

612 return tuple(getattr_rep(k) for k in Sensor.T.propnames) 

613 

614 @classmethod 

615 def from_channels(cls, channels): 

616 groups = defaultdict(list) 

617 for channel in channels: 

618 groups[channel._get_sensor_args()].append(channel) 

619 

620 return [ 

621 cls(**dict((k, v) for (k, v) in zip(cls.T.propnames, args))) 

622 for args, _ in groups.items()] 

623 

624 def _get_pyrocko_station_args(self): 

625 return ( 

626 self.codes.network, 

627 self.codes.station, 

628 self.codes.location, 

629 self.lat, 

630 self.lon, 

631 self.elevation, 

632 self.depth, 

633 self.north_shift, 

634 self.east_shift) 

635 

636 

637@squirrel_content 

638class Channel(Sensor): 

639 ''' 

640 A channel of a seismic station. 

641 ''' 

642 

643 dip = Float.T(optional=True) 

644 azimuth = Float.T(optional=True) 

645 

646 @classmethod 

647 def from_channels(cls, channels): 

648 raise NotImplementedError() 

649 

650 def get_pyrocko_channel(self): 

651 from pyrocko import model 

652 return model.Channel(*self._get_pyrocko_channel_args()) 

653 

654 def _get_pyrocko_channel_args(self): 

655 return ( 

656 self.codes.channel, 

657 self.azimuth, 

658 self.dip) 

659 

660 

661observational_quantities = [ 

662 'acceleration', 'velocity', 'displacement', 'pressure', 'rotation', 

663 'temperature'] 

664 

665 

666technical_quantities = [ 

667 'voltage', 'counts'] 

668 

669 

670class QuantityType(StringChoice): 

671 ''' 

672 Choice of observational or technical quantity. 

673 

674 SI units are used for all quantities, where applicable. 

675 ''' 

676 choices = observational_quantities + technical_quantities 

677 

678 

679class ResponseStage(Object): 

680 ''' 

681 Representation of a response stage. 

682 

683 Components of a seismic recording system are represented as a sequence of 

684 response stages, e.g. sensor, pre-amplifier, digitizer, digital 

685 downsampling. 

686 ''' 

687 input_quantity = QuantityType.T(optional=True) 

688 input_sample_rate = Float.T(optional=True) 

689 output_quantity = QuantityType.T(optional=True) 

690 output_sample_rate = Float.T(optional=True) 

691 elements = List.T(FrequencyResponse.T()) 

692 log = List.T(Tuple.T(3, String.T())) 

693 

694 @property 

695 def stage_type(self): 

696 if self.input_quantity in observational_quantities \ 

697 and self.output_quantity in observational_quantities: 

698 return 'conversion' 

699 

700 if self.input_quantity in observational_quantities \ 

701 and self.output_quantity == 'voltage': 

702 return 'sensor' 

703 

704 elif self.input_quantity == 'voltage' \ 

705 and self.output_quantity == 'voltage': 

706 return 'electronics' 

707 

708 elif self.input_quantity == 'voltage' \ 

709 and self.output_quantity == 'counts': 

710 return 'digitizer' 

711 

712 elif self.input_quantity == 'counts' \ 

713 and self.output_quantity == 'counts' \ 

714 and self.input_sample_rate != self.output_sample_rate: 

715 return 'decimation' 

716 

717 elif self.input_quantity in observational_quantities \ 

718 and self.output_quantity == 'counts': 

719 return 'combination' 

720 

721 else: 

722 return 'unknown' 

723 

724 @property 

725 def summary(self): 

726 irate = self.input_sample_rate 

727 orate = self.output_sample_rate 

728 factor = None 

729 if irate and orate: 

730 factor = irate / orate 

731 return 'ResponseStage, ' + ( 

732 '%s%s => %s%s%s' % ( 

733 self.input_quantity or '?', 

734 ' @ %g Hz' % irate if irate else '', 

735 self.output_quantity or '?', 

736 ' @ %g Hz' % orate if orate else '', 

737 ' :%g' % factor if factor else '') 

738 ) 

739 

740 def get_effective(self): 

741 return MultiplyResponse(responses=list(self.elements)) 

742 

743 

744D = 'displacement' 

745V = 'velocity' 

746A = 'acceleration' 

747 

748g_converters = { 

749 (V, D): IntegrationResponse(1), 

750 (A, D): IntegrationResponse(2), 

751 (D, V): DifferentiationResponse(1), 

752 (A, V): IntegrationResponse(1), 

753 (D, A): DifferentiationResponse(2), 

754 (V, A): DifferentiationResponse(1)} 

755 

756 

757def response_converters(input_quantity, output_quantity): 

758 if input_quantity is None or input_quantity == output_quantity: 

759 return [] 

760 

761 if output_quantity is None: 

762 raise ConversionError('Unspecified target quantity.') 

763 

764 try: 

765 return [g_converters[input_quantity, output_quantity]] 

766 

767 except KeyError: 

768 raise ConversionError('No rule to convert from "%s" to "%s".' % ( 

769 input_quantity, output_quantity)) 

770 

771 

772@squirrel_content 

773class Response(Object): 

774 ''' 

775 The instrument response of a seismic station channel. 

776 ''' 

777 

778 codes = CodesNSLCE.T() 

779 tmin = Timestamp.T(optional=True) 

780 tmax = Timestamp.T(optional=True) 

781 

782 stages = List.T(ResponseStage.T()) 

783 checkpoints = List.T(FrequencyResponseCheckpoint.T()) 

784 

785 deltat = Float.T(optional=True) 

786 log = List.T(Tuple.T(3, String.T())) 

787 

788 def __init__(self, **kwargs): 

789 kwargs['codes'] = CodesNSLCE(kwargs['codes']) 

790 Object.__init__(self, **kwargs) 

791 

792 @property 

793 def time_span(self): 

794 return (self.tmin, self.tmax) 

795 

796 @property 

797 def nstages(self): 

798 return len(self.stages) 

799 

800 @property 

801 def input_quantity(self): 

802 return self.stages[0].input_quantity if self.stages else None 

803 

804 @property 

805 def output_quantity(self): 

806 return self.stages[-1].input_quantity if self.stages else None 

807 

808 @property 

809 def output_sample_rate(self): 

810 return self.stages[-1].output_sample_rate if self.stages else None 

811 

812 @property 

813 def stages_summary(self): 

814 def grouped(xs): 

815 xs = list(xs) 

816 g = [] 

817 for i in range(len(xs)): 

818 g.append(xs[i]) 

819 if i+1 < len(xs) and xs[i+1] != xs[i]: 

820 yield g 

821 g = [] 

822 

823 if g: 

824 yield g 

825 

826 return '+'.join( 

827 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '') 

828 for g in grouped(stage.stage_type for stage in self.stages)) 

829 

830 @property 

831 def summary(self): 

832 orate = self.output_sample_rate 

833 return '%s %-16s %s' % ( 

834 self.__class__.__name__, self.str_codes, self.str_time_span) \ 

835 + ', ' + ', '.join(( 

836 '%s => %s' % ( 

837 self.input_quantity or '?', self.output_quantity or '?') 

838 + (' @ %g Hz' % orate if orate else ''), 

839 self.stages_summary, 

840 )) 

841 

842 def get_effective(self, input_quantity=None): 

843 elements = response_converters(input_quantity, self.input_quantity) 

844 

845 elements.extend( 

846 stage.get_effective() for stage in self.stages) 

847 

848 return MultiplyResponse(responses=simplify_responses(elements)) 

849 

850 

851@squirrel_content 

852class Event(Object): 

853 ''' 

854 A seismic event. 

855 ''' 

856 

857 name = String.T(optional=True) 

858 time = Timestamp.T() 

859 duration = Float.T(optional=True) 

860 

861 lat = Float.T() 

862 lon = Float.T() 

863 depth = Float.T(optional=True) 

864 

865 magnitude = Float.T(optional=True) 

866 

867 def get_pyrocko_event(self): 

868 from pyrocko import model 

869 model.Event( 

870 name=self.name, 

871 time=self.time, 

872 lat=self.lat, 

873 lon=self.lon, 

874 depth=self.depth, 

875 magnitude=self.magnitude, 

876 duration=self.duration) 

877 

878 @property 

879 def time_span(self): 

880 return (self.time, self.time) 

881 

882 

883def ehash(s): 

884 return hashlib.sha1(s.encode('utf8')).hexdigest() 

885 

886 

887def random_name(n=8): 

888 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii') 

889 

890 

891@squirrel_content 

892class WaveformPromise(Object): 

893 ''' 

894 Information about a waveform potentially downloadable from a remote site. 

895 

896 In the Squirrel framework, waveform promises are used to permit download of 

897 selected waveforms from a remote site. They are typically generated by 

898 calls to 

899 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`. 

900 Waveform promises are inserted and indexed in the database similar to 

901 normal waveforms. When processing a waveform query, e.g. from 

902 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local 

903 waveform is available for the queried time span, a matching promise can be 

904 resolved, i.e. an attempt is made to download the waveform from the remote 

905 site. The promise is removed after the download attempt (except when a 

906 network error occurs). This prevents Squirrel from making unnecessary 

907 queries for waveforms missing at the remote site. 

908 ''' 

909 

910 codes = CodesNSLCE.T() 

911 tmin = Timestamp.T() 

912 tmax = Timestamp.T() 

913 

914 deltat = Float.T(optional=True) 

915 

916 source_hash = String.T() 

917 

918 def __init__(self, **kwargs): 

919 kwargs['codes'] = CodesNSLCE(kwargs['codes']) 

920 Object.__init__(self, **kwargs) 

921 

922 @property 

923 def time_span(self): 

924 return (self.tmin, self.tmax) 

925 

926 

927class InvalidWaveform(Exception): 

928 pass 

929 

930 

931class WaveformOrder(Object): 

932 ''' 

933 Waveform request information. 

934 ''' 

935 

936 source_id = String.T() 

937 codes = CodesNSLCE.T() 

938 deltat = Float.T() 

939 tmin = Timestamp.T() 

940 tmax = Timestamp.T() 

941 gaps = List.T(Tuple.T(2, Timestamp.T())) 

942 

943 @property 

944 def client(self): 

945 return self.source_id.split(':')[1] 

946 

947 def describe(self, site='?'): 

948 return '%s:%s %s [%s - %s]' % ( 

949 self.client, site, str(self.codes), 

950 util.time_to_str(self.tmin), util.time_to_str(self.tmax)) 

951 

952 def validate(self, tr): 

953 if tr.ydata.size == 0: 

954 raise InvalidWaveform( 

955 'waveform with zero data samples') 

956 

957 if tr.deltat != self.deltat: 

958 raise InvalidWaveform( 

959 'incorrect sampling interval - waveform: %g s, ' 

960 'meta-data: %g s' % ( 

961 tr.deltat, self.deltat)) 

962 

963 if not num.all(num.isfinite(tr.ydata)): 

964 raise InvalidWaveform('waveform has NaN values') 

965 

966 

967def order_summary(orders): 

968 codes_list = sorted(set(order.codes for order in orders)) 

969 if len(codes_list) > 3: 

970 return '%i order%s: %s - %s' % ( 

971 len(orders), 

972 util.plural_s(orders), 

973 str(codes_list[0]), 

974 str(codes_list[1])) 

975 

976 else: 

977 return '%i order%s: %s' % ( 

978 len(orders), 

979 util.plural_s(orders), 

980 ', '.join(str(codes) for codes in codes_list)) 

981 

982 

983class Nut(Object): 

984 ''' 

985 Index entry referencing an elementary piece of content. 

986 

987 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common 

988 meta-information about individual pieces of waveforms, stations, channels, 

989 etc. together with the information where it was found or generated. 

990 ''' 

991 

992 file_path = String.T(optional=True) 

993 file_format = String.T(optional=True) 

994 file_mtime = Timestamp.T(optional=True) 

995 file_size = Int.T(optional=True) 

996 

997 file_segment = Int.T(optional=True) 

998 file_element = Int.T(optional=True) 

999 

1000 kind_id = Int.T() 

1001 codes = Codes.T() 

1002 

1003 tmin_seconds = Int.T(default=0) 

1004 tmin_offset = Int.T(default=0, optional=True) 

1005 tmax_seconds = Int.T(default=0) 

1006 tmax_offset = Int.T(default=0, optional=True) 

1007 

1008 deltat = Float.T(default=0.0) 

1009 

1010 content = Any.T(optional=True) 

1011 raw_content = Dict.T(String.T(), Any.T()) 

1012 

1013 content_in_db = False 

1014 

1015 def __init__( 

1016 self, 

1017 file_path=None, 

1018 file_format=None, 

1019 file_mtime=None, 

1020 file_size=None, 

1021 file_segment=None, 

1022 file_element=None, 

1023 kind_id=0, 

1024 codes=CodesX(''), 

1025 tmin_seconds=None, 

1026 tmin_offset=0, 

1027 tmax_seconds=None, 

1028 tmax_offset=0, 

1029 deltat=None, 

1030 content=None, 

1031 raw_content=None, 

1032 tmin=None, 

1033 tmax=None, 

1034 values_nocheck=None): 

1035 

1036 if values_nocheck is not None: 

1037 (self.file_path, self.file_format, self.file_mtime, 

1038 self.file_size, 

1039 self.file_segment, self.file_element, 

1040 self.kind_id, codes_safe_str, 

1041 self.tmin_seconds, self.tmin_offset, 

1042 self.tmax_seconds, self.tmax_offset, 

1043 self.deltat) = values_nocheck 

1044 

1045 self.codes = to_codes_simple(self.kind_id, codes_safe_str) 

1046 self.raw_content = {} 

1047 self.content = None 

1048 else: 

1049 if tmin is not None: 

1050 tmin_seconds, tmin_offset = tsplit(tmin) 

1051 

1052 if tmax is not None: 

1053 tmax_seconds, tmax_offset = tsplit(tmax) 

1054 

1055 self.kind_id = int(kind_id) 

1056 self.codes = codes 

1057 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1058 self.tmin_offset = int(tmin_offset) 

1059 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1060 self.tmax_offset = int(tmax_offset) 

1061 self.deltat = float_or_none(deltat) 

1062 self.file_path = str_or_none(file_path) 

1063 self.file_segment = int_or_none(file_segment) 

1064 self.file_element = int_or_none(file_element) 

1065 self.file_format = str_or_none(file_format) 

1066 self.file_mtime = float_or_none(file_mtime) 

1067 self.file_size = int_or_none(file_size) 

1068 self.content = content 

1069 if raw_content is None: 

1070 self.raw_content = {} 

1071 else: 

1072 self.raw_content = raw_content 

1073 

1074 Object.__init__(self, init_props=False) 

1075 

1076 def __eq__(self, other): 

1077 return (isinstance(other, Nut) and 

1078 self.equality_values == other.equality_values) 

1079 

1080 def hash(self): 

1081 return ehash(','.join(str(x) for x in self.key)) 

1082 

1083 def __ne__(self, other): 

1084 return not (self == other) 

1085 

1086 def get_io_backend(self): 

1087 from . import io 

1088 return io.get_backend(self.file_format) 

1089 

1090 def file_modified(self): 

1091 return self.get_io_backend().get_stats(self.file_path) \ 

1092 != (self.file_mtime, self.file_size) 

1093 

1094 @property 

1095 def dkey(self): 

1096 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes) 

1097 

1098 @property 

1099 def key(self): 

1100 return ( 

1101 self.file_path, 

1102 self.file_segment, 

1103 self.file_element, 

1104 self.file_mtime) 

1105 

1106 @property 

1107 def equality_values(self): 

1108 return ( 

1109 self.file_segment, self.file_element, 

1110 self.kind_id, self.codes, 

1111 self.tmin_seconds, self.tmin_offset, 

1112 self.tmax_seconds, self.tmax_offset, self.deltat) 

1113 

1114 @property 

1115 def tmin(self): 

1116 return tjoin(self.tmin_seconds, self.tmin_offset) 

1117 

1118 @tmin.setter 

1119 def tmin(self, t): 

1120 self.tmin_seconds, self.tmin_offset = tsplit(t) 

1121 

1122 @property 

1123 def tmax(self): 

1124 return tjoin(self.tmax_seconds, self.tmax_offset) 

1125 

1126 @tmax.setter 

1127 def tmax(self, t): 

1128 self.tmax_seconds, self.tmax_offset = tsplit(t) 

1129 

1130 @property 

1131 def kscale(self): 

1132 if self.tmin_seconds is None or self.tmax_seconds is None: 

1133 return 0 

1134 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds) 

1135 

1136 @property 

1137 def waveform_kwargs(self): 

1138 network, station, location, channel, extra = self.codes 

1139 

1140 return dict( 

1141 network=network, 

1142 station=station, 

1143 location=location, 

1144 channel=channel, 

1145 extra=extra, 

1146 tmin=self.tmin, 

1147 tmax=self.tmax, 

1148 deltat=self.deltat) 

1149 

1150 @property 

1151 def waveform_promise_kwargs(self): 

1152 return dict( 

1153 codes=self.codes, 

1154 tmin=self.tmin, 

1155 tmax=self.tmax, 

1156 deltat=self.deltat) 

1157 

1158 @property 

1159 def station_kwargs(self): 

1160 network, station, location = self.codes 

1161 return dict( 

1162 codes=self.codes, 

1163 tmin=tmin_or_none(self.tmin), 

1164 tmax=tmax_or_none(self.tmax)) 

1165 

1166 @property 

1167 def channel_kwargs(self): 

1168 network, station, location, channel, extra = self.codes 

1169 return dict( 

1170 codes=self.codes, 

1171 tmin=tmin_or_none(self.tmin), 

1172 tmax=tmax_or_none(self.tmax), 

1173 deltat=self.deltat) 

1174 

1175 @property 

1176 def response_kwargs(self): 

1177 return dict( 

1178 codes=self.codes, 

1179 tmin=tmin_or_none(self.tmin), 

1180 tmax=tmax_or_none(self.tmax), 

1181 deltat=self.deltat) 

1182 

1183 @property 

1184 def event_kwargs(self): 

1185 return dict( 

1186 name=self.codes, 

1187 time=self.tmin, 

1188 duration=(self.tmax - self.tmin) or None) 

1189 

1190 @property 

1191 def trace_kwargs(self): 

1192 network, station, location, channel, extra = self.codes 

1193 

1194 return dict( 

1195 network=network, 

1196 station=station, 

1197 location=location, 

1198 channel=channel, 

1199 extra=extra, 

1200 tmin=self.tmin, 

1201 tmax=self.tmax-self.deltat, 

1202 deltat=self.deltat) 

1203 

1204 @property 

1205 def dummy_trace(self): 

1206 return DummyTrace(self) 

1207 

1208 @property 

1209 def summary(self): 

1210 if self.tmin == self.tmax: 

1211 ts = util.time_to_str(self.tmin) 

1212 else: 

1213 ts = '%s - %s' % ( 

1214 util.time_to_str(self.tmin), 

1215 util.time_to_str(self.tmax)) 

1216 

1217 return ' '.join(( 

1218 ('%s,' % to_kind(self.kind_id)).ljust(9), 

1219 ('%s,' % str(self.codes)).ljust(18), 

1220 ts)) 

1221 

1222 

1223def make_waveform_nut(**kwargs): 

1224 return Nut(kind_id=WAVEFORM, **kwargs) 

1225 

1226 

1227def make_waveform_promise_nut(**kwargs): 

1228 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs) 

1229 

1230 

1231def make_station_nut(**kwargs): 

1232 return Nut(kind_id=STATION, **kwargs) 

1233 

1234 

1235def make_channel_nut(**kwargs): 

1236 return Nut(kind_id=CHANNEL, **kwargs) 

1237 

1238 

1239def make_response_nut(**kwargs): 

1240 return Nut(kind_id=RESPONSE, **kwargs) 

1241 

1242 

1243def make_event_nut(**kwargs): 

1244 return Nut(kind_id=EVENT, **kwargs) 

1245 

1246 

1247def group_channels(nuts): 

1248 by_ansl = {} 

1249 for nut in nuts: 

1250 if nut.kind_id != CHANNEL: 

1251 continue 

1252 

1253 ansl = nut.codes[:4] 

1254 

1255 if ansl not in by_ansl: 

1256 by_ansl[ansl] = {} 

1257 

1258 group = by_ansl[ansl] 

1259 

1260 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax 

1261 

1262 if k not in group: 

1263 group[k] = set() 

1264 

1265 group.add(nut.codes[4]) 

1266 

1267 return by_ansl 

1268 

1269 

1270class DummyTrace(object): 

1271 

1272 def __init__(self, nut): 

1273 self.nut = nut 

1274 self.codes = nut.codes 

1275 self.meta = {} 

1276 

1277 @property 

1278 def tmin(self): 

1279 return self.nut.tmin 

1280 

1281 @property 

1282 def tmax(self): 

1283 return self.nut.tmax 

1284 

1285 @property 

1286 def deltat(self): 

1287 return self.nut.deltat 

1288 

1289 @property 

1290 def nslc_id(self): 

1291 return self.codes.nslc 

1292 

1293 @property 

1294 def network(self): 

1295 return self.codes.network 

1296 

1297 @property 

1298 def station(self): 

1299 return self.codes.station 

1300 

1301 @property 

1302 def location(self): 

1303 return self.codes.location 

1304 

1305 @property 

1306 def channel(self): 

1307 return self.codes.channel 

1308 

1309 @property 

1310 def extra(self): 

1311 return self.codes.extra 

1312 

1313 def overlaps(self, tmin, tmax): 

1314 return not (tmax < self.nut.tmin or self.nut.tmax < tmin) 

1315 

1316 

1317def duration_to_str(t): 

1318 if t > 24*3600: 

1319 return '%gd' % (t / (24.*3600.)) 

1320 elif t > 3600: 

1321 return '%gh' % (t / 3600.) 

1322 elif t > 60: 

1323 return '%gm' % (t / 60.) 

1324 else: 

1325 return '%gs' % t 

1326 

1327 

1328class Coverage(Object): 

1329 ''' 

1330 Information about times covered by a waveform or other time series data. 

1331 ''' 

1332 kind_id = Int.T( 

1333 help='Content type.') 

1334 pattern = Codes.T( 

1335 help='The codes pattern in the request, which caused this entry to ' 

1336 'match.') 

1337 codes = Codes.T( 

1338 help='NSLCE or NSL codes identifier of the time series.') 

1339 deltat = Float.T( 

1340 help='Sampling interval [s]', 

1341 optional=True) 

1342 tmin = Timestamp.T( 

1343 help='Global start time of time series.', 

1344 optional=True) 

1345 tmax = Timestamp.T( 

1346 help='Global end time of time series.', 

1347 optional=True) 

1348 changes = List.T( 

1349 Tuple.T(2, Any.T()), 

1350 help='List of change points, with entries of the form ' 

1351 '``(time, count)``, where a ``count`` of zero indicates start of ' 

1352 'a gap, a value of 1 start of normal data coverage and a higher ' 

1353 'value duplicate or redundant data coverage.') 

1354 

1355 @classmethod 

1356 def from_values(cls, args): 

1357 pattern, codes, deltat, tmin, tmax, changes, kind_id = args 

1358 return cls( 

1359 kind_id=kind_id, 

1360 pattern=pattern, 

1361 codes=codes, 

1362 deltat=deltat, 

1363 tmin=tmin, 

1364 tmax=tmax, 

1365 changes=changes) 

1366 

1367 @property 

1368 def summary(self): 

1369 ts = '%s - %s,' % ( 

1370 util.time_to_str(self.tmin), 

1371 util.time_to_str(self.tmax)) 

1372 

1373 srate = self.sample_rate 

1374 

1375 return ' '.join(( 

1376 ('%s,' % to_kind(self.kind_id)).ljust(9), 

1377 ('%s,' % str(self.codes)).ljust(18), 

1378 ts, 

1379 '%10.3g,' % srate if srate else '', 

1380 '%4i' % len(self.changes), 

1381 '%s' % duration_to_str(self.total))) 

1382 

1383 @property 

1384 def sample_rate(self): 

1385 if self.deltat is None: 

1386 return None 

1387 elif self.deltat == 0.0: 

1388 return 0.0 

1389 else: 

1390 return 1.0 / self.deltat 

1391 

1392 @property 

1393 def labels(self): 

1394 srate = self.sample_rate 

1395 return ( 

1396 ('%s' % str(self.codes)), 

1397 '%.3g' % srate if srate else '') 

1398 

1399 @property 

1400 def total(self): 

1401 total_t = None 

1402 for tmin, tmax, _ in self.iter_spans(): 

1403 total_t = (total_t or 0.0) + (tmax - tmin) 

1404 

1405 return total_t 

1406 

1407 def iter_spans(self): 

1408 last = None 

1409 for (t, count) in self.changes: 

1410 if last is not None: 

1411 last_t, last_count = last 

1412 if last_count > 0: 

1413 yield last_t, t, last_count 

1414 

1415 last = (t, count) 

1416 

1417 

1418__all__ = [ 

1419 'to_codes', 

1420 'to_codes_guess', 

1421 'to_codes_simple', 

1422 'to_kind', 

1423 'to_kinds', 

1424 'to_kind_id', 

1425 'to_kind_ids', 

1426 'CodesError', 

1427 'Codes', 

1428 'CodesNSLCE', 

1429 'CodesNSL', 

1430 'CodesX', 

1431 'Station', 

1432 'Channel', 

1433 'Sensor', 

1434 'Response', 

1435 'Nut', 

1436 'Coverage', 

1437 'WaveformPromise', 

1438]