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, 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 

1199 @property 

1200 def client(self): 

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

1202 

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

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

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

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

1207 

1208 def validate(self, tr): 

1209 if tr.ydata.size == 0: 

1210 raise InvalidWaveform( 

1211 'waveform with zero data samples') 

1212 

1213 if tr.deltat != self.deltat: 

1214 raise InvalidWaveform( 

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

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

1217 tr.deltat, self.deltat)) 

1218 

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

1220 raise InvalidWaveform('waveform has NaN values') 

1221 

1222 

1223def order_summary(orders): 

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

1225 if len(codes_list) > 3: 

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

1227 len(orders), 

1228 util.plural_s(orders), 

1229 str(codes_list[0]), 

1230 str(codes_list[1])) 

1231 

1232 else: 

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

1234 len(orders), 

1235 util.plural_s(orders), 

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

1237 

1238 

1239class Nut(Object): 

1240 ''' 

1241 Index entry referencing an elementary piece of content. 

1242 

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

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

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

1246 ''' 

1247 

1248 file_path = String.T(optional=True) 

1249 file_format = String.T(optional=True) 

1250 file_mtime = Timestamp.T(optional=True) 

1251 file_size = Int.T(optional=True) 

1252 

1253 file_segment = Int.T(optional=True) 

1254 file_element = Int.T(optional=True) 

1255 

1256 kind_id = Int.T() 

1257 codes = Codes.T() 

1258 

1259 tmin_seconds = Int.T(default=0) 

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

1261 tmax_seconds = Int.T(default=0) 

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

1263 

1264 deltat = Float.T(default=0.0) 

1265 

1266 content = Any.T(optional=True) 

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

1268 

1269 content_in_db = False 

1270 

1271 def __init__( 

1272 self, 

1273 file_path=None, 

1274 file_format=None, 

1275 file_mtime=None, 

1276 file_size=None, 

1277 file_segment=None, 

1278 file_element=None, 

1279 kind_id=0, 

1280 codes=CodesX(''), 

1281 tmin_seconds=None, 

1282 tmin_offset=0, 

1283 tmax_seconds=None, 

1284 tmax_offset=0, 

1285 deltat=None, 

1286 content=None, 

1287 raw_content=None, 

1288 tmin=None, 

1289 tmax=None, 

1290 values_nocheck=None): 

1291 

1292 if values_nocheck is not None: 

1293 (self.file_path, self.file_format, self.file_mtime, 

1294 self.file_size, 

1295 self.file_segment, self.file_element, 

1296 self.kind_id, codes_safe_str, 

1297 self.tmin_seconds, self.tmin_offset, 

1298 self.tmax_seconds, self.tmax_offset, 

1299 self.deltat) = values_nocheck 

1300 

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

1302 self.deltat = self.deltat or None 

1303 self.raw_content = {} 

1304 self.content = None 

1305 else: 

1306 if tmin is not None: 

1307 tmin_seconds, tmin_offset = tsplit(tmin) 

1308 

1309 if tmax is not None: 

1310 tmax_seconds, tmax_offset = tsplit(tmax) 

1311 

1312 self.kind_id = int(kind_id) 

1313 self.codes = codes 

1314 self.tmin_seconds = int_or_g_tmin(tmin_seconds) 

1315 self.tmin_offset = int(tmin_offset) 

1316 self.tmax_seconds = int_or_g_tmax(tmax_seconds) 

1317 self.tmax_offset = int(tmax_offset) 

1318 self.deltat = float_or_none(deltat) 

1319 self.file_path = str_or_none(file_path) 

1320 self.file_segment = int_or_none(file_segment) 

1321 self.file_element = int_or_none(file_element) 

1322 self.file_format = str_or_none(file_format) 

1323 self.file_mtime = float_or_none(file_mtime) 

1324 self.file_size = int_or_none(file_size) 

1325 self.content = content 

1326 if raw_content is None: 

1327 self.raw_content = {} 

1328 else: 

1329 self.raw_content = raw_content 

1330 

1331 Object.__init__(self, init_props=False) 

1332 

1333 def __eq__(self, other): 

1334 return (isinstance(other, Nut) and 

1335 self.equality_values == other.equality_values) 

1336 

1337 def hash(self): 

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

1339 

1340 def __ne__(self, other): 

1341 return not (self == other) 

1342 

1343 def get_io_backend(self): 

1344 from . import io 

1345 return io.get_backend(self.file_format) 

1346 

1347 def file_modified(self): 

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

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

1350 

1351 @property 

1352 def dkey(self): 

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

1354 

1355 @property 

1356 def key(self): 

1357 return ( 

1358 self.file_path, 

1359 self.file_segment, 

1360 self.file_element, 

1361 self.file_mtime) 

1362 

1363 @property 

