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

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

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 match_codes(pattern, codes):

367 spattern = pattern.safe_str

368 scodes = codes.safe_str

369 if spattern not in g_codes_patterns:

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

371 g_codes_patterns[spattern] = rpattern

372

373 rpattern = g_codes_patterns[spattern]

374 return bool(rpattern.match(scodes))

375

376

377g_content_kinds = [

378 'undefined',

379 'waveform',

380 'station',

381 'channel',

382 'response',

383 'event',

384 'waveform_promise',

385 'empty']

386

387

388g_codes_classes = [

389 CodesX,

390 CodesNSLCE,

391 CodesNSL,

392 CodesNSLCE,

393 CodesNSLCE,

394 CodesX,

395 CodesNSLCE,

396 CodesX]

397

398g_codes_classes_ndot = {

399 0: CodesX,

400 2: CodesNSL,

401 3: CodesNSLCE,

402 4: CodesNSLCE}

403

404

405def to_codes_simple(kind_id, codes_safe_str):

406 return g_codes_classes[kind_id](safe_str=codes_safe_str)

407

408

409def to_codes(kind_id, obj):

410 return g_codes_classes[kind_id](obj)

411

412

413def to_codes_guess(s):

414 try:

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

416 except KeyError:

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

418

419

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

421class codes_patterns_list(list):

422 pass

423

424

425def codes_patterns_for_kind(kind_id, codes):

426 if isinstance(codes, codes_patterns_list):

427 return codes

428

429 if isinstance(codes, list):

430 lcodes = codes_patterns_list()

431 for sc in codes:

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

433

434 return lcodes

435

436 codes = to_codes(kind_id, codes)

437

438 lcodes = codes_patterns_list()

439 lcodes.append(codes)

440 if kind_id == STATION:

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

442

443 return lcodes

444

445

