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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

4# ---|P------/S----------~Lg----------

6'''

7Data model and content types handled by the Squirrel framework.

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.

15Further content types are defined here to handle waveform orders, waveform

16promises, data coverage and sensor information.

17'''

19import re

20import fnmatch

21import math

22import hashlib

23from os import urandom

24from base64 import urlsafe_b64encode

25from collections import defaultdict, namedtuple

27import numpy as num

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

37from .error import ConversionError, SquirrelError

39d2r = num.pi / 180.

40r2d = 1.0 / d2r

43guts_prefix = 'squirrel'

46def mkvec(x, y, z):

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

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)]))

56g_codes_pool = {}

59class CodesError(SquirrelError):

60 pass

63class Codes(SObject):

64 pass

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

68 if args and kwargs:

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

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]))

79 nargs = len(args)

80 if nargs == 5:

81 t = args

83 elif nargs == 4:

84 t = args + ('',)

86 elif nargs == 0:

87 d = dict(

88 network='',

89 station='',

90 location='',

91 channel='',

92 extra='')

94 d.update(kwargs)

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

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

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 _is_exact(pat):

367 return not ('*' in pat or '?' in pat or ']' in pat or '[' in pat)

368

369

370def classify_patterns(patterns):

371 pats_exact = []

372 pats_nonexact = []

373 for pat in patterns:

374 spat = pat.safe_str

375 (pats_exact if _is_exact(spat) else pats_nonexact).append(spat)

376

377 return pats_exact, pats_nonexact

378

379

380def get_regex_pattern(spattern):

381 if spattern not in g_codes_patterns:

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

383 g_codes_patterns[spattern] = rpattern

384

385 return g_codes_patterns[spattern]

386

387

388def match_codes(pattern, codes):

389 spattern = pattern.safe_str

390 scodes = codes.safe_str

391 rpattern = get_regex_pattern(spattern)

392 return bool(rpattern.match(scodes))

393

394

395class CodesMatcher:

396 def __init__(self, patterns):

397 self._pats_exact, self._pats_nonexact = classify_patterns(patterns)

398 self._pats_exact = set(self._pats_exact)

399 self._pats_nonexact = [

400 get_regex_pattern(spat) for spat in self._pats_nonexact]

401

402 def match(self, codes):

403 scodes = codes.safe_str

404 if scodes in self._pats_exact:

405 return True

406

407 return any(rpat.match(scodes) for rpat in self._pats_nonexact)

408

409

410def match_codes_any(patterns, codes):

411

412 pats_exact, pats_nonexact = classify_patterns(patterns)

413

414 if codes.safe_str in pats_exact:

415 return True

416

417 return any(match_codes(pattern, codes) for pattern in patterns)

418

419

420g_content_kinds = [

421 'undefined',

422 'waveform',

423 'station',

424 'channel',

425 'response',

426 'event',

427 'waveform_promise',

428 'empty']

429

430

431g_codes_classes = [

432 CodesX,

433 CodesNSLCE,

434 CodesNSL,

435 CodesNSLCE,

436 CodesNSLCE,

437 CodesX,

438 CodesNSLCE,

439 CodesX]

440

441g_codes_classes_ndot = {

442 0: CodesX,

443 2: CodesNSL,

444 3: CodesNSLCE,

445 4: CodesNSLCE}

446

447

448def to_codes_simple(kind_id, codes_safe_str):

449 return g_codes_classes[kind_id](safe_str=codes_safe_str)

450

451

452def to_codes(kind_id, obj):

453 return g_codes_classes[kind_id](obj)

454

455

456def to_codes_guess(s):

457 try:

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

459 except KeyError:

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

461

462

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

464class codes_patterns_list(list):

465 pass

466

467

468def codes_patterns_for_kind(kind_id, codes):

469 if isinstance(codes, codes_patterns_list):

470 return codes

471

472 if isinstance(codes, list):

473 lcodes = codes_patterns_list()

474 for sc in codes:

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

476

477 return lcodes

478

479 codes = to_codes(kind_id, codes)

480

481 lcodes = codes_patterns_list()

482 lcodes.append(codes)

483 if kind_id == STATION:

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

485

486 return lcodes

487

488

