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 2023-09-28 11:43 +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 

386 

387g_codes_classes = [ 

388 CodesX, 

389 CodesNSLCE, 

390 CodesNSL, 

391 CodesNSLCE, 

392 CodesNSLCE, 

393 CodesX, 

394 CodesNSLCE] 

395 

396g_codes_classes_ndot = { 

397 0: CodesX, 

398 2: CodesNSL, 

399 3: CodesNSLCE, 

400 4: CodesNSLCE} 

401 

402 

403def to_codes_simple(kind_id, codes_safe_str): 

404 return g_codes_classes[kind_id](safe_str=codes_safe_str) 

405 

406 

407def to_codes(kind_id, obj): 

408 return g_codes_classes[kind_id](obj) 

409 

410 

411def to_codes_guess(s): 

412 try: 

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

414 except KeyError: 

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

416 

417 

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

419class codes_patterns_list(list): 

420 pass 

421 

422 

423def codes_patterns_for_kind(kind_id, codes): 

424 if isinstance(codes, codes_patterns_list): 

425 return codes 

426 

427 if isinstance(codes, list): 

428 lcodes = codes_patterns_list() 

429 for sc in codes: 

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

431 

432 return lcodes 

433 

434 codes = to_codes(kind_id, codes) 

435 

436 lcodes = codes_patterns_list() 

437 lcodes.append(codes) 

438 if kind_id == STATION: 

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

440 

441 return lcodes 

442 

443 

444g_content_kind_ids = ( 

445 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT, 

446 WAVEFORM_PROMISE) = range(len(g_content_kinds)) 

447 

448 

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

450 

451 

452try: 

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

454except Exception: 

455 g_tmin_queries = g_tmin 

456 

457 

458def to_kind(kind_id): 

459 return g_content_kinds[kind_id] 

460 

461 

462def to_kinds(kind_ids): 

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

464 

465 

466def to_kind_id(kind): 

467 return g_content_kinds.index(kind) 

468 

469 

470def to_kind_ids(kinds): 

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

472 

473 

474g_kind_mask_all = 0xff 

475 

476 

477def to_kind_mask(kinds): 

478 if kinds: 

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

480 else: 

481 return g_kind_mask_all 

482 

483 

484def str_or_none(x): 

485 if x is None: 

486 return None 

487 else: 

488 return str(x) 

489 

490 

491def float_or_none(x): 

492 if x is None: 

493 return None 

494 else: 

495 return float(x) 

496 

497 

498def int_or_none(x): 

499 if x is None: 

500 return None 

501 else: 

502 return int(x) 

503 

504 

505def int_or_g_tmin(x): 

506 if x is None: 

507 return g_tmin 

508 else: 

509 return int(x) 

510 

511 

512def int_or_g_tmax(x): 

513 if x is None: 

514 return g_tmax 

515 else: 

516 return int(x) 

517 

518 

519def tmin_or_none(x): 

520 if x == g_tmin: 

521 return None 

522 else: 

523 return x 

524 

525 

526def tmax_or_none(x): 

527 if x == g_tmax: 

528 return None 

529 else: 

530 return x 

531 

532 

533def time_or_none_to_str(x): 

534 if x is None: 

535 return '...'.ljust(17) 

536 else: 

537 return util.time_to_str(x) 

538 

539 

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

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

542 

543 if len(codes) > 20: 

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

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

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

547 else: 

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

549 if codes else '<none>' 

550 

551 return scodes 

552 

553 

554g_offset_time_unit_inv = 1000000000 

555g_offset_time_unit = 1.0 / g_offset_time_unit_inv 

556 

557 

558def tsplit(t): 

559 if t is None: 

560 return None, 0 

561 

562 t = util.to_time_float(t) 

563 if type(t) is float: 

564 t = round(t, 5) 

565 else: 

566 t = round(t, 9) 

567 

568 seconds = num.floor(t) 

569 offset = t - seconds 

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

571 

572 

573def tjoin(seconds, offset): 

574 if seconds is None: 

575 return None 

576 

577 return util.to_time_float(seconds) \ 

578 + util.to_time_float(offset*g_offset_time_unit) 

579 

580 

