Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/squirrel/model.py: 69%

913 statements  

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

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 

19import re 

20import fnmatch 

21import math 

22import hashlib 

23from os import urandom 

24from base64 import urlsafe_b64encode 

25from collections import defaultdict, namedtuple 

26 

27import numpy as num 

28 

29from pyrocko import util 

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

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

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 

39d2r = num.pi / 180. 

40r2d = 1.0 / d2r 

41 

42 

43guts_prefix = 'squirrel' 

44 

45 

46def mkvec(x, y, z): 

47 return num.array([x, y, z], dtype=float) 

48 

49 

50def are_orthogonal(vecs, eps=0.01): 

51 return all(abs( 

52 num.dot(vecs[i], vecs[j]) < eps 

53 for (i, j) in [(0, 1), (1, 2), (2, 0)])) 

54 

55 

56g_codes_pool = {} 

57 

58 

59class CodesError(SquirrelError): 

60 pass 

61 

62 

63class Codes(SObject): 

64 pass 

65 

66 

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

68 if args and kwargs: 

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

70 

71 if len(args) == 1: 

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

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

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

75 args = args[0] 

76 else: 

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

78 

79 nargs = len(args) 

80 if nargs == 5: 

81 t = args 

82 

83 elif nargs == 4: 

84 t = args + ('',) 

85 

86 elif nargs == 0: 

87 d = dict( 

88 network='', 

89 station='', 

90 location='', 

91 channel='', 

92 extra='') 

93 

94 d.update(kwargs) 

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

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

97 

98 else: 