489g_content_kind_ids = (

490 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT,

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

492

493

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

495

496

497try:

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

499except Exception:

500 g_tmin_queries = g_tmin

501

502

503def to_kind(kind_id):

504 return g_content_kinds[kind_id]

505

506

507def to_kinds(kind_ids):

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

509

510

511def to_kind_id(kind):

512 return g_content_kinds.index(kind)

513

514

515def to_kind_ids(kinds):

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

517

518

519g_kind_mask_all = 0xff

520

521

522def to_kind_mask(kinds):

523 if kinds:

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

525 else:

526 return g_kind_mask_all

527

528

529def str_or_none(x):

530 if x is None:

531 return None

532 else:

533 return str(x)

534

535

536def float_or_none(x):

537 if x is None:

538 return None

539 else:

540 return float(x)

541

542

543def int_or_none(x):

544 if x is None:

545 return None

546 else:

547 return int(x)

548

549

550def int_or_g_tmin(x):

551 if x is None:

552 return g_tmin

553 else:

554 return int(x)

555

556

557def int_or_g_tmax(x):

558 if x is None:

559 return g_tmax

560 else:

561 return int(x)

562

563

564def tmin_or_none(x):

565 if x == g_tmin:

566 return None

567 else:

568 return x

569

570

571def tmax_or_none(x):

572 if x == g_tmax:

573 return None

574 else:

575 return x

576

577

578def time_or_none_to_str(x):

579 if x is None:

580 return '...'.ljust(17)

581 else:

582 return util.time_to_str(x)

583

584

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

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

587

588 if len(codes) > 20:

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

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

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

592 else:

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

594 if codes else '<none>'

595

596 return scodes

597

598

599g_offset_time_unit_inv = 1000000000

600g_offset_time_unit = 1.0 / g_offset_time_unit_inv

601

602

603def tsplit(t):

604 if t is None:

605 return None, 0

606

607 t = util.to_time_float(t)

608 if type(t) is float:

609 t = round(t, 5)

610 else:

611 t = round(t, 9)

612

613 seconds = num.floor(t)

614 offset = t - seconds

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

616

617

618def tjoin(seconds, offset):

619 if seconds is None:

620 return None

621

622 return util.to_time_float(seconds) \

623 + util.to_time_float(offset*g_offset_time_unit)

624

625

626tscale_min = 1

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

628tscale_logbase = 20

629

630tscale_edges = [tscale_min]

631while True:

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

633 if tscale_edges[-1] >= tscale_max:

634 break

635

636

637tscale_edges = num.array(tscale_edges)

638

639

640def tscale_to_kscale(tscale):

641

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

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

644 # ...

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

646

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

648

649

650@squirrel_content

651class Station(Location):

652 '''

653 A seismic station.

654 '''

655

656 codes = CodesNSL.T()

657

658 tmin = Timestamp.T(optional=True)

659 tmax = Timestamp.T(optional=True)

660

661 description = Unicode.T(optional=True)

662

663 def __init__(self, **kwargs):

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

665 Location.__init__(self, **kwargs)

666

667 @property

668 def time_span(self):

669 return (self.tmin, self.tmax)

670

671 def get_pyrocko_station(self):

672 '''

673 Get station as a classic Pyrocko station object.

674

675 :returns:

676 Converted station object.

677 :rtype:

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

679 '''

680 from pyrocko import model

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

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

696@squirrel_content

697class ChannelBase(Location):

698 codes = CodesNSLCE.T()

699

700 tmin = Timestamp.T(optional=True)

701 tmax = Timestamp.T(optional=True)

702

703 deltat = Float.T(optional=True)

704

705 def __init__(self, **kwargs):

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

707 Location.__init__(self, **kwargs)

708

709 @property

710 def time_span(self):

711 return (self.tmin, self.tmax)

712

713 def _get_sensor_codes(self):

714 return self.codes.replace(

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

716

717 def _get_sensor_args(self):

718 def getattr_rep(k):

719 if k == 'codes':

720 return self._get_sensor_codes()

721 else:

722 return getattr(self, k)

723

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

725

726 def _get_channel_args(self, component):

727 def getattr_rep(k):

728 if k == 'codes':

729 return self.codes.replace(

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

731 else:

732 return getattr(self, k)

733

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

735

736 def _get_pyrocko_station_args(self):

737 return (

738 self.codes.network,

739 self.codes.station,

740 self.codes.location,

741 self.lat,

742 self.lon,

743 self.elevation,

744 self.depth,

745 self.north_shift,

746 self.east_shift)

747

748

749class Channel(ChannelBase):

750 '''

751 A channel of a seismic station.

752 '''

753

754 dip = Float.T(optional=True)

755 azimuth = Float.T(optional=True)

756

757 def get_pyrocko_channel(self):

758 from pyrocko import model

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

760

761 def _get_pyrocko_channel_args(self):

762 return (

763 self.codes.channel,

764 self.azimuth,

765 self.dip)

766

767 @property

768 def orientation_enz(self):

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

770 return None

771

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

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

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

775 return mkvec(e, n, -d)

776

777

778def cut_intervals(channels):

779 channels = list(channels)

780 times = set()

781 for channel in channels:

782 if channel.tmin is not None:

783 times.add(channel.tmin)

784 if channel.tmax is not None:

785 times.add(channel.tmax)

786

787 times = sorted(times)

788

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

790 times[0:0] = [None]

791

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

793 times.append(None)

794

795 if len(times) <= 2:

796 return channels

797

798 channels_out = []

799 for channel in channels:

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

801 tmin = times[i]

802 tmax = times[i+1]

803 if ((channel.tmin is None or (

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

805 and (channel.tmax is None or (

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

807

808 channel_out = clone(channel)

809 channel_out.tmin = tmin

810 channel_out.tmax = tmax

811 channels_out.append(channel_out)

812

813 return channels_out

814

815

816class Sensor(ChannelBase):

817 '''

818 Representation of a channel group.

819 '''

820

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

822

823 @classmethod

824 def from_channels(cls, channels):

825 groups = defaultdict(list)

826 for channel in channels:

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

828

829 channels_cut = []

830 for group in groups.values():

831 channels_cut.extend(cut_intervals(group))

832

833 groups = defaultdict(list)

834 for channel in channels_cut:

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

836

837 return [

838 cls(channels=channels,

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

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

841

842 def channel_vectors(self):

843 return num.vstack(

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

845

846 def projected_channels(self, system, component_names):

847 return [

848 Channel(

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

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

851 **dict(zip(

852 ChannelBase.T.propnames,

853 self._get_channel_args(comp))))

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

855

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

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

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

859 return m

860

861 def projection_to(self, system, component_names):

862 return (

863 self.matrix_to(system),

864 self.channels,

865 self.projected_channels(system, component_names))

866

867 def projection_to_enz(self):

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

869

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

871 if azimuth is not None:

872 assert source is None

873 else:

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

875

876 return self.projection_to(num.array([

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

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

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

880

881 def project_to_enz(self, traces):

882 from pyrocko import trace

883

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

885

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

887

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

889 from pyrocko import trace

890

891 matrix, in_channels, out_channels = self.projection_to_trz(

892 source, azimuth=azimuth)

893

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

895

896

897observational_quantities = [

898 'acceleration', 'velocity', 'displacement', 'pressure',

899 'rotation_displacement', 'rotation_velocity', 'rotation_acceleration',

900 'temperature']

901

902

903technical_quantities = [

904 'voltage', 'counts']

905

906

907class QuantityType(StringChoice):

908 '''

909 Choice of observational or technical quantity.

910

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

912 '''

913 choices = observational_quantities + technical_quantities

914

915

916class ResponseStage(Object):

917 '''

918 Representation of a response stage.

919

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

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

922 downsampling.

923 '''

924 input_quantity = QuantityType.T(optional=True)

925 input_sample_rate = Float.T(optional=True)

926 output_quantity = QuantityType.T(optional=True)

927 output_sample_rate = Float.T(optional=True)

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

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

930

931 @property

932 def stage_type(self):

933 if self.input_quantity in observational_quantities \

934 and self.output_quantity in observational_quantities:

935 return 'conversion'

936

937 if self.input_quantity in observational_quantities \

938 and self.output_quantity == 'voltage':

939 return 'sensor'

940

941 elif self.input_quantity == 'voltage' \

942 and self.output_quantity == 'voltage':

943 return 'electronics'

944

945 elif self.input_quantity == 'voltage' \

946 and self.output_quantity == 'counts':

947 return 'digitizer'

948

949 elif self.decimation_factor is not None \

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

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

952 return 'decimation'

953

954 elif self.input_quantity in observational_quantities \

955 and self.output_quantity == 'counts':

956 return 'combination'

957

958 else:

959 return 'unknown'

960

961 @property

962 def decimation_factor(self):

963 irate = self.input_sample_rate

964 orate = self.output_sample_rate

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

966 and irate > orate and irate / orate > 1.0:

967

968 return irate / orate

969 else:

970 return None

971

972 @property

973 def summary_quantities(self):

974 return '%s -> %s' % (

975 self.input_quantity or '?',

976 self.output_quantity or '?')

977

978 @property

979 def summary_rates(self):

980 irate = self.input_sample_rate

981 orate = self.output_sample_rate

982 factor = self.decimation_factor

983

984 if irate and orate is None:

985 return '%g Hz' % irate

986

987 elif orate and irate is None:

988 return '%g Hz' % orate

989

990 elif irate and orate and irate == orate:

991 return '%g Hz' % irate

992

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

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

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

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

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

998 else:

999 return ''

1000

1001 @property

1002 def summary_elements(self):

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

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

1005

1006 @property

1007 def summary_log(self):

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

1009

1010 @property

1011 def summary_entries(self):

1012 return (

1013 self.__class__.__name__,

1014 self.stage_type,

1015 self.summary_log,

1016 self.summary_quantities,

1017 self.summary_rates,

1018 self.summary_elements)

1019

1020 @property

1021 def summary(self):

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

1023

1024 def get_effective(self):

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

1026

1027

1028D = 'displacement'

1029V = 'velocity'

1030A = 'acceleration'

1031

1032g_converters = {

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

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

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

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

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

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

1039

1040

1041def response_converters(input_quantity, output_quantity):

1042 if input_quantity is None or input_quantity == output_quantity:

1043 return []

1044

1045 if output_quantity is None:

1046 raise ConversionError('Unspecified target quantity.')

1047

1048 try:

1049 return [g_converters[input_quantity, output_quantity]]

1050

1051 except KeyError:

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

1053 input_quantity, output_quantity))

1054

1055

1056@squirrel_content

1057class Response(Object):

1058 '''

1059 The instrument response of a seismic station channel.

1060 '''

1061

1062 codes = CodesNSLCE.T()

1063 tmin = Timestamp.T(optional=True)

1064 tmax = Timestamp.T(optional=True)

1065

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

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

1068

1069 deltat = Float.T(optional=True)

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

1071

1072 def __init__(self, **kwargs):

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

1074 Object.__init__(self, **kwargs)

1075 self._effective_responses_cache = {}

1076

1077 @property

1078 def time_span(self):

1079 return (self.tmin, self.tmax)

1080

1081 @property

1082 def nstages(self):

1083 return len(self.stages)

1084

1085 @property

1086 def input_quantity(self):

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

1088

1089 @property

1090 def output_quantity(self):

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

1092

1093 @property

1094 def output_sample_rate(self):

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

1096

1097 @property

1098 def summary_stages(self):

1099 def grouped(xs):

1100 xs = list(xs)

1101 g = []

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

1103 g.append(xs[i])

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

1105 yield g

1106 g = []

1107

1108 if g:

1109 yield g

1110

1111 return '+'.join(

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

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

1114

1115 @property

1116 def summary_quantities(self):

1117 orate = self.output_sample_rate

1118 return '%s -> %s%s' % (

1119 self.input_quantity or '?',

1120 self.output_quantity or '?',

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

1122

1123 @property

1124 def summary_log(self):

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

1126

1127 @property

1128 def summary_entries(self):

1129 return (

1130 self.__class__.__name__,

1131 str(self.codes),

1132 self.str_time_span,

1133 self.summary_log,

1134 self.summary_quantities,

1135 self.summary_stages)

1136

1137 @property

1138 def summary(self):

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

1140

1141 def get_effective(

1142 self,

1143 input_quantity=None,

1144 stages=(None, None),

1145 mode='complete',

1146 gain_frequency=None):

1147

1148 assert mode in ('complete', 'sensor')

1149 assert not (mode == 'sensor' and gain_frequency is None)

1150

1151 cache_key = (input_quantity, stages)

1152 if cache_key in self._effective_responses_cache:

1153 return self._effective_responses_cache[cache_key]

1154

1155 try:

1156 elements = response_converters(input_quantity, self.input_quantity)

1157 except ConversionError as e:

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

1159

1160 for istage, stage in enumerate(self.stages):

1161 if (stages[0] is None or stages[0] <= istage) \

1162 and (stages[1] is None or istage < stages[1]):

1163

1164 is_sensor = istage == 0 or stage.stage_type == 'sensor'

1165

1166 if mode == 'complete' or (mode == 'sensor' and is_sensor):

1167 elements.append(stage.get_effective())

1168 else:

1169 resp = stage.get_effective()

1170 gain = Gain(constant=resp.evaluate1(gain_frequency))

1171 # maybe check derivative if response is flat enough?

1172 # or check frequency band?

1173 elements.append(gain)

1174

1175 if input_quantity is None \

1176 or input_quantity == self.input_quantity:

1177 checkpoints = self.checkpoints

1178 else:

1179 checkpoints = []

1180

1181 resp = MultiplyResponse(

1182 responses=simplify_responses(elements),

1183 checkpoints=checkpoints)

1184

1185 self._effective_responses_cache[cache_key] = resp

1186 return resp

1187

1188

1189@squirrel_content

1190class Event(Object):

1191 '''

1192 A seismic event.

1193 '''

1194

1195 name = String.T(optional=True)

1196 time = Timestamp.T()

1197 duration = Float.T(optional=True)

1198

1199 lat = Float.T()

1200 lon = Float.T()

1201 depth = Float.T(optional=True)

1202

1203 magnitude = Float.T(optional=True)

1204

1205 def get_pyrocko_event(self):

1206 from pyrocko import model

1207 model.Event(

1208 name=self.name,

1209 time=self.time,

1210 lat=self.lat,

1211 lon=self.lon,

1212 depth=self.depth,

1213 magnitude=self.magnitude,

1214 duration=self.duration)

1215

1216 @property

1217 def time_span(self):

1218 return (self.time, self.time)

1219

1220

1221def ehash(s):

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

1223

1224

1225def random_name(n=8):

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

1227

1228

1229@squirrel_content

1230class WaveformPromise(Object):

1231 '''

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

1233

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

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

1236 calls to

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

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

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

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

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

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

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

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

1245 queries for waveforms missing at the remote site.

1246 '''

1247

1248 codes = CodesNSLCE.T()

1249 tmin = Timestamp.T()

1250 tmax = Timestamp.T()

1251

1252 deltat = Float.T(optional=True)

1253

1254 source_hash = String.T()

1255

1256 def __init__(self, **kwargs):

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

1258 Object.__init__(self, **kwargs)

1259

1260 @property

1261 def time_span(self):

1262 return (self.tmin, self.tmax)

1263

1264

1265class InvalidWaveform(Exception):

1266 pass

1267

1268

1269class WaveformOrder(Object):

1270 '''

1271 Waveform request information.

1272 '''

1273

1274 source_id = String.T()

1275 codes = CodesNSLCE.T()

1276 deltat = Float.T()

1277 tmin = Timestamp.T()

1278 tmax = Timestamp.T()

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

1280 time_created = Timestamp.T()

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

1282

1283 @property

1284 def client(self):

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

1286

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

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

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

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

1291

1292 def validate(self, tr):

1293 if tr.ydata.size == 0:

1294 raise InvalidWaveform(

1295 'waveform with zero data samples')

1296

1297 if tr.deltat != self.deltat:

1298 raise InvalidWaveform(

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

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

1301 tr.deltat, self.deltat))

1302

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

1304 raise InvalidWaveform('waveform has NaN values')

1305

1306 def estimate_nsamples(self):

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

1308

1309 def is_near_real_time(self):

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

1311

1312

1313def order_summary(orders):

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

1315 if len(codes_list) > 3:

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

1317 len(orders),

1318 util.plural_s(orders),

1319 str(codes_list[0]),

1320 str(codes_list[1]))

1321

1322 else:

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

1324 len(orders),

1325 util.plural_s(orders),

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

1327

1328

1329class Nut(Object):

1330 '''

1331 Index entry referencing an elementary piece of content.

1332

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

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

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

1336 '''

1337

1338 file_path = String.T(optional=True)

1339 file_format = String.T(optional=True)

1340 file_mtime = Timestamp.T(optional=True)

1341 file_size = Int.T(optional=True)

1342

1343 file_segment = Int.T(optional=True)

1344 file_element = Int.T(optional=True)

1345

1346 kind_id = Int.T()

1347 codes = Codes.T()

1348

1349 tmin_seconds = Int.T(default=0)

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

1351 tmax_seconds = Int.T(default=0)

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

1353

1354 deltat = Float.T(default=0.0)

1355

1356 content = Any.T(optional=True)

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

1358

1359 content_in_db = False

1360

1361 def __init__(

1362 self,

1363 file_path=None,

1364 file_format=None,

1365 file_mtime=None,

1366 file_size=None,

1367 file_segment=None,

1368 file_element=None,

1369 kind_id=0,

1370 codes=CodesX(''),

1371 tmin_seconds=None,

1372 tmin_offset=0,

1373 tmax_seconds=None,

1374 tmax_offset=0,

1375 deltat=None,

1376 content=None,

1377 raw_content=None,

1378 tmin=None,

1379 tmax=None,

1380 values_nocheck=None):

1381

1382 if values_nocheck is not None:

1383 (self.file_path, self.file_format, self.file_mtime,

1384 self.file_size,

1385 self.file_segment, self.file_element,

1386 self.kind_id, codes_safe_str,

1387 self.tmin_seconds, self.tmin_offset,

1388 self.tmax_seconds, self.tmax_offset,

1389 self.deltat) = values_nocheck

1390

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

1392 self.deltat = self.deltat or None

1393 self.raw_content = {}

1394 self.content = None

1395 else:

1396 if tmin is not None:

1397 tmin_seconds, tmin_offset = tsplit(tmin)

1398

1399 if tmax is not None:

1400 tmax_seconds, tmax_offset = tsplit(tmax)

1401

1402 self.kind_id = int(kind_id)

1403 self.codes = codes

1404 self.tmin_seconds = int_or_g_tmin(tmin_seconds)

1405 self.tmin_offset = int(tmin_offset)

1406 self.tmax_seconds = int_or_g_tmax(tmax_seconds)

1407 self.tmax_offset = int(tmax_offset)

1408 self.deltat = float_or_none(deltat)

1409 self.file_path = str_or_none(file_path)

1410 self.file_segment = int_or_none(file_segment)

1411 self.file_element = int_or_none(file_element)

1412 self.file_format = str_or_none(file_format)

1413 self.file_mtime = float_or_none(file_mtime)

1414 self.file_size = int_or_none(file_size)

1415 self.content = content

1416 if raw_content is None:

1417 self.raw_content = {}

1418 else:

1419 self.raw_content = raw_content

1420

1421 Object.__init__(self, init_props=False)

1422

1423 def __eq__(self, other):

1424 return (isinstance(other, Nut) and

1425 self.equality_values == other.equality_values)

1426

1427 def hash(self):

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

1429

1430 def __ne__(self, other):

1431 return not (self == other)

1432

1433 def get_io_backend(self):

1434 from . import io

1435 return io.get_backend(self.file_format)

1436

1437 def file_modified(self):

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

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

1440

1441 @property

1442 def dkey(self):

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

1444

1445 @property

1446 def key(self):

1447 return (

1448 self.file_path,

1449 self.file_segment,

1450 self.file_element,

1451 self.file_mtime)

1452

1453 @property

1454 def equality_values(self):

1455 return (

1456 self.file_segment, self.file_element,

1457 self.kind_id, self.codes,

1458 self.tmin_seconds, self.tmin_offset,

1459 self.tmax_seconds, self.tmax_offset, self.deltat)

1460

1461 def diff(self, other):

1462 names = [

1463 'file_segment', 'file_element', 'kind_id', 'codes',

1464 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',

1465 'deltat']

1466

1467 d = []

1468 for name, a, b in zip(

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

1470

1471 if a != b:

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

1473

1474 return d

1475

1476 @property

1477 def tmin(self):

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

1479

1480 @tmin.setter

1481 def tmin(self, t):

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

1483

1484 @property

1485 def tmax(self):

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

1487

1488 @tmax.setter

1489 def tmax(self, t):

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

1491

1492 @property

1493 def kscale(self):

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

1495 return 0

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

1497

1498 @property

1499 def waveform_kwargs(self):

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

1501

1502 return dict(

1503 network=network,

1504 station=station,

1505 location=location,

1506 channel=channel,

1507 extra=extra,

1508 tmin=self.tmin,

1509 tmax=self.tmax,

1510 deltat=self.deltat)

1511

1512 @property

1513 def waveform_promise_kwargs(self):

1514 return dict(

1515 codes=self.codes,

1516 tmin=self.tmin,

1517 tmax=self.tmax,

1518 deltat=self.deltat)

1519

1520 @property

1521 def station_kwargs(self):

1522 network, station, location = self.codes

1523 return dict(

1524 codes=self.codes,

1525 tmin=tmin_or_none(self.tmin),

1526 tmax=tmax_or_none(self.tmax))

1527

1528 @property

1529 def channel_kwargs(self):

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

1531 return dict(

1532 codes=self.codes,

1533 tmin=tmin_or_none(self.tmin),

1534 tmax=tmax_or_none(self.tmax),

1535 deltat=self.deltat)

1536

1537 @property

1538 def response_kwargs(self):

1539 return dict(

1540 codes=self.codes,

1541 tmin=tmin_or_none(self.tmin),

1542 tmax=tmax_or_none(self.tmax),

1543 deltat=self.deltat)

1544

1545 @property

1546 def event_kwargs(self):

1547 return dict(

1548 name=self.codes,

1549 time=self.tmin,

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

1551

1552 @property

1553 def trace_kwargs(self):

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

1555

1556 return dict(

1557 network=network,

1558 station=station,

1559 location=location,

1560 channel=channel,

1561 extra=extra,

1562 tmin=self.tmin,

1563 tmax=self.tmax-self.deltat,

1564 deltat=self.deltat)

1565

1566 @property

1567 def dummy_trace(self):

1568 return DummyTrace(self)

1569

1570 @property

1571 def summary_entries(self):

1572 if self.tmin == self.tmax:

1573 ts = util.time_to_str(self.tmin)

1574 else:

1575 ts = '%s - %s' % (

1576 util.time_to_str(self.tmin),

1577 util.time_to_str(self.tmax))

1578

1579 return (

1580 self.__class__.__name__,

1581 to_kind(self.kind_id),

1582 str(self.codes),

1583 ts)

1584

1585 @property

1586 def summary(self):

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

1588

1589

1590def make_waveform_nut(**kwargs):

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

1592

1593

1594def make_waveform_promise_nut(**kwargs):

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

1596

1597

1598def make_station_nut(**kwargs):

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

1600

1601

1602def make_channel_nut(**kwargs):

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

1604

1605

1606def make_response_nut(**kwargs):

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

1608

1609

1610def make_event_nut(**kwargs):

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

1612

1613

1614def group_channels(nuts):

1615 by_ansl = {}

1616 for nut in nuts:

1617 if nut.kind_id != CHANNEL:

1618 continue

1619

1620 ansl = nut.codes[:4]

1621

1622 if ansl not in by_ansl:

1623 by_ansl[ansl] = {}

1624

1625 group = by_ansl[ansl]

1626

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

1628

1629 if k not in group:

1630 group[k] = set()

1631

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

1633

1634 return by_ansl

1635

1636

1637class DummyTrace(object):

1638

1639 def __init__(self, nut):

1640 self.nut = nut

1641 self.codes = nut.codes

1642 self.meta = {}

1643

1644 @property

1645 def tmin(self):

1646 return self.nut.tmin

1647

1648 @property

1649 def tmax(self):

1650 return self.nut.tmax

1651

1652 @property

1653 def deltat(self):

1654 return self.nut.deltat

1655

1656 @property

1657 def nslc_id(self):

1658 return self.codes.nslc

1659

1660 @property

1661 def network(self):

1662 return self.codes.network

1663

1664 @property

1665 def station(self):

1666 return self.codes.station

1667

1668 @property

1669 def location(self):

1670 return self.codes.location

1671

1672 @property

1673 def channel(self):

1674 return self.codes.channel

1675

1676 @property

1677 def extra(self):

1678 return self.codes.extra

1679

1680 def overlaps(self, tmin, tmax):

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

1682

1683

1684def duration_to_str(t):

1685 if t > 24*3600:

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

1687 elif t > 3600:

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

1689 elif t > 60:

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

1691 else:

1692 return '%gs' % t

1693

1694

1695class Coverage(Object):

1696 '''

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

1698 '''

1699 kind_id = Int.T(

1700 help='Content type.')

1701 pattern = Codes.T(

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

1703 'match.')

1704 codes = Codes.T(

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

1706 deltat = Float.T(

1707 help='Sampling interval [s]',

1708 optional=True)

1709 tmin = Timestamp.T(

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

1711 optional=True)

1712 tmax = Timestamp.T(

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

1714 optional=True)

1715 changes = List.T(

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

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

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

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

1720 'value duplicate or redundant data coverage.')

1721

1722 @classmethod

1723 def from_values(cls, args):

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

1725 return cls(

1726 kind_id=kind_id,

1727 pattern=pattern,

1728 codes=codes,

1729 deltat=deltat,

1730 tmin=tmin,

1731 tmax=tmax,

1732 changes=changes)

1733

1734 @property

1735 def summary_entries(self):

1736 ts = '%s - %s' % (

1737 util.time_to_str(self.tmin),

1738 util.time_to_str(self.tmax))

1739

1740 srate = self.sample_rate

1741

1742 total = self.total

1743

1744 return (

1745 self.__class__.__name__,

1746 to_kind(self.kind_id),

1747 str(self.codes),

1748 ts,

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

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

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

1752

1753 @property

1754 def summary(self):

1755 return util.fmt_summary(

1756 self.summary_entries,

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

1758

1759 @property

1760 def sample_rate(self):

1761 if self.deltat is None:

1762 return None

1763 elif self.deltat == 0.0:

1764 return 0.0

1765 else:

1766 return 1.0 / self.deltat

1767

1768 @property

1769 def labels(self):

1770 srate = self.sample_rate

1771 return (

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

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

1774

1775 @property

1776 def total(self):

1777 total_t = None

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

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

1780

1781 return total_t

1782

1783 def iter_spans(self):

1784 last = None

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

1786 if last is not None:

1787 last_t, last_count = last

1788 if last_count > 0:

1789 yield last_t, t, last_count

1790

1791 last = (t, count)

1792

1793 def iter_uncovered_by(self, other):

1794 a = self

1795 b = other

1796 ia = ib = -1

1797 ca = cb = 0

1798 last = None

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

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

1801 ia += 1

1802 t, ca = a.changes[ia]

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

1804 ib += 1

1805 t, cb = b.changes[ib]

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

1807 ia += 1

1808 t, ca = a.changes[ia]

1809 else:

1810 ib += 1

1811 t, cb = b.changes[ib]

1812

1813 if last is not None:

1814 tl, cal, cbl = last

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

1816 yield tl, t, ia, ib

1817

1818 last = (t, ca, cb)

1819

1820 def iter_uncovered_by_combined(self, other):

1821 ib_last = None

1822 group = None

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

1824 if ib_last is None or ib != ib_last:

1825 if group:

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

1827

1828 group = []

1829

1830 group.append((tmin, tmax))

1831 ib_last = ib

1832

1833 if group:

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

1835

1836

1837__all__ = [

1838 'UNDEFINED',

1839 'WAVEFORM',

1840 'STATION',

1841 'CHANNEL',

1842 'RESPONSE',

1843 'EVENT',

1844 'WAVEFORM_PROMISE',

1845 'EMPTY',

1846 'to_codes',

1847 'to_codes_guess',

1848 'to_codes_simple',

1849 'to_kind',

1850 'to_kinds',

1851 'to_kind_id',

1852 'to_kind_ids',

1853 'match_codes',

1854 'codes_patterns_for_kind',

1855 'CodesMatcher',

1856 'CodesError',

1857 'Codes',

1858 'CodesNSLCE',

1859 'CodesNSL',

1860 'CodesX',

1861 'Station',

1862 'Channel',

1863 'Sensor',

1864 'Response',

1865 'Nut',

1866 'Coverage',

1867 'WaveformPromise',

1868]