581tscale_min = 1 

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

583tscale_logbase = 20 

584 

585tscale_edges = [tscale_min] 

586while True: 

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

588 if tscale_edges[-1] >= tscale_max: 

589 break 

590 

591 

592tscale_edges = num.array(tscale_edges) 

593 

594 

595def tscale_to_kscale(tscale): 

596 

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

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

599 # ... 

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

601 

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

603 

604 

605@squirrel_content 

606class Station(Location): 

607 ''' 

608 A seismic station. 

609 ''' 

610 

611 codes = CodesNSL.T() 

612 

613 tmin = Timestamp.T(optional=True) 

614 tmax = Timestamp.T(optional=True) 

615 

616 description = Unicode.T(optional=True) 

617 

618 def __init__(self, **kwargs): 

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

620 Location.__init__(self, **kwargs) 

621 

622 @property 

623 def time_span(self): 

624 return (self.tmin, self.tmax) 

625 

626 def get_pyrocko_station(self): 

627 ''' 

628 Get station as a classic Pyrocko station object. 

629 

630 :returns: 

631 Converted station object. 

632 :rtype: 

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

634 ''' 

635 from pyrocko import model 

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

637 

638 def _get_pyrocko_station_args(self): 

639 return ( 

640 self.codes.network, 

641 self.codes.station, 

642 self.codes.location, 

643 self.lat, 

644 self.lon, 

645 self.elevation, 

646 self.depth, 

647 self.north_shift, 

648 self.east_shift) 

649 

650 

651@squirrel_content 

652class ChannelBase(Location): 

653 codes = CodesNSLCE.T() 

654 

655 tmin = Timestamp.T(optional=True) 

656 tmax = Timestamp.T(optional=True) 

657 

658 deltat = Float.T(optional=True) 

659 

660 def __init__(self, **kwargs): 

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

662 Location.__init__(self, **kwargs) 

663 

664 @property 

665 def time_span(self): 

666 return (self.tmin, self.tmax) 

667 

668 def _get_sensor_codes(self): 

