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 

385# derived list class to enable detection of already preprocessed codes patterns 

386class codes_patterns_list(list): 

387 pass 

388 

389 

390def codes_patterns_for_kind(kind_id, codes): 

391 if isinstance(codes, codes_patterns_list): 

392 return codes 

393 

394 if isinstance(codes, list): 

395 lcodes = codes_patterns_list() 

396 for sc in codes: 

397 lcodes.extend(codes_patterns_for_kind(kind_id, sc)) 

398 

399 return lcodes 

400 

401 codes = to_codes(kind_id, codes) 

402 

403 lcodes = codes_patterns_list() 

404 lcodes.append(codes) 

405 if kind_id == STATION: 

406 lcodes.append(codes.replace(location='[*]')) 

407 

408 return lcodes 

409 

410 

411g_content_kind_ids = ( 

412 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT, 

413 WAVEFORM_PROMISE) = range(len(g_content_kinds)) 

414 

415 

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

417 

418 

419try: 

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

421except Exception: 

422 g_tmin_queries = g_tmin 

423 

424 

425def to_kind(kind_id): 

426 return g_content_kinds[kind_id] 

427 

428 

429def to_kinds(kind_ids): 

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

431 

432 

433def to_kind_id(kind): 

434 return g_content_kinds.index(kind) 

435 

436 

437def to_kind_ids(kinds): 

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

439 

440 

441g_kind_mask_all = 0xff 

442 

443 

444def to_kind_mask(kinds): 

445 if kinds: 

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

447 else: 

448 return g_kind_mask_all 

449 

450 

451def str_or_none(x): 

452 if x is None: 

453 return None 

454 else: 

455 return str(x) 

456 

457 

458def float_or_none(x): 

459 if x is None: 

460 return None 

461 else: 

462 return float(x) 

463 

464 

465def int_or_none(x): 

466 if x is None: 

467 return None 

468 else: 

469 return int(x) 

470 

471 

472def int_or_g_tmin(x): 

473 if x is None: 

474 return g_tmin 

475 else: 

476 return int(x) 

477 

478 

479def int_or_g_tmax(x): 

480 if x is None: 

481 return g_tmax 

482 else: 

483 return int(x) 

484 

485 

486def tmin_or_none(x): 

487 if x == g_tmin: 

488 return None 

489 else: 

490 return x 

491 

492 

493def tmax_or_none(x): 

494 if x == g_tmax: 

495 return None 

496 else: 

497 return x 

498 

499 

500def time_or_none_to_str(x): 

501 if x is None: 

502 return '...'.ljust(17) 

503 else: 

504 return util.time_to_str(x) 

505 

506 

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

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

509 

510 if len(codes) > 20: 

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

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

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

514 else: 

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

516 if codes else '<none>' 

517 

518 return scodes 

519 

520 

521g_offset_time_unit_inv = 1000000000 

522g_offset_time_unit = 1.0 / g_offset_time_unit_inv 

523 

524 

525def tsplit(t): 

526 if t is None: 

527 return None, 0 

528 

529 t = util.to_time_float(t) 

530 if type(t) is float: 

531 t = round(t, 5) 

532 else: 

533 t = round(t, 9) 

534 

535 seconds = num.floor(t) 

536 offset = t - seconds 

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

538 

539 

540def tjoin(seconds, offset): 

541 if seconds is None: 

542 return None 

543 

544 return util.to_time_float(seconds) \ 

545 + util.to_time_float(offset*g_offset_time_unit) 

546 

547 

548tscale_min = 1 

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

550tscale_logbase = 20 

551 

552tscale_edges = [tscale_min] 

553while True: 

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

555 if tscale_edges[-1] >= tscale_max: 

556 break 

557 

558 

559tscale_edges = num.array(tscale_edges) 

560 

561 

562def tscale_to_kscale(tscale): 

563 

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

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

566 # ... 

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

568 

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

570 

571 

572@squirrel_content 

573class Station(Location): 

574 ''' 

575 A seismic station. 

576 ''' 

577 

578 codes = CodesNSL.T() 

579 

580 tmin = Timestamp.T(optional=True) 

581 tmax = Timestamp.T(optional=True) 

582 

583 description = Unicode.T(optional=True) 

584 

585 def __init__(self, **kwargs): 

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

587 Location.__init__(self, **kwargs) 

588 

589 @property 

590 def time_span(self): 