446g_content_kind_ids = (

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

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

449

450

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

452

453

454try:

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

456except Exception:

457 g_tmin_queries = g_tmin

458

459

460def to_kind(kind_id):

461 return g_content_kinds[kind_id]

462

463

464def to_kinds(kind_ids):

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

466

467

468def to_kind_id(kind):

469 return g_content_kinds.index(kind)

470

471

472def to_kind_ids(kinds):

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

474

475

476g_kind_mask_all = 0xff

477

478

479def to_kind_mask(kinds):

480 if kinds:

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

482 else:

483 return g_kind_mask_all

484

485

486def str_or_none(x):

487 if x is None:

488 return None

489 else:

490 return str(x)

491

492

493def float_or_none(x):

494 if x is None:

495 return None

496 else:

497 return float(x)

498

499

500def int_or_none(x):

501 if x is None:

502 return None

503 else:

504 return int(x)

505

506

507def int_or_g_tmin(x):

508 if x is None:

509 return g_tmin

510 else:

511 return int(x)

512

513

514def int_or_g_tmax(x):

515 if x is None:

516 return g_tmax

517 else:

518 return int(x)

519

520

521def tmin_or_none(x):

522 if x == g_tmin:

523 return None

524 else:

525 return x

526

527

528def tmax_or_none(x):

529 if x == g_tmax:

530 return None

531 else:

532 return x

533

534

535def time_or_none_to_str(x):

536 if x is None:

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

538 else:

539 return util.time_to_str(x)

540

541

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

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

544

545 if len(codes) > 20:

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

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

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

549 else:

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

551 if codes else '<none>'

552

553 return scodes

554

555

556g_offset_time_unit_inv = 1000000000

557g_offset_time_unit = 1.0 / g_offset_time_unit_inv

558

559

560def tsplit(t):

561 if t is None:

562 return None, 0

563

564 t = util.to_time_float(t)

565 if type(t) is float:

566 t = round(t, 5)

567 else:

568 t = round(t, 9)

569

570 seconds = num.floor(t)

571 offset = t - seconds

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

573

574

575def tjoin(seconds, offset):

576 if seconds is None:

577 return None

578

579 return util.to_time_float(seconds) \

580 + util.to_time_float(offset*g_offset_time_unit)

581

582

583tscale_min = 1

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

585tscale_logbase = 20

586

587tscale_edges = [tscale_min]

588while True:

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

590 if tscale_edges[-1] >= tscale_max:

591 break

592

593

594tscale_edges = num.array(tscale_edges)

595

596

597def tscale_to_kscale(tscale):

598

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

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

601 # ...

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

603

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

605

606

607@squirrel_content

608class Station(Location):

609 '''

610 A seismic station.

611 '''

612

613 codes = CodesNSL.T()

614

615 tmin = Timestamp.T(optional=True)

616 tmax = Timestamp.T(optional=True)

617

618 description = Unicode.T(optional=True)

619

620 def __init__(self, **kwargs):

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

622 Location.__init__(self, **kwargs)

623

624 @property

625 def time_span(self):

626 return (self.tmin, self.tmax)

627

628 def get_pyrocko_station(self):

629 '''

630 Get station as a classic Pyrocko station object.

631

632 :returns:

633 Converted station object.

634 :rtype:

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

636 '''

637 from pyrocko import model

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

639

640 def _get_pyrocko_station_args(self):

641 return (

642 self.codes.network,

643 self.codes.station,

644 self.codes.location,

645 self.lat,

646 self.lon,

647 self.elevation,

648 self.depth,

649 self.north_shift,

650 self.east_shift)

651

652

653@squirrel_content

654class ChannelBase(Location):

655 codes = CodesNSLCE.T()

656

657 tmin = Timestamp.T(optional=True)

658 tmax = Timestamp.T(optional=True)

659

660 deltat = Float.T(optional=True)

661

662 def __init__(self, **kwargs):

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

664 Location.__init__(self, **kwargs)

665

666 @property

667 def time_span(self):

668 return (self.tmin, self.tmax)

669

670 def _get_sensor_codes(self):

671 return self.codes.replace(

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

673

674 def _get_sensor_args(self):

675 def getattr_rep(k):

676 if k == 'codes':

677 return self._get_sensor_codes()

678 else:

679 return getattr(self, k)

680

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

682

683 def _get_channel_args(self, component):

684 def getattr_rep(k):

685 if k == 'codes':

686 return self.codes.replace(

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

688 else:

689 return getattr(self, k)

690

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

692

693 def _get_pyrocko_station_args(self):

694 return (

695 self.codes.network,

696 self.codes.station,

697 self.codes.location,

698 self.lat,

699 self.lon,

700 self.elevation,

701 self.depth,

702 self.north_shift,

703 self.east_shift)

704

705

706class Channel(ChannelBase):

707 '''

708 A channel of a seismic station.

709 '''

710

711 dip = Float.T(optional=True)

712 azimuth = Float.T(optional=True)

713

714 def get_pyrocko_channel(self):

715 from pyrocko import model

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

717

718 def _get_pyrocko_channel_args(self):

719 return (

720 self.codes.channel,

721 self.azimuth,

722 self.dip)

723

724 @property

725 def orientation_enz(self):

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

727 return None

728

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

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

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

732 return mkvec(e, n, -d)

733

734

735def cut_intervals(channels):

736 channels = list(channels)

737 times = set()

738 for channel in channels:

739 if channel.tmin is not None:

740 times.add(channel.tmin)

741 if channel.tmax is not None:

742 times.add(channel.tmax)

743

744 times = sorted(times)

745

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

747 times[0:0] = [None]

748

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

750 times.append(None)

751

752 if len(times) <= 2:

753 return channels

754

755 channels_out = []

756 for channel in channels:

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

758 tmin = times[i]

759 tmax = times[i+1]

760 if ((channel.tmin is None or (

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

762 and (channel.tmax is None or (

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

764

765 channel_out = clone(channel)

766 channel_out.tmin = tmin

767 channel_out.tmax = tmax

768 channels_out.append(channel_out)

769

770 return channels_out

771

772

773class Sensor(ChannelBase):

774 '''

775 Representation of a channel group.

776 '''

777

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

779

780 @classmethod

781 def from_channels(cls, channels):

782 groups = defaultdict(list)

783 for channel in channels:

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

785

786 channels_cut = []

787 for group in groups.values():

788 channels_cut.extend(cut_intervals(group))

789

790 groups = defaultdict(list)

791 for channel in channels_cut:

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

793

794 return [

795 cls(channels=channels,

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

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

798

799 def channel_vectors(self):

800 return num.vstack(

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

802

803 def projected_channels(self, system, component_names):

804 return [

805 Channel(

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

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

808 **dict(zip(

809 ChannelBase.T.propnames,

810 self._get_channel_args(comp))))

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

812

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

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

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

816 return m

817

818 def projection_to(self, system, component_names):

819 return (

820 self.matrix_to(system),

821 self.channels,

822 self.projected_channels(system, component_names))

823

824 def projection_to_enz(self):

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

826

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

828 if azimuth is not None:

829 assert source is None

830 else:

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

832

833 return self.projection_to(num.array([

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

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

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

837

838 def project_to_enz(self, traces):

839 from pyrocko import trace

840

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

842

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

844

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

846 from pyrocko import trace

847

848 matrix, in_channels, out_channels = self.projection_to_trz(

849 source, azimuth=azimuth)

850

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

852

853

854observational_quantities = [

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

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

857 'temperature']

858

859

860technical_quantities = [

861 'voltage', 'counts']

862

863

864class QuantityType(StringChoice):

865 '''

866 Choice of observational or technical quantity.

867

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

869 '''

870 choices = observational_quantities + technical_quantities

871

872

873class ResponseStage(Object):

874 '''

875 Representation of a response stage.

876

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

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

879 downsampling.

880 '''

881 input_quantity = QuantityType.T(optional=True)

882 input_sample_rate = Float.T(optional=True)

883 output_quantity = QuantityType.T(optional=True)

884 output_sample_rate = Float.T(optional=True)

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

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

887

888 @property

889 def stage_type(self):

890 if self.input_quantity in observational_quantities \

891 and self.output_quantity in observational_quantities:

892 return 'conversion'

893

894 if self.input_quantity in observational_quantities \

895 and self.output_quantity == 'voltage':

896 return 'sensor'

897

898 elif self.input_quantity == 'voltage' \

899 and self.output_quantity == 'voltage':

900 return 'electronics'

901

902 elif self.input_quantity == 'voltage' \

903 and self.output_quantity == 'counts':

904 return 'digitizer'

905

906 elif self.decimation_factor is not None \

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

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

909 return 'decimation'

910

911 elif self.input_quantity in observational_quantities \

912 and self.output_quantity == 'counts':

913 return 'combination'

914

915 else:

916 return 'unknown'

917

918 @property

919 def decimation_factor(self):

920 irate = self.input_sample_rate

921 orate = self.output_sample_rate

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

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

924

925 return irate / orate

926 else:

927 return None

928

929 @property

930 def summary_quantities(self):

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

932 self.input_quantity or '?',

933 self.output_quantity or '?')

934

935 @property

936 def summary_rates(self):

937 irate = self.input_sample_rate

938 orate = self.output_sample_rate

939 factor = self.decimation_factor

940

941 if irate and orate is None:

942 return '%g Hz' % irate

943

944 elif orate and irate is None:

945 return '%g Hz' % orate

946

947 elif irate and orate and irate == orate:

948 return '%g Hz' % irate

949

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

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

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

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

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

955 else:

956 return ''

957

958 @property

959 def summary_elements(self):

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

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

962

963 @property

964 def summary_log(self):

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

966

967 @property

968 def summary_entries(self):

969 return (

970 self.__class__.__name__,

971 self.stage_type,

972 self.summary_log,

973 self.summary_quantities,

974 self.summary_rates,

975 self.summary_elements)

976

977 @property

978 def summary(self):

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

980

981 def get_effective(self):

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

983

984

985D = 'displacement'

986V = 'velocity'

987A = 'acceleration'

988

989g_converters = {

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

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

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

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

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

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

996

997

998def response_converters(input_quantity, output_quantity):

999 if input_quantity is None or input_quantity == output_quantity:

1000 return []

1001

1002 if output_quantity is None:

1003 raise ConversionError('Unspecified target quantity.')

1004

1005 try:

1006 return [g_converters[input_quantity, output_quantity]]

1007

1008 except KeyError:

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

1010 input_quantity, output_quantity))

1011

1012

1013@squirrel_content

1014class Response(Object):

1015 '''

1016 The instrument response of a seismic station channel.

1017 '''

1018

1019 codes = CodesNSLCE.T()

1020 tmin = Timestamp.T(optional=True)

1021 tmax = Timestamp.T(optional=True)

1022

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

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

1025

1026 deltat = Float.T(optional=True)

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

1028

1029 def __init__(self, **kwargs):

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

1031 Object.__init__(self, **kwargs)

1032 self._effective_responses_cache = {}

1033

1034 @property

1035 def time_span(self):

1036 return (self.tmin, self.tmax)

1037

1038 @property

1039 def nstages(self):

1040 return len(self.stages)

1041

1042 @property

1043 def input_quantity(self):

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

1045

1046 @property

1047 def output_quantity(self):

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

1049

1050 @property

1051 def output_sample_rate(self):

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

1053

1054 @property

1055 def summary_stages(self):

1056 def grouped(xs):

1057 xs = list(xs)

1058 g = []

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

1060 g.append(xs[i])

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

1062 yield g

1063 g = []

1064

1065 if g:

1066 yield g

1067

1068 return '+'.join(

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

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

1071

1072 @property

1073 def summary_quantities(self):

1074 orate = self.output_sample_rate

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

1076 self.input_quantity or '?',

1077 self.output_quantity or '?',

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

1079

1080 @property

1081 def summary_log(self):

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

1083

1084 @property

1085 def summary_entries(self):

1086 return (

1087 self.__class__.__name__,

1088 str(self.codes),

1089 self.str_time_span,

1090 self.summary_log,

1091 self.summary_quantities,

1092 self.summary_stages)

1093

1094 @property

1095 def summary(self):

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

1097

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

1099 cache_key = (input_quantity, stages)

1100 if cache_key in self._effective_responses_cache:

1101 return self._effective_responses_cache[cache_key]

1102

1103 try:

1104 elements = response_converters(input_quantity, self.input_quantity)

1105 except ConversionError as e:

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

1107

1108 elements.extend(

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

1110

1111 if input_quantity is None \

1112 or input_quantity == self.input_quantity:

1113 checkpoints = self.checkpoints

1114 else:

1115 checkpoints = []

1116

1117 resp = MultiplyResponse(

1118 responses=simplify_responses(elements),

1119 checkpoints=checkpoints)

1120

1121 self._effective_responses_cache[cache_key] = resp

1122 return resp

1123

1124

1125@squirrel_content

1126class Event(Object):

1127 '''

1128 A seismic event.

1129 '''

1130

1131 name = String.T(optional=True)

1132 time = Timestamp.T()

1133 duration = Float.T(optional=True)

1134

1135 lat = Float.T()

1136 lon = Float.T()

1137 depth = Float.T(optional=True)

1138

1139 magnitude = Float.T(optional=True)

1140

1141 def get_pyrocko_event(self):

1142 from pyrocko import model

1143 model.Event(

1144 name=self.name,

1145 time=self.time,

1146 lat=self.lat,

1147 lon=self.lon,

1148 depth=self.depth,

1149 magnitude=self.magnitude,

1150 duration=self.duration)

1151

1152 @property

1153 def time_span(self):

1154 return (self.time, self.time)

1155

1156

1157def ehash(s):

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

1159

1160

1161def random_name(n=8):

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

1163

1164

1165@squirrel_content

1166class WaveformPromise(Object):

1167 '''

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

1169

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

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

1172 calls to

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

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

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

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

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

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

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

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

1181 queries for waveforms missing at the remote site.

1182 '''

1183

1184 codes = CodesNSLCE.T()

1185 tmin = Timestamp.T()

1186 tmax = Timestamp.T()

1187

1188 deltat = Float.T(optional=True)

1189

1190 source_hash = String.T()

1191

1192 def __init__(self, **kwargs):

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

1194 Object.__init__(self, **kwargs)

1195

1196 @property

1197 def time_span(self):

1198 return (self.tmin, self.tmax)

1199

1200

1201class InvalidWaveform(Exception):

1202 pass

1203

1204

1205class WaveformOrder(Object):

1206 '''

1207 Waveform request information.

1208 '''

1209

1210 source_id = String.T()

1211 codes = CodesNSLCE.T()

1212 deltat = Float.T()

1213 tmin = Timestamp.T()

1214 tmax = Timestamp.T()

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

1216 time_created = Timestamp.T()

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

1218

1219 @property

1220 def client(self):

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

1222

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

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

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

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

1227

1228 def validate(self, tr):

1229 if tr.ydata.size == 0:

1230 raise InvalidWaveform(

1231 'waveform with zero data samples')

1232

1233 if tr.deltat != self.deltat:

1234 raise InvalidWaveform(

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

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

1237 tr.deltat, self.deltat))

1238

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

1240 raise InvalidWaveform('waveform has NaN values')

1241

1242 def estimate_nsamples(self):

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

1244

1245 def is_near_real_time(self):

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

1247

1248

1249def order_summary(orders):

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

1251 if len(codes_list) > 3:

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

1253 len(orders),

1254 util.plural_s(orders),

1255 str(codes_list[0]),

1256 str(codes_list[1]))

1257

1258 else:

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

1260 len(orders),

1261 util.plural_s(orders),

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

1263

1264

1265class Nut(Object):

1266 '''

1267 Index entry referencing an elementary piece of content.

1268

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

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

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

1272 '''

1273

1274 file_path = String.T(optional=True)

1275 file_format = String.T(optional=True)

1276 file_mtime = Timestamp.T(optional=True)

1277 file_size = Int.T(optional=True)

1278

1279 file_segment = Int.T(optional=True)

1280 file_element = Int.T(optional=True)

1281

1282 kind_id = Int.T()

1283 codes = Codes.T()

1284

1285 tmin_seconds = Int.T(default=0)

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

1287 tmax_seconds = Int.T(default=0)

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

1289

1290 deltat = Float.T(default=0.0)

1291

1292 content = Any.T(optional=True)

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

1294

1295 content_in_db = False

1296

1297 def __init__(

1298 self,

1299 file_path=None,

1300 file_format=None,

1301 file_mtime=None,

1302 file_size=None,

1303 file_segment=None,

1304 file_element=None,

1305 kind_id=0,

1306 codes=CodesX(''),

1307 tmin_seconds=None,

1308 tmin_offset=0,

1309 tmax_seconds=None,

1310 tmax_offset=0,

1311 deltat=None,

1312 content=None,

1313 raw_content=None,

1314 tmin=None,

1315 tmax=None,

1316 values_nocheck=None):

1317

1318 if values_nocheck is not None:

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

1320 self.file_size,

1321 self.file_segment, self.file_element,

1322 self.kind_id, codes_safe_str,

1323 self.tmin_seconds, self.tmin_offset,

1324 self.tmax_seconds, self.tmax_offset,

1325 self.deltat) = values_nocheck

1326

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

1328 self.deltat = self.deltat or None

1329 self.raw_content = {}

1330 self.content = None

1331 else:

1332 if tmin is not None:

1333 tmin_seconds, tmin_offset = tsplit(tmin)

1334

1335 if tmax is not None:

1336 tmax_seconds, tmax_offset = tsplit(tmax)

1337

1338 self.kind_id = int(kind_id)

1339 self.codes = codes

1340 self.tmin_seconds = int_or_g_tmin(tmin_seconds)

1341 self.tmin_offset = int(tmin_offset)

1342 self.tmax_seconds = int_or_g_tmax(tmax_seconds)

1343 self.tmax_offset = int(tmax_offset)

1344 self.deltat = float_or_none(deltat)

1345 self.file_path = str_or_none(file_path)

1346 self.file_segment = int_or_none(file_segment)

1347 self.file_element = int_or_none(file_element)

1348 self.file_format = str_or_none(file_format)

1349 self.file_mtime = float_or_none(file_mtime)

1350 self.file_size = int_or_none(file_size)

1351 self.content = content

1352 if raw_content is None:

1353 self.raw_content = {}

1354 else:

1355 self.raw_content = raw_content

1356

1357 Object.__init__(self, init_props=False)

1358

1359 def __eq__(self, other):

1360 return (isinstance(other, Nut) and

1361 self.equality_values == other.equality_values)

1362

1363 def hash(self):

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

1365

1366 def __ne__(self, other):

1367 return not (self == other)

1368

1369 def get_io_backend(self):

1370 from . import io

1371 return io.get_backend(self.file_format)

1372

1373 def file_modified(self):

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

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

1376

1377 @property

1378 def dkey(self):

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

1380

1381 @property

1382 def key(self):

1383 return (

1384 self.file_path,

1385 self.file_segment,

1386 self.file_element,

1387 self.file_mtime)

1388

1389 @property

1390 def equality_values(self):

1391 return (

1392 self.file_segment, self.file_element,

1393 self.kind_id, self.codes,

1394 self.tmin_seconds, self.tmin_offset,

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

1396

1397 def diff(self, other):

1398 names = [

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

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

1401 'deltat']

1402

1403 d = []

1404 for name, a, b in zip(

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

1406

1407 if a != b:

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

1409

1410 return d

1411

1412 @property

1413 def tmin(self):

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

1415

1416 @tmin.setter

1417 def tmin(self, t):

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

1419

1420 @property

1421 def tmax(self):

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

1423

1424 @tmax.setter

1425 def tmax(self, t):

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

1427

1428 @property

1429 def kscale(self):

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

1431 return 0

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

1433

1434 @property

1435 def waveform_kwargs(self):

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

1437

1438 return dict(

1439 network=network,

1440 station=station,

1441 location=location,

1442 channel=channel,

1443 extra=extra,

1444 tmin=self.tmin,

1445 tmax=self.tmax,

1446 deltat=self.deltat)

1447

1448 @property

1449 def waveform_promise_kwargs(self):

1450 return dict(

1451 codes=self.codes,

1452 tmin=self.tmin,

1453 tmax=self.tmax,

1454 deltat=self.deltat)

1455

1456 @property

1457 def station_kwargs(self):

1458 network, station, location = self.codes

1459 return dict(

1460 codes=self.codes,

1461 tmin=tmin_or_none(self.tmin),

1462 tmax=tmax_or_none(self.tmax))

1463

1464 @property

1465 def channel_kwargs(self):

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

1467 return dict(

1468 codes=self.codes,

1469 tmin=tmin_or_none(self.tmin),

1470 tmax=tmax_or_none(self.tmax),

1471 deltat=self.deltat)

1472

1473 @property

1474 def response_kwargs(self):

1475 return dict(

1476 codes=self.codes,

1477 tmin=tmin_or_none(self.tmin),

1478 tmax=tmax_or_none(self.tmax),

1479 deltat=self.deltat)

1480

1481 @property

1482 def event_kwargs(self):

1483 return dict(

1484 name=self.codes,

1485 time=self.tmin,

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

1487

1488 @property

1489 def trace_kwargs(self):

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

1491

1492 return dict(

1493 network=network,

1494 station=station,

1495 location=location,

1496 channel=channel,

1497 extra=extra,

1498 tmin=self.tmin,

1499 tmax=self.tmax-self.deltat,

1500 deltat=self.deltat)

1501

1502 @property

1503 def dummy_trace(self):

1504 return DummyTrace(self)

1505

1506 @property

1507 def summary_entries(self):

1508 if self.tmin == self.tmax:

1509 ts = util.time_to_str(self.tmin)

1510 else:

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

1512 util.time_to_str(self.tmin),

1513 util.time_to_str(self.tmax))

1514

1515 return (

1516 self.__class__.__name__,

1517 to_kind(self.kind_id),

1518 str(self.codes),

1519 ts)

1520

1521 @property

1522 def summary(self):

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

1524

1525

1526def make_waveform_nut(**kwargs):

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

1528

1529

1530def make_waveform_promise_nut(**kwargs):

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

1532

1533

1534def make_station_nut(**kwargs):

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

1536

1537

1538def make_channel_nut(**kwargs):

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

1540

1541

1542def make_response_nut(**kwargs):

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

1544

1545

1546def make_event_nut(**kwargs):

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

1548

1549

1550def group_channels(nuts):

1551 by_ansl = {}

1552 for nut in nuts:

1553 if nut.kind_id != CHANNEL:

1554 continue

1555

1556 ansl = nut.codes[:4]

1557

1558 if ansl not in by_ansl:

1559 by_ansl[ansl] = {}

1560

1561 group = by_ansl[ansl]

1562

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

1564

1565 if k not in group:

1566 group[k] = set()

1567

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

1569

1570 return by_ansl

1571

1572

1573class DummyTrace(object):

1574

1575 def __init__(self, nut):

1576 self.nut = nut

1577 self.codes = nut.codes

1578 self.meta = {}

1579

1580 @property

1581 def tmin(self):

1582 return self.nut.tmin

1583

1584 @property

1585 def tmax(self):

1586 return self.nut.tmax

1587

1588 @property

1589 def deltat(self):

1590 return self.nut.deltat

1591

1592 @property

1593 def nslc_id(self):

1594 return self.codes.nslc

1595

1596 @property

1597 def network(self):

1598 return self.codes.network

1599

1600 @property

1601 def station(self):

1602 return self.codes.station

1603

1604 @property

1605 def location(self):

1606 return self.codes.location

1607

1608 @property

1609 def channel(self):

1610 return self.codes.channel

1611

1612 @property

1613 def extra(self):

1614 return self.codes.extra

1615

1616 def overlaps(self, tmin, tmax):

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

1618

1619

1620def duration_to_str(t):

1621 if t > 24*3600:

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

1623 elif t > 3600:

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

1625 elif t > 60:

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

1627 else:

1628 return '%gs' % t

1629

1630

1631class Coverage(Object):

1632 '''

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

1634 '''

1635 kind_id = Int.T(

1636 help='Content type.')

1637 pattern = Codes.T(

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

1639 'match.')

1640 codes = Codes.T(

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

1642 deltat = Float.T(

1643 help='Sampling interval [s]',

1644 optional=True)

1645 tmin = Timestamp.T(

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

1647 optional=True)

1648 tmax = Timestamp.T(

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

1650 optional=True)

1651 changes = List.T(

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

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

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

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

1656 'value duplicate or redundant data coverage.')

1657

1658 @classmethod

1659 def from_values(cls, args):

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

1661 return cls(

1662 kind_id=kind_id,

1663 pattern=pattern,

1664 codes=codes,

1665 deltat=deltat,

1666 tmin=tmin,

1667 tmax=tmax,

1668 changes=changes)

1669

1670 @property

1671 def summary_entries(self):

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

1673 util.time_to_str(self.tmin),

1674 util.time_to_str(self.tmax))

1675

1676 srate = self.sample_rate

1677

1678 total = self.total

1679

1680 return (

1681 self.__class__.__name__,

1682 to_kind(self.kind_id),

1683 str(self.codes),

1684 ts,

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

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

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

1688

1689 @property

1690 def summary(self):

1691 return util.fmt_summary(

1692 self.summary_entries,

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

1694

1695 @property

1696 def sample_rate(self):

1697 if self.deltat is None:

1698 return None

1699 elif self.deltat == 0.0:

1700 return 0.0

1701 else:

1702 return 1.0 / self.deltat

1703

1704 @property

1705 def labels(self):

1706 srate = self.sample_rate

1707 return (

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

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

1710

1711 @property

1712 def total(self):

1713 total_t = None

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

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

1716

1717 return total_t

1718

1719 def iter_spans(self):

1720 last = None

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

1722 if last is not None:

1723 last_t, last_count = last

1724 if last_count > 0:

1725 yield last_t, t, last_count

1726

1727 last = (t, count)

1728

1729 def iter_uncovered_by(self, other):

1730 a = self

1731 b = other

1732 ia = ib = -1

1733 ca = cb = 0

1734 last = None

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

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

1737 ia += 1

1738 t, ca = a.changes[ia]

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

1740 ib += 1

1741 t, cb = b.changes[ib]

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

1743 ia += 1

1744 t, ca = a.changes[ia]

1745 else:

1746 ib += 1

1747 t, cb = b.changes[ib]

1748

1749 if last is not None:

1750 tl, cal, cbl = last

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

1752 yield tl, t, ia, ib

1753

1754 last = (t, ca, cb)

1755

1756 def iter_uncovered_by_combined(self, other):

1757 ib_last = None

1758 group = None

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

1760 if ib_last is None or ib != ib_last:

1761 if group:

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

1763

1764 group = []

1765

1766 group.append((tmin, tmax))

1767 ib_last = ib

1768

1769 if group:

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

1771

1772

1773__all__ = [

1774 'to_codes',

1775 'to_codes_guess',

1776 'to_codes_simple',

1777 'to_kind',

1778 'to_kinds',

1779 'to_kind_id',

1780 'to_kind_ids',

1781 'CodesError',

1782 'Codes',

1783 'CodesNSLCE',

1784 'CodesNSL',

1785 'CodesX',

1786 'Station',

1787 'Channel',

1788 'Sensor',

1789 'Response',

1790 'Nut',

1791 'Coverage',

1792 'WaveformPromise',

1793]