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_args(self):
655 def getattr_rep(k):
656 if k == 'codes':
657 return self.codes.replace(
658 channel=self.codes.channel[:-1] + '?')
659 else:
660 return getattr(self, k)
662 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
664 def _get_channel_args(self, component):
665 def getattr_rep(k):
666 if k == 'codes':
667 return self.codes.replace(
668 channel=self.codes.channel[:-1] + component)
669 else:
670 return getattr(self, k)
672 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
674 def _get_pyrocko_station_args(self):
675 return (
676 self.codes.network,
677 self.codes.station,
678 self.codes.location,
679 self.lat,
680 self.lon,
681 self.elevation,
682 self.depth,
683 self.north_shift,
684 self.east_shift)
687class Channel(ChannelBase):
688 '''
689 A channel of a seismic station.
690 '''
692 dip = Float.T(optional=True)
693 azimuth = Float.T(optional=True)
695 def get_pyrocko_channel(self):
696 from pyrocko import model
697 return model.Channel(*self._get_pyrocko_channel_args())
699 def _get_pyrocko_channel_args(self):
700 return (
701 self.codes.channel,
702 self.azimuth,
703 self.dip)
705 @property
706 def orientation_enz(self):
707 if None in (self.azimuth, self.dip):
708 return None
710 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
711 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
712 d = math.sin(self.dip*d2r)
713 return mkvec(e, n, -d)
716def cut_intervals(channels):
717 channels = list(channels)
718 times = set()
719 for channel in channels:
720 if channel.tmin is not None:
721 times.add(channel.tmin)
722 if channel.tmax is not None:
723 times.add(channel.tmax)
725 times = sorted(times)
727 if len(times) < 2:
728 return channels
730 channels_out = []
731 for channel in channels:
732 tstart = channel.tmin
733 for t in times:
734 if (channel.tmin is None or channel.tmin < t) \
735 and (channel.tmax is None or t < channel.tmax):
737 channel_out = clone(channel)
738 channel_out.tmin = tstart
739 channel_out.tmax = t
740 channels_out.append(channel_out)
741 tstart = t
743 if channel.tmax is None:
744 channel_out = clone(channel)
745 channel_out.tmin = tstart
746 channels_out.append(channel_out)
748 return channels_out
751class Sensor(ChannelBase):
752 '''
753 Representation of a channel group.
754 '''
756 channels = List.T(Channel.T())
758 @classmethod
759 def from_channels(cls, channels):
760 channels = cut_intervals(channels)
762 groups = defaultdict(list)
763 for channel in channels:
764 groups[channel._get_sensor_args()].append(channel)
766 return [
767 cls(channels=channels,
768 **dict(zip(ChannelBase.T.propnames, args)))
769 for args, channels in groups.items()]
771 def channel_vectors(self):
772 return num.vstack(
773 [channel.orientation_enz for channel in self.channels])
775 def projected_channels(self, system, component_names):
776 return [
777 Channel(
778 azimuth=math.atan2(e, n) * r2d,
779 dip=-math.asin(z) * r2d,
780 **dict(zip(
781 ChannelBase.T.propnames,
782 self._get_channel_args(comp))))
783 for comp, (e, n, z) in zip(component_names, system)]
785 def matrix_to(self, system, epsilon=1e-7):
786 m = num.dot(system, self.channel_vectors().T)
787 m[num.abs(m) < epsilon] = 0.0
788 return m
790 def projection_to(self, system, component_names):
791 return (
792 self.matrix_to(system),
793 self.channels,
794 self.projected_channels(system, component_names))
796 def projection_to_enz(self):
797 return self.projection_to(num.identity(3), 'ENZ')
799 def projection_to_trz(self, source, azimuth=None):
800 if azimuth is not None:
801 assert source is None
802 else:
803 azimuth = source.azibazi_to(self)[1] + 180.
805 return self.projection_to(num.array([
806 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
807 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
808 [0., 0., 1.]], dtype=float), 'TRZ')
810 def project_to_enz(self, traces):
811 from pyrocko import trace
813 matrix, in_channels, out_channels = self.projection_to_enz()
815 return trace.project(traces, matrix, in_channels, out_channels)
817 def project_to_trz(self, source, traces, azimuth=None):
818 from pyrocko import trace
820 matrix, in_channels, out_channels = self.projection_to_trz(
821 source, azimuth=azimuth)
823 return trace.project(traces, matrix, in_channels, out_channels)
826observational_quantities = [
827 'acceleration', 'velocity', 'displacement', 'pressure',
828 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
829 'temperature']
832technical_quantities = [
833 'voltage', 'counts']
836class QuantityType(StringChoice):
837 '''
838 Choice of observational or technical quantity.
840 SI units are used for all quantities, where applicable.
841 '''
842 choices = observational_quantities + technical_quantities
845class ResponseStage(Object):
846 '''
847 Representation of a response stage.
849 Components of a seismic recording system are represented as a sequence of
850 response stages, e.g. sensor, pre-amplifier, digitizer, digital
851 downsampling.
852 '''
853 input_quantity = QuantityType.T(optional=True)
854 input_sample_rate = Float.T(optional=True)
855 output_quantity = QuantityType.T(optional=True)
856 output_sample_rate = Float.T(optional=True)
857 elements = List.T(FrequencyResponse.T())
858 log = List.T(Tuple.T(3, String.T()))
860 @property
861 def stage_type(self):
862 if self.input_quantity in observational_quantities \
863 and self.output_quantity in observational_quantities:
864 return 'conversion'
866 if self.input_quantity in observational_quantities \
867 and self.output_quantity == 'voltage':
868 return 'sensor'
870 elif self.input_quantity == 'voltage' \
871 and self.output_quantity == 'voltage':
872 return 'electronics'
874 elif self.input_quantity == 'voltage' \
875 and self.output_quantity == 'counts':
876 return 'digitizer'
878 elif self.input_quantity == 'counts' \
879 and self.output_quantity == 'counts' \
880 and self.input_sample_rate != self.output_sample_rate:
881 return 'decimation'
883 elif self.input_quantity in observational_quantities \
884 and self.output_quantity == 'counts':
885 return 'combination'
887 else:
888 return 'unknown'
890 @property
891 def summary(self):
892 irate = self.input_sample_rate
893 orate = self.output_sample_rate
894 factor = None
895 if irate and orate:
896 factor = irate / orate
897 return 'ResponseStage, ' + (
898 '%s%s => %s%s%s' % (
899 self.input_quantity or '?',
900 ' @ %g Hz' % irate if irate else '',
901 self.output_quantity or '?',
902 ' @ %g Hz' % orate if orate else '',
903 ' :%g' % factor if factor else '')
904 )
906 def get_effective(self):
907 return MultiplyResponse(responses=list(self.elements))
910D = 'displacement'
911V = 'velocity'
912A = 'acceleration'
914g_converters = {
915 (V, D): IntegrationResponse(1),
916 (A, D): IntegrationResponse(2),
917 (D, V): DifferentiationResponse(1),
918 (A, V): IntegrationResponse(1),
919 (D, A): DifferentiationResponse(2),
920 (V, A): DifferentiationResponse(1)}
923def response_converters(input_quantity, output_quantity):
924 if input_quantity is None or input_quantity == output_quantity:
925 return []
927 if output_quantity is None:
928 raise ConversionError('Unspecified target quantity.')
930 try:
931 return [g_converters[input_quantity, output_quantity]]
933 except KeyError:
934 raise ConversionError('No rule to convert from "%s" to "%s".' % (
935 input_quantity, output_quantity))
938@squirrel_content
939class Response(Object):
940 '''
941 The instrument response of a seismic station channel.
942 '''
944 codes = CodesNSLCE.T()
945 tmin = Timestamp.T(optional=True)
946 tmax = Timestamp.T(optional=True)
948 stages = List.T(ResponseStage.T())
949 checkpoints = List.T(FrequencyResponseCheckpoint.T())
951 deltat = Float.T(optional=True)
952 log = List.T(Tuple.T(3, String.T()))
954 def __init__(self, **kwargs):
955 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
956 Object.__init__(self, **kwargs)
958 @property
959 def time_span(self):
960 return (self.tmin, self.tmax)
962 @property
963 def nstages(self):
964 return len(self.stages)
966 @property
967 def input_quantity(self):
968 return self.stages[0].input_quantity if self.stages else None
970 @property
971 def output_quantity(self):
972 return self.stages[-1].input_quantity if self.stages else None
974 @property
975 def output_sample_rate(self):
976 return self.stages[-1].output_sample_rate if self.stages else None
978 @property
979 def stages_summary(self):
980 def grouped(xs):
981 xs = list(xs)
982 g = []
983 for i in range(len(xs)):
984 g.append(xs[i])
985 if i+1 < len(xs) and xs[i+1] != xs[i]:
986 yield g
987 g = []
989 if g:
990 yield g
992 return '+'.join(
993 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
994 for g in grouped(stage.stage_type for stage in self.stages))
996 @property
997 def summary(self):
998 orate = self.output_sample_rate
999 return '%s %-16s %s' % (
1000 self.__class__.__name__, self.str_codes, self.str_time_span) \
1001 + ', ' + ', '.join((
1002 '%s => %s' % (
1003 self.input_quantity or '?', self.output_quantity or '?')
1004 + (' @ %g Hz' % orate if orate else ''),
1005 self.stages_summary,
1006 ))
1008 def get_effective(self, input_quantity=None):
1009 try:
1010 elements = response_converters(input_quantity, self.input_quantity)
1011 except ConversionError as e:
1012 raise ConversionError(str(e) + ' (%s)' % self.summary)
1014 elements.extend(
1015 stage.get_effective() for stage in self.stages)
1017 return MultiplyResponse(responses=simplify_responses(elements))
1020@squirrel_content
1021class Event(Object):
1022 '''
1023 A seismic event.
1024 '''
1026 name = String.T(optional=True)
1027 time = Timestamp.T()
1028 duration = Float.T(optional=True)
1030 lat = Float.T()
1031 lon = Float.T()
1032 depth = Float.T(optional=True)
1034 magnitude = Float.T(optional=True)
1036 def get_pyrocko_event(self):
1037 from pyrocko import model
1038 model.Event(
1039 name=self.name,
1040 time=self.time,
1041 lat=self.lat,
1042 lon=self.lon,
1043 depth=self.depth,
1044 magnitude=self.magnitude,
1045 duration=self.duration)
1047 @property
1048 def time_span(self):
1049 return (self.time, self.time)
1052def ehash(s):
1053 return hashlib.sha1(s.encode('utf8')).hexdigest()
1056def random_name(n=8):
1057 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1060@squirrel_content
1061class WaveformPromise(Object):
1062 '''
1063 Information about a waveform potentially downloadable from a remote site.
1065 In the Squirrel framework, waveform promises are used to permit download of
1066 selected waveforms from a remote site. They are typically generated by
1067 calls to
1068 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1069 Waveform promises are inserted and indexed in the database similar to
1070 normal waveforms. When processing a waveform query, e.g. from
1071 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1072 waveform is available for the queried time span, a matching promise can be
1073 resolved, i.e. an attempt is made to download the waveform from the remote
1074 site. The promise is removed after the download attempt (except when a
1075 network error occurs). This prevents Squirrel from making unnecessary
1076 queries for waveforms missing at the remote site.
1077 '''
1079 codes = CodesNSLCE.T()
1080 tmin = Timestamp.T()
1081 tmax = Timestamp.T()
1083 deltat = Float.T(optional=True)
1085 source_hash = String.T()
1087 def __init__(self, **kwargs):
1088 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1089 Object.__init__(self, **kwargs)
1091 @property
1092 def time_span(self):
1093 return (self.tmin, self.tmax)
1096class InvalidWaveform(Exception):
1097 pass
1100class WaveformOrder(Object):
1101 '''
1102 Waveform request information.
1103 '''
1105 source_id = String.T()
1106 codes = CodesNSLCE.T()
1107 deltat = Float.T()
1108 tmin = Timestamp.T()
1109 tmax = Timestamp.T()
1110 gaps = List.T(Tuple.T(2, Timestamp.T()))
1112 @property
1113 def client(self):
1114 return self.source_id.split(':')[1]
1116 def describe(self, site='?'):
1117 return '%s:%s %s [%s - %s]' % (
1118 self.client, site, str(self.codes),
1119 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1121 def validate(self, tr):
1122 if tr.ydata.size == 0:
1123 raise InvalidWaveform(
1124 'waveform with zero data samples')
1126 if tr.deltat != self.deltat:
1127 raise InvalidWaveform(
1128 'incorrect sampling interval - waveform: %g s, '
1129 'meta-data: %g s' % (
1130 tr.deltat, self.deltat))
1132 if not num.all(num.isfinite(tr.ydata)):
1133 raise InvalidWaveform('waveform has NaN values')
1136def order_summary(orders):
1137 codes_list = sorted(set(order.codes for order in orders))
1138 if len(codes_list) > 3:
1139 return '%i order%s: %s - %s' % (
1140 len(orders),
1141 util.plural_s(orders),
1142 str(codes_list[0]),
1143 str(codes_list[1]))
1145 else:
1146 return '%i order%s: %s' % (
1147 len(orders),
1148 util.plural_s(orders),
1149 ', '.join(str(codes) for codes in codes_list))
1152class Nut(Object):
1153 '''
1154 Index entry referencing an elementary piece of content.
1156 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1157 meta-information about individual pieces of waveforms, stations, channels,
1158 etc. together with the information where it was found or generated.
1159 '''
1161 file_path = String.T(optional=True)
1162 file_format = String.T(optional=True)
1163 file_mtime = Timestamp.T(optional=True)
1164 file_size = Int.T(optional=True)
1166 file_segment = Int.T(optional=True)
1167 file_element = Int.T(optional=True)
1169 kind_id = Int.T()
1170 codes = Codes.T()
1172 tmin_seconds = Int.T(default=0)
1173 tmin_offset = Int.T(default=0, optional=True)
1174 tmax_seconds = Int.T(default=0)
1175 tmax_offset = Int.T(default=0, optional=True)
1177 deltat = Float.T(default=0.0)
1179 content = Any.T(optional=True)
1180 raw_content = Dict.T(String.T(), Any.T())
1182 content_in_db = False
1184 def __init__(
1185 self,
1186 file_path=None,
1187 file_format=None,
1188 file_mtime=None,
1189 file_size=None,
1190 file_segment=None,
1191 file_element=None,
1192 kind_id=0,
1193 codes=CodesX(''),
1194 tmin_seconds=None,
1195 tmin_offset=0,
1196 tmax_seconds=None,
1197 tmax_offset=0,
1198 deltat=None,
1199 content=None,
1200 raw_content=None,
1201 tmin=None,
1202 tmax=None,
1203 values_nocheck=None):
1205 if values_nocheck is not None:
1206 (self.file_path, self.file_format, self.file_mtime,
1207 self.file_size,
1208 self.file_segment, self.file_element,
1209 self.kind_id, codes_safe_str,
1210 self.tmin_seconds, self.tmin_offset,
1211 self.tmax_seconds, self.tmax_offset,
1212 self.deltat) = values_nocheck
1214 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1215 self.deltat = self.deltat or None
1216 self.raw_content = {}
1217 self.content = None
1218 else:
1219 if tmin is not None:
1220 tmin_seconds, tmin_offset = tsplit(tmin)
1222 if tmax is not None:
1223 tmax_seconds, tmax_offset = tsplit(tmax)
1225 self.kind_id = int(kind_id)
1226 self.codes = codes
1227 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1228 self.tmin_offset = int(tmin_offset)
1229 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1230 self.tmax_offset = int(tmax_offset)
1231 self.deltat = float_or_none(deltat)
1232 self.file_path = str_or_none(file_path)
1233 self.file_segment = int_or_none(file_segment)
1234 self.file_element = int_or_none(file_element)
1235 self.file_format = str_or_none(file_format)
1236 self.file_mtime = float_or_none(file_mtime)
1237 self.file_size = int_or_none(file_size)
1238 self.content = content
1239 if raw_content is None:
1240 self.raw_content = {}
1241 else:
1242 self.raw_content = raw_content
1244 Object.__init__(self, init_props=False)
1246 def __eq__(self, other):
1247 return (isinstance(other, Nut) and
1248 self.equality_values == other.equality_values)
1250 def hash(self):
1251 return ehash(','.join(str(x) for x in self.key))
1253 def __ne__(self, other):
1254 return not (self == other)
1256 def get_io_backend(self):
1257 from . import io
1258 return io.get_backend(self.file_format)
1260 def file_modified(self):
1261 return self.get_io_backend().get_stats(self.file_path) \
1262 != (self.file_mtime, self.file_size)
1264 @property
1265 def dkey(self):
1266 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1268 @property
1269 def key(self):
1270 return (
1271 self.file_path,
1272 self.file_segment,
1273 self.file_element,
1274 self.file_mtime)
1276 @property
1277 def equality_values(self):
1278 return (
1279 self.file_segment, self.file_element,
1280 self.kind_id, self.codes,
1281 self.tmin_seconds, self.tmin_offset,
1282 self.tmax_seconds, self.tmax_offset, self.deltat)
1284 def diff(self, other):
1285 names = [
1286 'file_segment', 'file_element', 'kind_id', 'codes',
1287 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1288 'deltat']
1290 d = []
1291 for name, a, b in zip(
1292 names, self.equality_values, other.equality_values):
1294 if a != b:
1295 d.append((name, a, b))
1297 return d
1299 @property
1300 def tmin(self):
1301 return tjoin(self.tmin_seconds, self.tmin_offset)
1303 @tmin.setter
1304 def tmin(self, t):
1305 self.tmin_seconds, self.tmin_offset = tsplit(t)
1307 @property
1308 def tmax(self):
1309 return tjoin(self.tmax_seconds, self.tmax_offset)
1311 @tmax.setter
1312 def tmax(self, t):
1313 self.tmax_seconds, self.tmax_offset = tsplit(t)
1315 @property
1316 def kscale(self):
1317 if self.tmin_seconds is None or self.tmax_seconds is None:
1318 return 0
1319 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1321 @property
1322 def waveform_kwargs(self):
1323 network, station, location, channel, extra = self.codes
1325 return dict(
1326 network=network,
1327 station=station,
1328 location=location,
1329 channel=channel,
1330 extra=extra,
1331 tmin=self.tmin,
1332 tmax=self.tmax,
1333 deltat=self.deltat)
1335 @property
1336 def waveform_promise_kwargs(self):
1337 return dict(
1338 codes=self.codes,
1339 tmin=self.tmin,
1340 tmax=self.tmax,
1341 deltat=self.deltat)
1343 @property
1344 def station_kwargs(self):
1345 network, station, location = self.codes
1346 return dict(
1347 codes=self.codes,
1348 tmin=tmin_or_none(self.tmin),
1349 tmax=tmax_or_none(self.tmax))
1351 @property
1352 def channel_kwargs(self):
1353 network, station, location, channel, extra = self.codes
1354 return dict(
1355 codes=self.codes,
1356 tmin=tmin_or_none(self.tmin),
1357 tmax=tmax_or_none(self.tmax),
1358 deltat=self.deltat)
1360 @property
1361 def response_kwargs(self):
1362 return dict(
1363 codes=self.codes,
1364 tmin=tmin_or_none(self.tmin),
1365 tmax=tmax_or_none(self.tmax),
1366 deltat=self.deltat)
1368 @property
1369 def event_kwargs(self):
1370 return dict(
1371 name=self.codes,
1372 time=self.tmin,
1373 duration=(self.tmax - self.tmin) or None)
1375 @property
1376 def trace_kwargs(self):
1377 network, station, location, channel, extra = self.codes
1379 return dict(
1380 network=network,
1381 station=station,
1382 location=location,
1383 channel=channel,
1384 extra=extra,
1385 tmin=self.tmin,
1386 tmax=self.tmax-self.deltat,
1387 deltat=self.deltat)
1389 @property
1390 def dummy_trace(self):
1391 return DummyTrace(self)
1393 @property
1394 def summary(self):
1395 if self.tmin == self.tmax:
1396 ts = util.time_to_str(self.tmin)
1397 else:
1398 ts = '%s - %s' % (
1399 util.time_to_str(self.tmin),
1400 util.time_to_str(self.tmax))
1402 return ' '.join((
1403 ('%s,' % to_kind(self.kind_id)).ljust(9),
1404 ('%s,' % str(self.codes)).ljust(18),
1405 ts))
1408def make_waveform_nut(**kwargs):
1409 return Nut(kind_id=WAVEFORM, **kwargs)
1412def make_waveform_promise_nut(**kwargs):
1413 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1416def make_station_nut(**kwargs):
1417 return Nut(kind_id=STATION, **kwargs)
1420def make_channel_nut(**kwargs):
1421 return Nut(kind_id=CHANNEL, **kwargs)
1424def make_response_nut(**kwargs):
1425 return Nut(kind_id=RESPONSE, **kwargs)
1428def make_event_nut(**kwargs):
1429 return Nut(kind_id=EVENT, **kwargs)
1432def group_channels(nuts):
1433 by_ansl = {}
1434 for nut in nuts:
1435 if nut.kind_id != CHANNEL:
1436 continue
1438 ansl = nut.codes[:4]
1440 if ansl not in by_ansl:
1441 by_ansl[ansl] = {}
1443 group = by_ansl[ansl]
1445 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1447 if k not in group:
1448 group[k] = set()
1450 group.add(nut.codes[4])
1452 return by_ansl
1455class DummyTrace(object):
1457 def __init__(self, nut):
1458 self.nut = nut
1459 self.codes = nut.codes
1460 self.meta = {}
1462 @property
1463 def tmin(self):
1464 return self.nut.tmin
1466 @property
1467 def tmax(self):
1468 return self.nut.tmax
1470 @property
1471 def deltat(self):
1472 return self.nut.deltat
1474 @property
1475 def nslc_id(self):
1476 return self.codes.nslc
1478 @property
1479 def network(self):
1480 return self.codes.network
1482 @property
1483 def station(self):
1484 return self.codes.station
1486 @property
1487 def location(self):
1488 return self.codes.location
1490 @property
1491 def channel(self):
1492 return self.codes.channel
1494 @property
1495 def extra(self):
1496 return self.codes.extra
1498 def overlaps(self, tmin, tmax):
1499 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1502def duration_to_str(t):
1503 if t > 24*3600:
1504 return '%gd' % (t / (24.*3600.))
1505 elif t > 3600:
1506 return '%gh' % (t / 3600.)
1507 elif t > 60:
1508 return '%gm' % (t / 60.)
1509 else:
1510 return '%gs' % t
1513class Coverage(Object):
1514 '''
1515 Information about times covered by a waveform or other time series data.
1516 '''
1517 kind_id = Int.T(
1518 help='Content type.')
1519 pattern = Codes.T(
1520 help='The codes pattern in the request, which caused this entry to '
1521 'match.')
1522 codes = Codes.T(
1523 help='NSLCE or NSL codes identifier of the time series.')
1524 deltat = Float.T(
1525 help='Sampling interval [s]',
1526 optional=True)
1527 tmin = Timestamp.T(
1528 help='Global start time of time series.',
1529 optional=True)
1530 tmax = Timestamp.T(
1531 help='Global end time of time series.',
1532 optional=True)
1533 changes = List.T(
1534 Tuple.T(2, Any.T()),
1535 help='List of change points, with entries of the form '
1536 '``(time, count)``, where a ``count`` of zero indicates start of '
1537 'a gap, a value of 1 start of normal data coverage and a higher '
1538 'value duplicate or redundant data coverage.')
1540 @classmethod
1541 def from_values(cls, args):
1542 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1543 return cls(
1544 kind_id=kind_id,
1545 pattern=pattern,
1546 codes=codes,
1547 deltat=deltat,
1548 tmin=tmin,
1549 tmax=tmax,
1550 changes=changes)
1552 @property
1553 def summary(self):
1554 ts = '%s - %s,' % (
1555 util.time_to_str(self.tmin),
1556 util.time_to_str(self.tmax))
1558 srate = self.sample_rate
1560 return ' '.join((
1561 ('%s,' % to_kind(self.kind_id)).ljust(9),
1562 ('%s,' % str(self.codes)).ljust(18),
1563 ts,
1564 '%10.3g,' % srate if srate else '',
1565 '%4i' % len(self.changes),
1566 '%s' % duration_to_str(self.total)))
1568 @property
1569 def sample_rate(self):
1570 if self.deltat is None:
1571 return None
1572 elif self.deltat == 0.0:
1573 return 0.0
1574 else:
1575 return 1.0 / self.deltat
1577 @property
1578 def labels(self):
1579 srate = self.sample_rate
1580 return (
1581 ('%s' % str(self.codes)),
1582 '%.3g' % srate if srate else '')
1584 @property
1585 def total(self):
1586 total_t = None
1587 for tmin, tmax, _ in self.iter_spans():
1588 total_t = (total_t or 0.0) + (tmax - tmin)
1590 return total_t
1592 def iter_spans(self):
1593 last = None
1594 for (t, count) in self.changes:
1595 if last is not None:
1596 last_t, last_count = last
1597 if last_count > 0:
1598 yield last_t, t, last_count
1600 last = (t, count)
1603__all__ = [
1604 'to_codes',
1605 'to_codes_guess',
1606 'to_codes_simple',
1607 'to_kind',
1608 'to_kinds',
1609 'to_kind_id',
1610 'to_kind_ids',
1611 'CodesError',
1612 'Codes',
1613 'CodesNSLCE',
1614 'CodesNSL',
1615 'CodesX',
1616 'Station',
1617 'Channel',
1618 'Sensor',
1619 'Response',
1620 'Nut',
1621 'Coverage',
1622 'WaveformPromise',
1623]