99 raise CodesError( 

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

101 

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

103 raise CodesError( 

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

105 

106 return t 

107 

108 

109CodesNSLCEBase = namedtuple( 

110 'CodesNSLCEBase', [ 

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

112 

113 

114class CodesNSLCE(CodesNSLCEBase, Codes): 

115 ''' 

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

117 

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

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

120 processing applied to a channel. 

121 ''' 

122 

123 __slots__ = () 

124 __hash__ = CodesNSLCEBase.__hash__ 

125 

126 as_dict = CodesNSLCEBase._asdict 

127 

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

129 nargs = len(args) 

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

131 return args[0] 

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

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

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

135 t = ('*', '*', '*', '*', '*') 

136 elif safe_str is not None: 

137 t = safe_str.split('.') 

138 else: 

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

140 

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

142 return g_codes_pool.setdefault(x, x) 

143 

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

145 Codes.__init__(self) 

146 

147 def __str__(self): 

148 if self.extra == '': 

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

150 else: 

151 return '.'.join(self) 

152 

153 def __eq__(self, other): 

154 if not isinstance(other, CodesNSLCE): 

155 other = CodesNSLCE(other) 

156 

157 return CodesNSLCEBase.__eq__(self, other) 

158 

159 def matches(self, pattern): 

160 if not isinstance(pattern, CodesNSLCE): 

161 pattern = CodesNSLCE(pattern) 

162 

163 return match_codes(pattern, self) 

164 

165 @property 

166 def safe_str(self): 

167 return '.'.join(self) 

168 

169 @property 

170 def nslce(self): 

171 return self[:5] 

172 

173 @property 

174 def nslc(self): 

175 return self[:4] 

176 

177 @property 

178 def nsl(self): 

179 return self[:3] 

180 

181 @property 

182 def ns(self): 

183 return self[:2] 

184 

185 @property 

186 def codes_nsl(self): 

187 return CodesNSL(self) 

188 

189 @property 

190 def codes_nsl_star(self): 

191 return CodesNSL(self.network, self.station, '*') 

192 

193 def as_tuple(self): 

194 return tuple(self) 

195 

196 def replace(self, **kwargs): 

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

198 return g_codes_pool.setdefault(x, x) 

199 

200 

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

202 if args and kwargs: 

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

204 

205 if len(args) == 1: 

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

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

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

209 args = args[0] 

210 else: 

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

212 

213 nargs = len(args) 

214 if nargs == 3: 

215 t = args 

216 

217 elif nargs == 0: 

218 d = dict( 

219 network='', 

220 station='', 

221 location='') 

222 

223 d.update(kwargs) 

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

225 'network', 'station', 'location')) 

226 

227 else: 

228 raise CodesError( 

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

230 

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

232 raise CodesError( 

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

234 

235 return t 

236 

237 

238CodesNSLBase = namedtuple( 

239 'CodesNSLBase', [ 

240 'network', 'station', 'location']) 

241 

242 

243class CodesNSL(CodesNSLBase, Codes): 

244 ''' 

245 Codes denominating a seismic station (NSL). 

246 

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

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

249 compatible behaviour in most cases. 

250 ''' 

251 

252 __slots__ = () 

253 __hash__ = CodesNSLBase.__hash__ 

254 

255 as_dict = CodesNSLBase._asdict 

256 

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

258 nargs = len(args) 

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

260 return args[0] 

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

262 t = args[0].nsl 

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

264 t = ('*', '*', '*') 

265 elif safe_str is not None: 

266 t = safe_str.split('.') 

267 else: 

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

269 

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

271 return g_codes_pool.setdefault(x, x) 

272 

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

274 Codes.__init__(self) 

275 

276 def __str__(self): 

277 return '.'.join(self) 

278 

279 def __eq__(self, other): 

280 if not isinstance(other, CodesNSL): 

281 other = CodesNSL(other) 

282 

283 return CodesNSLBase.__eq__(self, other) 

284 

285 def matches(self, pattern): 

286 if not isinstance(pattern, CodesNSL): 

287 pattern = CodesNSL(pattern) 

288 

289 return match_codes(pattern, self) 

290 

291 @property 

292 def safe_str(self): 

293 return '.'.join(self) 

294 

295 @property 

296 def ns(self): 

297 return self[:2] 

298 

299 @property 

300 def nsl(self): 

301 return self[:3] 

302 

303 def as_tuple(self): 

304 return tuple(self) 

305 

306 def replace(self, **kwargs): 

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

308 return g_codes_pool.setdefault(x, x) 

309 

310 

311CodesXBase = namedtuple( 

312 'CodesXBase', [ 

313 'name']) 

314 

315 

316class CodesX(CodesXBase, Codes): 

317 ''' 

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

319 ''' 

320 

321 __slots__ = () 

322 __hash__ = CodesXBase.__hash__ 

323 __eq__ = CodesXBase.__eq__ 

324 

325 as_dict = CodesXBase._asdict 

326 

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

328 if isinstance(name, CodesX): 

329 return name 

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

331 name = '*' 

332 elif safe_str is not None: 

333 name = safe_str 

334 else: 

335 if '.' in name: 

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

337 

338 x = CodesXBase.__new__(cls, name) 

339 return g_codes_pool.setdefault(x, x) 

340 

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

342 Codes.__init__(self) 

343 

344 def __str__(self): 

345 return '.'.join(self) 

346 

347 @property 

348 def safe_str(self): 

349 return '.'.join(self) 

350 

351 @property 

352 def ns(self): 

353 return self[:2] 

354 

355 def as_tuple(self): 

356 return tuple(self) 

357 

358 def replace(self, **kwargs): 

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

360 return g_codes_pool.setdefault(x, x) 

361 

362 

363g_codes_patterns = {} 

364 

365 

366def match_codes(pattern, codes): 

367 spattern = pattern.safe_str 

368 scodes = codes.safe_str 

369 if spattern not in g_codes_patterns: 

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

371 g_codes_patterns[spattern] = rpattern 

372 

373 rpattern = g_codes_patterns[spattern] 

374 return bool(rpattern.match(scodes)) 

375 

376 

377g_content_kinds = [ 

378 'undefined', 

379 'waveform', 

380 'station', 

381 'channel', 

382 'response', 

383 'event', 

384 'waveform_promise', 

385 'empty'] 

386 

387 

388g_codes_classes = [ 

389 CodesX, 

390 CodesNSLCE, 

391 CodesNSL, 

392 CodesNSLCE, 

393 CodesNSLCE, 

394 CodesX, 

395 CodesNSLCE, 

396 CodesX] 

397 

398g_codes_classes_ndot = { 

399 0: CodesX, 

400 2: CodesNSL, 

401 3: CodesNSLCE, 

402 4: CodesNSLCE} 

403 

404 

405def to_codes_simple(kind_id, codes_safe_str): 

406 return g_codes_classes[kind_id](safe_str=codes_safe_str) 

407 

408 

409def to_codes(kind_id, obj): 

410 return g_codes_classes[kind_id](obj) 

411 

412 

413def to_codes_guess(s): 

414 try: 

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

416 except KeyError: 

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

418 

419 

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

421class codes_patterns_list(list): 

422 pass 

423 

424 

425def codes_patterns_for_kind(kind_id, codes): 

426 if isinstance(codes, codes_patterns_list): 

427 return codes 

428 

429 if isinstance(codes, list): 

430 lcodes = codes_patterns_list() 

431 for sc in codes: 

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

433 

434 return lcodes 

435 

436 codes = to_codes(kind_id, codes) 

437 

438 lcodes = codes_patterns_list() 

439 lcodes.append(codes) 

440 if kind_id == STATION: 

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

442 

443 return lcodes 

444 

445 

446g_content_kind_ids = ( 

447 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT, 

448 WAVEFORM_PROMISE, EMPTY) = range(len(g_content_kinds)) 

449 

450 

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

452 

453 

454try: 

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

456except Exception: 

457 g_tmin_queries = g_tmin 

458 

459 

460def to_kind(kind_id): 

461 return g_content_kinds[kind_id] 

462 

463 

464def to_kinds(kind_ids): 

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

466 

467 

468def to_kind_id(kind): 

469 return g_content_kinds.index(kind) 

470 

471 

472def to_kind_ids(kinds): 

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

474 

475 

476g_kind_mask_all = 0xff 

477 

478 

479def to_kind_mask(kinds): 

480 if kinds: 

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

482 else: 

483 return g_kind_mask_all 

484 

485 

486def str_or_none(x): 

487 if x is None: 

488 return None 

489 else: 

490 return str(x) 

491 

492 

493def float_or_none(x): 

494 if x is None: 

495 return None 

496 else: 

497 return float(x) 

498 

499 

500def int_or_none(x): 

501 if x is None: 

502 return None 

503 else: 

504 return int(x) 

505 

506 

507def int_or_g_tmin(x): 

508 if x is None: 

509 return g_tmin 

510 else: 

511 return int(x) 

512 

513 

514def int_or_g_tmax(x): 

515 if x is None: 

516 return g_tmax 

517 else: 

518 return int(x) 

519 

520 

521def tmin_or_none(x): 

522 if x == g_tmin: 

523 return None 

524 else: 

525 return x 

526 

527 

528def tmax_or_none(x): 

529 if x == g_tmax: 

530 return None 

531 else: 

532 return x 

533 

534 

535def time_or_none_to_str(x): 

536 if x is None: 

537 return '...'.ljust(17) 

538 else: 

539 return util.time_to_str(x) 

540 

541 

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

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

544 

545 if len(codes) > 20: 

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

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

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

549 else: 

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

551 if codes else '<none>' 

552 

553 return scodes 

554 

555 

556g_offset_time_unit_inv = 1000000000 

557g_offset_time_unit = 1.0 / g_offset_time_unit_inv 

558 

559 

560def tsplit(t): 

561 if t is None: 

562 return None, 0 

563 

564 t = util.to_time_float(t) 

565 if type(t) is float: 

566 t = round(t, 5) 

567 else: 

568 t = round(t, 9) 

569 

570 seconds = num.floor(t) 

571 offset = t - seconds 

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

573 

574 

575def tjoin(seconds, offset): 

576 if seconds is None: 

577 return None 

578 

579 return util.to_time_float(seconds) \ 

580 + util.to_time_float(offset*g_offset_time_unit) 

581 

582 

583tscale_min = 1 

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

585tscale_logbase = 20 

586 

587tscale_edges = [tscale_min] 

588while True: 

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

590 if tscale_edges[-1] >= tscale_max: 

591 break 

592 

593 

594tscale_edges = num.array(tscale_edges) 

595 

596 

597def tscale_to_kscale(tscale): 

598 

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

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

601 # ... 

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

603 

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

605 

606 

607@squirrel_content 

608class Station(Location): 

609 ''' 

610 A seismic station. 

611 ''' 

612 

613 codes = CodesNSL.T() 

614 

615 tmin = Timestamp.T(optional=True) 

616 tmax = Timestamp.T(optional=True) 

617 

618 description = Unicode.T(optional=True) 

619 

620 def __init__(self, **kwargs): 

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

622 Location.__init__(self, **kwargs) 

623 

624 @property 

625 def time_span(self): 

626 return (self.tmin, self.tmax) 

627 

628 def get_pyrocko_station(self): 

629 ''' 

630 Get station as a classic Pyrocko station object. 

631 

632 :returns: 

633 Converted station object. 

634 :rtype: 

635 :py:class:`pyrocko.model.station.Station` 

636 ''' 

637 from pyrocko import model 

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

639 

640 def _get_pyrocko_station_args(self): 

641 return ( 

642 self.codes.network, 

643 self.codes.station, 

644 self.codes.location, 

645 self.lat, 

646 self.lon, 

647 self.elevation, 

648 self.depth, 

649 self.north_shift, 

650 self.east_shift) 

651 

652 

653@squirrel_content 

654class ChannelBase(Location): 

655 codes = CodesNSLCE.T() 

656 

657 tmin = Timestamp.T(optional=True) 

658 tmax = Timestamp.T(optional=True) 

659 

660 deltat = Float.T(optional=True) 

661 

662 def __init__(self, **kwargs): 

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

664 Location.__init__(self, **kwargs) 

665 

666 @property 

667 def time_span(self): 

668 return (self.tmin, self.tmax) 

669 

670 def _get_sensor_codes(self): 

671 return self.codes.replace( 

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

673 

674 def _get_sensor_args(self): 

675 def getattr_rep(k): 

676 if k == 'codes': 

677 return self._get_sensor_codes() 

678 else: 

679 return getattr(self, k) 

680 

681 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames) 

682 

683 def _get_channel_args(self, component): 

684 def getattr_rep(k): 

685 if k == 'codes': 

686 return self.codes.replace( 

687 channel=self.codes.channel[:-1] + component) 

688 else: 

689 return getattr(self, k) 

690 

691 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames) 

692 

693 def _get_pyrocko_station_args(self): 

694 return ( 

695 self.codes.network, 

696 self.codes.station, 

697 self.codes.location, 

698 self.lat, 

699 self.lon, 

700 self.elevation, 

701 self.depth, 

702 self.north_shift, 

703 self.east_shift) 

704 

705 

706class Channel(ChannelBase): 

707 ''' 

708 A channel of a seismic station. 

709 ''' 

710 

711 dip = Float.T(optional=True) 

712 azimuth = Float.T(optional=True) 

713 

714 def get_pyrocko_channel(self): 

715 from pyrocko import model 

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

717 

718 def _get_pyrocko_channel_args(self): 

719 return ( 

720 self.codes.channel, 

721 self.azimuth, 

722 self.dip) 

723 

724 @property 

725 def orientation_enz(self): 

726 if None in (self.azimuth, self.dip): 

727 return None 

728 

729 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r) 

730 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r) 

731 d = math.sin(self.dip*d2r) 

732 return mkvec(e, n, -d) 

733 

734 

735def cut_intervals(channels): 

736 channels = list(channels) 

737 times = set() 

738 for channel in channels: 

739 if channel.tmin is not None: 

740 times.add(channel.tmin) 

741 if channel.tmax is not None: 

742 times.add(channel.tmax) 

743 

744 times = sorted(times) 

745 

746 if any(channel.tmin is None for channel in channels): 

747 times[0:0] = [None] 

748 

749 if any(channel.tmax is None for channel in channels): 

750 times.append(None) 

751 

752 if len(times) <= 2: 

753 return channels 

754 

755 channels_out = [] 

756 for channel in channels: 

757 for i in range(len(times)-1): 

758 tmin = times[i] 

759 tmax = times[i+1] 

760 if ((channel.tmin is None or ( 

761 tmin is not None and channel.tmin <= tmin)) 

762 and (channel.tmax is None or ( 

763 tmax is not None and tmax <= channel.tmax))): 

764 

765 channel_out = clone(channel) 

766 channel_out.tmin = tmin 

767 channel_out.tmax = tmax 

768 channels_out.append(channel_out) 

769 

770 return channels_out 

771 

772 

773class Sensor(ChannelBase): 

774 ''' 

775 Representation of a channel group. 

776 ''' 

777 

778 channels = List.T(Channel.T()) 

779 

780 @classmethod 

781 def from_channels(cls, channels): 

782 groups = defaultdict(list) 

783 for channel in channels: 

784 groups[channel._get_sensor_codes()].append(channel) 

785 

786 channels_cut = [] 

787 for group in groups.values(): 

788 channels_cut.extend(cut_intervals(group)) 

789 

790 groups = defaultdict(list) 

791 for channel in channels_cut: 

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

793 

794 return [ 

795 cls(channels=channels, 

796 **dict(zip(ChannelBase.T.propnames, args))) 

797 for args, channels in groups.items()] 

798 

799 def channel_vectors(self): 

800 return num.vstack( 

801 [channel.orientation_enz for channel in self.channels]) 

802 

803 def projected_channels(self, system, component_names): 

804 return [ 

805 Channel( 

806 azimuth=math.atan2(e, n) * r2d, 

807 dip=-math.asin(z) * r2d, 

808 **dict(zip( 

809 ChannelBase.T.propnames, 

810 self._get_channel_args(comp)))) 

811 for comp, (e, n, z) in zip(component_names, system)] 

812 

813 def matrix_to(self, system, epsilon=1e-7): 

814 m = num.dot(system, self.channel_vectors().T) 

815 m[num.abs(m) < epsilon] = 0.0 

816 return m 

817 

818 def projection_to(self, system, component_names): 

819 return ( 

820 self.matrix_to(system), 

821 self.channels, 

822 self.projected_channels(system, component_names)) 

823 

824 def projection_to_enz(self): 

825 return self.projection_to(num.identity(3), 'ENZ') 

826 

827 def projection_to_trz(self, source, azimuth=None): 

828 if azimuth is not None: 

829 assert source is None 

830 else: 

831 azimuth = source.azibazi_to(self)[1] + 180. 

832 

833 return self.projection_to(num.array([ 

834 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.], 

835 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.], 

836 [0., 0., 1.]], dtype=float), 'TRZ') 