591 return (self.tmin, self.tmax) 

592 

593 def get_pyrocko_station(self): 

594 from pyrocko import model 

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

596 

597 def _get_pyrocko_station_args(self): 

598 return ( 

599 self.codes.network, 

600 self.codes.station, 

601 self.codes.location, 

602 self.lat, 

603 self.lon, 

604 self.elevation, 

605 self.depth, 

606 self.north_shift, 

607 self.east_shift) 

608 

609 

610class Sensor(Location): 

611 ''' 

612 Representation of a channel group. 

613 ''' 

614 

615 codes = CodesNSLCE.T() 

616 

617 tmin = Timestamp.T(optional=True) 

618 tmax = Timestamp.T(optional=True) 

619 

620 deltat = Float.T(optional=True) 

621 

622 @property 

623 def time_span(self): 

624 return (self.tmin, self.tmax) 

625 

626 def __init__(self, **kwargs): 

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

628 Location.__init__(self, **kwargs) 

629 

630 def _get_sensor_args(self): 

631 def getattr_rep(k): 

632 if k == 'codes': 

633 return self.codes.replace( 

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

635 else: 

636 return getattr(self, k) 

637 

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

639 

640 @classmethod 

641 def from_channels(cls, channels): 

642 groups = defaultdict(list) 

643 for channel in channels: 

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

645 

646 return [ 

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

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

649 

650 def _get_pyrocko_station_args(self): 

651 return ( 

652 self.codes.network, 

653 self.codes.station, 

654 self.codes.location, 

655 self.lat, 

656 self.lon, 

657 self.elevation, 

658 self.depth, 

659 self.north_shift, 

660 self.east_shift) 

661 

662 

663@squirrel_content 

664class Channel(Sensor): 

665 ''' 

666 A channel of a seismic station. 

667 ''' 

668 

669 dip = Float.T(optional=True) 

670 azimuth = Float.T(optional=True) 

671 

672 @classmethod 

673 def from_channels(cls, channels): 

674 raise NotImplementedError() 

675 

676 def get_pyrocko_channel(self): 

677 from pyrocko import model 

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

679 

680 def _get_pyrocko_channel_args(self): 

681 return ( 

682 self.codes.channel, 

683 self.azimuth, 

684 self.dip) 

685 

686 

687observational_quantities = [ 

688 'acceleration', 'velocity', 'displacement', 'pressure', 'rotation', 

689 'temperature'] 

690 

691 

692technical_quantities = [ 

693 'voltage', 'counts'] 

694 

695 

696class QuantityType(StringChoice): 

697 ''' 

698 Choice of observational or technical quantity. 

699 

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

701 ''' 

702 choices = observational_quantities + technical_quantities 

703 

704 

705class ResponseStage(Object): 

706 ''' 

707 Representation of a response stage. 

708 

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

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

711 downsampling. 

712 ''' 

713 input_quantity = QuantityType.T(optional=True) 

714 input_sample_rate = Float.T(optional=True) 

715 output_quantity = QuantityType.T(optional=True) 

716 output_sample_rate = Float.T(optional=True) 

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

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

719 

720 @property 

721 def stage_type(self): 

722 if self.input_quantity in observational_quantities \ 

723 and self.output_quantity in observational_quantities: 

724 return 'conversion' 

725 

726 if self.input_quantity in observational_quantities \ 

727 and self.output_quantity == 'voltage': 

728 return 'sensor' 

729 

730 elif self.input_quantity == 'voltage' \ 

731 and self.output_quantity == 'voltage': 

732 return 'electronics' 

733 

734 elif self.input_quantity == 'voltage' \ 

735 and self.output_quantity == 'counts': 

736 return 'digitizer' 

737 

738 elif self.input_quantity == 'counts' \ 

739 and self.output_quantity == 'counts' \ 

740 and self.input_sample_rate != self.output_sample_rate: 

741 return 'decimation' 

742 

743 elif self.input_quantity in observational_quantities \ 

744 and self.output_quantity == 'counts': 

745 return 'combination' 

746 

747 else: 

748 return 'unknown' 

749 

750 @property 

751 def summary(self): 

752 irate = self.input_sample_rate 

753 orate = self.output_sample_rate 

754 factor = None 

755 if irate and orate: 

756 factor = irate / orate 

757 return 'ResponseStage, ' + ( 

758 '%s%s => %s%s%s' % ( 

759 self.input_quantity or '?', 

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

761 self.output_quantity or '?', 

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

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

764 ) 

765 

766 def get_effective(self): 

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

768 

769 

770D = 'displacement' 

771V = 'velocity' 

772A = 'acceleration' 

773 

774g_converters = { 

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

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

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

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

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

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

781 

782 

783def response_converters(input_quantity, output_quantity): 

784 if input_quantity is None or input_quantity == output_quantity: 

785 return [] 

786 

787 if output_quantity is None: 

788 raise ConversionError('Unspecified target quantity.') 

789 

790 try: 

791 return [g_converters[input_quantity, output_quantity]] 

792 

793 except KeyError: 

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

795 input_quantity, output_quantity)) 

796 

797 

798@squirrel_content 

799class Response(Object): 

800 ''' 

801 The instrument response of a seismic station channel. 

802 ''' 

803 

804 codes = CodesNSLCE.T() 

805 tmin = Timestamp.T(optional=True) 

806 tmax = Timestamp.T(optional=True) 

807 

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

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

810 

811 deltat = Float.T(optional=True) 

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

813 

814 def __init__(self, **kwargs): 

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

816 Object.__init__(self, **kwargs) 

817 

818 @property 

819 def time_span(self): 

820 return (self.tmin, self.tmax) 

821 

822 @property 

823 def nstages(self): 

824 return len(self.stages) 

825 

826 @property 

827 def input_quantity(self): 

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

829 

830 @property 

831 def output_quantity(self): 

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

833 

834 @property 

835 def output_sample_rate(self): 

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

837 

838 @property 

839 def stages_summary(self): 

840 def grouped(xs): 

841 xs = list(xs) 

842 g = [] 

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

844 g.append(xs[i]) 

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

846 yield g 

847 g = [] 

848 

849 if g: 

850 yield g 

851 

852 return '+'.join( 

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

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

855 

856 @property 

857 def summary(self): 

858 orate = self.output_sample_rate 

859 return '%s %-16s %s' % ( 

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

861 + ', ' + ', '.join(( 

862 '%s => %s' % ( 

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

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

865 self.stages_summary, 

866 )) 

867 

868 def get_effective(self, input_quantity=None): 

869 try: 

870 elements = response_converters(input_quantity, self.input_quantity) 

871 except ConversionError as e: 

872 raise ConversionError(str(e) + ' (%s)' % self.summary) 

873 

874 elements.extend( 

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

876 

877 return MultiplyResponse(responses=simplify_responses(elements)) 

878 

879 

880@squirrel_content 

881class Event(Object): 

882 ''' 

883 A seismic event. 

884 ''' 

885 

886 name = String.T(optional=True) 

887 time = Timestamp.T() 

888 duration = Float.T(optional=True) 

889 

890 lat = Float.T() 

891 lon = Float.T() 

892 depth = Float.T(optional=True) 

893 

894 magnitude = Float.T(optional=True) 

895 

896 def get_pyrocko_event(self): 

897 from pyrocko import model 

898 model.Event( 

899 name=self.name, 

900 time=self.time, 

901 lat=self.lat, 

902 lon=self.lon, 

903 depth=self.depth, 

904 magnitude=self.magnitude, 

905 duration=self.duration) 

906 

907 @property 

908 def time_span(self): 

909 return (self.time, self.time) 

910 

911 

912def ehash(s): 

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

914 

915 

916def random_name(n=8): 

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

918 

919 

920@squirrel_content 

921class WaveformPromise(Object): 

922 ''' 

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

924 

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

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

927 calls to 

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

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

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

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

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

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

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

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

936 queries for waveforms missing at the remote site. 

937 ''' 

938 

939 codes = CodesNSLCE.T() 

940 tmin = Timestamp.T() 

941 tmax = Timestamp.T() 

942 

943 deltat = Float.T(optional=True) 

944 

945 source_hash = String.T() 

946 

947 def __init__(self, **kwargs): 

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

949 Object.__init__(self, **kwargs) 

950 

951 @property 

952 def time_span(self): 

953 return (self.tmin, self.tmax) 

954 

955 

956class InvalidWaveform(Exception): 

957 pass 

958 

959 

960class WaveformOrder(Object): 

961 ''' 

962 Waveform request information. 

963 ''' 

964 

965 source_id = String.T() 

966 codes = CodesNSLCE.T() 

967 deltat = Float.T() 

968 tmin = Timestamp.T() 

969 tmax = Timestamp.T() 

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

971 

972 @property 

973 def client(self): 

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

975 

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

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

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

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

980 

981 def validate(self, tr): 

982 if tr.ydata.size == 0: 

983 raise InvalidWaveform( 

984 'waveform with zero data samples') 

985 

986 if tr.deltat != self.deltat: 

987 raise InvalidWaveform( 

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

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

990 tr.deltat, self.deltat)) 

991 

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

993 raise InvalidWaveform('waveform has NaN values') 

994 

995 

996def order_summary(orders): 

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

998 if len(codes_list) > 3: 

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

1000 len(orders), 

1001 util.plural_s(orders), 

1002 str(codes_list[0]), 

1003 str(codes_list[1])) 

1004 

1005 else: 

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

1007 len(orders), 

1008 util.plural_s(orders), 

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

1010 

1011 

1012class Nut(Object): 

1013 ''' 

1014 Index entry referencing an elementary piece of content. 

1015 

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

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

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

1019 ''' 

1020 

1021 file_path = String.T(optional=True) 

1022 file_format = String.T(optional=True) 

1023 file_mtime = Timestamp.T(optional=True) 

1024 file_size = Int.T(optional=True) 

1025 

1026 file_segment = Int.T(optional=True) 

1027 file_element = Int.T(optional=True) 

1028 

1029 kind_id = Int.T() 

1030 codes = Codes.T() 

1031 

1032 tmin_seconds = Int.T(default=0) 

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

1034 tmax_seconds = Int.T(default=0) 

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

1036 

1037 deltat = Float.T(default=0.0) 

1038 

1039 content = Any.T(optional=True) 

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

1041 

1042 content_in_db = False 

1043 

1044 def __init__( 

1045 self, 

1046 file_path=None, 

1047 file_format=None, 

1048 file_mtime=None, 

1049 file_size=None, 

1050 file_segment=None, 

1051 file_element=None, 

1052 kind_id=0, 

1053 codes=CodesX(''), 

1054 tmin_seconds=None, 

1055 tmin_offset=0, 

1056 tmax_seconds=None, 

1057 tmax_offset=0, 

1058 deltat=None, 

1059 content=None, 

1060 raw_content=None, 

1061 tmin=None, 

1062 tmax=None, 

1063 values_nocheck=None): 

1064 

1065 if values_nocheck is not None: 

1066 (self.file_path, self.file_format, self.file_mtime, 

1067 self.file_size, 

1068 self.file_segment, self.file_element, 

1069 self.kind_id, codes_safe_str, 

1070 self.tmin_seconds, self.tmin_offset, 

1071 self.tmax_seconds, self.tmax_offset, 

1072 self.deltat) = values_nocheck 

1073 

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

1075 self.deltat = self.deltat or None 

1076 self.raw_content = {} 

1077 self.content = None 

1078 else: 

1079 if tmin is not None: 

1080 tmin_seconds, tmin_offset = tsplit(tmin) 

1081 

1082 if tmax is not None: 

1083 tmax_seconds, tmax_offset = tsplit(tmax) 

1084 

1085 self.kind_id = int(kind_id) 

1086 self.codes = codes 

1087 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1088 self.tmin_offset = int(tmin_offset) 

1089 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1090 self.tmax_offset = int(tmax_offset) 

1091 self.deltat = float_or_none(deltat) 

1092 self.file_path = str_or_none(file_path) 

1093 self.file_segment = int_or_none(file_segment) 

1094 self.file_element = int_or_none(file_element) 

1095 self.file_format = str_or_none(file_format) 

1096 self.file_mtime = float_or_none(file_mtime) 

1097 self.file_size = int_or_none(file_size) 

1098 self.content = content 

1099 if raw_content is None: 

1100 self.raw_content = {} 

1101 else: 

1102 self.raw_content = raw_content 

1103 

1104 Object.__init__(self, init_props=False) 

1105 

1106 def __eq__(self, other): 

1107 return (isinstance(other, Nut) and 

1108 self.equality_values == other.equality_values) 

1109 

1110 def hash(self): 

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

1112 

1113 def __ne__(self, other): 

1114 return not (self == other) 

1115 

1116 def get_io_backend(self): 

1117 from . import io 

1118 return io.get_backend(self.file_format) 

1119 

1120 def file_modified(self): 

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

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

1123 

1124 @property 

1125 def dkey(self): 

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

1127 

1128 @property 

1129 def key(self): 

1130 return ( 

1131 self.file_path, 

1132 self.file_segment, 

1133 self.file_element, 

1134 self.file_mtime) 

1135 

1136 @property 

1137 def equality_values(self): 

1138 return ( 

1139 self.file_segment, self.file_element, 

1140 self.kind_id, self.codes, 

1141 self.tmin_seconds, self.tmin_offset, 

1142 self.tmax_seconds, self.tmax_offset, self.deltat) 

1143 

1144 def diff(self, other): 

1145 names = [ 

1146 'file_segment', 'file_element', 'kind_id', 'codes', 

1147 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset', 

1148 'deltat'] 

1149 

1150 d = [] 

1151 for name, a, b in zip( 

1152 names, self.equality_values, other.equality_values): 

1153 

1154 if a != b: 

1155 d.append((name, a, b)) 

1156 

1157 return d 

1158 

1159 @property 

1160 def tmin(self): 

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

1162 

1163 @tmin.setter 

1164 def tmin(self, t): 

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

1166 

1167 @property 

1168 def tmax(self): 

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

1170 

1171 @tmax.setter 

1172 def tmax(self, t): 

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

1174 

1175 @property 

1176 def kscale(self): 

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

1178 return 0 

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

1180 

1181 @property 

1182 def waveform_kwargs(self): 

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

1184 

1185 return dict( 

1186 network=network, 

1187 station=station, 

1188 location=location, 

1189 channel=channel, 

1190 extra=extra, 

1191 tmin=self.tmin, 

1192 tmax=self.tmax, 

1193 deltat=self.deltat) 

1194 

1195 @property 

1196 def waveform_promise_kwargs(self): 

1197 return dict( 

1198 codes=self.codes, 

1199 tmin=self.tmin, 

1200 tmax=self.tmax, 

1201 deltat=self.deltat) 

1202 

1203 @property 

1204 def station_kwargs(self): 

1205 network, station, location = self.codes 

1206 return dict( 

1207 codes=self.codes, 

1208 tmin=tmin_or_none(self.tmin), 

1209 tmax=tmax_or_none(self.tmax)) 

1210 

1211 @property 

1212 def channel_kwargs(self): 

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

1214 return dict( 

1215 codes=self.codes, 

1216 tmin=tmin_or_none(self.tmin), 

1217 tmax=tmax_or_none(self.tmax), 

1218 deltat=self.deltat) 

1219 

1220 @property 

1221 def response_kwargs(self): 

1222 return dict( 

1223 codes=self.codes, 

1224 tmin=tmin_or_none(self.tmin), 

1225 tmax=tmax_or_none(self.tmax), 

1226 deltat=self.deltat) 

1227 

1228 @property 

1229 def event_kwargs(self): 

1230 return dict( 

1231 name=self.codes, 

1232 time=self.tmin, 

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

1234 

1235 @property 

1236 def trace_kwargs(self): 

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

1238 

1239 return dict( 

1240 network=network, 

1241 station=station, 

1242 location=location, 

1243 channel=channel, 

1244 extra=extra, 

1245 tmin=self.tmin, 

1246 tmax=self.tmax-self.deltat, 

1247 deltat=self.deltat) 

1248 

1249 @property 

1250 def dummy_trace(self): 

1251 return DummyTrace(self) 

1252 

1253 @property 

1254 def summary(self): 

1255 if self.tmin == self.tmax: 

1256 ts = util.time_to_str(self.tmin) 

1257 else: 

1258 ts = '%s - %s' % ( 

1259 util.time_to_str(self.tmin), 

1260 util.time_to_str(self.tmax)) 

1261 

1262 return ' '.join(( 

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

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

1265 ts)) 

1266 

1267 

1268def make_waveform_nut(**kwargs): 

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

1270 

1271 

1272def make_waveform_promise_nut(**kwargs): 

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

1274 

1275 

1276def make_station_nut(**kwargs): 

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

1278 

1279 

1280def make_channel_nut(**kwargs): 

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

1282 

1283 

1284def make_response_nut(**kwargs): 

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

1286 

1287 

1288def make_event_nut(**kwargs): 

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

1290 

1291 

1292def group_channels(nuts): 

1293 by_ansl = {} 

1294 for nut in nuts: 

1295 if nut.kind_id != CHANNEL: 

1296 continue 

1297 

1298 ansl = nut.codes[:4] 

1299 

1300 if ansl not in by_ansl: 

1301 by_ansl[ansl] = {} 

1302 

1303 group = by_ansl[ansl] 

1304 

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

1306 

1307 if k not in group: 

1308 group[k] = set() 

1309 

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

1311 

1312 return by_ansl 

1313 

1314 

1315class DummyTrace(object): 

1316 

1317 def __init__(self, nut): 

1318 self.nut = nut 

1319 self.codes = nut.codes 

1320 self.meta = {} 

1321 

1322 @property 

1323 def tmin(self): 

1324 return self.nut.tmin 

1325 

1326 @property 

1327 def tmax(self): 

1328 return self.nut.tmax 

1329 

1330 @property 

1331 def deltat(self): 

1332 return self.nut.deltat 

1333 

1334 @property 

1335 def nslc_id(self): 

1336 return self.codes.nslc 

1337 

1338 @property 

1339 def network(self): 

1340 return self.codes.network 

1341 

1342 @property 

1343 def station(self): 

1344 return self.codes.station 

1345 

1346 @property 

1347 def location(self): 

1348 return self.codes.location 

1349 

1350 @property 

1351 def channel(self): 

1352 return self.codes.channel 

1353 

1354 @property 

1355 def extra(self): 

1356 return self.codes.extra 

1357 

1358 def overlaps(self, tmin, tmax): 

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

1360 

1361 

1362def duration_to_str(t): 

1363 if t > 24*3600: 

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

1365 elif t > 3600: 

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

1367 elif t > 60: 

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

1369 else: 

1370 return '%gs' % t 

1371 

1372 

1373class Coverage(Object): 

1374 ''' 

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

1376 ''' 

1377 kind_id = Int.T( 

1378 help='Content type.') 

1379 pattern = Codes.T( 

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

1381 'match.') 

1382 codes = Codes.T( 

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

1384 deltat = Float.T( 

1385 help='Sampling interval [s]', 

1386 optional=True) 

1387 tmin = Timestamp.T( 

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

1389 optional=True) 

1390 tmax = Timestamp.T( 

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

1392 optional=True) 

1393 changes = List.T( 

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

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

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

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

1398 'value duplicate or redundant data coverage.') 

1399 

1400 @classmethod 

1401 def from_values(cls, args): 

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

1403 return cls( 

1404 kind_id=kind_id, 

1405 pattern=pattern, 

1406 codes=codes, 

1407 deltat=deltat, 

1408 tmin=tmin, 

1409 tmax=tmax, 

1410 changes=changes) 

1411 

1412 @property 

1413 def summary(self): 

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

1415 util.time_to_str(self.tmin), 

1416 util.time_to_str(self.tmax)) 

1417 

1418 srate = self.sample_rate 

1419 

1420 return ' '.join(( 

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

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

1423 ts, 

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

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

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

1427 

1428 @property 

1429 def sample_rate(self): 

1430 if self.deltat is None: 

1431 return None 

1432 elif self.deltat == 0.0: 

1433 return 0.0 

1434 else: 

1435 return 1.0 / self.deltat 

1436 

1437 @property 

1438 def labels(self): 

1439 srate = self.sample_rate 

1440 return ( 

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

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

1443 

1444 @property 

1445 def total(self): 

1446 total_t = None 

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

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

1449 

1450 return total_t 

1451 

1452 def iter_spans(self): 

1453 last = None 

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

1455 if last is not None: 

1456 last_t, last_count = last 

1457 if last_count > 0: 

1458 yield last_t, t, last_count 

1459 

1460 last = (t, count) 

1461 

1462 

1463__all__ = [ 

1464 'to_codes', 

1465 'to_codes_guess', 

1466 'to_codes_simple', 

1467 'to_kind', 

1468 'to_kinds', 

1469 'to_kind_id', 

1470 'to_kind_ids', 

1471 'CodesError', 

1472 'Codes', 

1473 'CodesNSLCE', 

1474 'CodesNSL', 

1475 'CodesX', 

1476 'Station', 

1477 'Channel', 

1478 'Sensor', 

1479 'Response', 

1480 'Nut', 

1481 'Coverage', 

1482 'WaveformPromise', 

1483]