1364 def equality_values(self): 

1365 return ( 

1366 self.file_segment, self.file_element, 

1367 self.kind_id, self.codes, 

1368 self.tmin_seconds, self.tmin_offset, 

1369 self.tmax_seconds, self.tmax_offset, self.deltat) 

1370 

1371 def diff(self, other): 

1372 names = [ 

1373 'file_segment', 'file_element', 'kind_id', 'codes', 

1374 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset', 

1375 'deltat'] 

1376 

1377 d = [] 

1378 for name, a, b in zip( 

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

1380 

1381 if a != b: 

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

1383 

1384 return d 

1385 

1386 @property 

1387 def tmin(self): 

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

1389 

1390 @tmin.setter 

1391 def tmin(self, t): 

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

1393 

1394 @property 

1395 def tmax(self): 

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

1397 

1398 @tmax.setter 

1399 def tmax(self, t): 

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

1401 

1402 @property 

1403 def kscale(self): 

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

1405 return 0 

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

1407 

1408 @property 

1409 def waveform_kwargs(self): 

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

1411 

1412 return dict( 

1413 network=network, 

1414 station=station, 

1415 location=location, 

1416 channel=channel, 

1417 extra=extra, 

1418 tmin=self.tmin, 

1419 tmax=self.tmax, 

1420 deltat=self.deltat) 

1421 

1422 @property 

1423 def waveform_promise_kwargs(self): 

1424 return dict( 

1425 codes=self.codes, 

1426 tmin=self.tmin, 

1427 tmax=self.tmax, 

1428 deltat=self.deltat) 

1429 

1430 @property 

1431 def station_kwargs(self): 

1432 network, station, location = self.codes 

1433 return dict( 

1434 codes=self.codes, 

1435 tmin=tmin_or_none(self.tmin), 

1436 tmax=tmax_or_none(self.tmax)) 

1437 

1438 @property 

1439 def channel_kwargs(self): 

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

1441 return dict( 

1442 codes=self.codes, 

1443 tmin=tmin_or_none(self.tmin), 

1444 tmax=tmax_or_none(self.tmax), 

1445 deltat=self.deltat) 

1446 

1447 @property 

1448 def response_kwargs(self): 

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 event_kwargs(self): 

1457 return dict( 

1458 name=self.codes, 

1459 time=self.tmin, 

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

1461 

1462 @property 

1463 def trace_kwargs(self): 

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

1465 

1466 return dict( 

1467 network=network, 

1468 station=station, 

1469 location=location, 

1470 channel=channel, 

1471 extra=extra, 

1472 tmin=self.tmin, 

1473 tmax=self.tmax-self.deltat, 

1474 deltat=self.deltat) 

1475 

1476 @property 

1477 def dummy_trace(self): 

1478 return DummyTrace(self) 

1479 

1480 @property 

1481 def summary_entries(self): 

1482 if self.tmin == self.tmax: 

1483 ts = util.time_to_str(self.tmin) 

1484 else: 

1485 ts = '%s - %s' % ( 

1486 util.time_to_str(self.tmin), 

1487 util.time_to_str(self.tmax)) 

1488 

1489 return ( 

1490 self.__class__.__name__, 

1491 to_kind(self.kind_id), 

1492 str(self.codes), 

1493 ts) 

1494 

1495 @property 

1496 def summary(self): 

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

1498 

1499 

1500def make_waveform_nut(**kwargs): 

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

1502 

1503 

1504def make_waveform_promise_nut(**kwargs): 

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

1506 

1507 

1508def make_station_nut(**kwargs): 

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

1510 

1511 

1512def make_channel_nut(**kwargs): 

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

1514 

1515 

1516def make_response_nut(**kwargs): 

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

1518 

1519 

1520def make_event_nut(**kwargs): 

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

1522 

1523 

1524def group_channels(nuts): 

1525 by_ansl = {} 

1526 for nut in nuts: 

1527 if nut.kind_id != CHANNEL: 

1528 continue 

1529 

1530 ansl = nut.codes[:4] 

1531 

1532 if ansl not in by_ansl: 

1533 by_ansl[ansl] = {} 

1534 

1535 group = by_ansl[ansl] 

1536 

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

1538 

1539 if k not in group: 

1540 group[k] = set() 

1541 

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

1543 

1544 return by_ansl 

1545 

1546 

1547class DummyTrace(object): 

1548 

1549 def __init__(self, nut): 

1550 self.nut = nut 

1551 self.codes = nut.codes 

1552 self.meta = {} 

1553 

1554 @property 

1555 def tmin(self): 

1556 return self.nut.tmin 

1557 

1558 @property 

1559 def tmax(self): 

1560 return self.nut.tmax 

1561 

1562 @property 

1563 def deltat(self): 

1564 return self.nut.deltat 

1565 

1566 @property 

1567 def nslc_id(self): 

1568 return self.codes.nslc 

1569 

1570 @property 

1571 def network(self): 

1572 return self.codes.network 

1573 

1574 @property 

1575 def station(self): 

1576 return self.codes.station 

1577 

1578 @property 

1579 def location(self): 

1580 return self.codes.location 

1581 

1582 @property 

1583 def channel(self): 

1584 return self.codes.channel 

1585 

1586 @property 

1587 def extra(self): 

1588 return self.codes.extra 

1589 

1590 def overlaps(self, tmin, tmax): 

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

1592 

1593 

1594def duration_to_str(t): 

1595 if t > 24*3600: 

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

1597 elif t > 3600: 

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

1599 elif t > 60: 

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

1601 else: 

1602 return '%gs' % t 

1603 

1604 

1605class Coverage(Object): 

1606 ''' 

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

1608 ''' 

1609 kind_id = Int.T( 

1610 help='Content type.') 

1611 pattern = Codes.T( 

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

1613 'match.') 

1614 codes = Codes.T( 

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

1616 deltat = Float.T( 

1617 help='Sampling interval [s]', 

1618 optional=True) 

1619 tmin = Timestamp.T( 

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

1621 optional=True) 

1622 tmax = Timestamp.T( 

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

1624 optional=True) 

1625 changes = List.T( 

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

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

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

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

1630 'value duplicate or redundant data coverage.') 

1631 

1632 @classmethod 

1633 def from_values(cls, args): 

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

1635 return cls( 

1636 kind_id=kind_id, 

1637 pattern=pattern, 

1638 codes=codes, 

1639 deltat=deltat, 

1640 tmin=tmin, 

1641 tmax=tmax, 

1642 changes=changes) 

1643 

1644 @property 

1645 def summary_entries(self): 

1646 ts = '%s - %s' % ( 

1647 util.time_to_str(self.tmin), 

1648 util.time_to_str(self.tmax)) 

1649 

1650 srate = self.sample_rate 

1651 

1652 total = self.total 

1653 

1654 return ( 

1655 self.__class__.__name__, 

1656 to_kind(self.kind_id), 

1657 str(self.codes), 

1658 ts, 

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

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

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

1662 

1663 @property 

1664 def summary(self): 

1665 return util.fmt_summary( 

1666 self.summary_entries, 

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

1668 

1669 @property 

1670 def sample_rate(self): 

1671 if self.deltat is None: 

1672 return None 

1673 elif self.deltat == 0.0: 

1674 return 0.0 

1675 else: 

1676 return 1.0 / self.deltat 

1677 

1678 @property 

1679 def labels(self): 

1680 srate = self.sample_rate 

1681 return ( 

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

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

1684 

1685 @property 

1686 def total(self): 

1687 total_t = None 

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

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

1690 

1691 return total_t 

1692 

1693 def iter_spans(self): 

1694 last = None 

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

1696 if last is not None: 

1697 last_t, last_count = last 

1698 if last_count > 0: 

1699 yield last_t, t, last_count 

1700 

1701 last = (t, count) 

1702 

1703 def iter_uncovered_by(self, other): 

1704 a = self 

1705 b = other 

1706 ia = ib = -1 

1707 ca = cb = 0 

1708 last = None 

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

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

1711 ia += 1 

1712 t, ca = a.changes[ia] 

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

1714 ib += 1 

1715 t, cb = b.changes[ib] 

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

1717 ia += 1 

1718 t, ca = a.changes[ia] 

1719 else: 

1720 ib += 1 

1721 t, cb = b.changes[ib] 

1722 

1723 if last is not None: 

1724 tl, cal, cbl = last 

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

1726 yield tl, t, ia, ib 

1727 

1728 last = (t, ca, cb) 

1729 

1730 def iter_uncovered_by_combined(self, other): 

1731 ib_last = None 

1732 group = None 

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

1734 if ib_last is None or ib != ib_last: 

1735 if group: 

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

1737 

1738 group = [] 

1739 

1740 group.append((tmin, tmax)) 

1741 ib_last = ib 

1742 

1743 if group: 

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

1745 

1746 

1747__all__ = [ 

1748 'to_codes', 

1749 'to_codes_guess', 

1750 'to_codes_simple', 

1751 'to_kind', 

1752 'to_kinds', 

1753 'to_kind_id', 

1754 'to_kind_ids', 

1755 'CodesError', 

1756 'Codes', 

1757 'CodesNSLCE', 

1758 'CodesNSL', 

1759 'CodesX', 

1760 'Station', 

1761 'Channel', 

1762 'Sensor', 

1763 'Response', 

1764 'Nut', 

1765 'Coverage', 

1766 'WaveformPromise', 

1767]