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[:5]
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 @property
188 def codes_nsl(self):
189 return CodesNSL(self)
191 @property
192 def codes_nsl_star(self):
193 return CodesNSL(self.network, self.station, '*')
195 def as_tuple(self):
196 return tuple(self)
198 def replace(self, **kwargs):
199 x = CodesNSLCEBase._replace(self, **kwargs)
200 return g_codes_pool.setdefault(x, x)
203def normalize_nsl(*args, **kwargs):
204 if args and kwargs:
205 raise ValueError('Either *args or **kwargs accepted, not both.')
207 if len(args) == 1:
208 if isinstance(args[0], str):
209 args = tuple(args[0].split('.'))
210 elif isinstance(args[0], tuple):
211 args = args[0]
212 else:
213 raise ValueError('Invalid argument type: %s' % type(args[0]))
215 nargs = len(args)
216 if nargs == 3:
217 t = args
219 elif nargs == 0:
220 d = dict(
221 network='',
222 station='',
223 location='')
225 d.update(kwargs)
226 t = tuple(kwargs.get(k, '') for k in (
227 'network', 'station', 'location'))
229 else:
230 raise CodesError(
231 'Does not match NSL codes pattern: %s' % '.'.join(args))
233 if '.'.join(t).count('.') != 2:
234 raise CodesError(
235 'Codes may not contain a ".": "%s", "%s", "%s"' % t)
237 return t
240CodesNSLBase = namedtuple(
241 'CodesNSLBase', [
242 'network', 'station', 'location'])
245class CodesNSL(CodesNSLBase, Codes):
246 '''
247 Codes denominating a seismic station (NSL).
249 NET.STA.LOC is accepted, slightly different from SEED/StationXML, where
250 LOC is part of the channel. By setting location='*' is possible to get
251 compatible behaviour in most cases.
252 '''
254 __slots__ = ()
255 __hash__ = CodesNSLBase.__hash__
257 as_dict = CodesNSLBase._asdict
259 def __new__(cls, *args, safe_str=None, **kwargs):
260 nargs = len(args)
261 if nargs == 1 and isinstance(args[0], CodesNSL):
262 return args[0]
263 elif nargs == 1 and isinstance(args[0], CodesNSLCE):
264 t = args[0].nsl
265 elif nargs == 1 and isinstance(args[0], CodesX):
266 t = ('*', '*', '*')
267 elif safe_str is not None:
268 t = safe_str.split('.')
269 else:
270 t = normalize_nsl(*args, **kwargs)
272 x = CodesNSLBase.__new__(cls, *t)
273 return g_codes_pool.setdefault(x, x)
275 def __init__(self, *args, **kwargs):
276 Codes.__init__(self)
278 def __str__(self):
279 return '.'.join(self)
281 def __eq__(self, other):
282 if not isinstance(other, CodesNSL):
283 other = CodesNSL(other)
285 return CodesNSLBase.__eq__(self, other)
287 def matches(self, pattern):
288 if not isinstance(pattern, CodesNSL):
289 pattern = CodesNSL(pattern)
291 return match_codes(pattern, self)
293 @property
294 def safe_str(self):
295 return '.'.join(self)
297 @property
298 def ns(self):
299 return self[:2]
301 @property
302 def nsl(self):
303 return self[:3]
305 def as_tuple(self):
306 return tuple(self)
308 def replace(self, **kwargs):
309 x = CodesNSLBase._replace(self, **kwargs)
310 return g_codes_pool.setdefault(x, x)
313CodesXBase = namedtuple(
314 'CodesXBase', [
315 'name'])
318class CodesX(CodesXBase, Codes):
319 '''
320 General purpose codes for anything other than channels or stations.
321 '''
323 __slots__ = ()
324 __hash__ = CodesXBase.__hash__
325 __eq__ = CodesXBase.__eq__
327 as_dict = CodesXBase._asdict
329 def __new__(cls, name='', safe_str=None):
330 if isinstance(name, CodesX):
331 return name
332 elif isinstance(name, (CodesNSLCE, CodesNSL)):
333 name = '*'
334 elif safe_str is not None:
335 name = safe_str
336 else:
337 if '.' in name:
338 raise CodesError('Code may not contain a ".": %s' % name)
340 x = CodesXBase.__new__(cls, name)
341 return g_codes_pool.setdefault(x, x)
343 def __init__(self, *args, **kwargs):
344 Codes.__init__(self)
346 def __str__(self):
347 return '.'.join(self)
349 @property
350 def safe_str(self):
351 return '.'.join(self)
353 @property
354 def ns(self):
355 return self[:2]
357 def as_tuple(self):
358 return tuple(self)
360 def replace(self, **kwargs):
361 x = CodesXBase._replace(self, **kwargs)
362 return g_codes_pool.setdefault(x, x)
365g_codes_patterns = {}
368def match_codes(pattern, codes):
369 spattern = pattern.safe_str
370 scodes = codes.safe_str
371 if spattern not in g_codes_patterns:
372 rpattern = re.compile(fnmatch.translate(spattern), re.I)
373 g_codes_patterns[spattern] = rpattern
375 rpattern = g_codes_patterns[spattern]
376 return bool(rpattern.match(scodes))
379g_content_kinds = [
380 'undefined',
381 'waveform',
382 'station',
383 'channel',
384 'response',
385 'event',
386 'waveform_promise']
389g_codes_classes = [
390 CodesX,
391 CodesNSLCE,
392 CodesNSL,
393 CodesNSLCE,
394 CodesNSLCE,
395 CodesX,
396 CodesNSLCE]
398g_codes_classes_ndot = {
399 0: CodesX,
400 2: CodesNSL,
401 3: CodesNSLCE,
402 4: CodesNSLCE}
405def to_codes_simple(kind_id, codes_safe_str):
406 return g_codes_classes[kind_id](safe_str=codes_safe_str)
409def to_codes(kind_id, obj):
410 return g_codes_classes[kind_id](obj)
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)
420# derived list class to enable detection of already preprocessed codes patterns
421class codes_patterns_list(list):
422 pass
425def codes_patterns_for_kind(kind_id, codes):
426 if isinstance(codes, codes_patterns_list):
427 return codes
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))
434 return lcodes
436 codes = to_codes(kind_id, codes)
438 lcodes = codes_patterns_list()
439 lcodes.append(codes)
440 if kind_id == STATION:
441 lcodes.append(codes.replace(location='[*]'))
443 return lcodes
446g_content_kind_ids = (
447 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT,
448 WAVEFORM_PROMISE) = range(len(g_content_kinds))
451g_tmin, g_tmax = util.get_working_system_time_range()[:2]
454try:
455 g_tmin_queries = max(g_tmin, util.str_to_time_fillup('1900-01-01'))
456except Exception:
457 g_tmin_queries = g_tmin
460def to_kind(kind_id):
461 return g_content_kinds[kind_id]
464def to_kinds(kind_ids):
465 return [g_content_kinds[kind_id] for kind_id in kind_ids]
468def to_kind_id(kind):
469 return g_content_kinds.index(kind)
472def to_kind_ids(kinds):
473 return [g_content_kinds.index(kind) for kind in kinds]
476g_kind_mask_all = 0xff
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
486def str_or_none(x):
487 if x is None:
488 return None
489 else:
490 return str(x)
493def float_or_none(x):
494 if x is None:
495 return None
496 else:
497 return float(x)
500def int_or_none(x):
501 if x is None:
502 return None
503 else:
504 return int(x)
507def int_or_g_tmin(x):
508 if x is None:
509 return g_tmin
510 else:
511 return int(x)
514def int_or_g_tmax(x):
515 if x is None:
516 return g_tmax
517 else:
518 return int(x)
521def tmin_or_none(x):
522 if x == g_tmin:
523 return None
524 else:
525 return x
528def tmax_or_none(x):
529 if x == g_tmax:
530 return None
531 else:
532 return x
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)
542def codes_to_str_abbreviated(codes, indent=' '):
543 codes = [str(x) for x in codes]
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>'
553 return scodes
556g_offset_time_unit_inv = 1000000000
557g_offset_time_unit = 1.0 / g_offset_time_unit_inv
560def tsplit(t):
561 if t is None:
562 return None, 0
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)
570 seconds = num.floor(t)
571 offset = t - seconds
572 return int(seconds), int(round(offset * g_offset_time_unit_inv))
575def tjoin(seconds, offset):
576 if seconds is None:
577 return None
579 return util.to_time_float(seconds) \
580 + util.to_time_float(offset*g_offset_time_unit)
583tscale_min = 1
584tscale_max = 365 * 24 * 3600 # last edge is one above
585tscale_logbase = 20
587tscale_edges = [tscale_min]
588while True:
589 tscale_edges.append(tscale_edges[-1]*tscale_logbase)
590 if tscale_edges[-1] >= tscale_max:
591 break
594tscale_edges = num.array(tscale_edges)
597def tscale_to_kscale(tscale):
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)
604 return int(num.searchsorted(tscale_edges, tscale))
607@squirrel_content
608class Station(Location):
609 '''
610 A seismic station.
611 '''
613 codes = CodesNSL.T()
615 tmin = Timestamp.T(optional=True)
616 tmax = Timestamp.T(optional=True)
618 description = Unicode.T(optional=True)
620 def __init__(self, **kwargs):
621 kwargs['codes'] = CodesNSL(kwargs['codes'])
622 Location.__init__(self, **kwargs)
624 @property
625 def time_span(self):
626 return (self.tmin, self.tmax)
628 def get_pyrocko_station(self):
629 from pyrocko import model
630 return model.Station(*self._get_pyrocko_station_args())
632 def _get_pyrocko_station_args(self):
633 return (
634 self.codes.network,
635 self.codes.station,
636 self.codes.location,
637 self.lat,
638 self.lon,
639 self.elevation,
640 self.depth,
641 self.north_shift,
642 self.east_shift)
645@squirrel_content
646class ChannelBase(Location):
647 codes = CodesNSLCE.T()
649 tmin = Timestamp.T(optional=True)
650 tmax = Timestamp.T(optional=True)
652 deltat = Float.T(optional=True)
654 def __init__(self, **kwargs):
655 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
656 Location.__init__(self, **kwargs)
658 @property
659 def time_span(self):
660 return (self.tmin, self.tmax)
662 def _get_sensor_codes(self):
663 return self.codes.replace(
664 channel=self.codes.channel[:-1] + '?')
666 def _get_sensor_args(self):
667 def getattr_rep(k):
668 if k == 'codes':
669 return self._get_sensor_codes()
670 else:
671 return getattr(self, k)
673 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
675 def _get_channel_args(self, component):
676 def getattr_rep(k):
677 if k == 'codes':
678 return self.codes.replace(
679 channel=self.codes.channel[:-1] + component)
680 else:
681 return getattr(self, k)
683 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
685 def _get_pyrocko_station_args(self):
686 return (
687 self.codes.network,
688 self.codes.station,
689 self.codes.location,
690 self.lat,
691 self.lon,
692 self.elevation,
693 self.depth,
694 self.north_shift,
695 self.east_shift)
698class Channel(ChannelBase):
699 '''
700 A channel of a seismic station.
701 '''
703 dip = Float.T(optional=True)
704 azimuth = Float.T(optional=True)
706 def get_pyrocko_channel(self):
707 from pyrocko import model
708 return model.Channel(*self._get_pyrocko_channel_args())
710 def _get_pyrocko_channel_args(self):
711 return (
712 self.codes.channel,
713 self.azimuth,
714 self.dip)
716 @property
717 def orientation_enz(self):
718 if None in (self.azimuth, self.dip):
719 return None
721 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
722 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
723 d = math.sin(self.dip*d2r)
724 return mkvec(e, n, -d)
727def cut_intervals(channels):
728 channels = list(channels)
729 times = set()
730 for channel in channels:
731 if channel.tmin is not None:
732 times.add(channel.tmin)
733 if channel.tmax is not None:
734 times.add(channel.tmax)
736 times = sorted(times)
738 if any(channel.tmin is None for channel in channels):
739 times[0:0] = [None]
741 if any(channel.tmax is None for channel in channels):
742 times.append(None)
744 if len(times) <= 2:
745 return channels
747 channels_out = []
748 for channel in channels:
749 for i in range(len(times)-1):
750 tmin = times[i]
751 tmax = times[i+1]
752 if ((channel.tmin is None or (
753 tmin is not None and channel.tmin <= tmin))
754 and (channel.tmax is None or (
755 tmax is not None and tmax <= channel.tmax))):
757 channel_out = clone(channel)
758 channel_out.tmin = tmin
759 channel_out.tmax = tmax
760 channels_out.append(channel_out)
762 return channels_out
765class Sensor(ChannelBase):
766 '''
767 Representation of a channel group.
768 '''
770 channels = List.T(Channel.T())
772 @classmethod
773 def from_channels(cls, channels):
774 groups = defaultdict(list)
775 for channel in channels:
776 groups[channel._get_sensor_codes()].append(channel)
778 channels_cut = []
779 for group in groups.values():
780 channels_cut.extend(cut_intervals(group))
782 groups = defaultdict(list)
783 for channel in channels_cut:
784 groups[channel._get_sensor_args()].append(channel)
786 return [
787 cls(channels=channels,
788 **dict(zip(ChannelBase.T.propnames, args)))
789 for args, channels in groups.items()]
791 def channel_vectors(self):
792 return num.vstack(
793 [channel.orientation_enz for channel in self.channels])
795 def projected_channels(self, system, component_names):
796 return [
797 Channel(
798 azimuth=math.atan2(e, n) * r2d,
799 dip=-math.asin(z) * r2d,
800 **dict(zip(
801 ChannelBase.T.propnames,
802 self._get_channel_args(comp))))
803 for comp, (e, n, z) in zip(component_names, system)]
805 def matrix_to(self, system, epsilon=1e-7):
806 m = num.dot(system, self.channel_vectors().T)
807 m[num.abs(m) < epsilon] = 0.0
808 return m
810 def projection_to(self, system, component_names):
811 return (
812 self.matrix_to(system),
813 self.channels,
814 self.projected_channels(system, component_names))
816 def projection_to_enz(self):
817 return self.projection_to(num.identity(3), 'ENZ')
819 def projection_to_trz(self, source, azimuth=None):
820 if azimuth is not None:
821 assert source is None
822 else:
823 azimuth = source.azibazi_to(self)[1] + 180.
825 return self.projection_to(num.array([
826 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
827 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
828 [0., 0., 1.]], dtype=float), 'TRZ')
830 def project_to_enz(self, traces):
831 from pyrocko import trace
833 matrix, in_channels, out_channels = self.projection_to_enz()
835 return trace.project(traces, matrix, in_channels, out_channels)
837 def project_to_trz(self, source, traces, azimuth=None):
838 from pyrocko import trace
840 matrix, in_channels, out_channels = self.projection_to_trz(
841 source, azimuth=azimuth)
843 return trace.project(traces, matrix, in_channels, out_channels)
846observational_quantities = [
847 'acceleration', 'velocity', 'displacement', 'pressure',
848 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
849 'temperature']
852technical_quantities = [
853 'voltage', 'counts']
856class QuantityType(StringChoice):
857 '''
858 Choice of observational or technical quantity.
860 SI units are used for all quantities, where applicable.
861 '''
862 choices = observational_quantities + technical_quantities
865class ResponseStage(Object):
866 '''
867 Representation of a response stage.
869 Components of a seismic recording system are represented as a sequence of
870 response stages, e.g. sensor, pre-amplifier, digitizer, digital
871 downsampling.
872 '''
873 input_quantity = QuantityType.T(optional=True)
874 input_sample_rate = Float.T(optional=True)
875 output_quantity = QuantityType.T(optional=True)
876 output_sample_rate = Float.T(optional=True)
877 elements = List.T(FrequencyResponse.T())
878 log = List.T(Tuple.T(3, String.T()))
880 @property
881 def stage_type(self):
882 if self.input_quantity in observational_quantities \
883 and self.output_quantity in observational_quantities:
884 return 'conversion'
886 if self.input_quantity in observational_quantities \
887 and self.output_quantity == 'voltage':
888 return 'sensor'
890 elif self.input_quantity == 'voltage' \
891 and self.output_quantity == 'voltage':
892 return 'electronics'
894 elif self.input_quantity == 'voltage' \
895 and self.output_quantity == 'counts':
896 return 'digitizer'
898 elif self.input_quantity == 'counts' \
899 and self.output_quantity == 'counts' \
900 and self.input_sample_rate != self.output_sample_rate:
901 return 'decimation'
903 elif self.input_quantity in observational_quantities \
904 and self.output_quantity == 'counts':
905 return 'combination'
907 else:
908 return 'unknown'
910 @property
911 def summary(self):
912 irate = self.input_sample_rate
913 orate = self.output_sample_rate
914 factor = None
915 if irate and orate:
916 factor = irate / orate
917 return 'ResponseStage, ' + (
918 '%s%s => %s%s%s' % (
919 self.input_quantity or '?',
920 ' @ %g Hz' % irate if irate else '',
921 self.output_quantity or '?',
922 ' @ %g Hz' % orate if orate else '',
923 ' :%g' % factor if factor else '')
924 )
926 def get_effective(self):
927 return MultiplyResponse(responses=list(self.elements))
930D = 'displacement'
931V = 'velocity'
932A = 'acceleration'
934g_converters = {
935 (V, D): IntegrationResponse(1),
936 (A, D): IntegrationResponse(2),
937 (D, V): DifferentiationResponse(1),
938 (A, V): IntegrationResponse(1),
939 (D, A): DifferentiationResponse(2),
940 (V, A): DifferentiationResponse(1)}
943def response_converters(input_quantity, output_quantity):
944 if input_quantity is None or input_quantity == output_quantity:
945 return []
947 if output_quantity is None:
948 raise ConversionError('Unspecified target quantity.')
950 try:
951 return [g_converters[input_quantity, output_quantity]]
953 except KeyError:
954 raise ConversionError('No rule to convert from "%s" to "%s".' % (
955 input_quantity, output_quantity))
958@squirrel_content
959class Response(Object):
960 '''
961 The instrument response of a seismic station channel.
962 '''
964 codes = CodesNSLCE.T()
965 tmin = Timestamp.T(optional=True)
966 tmax = Timestamp.T(optional=True)
968 stages = List.T(ResponseStage.T())
969 checkpoints = List.T(FrequencyResponseCheckpoint.T())
971 deltat = Float.T(optional=True)
972 log = List.T(Tuple.T(3, String.T()))
974 def __init__(self, **kwargs):
975 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
976 Object.__init__(self, **kwargs)
978 @property
979 def time_span(self):
980 return (self.tmin, self.tmax)
982 @property
983 def nstages(self):
984 return len(self.stages)
986 @property
987 def input_quantity(self):
988 return self.stages[0].input_quantity if self.stages else None
990 @property
991 def output_quantity(self):
992 return self.stages[-1].input_quantity if self.stages else None
994 @property
995 def output_sample_rate(self):
996 return self.stages[-1].output_sample_rate if self.stages else None
998 @property
999 def stages_summary(self):
1000 def grouped(xs):
1001 xs = list(xs)
1002 g = []
1003 for i in range(len(xs)):
1004 g.append(xs[i])
1005 if i+1 < len(xs) and xs[i+1] != xs[i]:
1006 yield g
1007 g = []
1009 if g:
1010 yield g
1012 return '+'.join(
1013 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1014 for g in grouped(stage.stage_type for stage in self.stages))
1016 @property
1017 def summary(self):
1018 orate = self.output_sample_rate
1019 return '%s %-16s %s' % (
1020 self.__class__.__name__, self.str_codes, self.str_time_span) \
1021 + ', ' + ', '.join((
1022 '%s => %s' % (
1023 self.input_quantity or '?', self.output_quantity or '?')
1024 + (' @ %g Hz' % orate if orate else ''),
1025 self.stages_summary,
1026 ))
1028 def get_effective(self, input_quantity=None):
1029 try:
1030 elements = response_converters(input_quantity, self.input_quantity)
1031 except ConversionError as e:
1032 raise ConversionError(str(e) + ' (%s)' % self.summary)
1034 elements.extend(
1035 stage.get_effective() for stage in self.stages)
1037 return MultiplyResponse(responses=simplify_responses(elements))
1040@squirrel_content
1041class Event(Object):
1042 '''
1043 A seismic event.
1044 '''
1046 name = String.T(optional=True)
1047 time = Timestamp.T()
1048 duration = Float.T(optional=True)
1050 lat = Float.T()
1051 lon = Float.T()
1052 depth = Float.T(optional=True)
1054 magnitude = Float.T(optional=True)
1056 def get_pyrocko_event(self):
1057 from pyrocko import model
1058 model.Event(
1059 name=self.name,
1060 time=self.time,
1061 lat=self.lat,
1062 lon=self.lon,
1063 depth=self.depth,
1064 magnitude=self.magnitude,
1065 duration=self.duration)
1067 @property
1068 def time_span(self):
1069 return (self.time, self.time)
1072def ehash(s):
1073 return hashlib.sha1(s.encode('utf8')).hexdigest()
1076def random_name(n=8):
1077 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1080@squirrel_content
1081class WaveformPromise(Object):
1082 '''
1083 Information about a waveform potentially downloadable from a remote site.
1085 In the Squirrel framework, waveform promises are used to permit download of
1086 selected waveforms from a remote site. They are typically generated by
1087 calls to
1088 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1089 Waveform promises are inserted and indexed in the database similar to
1090 normal waveforms. When processing a waveform query, e.g. from
1091 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1092 waveform is available for the queried time span, a matching promise can be
1093 resolved, i.e. an attempt is made to download the waveform from the remote
1094 site. The promise is removed after the download attempt (except when a
1095 network error occurs). This prevents Squirrel from making unnecessary
1096 queries for waveforms missing at the remote site.
1097 '''
1099 codes = CodesNSLCE.T()
1100 tmin = Timestamp.T()
1101 tmax = Timestamp.T()
1103 deltat = Float.T(optional=True)
1105 source_hash = String.T()
1107 def __init__(self, **kwargs):
1108 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1109 Object.__init__(self, **kwargs)
1111 @property
1112 def time_span(self):
1113 return (self.tmin, self.tmax)
1116class InvalidWaveform(Exception):
1117 pass
1120class WaveformOrder(Object):
1121 '''
1122 Waveform request information.
1123 '''
1125 source_id = String.T()
1126 codes = CodesNSLCE.T()
1127 deltat = Float.T()
1128 tmin = Timestamp.T()
1129 tmax = Timestamp.T()
1130 gaps = List.T(Tuple.T(2, Timestamp.T()))
1132 @property
1133 def client(self):
1134 return self.source_id.split(':')[1]
1136 def describe(self, site='?'):
1137 return '%s:%s %s [%s - %s]' % (
1138 self.client, site, str(self.codes),
1139 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1141 def validate(self, tr):
1142 if tr.ydata.size == 0:
1143 raise InvalidWaveform(
1144 'waveform with zero data samples')
1146 if tr.deltat != self.deltat:
1147 raise InvalidWaveform(
1148 'incorrect sampling interval - waveform: %g s, '
1149 'meta-data: %g s' % (
1150 tr.deltat, self.deltat))
1152 if not num.all(num.isfinite(tr.ydata)):
1153 raise InvalidWaveform('waveform has NaN values')
1156def order_summary(orders):
1157 codes_list = sorted(set(order.codes for order in orders))
1158 if len(codes_list) > 3:
1159 return '%i order%s: %s - %s' % (
1160 len(orders),
1161 util.plural_s(orders),
1162 str(codes_list[0]),
1163 str(codes_list[1]))
1165 else:
1166 return '%i order%s: %s' % (
1167 len(orders),
1168 util.plural_s(orders),
1169 ', '.join(str(codes) for codes in codes_list))
1172class Nut(Object):
1173 '''
1174 Index entry referencing an elementary piece of content.
1176 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1177 meta-information about individual pieces of waveforms, stations, channels,
1178 etc. together with the information where it was found or generated.
1179 '''
1181 file_path = String.T(optional=True)
1182 file_format = String.T(optional=True)
1183 file_mtime = Timestamp.T(optional=True)
1184 file_size = Int.T(optional=True)
1186 file_segment = Int.T(optional=True)
1187 file_element = Int.T(optional=True)
1189 kind_id = Int.T()
1190 codes = Codes.T()
1192 tmin_seconds = Int.T(default=0)
1193 tmin_offset = Int.T(default=0, optional=True)
1194 tmax_seconds = Int.T(default=0)
1195 tmax_offset = Int.T(default=0, optional=True)
1197 deltat = Float.T(default=0.0)
1199 content = Any.T(optional=True)
1200 raw_content = Dict.T(String.T(), Any.T())
1202 content_in_db = False
1204 def __init__(
1205 self,
1206 file_path=None,
1207 file_format=None,
1208 file_mtime=None,
1209 file_size=None,
1210 file_segment=None,
1211 file_element=None,
1212 kind_id=0,
1213 codes=CodesX(''),
1214 tmin_seconds=None,
1215 tmin_offset=0,
1216 tmax_seconds=None,
1217 tmax_offset=0,
1218 deltat=None,
1219 content=None,
1220 raw_content=None,
1221 tmin=None,
1222 tmax=None,
1223 values_nocheck=None):
1225 if values_nocheck is not None:
1226 (self.file_path, self.file_format, self.file_mtime,
1227 self.file_size,
1228 self.file_segment, self.file_element,
1229 self.kind_id, codes_safe_str,
1230 self.tmin_seconds, self.tmin_offset,
1231 self.tmax_seconds, self.tmax_offset,
1232 self.deltat) = values_nocheck
1234 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1235 self.deltat = self.deltat or None
1236 self.raw_content = {}
1237 self.content = None
1238 else:
1239 if tmin is not None:
1240 tmin_seconds, tmin_offset = tsplit(tmin)
1242 if tmax is not None:
1243 tmax_seconds, tmax_offset = tsplit(tmax)
1245 self.kind_id = int(kind_id)
1246 self.codes = codes
1247 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1248 self.tmin_offset = int(tmin_offset)
1249 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1250 self.tmax_offset = int(tmax_offset)
1251 self.deltat = float_or_none(deltat)
1252 self.file_path = str_or_none(file_path)
1253 self.file_segment = int_or_none(file_segment)
1254 self.file_element = int_or_none(file_element)
1255 self.file_format = str_or_none(file_format)
1256 self.file_mtime = float_or_none(file_mtime)
1257 self.file_size = int_or_none(file_size)
1258 self.content = content
1259 if raw_content is None:
1260 self.raw_content = {}
1261 else:
1262 self.raw_content = raw_content
1264 Object.__init__(self, init_props=False)
1266 def __eq__(self, other):
1267 return (isinstance(other, Nut) and
1268 self.equality_values == other.equality_values)
1270 def hash(self):
1271 return ehash(','.join(str(x) for x in self.key))
1273 def __ne__(self, other):
1274 return not (self == other)
1276 def get_io_backend(self):
1277 from . import io
1278 return io.get_backend(self.file_format)
1280 def file_modified(self):
1281 return self.get_io_backend().get_stats(self.file_path) \
1282 != (self.file_mtime, self.file_size)
1284 @property
1285 def dkey(self):
1286 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1288 @property
1289 def key(self):
1290 return (
1291 self.file_path,
1292 self.file_segment,
1293 self.file_element,
1294 self.file_mtime)
1296 @property
1297 def equality_values(self):
1298 return (
1299 self.file_segment, self.file_element,
1300 self.kind_id, self.codes,
1301 self.tmin_seconds, self.tmin_offset,
1302 self.tmax_seconds, self.tmax_offset, self.deltat)
1304 def diff(self, other):
1305 names = [
1306 'file_segment', 'file_element', 'kind_id', 'codes',
1307 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1308 'deltat']
1310 d = []
1311 for name, a, b in zip(
1312 names, self.equality_values, other.equality_values):
1314 if a != b:
1315 d.append((name, a, b))
1317 return d
1319 @property
1320 def tmin(self):
1321 return tjoin(self.tmin_seconds, self.tmin_offset)
1323 @tmin.setter
1324 def tmin(self, t):
1325 self.tmin_seconds, self.tmin_offset = tsplit(t)
1327 @property
1328 def tmax(self):
1329 return tjoin(self.tmax_seconds, self.tmax_offset)
1331 @tmax.setter
1332 def tmax(self, t):
1333 self.tmax_seconds, self.tmax_offset = tsplit(t)
1335 @property
1336 def kscale(self):
1337 if self.tmin_seconds is None or self.tmax_seconds is None:
1338 return 0
1339 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1341 @property
1342 def waveform_kwargs(self):
1343 network, station, location, channel, extra = self.codes
1345 return dict(
1346 network=network,
1347 station=station,
1348 location=location,
1349 channel=channel,
1350 extra=extra,
1351 tmin=self.tmin,
1352 tmax=self.tmax,
1353 deltat=self.deltat)
1355 @property
1356 def waveform_promise_kwargs(self):
1357 return dict(
1358 codes=self.codes,
1359 tmin=self.tmin,
1360 tmax=self.tmax,
1361 deltat=self.deltat)
1363 @property
1364 def station_kwargs(self):
1365 network, station, location = self.codes
1366 return dict(
1367 codes=self.codes,
1368 tmin=tmin_or_none(self.tmin),
1369 tmax=tmax_or_none(self.tmax))
1371 @property
1372 def channel_kwargs(self):
1373 network, station, location, channel, extra = self.codes
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 response_kwargs(self):
1382 return dict(
1383 codes=self.codes,
1384 tmin=tmin_or_none(self.tmin),
1385 tmax=tmax_or_none(self.tmax),
1386 deltat=self.deltat)
1388 @property
1389 def event_kwargs(self):
1390 return dict(
1391 name=self.codes,
1392 time=self.tmin,
1393 duration=(self.tmax - self.tmin) or None)
1395 @property
1396 def trace_kwargs(self):
1397 network, station, location, channel, extra = self.codes
1399 return dict(
1400 network=network,
1401 station=station,
1402 location=location,
1403 channel=channel,
1404 extra=extra,
1405 tmin=self.tmin,
1406 tmax=self.tmax-self.deltat,
1407 deltat=self.deltat)
1409 @property
1410 def dummy_trace(self):
1411 return DummyTrace(self)
1413 @property
1414 def summary(self):
1415 if self.tmin == self.tmax:
1416 ts = util.time_to_str(self.tmin)
1417 else:
1418 ts = '%s - %s' % (
1419 util.time_to_str(self.tmin),
1420 util.time_to_str(self.tmax))
1422 return ' '.join((
1423 ('%s,' % to_kind(self.kind_id)).ljust(9),
1424 ('%s,' % str(self.codes)).ljust(18),
1425 ts))
1428def make_waveform_nut(**kwargs):
1429 return Nut(kind_id=WAVEFORM, **kwargs)
1432def make_waveform_promise_nut(**kwargs):
1433 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1436def make_station_nut(**kwargs):
1437 return Nut(kind_id=STATION, **kwargs)
1440def make_channel_nut(**kwargs):
1441 return Nut(kind_id=CHANNEL, **kwargs)
1444def make_response_nut(**kwargs):
1445 return Nut(kind_id=RESPONSE, **kwargs)
1448def make_event_nut(**kwargs):
1449 return Nut(kind_id=EVENT, **kwargs)
1452def group_channels(nuts):
1453 by_ansl = {}
1454 for nut in nuts:
1455 if nut.kind_id != CHANNEL:
1456 continue
1458 ansl = nut.codes[:4]
1460 if ansl not in by_ansl:
1461 by_ansl[ansl] = {}
1463 group = by_ansl[ansl]
1465 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1467 if k not in group:
1468 group[k] = set()
1470 group.add(nut.codes[4])
1472 return by_ansl
1475class DummyTrace(object):
1477 def __init__(self, nut):
1478 self.nut = nut
1479 self.codes = nut.codes
1480 self.meta = {}
1482 @property
1483 def tmin(self):
1484 return self.nut.tmin
1486 @property
1487 def tmax(self):
1488 return self.nut.tmax
1490 @property
1491 def deltat(self):
1492 return self.nut.deltat
1494 @property
1495 def nslc_id(self):
1496 return self.codes.nslc
1498 @property
1499 def network(self):
1500 return self.codes.network
1502 @property
1503 def station(self):
1504 return self.codes.station
1506 @property
1507 def location(self):
1508 return self.codes.location
1510 @property
1511 def channel(self):
1512 return self.codes.channel
1514 @property
1515 def extra(self):
1516 return self.codes.extra
1518 def overlaps(self, tmin, tmax):
1519 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1522def duration_to_str(t):
1523 if t > 24*3600:
1524 return '%gd' % (t / (24.*3600.))
1525 elif t > 3600:
1526 return '%gh' % (t / 3600.)
1527 elif t > 60:
1528 return '%gm' % (t / 60.)
1529 else:
1530 return '%gs' % t
1533class Coverage(Object):
1534 '''
1535 Information about times covered by a waveform or other time series data.
1536 '''
1537 kind_id = Int.T(
1538 help='Content type.')
1539 pattern = Codes.T(
1540 help='The codes pattern in the request, which caused this entry to '
1541 'match.')
1542 codes = Codes.T(
1543 help='NSLCE or NSL codes identifier of the time series.')
1544 deltat = Float.T(
1545 help='Sampling interval [s]',
1546 optional=True)
1547 tmin = Timestamp.T(
1548 help='Global start time of time series.',
1549 optional=True)
1550 tmax = Timestamp.T(
1551 help='Global end time of time series.',
1552 optional=True)
1553 changes = List.T(
1554 Tuple.T(2, Any.T()),
1555 help='List of change points, with entries of the form '
1556 '``(time, count)``, where a ``count`` of zero indicates start of '
1557 'a gap, a value of 1 start of normal data coverage and a higher '
1558 'value duplicate or redundant data coverage.')
1560 @classmethod
1561 def from_values(cls, args):
1562 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1563 return cls(
1564 kind_id=kind_id,
1565 pattern=pattern,
1566 codes=codes,
1567 deltat=deltat,
1568 tmin=tmin,
1569 tmax=tmax,
1570 changes=changes)
1572 @property
1573 def summary(self):
1574 ts = '%s - %s,' % (
1575 util.time_to_str(self.tmin),
1576 util.time_to_str(self.tmax))
1578 srate = self.sample_rate
1580 return ' '.join((
1581 ('%s,' % to_kind(self.kind_id)).ljust(9),
1582 ('%s,' % str(self.codes)).ljust(18),
1583 ts,
1584 '%10.3g,' % srate if srate else '',
1585 '%4i' % len(self.changes),
1586 '%s' % duration_to_str(self.total)))
1588 @property
1589 def sample_rate(self):
1590 if self.deltat is None:
1591 return None
1592 elif self.deltat == 0.0:
1593 return 0.0
1594 else:
1595 return 1.0 / self.deltat
1597 @property
1598 def labels(self):
1599 srate = self.sample_rate
1600 return (
1601 ('%s' % str(self.codes)),
1602 '%.3g' % srate if srate else '')
1604 @property
1605 def total(self):
1606 total_t = None
1607 for tmin, tmax, _ in self.iter_spans():
1608 total_t = (total_t or 0.0) + (tmax - tmin)
1610 return total_t
1612 def iter_spans(self):
1613 last = None
1614 for (t, count) in self.changes:
1615 if last is not None:
1616 last_t, last_count = last
1617 if last_count > 0:
1618 yield last_t, t, last_count
1620 last = (t, count)
1622 def iter_uncovered_by(self, other):
1623 a = self
1624 b = other
1625 ia = ib = -1
1626 ca = cb = 0
1627 last = None
1628 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1629 if ib + 1 == len(b.changes):
1630 ia += 1
1631 t, ca = a.changes[ia]
1632 elif ia + 1 == len(a.changes):
1633 ib += 1
1634 t, cb = b.changes[ib]
1635 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1636 ia += 1
1637 t, ca = a.changes[ia]
1638 else:
1639 ib += 1
1640 t, cb = b.changes[ib]
1642 if last is not None:
1643 tl, cal, cbl = last
1644 if tl < t and cal > 0 and cbl == 0:
1645 yield tl, t, ia, ib
1647 last = (t, ca, cb)
1649 def iter_uncovered_by_combined(self, other):
1650 ib_last = None
1651 group = None
1652 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1653 if ib_last is None or ib != ib_last:
1654 if group:
1655 yield (group[0][0], group[-1][1])
1657 group = []
1659 group.append((tmin, tmax))
1660 ib_last = ib
1662 if group:
1663 yield (group[0][0], group[-1][1])
1666class LocationPool(object):
1668 def __init__(self, squirrel, tmin, tmax):
1670 locations = {}
1671 for station in squirrel.get_stations(tmin=tmin, tmax=tmax):
1672 c = station.codes
1673 if c not in locations:
1674 locations[c] = station
1675 else:
1676 if locations[c] is not None \
1677 and not locations[c].same_location(station):
1679 locations[c] = None
1681 for channel in squirrel.get_channels(tmin=tmin, tmax=tmax):
1682 c = channel.codes
1683 if c not in locations:
1684 locations[c] = channel
1685 else:
1686 if locations[c] is not None \
1687 and not locations[c].same_location(channel):
1689 locations[c] = None
1691 self._locations = locations
1693 def get(self, codes):
1694 try:
1695 return self._locations[codes]
1696 except KeyError:
1697 pass
1699 try:
1700 return self._locations[codes.codes_nsl]
1701 except KeyError:
1702 pass
1704 try:
1705 return self._locations[codes.codes_nsl_star]
1706 except KeyError:
1707 return None
1710__all__ = [
1711 'to_codes',
1712 'to_codes_guess',
1713 'to_codes_simple',
1714 'to_kind',
1715 'to_kinds',
1716 'to_kind_id',
1717 'to_kind_ids',
1718 'CodesError',
1719 'Codes',
1720 'CodesNSLCE',
1721 'CodesNSL',
1722 'CodesX',
1723 'Station',
1724 'Channel',
1725 'Sensor',
1726 'Response',
1727 'Nut',
1728 'Coverage',
1729 'WaveformPromise',
1730]