837 

838 def project_to_enz(self, traces): 

839 from pyrocko import trace 

840 

841 matrix, in_channels, out_channels = self.projection_to_enz() 

842 

843 return trace.project(traces, matrix, in_channels, out_channels) 

844 

845 def project_to_trz(self, source, traces, azimuth=None): 

846 from pyrocko import trace 

847 

848 matrix, in_channels, out_channels = self.projection_to_trz( 

849 source, azimuth=azimuth) 

850 

851 return trace.project(traces, matrix, in_channels, out_channels) 

852 

853 

854observational_quantities = [ 

855 'acceleration', 'velocity', 'displacement', 'pressure', 

856 'rotation_displacement', 'rotation_velocity', 'rotation_acceleration', 

857 'temperature'] 

858 

859 

860technical_quantities = [ 

861 'voltage', 'counts'] 

862 

863 

864class QuantityType(StringChoice): 

865 ''' 

866 Choice of observational or technical quantity. 

867 

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

869 ''' 

870 choices = observational_quantities + technical_quantities 

871 

872 

873class ResponseStage(Object): 

874 ''' 

875 Representation of a response stage. 

876 

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

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

879 downsampling. 

880 ''' 

881 input_quantity = QuantityType.T(optional=True) 

882 input_sample_rate = Float.T(optional=True) 

