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 from pyrocko import model 

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

629 

630 def _get_pyrocko_station_args(self): 

631 return ( 

632 self.codes.network, 

633 self.codes.station, 

634 self.codes.location, 

635 self.lat, 

636 self.lon, 

637 self.elevation, 

638 self.depth, 

639 self.north_shift, 

640 self.east_shift) 

641 

642 

643@squirrel_content 

644class ChannelBase(Location): 

645 codes = CodesNSLCE.T() 

646 

647 tmin = Timestamp.T(optional=True) 

648 tmax = Timestamp.T(optional=True) 

649 

650 deltat = Float.T(optional=True) 

651 

652 def __init__(self, **kwargs): 

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

654 Location.__init__(self, **kwargs) 

655 

656 @property 

657 def time_span(self): 

658 return (self.tmin, self.tmax) 

659 

660 def _get_sensor_codes(self): 

661 return self.codes.replace( 

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

663 

664 def _get_sensor_args(self): 

665 def getattr_rep(k): 

666 if k == 'codes': 

667 return self._get_sensor_codes() 

668 else: 

669 return getattr(self, k) 

670 

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

672 

673 def _get_channel_args(self, component): 

674 def getattr_rep(k): 

675 if k == 'codes': 

676 return self.codes.replace( 

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

678 else: 

679 return getattr(self, k) 

680 

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

682 

683 def _get_pyrocko_station_args(self): 

684 return ( 

685 self.codes.network, 

686 self.codes.station, 

687 self.codes.location, 

688 self.lat, 

689 self.lon, 

690 self.elevation, 

691 self.depth, 

692 self.north_shift, 

693 self.east_shift) 

694 

695 

696class Channel(ChannelBase): 

697 ''' 

698 A channel of a seismic station. 

699 ''' 

700 

701 dip = Float.T(optional=True) 

702 azimuth = Float.T(optional=True) 

703 

704 def get_pyrocko_channel(self): 

705 from pyrocko import model 

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

707 

708 def _get_pyrocko_channel_args(self): 

709 return ( 

710 self.codes.channel, 

711 self.azimuth, 

712 self.dip) 

713 

714 @property 

715 def orientation_enz(self): 

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

717 return None 

718 

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

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

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

722 return mkvec(e, n, -d) 

723 

724 

725def cut_intervals(channels): 

726 channels = list(channels) 

727 times = set() 

728 for channel in channels: 

729 if channel.tmin is not None: 

730 times.add(channel.tmin) 

731 if channel.tmax is not None: 

732 times.add(channel.tmax) 

733 

734 times = sorted(times) 

735 

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

737 times[0:0] = [None] 

738 

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

740 times.append(None) 

741 

742 if len(times) <= 2: 

743 return channels 

744 

745 channels_out = [] 

746 for channel in channels: 

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

748 tmin = times[i] 

749 tmax = times[i+1] 

750 if ((channel.tmin is None or ( 

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

752 and (channel.tmax is None or ( 

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

754 

755 channel_out = clone(channel) 

756 channel_out.tmin = tmin 

757 channel_out.tmax = tmax 

758 channels_out.append(channel_out) 

759 

760 return channels_out 

761 

762 

763class Sensor(ChannelBase): 

764 ''' 

765 Representation of a channel group. 

766 ''' 

767 

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

769 

770 @classmethod 

771 def from_channels(cls, channels): 

772 groups = defaultdict(list) 

773 for channel in channels: 

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

775 

776 channels_cut = [] 

777 for group in groups.values(): 

778 channels_cut.extend(cut_intervals(group)) 

779 

780 groups = defaultdict(list) 

781 for channel in channels_cut: 

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

783 

784 return [ 

785 cls(channels=channels, 

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

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

788 

789 def channel_vectors(self): 

790 return num.vstack( 

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

792 

793 def projected_channels(self, system, component_names): 

794 return [ 

795 Channel( 

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

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

798 **dict(zip( 

799 ChannelBase.T.propnames, 

800 self._get_channel_args(comp)))) 

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

802 

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

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

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

806 return m 

807 

808 def projection_to(self, system, component_names): 

809 return ( 

810 self.matrix_to(system), 

811 self.channels, 

812 self.projected_channels(system, component_names)) 

813 

814 def projection_to_enz(self): 

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

816 

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

818 if azimuth is not None: 

819 assert source is None 

820 else: 

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

822 

823 return self.projection_to(num.array([ 

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

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

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

827 

828 def project_to_enz(self, traces): 

829 from pyrocko import trace 

830 

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

832 

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

834 

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

836 from pyrocko import trace 

837 

838 matrix, in_channels, out_channels = self.projection_to_trz( 

839 source, azimuth=azimuth) 

840 

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

842 

843 

844observational_quantities = [ 

845 'acceleration', 'velocity', 'displacement', 'pressure', 

846 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration', 

847 'temperature'] 

848 

849 

850technical_quantities = [ 

851 'voltage', 'counts'] 

852 

853 

854class QuantityType(StringChoice): 

855 ''' 

856 Choice of observational or technical quantity. 

857 

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

859 ''' 

860 choices = observational_quantities + technical_quantities 

861 

862 

863class ResponseStage(Object): 

864 ''' 

865 Representation of a response stage. 

866 

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

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

869 downsampling. 

870 ''' 

871 input_quantity = QuantityType.T(optional=True) 

872 input_sample_rate = Float.T(optional=True) 

873 output_quantity = QuantityType.T(optional=True) 

874 output_sample_rate = Float.T(optional=True) 

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

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

877 

878 @property 

879 def stage_type(self): 

880 if self.input_quantity in observational_quantities \ 

881 and self.output_quantity in observational_quantities: 

882 return 'conversion' 

883 

884 if self.input_quantity in observational_quantities \ 

885 and self.output_quantity == 'voltage': 

886 return 'sensor' 

887 

888 elif self.input_quantity == 'voltage' \ 

889 and self.output_quantity == 'voltage': 

890 return 'electronics' 

891 

892 elif self.input_quantity == 'voltage' \ 

893 and self.output_quantity == 'counts': 

894 return 'digitizer' 

895 

896 elif self.decimation_factor is not None \ 

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

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

899 return 'decimation' 

900 

901 elif self.input_quantity in observational_quantities \ 

902 and self.output_quantity == 'counts': 

903 return 'combination' 

904 

905 else: 

906 return 'unknown' 

907 

908 @property 

909 def decimation_factor(self): 

910 irate = self.input_sample_rate 

911 orate = self.output_sample_rate 

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

913 and irate > orate and irate / orate > 1.0: 

914 

915 return irate / orate 

916 else: 

917 return None 

918 

919 @property 

920 def summary_quantities(self): 

921 return '%s -> %s' % ( 

922 self.input_quantity or '?', 

923 self.output_quantity or '?') 

924 

925 @property 

926 def summary_rates(self): 

927 irate = self.input_sample_rate 

928 orate = self.output_sample_rate 

929 factor = self.decimation_factor 

930 

931 if irate and orate is None: 

932 return '%g Hz' % irate 

933 

934 elif orate and irate is None: 

935 return '%g Hz' % orate 

936 

937 elif irate and orate and irate == orate: 

938 return '%g Hz' % irate 

939 

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

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

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

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

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

945 else: 

946 return '' 

947 

948 @property 

949 def summary_elements(self): 

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

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

952 

953 @property 

954 def summary_log(self): 

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

956 

957 @property 

958 def summary_entries(self): 

959 return ( 

960 self.__class__.__name__, 

961 self.stage_type, 

962 self.summary_log, 

963 self.summary_quantities, 

964 self.summary_rates, 

965 self.summary_elements) 

966 

967 @property 

968 def summary(self): 

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

970 

971 def get_effective(self): 

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

973 

974 

975D = 'displacement' 

976V = 'velocity' 

977A = 'acceleration' 

978 

979g_converters = { 

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

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

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

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

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

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

986 

987 

988def response_converters(input_quantity, output_quantity): 

989 if input_quantity is None or input_quantity == output_quantity: 

990 return [] 

991 

992 if output_quantity is None: 

993 raise ConversionError('Unspecified target quantity.') 

994 

995 try: 

996 return [g_converters[input_quantity, output_quantity]] 

997 

998 except KeyError: 

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

1000 input_quantity, output_quantity)) 

1001 

1002 

1003@squirrel_content 

1004class Response(Object): 

1005 ''' 

1006 The instrument response of a seismic station channel. 

1007 ''' 

1008 

1009 codes = CodesNSLCE.T() 

1010 tmin = Timestamp.T(optional=True) 

1011 tmax = Timestamp.T(optional=True) 

1012 

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

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

1015 

1016 deltat = Float.T(optional=True) 

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

1018 

1019 def __init__(self, **kwargs): 

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

1021 Object.__init__(self, **kwargs) 

1022 

1023 @property 

1024 def time_span(self): 

1025 return (self.tmin, self.tmax) 

1026 

1027 @property 

1028 def nstages(self): 

1029 return len(self.stages) 

1030 

1031 @property 

1032 def input_quantity(self): 

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

1034 

1035 @property 

1036 def output_quantity(self): 

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

1038 

1039 @property 

1040 def output_sample_rate(self): 

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

1042 

1043 @property 

1044 def summary_stages(self): 

1045 def grouped(xs): 

1046 xs = list(xs) 

1047 g = [] 

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

1049 g.append(xs[i]) 

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

1051 yield g 

1052 g = [] 

1053 

1054 if g: 

1055 yield g 

1056 

1057 return '+'.join( 

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

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

1060 

1061 @property 

1062 def summary_quantities(self): 

1063 orate = self.output_sample_rate 

1064 return '%s -> %s%s' % ( 

1065 self.input_quantity or '?', 

1066 self.output_quantity or '?', 

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

1068 

1069 @property 

1070 def summary_log(self): 

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

1072 

1073 @property 

1074 def summary_entries(self): 

1075 return ( 

1076 self.__class__.__name__, 

1077 str(self.codes), 

1078 self.str_time_span, 

1079 self.summary_log, 

1080 self.summary_quantities, 

1081 self.summary_stages) 

1082 

1083 @property 

1084 def summary(self): 

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

1086 

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

1088 try: 

1089 elements = response_converters(input_quantity, self.input_quantity) 

1090 except ConversionError as e: 

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

1092 

1093 elements.extend( 

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

1095 

1096 if input_quantity is None \ 

1097 or input_quantity == self.input_quantity: 

1098 checkpoints = self.checkpoints 

1099 else: 

1100 checkpoints = [] 

1101 

1102 return MultiplyResponse( 

1103 responses=simplify_responses(elements), 

1104 checkpoints=checkpoints) 

1105 

1106 

1107@squirrel_content 

1108class Event(Object): 

1109 ''' 

1110 A seismic event. 

1111 ''' 

1112 

1113 name = String.T(optional=True) 

1114 time = Timestamp.T() 

1115 duration = Float.T(optional=True) 

1116 

1117 lat = Float.T() 

1118 lon = Float.T() 

1119 depth = Float.T(optional=True) 

1120 

1121 magnitude = Float.T(optional=True) 

1122 

1123 def get_pyrocko_event(self): 

1124 from pyrocko import model 

1125 model.Event( 

1126 name=self.name, 

1127 time=self.time, 

1128 lat=self.lat, 

1129 lon=self.lon, 

1130 depth=self.depth, 

1131 magnitude=self.magnitude, 

1132 duration=self.duration) 

1133 

1134 @property 

1135 def time_span(self): 

1136 return (self.time, self.time) 

1137 

1138 

1139def ehash(s): 

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

1141 

1142 

1143def random_name(n=8): 

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

1145 

1146 

1147@squirrel_content 

1148class WaveformPromise(Object): 

1149 ''' 

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

1151 

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

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

1154 calls to 

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

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

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

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

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

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

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

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

1163 queries for waveforms missing at the remote site. 

1164 ''' 

1165 

1166 codes = CodesNSLCE.T() 

1167 tmin = Timestamp.T() 

1168 tmax = Timestamp.T() 

1169 

1170 deltat = Float.T(optional=True) 

1171 

1172 source_hash = String.T() 

1173 

1174 def __init__(self, **kwargs): 

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

1176 Object.__init__(self, **kwargs) 

1177 

1178 @property 

1179 def time_span(self): 

1180 return (self.tmin, self.tmax) 

1181 

1182 

1183class InvalidWaveform(Exception): 

1184 pass 

1185 

1186 

1187class WaveformOrder(Object): 

1188 ''' 

1189 Waveform request information. 

1190 ''' 

1191 

1192 source_id = String.T() 

1193 codes = CodesNSLCE.T() 

1194 deltat = Float.T() 

1195 tmin = Timestamp.T() 

1196 tmax = Timestamp.T() 

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

1198 time_created = Timestamp.T() 

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

1200 

1201 @property 

1202 def client(self): 

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

1204 

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

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

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

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

1209 

1210 def validate(self, tr): 

1211 if tr.ydata.size == 0: 

1212 raise InvalidWaveform( 

1213 'waveform with zero data samples') 

1214 

1215 if tr.deltat != self.deltat: 

1216 raise InvalidWaveform( 

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

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

1219 tr.deltat, self.deltat)) 

1220 

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

1222 raise InvalidWaveform('waveform has NaN values') 

1223 

1224 def estimate_nsamples(self): 

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

1226 

1227 def is_near_real_time(self): 

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

1229 

1230 

1231def order_summary(orders): 

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

1233 if len(codes_list) > 3: 

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

1235 len(orders), 

1236 util.plural_s(orders), 

1237 str(codes_list[0]), 

1238 str(codes_list[1])) 

1239 

1240 else: 

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

1242 len(orders), 

1243 util.plural_s(orders), 

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

1245 

1246 

1247class Nut(Object): 

1248 ''' 

1249 Index entry referencing an elementary piece of content. 

1250 

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

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

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

1254 ''' 

1255 

1256 file_path = String.T(optional=True) 

1257 file_format = String.T(optional=True) 

1258 file_mtime = Timestamp.T(optional=True) 

1259 file_size = Int.T(optional=True) 

1260 

1261 file_segment = Int.T(optional=True) 

1262 file_element = Int.T(optional=True) 

1263 

1264 kind_id = Int.T() 

1265 codes = Codes.T() 

1266 

1267 tmin_seconds = Int.T(default=0) 

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

1269 tmax_seconds = Int.T(default=0) 

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

1271 

1272 deltat = Float.T(default=0.0) 

1273 

1274 content = Any.T(optional=True) 

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

1276 

1277 content_in_db = False 

1278 

1279 def __init__( 

1280 self, 

1281 file_path=None, 

1282 file_format=None, 

1283 file_mtime=None, 

1284 file_size=None, 

1285 file_segment=None, 

1286 file_element=None, 

1287 kind_id=0, 

1288 codes=CodesX(''), 

1289 tmin_seconds=None, 

1290 tmin_offset=0, 

1291 tmax_seconds=None, 

1292 tmax_offset=0, 

1293 deltat=None, 

1294 content=None, 

1295 raw_content=None, 

1296 tmin=None, 

1297 tmax=None, 

1298 values_nocheck=None): 

1299 

1300 if values_nocheck is not None: 

1301 (self.file_path, self.file_format, self.file_mtime, 

1302 self.file_size, 

1303 self.file_segment, self.file_element, 

1304 self.kind_id, codes_safe_str, 

1305 self.tmin_seconds, self.tmin_offset, 

1306 self.tmax_seconds, self.tmax_offset, 

1307 self.deltat) = values_nocheck 

1308 

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

1310 self.deltat = self.deltat or None 

1311 self.raw_content = {} 

1312 self.content = None 

1313 else: 

1314 if tmin is not None: 

1315 tmin_seconds, tmin_offset = tsplit(tmin) 

1316 

1317 if tmax is not None: 

1318 tmax_seconds, tmax_offset = tsplit(tmax) 

1319 

1320 self.kind_id = int(kind_id) 

1321 self.codes = codes 

1322 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1323 self.tmin_offset = int(tmin_offset) 

1324 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1325 self.tmax_offset = int(tmax_offset) 

1326 self.deltat = float_or_none(deltat) 

1327 self.file_path = str_or_none(file_path) 

1328 self.file_segment = int_or_none(file_segment) 

1329 self.file_element = int_or_none(file_element) 

1330 self.file_format = str_or_none(file_format) 

1331 self.file_mtime = float_or_none(file_mtime) 

1332 self.file_size = int_or_none(file_size) 

1333 self.content = content 

1334 if raw_content is None: 

1335 self.raw_content = {} 

1336 else: 

1337 self.raw_content = raw_content 

1338 

1339 Object.__init__(self, init_props=False) 

1340 

1341 def __eq__(self, other): 

1342 return (isinstance(other, Nut) and 

1343 self.equality_values == other.equality_values) 

1344 

1345 def hash(self): 

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

1347 

1348 def __ne__(self, other): 

1349 return not (self == other) 

1350 

1351 def get_io_backend(self): 

1352 from . import io 

1353 return io.get_backend(self.file_format) 

1354 

1355 def file_modified(self): 

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

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

1358 

1359 @property 

1360 def dkey(self): 

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

1362 

1363 @property 

1364 def key(self): 

1365 return ( 

1366 self.file_path, 

1367 self.file_segment, 

1368 self.file_element, 

1369 self.file_mtime) 

1370 

1371 @property 

1372 def equality_values(self): 

1373 return ( 

1374 self.file_segment, self.file_element, 

1375 self.kind_id, self.codes, 

1376 self.tmin_seconds, self.tmin_offset, 

1377 self.tmax_seconds, self.tmax_offset, self.deltat) 

1378 

1379 def diff(self, other): 

1380 names = [ 

1381 'file_segment', 'file_element', 'kind_id', 'codes', 

1382 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset', 

1383 'deltat'] 

1384 

1385 d = [] 

1386 for name, a, b in zip( 

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

1388 

1389 if a != b: 

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

1391 

1392 return d 

1393 

1394 @property 

1395 def tmin(self): 

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

1397 

1398 @tmin.setter 

1399 def tmin(self, t): 

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

1401 

1402 @property 

1403 def tmax(self): 

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

1405 

1406 @tmax.setter 

1407 def tmax(self, t): 

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

1409 

1410 @property 

1411 def kscale(self): 

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

1413 return 0 

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

1415 

1416 @property 

1417 def waveform_kwargs(self): 

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

1419 

1420 return dict( 

1421 network=network, 

1422 station=station, 

1423 location=location, 

1424 channel=channel, 

1425 extra=extra, 

1426 tmin=self.tmin, 

1427 tmax=self.tmax, 

1428 deltat=self.deltat) 

1429 

1430 @property 

1431 def waveform_promise_kwargs(self): 

1432 return dict( 

1433 codes=self.codes, 

1434 tmin=self.tmin, 

1435 tmax=self.tmax, 

1436 deltat=self.deltat) 

1437 

1438 @property 

1439 def station_kwargs(self): 

1440 network, station, location = self.codes 

1441 return dict( 

1442 codes=self.codes, 

1443 tmin=tmin_or_none(self.tmin), 

1444 tmax=tmax_or_none(self.tmax)) 

1445 

1446 @property 

1447 def channel_kwargs(self): 

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

1449 return dict( 

1450 codes=self.codes, 

1451 tmin=tmin_or_none(self.tmin), 

1452 tmax=tmax_or_none(self.tmax), 

1453 deltat=self.deltat) 

1454 

1455 @property 

1456 def response_kwargs(self): 

1457 return dict( 

1458 codes=self.codes, 

1459 tmin=tmin_or_none(self.tmin), 

1460 tmax=tmax_or_none(self.tmax), 

1461 deltat=self.deltat) 

1462 

1463 @property 

1464 def event_kwargs(self): 

1465 return dict( 

1466 name=self.codes, 

1467 time=self.tmin, 

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

1469 

1470 @property 

1471 def trace_kwargs(self): 

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

1473 

1474 return dict( 

1475 network=network, 

1476 station=station, 

1477 location=location, 

1478 channel=channel, 

1479 extra=extra, 

1480 tmin=self.tmin, 

1481 tmax=self.tmax-self.deltat, 

1482 deltat=self.deltat) 

1483 

1484 @property 

1485 def dummy_trace(self): 

1486 return DummyTrace(self) 

1487 

1488 @property 

1489 def summary_entries(self): 

1490 if self.tmin == self.tmax: 

1491 ts = util.time_to_str(self.tmin) 

1492 else: 

1493 ts = '%s - %s' % ( 

1494 util.time_to_str(self.tmin), 

1495 util.time_to_str(self.tmax)) 

1496 

1497 return ( 

1498 self.__class__.__name__, 

1499 to_kind(self.kind_id), 

1500 str(self.codes), 

1501 ts) 

1502 

1503 @property 

1504 def summary(self): 

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

1506 

1507 

1508def make_waveform_nut(**kwargs): 

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

1510 

1511 

1512def make_waveform_promise_nut(**kwargs): 

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

1514 

1515 

1516def make_station_nut(**kwargs): 

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

1518 

1519 

1520def make_channel_nut(**kwargs): 

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

1522 

1523 

1524def make_response_nut(**kwargs): 

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

1526 

1527 

1528def make_event_nut(**kwargs): 

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

1530 

1531 

1532def group_channels(nuts): 

1533 by_ansl = {} 

1534 for nut in nuts: 

1535 if nut.kind_id != CHANNEL: 

1536 continue 

1537 

1538 ansl = nut.codes[:4] 

1539 

1540 if ansl not in by_ansl: 

1541 by_ansl[ansl] = {} 

1542 

1543 group = by_ansl[ansl] 

1544 

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

1546 

1547 if k not in group: 

1548 group[k] = set() 

1549 

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

1551 

1552 return by_ansl 

1553 

1554 

1555class DummyTrace(object): 

1556 

1557 def __init__(self, nut): 

1558 self.nut = nut 

1559 self.codes = nut.codes 

1560 self.meta = {} 

1561 

1562 @property 

1563 def tmin(self): 

1564 return self.nut.tmin 

1565 

1566 @property 

1567 def tmax(self): 

1568 return self.nut.tmax 

1569 

1570 @property 

1571 def deltat(self): 

1572 return self.nut.deltat 

1573 

1574 @property 

1575 def nslc_id(self): 

1576 return self.codes.nslc 

1577 

1578 @property 

1579 def network(self): 

1580 return self.codes.network 

1581 

1582 @property 

1583 def station(self): 

1584 return self.codes.station 

1585 

1586 @property 

1587 def location(self): 

1588 return self.codes.location 

1589 

1590 @property 

1591 def channel(self): 

1592 return self.codes.channel 

1593 

1594 @property 

1595 def extra(self): 

1596 return self.codes.extra 

1597 

1598 def overlaps(self, tmin, tmax): 

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

1600 

1601 

1602def duration_to_str(t): 

1603 if t > 24*3600: 

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

1605 elif t > 3600: 

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

1607 elif t > 60: 

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

1609 else: 

1610 return '%gs' % t 

1611 

1612 

1613class Coverage(Object): 

1614 ''' 

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

1616 ''' 

1617 kind_id = Int.T( 

1618 help='Content type.') 

1619 pattern = Codes.T( 

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

1621 'match.') 

1622 codes = Codes.T( 

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

1624 deltat = Float.T( 

1625 help='Sampling interval [s]', 

1626 optional=True) 

1627 tmin = Timestamp.T( 

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

1629 optional=True) 

1630 tmax = Timestamp.T( 

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

1632 optional=True) 

1633 changes = List.T( 

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

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

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

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

1638 'value duplicate or redundant data coverage.') 

1639 

1640 @classmethod 

1641 def from_values(cls, args): 

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

1643 return cls( 

1644 kind_id=kind_id, 

1645 pattern=pattern, 

1646 codes=codes, 

1647 deltat=deltat, 

1648 tmin=tmin, 

1649 tmax=tmax, 

1650 changes=changes) 

1651 

1652 @property 

1653 def summary_entries(self): 

1654 ts = '%s - %s' % ( 

1655 util.time_to_str(self.tmin), 

1656 util.time_to_str(self.tmax)) 

1657 

1658 srate = self.sample_rate 

1659 

1660 total = self.total 

1661 

1662 return ( 

1663 self.__class__.__name__, 

1664 to_kind(self.kind_id), 

1665 str(self.codes), 

1666 ts, 

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

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

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

1670 

1671 @property 

1672 def summary(self): 

1673 return util.fmt_summary( 

1674 self.summary_entries, 

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

1676 

1677 @property 

1678 def sample_rate(self): 

1679 if self.deltat is None: 

1680 return None 

1681 elif self.deltat == 0.0: 

1682 return 0.0 

1683 else: 

1684 return 1.0 / self.deltat 

1685 

1686 @property 

1687 def labels(self): 

1688 srate = self.sample_rate 

1689 return ( 

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

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

1692 

1693 @property 

1694 def total(self): 

1695 total_t = None 

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

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

1698 

1699 return total_t 

1700 

1701 def iter_spans(self): 

1702 last = None 

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

1704 if last is not None: 

1705 last_t, last_count = last 

1706 if last_count > 0: 

1707 yield last_t, t, last_count 

1708 

1709 last = (t, count) 

1710 

1711 def iter_uncovered_by(self, other): 

1712 a = self 

1713 b = other 

1714 ia = ib = -1 

1715 ca = cb = 0 

1716 last = None 

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

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

1719 ia += 1 

1720 t, ca = a.changes[ia] 

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

1722 ib += 1 

1723 t, cb = b.changes[ib] 

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

1725 ia += 1 

1726 t, ca = a.changes[ia] 

1727 else: 

1728 ib += 1 

1729 t, cb = b.changes[ib] 

1730 

1731 if last is not None: 

1732 tl, cal, cbl = last 

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

1734 yield tl, t, ia, ib 

1735 

1736 last = (t, ca, cb) 

1737 

1738 def iter_uncovered_by_combined(self, other): 

1739 ib_last = None 

1740 group = None 

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

1742 if ib_last is None or ib != ib_last: 

1743 if group: 

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

1745 

1746 group = [] 

1747 

1748 group.append((tmin, tmax)) 

1749 ib_last = ib 

1750 

1751 if group: 

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

1753 

1754 

1755__all__ = [ 

1756 'to_codes', 

1757 'to_codes_guess', 

1758 'to_codes_simple', 

1759 'to_kind', 

1760 'to_kinds', 

1761 'to_kind_id', 

1762 'to_kind_ids', 

1763 'CodesError', 

1764 'Codes', 

1765 'CodesNSLCE', 

1766 'CodesNSL', 

1767 'CodesX', 

1768 'Station', 

1769 'Channel', 

1770 'Sensor', 

1771 'Response', 

1772 'Nut', 

1773 'Coverage', 

1774 'WaveformPromise', 

1775]