1# http://pyrocko.org - GPLv3
2#
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'''
19from __future__ import absolute_import, print_function
21import re
22import fnmatch
23import math
24import hashlib
25from os import urandom
26from base64 import urlsafe_b64encode
27from collections import defaultdict, namedtuple
29import numpy as num
31from pyrocko import util
32from pyrocko.guts import Object, SObject, String, Timestamp, Float, Int, \
33 Unicode, Tuple, List, StringChoice, Any, Dict, clone
34from pyrocko.model import squirrel_content, Location
35from pyrocko.response import FrequencyResponse, MultiplyResponse, \
36 IntegrationResponse, DifferentiationResponse, simplify_responses, \
37 FrequencyResponseCheckpoint
39from .error import ConversionError, SquirrelError
41d2r = num.pi / 180.
42r2d = 1.0 / d2r
45guts_prefix = 'squirrel'
48def mkvec(x, y, z):
49 return num.array([x, y, z], dtype=float)
52def are_orthogonal(vecs, eps=0.01):
53 return all(abs(
54 num.dot(vecs[i], vecs[j]) < eps
55 for (i, j) in [(0, 1), (1, 2), (2, 0)]))
58g_codes_pool = {}
61class CodesError(SquirrelError):
62 pass
65class Codes(SObject):
66 pass
69def normalize_nslce(*args, **kwargs):
70 if args and kwargs:
71 raise ValueError('Either *args or **kwargs accepted, not both.')
73 if len(args) == 1:
74 if isinstance(args[0], str):
75 args = tuple(args[0].split('.'))
76 elif isinstance(args[0], tuple):
77 args = args[0]
78 else:
79 raise ValueError('Invalid argument type: %s' % type(args[0]))
81 nargs = len(args)
82 if nargs == 5:
83 t = args
85 elif nargs == 4:
86 t = args + ('',)
88 elif nargs == 0:
89 d = dict(
90 network='',
91 station='',
92 location='',
93 channel='',
94 extra='')
96 d.update(kwargs)
97 t = tuple(kwargs.get(k, '') for k in (
98 'network', 'station', 'location', 'channel', 'extra'))
100 else:
101 raise CodesError(
102 'Does not match NSLC or NSLCE codes pattern: %s' % '.'.join(args))
104 if '.'.join(t).count('.') != 4:
105 raise CodesError(
106 'Codes may not contain a ".": "%s", "%s", "%s", "%s", "%s"' % t)
108 return t
111CodesNSLCEBase = namedtuple(
112 'CodesNSLCEBase', [
113 'network', 'station', 'location', 'channel', 'extra'])
116class CodesNSLCE(CodesNSLCEBase, Codes):
117 '''
118 Codes denominating a seismic channel (NSLC or NSLCE).
120 FDSN/SEED style NET.STA.LOC.CHA is accepted or NET.STA.LOC.CHA.EXTRA, where
121 the EXTRA part in the latter form can be used to identify a custom
122 processing applied to a channel.
123 '''
125 __slots__ = ()
126 __hash__ = CodesNSLCEBase.__hash__
128 as_dict = CodesNSLCEBase._asdict
130 def __new__(cls, *args, safe_str=None, **kwargs):
131 nargs = len(args)
132 if nargs == 1 and isinstance(args[0], CodesNSLCE):
133 return args[0]
134 elif nargs == 1 and isinstance(args[0], CodesNSL):
135 t = (args[0].nsl) + ('*', '*')
136 elif nargs == 1 and isinstance(args[0], CodesX):
137 t = ('*', '*', '*', '*', '*')
138 elif safe_str is not None:
139 t = safe_str.split('.')
140 else:
141 t = normalize_nslce(*args, **kwargs)
143 x = CodesNSLCEBase.__new__(cls, *t)
144 return g_codes_pool.setdefault(x, x)
146 def __init__(self, *args, **kwargs):
147 Codes.__init__(self)
149 def __str__(self):
150 if self.extra == '':
151 return '.'.join(self[:-1])
152 else:
153 return '.'.join(self)
155 def __eq__(self, other):
156 if not isinstance(other, CodesNSLCE):
157 other = CodesNSLCE(other)
159 return CodesNSLCEBase.__eq__(self, other)
161 def matches(self, pattern):
162 if not isinstance(pattern, CodesNSLCE):
163 pattern = CodesNSLCE(pattern)
165 return match_codes(pattern, self)
167 @property
168 def safe_str(self):
169 return '.'.join(self)
171 @property
172 def nslce(self):
173 return self[:4]
175 @property
176 def nslc(self):
177 return self[:4]
179 @property
180 def nsl(self):
181 return self[:3]
183 @property
184 def ns(self):
185 return self[:2]
187 def as_tuple(self):
188 return tuple(self)
190 def replace(self, **kwargs):
191 x = CodesNSLCEBase._replace(self, **kwargs)
192 return g_codes_pool.setdefault(x, x)
195def normalize_nsl(*args, **kwargs):
196 if args and kwargs:
197 raise ValueError('Either *args or **kwargs accepted, not both.')
199 if len(args) == 1:
200 if isinstance(args[0], str):
201 args = tuple(args[0].split('.'))
202 elif isinstance(args[0], tuple):
203 args = args[0]
204 else:
205 raise ValueError('Invalid argument type: %s' % type(args[0]))
207 nargs = len(args)
208 if nargs == 3:
209 t = args
211 elif nargs == 0:
212 d = dict(
213 network='',
214 station='',
215 location='')
217 d.update(kwargs)
218 t = tuple(kwargs.get(k, '') for k in (
219 'network', 'station', 'location'))
221 else:
222 raise CodesError(
223 'Does not match NSL codes pattern: %s' % '.'.join(args))
225 if '.'.join(t).count('.') != 2:
226 raise CodesError(
227 'Codes may not contain a ".": "%s", "%s", "%s"' % t)
229 return t
232CodesNSLBase = namedtuple(
233 'CodesNSLBase', [
234 'network', 'station', 'location'])
237class CodesNSL(CodesNSLBase, Codes):
238 '''
239 Codes denominating a seismic station (NSL).
241 NET.STA.LOC is accepted, slightly different from SEED/StationXML, where
242 LOC is part of the channel. By setting location='*' is possible to get
243 compatible behaviour in most cases.
244 '''
246 __slots__ = ()
247 __hash__ = CodesNSLBase.__hash__
249 as_dict = CodesNSLBase._asdict
251 def __new__(cls, *args, safe_str=None, **kwargs):
252 nargs = len(args)
253 if nargs == 1 and isinstance(args[0], CodesNSL):
254 return args[0]
255 elif nargs == 1 and isinstance(args[0], CodesNSLCE):
256 t = args[0].nsl
257 elif nargs == 1 and isinstance(args[0], CodesX):
258 t = ('*', '*', '*')
259 elif safe_str is not None:
260 t = safe_str.split('.')
261 else:
262 t = normalize_nsl(*args, **kwargs)
264 x = CodesNSLBase.__new__(cls, *t)
265 return g_codes_pool.setdefault(x, x)
267 def __init__(self, *args, **kwargs):
268 Codes.__init__(self)
270 def __str__(self):
271 return '.'.join(self)
273 def __eq__(self, other):
274 if not isinstance(other, CodesNSL):
275 other = CodesNSL(other)
277 return CodesNSLBase.__eq__(self, other)
279 def matches(self, pattern):
280 if not isinstance(pattern, CodesNSL):
281 pattern = CodesNSL(pattern)
283 return match_codes(pattern, self)
285 @property
286 def safe_str(self):
287 return '.'.join(self)
289 @property
290 def ns(self):
291 return self[:2]
293 @property
294 def nsl(self):
295 return self[:3]
297 def as_tuple(self):
298 return tuple(self)
300 def replace(self, **kwargs):
301 x = CodesNSLBase._replace(self, **kwargs)
302 return g_codes_pool.setdefault(x, x)
305CodesXBase = namedtuple(
306 'CodesXBase', [
307 'name'])
310class CodesX(CodesXBase, Codes):
311 '''
312 General purpose codes for anything other than channels or stations.
313 '''
315 __slots__ = ()
316 __hash__ = CodesXBase.__hash__
317 __eq__ = CodesXBase.__eq__
319 as_dict = CodesXBase._asdict
321 def __new__(cls, name='', safe_str=None):
322 if isinstance(name, CodesX):
323 return name
324 elif isinstance(name, (CodesNSLCE, CodesNSL)):
325 name = '*'
326 elif safe_str is not None:
327 name = safe_str
328 else:
329 if '.' in name:
330 raise CodesError('Code may not contain a ".": %s' % name)
332 x = CodesXBase.__new__(cls, name)
333 return g_codes_pool.setdefault(x, x)
335 def __init__(self, *args, **kwargs):
336 Codes.__init__(self)
338 def __str__(self):
339 return '.'.join(self)
341 @property
342 def safe_str(self):
343 return '.'.join(self)
345 @property
346 def ns(self):
347 return self[:2]
349 def as_tuple(self):
350 return tuple(self)
352 def replace(self, **kwargs):
353 x = CodesXBase._replace(self, **kwargs)
354 return g_codes_pool.setdefault(x, x)
357g_codes_patterns = {}
360def match_codes(pattern, codes):
361 spattern = pattern.safe_str
362 scodes = codes.safe_str
363 if spattern not in g_codes_patterns:
364 rpattern = re.compile(fnmatch.translate(spattern), re.I)
365 g_codes_patterns[spattern] = rpattern
367 rpattern = g_codes_patterns[spattern]
368 return bool(rpattern.match(scodes))
371g_content_kinds = [
372 'undefined',
373 'waveform',
374 'station',
375 'channel',
376 'response',
377 'event',
378 'waveform_promise']
381g_codes_classes = [
382 CodesX,
383 CodesNSLCE,
384 CodesNSL,
385 CodesNSLCE,
386 CodesNSLCE,
387 CodesX,
388 CodesNSLCE]
390g_codes_classes_ndot = {
391 0: CodesX,
392 2: CodesNSL,
393 3: CodesNSLCE,
394 4: CodesNSLCE}
397def to_codes_simple(kind_id, codes_safe_str):
398 return g_codes_classes[kind_id](safe_str=codes_safe_str)
401def to_codes(kind_id, obj):
402 return g_codes_classes[kind_id](obj)
405def to_codes_guess(s):
406 try:
407 return g_codes_classes_ndot[s.count('.')](s)
408 except KeyError:
409 raise CodesError('Cannot guess codes type: %s' % s)
412# derived list class to enable detection of already preprocessed codes patterns
413class codes_patterns_list(list):
414 pass
417def codes_patterns_for_kind(kind_id, codes):
418 if isinstance(codes, codes_patterns_list):
419 return codes
421 if isinstance(codes, list):
422 lcodes = codes_patterns_list()
423 for sc in codes:
424 lcodes.extend(codes_patterns_for_kind(kind_id, sc))
426 return lcodes
428 codes = to_codes(kind_id, codes)
430 lcodes = codes_patterns_list()
431 lcodes.append(codes)
432 if kind_id == STATION:
433 lcodes.append(codes.replace(location='[*]'))
435 return lcodes
438g_content_kind_ids = (
439 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT,
440 WAVEFORM_PROMISE) = range(len(g_content_kinds))
443g_tmin, g_tmax = util.get_working_system_time_range()[:2]
446try:
447 g_tmin_queries = max(g_tmin, util.str_to_time_fillup('1900-01-01'))
448except Exception:
449 g_tmin_queries = g_tmin
452def to_kind(kind_id):
453 return g_content_kinds[kind_id]
456def to_kinds(kind_ids):
457 return [g_content_kinds[kind_id] for kind_id in kind_ids]
460def to_kind_id(kind):
461 return g_content_kinds.index(kind)
464def to_kind_ids(kinds):
465 return [g_content_kinds.index(kind) for kind in kinds]
468g_kind_mask_all = 0xff
471def to_kind_mask(kinds):
472 if kinds:
473 return sum(1 << kind_id for kind_id in to_kind_ids(kinds))
474 else:
475 return g_kind_mask_all
478def str_or_none(x):
479 if x is None:
480 return None
481 else:
482 return str(x)
485def float_or_none(x):
486 if x is None:
487 return None
488 else:
489 return float(x)
492def int_or_none(x):
493 if x is None:
494 return None
495 else:
496 return int(x)
499def int_or_g_tmin(x):
500 if x is None:
501 return g_tmin
502 else:
503 return int(x)
506def int_or_g_tmax(x):
507 if x is None:
508 return g_tmax
509 else:
510 return int(x)
513def tmin_or_none(x):
514 if x == g_tmin:
515 return None
516 else:
517 return x
520def tmax_or_none(x):
521 if x == g_tmax:
522 return None
523 else:
524 return x
527def time_or_none_to_str(x):
528 if x is None:
529 return '...'.ljust(17)
530 else:
531 return util.time_to_str(x)
534def codes_to_str_abbreviated(codes, indent=' '):
535 codes = [str(x) for x in codes]
537 if len(codes) > 20:
538 scodes = '\n' + util.ewrap(codes[:10], indent=indent) \
539 + '\n%s[%i more]\n' % (indent, len(codes) - 20) \
540 + util.ewrap(codes[-10:], indent=' ')
541 else:
542 scodes = '\n' + util.ewrap(codes, indent=indent) \
543 if codes else '<none>'
545 return scodes
548g_offset_time_unit_inv = 1000000000
549g_offset_time_unit = 1.0 / g_offset_time_unit_inv
552def tsplit(t):
553 if t is None:
554 return None, 0
556 t = util.to_time_float(t)
557 if type(t) is float:
558 t = round(t, 5)
559 else:
560 t = round(t, 9)
562 seconds = num.floor(t)
563 offset = t - seconds
564 return int(seconds), int(round(offset * g_offset_time_unit_inv))
567def tjoin(seconds, offset):
568 if seconds is None:
569 return None
571 return util.to_time_float(seconds) \
572 + util.to_time_float(offset*g_offset_time_unit)
575tscale_min = 1
576tscale_max = 365 * 24 * 3600 # last edge is one above
577tscale_logbase = 20
579tscale_edges = [tscale_min]
580while True:
581 tscale_edges.append(tscale_edges[-1]*tscale_logbase)
582 if tscale_edges[-1] >= tscale_max:
583 break
586tscale_edges = num.array(tscale_edges)
589def tscale_to_kscale(tscale):
591 # 0 <= x < tscale_edges[1]: 0
592 # tscale_edges[1] <= x < tscale_edges[2]: 1
593 # ...
594 # tscale_edges[len(tscale_edges)-1] <= x: len(tscale_edges)
596 return int(num.searchsorted(tscale_edges, tscale))
599@squirrel_content
600class Station(Location):
601 '''
602 A seismic station.
603 '''
605 codes = CodesNSL.T()
607 tmin = Timestamp.T(optional=True)
608 tmax = Timestamp.T(optional=True)
610 description = Unicode.T(optional=True)
612 def __init__(self, **kwargs):
613 kwargs['codes'] = CodesNSL(kwargs['codes'])
614 Location.__init__(self, **kwargs)
616 @property
617 def time_span(self):
618 return (self.tmin, self.tmax)
620 def get_pyrocko_station(self):
621 from pyrocko import model
622 return model.Station(*self._get_pyrocko_station_args())
624 def _get_pyrocko_station_args(self):
625 return (
626 self.codes.network,
627 self.codes.station,
628 self.codes.location,
629 self.lat,
630 self.lon,
631 self.elevation,
632 self.depth,
633 self.north_shift,
634 self.east_shift)
637@squirrel_content
638class ChannelBase(Location):
639 codes = CodesNSLCE.T()
641 tmin = Timestamp.T(optional=True)
642 tmax = Timestamp.T(optional=True)
644 deltat = Float.T(optional=True)
646 def __init__(self, **kwargs):
647 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
648 Location.__init__(self, **kwargs)
650 @property
651 def time_span(self):
652 return (self.tmin, self.tmax)
654 def _get_sensor_codes(self):
655 return self.codes.replace(
656 channel=self.codes.channel[:-1] + '?')
658 def _get_sensor_args(self):
659 def getattr_rep(k):
660 if k == 'codes':
661 return self._get_sensor_codes()
662 else:
663 return getattr(self, k)
665 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
667 def _get_channel_args(self, component):
668 def getattr_rep(k):
669 if k == 'codes':
670 return self.codes.replace(
671 channel=self.codes.channel[:-1] + component)
672 else:
673 return getattr(self, k)
675 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
677 def _get_pyrocko_station_args(self):
678 return (
679 self.codes.network,
680 self.codes.station,
681 self.codes.location,
682 self.lat,
683 self.lon,
684 self.elevation,
685 self.depth,
686 self.north_shift,
687 self.east_shift)
690class Channel(ChannelBase):
691 '''
692 A channel of a seismic station.
693 '''
695 dip = Float.T(optional=True)
696 azimuth = Float.T(optional=True)
698 def get_pyrocko_channel(self):
699 from pyrocko import model
700 return model.Channel(*self._get_pyrocko_channel_args())
702 def _get_pyrocko_channel_args(self):
703 return (
704 self.codes.channel,
705 self.azimuth,
706 self.dip)
708 @property
709 def orientation_enz(self):
710 if None in (self.azimuth, self.dip):
711 return None
713 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
714 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
715 d = math.sin(self.dip*d2r)
716 return mkvec(e, n, -d)
719def cut_intervals(channels):
720 channels = list(channels)
721 times = set()
722 for channel in channels:
723 if channel.tmin is not None:
724 times.add(channel.tmin)
725 if channel.tmax is not None:
726 times.add(channel.tmax)
728 times = sorted(times)
730 if any(channel.tmin is None for channel in channels):
731 times[0:0] = [None]
733 if any(channel.tmax is None for channel in channels):
734 times.append(None)
736 if len(times) <= 2:
737 return channels
739 channels_out = []
740 for channel in channels:
741 for i in range(len(times)-1):
742 tmin = times[i]
743 tmax = times[i+1]
744 if ((channel.tmin is None or (
745 tmin is not None and channel.tmin <= tmin))
746 and (channel.tmax is None or (
747 tmax is not None and tmax <= channel.tmax))):
749 channel_out = clone(channel)
750 channel_out.tmin = tmin
751 channel_out.tmax = tmax
752 channels_out.append(channel_out)
754 return channels_out
757class Sensor(ChannelBase):
758 '''
759 Representation of a channel group.
760 '''
762 channels = List.T(Channel.T())
764 @classmethod
765 def from_channels(cls, channels):
766 groups = defaultdict(list)
767 for channel in channels:
768 groups[channel._get_sensor_codes()].append(channel)
770 channels_cut = []
771 for group in groups.values():
772 channels_cut.extend(cut_intervals(group))
774 groups = defaultdict(list)
775 for channel in channels_cut:
776 groups[channel._get_sensor_args()].append(channel)
778 return [
779 cls(channels=channels,
780 **dict(zip(ChannelBase.T.propnames, args)))
781 for args, channels in groups.items()]
783 def channel_vectors(self):
784 return num.vstack(
785 [channel.orientation_enz for channel in self.channels])
787 def projected_channels(self, system, component_names):
788 return [
789 Channel(
790 azimuth=math.atan2(e, n) * r2d,
791 dip=-math.asin(z) * r2d,
792 **dict(zip(
793 ChannelBase.T.propnames,
794 self._get_channel_args(comp))))
795 for comp, (e, n, z) in zip(component_names, system)]
797 def matrix_to(self, system, epsilon=1e-7):
798 m = num.dot(system, self.channel_vectors().T)
799 m[num.abs(m) < epsilon] = 0.0
800 return m
802 def projection_to(self, system, component_names):
803 return (
804 self.matrix_to(system),
805 self.channels,
806 self.projected_channels(system, component_names))
808 def projection_to_enz(self):
809 return self.projection_to(num.identity(3), 'ENZ')
811 def projection_to_trz(self, source, azimuth=None):
812 if azimuth is not None:
813 assert source is None
814 else:
815 azimuth = source.azibazi_to(self)[1] + 180.
817 return self.projection_to(num.array([
818 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
819 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
820 [0., 0., 1.]], dtype=float), 'TRZ')
822 def project_to_enz(self, traces):
823 from pyrocko import trace
825 matrix, in_channels, out_channels = self.projection_to_enz()
827 return trace.project(traces, matrix, in_channels, out_channels)
829 def project_to_trz(self, source, traces, azimuth=None):
830 from pyrocko import trace
832 matrix, in_channels, out_channels = self.projection_to_trz(
833 source, azimuth=azimuth)
835 return trace.project(traces, matrix, in_channels, out_channels)
838observational_quantities = [
839 'acceleration', 'velocity', 'displacement', 'pressure',
840 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
841 'temperature']
844technical_quantities = [
845 'voltage', 'counts']
848class QuantityType(StringChoice):
849 '''
850 Choice of observational or technical quantity.
852 SI units are used for all quantities, where applicable.
853 '''
854 choices = observational_quantities + technical_quantities
857class ResponseStage(Object):
858 '''
859 Representation of a response stage.
861 Components of a seismic recording system are represented as a sequence of
862 response stages, e.g. sensor, pre-amplifier, digitizer, digital
863 downsampling.
864 '''
865 input_quantity = QuantityType.T(optional=True)
866 input_sample_rate = Float.T(optional=True)
867 output_quantity = QuantityType.T(optional=True)
868 output_sample_rate = Float.T(optional=True)
869 elements = List.T(FrequencyResponse.T())
870 log = List.T(Tuple.T(3, String.T()))
872 @property
873 def stage_type(self):
874 if self.input_quantity in observational_quantities \
875 and self.output_quantity in observational_quantities:
876 return 'conversion'
878 if self.input_quantity in observational_quantities \
879 and self.output_quantity == 'voltage':
880 return 'sensor'
882 elif self.input_quantity == 'voltage' \
883 and self.output_quantity == 'voltage':
884 return 'electronics'
886 elif self.input_quantity == 'voltage' \
887 and self.output_quantity == 'counts':
888 return 'digitizer'
890 elif self.input_quantity == 'counts' \
891 and self.output_quantity == 'counts' \
892 and self.input_sample_rate != self.output_sample_rate:
893 return 'decimation'
895 elif self.input_quantity in observational_quantities \
896 and self.output_quantity == 'counts':
897 return 'combination'
899 else:
900 return 'unknown'
902 @property
903 def summary(self):
904 irate = self.input_sample_rate
905 orate = self.output_sample_rate
906 factor = None
907 if irate and orate:
908 factor = irate / orate
909 return 'ResponseStage, ' + (
910 '%s%s => %s%s%s' % (
911 self.input_quantity or '?',
912 ' @ %g Hz' % irate if irate else '',
913 self.output_quantity or '?',
914 ' @ %g Hz' % orate if orate else '',
915 ' :%g' % factor if factor else '')
916 )
918 def get_effective(self):
919 return MultiplyResponse(responses=list(self.elements))
922D = 'displacement'
923V = 'velocity'
924A = 'acceleration'
926g_converters = {
927 (V, D): IntegrationResponse(1),
928 (A, D): IntegrationResponse(2),
929 (D, V): DifferentiationResponse(1),
930 (A, V): IntegrationResponse(1),
931 (D, A): DifferentiationResponse(2),
932 (V, A): DifferentiationResponse(1)}
935def response_converters(input_quantity, output_quantity):
936 if input_quantity is None or input_quantity == output_quantity:
937 return []
939 if output_quantity is None:
940 raise ConversionError('Unspecified target quantity.')
942 try:
943 return [g_converters[input_quantity, output_quantity]]
945 except KeyError:
946 raise ConversionError('No rule to convert from "%s" to "%s".' % (
947 input_quantity, output_quantity))
950@squirrel_content
951class Response(Object):
952 '''
953 The instrument response of a seismic station channel.
954 '''
956 codes = CodesNSLCE.T()
957 tmin = Timestamp.T(optional=True)
958 tmax = Timestamp.T(optional=True)
960 stages = List.T(ResponseStage.T())
961 checkpoints = List.T(FrequencyResponseCheckpoint.T())
963 deltat = Float.T(optional=True)
964 log = List.T(Tuple.T(3, String.T()))
966 def __init__(self, **kwargs):
967 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
968 Object.__init__(self, **kwargs)
970 @property
971 def time_span(self):
972 return (self.tmin, self.tmax)
974 @property
975 def nstages(self):
976 return len(self.stages)
978 @property
979 def input_quantity(self):
980 return self.stages[0].input_quantity if self.stages else None
982 @property
983 def output_quantity(self):
984 return self.stages[-1].input_quantity if self.stages else None
986 @property
987 def output_sample_rate(self):
988 return self.stages[-1].output_sample_rate if self.stages else None
990 @property
991 def stages_summary(self):
992 def grouped(xs):
993 xs = list(xs)
994 g = []
995 for i in range(len(xs)):
996 g.append(xs[i])
997 if i+1 < len(xs) and xs[i+1] != xs[i]:
998 yield g
999 g = []
1001 if g:
1002 yield g
1004 return '+'.join(
1005 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1006 for g in grouped(stage.stage_type for stage in self.stages))
1008 @property
1009 def summary(self):
1010 orate = self.output_sample_rate
1011 return '%s %-16s %s' % (
1012 self.__class__.__name__, self.str_codes, self.str_time_span) \
1013 + ', ' + ', '.join((
1014 '%s => %s' % (
1015 self.input_quantity or '?', self.output_quantity or '?')
1016 + (' @ %g Hz' % orate if orate else ''),
1017 self.stages_summary,
1018 ))
1020 def get_effective(self, input_quantity=None):
1021 try:
1022 elements = response_converters(input_quantity, self.input_quantity)
1023 except ConversionError as e:
1024 raise ConversionError(str(e) + ' (%s)' % self.summary)
1026 elements.extend(
1027 stage.get_effective() for stage in self.stages)
1029 return MultiplyResponse(responses=simplify_responses(elements))
1032@squirrel_content
1033class Event(Object):
1034 '''
1035 A seismic event.
1036 '''
1038 name = String.T(optional=True)
1039 time = Timestamp.T()
1040 duration = Float.T(optional=True)
1042 lat = Float.T()
1043 lon = Float.T()
1044 depth = Float.T(optional=True)
1046 magnitude = Float.T(optional=True)
1048 def get_pyrocko_event(self):
1049 from pyrocko import model
1050 model.Event(
1051 name=self.name,
1052 time=self.time,
1053 lat=self.lat,
1054 lon=self.lon,
1055 depth=self.depth,
1056 magnitude=self.magnitude,
1057 duration=self.duration)
1059 @property
1060 def time_span(self):
1061 return (self.time, self.time)
1064def ehash(s):
1065 return hashlib.sha1(s.encode('utf8')).hexdigest()
1068def random_name(n=8):
1069 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1072@squirrel_content
1073class WaveformPromise(Object):
1074 '''
1075 Information about a waveform potentially downloadable from a remote site.
1077 In the Squirrel framework, waveform promises are used to permit download of
1078 selected waveforms from a remote site. They are typically generated by
1079 calls to
1080 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1081 Waveform promises are inserted and indexed in the database similar to
1082 normal waveforms. When processing a waveform query, e.g. from
1083 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1084 waveform is available for the queried time span, a matching promise can be
1085 resolved, i.e. an attempt is made to download the waveform from the remote
1086 site. The promise is removed after the download attempt (except when a
1087 network error occurs). This prevents Squirrel from making unnecessary
1088 queries for waveforms missing at the remote site.
1089 '''
1091 codes = CodesNSLCE.T()
1092 tmin = Timestamp.T()
1093 tmax = Timestamp.T()
1095 deltat = Float.T(optional=True)
1097 source_hash = String.T()
1099 def __init__(self, **kwargs):
1100 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1101 Object.__init__(self, **kwargs)
1103 @property
1104 def time_span(self):
1105 return (self.tmin, self.tmax)
1108class InvalidWaveform(Exception):
1109 pass
1112class WaveformOrder(Object):
1113 '''
1114 Waveform request information.
1115 '''
1117 source_id = String.T()
1118 codes = CodesNSLCE.T()
1119 deltat = Float.T()
1120 tmin = Timestamp.T()
1121 tmax = Timestamp.T()
1122 gaps = List.T(Tuple.T(2, Timestamp.T()))
1124 @property
1125 def client(self):
1126 return self.source_id.split(':')[1]
1128 def describe(self, site='?'):
1129 return '%s:%s %s [%s - %s]' % (
1130 self.client, site, str(self.codes),
1131 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1133 def validate(self, tr):
1134 if tr.ydata.size == 0:
1135 raise InvalidWaveform(
1136 'waveform with zero data samples')
1138 if tr.deltat != self.deltat:
1139 raise InvalidWaveform(
1140 'incorrect sampling interval - waveform: %g s, '
1141 'meta-data: %g s' % (
1142 tr.deltat, self.deltat))
1144 if not num.all(num.isfinite(tr.ydata)):
1145 raise InvalidWaveform('waveform has NaN values')
1148def order_summary(orders):
1149 codes_list = sorted(set(order.codes for order in orders))
1150 if len(codes_list) > 3:
1151 return '%i order%s: %s - %s' % (
1152 len(orders),
1153 util.plural_s(orders),
1154 str(codes_list[0]),
1155 str(codes_list[1]))
1157 else:
1158 return '%i order%s: %s' % (
1159 len(orders),
1160 util.plural_s(orders),
1161 ', '.join(str(codes) for codes in codes_list))
1164class Nut(Object):
1165 '''
1166 Index entry referencing an elementary piece of content.
1168 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1169 meta-information about individual pieces of waveforms, stations, channels,
1170 etc. together with the information where it was found or generated.
1171 '''
1173 file_path = String.T(optional=True)
1174 file_format = String.T(optional=True)
1175 file_mtime = Timestamp.T(optional=True)
1176 file_size = Int.T(optional=True)
1178 file_segment = Int.T(optional=True)
1179 file_element = Int.T(optional=True)
1181 kind_id = Int.T()
1182 codes = Codes.T()
1184 tmin_seconds = Int.T(default=0)
1185 tmin_offset = Int.T(default=0, optional=True)
1186 tmax_seconds = Int.T(default=0)
1187 tmax_offset = Int.T(default=0, optional=True)
1189 deltat = Float.T(default=0.0)
1191 content = Any.T(optional=True)
1192 raw_content = Dict.T(String.T(), Any.T())
1194 content_in_db = False
1196 def __init__(
1197 self,
1198 file_path=None,
1199 file_format=None,
1200 file_mtime=None,
1201 file_size=None,
1202 file_segment=None,
1203 file_element=None,
1204 kind_id=0,
1205 codes=CodesX(''),
1206 tmin_seconds=None,
1207 tmin_offset=0,
1208 tmax_seconds=None,
1209 tmax_offset=0,
1210 deltat=None,
1211 content=None,
1212 raw_content=None,
1213 tmin=None,
1214 tmax=None,
1215 values_nocheck=None):
1217 if values_nocheck is not None:
1218 (self.file_path, self.file_format, self.file_mtime,
1219 self.file_size,
1220 self.file_segment, self.file_element,
1221 self.kind_id, codes_safe_str,
1222 self.tmin_seconds, self.tmin_offset,
1223 self.tmax_seconds, self.tmax_offset,
1224 self.deltat) = values_nocheck
1226 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1227 self.deltat = self.deltat or None
1228 self.raw_content = {}
1229 self.content = None
1230 else:
1231 if tmin is not None:
1232 tmin_seconds, tmin_offset = tsplit(tmin)
1234 if tmax is not None:
1235 tmax_seconds, tmax_offset = tsplit(tmax)
1237 self.kind_id = int(kind_id)
1238 self.codes = codes
1239 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1240 self.tmin_offset = int(tmin_offset)
1241 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1242 self.tmax_offset = int(tmax_offset)
1243 self.deltat = float_or_none(deltat)
1244 self.file_path = str_or_none(file_path)
1245 self.file_segment = int_or_none(file_segment)
1246 self.file_element = int_or_none(file_element)
1247 self.file_format = str_or_none(file_format)
1248 self.file_mtime = float_or_none(file_mtime)
1249 self.file_size = int_or_none(file_size)
1250 self.content = content
1251 if raw_content is None:
1252 self.raw_content = {}
1253 else:
1254 self.raw_content = raw_content
1256 Object.__init__(self, init_props=False)
1258 def __eq__(self, other):
1259 return (isinstance(other, Nut) and
1260 self.equality_values == other.equality_values)
1262 def hash(self):
1263 return ehash(','.join(str(x) for x in self.key))
1265 def __ne__(self, other):
1266 return not (self == other)
1268 def get_io_backend(self):
1269 from . import io
1270 return io.get_backend(self.file_format)
1272 def file_modified(self):
1273 return self.get_io_backend().get_stats(self.file_path) \
1274 != (self.file_mtime, self.file_size)
1276 @property
1277 def dkey(self):
1278 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1280 @property
1281 def key(self):
1282 return (
1283 self.file_path,
1284 self.file_segment,
1285 self.file_element,
1286 self.file_mtime)
1288 @property
1289 def equality_values(self):
1290 return (
1291 self.file_segment, self.file_element,
1292 self.kind_id, self.codes,
1293 self.tmin_seconds, self.tmin_offset,
1294 self.tmax_seconds, self.tmax_offset, self.deltat)
1296 def diff(self, other):
1297 names = [
1298 'file_segment', 'file_element', 'kind_id', 'codes',
1299 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1300 'deltat']
1302 d = []
1303 for name, a, b in zip(
1304 names, self.equality_values, other.equality_values):
1306 if a != b:
1307 d.append((name, a, b))
1309 return d
1311 @property
1312 def tmin(self):
1313 return tjoin(self.tmin_seconds, self.tmin_offset)
1315 @tmin.setter
1316 def tmin(self, t):
1317 self.tmin_seconds, self.tmin_offset = tsplit(t)
1319 @property
1320 def tmax(self):
1321 return tjoin(self.tmax_seconds, self.tmax_offset)
1323 @tmax.setter
1324 def tmax(self, t):
1325 self.tmax_seconds, self.tmax_offset = tsplit(t)
1327 @property
1328 def kscale(self):
1329 if self.tmin_seconds is None or self.tmax_seconds is None:
1330 return 0
1331 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1333 @property
1334 def waveform_kwargs(self):
1335 network, station, location, channel, extra = self.codes
1337 return dict(
1338 network=network,
1339 station=station,
1340 location=location,
1341 channel=channel,
1342 extra=extra,
1343 tmin=self.tmin,
1344 tmax=self.tmax,
1345 deltat=self.deltat)
1347 @property
1348 def waveform_promise_kwargs(self):
1349 return dict(
1350 codes=self.codes,
1351 tmin=self.tmin,
1352 tmax=self.tmax,
1353 deltat=self.deltat)
1355 @property
1356 def station_kwargs(self):
1357 network, station, location = self.codes
1358 return dict(
1359 codes=self.codes,
1360 tmin=tmin_or_none(self.tmin),
1361 tmax=tmax_or_none(self.tmax))
1363 @property
1364 def channel_kwargs(self):
1365 network, station, location, channel, extra = self.codes
1366 return dict(
1367 codes=self.codes,
1368 tmin=tmin_or_none(self.tmin),
1369 tmax=tmax_or_none(self.tmax),
1370 deltat=self.deltat)
1372 @property
1373 def response_kwargs(self):
1374 return dict(
1375 codes=self.codes,
1376 tmin=tmin_or_none(self.tmin),
1377 tmax=tmax_or_none(self.tmax),
1378 deltat=self.deltat)
1380 @property
1381 def event_kwargs(self):
1382 return dict(
1383 name=self.codes,
1384 time=self.tmin,
1385 duration=(self.tmax - self.tmin) or None)
1387 @property
1388 def trace_kwargs(self):
1389 network, station, location, channel, extra = self.codes
1391 return dict(
1392 network=network,
1393 station=station,
1394 location=location,
1395 channel=channel,
1396 extra=extra,
1397 tmin=self.tmin,
1398 tmax=self.tmax-self.deltat,
1399 deltat=self.deltat)
1401 @property
1402 def dummy_trace(self):
1403 return DummyTrace(self)
1405 @property
1406 def summary(self):
1407 if self.tmin == self.tmax:
1408 ts = util.time_to_str(self.tmin)
1409 else:
1410 ts = '%s - %s' % (
1411 util.time_to_str(self.tmin),
1412 util.time_to_str(self.tmax))
1414 return ' '.join((
1415 ('%s,' % to_kind(self.kind_id)).ljust(9),
1416 ('%s,' % str(self.codes)).ljust(18),
1417 ts))
1420def make_waveform_nut(**kwargs):
1421 return Nut(kind_id=WAVEFORM, **kwargs)
1424def make_waveform_promise_nut(**kwargs):
1425 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1428def make_station_nut(**kwargs):
1429 return Nut(kind_id=STATION, **kwargs)
1432def make_channel_nut(**kwargs):
1433 return Nut(kind_id=CHANNEL, **kwargs)
1436def make_response_nut(**kwargs):
1437 return Nut(kind_id=RESPONSE, **kwargs)
1440def make_event_nut(**kwargs):
1441 return Nut(kind_id=EVENT, **kwargs)
1444def group_channels(nuts):
1445 by_ansl = {}
1446 for nut in nuts:
1447 if nut.kind_id != CHANNEL:
1448 continue
1450 ansl = nut.codes[:4]
1452 if ansl not in by_ansl:
1453 by_ansl[ansl] = {}
1455 group = by_ansl[ansl]
1457 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1459 if k not in group:
1460 group[k] = set()
1462 group.add(nut.codes[4])
1464 return by_ansl
1467class DummyTrace(object):
1469 def __init__(self, nut):
1470 self.nut = nut
1471 self.codes = nut.codes
1472 self.meta = {}
1474 @property
1475 def tmin(self):
1476 return self.nut.tmin
1478 @property
1479 def tmax(self):
1480 return self.nut.tmax
1482 @property
1483 def deltat(self):
1484 return self.nut.deltat
1486 @property
1487 def nslc_id(self):
1488 return self.codes.nslc
1490 @property
1491 def network(self):
1492 return self.codes.network
1494 @property
1495 def station(self):
1496 return self.codes.station
1498 @property
1499 def location(self):
1500 return self.codes.location
1502 @property
1503 def channel(self):
1504 return self.codes.channel
1506 @property
1507 def extra(self):
1508 return self.codes.extra
1510 def overlaps(self, tmin, tmax):
1511 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1514def duration_to_str(t):
1515 if t > 24*3600:
1516 return '%gd' % (t / (24.*3600.))
1517 elif t > 3600:
1518 return '%gh' % (t / 3600.)
1519 elif t > 60:
1520 return '%gm' % (t / 60.)
1521 else:
1522 return '%gs' % t
1525class Coverage(Object):
1526 '''
1527 Information about times covered by a waveform or other time series data.
1528 '''
1529 kind_id = Int.T(
1530 help='Content type.')
1531 pattern = Codes.T(
1532 help='The codes pattern in the request, which caused this entry to '
1533 'match.')
1534 codes = Codes.T(
1535 help='NSLCE or NSL codes identifier of the time series.')
1536 deltat = Float.T(
1537 help='Sampling interval [s]',
1538 optional=True)
1539 tmin = Timestamp.T(
1540 help='Global start time of time series.',
1541 optional=True)
1542 tmax = Timestamp.T(
1543 help='Global end time of time series.',
1544 optional=True)
1545 changes = List.T(
1546 Tuple.T(2, Any.T()),
1547 help='List of change points, with entries of the form '
1548 '``(time, count)``, where a ``count`` of zero indicates start of '
1549 'a gap, a value of 1 start of normal data coverage and a higher '
1550 'value duplicate or redundant data coverage.')
1552 @classmethod
1553 def from_values(cls, args):
1554 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1555 return cls(
1556 kind_id=kind_id,
1557 pattern=pattern,
1558 codes=codes,
1559 deltat=deltat,
1560 tmin=tmin,
1561 tmax=tmax,
1562 changes=changes)
1564 @property
1565 def summary(self):
1566 ts = '%s - %s,' % (
1567 util.time_to_str(self.tmin),
1568 util.time_to_str(self.tmax))
1570 srate = self.sample_rate
1572 total = self.total
1574 return ' '.join((
1575 ('%s,' % to_kind(self.kind_id)).ljust(9),
1576 ('%s,' % str(self.codes)).ljust(18),
1577 ts,
1578 '%10.3g,' % srate if srate else '',
1579 '%4i' % len(self.changes),
1580 '%s' % duration_to_str(total) if total else 'none'))
1582 @property
1583 def sample_rate(self):
1584 if self.deltat is None:
1585 return None
1586 elif self.deltat == 0.0:
1587 return 0.0
1588 else:
1589 return 1.0 / self.deltat
1591 @property
1592 def labels(self):
1593 srate = self.sample_rate
1594 return (
1595 ('%s' % str(self.codes)),
1596 '%.3g' % srate if srate else '')
1598 @property
1599 def total(self):
1600 total_t = None
1601 for tmin, tmax, _ in self.iter_spans():
1602 total_t = (total_t or 0.0) + (tmax - tmin)
1604 return total_t
1606 def iter_spans(self):
1607 last = None
1608 for (t, count) in self.changes:
1609 if last is not None:
1610 last_t, last_count = last
1611 if last_count > 0:
1612 yield last_t, t, last_count
1614 last = (t, count)
1616 def iter_uncovered_by(self, other):
1617 a = self
1618 b = other
1619 ia = ib = -1
1620 ca = cb = 0
1621 last = None
1622 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1623 if ib + 1 == len(b.changes):
1624 ia += 1
1625 t, ca = a.changes[ia]
1626 elif ia + 1 == len(a.changes):
1627 ib += 1
1628 t, cb = b.changes[ib]
1629 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1630 ia += 1
1631 t, ca = a.changes[ia]
1632 else:
1633 ib += 1
1634 t, cb = b.changes[ib]
1636 if last is not None:
1637 tl, cal, cbl = last
1638 if tl < t and cal > 0 and cbl == 0:
1639 yield tl, t, ia, ib
1641 last = (t, ca, cb)
1643 def iter_uncovered_by_combined(self, other):
1644 ib_last = None
1645 group = None
1646 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1647 if ib_last is None or ib != ib_last:
1648 if group:
1649 yield (group[0][0], group[-1][1])
1651 group = []
1653 group.append((tmin, tmax))
1654 ib_last = ib
1656 if group:
1657 yield (group[0][0], group[-1][1])
1660__all__ = [
1661 'to_codes',
1662 'to_codes_guess',
1663 'to_codes_simple',
1664 'to_kind',
1665 'to_kinds',
1666 'to_kind_id',
1667 'to_kind_ids',
1668 'CodesError',
1669 'Codes',
1670 'CodesNSLCE',
1671 'CodesNSL',
1672 'CodesX',
1673 'Station',
1674 'Channel',
1675 'Sensor',
1676 'Response',
1677 'Nut',
1678 'Coverage',
1679 'WaveformPromise',
1680]