883 output_quantity = QuantityType.T(optional=True) 

884 output_sample_rate = Float.T(optional=True) 

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

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

887 

888 @property 

889 def stage_type(self): 

890 if self.input_quantity in observational_quantities \ 

891 and self.output_quantity in observational_quantities: 

892 return 'conversion' 

893 

894 if self.input_quantity in observational_quantities \ 

895 and self.output_quantity == 'voltage': 

896 return 'sensor' 

897 

898 elif self.input_quantity == 'voltage' \ 

899 and self.output_quantity == 'voltage': 

900 return 'electronics' 

901 

902 elif self.input_quantity == 'voltage' \ 

903 and self.output_quantity == 'counts': 

904 return 'digitizer' 

905 

906 elif self.decimation_factor is not None \ 

907 and (self.input_quantity is None or self.input_quantity == 'counts') \ 

908 and (self.output_quantity is None or self.output_quantity == 'counts'): # noqa 

909 return 'decimation' 

910 

911 elif self.input_quantity in observational_quantities \ 

912 and self.output_quantity == 'counts': 

913 return 'combination' 

914 

915 else: 

916 return 'unknown' 

917 

918 @property 

919 def decimation_factor(self): 

920 irate = self.input_sample_rate 

921 orate = self.output_sample_rate 

922 if irate is not None and orate is not None \ 

923 and irate > orate and irate / orate > 1.0: 

924 

925 return irate / orate 

926 else: 

927 return None 

928 

929 @property 

930 def summary_quantities(self): 

931 return '%s -> %s' % ( 

932 self.input_quantity or '?', 

933 self.output_quantity or '?') 

934 

935 @property 

936 def summary_rates(self): 

937 irate = self.input_sample_rate 

938 orate = self.output_sample_rate 

939 factor = self.decimation_factor 

940 

941 if irate and orate is None: 

942 return '%g Hz' % irate 

943 

944 elif orate and irate is None: 

945 return '%g Hz' % orate 

946 

947 elif irate and orate and irate == orate: 

948 return '%g Hz' % irate 

949 

950 elif any(x for x in (irate, orate, factor)): 

