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_args(self):
663 def getattr_rep(k):
664 if k == 'codes':
665 return self.codes.replace(
666 channel=self.codes.channel[:-1] + '?')
667 else:
668 return getattr(self, k)
670 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
672 def _get_channel_args(self, component):
673 def getattr_rep(k):
674 if k == 'codes':
675 return self.codes.replace(
676 channel=self.codes.channel[:-1] + component)
677 else:
678 return getattr(self, k)
680 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
682 def _get_pyrocko_station_args(self):
683 return (
684 self.codes.network,
685 self.codes.station,
686 self.codes.location,
687 self.lat,
688 self.lon,
689 self.elevation,
690 self.depth,
691 self.north_shift,
692 self.east_shift)
695class Channel(ChannelBase):
696 '''
697 A channel of a seismic station.
698 '''
700 dip = Float.T(optional=True)
701 azimuth = Float.T(optional=True)
703 def get_pyrocko_channel(self):
704 from pyrocko import model
705 return model.Channel(*self._get_pyrocko_channel_args())
707 def _get_pyrocko_channel_args(self):
708 return (
709 self.codes.channel,
710 self.azimuth,
711 self.dip)
713 @property
714 def orientation_enz(self):
715 if None in (self.azimuth, self.dip):
716 return None
718 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
719 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
720 d = math.sin(self.dip*d2r)
721 return mkvec(e, n, -d)
724def cut_intervals(channels):
725 channels = list(channels)
726 times = set()
727 for channel in channels:
728 if channel.tmin is not None:
729 times.add(channel.tmin)
730 if channel.tmax is not None:
731 times.add(channel.tmax)
733 times = sorted(times)
735 if len(times) < 2:
736 return channels
738 channels_out = []
739 for channel in channels:
740 tstart = channel.tmin
741 for t in times:
742 if (channel.tmin is None or channel.tmin < t) \
743 and (channel.tmax is None or t < channel.tmax):
745 channel_out = clone(channel)
746 channel_out.tmin = tstart
747 channel_out.tmax = t
748 channels_out.append(channel_out)
749 tstart = t
751 if channel.tmax is None:
752 channel_out = clone(channel)
753 channel_out.tmin = tstart
754 channels_out.append(channel_out)
756 return channels_out
759class Sensor(ChannelBase):
760 '''
761 Representation of a channel group.
762 '''
764 channels = List.T(Channel.T())
766 @classmethod
767 def from_channels(cls, channels):
768 channels = cut_intervals(channels)
770 groups = defaultdict(list)
771 for channel in channels:
772 groups[channel._get_sensor_args()].append(channel)
774 return [
775 cls(channels=channels,
776 **dict(zip(ChannelBase.T.propnames, args)))
777 for args, channels in groups.items()]
779 def channel_vectors(self):
780 return num.vstack(
781 [channel.orientation_enz for channel in self.channels])
783 def projected_channels(self, system, component_names):
784 return [
785 Channel(
786 azimuth=math.atan2(e, n) * r2d,
787 dip=-math.asin(z) * r2d,
788 **dict(zip(
789 ChannelBase.T.propnames,
790 self._get_channel_args(comp))))
791 for comp, (e, n, z) in zip(component_names, system)]
793 def matrix_to(self, system, epsilon=1e-7):
794 m = num.dot(system, self.channel_vectors().T)
795 m[num.abs(m) < epsilon] = 0.0
796 return m
798 def projection_to(self, system, component_names):
799 return (
800 self.matrix_to(system),
801 self.channels,
802 self.projected_channels(system, component_names))
804 def projection_to_enz(self):
805 return self.projection_to(num.identity(3), 'ENZ')
807 def projection_to_trz(self, source, azimuth=None):
808 if azimuth is not None:
809 assert source is None
810 else:
811 azimuth = source.azibazi_to(self)[1] + 180.
813 return self.projection_to(num.array([
814 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
815 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
816 [0., 0., 1.]], dtype=float), 'TRZ')
818 def project_to_enz(self, traces):
819 from pyrocko import trace
821 matrix, in_channels, out_channels = self.projection_to_enz()
823 return trace.project(traces, matrix, in_channels, out_channels)
825 def project_to_trz(self, source, traces, azimuth=None):
826 from pyrocko import trace
828 matrix, in_channels, out_channels = self.projection_to_trz(
829 source, azimuth=azimuth)
831 return trace.project(traces, matrix, in_channels, out_channels)
834observational_quantities = [
835 'acceleration', 'velocity', 'displacement', 'pressure',
836 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
837 'temperature']
840technical_quantities = [
841 'voltage', 'counts']
844class QuantityType(StringChoice):
845 '''
846 Choice of observational or technical quantity.
848 SI units are used for all quantities, where applicable.
849 '''
850 choices = observational_quantities + technical_quantities
853class ResponseStage(Object):
854 '''
855 Representation of a response stage.
857 Components of a seismic recording system are represented as a sequence of
858 response stages, e.g. sensor, pre-amplifier, digitizer, digital
859 downsampling.
860 '''
861 input_quantity = QuantityType.T(optional=True)
862 input_sample_rate = Float.T(optional=True)
863 output_quantity = QuantityType.T(optional=True)
864 output_sample_rate = Float.T(optional=True)
865 elements = List.T(FrequencyResponse.T())
866 log = List.T(Tuple.T(3, String.T()))
868 @property
869 def stage_type(self):
870 if self.input_quantity in observational_quantities \
871 and self.output_quantity in observational_quantities:
872 return 'conversion'
874 if self.input_quantity in observational_quantities \
875 and self.output_quantity == 'voltage':
876 return 'sensor'
878 elif self.input_quantity == 'voltage' \
879 and self.output_quantity == 'voltage':
880 return 'electronics'
882 elif self.input_quantity == 'voltage' \
883 and self.output_quantity == 'counts':
884 return 'digitizer'
886 elif self.input_quantity == 'counts' \
887 and self.output_quantity == 'counts' \
888 and self.input_sample_rate != self.output_sample_rate:
889 return 'decimation'
891 elif self.input_quantity in observational_quantities \
892 and self.output_quantity == 'counts':
893 return 'combination'
895 else:
896 return 'unknown'
898 @property
899 def summary(self):
900 irate = self.input_sample_rate
901 orate = self.output_sample_rate
902 factor = None
903 if irate and orate:
904 factor = irate / orate
905 return 'ResponseStage, ' + (
906 '%s%s => %s%s%s' % (
907 self.input_quantity or '?',
908 ' @ %g Hz' % irate if irate else '',
909 self.output_quantity or '?',
910 ' @ %g Hz' % orate if orate else '',
911 ' :%g' % factor if factor else '')
912 )
914 def get_effective(self):
915 return MultiplyResponse(responses=list(self.elements))
918D = 'displacement'
919V = 'velocity'
920A = 'acceleration'
922g_converters = {
923 (V, D): IntegrationResponse(1),
924 (A, D): IntegrationResponse(2),
925 (D, V): DifferentiationResponse(1),
926 (A, V): IntegrationResponse(1),
927 (D, A): DifferentiationResponse(2),
928 (V, A): DifferentiationResponse(1)}
931def response_converters(input_quantity, output_quantity):
932 if input_quantity is None or input_quantity == output_quantity:
933 return []
935 if output_quantity is None:
936 raise ConversionError('Unspecified target quantity.')
938 try:
939 return [g_converters[input_quantity, output_quantity]]
941 except KeyError:
942 raise ConversionError('No rule to convert from "%s" to "%s".' % (
943 input_quantity, output_quantity))
946@squirrel_content
947class Response(Object):
948 '''
949 The instrument response of a seismic station channel.
950 '''
952 codes = CodesNSLCE.T()
953 tmin = Timestamp.T(optional=True)
954 tmax = Timestamp.T(optional=True)
956 stages = List.T(ResponseStage.T())
957 checkpoints = List.T(FrequencyResponseCheckpoint.T())
959 deltat = Float.T(optional=True)
960 log = List.T(Tuple.T(3, String.T()))
962 def __init__(self, **kwargs):
963 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
964 Object.__init__(self, **kwargs)
966 @property
967 def time_span(self):
968 return (self.tmin, self.tmax)
970 @property
971 def nstages(self):
972 return len(self.stages)
974 @property
975 def input_quantity(self):
976 return self.stages[0].input_quantity if self.stages else None
978 @property
979 def output_quantity(self):
980 return self.stages[-1].input_quantity if self.stages else None
982 @property
983 def output_sample_rate(self):
984 return self.stages[-1].output_sample_rate if self.stages else None
986 @property
987 def stages_summary(self):
988 def grouped(xs):
989 xs = list(xs)
990 g = []
991 for i in range(len(xs)):
992 g.append(xs[i])
993 if i+1 < len(xs) and xs[i+1] != xs[i]:
994 yield g
995 g = []
997 if g:
998 yield g
1000 return '+'.join(
1001 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1002 for g in grouped(stage.stage_type for stage in self.stages))
1004 @property
1005 def summary(self):
1006 orate = self.output_sample_rate
1007 return '%s %-16s %s' % (
1008 self.__class__.__name__, self.str_codes, self.str_time_span) \
1009 + ', ' + ', '.join((
1010 '%s => %s' % (
1011 self.input_quantity or '?', self.output_quantity or '?')
1012 + (' @ %g Hz' % orate if orate else ''),
1013 self.stages_summary,
1014 ))
1016 def get_effective(self, input_quantity=None):
1017 try:
1018 elements = response_converters(input_quantity, self.input_quantity)
1019 except ConversionError as e:
1020 raise ConversionError(str(e) + ' (%s)' % self.summary)
1022 elements.extend(
1023 stage.get_effective() for stage in self.stages)
1025 return MultiplyResponse(responses=simplify_responses(elements))
1028@squirrel_content
1029class Event(Object):
1030 '''
1031 A seismic event.
1032 '''
1034 name = String.T(optional=True)
1035 time = Timestamp.T()
1036 duration = Float.T(optional=True)
1038 lat = Float.T()
1039 lon = Float.T()
1040 depth = Float.T(optional=True)
1042 magnitude = Float.T(optional=True)
1044 def get_pyrocko_event(self):
1045 from pyrocko import model
1046 model.Event(
1047 name=self.name,
1048 time=self.time,
1049 lat=self.lat,
1050 lon=self.lon,
1051 depth=self.depth,
1052 magnitude=self.magnitude,
1053 duration=self.duration)
1055 @property
1056 def time_span(self):
1057 return (self.time, self.time)
1060def ehash(s):
1061 return hashlib.sha1(s.encode('utf8')).hexdigest()
1064def random_name(n=8):
1065 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1068@squirrel_content
1069class WaveformPromise(Object):
1070 '''
1071 Information about a waveform potentially downloadable from a remote site.
1073 In the Squirrel framework, waveform promises are used to permit download of
1074 selected waveforms from a remote site. They are typically generated by
1075 calls to
1076 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1077 Waveform promises are inserted and indexed in the database similar to
1078 normal waveforms. When processing a waveform query, e.g. from
1079 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1080 waveform is available for the queried time span, a matching promise can be
1081 resolved, i.e. an attempt is made to download the waveform from the remote
1082 site. The promise is removed after the download attempt (except when a
1083 network error occurs). This prevents Squirrel from making unnecessary
1084 queries for waveforms missing at the remote site.
1085 '''
1087 codes = CodesNSLCE.T()
1088 tmin = Timestamp.T()
1089 tmax = Timestamp.T()
1091 deltat = Float.T(optional=True)
1093 source_hash = String.T()
1095 def __init__(self, **kwargs):
1096 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1097 Object.__init__(self, **kwargs)
1099 @property
1100 def time_span(self):
1101 return (self.tmin, self.tmax)
1104class InvalidWaveform(Exception):
1105 pass
1108class WaveformOrder(Object):
1109 '''
1110 Waveform request information.
1111 '''
1113 source_id = String.T()
1114 codes = CodesNSLCE.T()
1115 deltat = Float.T()
1116 tmin = Timestamp.T()
1117 tmax = Timestamp.T()
1118 gaps = List.T(Tuple.T(2, Timestamp.T()))
1120 @property
1121 def client(self):
1122 return self.source_id.split(':')[1]
1124 def describe(self, site='?'):
1125 return '%s:%s %s [%s - %s]' % (
1126 self.client, site, str(self.codes),
1127 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1129 def validate(self, tr):
1130 if tr.ydata.size == 0:
1131 raise InvalidWaveform(
1132 'waveform with zero data samples')
1134 if tr.deltat != self.deltat:
1135 raise InvalidWaveform(
1136 'incorrect sampling interval - waveform: %g s, '
1137 'meta-data: %g s' % (
1138 tr.deltat, self.deltat))
1140 if not num.all(num.isfinite(tr.ydata)):
1141 raise InvalidWaveform('waveform has NaN values')
1144def order_summary(orders):
1145 codes_list = sorted(set(order.codes for order in orders))
1146 if len(codes_list) > 3:
1147 return '%i order%s: %s - %s' % (
1148 len(orders),
1149 util.plural_s(orders),
1150 str(codes_list[0]),
1151 str(codes_list[1]))
1153 else:
1154 return '%i order%s: %s' % (
1155 len(orders),
1156 util.plural_s(orders),
1157 ', '.join(str(codes) for codes in codes_list))
1160class Nut(Object):
1161 '''
1162 Index entry referencing an elementary piece of content.
1164 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1165 meta-information about individual pieces of waveforms, stations, channels,
1166 etc. together with the information where it was found or generated.
1167 '''
1169 file_path = String.T(optional=True)
1170 file_format = String.T(optional=True)
1171 file_mtime = Timestamp.T(optional=True)
1172 file_size = Int.T(optional=True)
1174 file_segment = Int.T(optional=True)
1175 file_element = Int.T(optional=True)
1177 kind_id = Int.T()
1178 codes = Codes.T()
1180 tmin_seconds = Int.T(default=0)
1181 tmin_offset = Int.T(default=0, optional=True)
1182 tmax_seconds = Int.T(default=0)
1183 tmax_offset = Int.T(default=0, optional=True)
1185 deltat = Float.T(default=0.0)
1187 content = Any.T(optional=True)
1188 raw_content = Dict.T(String.T(), Any.T())
1190 content_in_db = False
1192 def __init__(
1193 self,
1194 file_path=None,
1195 file_format=None,
1196 file_mtime=None,
1197 file_size=None,
1198 file_segment=None,
1199 file_element=None,
1200 kind_id=0,
1201 codes=CodesX(''),
1202 tmin_seconds=None,
1203 tmin_offset=0,
1204 tmax_seconds=None,
1205 tmax_offset=0,
1206 deltat=None,
1207 content=None,
1208 raw_content=None,
1209 tmin=None,
1210 tmax=None,
1211 values_nocheck=None):
1213 if values_nocheck is not None:
1214 (self.file_path, self.file_format, self.file_mtime,
1215 self.file_size,
1216 self.file_segment, self.file_element,
1217 self.kind_id, codes_safe_str,
1218 self.tmin_seconds, self.tmin_offset,
1219 self.tmax_seconds, self.tmax_offset,
1220 self.deltat) = values_nocheck
1222 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1223 self.deltat = self.deltat or None
1224 self.raw_content = {}
1225 self.content = None
1226 else:
1227 if tmin is not None:
1228 tmin_seconds, tmin_offset = tsplit(tmin)
1230 if tmax is not None:
1231 tmax_seconds, tmax_offset = tsplit(tmax)
1233 self.kind_id = int(kind_id)
1234 self.codes = codes
1235 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1236 self.tmin_offset = int(tmin_offset)
1237 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1238 self.tmax_offset = int(tmax_offset)
1239 self.deltat = float_or_none(deltat)
1240 self.file_path = str_or_none(file_path)
1241 self.file_segment = int_or_none(file_segment)
1242 self.file_element = int_or_none(file_element)
1243 self.file_format = str_or_none(file_format)
1244 self.file_mtime = float_or_none(file_mtime)
1245 self.file_size = int_or_none(file_size)
1246 self.content = content
1247 if raw_content is None:
1248 self.raw_content = {}
1249 else:
1250 self.raw_content = raw_content
1252 Object.__init__(self, init_props=False)
1254 def __eq__(self, other):
1255 return (isinstance(other, Nut) and
1256 self.equality_values == other.equality_values)
1258 def hash(self):
1259 return ehash(','.join(str(x) for x in self.key))
1261 def __ne__(self, other):
1262 return not (self == other)
1264 def get_io_backend(self):
1265 from . import io
1266 return io.get_backend(self.file_format)
1268 def file_modified(self):
1269 return self.get_io_backend().get_stats(self.file_path) \
1270 != (self.file_mtime, self.file_size)
1272 @property
1273 def dkey(self):
1274 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1276 @property
1277 def key(self):
1278 return (
1279 self.file_path,
1280 self.file_segment,
1281 self.file_element,
1282 self.file_mtime)
1284 @property
1285 def equality_values(self):
1286 return (
1287 self.file_segment, self.file_element,
1288 self.kind_id, self.codes,
1289 self.tmin_seconds, self.tmin_offset,
1290 self.tmax_seconds, self.tmax_offset, self.deltat)
1292 def diff(self, other):
1293 names = [
1294 'file_segment', 'file_element', 'kind_id', 'codes',
1295 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1296 'deltat']
1298 d = []
1299 for name, a, b in zip(
1300 names, self.equality_values, other.equality_values):
1302 if a != b:
1303 d.append((name, a, b))
1305 return d
1307 @property
1308 def tmin(self):
1309 return tjoin(self.tmin_seconds, self.tmin_offset)
1311 @tmin.setter
1312 def tmin(self, t):
1313 self.tmin_seconds, self.tmin_offset = tsplit(t)
1315 @property
1316 def tmax(self):
1317 return tjoin(self.tmax_seconds, self.tmax_offset)
1319 @tmax.setter
1320 def tmax(self, t):
1321 self.tmax_seconds, self.tmax_offset = tsplit(t)
1323 @property
1324 def kscale(self):
1325 if self.tmin_seconds is None or self.tmax_seconds is None:
1326 return 0
1327 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1329 @property
1330 def waveform_kwargs(self):
1331 network, station, location, channel, extra = self.codes
1333 return dict(
1334 network=network,
1335 station=station,
1336 location=location,
1337 channel=channel,
1338 extra=extra,
1339 tmin=self.tmin,
1340 tmax=self.tmax,
1341 deltat=self.deltat)
1343 @property
1344 def waveform_promise_kwargs(self):
1345 return dict(
1346 codes=self.codes,
1347 tmin=self.tmin,
1348 tmax=self.tmax,
1349 deltat=self.deltat)
1351 @property
1352 def station_kwargs(self):
1353 network, station, location = self.codes
1354 return dict(
1355 codes=self.codes,
1356 tmin=tmin_or_none(self.tmin),
1357 tmax=tmax_or_none(self.tmax))
1359 @property
1360 def channel_kwargs(self):
1361 network, station, location, channel, extra = self.codes
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 response_kwargs(self):
1370 return dict(
1371 codes=self.codes,
1372 tmin=tmin_or_none(self.tmin),
1373 tmax=tmax_or_none(self.tmax),
1374 deltat=self.deltat)
1376 @property
1377 def event_kwargs(self):
1378 return dict(
1379 name=self.codes,
1380 time=self.tmin,
1381 duration=(self.tmax - self.tmin) or None)
1383 @property
1384 def trace_kwargs(self):
1385 network, station, location, channel, extra = self.codes
1387 return dict(
1388 network=network,
1389 station=station,
1390 location=location,
1391 channel=channel,
1392 extra=extra,
1393 tmin=self.tmin,
1394 tmax=self.tmax-self.deltat,
1395 deltat=self.deltat)
1397 @property
1398 def dummy_trace(self):
1399 return DummyTrace(self)
1401 @property
1402 def summary(self):
1403 if self.tmin == self.tmax:
1404 ts = util.time_to_str(self.tmin)
1405 else:
1406 ts = '%s - %s' % (
1407 util.time_to_str(self.tmin),
1408 util.time_to_str(self.tmax))
1410 return ' '.join((
1411 ('%s,' % to_kind(self.kind_id)).ljust(9),
1412 ('%s,' % str(self.codes)).ljust(18),
1413 ts))
1416def make_waveform_nut(**kwargs):
1417 return Nut(kind_id=WAVEFORM, **kwargs)
1420def make_waveform_promise_nut(**kwargs):
1421 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1424def make_station_nut(**kwargs):
1425 return Nut(kind_id=STATION, **kwargs)
1428def make_channel_nut(**kwargs):
1429 return Nut(kind_id=CHANNEL, **kwargs)
1432def make_response_nut(**kwargs):
1433 return Nut(kind_id=RESPONSE, **kwargs)
1436def make_event_nut(**kwargs):
1437 return Nut(kind_id=EVENT, **kwargs)
1440def group_channels(nuts):
1441 by_ansl = {}
1442 for nut in nuts:
1443 if nut.kind_id != CHANNEL:
1444 continue
1446 ansl = nut.codes[:4]
1448 if ansl not in by_ansl:
1449 by_ansl[ansl] = {}
1451 group = by_ansl[ansl]
1453 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1455 if k not in group:
1456 group[k] = set()
1458 group.add(nut.codes[4])
1460 return by_ansl
1463class DummyTrace(object):
1465 def __init__(self, nut):
1466 self.nut = nut
1467 self.codes = nut.codes
1468 self.meta = {}
1470 @property
1471 def tmin(self):
1472 return self.nut.tmin
1474 @property
1475 def tmax(self):
1476 return self.nut.tmax
1478 @property
1479 def deltat(self):
1480 return self.nut.deltat
1482 @property
1483 def nslc_id(self):
1484 return self.codes.nslc
1486 @property
1487 def network(self):
1488 return self.codes.network
1490 @property
1491 def station(self):
1492 return self.codes.station
1494 @property
1495 def location(self):
1496 return self.codes.location
1498 @property
1499 def channel(self):
1500 return self.codes.channel
1502 @property
1503 def extra(self):
1504 return self.codes.extra
1506 def overlaps(self, tmin, tmax):
1507 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1510def duration_to_str(t):
1511 if t > 24*3600:
1512 return '%gd' % (t / (24.*3600.))
1513 elif t > 3600:
1514 return '%gh' % (t / 3600.)
1515 elif t > 60:
1516 return '%gm' % (t / 60.)
1517 else:
1518 return '%gs' % t
1521class Coverage(Object):
1522 '''
1523 Information about times covered by a waveform or other time series data.
1524 '''
1525 kind_id = Int.T(
1526 help='Content type.')
1527 pattern = Codes.T(
1528 help='The codes pattern in the request, which caused this entry to '
1529 'match.')
1530 codes = Codes.T(
1531 help='NSLCE or NSL codes identifier of the time series.')
1532 deltat = Float.T(
1533 help='Sampling interval [s]',
1534 optional=True)
1535 tmin = Timestamp.T(
1536 help='Global start time of time series.',
1537 optional=True)
1538 tmax = Timestamp.T(
1539 help='Global end time of time series.',
1540 optional=True)
1541 changes = List.T(
1542 Tuple.T(2, Any.T()),
1543 help='List of change points, with entries of the form '
1544 '``(time, count)``, where a ``count`` of zero indicates start of '
1545 'a gap, a value of 1 start of normal data coverage and a higher '
1546 'value duplicate or redundant data coverage.')
1548 @classmethod
1549 def from_values(cls, args):
1550 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1551 return cls(
1552 kind_id=kind_id,
1553 pattern=pattern,
1554 codes=codes,
1555 deltat=deltat,
1556 tmin=tmin,
1557 tmax=tmax,
1558 changes=changes)
1560 @property
1561 def summary(self):
1562 ts = '%s - %s,' % (
1563 util.time_to_str(self.tmin),
1564 util.time_to_str(self.tmax))
1566 srate = self.sample_rate
1568 return ' '.join((
1569 ('%s,' % to_kind(self.kind_id)).ljust(9),
1570 ('%s,' % str(self.codes)).ljust(18),
1571 ts,
1572 '%10.3g,' % srate if srate else '',
1573 '%4i' % len(self.changes),
1574 '%s' % duration_to_str(self.total)))
1576 @property
1577 def sample_rate(self):
1578 if self.deltat is None:
1579 return None
1580 elif self.deltat == 0.0:
1581 return 0.0
1582 else:
1583 return 1.0 / self.deltat
1585 @property
1586 def labels(self):
1587 srate = self.sample_rate
1588 return (
1589 ('%s' % str(self.codes)),
1590 '%.3g' % srate if srate else '')
1592 @property
1593 def total(self):
1594 total_t = None
1595 for tmin, tmax, _ in self.iter_spans():
1596 total_t = (total_t or 0.0) + (tmax - tmin)
1598 return total_t
1600 def iter_spans(self):
1601 last = None
1602 for (t, count) in self.changes:
1603 if last is not None:
1604 last_t, last_count = last
1605 if last_count > 0:
1606 yield last_t, t, last_count
1608 last = (t, count)
1610 def iter_uncovered_by(self, other):
1611 a = self
1612 b = other
1613 ia = ib = -1
1614 ca = cb = 0
1615 last = None
1616 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1617 if ib + 1 == len(b.changes):
1618 ia += 1
1619 t, ca = a.changes[ia]
1620 elif ia + 1 == len(a.changes):
1621 ib += 1
1622 t, cb = b.changes[ib]
1623 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1624 ia += 1
1625 t, ca = a.changes[ia]
1626 else:
1627 ib += 1
1628 t, cb = b.changes[ib]
1630 if last is not None:
1631 tl, cal, cbl = last
1632 if tl < t and cal > 0 and cbl == 0:
1633 yield tl, t, ia, ib
1635 last = (t, ca, cb)
1637 def iter_uncovered_by_combined(self, other):
1638 ib_last = None
1639 group = None
1640 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1641 if ib_last is None or ib != ib_last:
1642 if group:
1643 yield (group[0][0], group[-1][1])
1645 group = []
1647 group.append((tmin, tmax))
1648 ib_last = ib
1650 if group:
1651 yield (group[0][0], group[-1][1])
1654class LocationPool(object):
1656 def __init__(self, squirrel, tmin, tmax):
1658 locations = {}
1659 for station in squirrel.get_stations(tmin=tmin, tmax=tmax):
1660 c = station.codes
1661 if c not in locations:
1662 locations[c] = station
1663 else:
1664 if locations[c] is not None \
1665 and not locations[c].same_location(station):
1667 locations[c] = None
1669 for channel in squirrel.get_channels(tmin=tmin, tmax=tmax):
1670 c = channel.codes
1671 if c not in locations:
1672 locations[c] = channel
1673 else:
1674 if locations[c] is not None \
1675 and not locations[c].same_location(channel):
1677 locations[c] = None
1679 self._locations = locations
1681 def get(self, codes):
1682 try:
1683 return self._locations[codes]
1684 except KeyError:
1685 pass
1687 try:
1688 return self._locations[codes.codes_nsl]
1689 except KeyError:
1690 pass
1692 try:
1693 return self._locations[codes.codes_nsl_star]
1694 except KeyError:
1695 return None
1698__all__ = [
1699 'to_codes',
1700 'to_codes_guess',
1701 'to_codes_simple',
1702 'to_kind',
1703 'to_kinds',
1704 'to_kind_id',
1705 'to_kind_ids',
1706 'CodesError',
1707 'Codes',
1708 'CodesNSLCE',
1709 'CodesNSL',
1710 'CodesX',
1711 'Station',
1712 'Channel',
1713 'Sensor',
1714 'Response',
1715 'Nut',
1716 'Coverage',
1717 'WaveformPromise',
1718]