669 return self.codes.replace( 

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

671 

672 def _get_sensor_args(self): 

673 def getattr_rep(k): 

674 if k == 'codes': 

675 return self._get_sensor_codes() 

676 else: 

677 return getattr(self, k) 

678 

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

680 

681 def _get_channel_args(self, component): 

682 def getattr_rep(k): 

683 if k == 'codes': 

684 return self.codes.replace( 

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

686 else: 

687 return getattr(self, k) 

688 

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

690 

691 def _get_pyrocko_station_args(self): 

692 return ( 

693 self.codes.network, 

694 self.codes.station, 

695 self.codes.location, 

696 self.lat, 

697 self.lon, 

698 self.elevation, 

699 self.depth, 

700 self.north_shift, 

701 self.east_shift) 

702 

703 

704class Channel(ChannelBase): 

705 ''' 

706 A channel of a seismic station. 

707 ''' 

708 

709 dip = Float.T(optional=True) 

710 azimuth = Float.T(optional=True) 

711 

712 def get_pyrocko_channel(self): 

713 from pyrocko import model 

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

715 

716 def _get_pyrocko_channel_args(self): 

717 return ( 

718 self.codes.channel, 

719 self.azimuth, 

720 self.dip) 

721 

722 @property 

723 def orientation_enz(self): 

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

725 return None 

726 

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

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

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

730 return mkvec(e, n, -d) 

731 

732 

733def cut_intervals(channels): 

734 channels = list(channels) 

735 times = set() 

736 for channel in channels: 

737 if channel.tmin is not None: 

738 times.add(channel.tmin) 

739 if channel.tmax is not None: 

740 times.add(channel.tmax) 

741 

742 times = sorted(times) 

743 

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

745 times[0:0] = [None] 

746 

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

748 times.append(None) 

749 

750 if len(times) <= 2: 

751 return channels 

752 

753 channels_out = [] 

754 for channel in channels: 

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

756 tmin = times[i] 

757 tmax = times[i+1] 

758 if ((channel.tmin is None or ( 

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

760 and (channel.tmax is None or ( 

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

762 

763 channel_out = clone(channel) 

764 channel_out.tmin = tmin 

765 channel_out.tmax = tmax 

766 channels_out.append(channel_out) 

767 

768 return channels_out 

769 

770 

771class Sensor(ChannelBase): 

772 ''' 

773 Representation of a channel group. 

774 ''' 

775 

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

777 

778 @classmethod 

779 def from_channels(cls, channels): 

780 groups = defaultdict(list) 

781 for channel in channels: 

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

783 

784 channels_cut = [] 

785 for group in groups.values(): 

786 channels_cut.extend(cut_intervals(group)) 

787 

788 groups = defaultdict(list) 

789 for channel in channels_cut: 

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

791 

792 return [ 

793 cls(channels=channels, 

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

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

796 

797 def channel_vectors(self): 

798 return num.vstack( 

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

800 

801 def projected_channels(self, system, component_names): 

802 return [ 

803 Channel( 

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

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

806 **dict(zip( 

807 ChannelBase.T.propnames, 

808 self._get_channel_args(comp)))) 

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

810 

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

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

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

814 return m 

815 

816 def projection_to(self, system, component_names): 

817 return ( 

818 self.matrix_to(system), 

819 self.channels, 

820 self.projected_channels(system, component_names)) 

821 

822 def projection_to_enz(self): 

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

824 

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

826 if azimuth is not None: 

827 assert source is None 

828 else: 

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

830 

831 return self.projection_to(num.array([ 

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

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

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

835 

836 def project_to_enz(self, traces): 

837 from pyrocko import trace 

838 

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

840 

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

842 

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

844 from pyrocko import trace 

845 

846 matrix, in_channels, out_channels = self.projection_to_trz( 

847 source, azimuth=azimuth) 

848 

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

850 

851 

852observational_quantities = [ 

853 'acceleration', 'velocity', 'displacement', 'pressure', 

854 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration', 

855 'temperature'] 

856 

857 

858technical_quantities = [ 

859 'voltage', 'counts'] 

860 

861 

862class QuantityType(StringChoice): 

863 ''' 

864 Choice of observational or technical quantity. 

865 

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

867 ''' 

868 choices = observational_quantities + technical_quantities 

869 

870 

871class ResponseStage(Object): 

872 ''' 

873 Representation of a response stage. 

874 

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

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

877 downsampling. 

878 ''' 

879 input_quantity = QuantityType.T(optional=True) 

880 input_sample_rate = Float.T(optional=True) 

881 output_quantity = QuantityType.T(optional=True) 

882 output_sample_rate = Float.T(optional=True) 

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

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

885 

886 @property 

887 def stage_type(self): 

888 if self.input_quantity in observational_quantities \ 

889 and self.output_quantity in observational_quantities: 

890 return 'conversion' 

891 

892 if self.input_quantity in observational_quantities \ 

893 and self.output_quantity == 'voltage': 

894 return 'sensor' 

895 

896 elif self.input_quantity == 'voltage' \ 

897 and self.output_quantity == 'voltage': 

898 return 'electronics' 

899 

900 elif self.input_quantity == 'voltage' \ 

901 and self.output_quantity == 'counts': 

902 return 'digitizer' 

903 

904 elif self.decimation_factor is not None \ 

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

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

907 return 'decimation' 

908 

909 elif self.input_quantity in observational_quantities \ 

910 and self.output_quantity == 'counts': 

911 return 'combination' 

912 

913 else: 

914 return 'unknown' 

915 

916 @property 

917 def decimation_factor(self): 

918 irate = self.input_sample_rate 

919 orate = self.output_sample_rate 

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

921 and irate > orate and irate / orate > 1.0: 

922 

923 return irate / orate 

924 else: 

925 return None 

926 

927 @property 

928 def summary_quantities(self): 

929 return '%s -> %s' % ( 

930 self.input_quantity or '?', 

931 self.output_quantity or '?') 

932 

933 @property 

934 def summary_rates(self): 

935 irate = self.input_sample_rate 

936 orate = self.output_sample_rate 

937 factor = self.decimation_factor 

938 

939 if irate and orate is None: 

940 return '%g Hz' % irate 

941 

942 elif orate and irate is None: 

943 return '%g Hz' % orate 

944 

945 elif irate and orate and irate == orate: 

946 return '%g Hz' % irate 

947 

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

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

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

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

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

953 else: 

954 return '' 

955 

956 @property 

957 def summary_elements(self): 

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

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

960 

961 @property 

962 def summary_log(self): 

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

964 

965 @property 

966 def summary_entries(self): 

967 return ( 

968 self.__class__.__name__, 

969 self.stage_type, 

970 self.summary_log, 

971 self.summary_quantities, 

972 self.summary_rates, 

973 self.summary_elements) 

974 

975 @property 

976 def summary(self): 

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

978 

979 def get_effective(self): 

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

981 

982 

983D = 'displacement' 

984V = 'velocity' 

985A = 'acceleration' 

986 

987g_converters = { 

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

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

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

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

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

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

994 

995 

996def response_converters(input_quantity, output_quantity): 

997 if input_quantity is None or input_quantity == output_quantity: 

998 return [] 

999 

1000 if output_quantity is None: 

1001 raise ConversionError('Unspecified target quantity.') 

1002 

1003 try: 

1004 return [g_converters[input_quantity, output_quantity]] 

1005 

1006 except KeyError: 

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

1008 input_quantity, output_quantity)) 

1009 

1010 

1011@squirrel_content 

1012class Response(Object): 

1013 ''' 

1014 The instrument response of a seismic station channel. 

1015 ''' 

1016 

1017 codes = CodesNSLCE.T() 

1018 tmin = Timestamp.T(optional=True) 

1019 tmax = Timestamp.T(optional=True) 

1020 

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

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

1023 

1024 deltat = Float.T(optional=True) 

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

1026 

1027 def __init__(self, **kwargs): 

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

1029 Object.__init__(self, **kwargs) 

1030 self._effective_responses_cache = {} 

1031 

1032 @property 

1033 def time_span(self): 

1034 return (self.tmin, self.tmax) 

1035 

1036 @property 

1037 def nstages(self): 

1038 return len(self.stages) 

1039 

1040 @property 

1041 def input_quantity(self): 

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

1043 

1044 @property 

1045 def output_quantity(self): 

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

1047 

1048 @property 

1049 def output_sample_rate(self): 

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

1051 

1052 @property 

1053 def summary_stages(self): 

1054 def grouped(xs): 

1055 xs = list(xs) 

1056 g = [] 

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

1058 g.append(xs[i]) 

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

1060 yield g 

1061 g = [] 

1062 

1063 if g: 

1064 yield g 

1065 

1066 return '+'.join( 

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

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

1069 

1070 @property 

1071 def summary_quantities(self): 

1072 orate = self.output_sample_rate 

1073 return '%s -> %s%s' % ( 

1074 self.input_quantity or '?', 

1075 self.output_quantity or '?', 

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

1077 

1078 @property 

1079 def summary_log(self): 

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

1081 

1082 @property 

1083 def summary_entries(self): 

1084 return ( 

1085 self.__class__.__name__, 

1086 str(self.codes), 

1087 self.str_time_span, 

1088 self.summary_log, 

1089 self.summary_quantities, 

1090 self.summary_stages) 

1091 

1092 @property 

1093 def summary(self): 

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

1095 

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

1097 cache_key = (input_quantity, stages) 

1098 if cache_key in self._effective_responses_cache: 

1099 return self._effective_responses_cache[cache_key] 

1100 

1101 try: 

1102 elements = response_converters(input_quantity, self.input_quantity) 

1103 except ConversionError as e: 

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

1105 

1106 elements.extend( 

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

1108 

1109 if input_quantity is None \ 

1110 or input_quantity == self.input_quantity: 

1111 checkpoints = self.checkpoints 

1112 else: 

1113 checkpoints = [] 

1114 

1115 resp = MultiplyResponse( 

1116 responses=simplify_responses(elements), 

1117 checkpoints=checkpoints) 

1118 

1119 self._effective_responses_cache[cache_key] = resp 

1120 return resp 

1121 

1122 

1123@squirrel_content 

1124class Event(Object): 

1125 ''' 

1126 A seismic event. 

1127 ''' 

1128 

1129 name = String.T(optional=True) 

1130 time = Timestamp.T() 

1131 duration = Float.T(optional=True) 

1132 

1133 lat = Float.T() 

1134 lon = Float.T() 

1135 depth = Float.T(optional=True) 

1136 

1137 magnitude = Float.T(optional=True) 

1138 

1139 def get_pyrocko_event(self): 

1140 from pyrocko import model 

1141 model.Event( 

1142 name=self.name, 

1143 time=self.time, 

1144 lat=self.lat, 

1145 lon=self.lon, 

1146 depth=self.depth, 

1147 magnitude=self.magnitude, 

1148 duration=self.duration) 

1149 

1150 @property 

1151 def time_span(self): 

1152 return (self.time, self.time) 

1153 

1154 

1155def ehash(s): 

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

1157 

1158 

1159def random_name(n=8): 

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

1161 

1162 

1163@squirrel_content 

1164class WaveformPromise(Object): 

1165 ''' 

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

1167 

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

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

1170 calls to 

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

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

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

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

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

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

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

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

1179 queries for waveforms missing at the remote site. 

1180 ''' 

1181 

1182 codes = CodesNSLCE.T() 

1183 tmin = Timestamp.T() 

1184 tmax = Timestamp.T() 

1185 

1186 deltat = Float.T(optional=True) 

1187 

1188 source_hash = String.T() 

1189 

1190 def __init__(self, **kwargs): 

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

1192 Object.__init__(self, **kwargs) 

1193 

1194 @property 

1195 def time_span(self): 

1196 return (self.tmin, self.tmax) 

1197 

1198 

1199class InvalidWaveform(Exception): 

1200 pass 

1201 

1202 

1203class WaveformOrder(Object): 

1204 ''' 

1205 Waveform request information. 

1206 ''' 

1207 

1208 source_id = String.T() 

1209 codes = CodesNSLCE.T() 

1210 deltat = Float.T() 

1211 tmin = Timestamp.T() 

1212 tmax = Timestamp.T() 

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

1214 time_created = Timestamp.T() 

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

1216 

1217 @property 

1218 def client(self): 

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

1220 

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

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

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

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

1225 

1226 def validate(self, tr): 

1227 if tr.ydata.size == 0: 

1228 raise InvalidWaveform( 

1229 'waveform with zero data samples') 

1230 

1231 if tr.deltat != self.deltat: 

1232 raise InvalidWaveform( 

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

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

1235 tr.deltat, self.deltat)) 

1236 

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

1238 raise InvalidWaveform('waveform has NaN values') 

1239 

1240 def estimate_nsamples(self): 

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

1242 

1243 def is_near_real_time(self): 

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

1245 

1246 

1247def order_summary(orders): 

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

1249 if len(codes_list) > 3: 

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

1251 len(orders), 

1252 util.plural_s(orders), 

1253 str(codes_list[0]), 

1254 str(codes_list[1])) 

1255 

1256 else: 

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

1258 len(orders), 

1259 util.plural_s(orders), 

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

1261 

1262 

1263class Nut(Object): 

1264 ''' 

1265 Index entry referencing an elementary piece of content. 

1266 

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

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

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

1270 ''' 

1271 

1272 file_path = String.T(optional=True) 

1273 file_format = String.T(optional=True) 

1274 file_mtime = Timestamp.T(optional=True) 

1275 file_size = Int.T(optional=True) 

1276 

1277 file_segment = Int.T(optional=True) 

1278 file_element = Int.T(optional=True) 

1279 

1280 kind_id = Int.T() 

1281 codes = Codes.T() 

1282 

1283 tmin_seconds = Int.T(default=0) 

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

1285 tmax_seconds = Int.T(default=0) 

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

1287 

1288 deltat = Float.T(default=0.0) 

1289 

1290 content = Any.T(optional=True) 

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

1292 

1293 content_in_db = False 

1294 

1295 def __init__( 

1296 self, 

1297 file_path=None, 

1298 file_format=None, 

1299 file_mtime=None, 

1300 file_size=None, 

1301 file_segment=None, 

1302 file_element=None, 

1303 kind_id=0, 

1304 codes=CodesX(''), 

1305 tmin_seconds=None, 

1306 tmin_offset=0, 

1307 tmax_seconds=None, 

1308 tmax_offset=0, 

1309 deltat=None, 

1310 content=None, 

1311 raw_content=None, 

1312 tmin=None, 

1313 tmax=None, 

1314 values_nocheck=None): 

1315 

1316 if values_nocheck is not None: 

1317 (self.file_path, self.file_format, self.file_mtime, 

1318 self.file_size, 

1319 self.file_segment, self.file_element, 

1320 self.kind_id, codes_safe_str, 

1321 self.tmin_seconds, self.tmin_offset, 

1322 self.tmax_seconds, self.tmax_offset, 

1323 self.deltat) = values_nocheck 

1324 

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

1326 self.deltat = self.deltat or None 

1327 self.raw_content = {} 

1328 self.content = None 

1329 else: 

1330 if tmin is not None: 

1331 tmin_seconds, tmin_offset = tsplit(tmin) 

1332 

1333 if tmax is not None: 

1334 tmax_seconds, tmax_offset = tsplit(tmax) 

1335 

1336 self.kind_id = int(kind_id) 

1337 self.codes = codes 

1338 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1339 self.tmin_offset = int(tmin_offset) 

1340 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1341 self.tmax_offset = int(tmax_offset) 

1342 self.deltat = float_or_none(deltat) 

1343 self.file_path = str_or_none(file_path) 

1344 self.file_segment = int_or_none(file_segment) 

1345 self.file_element = int_or_none(file_element) 

1346 self.file_format = str_or_none(file_format) 

1347 self.file_mtime = float_or_none(file_mtime) 

1348 self.file_size = int_or_none(file_size) 

1349 self.content = content 

1350 if raw_content is None: 

1351 self.raw_content = {} 

1352 else: 

1353 self.raw_content = raw_content 

1354 

1355 Object.__init__(self, init_props=False) 

1356 

1357 def __eq__(self, other): 

1358 return (isinstance(other, Nut) and 

1359 self.equality_values == other.equality_values) 

1360 

1361 def hash(self): 

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

1363 

1364 def __ne__(self, other): 

1365 return not (self == other) 

1366 

1367 def get_io_backend(self): 

1368 from . import io 

1369 return io.get_backend(self.file_format) 

1370 

1371 def file_modified(self): 

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

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

1374 

1375 @property 

1376 def dkey(self): 

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

1378 

1379 @property 

1380 def key(self): 

1381 return ( 

1382 self.file_path, 

1383 self.file_segment, 

1384 self.file_element, 

1385 self.file_mtime) 

1386 

1387 @property 

1388 def equality_values(self): 

1389 return ( 

1390 self.file_segment, self.file_element, 

1391 self.kind_id, self.codes, 

1392 self.tmin_seconds, self.tmin_offset, 

1393 self.tmax_seconds, self.tmax_offset, self.deltat) 

1394 

1395 def diff(self, other): 

1396 names = [ 

1397 'file_segment', 'file_element', 'kind_id', 'codes', 

1398 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset', 

1399 'deltat'] 

1400 

1401 d = [] 

1402 for name, a, b in zip( 

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

1404 

1405 if a != b: 

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

1407 

1408 return d 

1409 

1410 @property 

1411 def tmin(self): 

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

1413 

1414 @tmin.setter 

1415 def tmin(self, t): 

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

1417 

1418 @property 

1419 def tmax(self): 

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

1421 

1422 @tmax.setter 

1423 def tmax(self, t): 

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

1425 

1426 @property 

1427 def kscale(self): 

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

1429 return 0 

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

1431 

1432 @property 

1433 def waveform_kwargs(self): 

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

1435 

1436 return dict( 

1437 network=network, 

1438 station=station, 

1439 location=location, 

1440 channel=channel, 

1441 extra=extra, 

1442 tmin=self.tmin, 

1443 tmax=self.tmax, 

1444 deltat=self.deltat) 

1445 

1446 @property 

1447 def waveform_promise_kwargs(self): 

1448 return dict( 

1449 codes=self.codes, 

1450 tmin=self.tmin, 

1451 tmax=self.tmax, 

1452 deltat=self.deltat) 

1453 

1454 @property 

1455 def station_kwargs(self): 

1456 network, station, location = self.codes 

1457 return dict( 

1458 codes=self.codes, 

1459 tmin=tmin_or_none(self.tmin), 

1460 tmax=tmax_or_none(self.tmax)) 

1461 

1462 @property 

1463 def channel_kwargs(self): 

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

1465 return dict( 

1466 codes=self.codes, 

1467 tmin=tmin_or_none(self.tmin), 

1468 tmax=tmax_or_none(self.tmax), 

1469 deltat=self.deltat) 

1470 

1471 @property 

1472 def response_kwargs(self): 

1473 return dict( 

1474 codes=self.codes, 

1475 tmin=tmin_or_none(self.tmin), 

1476 tmax=tmax_or_none(self.tmax), 

1477 deltat=self.deltat) 

1478 

1479 @property 

1480 def event_kwargs(self): 

1481 return dict( 

1482 name=self.codes, 

1483 time=self.tmin, 

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

1485 

1486 @property 

1487 def trace_kwargs(self): 

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

1489 

1490 return dict( 

1491 network=network, 

1492 station=station, 

1493 location=location, 

1494 channel=channel, 

1495 extra=extra, 

1496 tmin=self.tmin, 

1497 tmax=self.tmax-self.deltat, 

1498 deltat=self.deltat) 

1499 

1500 @property 

1501 def dummy_trace(self): 

1502 return DummyTrace(self) 

1503 

1504 @property 

1505 def summary_entries(self): 

1506 if self.tmin == self.tmax: 

1507 ts = util.time_to_str(self.tmin) 

1508 else: 

1509 ts = '%s - %s' % ( 

1510 util.time_to_str(self.tmin), 

1511 util.time_to_str(self.tmax)) 

1512 

1513 return ( 

1514 self.__class__.__name__, 

1515 to_kind(self.kind_id), 

1516 str(self.codes), 

1517 ts) 

1518 

1519 @property 

1520 def summary(self): 

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

1522 

1523 

1524def make_waveform_nut(**kwargs): 

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

1526 

1527 

1528def make_waveform_promise_nut(**kwargs): 

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

1530 

1531 

1532def make_station_nut(**kwargs): 

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

1534 

1535 

1536def make_channel_nut(**kwargs): 

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

1538 

1539 

1540def make_response_nut(**kwargs): 

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

1542 

1543 

1544def make_event_nut(**kwargs): 

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

1546 

1547 

1548def group_channels(nuts): 

1549 by_ansl = {} 

1550 for nut in nuts: 

1551 if nut.kind_id != CHANNEL: 

1552 continue 

1553 

1554 ansl = nut.codes[:4] 

1555 

1556 if ansl not in by_ansl: 

1557 by_ansl[ansl] = {} 

1558 

1559 group = by_ansl[ansl] 

1560 

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

1562 

1563 if k not in group: 

1564 group[k] = set() 

1565 

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

1567 

1568 return by_ansl 

1569 

1570 

1571class DummyTrace(object): 

1572 

1573 def __init__(self, nut): 

1574 self.nut = nut 

1575 self.codes = nut.codes 

1576 self.meta = {} 

1577 

1578 @property 

1579 def tmin(self): 

1580 return self.nut.tmin 

1581 

1582 @property 

1583 def tmax(self): 

1584 return self.nut.tmax 

1585 

1586 @property 

1587 def deltat(self): 

1588 return self.nut.deltat 

1589 

1590 @property 

1591 def nslc_id(self): 

1592 return self.codes.nslc 

1593 

1594 @property 

1595 def network(self): 

1596 return self.codes.network 

1597 

1598 @property 

1599 def station(self): 

1600 return self.codes.station 

1601 

1602 @property 

1603 def location(self): 

1604 return self.codes.location 

1605 

1606 @property 

1607 def channel(self): 

1608 return self.codes.channel 

1609 

1610 @property 

1611 def extra(self): 

1612 return self.codes.extra 

1613 

1614 def overlaps(self, tmin, tmax): 

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

1616 

1617 

1618def duration_to_str(t): 

1619 if t > 24*3600: 

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

1621 elif t > 3600: 

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

1623 elif t > 60: 

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

1625 else: 

1626 return '%gs' % t 

1627 

1628 

1629class Coverage(Object): 

1630 ''' 

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

1632 ''' 

1633 kind_id = Int.T( 

1634 help='Content type.') 

1635 pattern = Codes.T( 

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

1637 'match.') 

1638 codes = Codes.T( 

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

1640 deltat = Float.T( 

1641 help='Sampling interval [s]', 

1642 optional=True) 

1643 tmin = Timestamp.T( 

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

1645 optional=True) 

1646 tmax = Timestamp.T( 

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

1648 optional=True) 

1649 changes = List.T( 

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

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

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

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

1654 'value duplicate or redundant data coverage.') 

1655 

1656 @classmethod 

1657 def from_values(cls, args): 

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

1659 return cls( 

1660 kind_id=kind_id, 

1661 pattern=pattern, 

1662 codes=codes, 

1663 deltat=deltat, 

1664 tmin=tmin, 

1665 tmax=tmax, 

1666 changes=changes) 

1667 

1668 @property 

1669 def summary_entries(self): 

1670 ts = '%s - %s' % ( 

1671 util.time_to_str(self.tmin), 

1672 util.time_to_str(self.tmax)) 

1673 

1674 srate = self.sample_rate 

1675 

1676 total = self.total 

1677 

1678 return ( 

1679 self.__class__.__name__, 

1680 to_kind(self.kind_id), 

1681 str(self.codes), 

1682 ts, 

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

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

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

1686 

1687 @property 

1688 def summary(self): 

1689 return util.fmt_summary( 

1690 self.summary_entries, 

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

1692 

1693 @property 

1694 def sample_rate(self): 

1695 if self.deltat is None: 

1696 return None 

1697 elif self.deltat == 0.0: 

1698 return 0.0 

1699 else: 

1700 return 1.0 / self.deltat 

1701 

1702 @property 

1703 def labels(self): 

1704 srate = self.sample_rate 

1705 return ( 

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

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

1708 

1709 @property 

1710 def total(self): 

1711 total_t = None 

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

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

1714 

1715 return total_t 

1716 

1717 def iter_spans(self): 

1718 last = None 

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

1720 if last is not None: 

1721 last_t, last_count = last 

1722 if last_count > 0: 

1723 yield last_t, t, last_count 

1724 

1725 last = (t, count) 

1726 

1727 def iter_uncovered_by(self, other): 

1728 a = self 

1729 b = other 

1730 ia = ib = -1 

1731 ca = cb = 0 

1732 last = None 

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

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

1735 ia += 1 

1736 t, ca = a.changes[ia] 

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

1738 ib += 1 

1739 t, cb = b.changes[ib] 

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

1741 ia += 1 

1742 t, ca = a.changes[ia] 

1743 else: 

1744 ib += 1 

1745 t, cb = b.changes[ib] 

1746 

1747 if last is not None: 

1748 tl, cal, cbl = last 

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

1750 yield tl, t, ia, ib 

1751 

1752 last = (t, ca, cb) 

1753 

1754 def iter_uncovered_by_combined(self, other): 

1755 ib_last = None 

1756 group = None 

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

1758 if ib_last is None or ib != ib_last: 

1759 if group: 

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

1761 

1762 group = [] 

1763 

1764 group.append((tmin, tmax)) 

1765 ib_last = ib 

1766 

1767 if group: 

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

1769 

1770 

1771__all__ = [ 

1772 'to_codes', 

1773 'to_codes_guess', 

1774 'to_codes_simple', 

1775 'to_kind', 

1776 'to_kinds', 

1777 'to_kind_id', 

1778 'to_kind_ids', 

1779 'CodesError', 

1780 'Codes', 

1781 'CodesNSLCE', 

1782 'CodesNSL', 

1783 'CodesX', 

1784 'Station', 

1785 'Channel', 

1786 'Sensor', 

1787 'Response', 

1788 'Nut', 

1789 'Coverage', 

1790 'WaveformPromise', 

1791]