951 return '%s -> %s Hz (%s)' % ( 

952 '%g' % irate if irate else '?', 

953 '%g' % orate if orate else '?', 

954 ':%g' % factor if factor else '?') 

955 else: 

956 return '' 

957 

958 @property 

959 def summary_elements(self): 

960 xs = [x.summary for x in self.elements] 

961 return '%s' % ('*'.join(x for x in xs if x != 'one') or 'one') 

962 

963 @property 

964 def summary_log(self): 

965 return ''.join(sorted(set(x[0][0].upper() for x in self.log))) 

966 

967 @property 

968 def summary_entries(self): 

969 return ( 

970 self.__class__.__name__, 

971 self.stage_type, 

972 self.summary_log, 

973 self.summary_quantities, 

974 self.summary_rates, 

975 self.summary_elements) 

976 

977 @property 

978 def summary(self): 

979 return util.fmt_summary(self.summary_entries, (10, 15, 3, 30, 30, 0)) 

980 

981 def get_effective(self): 

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

983 

984 

985D = 'displacement' 

986V = 'velocity' 

987A = 'acceleration' 

988 

989g_converters = { 

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

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

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

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

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

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

996 

997 

998def response_converters(input_quantity, output_quantity): 

999 if input_quantity is None or input_quantity == output_quantity: 

1000 return [] 

1001 

1002 if output_quantity is None: 

1003 raise ConversionError('Unspecified target quantity.') 

1004 

1005 try: 

1006 return [g_converters[input_quantity, output_quantity]] 

1007 

1008 except KeyError: 

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

1010 input_quantity, output_quantity)) 

1011 

1012 

1013@squirrel_content 

1014class Response(Object): 

1015 ''' 

1016 The instrument response of a seismic station channel. 

1017 ''' 

1018 

1019 codes = CodesNSLCE.T() 

1020 tmin = Timestamp.T(optional=True) 

1021 tmax = Timestamp.T(optional=True) 

1022 

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

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

1025 

1026 deltat = Float.T(optional=True) 

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

1028 

1029 def __init__(self, **kwargs): 

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

1031 Object.__init__(self, **kwargs) 

1032 self._effective_responses_cache = {} 

1033 

1034 @property 

1035 def time_span(self): 

1036 return (self.tmin, self.tmax) 

1037 

1038 @property 

1039 def nstages(self): 

1040 return len(self.stages) 

1041 

1042 @property 

1043 def input_quantity(self): 

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

1045 

1046 @property 

1047 def output_quantity(self): 

1048 return self.stages[-1].output_quantity if self.stages else None 

1049 

1050 @property 

1051 def output_sample_rate(self): 

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

1053 

1054 @property 

1055 def summary_stages(self): 

1056 def grouped(xs): 

1057 xs = list(xs) 

1058 g = [] 

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

1060 g.append(xs[i]) 

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

1062 yield g 

1063 g = [] 

1064 

1065 if g: 

1066 yield g 

1067 

1068 return '+'.join( 

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

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

1071 

1072 @property 

1073 def summary_quantities(self): 

1074 orate = self.output_sample_rate 

1075 return '%s -> %s%s' % ( 

1076 self.input_quantity or '?', 

1077 self.output_quantity or '?', 

1078 ' @ %g Hz' % orate if orate else '') 

1079 

1080 @property 

1081 def summary_log(self): 

1082 return ''.join(sorted(set(x[0][0].upper() for x in self.log))) 

1083 

1084 @property 

1085 def summary_entries(self): 

1086 return ( 

1087 self.__class__.__name__, 

1088 str(self.codes), 

1089 self.str_time_span, 

1090 self.summary_log, 

1091 self.summary_quantities, 

1092 self.summary_stages) 

1093 

1094 @property 

1095 def summary(self): 

1096 return util.fmt_summary(self.summary_entries, (10, 20, 55, 3, 35, 0)) 

1097 

1098 def get_effective(self, input_quantity=None, stages=(None, None)): 

1099 cache_key = (input_quantity, stages) 

1100 if cache_key in self._effective_responses_cache: 

1101 return self._effective_responses_cache[cache_key] 

1102 

1103 try: 

1104 elements = response_converters(input_quantity, self.input_quantity) 

1105 except ConversionError as e: 

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

1107 

1108 elements.extend( 

1109 stage.get_effective() for stage in self.stages[slice(*stages)]) 

1110 

1111 if input_quantity is None \ 

1112 or input_quantity == self.input_quantity: 

1113 checkpoints = self.checkpoints 

1114 else: 

1115 checkpoints = [] 

1116 

1117 resp = MultiplyResponse( 

1118 responses=simplify_responses(elements), 

1119 checkpoints=checkpoints) 

1120 

1121 self._effective_responses_cache[cache_key] = resp 

1122 return resp 

1123 

1124 

1125@squirrel_content 

1126class Event(Object): 

1127 ''' 

1128 A seismic event. 

1129 ''' 

1130 

1131 name = String.T(optional=True) 

1132 time = Timestamp.T() 

1133 duration = Float.T(optional=True) 

1134 

1135 lat = Float.T() 

1136 lon = Float.T() 

1137 depth = Float.T(optional=True) 

1138 

1139 magnitude = Float.T(optional=True) 

1140 

1141 def get_pyrocko_event(self): 

1142 from pyrocko import model 

1143 model.Event( 

1144 name=self.name, 

1145 time=self.time, 

1146 lat=self.lat, 

1147 lon=self.lon, 

1148 depth=self.depth, 

1149 magnitude=self.magnitude, 

1150 duration=self.duration) 

1151 

1152 @property 

1153 def time_span(self): 

1154 return (self.time, self.time) 

1155 

1156 

1157def ehash(s): 

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

1159 

1160 

1161def random_name(n=8): 

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

1163 

1164 

1165@squirrel_content 

1166class WaveformPromise(Object): 

1167 ''' 

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

1169 

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

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

1172 calls to 

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

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

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

1176 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveforms`, and no local 

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

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

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

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

1181 queries for waveforms missing at the remote site. 

1182 ''' 

1183 

1184 codes = CodesNSLCE.T() 

1185 tmin = Timestamp.T() 

1186 tmax = Timestamp.T() 

1187 

1188 deltat = Float.T(optional=True) 

1189 

1190 source_hash = String.T() 

1191 

1192 def __init__(self, **kwargs): 

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

1194 Object.__init__(self, **kwargs) 

1195 

1196 @property 

1197 def time_span(self): 

1198 return (self.tmin, self.tmax) 

1199 

1200 

1201class InvalidWaveform(Exception): 

1202 pass 

1203 

1204 

1205class WaveformOrder(Object): 

1206 ''' 

1207 Waveform request information. 

1208 ''' 

1209 

1210 source_id = String.T() 

1211 codes = CodesNSLCE.T() 

1212 deltat = Float.T() 

1213 tmin = Timestamp.T() 

1214 tmax = Timestamp.T() 

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

1216 time_created = Timestamp.T() 

1217 anxious = Duration.T(default=600.) 

1218 

1219 @property 

1220 def client(self): 

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

1222 

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

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

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

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

1227 

1228 def validate(self, tr): 

1229 if tr.ydata.size == 0: 

1230 raise InvalidWaveform( 

1231 'waveform with zero data samples') 

1232 

1233 if tr.deltat != self.deltat: 

1234 raise InvalidWaveform( 

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

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

1237 tr.deltat, self.deltat)) 

1238 

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

1240 raise InvalidWaveform('waveform has NaN values') 

1241 

1242 def estimate_nsamples(self): 

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

1244 

1245 def is_near_real_time(self): 

1246 return self.tmax > self.time_created - self.anxious 

1247 

1248 

1249def order_summary(orders): 

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

1251 if len(codes_list) > 3: 

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

1253 len(orders), 

1254 util.plural_s(orders), 

1255 str(codes_list[0]), 

1256 str(codes_list[1])) 

1257 

1258 else: 

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

1260 len(orders), 

1261 util.plural_s(orders), 

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

1263 

1264 

1265class Nut(Object): 

1266 ''' 

1267 Index entry referencing an elementary piece of content. 

1268 

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

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

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

1272 ''' 

1273 

1274 file_path = String.T(optional=True) 

1275 file_format = String.T(optional=True) 

1276 file_mtime = Timestamp.T(optional=True) 

1277 file_size = Int.T(optional=True) 

1278 

1279 file_segment = Int.T(optional=True) 

1280 file_element = Int.T(optional=True) 

1281 

1282 kind_id = Int.T() 

1283 codes = Codes.T() 

1284 

1285 tmin_seconds = Int.T(default=0) 

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

1287 tmax_seconds = Int.T(default=0) 

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

1289 

1290 deltat = Float.T(default=0.0) 

1291 

1292 content = Any.T(optional=True) 

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

1294 

1295 content_in_db = False 

1296 

1297 def __init__( 

1298 self, 

1299 file_path=None, 

1300 file_format=None, 

1301 file_mtime=None, 

1302 file_size=None, 

1303 file_segment=None, 

1304 file_element=None, 

1305 kind_id=0, 

1306 codes=CodesX(''), 

1307 tmin_seconds=None, 

1308 tmin_offset=0, 

1309 tmax_seconds=None, 

1310 tmax_offset=0, 

1311 deltat=None, 

1312 content=None, 

1313 raw_content=None, 

1314 tmin=None, 

1315 tmax=None, 

1316 values_nocheck=None): 

1317 

1318 if values_nocheck is not None: 

1319 (self.file_path, self.file_format, self.file_mtime, 

1320 self.file_size, 

1321 self.file_segment, self.file_element, 

1322 self.kind_id, codes_safe_str, 

1323 self.tmin_seconds, self.tmin_offset, 

1324 self.tmax_seconds, self.tmax_offset, 

1325 self.deltat) = values_nocheck 

1326 

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

1328 self.deltat = self.deltat or None 

1329 self.raw_content = {} 

1330 self.content = None 

1331 else: 

1332 if tmin is not None: 

1333 tmin_seconds, tmin_offset = tsplit(tmin) 

1334 

1335 if tmax is not None: 

1336 tmax_seconds, tmax_offset = tsplit(tmax) 

1337 

1338 self.kind_id = int(kind_id) 

1339 self.codes = codes 

1340 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1341 self.tmin_offset = int(tmin_offset) 

1342 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1343 self.tmax_offset = int(tmax_offset) 

1344 self.deltat = float_or_none(deltat) 

1345 self.file_path = str_or_none(file_path) 

1346 self.file_segment = int_or_none(file_segment) 

1347 self.file_element = int_or_none(file_element) 

1348 self.file_format = str_or_none(file_format) 

1349 self.file_mtime = float_or_none(file_mtime) 

1350 self.file_size = int_or_none(file_size) 

1351 self.content = content 

1352 if raw_content is None: 

1353 self.raw_content = {} 

1354 else: 

1355 self.raw_content = raw_content 

1356 

1357 Object.__init__(self, init_props=False) 

1358 

1359 def __eq__(self, other): 

1360 return (isinstance(other, Nut) and 

1361 self.equality_values == other.equality_values) 

1362 

1363 def hash(self): 

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

1365 

1366 def __ne__(self, other): 

1367 return not (self == other) 

1368 

1369 def get_io_backend(self): 

1370 from . import io 

1371 return io.get_backend(self.file_format) 

1372 

1373 def file_modified(self): 

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

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

1376 

1377 @property 

1378 def dkey(self): 

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

1380 

1381 @property 

1382 def key(self): 

1383 return ( 

1384 self.file_path, 

1385 self.file_segment, 

1386 self.file_element, 

1387 self.file_mtime) 

1388 

1389 @property 

1390 def equality_values(self): 

1391 return ( 

1392 self.file_segment, self.file_element, 

1393 self.kind_id, self.codes, 

1394 self.tmin_seconds, self.tmin_offset, 

1395 self.tmax_seconds, self.tmax_offset, self.deltat) 

1396 

1397 def diff(self, other): 

1398 names = [ 

1399 'file_segment', 'file_element', 'kind_id', 'codes', 

1400 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset', 

1401 'deltat'] 

1402 

1403 d = [] 

1404 for name, a, b in zip( 

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

1406 

1407 if a != b: 

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

1409 

1410 return d 

1411 

1412 @property 

1413 def tmin(self): 

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

1415 

1416 @tmin.setter 

1417 def tmin(self, t): 

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

1419 

1420 @property 

1421 def tmax(self): 

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

1423 

1424 @tmax.setter 

1425 def tmax(self, t): 

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

1427 

1428 @property 

1429 def kscale(self): 

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

1431 return 0 

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

1433 

1434 @property 

1435 def waveform_kwargs(self): 

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

1437 

1438 return dict( 

1439 network=network, 

1440 station=station, 

1441 location=location, 

1442 channel=channel, 

1443 extra=extra, 

1444 tmin=self.tmin, 

1445 tmax=self.tmax, 

1446 deltat=self.deltat) 

1447 

1448 @property 

1449 def waveform_promise_kwargs(self): 

1450 return dict( 

1451 codes=self.codes, 

1452 tmin=self.tmin, 

1453 tmax=self.tmax, 

1454 deltat=self.deltat) 

1455 

1456 @property 

1457 def station_kwargs(self): 

1458 network, station, location = self.codes 

1459 return dict( 

1460 codes=self.codes, 

1461 tmin=tmin_or_none(self.tmin), 

1462 tmax=tmax_or_none(self.tmax)) 

1463 

1464 @property 

1465 def channel_kwargs(self): 

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

1467 return dict( 

1468 codes=self.codes, 

1469 tmin=tmin_or_none(self.tmin), 

1470 tmax=tmax_or_none(self.tmax), 

1471 deltat=self.deltat) 

1472 

1473 @property 

1474 def response_kwargs(self): 

1475 return dict( 

1476 codes=self.codes, 

1477 tmin=tmin_or_none(self.tmin), 

1478 tmax=tmax_or_none(self.tmax), 

1479 deltat=self.deltat) 

1480 

1481 @property 

1482 def event_kwargs(self): 

1483 return dict( 

1484 name=self.codes, 

1485 time=self.tmin, 

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

1487 

1488 @property 

1489 def trace_kwargs(self): 

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

1491 

1492 return dict( 

1493 network=network, 

1494 station=station, 

1495 location=location, 

1496 channel=channel, 

1497 extra=extra, 

1498 tmin=self.tmin, 

1499 tmax=self.tmax-self.deltat, 

1500 deltat=self.deltat) 

1501 

1502 @property 

1503 def dummy_trace(self): 

1504 return DummyTrace(self) 

1505 

1506 @property 

1507 def summary_entries(self): 

1508 if self.tmin == self.tmax: 

1509 ts = util.time_to_str(self.tmin) 

1510 else: 

1511 ts = '%s - %s' % ( 

1512 util.time_to_str(self.tmin), 

1513 util.time_to_str(self.tmax)) 

1514 

1515 return ( 

1516 self.__class__.__name__, 

1517 to_kind(self.kind_id), 

1518 str(self.codes), 

1519 ts) 

1520 

1521 @property 

1522 def summary(self): 

1523 return util.fmt_summary(self.summary_entries, (10, 16, 20, 0)) 

1524 

1525 

1526def make_waveform_nut(**kwargs): 

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

1528 

1529 

1530def make_waveform_promise_nut(**kwargs): 

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

1532 

1533 

1534def make_station_nut(**kwargs): 

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

1536 

1537 

1538def make_channel_nut(**kwargs): 

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

1540 

1541 

1542def make_response_nut(**kwargs): 

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

1544 

1545 

1546def make_event_nut(**kwargs): 

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

1548 

1549 

1550def group_channels(nuts): 

1551 by_ansl = {} 

1552 for nut in nuts: 

1553 if nut.kind_id != CHANNEL: 

1554 continue 

1555 

1556 ansl = nut.codes[:4] 

1557 

1558 if ansl not in by_ansl: 

1559 by_ansl[ansl] = {} 

1560 

1561 group = by_ansl[ansl] 

1562 

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

1564 

1565 if k not in group: 

1566 group[k] = set() 

1567 

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

1569 

1570 return by_ansl 

1571 

1572 

1573class DummyTrace(object): 

1574 

1575 def __init__(self, nut): 

1576 self.nut = nut 

1577 self.codes = nut.codes 

1578 self.meta = {} 

1579 

1580 @property 

1581 def tmin(self): 

1582 return self.nut.tmin 

1583 

1584 @property 

1585 def tmax(self): 

1586 return self.nut.tmax 

1587 

1588 @property 

1589 def deltat(self): 

1590 return self.nut.deltat 

1591 

1592 @property 

1593 def nslc_id(self): 

1594 return self.codes.nslc 

1595 

1596 @property 

1597 def network(self): 

1598 return self.codes.network 

1599 

1600 @property 

1601 def station(self): 

1602 return self.codes.station 

1603 

1604 @property 

1605 def location(self): 

1606 return self.codes.location 

1607 

1608 @property 

1609 def channel(self): 

1610 return self.codes.channel 

1611 

1612 @property 

1613 def extra(self): 

1614 return self.codes.extra 

1615 

1616 def overlaps(self, tmin, tmax): 

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

1618 

1619 

1620def duration_to_str(t): 

1621 if t > 24*3600: 

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

1623 elif t > 3600: 

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

1625 elif t > 60: 

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

1627 else: 

1628 return '%gs' % t 

1629 

1630 

1631class Coverage(Object): 

1632 ''' 

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

1634 ''' 

1635 kind_id = Int.T( 

1636 help='Content type.') 

1637 pattern = Codes.T( 

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

1639 'match.') 

1640 codes = Codes.T( 

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

1642 deltat = Float.T( 

1643 help='Sampling interval [s]', 

1644 optional=True) 

1645 tmin = Timestamp.T( 

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

1647 optional=True) 

1648 tmax = Timestamp.T( 

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

1650 optional=True) 

1651 changes = List.T( 

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

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

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

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

1656 'value duplicate or redundant data coverage.') 

1657 

1658 @classmethod 

1659 def from_values(cls, args): 

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

1661 return cls( 

1662 kind_id=kind_id, 

1663 pattern=pattern, 

1664 codes=codes, 

1665 deltat=deltat, 

1666 tmin=tmin, 

1667 tmax=tmax, 

1668 changes=changes) 

1669 

1670 @property 

1671 def summary_entries(self): 

1672 ts = '%s - %s' % ( 

1673 util.time_to_str(self.tmin), 

1674 util.time_to_str(self.tmax)) 

1675 

1676 srate = self.sample_rate 

1677 

1678 total = self.total 

1679 

1680 return ( 

1681 self.__class__.__name__, 

1682 to_kind(self.kind_id), 

1683 str(self.codes), 

1684 ts, 

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

1686 '%i' % len(self.changes), 

1687 '%s' % duration_to_str(total) if total else 'none') 

1688 

1689 @property 

1690 def summary(self): 

1691 return util.fmt_summary( 

1692 self.summary_entries, 

1693 (10, 16, 20, 55, 10, 4, 0)) 

1694 

1695 @property 

1696 def sample_rate(self): 

1697 if self.deltat is None: 

1698 return None 

1699 elif self.deltat == 0.0: 

1700 return 0.0 

1701 else: 

1702 return 1.0 / self.deltat 

1703 

1704 @property 

1705 def labels(self): 

1706 srate = self.sample_rate 

1707 return ( 

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

1709 '%.4g' % srate if srate else '') 

1710 

1711 @property 

1712 def total(self): 

1713 total_t = None 

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

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

1716 

1717 return total_t 

1718 

1719 def iter_spans(self): 

1720 last = None 

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

1722 if last is not None: 

1723 last_t, last_count = last 

1724 if last_count > 0: 

1725 yield last_t, t, last_count 

1726 

1727 last = (t, count) 

1728 

1729 def iter_uncovered_by(self, other): 

1730 a = self 

1731 b = other 

1732 ia = ib = -1 

1733 ca = cb = 0 

1734 last = None 

1735 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)): 

1736 if ib + 1 == len(b.changes): 

1737 ia += 1 

1738 t, ca = a.changes[ia] 

1739 elif ia + 1 == len(a.changes): 

1740 ib += 1 

1741 t, cb = b.changes[ib] 

1742 elif a.changes[ia+1][0] < b.changes[ib+1][0]: 

1743 ia += 1 

1744 t, ca = a.changes[ia] 

1745 else: 

1746 ib += 1 

1747 t, cb = b.changes[ib] 

1748 

1749 if last is not None: 

1750 tl, cal, cbl = last 

1751 if tl < t and cal > 0 and cbl == 0: 

1752 yield tl, t, ia, ib 

1753 

1754 last = (t, ca, cb) 

1755 

1756 def iter_uncovered_by_combined(self, other): 

1757 ib_last = None 

1758 group = None 

1759 for tmin, tmax, _, ib in self.iter_uncovered_by(other): 

1760 if ib_last is None or ib != ib_last: 

1761 if group: 

1762 yield (group[0][0], group[-1][1]) 

1763 

1764 group = [] 

1765 

1766 group.append((tmin, tmax)) 

1767 ib_last = ib 

1768 

1769 if group: 

1770 yield (group[0][0], group[-1][1]) 

1771 

1772 

1773__all__ = [ 

1774 'to_codes', 

1775 'to_codes_guess', 

1776 'to_codes_simple', 

1777 'to_kind', 

1778 'to_kinds', 

1779 'to_kind_id', 

1780 'to_kind_ids', 

1781 'CodesError', 

1782 'Codes', 

1783 'CodesNSLCE', 

1784 'CodesNSL', 

1785 'CodesX', 

1786 'Station', 

1787 'Channel', 

1788 'Sensor', 

1789 'Response', 

1790 'Nut', 

1791 'Coverage', 

1792 'WaveformPromise', 

1793]