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'''
19import re
20import fnmatch
21import math
22import hashlib
23from os import urandom
24from base64 import urlsafe_b64encode
25from collections import defaultdict, namedtuple
27import numpy as num
29from pyrocko import util
30from pyrocko.guts import Object, SObject, String, Timestamp, Float, Int, \
31 Unicode, Tuple, List, StringChoice, Any, Dict, clone
32from pyrocko.model import squirrel_content, Location
33from pyrocko.response import FrequencyResponse, MultiplyResponse, \
34 IntegrationResponse, DifferentiationResponse, simplify_responses, \
35 FrequencyResponseCheckpoint
37from .error import ConversionError, SquirrelError
39d2r = num.pi / 180.
40r2d = 1.0 / d2r
43guts_prefix = 'squirrel'
46def mkvec(x, y, z):
47 return num.array([x, y, z], dtype=float)
50def are_orthogonal(vecs, eps=0.01):
51 return all(abs(
52 num.dot(vecs[i], vecs[j]) < eps
53 for (i, j) in [(0, 1), (1, 2), (2, 0)]))
56g_codes_pool = {}
59class CodesError(SquirrelError):
60 pass
63class Codes(SObject):
64 pass
67def normalize_nslce(*args, **kwargs):
68 if args and kwargs:
69 raise ValueError('Either *args or **kwargs accepted, not both.')
71 if len(args) == 1:
72 if isinstance(args[0], str):
73 args = tuple(args[0].split('.'))
74 elif isinstance(args[0], tuple):
75 args = args[0]
76 else:
77 raise ValueError('Invalid argument type: %s' % type(args[0]))
79 nargs = len(args)
80 if nargs == 5:
81 t = args
83 elif nargs == 4:
84 t = args + ('',)
86 elif nargs == 0:
87 d = dict(
88 network='',
89 station='',
90 location='',
91 channel='',
92 extra='')
94 d.update(kwargs)
95 t = tuple(kwargs.get(k, '') for k in (
96 'network', 'station', 'location', 'channel', 'extra'))
98 else:
99 raise CodesError(
100 'Does not match NSLC or NSLCE codes pattern: %s' % '.'.join(args))
102 if '.'.join(t).count('.') != 4:
103 raise CodesError(
104 'Codes may not contain a ".": "%s", "%s", "%s", "%s", "%s"' % t)
106 return t
109CodesNSLCEBase = namedtuple(
110 'CodesNSLCEBase', [
111 'network', 'station', 'location', 'channel', 'extra'])
114class CodesNSLCE(CodesNSLCEBase, Codes):
115 '''
116 Codes denominating a seismic channel (NSLC or NSLCE).
118 FDSN/SEED style NET.STA.LOC.CHA is accepted or NET.STA.LOC.CHA.EXTRA, where
119 the EXTRA part in the latter form can be used to identify a custom
120 processing applied to a channel.
121 '''
123 __slots__ = ()
124 __hash__ = CodesNSLCEBase.__hash__
126 as_dict = CodesNSLCEBase._asdict
128 def __new__(cls, *args, safe_str=None, **kwargs):
129 nargs = len(args)
130 if nargs == 1 and isinstance(args[0], CodesNSLCE):
131 return args[0]
132 elif nargs == 1 and isinstance(args[0], CodesNSL):
133 t = (args[0].nsl) + ('*', '*')
134 elif nargs == 1 and isinstance(args[0], CodesX):
135 t = ('*', '*', '*', '*', '*')
136 elif safe_str is not None:
137 t = safe_str.split('.')
138 else:
139 t = normalize_nslce(*args, **kwargs)
141 x = CodesNSLCEBase.__new__(cls, *t)
142 return g_codes_pool.setdefault(x, x)
144 def __init__(self, *args, **kwargs):
145 Codes.__init__(self)
147 def __str__(self):
148 if self.extra == '':
149 return '.'.join(self[:-1])
150 else:
151 return '.'.join(self)
153 def __eq__(self, other):
154 if not isinstance(other, CodesNSLCE):
155 other = CodesNSLCE(other)
157 return CodesNSLCEBase.__eq__(self, other)
159 def matches(self, pattern):
160 if not isinstance(pattern, CodesNSLCE):
161 pattern = CodesNSLCE(pattern)
163 return match_codes(pattern, self)
165 @property
166 def safe_str(self):
167 return '.'.join(self)
169 @property
170 def nslce(self):
171 return self[:5]
173 @property
174 def nslc(self):
175 return self[:4]
177 @property
178 def nsl(self):
179 return self[:3]
181 @property
182 def ns(self):
183 return self[:2]
185 @property
186 def codes_nsl(self):
187 return CodesNSL(self)
189 @property
190 def codes_nsl_star(self):
191 return CodesNSL(self.network, self.station, '*')
193 def as_tuple(self):
194 return tuple(self)
196 def replace(self, **kwargs):
197 x = CodesNSLCEBase._replace(self, **kwargs)
198 return g_codes_pool.setdefault(x, x)
201def normalize_nsl(*args, **kwargs):
202 if args and kwargs:
203 raise ValueError('Either *args or **kwargs accepted, not both.')
205 if len(args) == 1:
206 if isinstance(args[0], str):
207 args = tuple(args[0].split('.'))
208 elif isinstance(args[0], tuple):
209 args = args[0]
210 else:
211 raise ValueError('Invalid argument type: %s' % type(args[0]))
213 nargs = len(args)
214 if nargs == 3:
215 t = args
217 elif nargs == 0:
218 d = dict(
219 network='',
220 station='',
221 location='')
223 d.update(kwargs)
224 t = tuple(kwargs.get(k, '') for k in (
225 'network', 'station', 'location'))
227 else:
228 raise CodesError(
229 'Does not match NSL codes pattern: %s' % '.'.join(args))
231 if '.'.join(t).count('.') != 2:
232 raise CodesError(
233 'Codes may not contain a ".": "%s", "%s", "%s"' % t)
235 return t
238CodesNSLBase = namedtuple(
239 'CodesNSLBase', [
240 'network', 'station', 'location'])
243class CodesNSL(CodesNSLBase, Codes):
244 '''
245 Codes denominating a seismic station (NSL).
247 NET.STA.LOC is accepted, slightly different from SEED/StationXML, where
248 LOC is part of the channel. By setting location='*' is possible to get
249 compatible behaviour in most cases.
250 '''
252 __slots__ = ()
253 __hash__ = CodesNSLBase.__hash__
255 as_dict = CodesNSLBase._asdict
257 def __new__(cls, *args, safe_str=None, **kwargs):
258 nargs = len(args)
259 if nargs == 1 and isinstance(args[0], CodesNSL):
260 return args[0]
261 elif nargs == 1 and isinstance(args[0], CodesNSLCE):
262 t = args[0].nsl
263 elif nargs == 1 and isinstance(args[0], CodesX):
264 t = ('*', '*', '*')
265 elif safe_str is not None:
266 t = safe_str.split('.')
267 else:
268 t = normalize_nsl(*args, **kwargs)
270 x = CodesNSLBase.__new__(cls, *t)
271 return g_codes_pool.setdefault(x, x)
273 def __init__(self, *args, **kwargs):
274 Codes.__init__(self)
276 def __str__(self):
277 return '.'.join(self)
279 def __eq__(self, other):
280 if not isinstance(other, CodesNSL):
281 other = CodesNSL(other)
283 return CodesNSLBase.__eq__(self, other)
285 def matches(self, pattern):
286 if not isinstance(pattern, CodesNSL):
287 pattern = CodesNSL(pattern)
289 return match_codes(pattern, self)
291 @property
292 def safe_str(self):
293 return '.'.join(self)
295 @property
296 def ns(self):
297 return self[:2]
299 @property
300 def nsl(self):
301 return self[:3]
303 def as_tuple(self):
304 return tuple(self)
306 def replace(self, **kwargs):
307 x = CodesNSLBase._replace(self, **kwargs)
308 return g_codes_pool.setdefault(x, x)
311CodesXBase = namedtuple(
312 'CodesXBase', [
313 'name'])
316class CodesX(CodesXBase, Codes):
317 '''
318 General purpose codes for anything other than channels or stations.
319 '''
321 __slots__ = ()
322 __hash__ = CodesXBase.__hash__
323 __eq__ = CodesXBase.__eq__
325 as_dict = CodesXBase._asdict
327 def __new__(cls, name='', safe_str=None):
328 if isinstance(name, CodesX):
329 return name
330 elif isinstance(name, (CodesNSLCE, CodesNSL)):
331 name = '*'
332 elif safe_str is not None:
333 name = safe_str
334 else:
335 if '.' in name:
336 raise CodesError('Code may not contain a ".": %s' % name)
338 x = CodesXBase.__new__(cls, name)
339 return g_codes_pool.setdefault(x, x)
341 def __init__(self, *args, **kwargs):
342 Codes.__init__(self)
344 def __str__(self):
345 return '.'.join(self)
347 @property
348 def safe_str(self):
349 return '.'.join(self)
351 @property
352 def ns(self):
353 return self[:2]
355 def as_tuple(self):
356 return tuple(self)
358 def replace(self, **kwargs):
359 x = CodesXBase._replace(self, **kwargs)
360 return g_codes_pool.setdefault(x, x)
363g_codes_patterns = {}
366def match_codes(pattern, codes):
367 spattern = pattern.safe_str
368 scodes = codes.safe_str
369 if spattern not in g_codes_patterns:
370 rpattern = re.compile(fnmatch.translate(spattern), re.I)
371 g_codes_patterns[spattern] = rpattern
373 rpattern = g_codes_patterns[spattern]
374 return bool(rpattern.match(scodes))
377g_content_kinds = [
378 'undefined',
379 'waveform',
380 'station',
381 'channel',
382 'response',
383 'event',
384 'waveform_promise']
387g_codes_classes = [
388 CodesX,
389 CodesNSLCE,
390 CodesNSL,
391 CodesNSLCE,
392 CodesNSLCE,
393 CodesX,
394 CodesNSLCE]
396g_codes_classes_ndot = {
397 0: CodesX,
398 2: CodesNSL,
399 3: CodesNSLCE,
400 4: CodesNSLCE}
403def to_codes_simple(kind_id, codes_safe_str):
404 return g_codes_classes[kind_id](safe_str=codes_safe_str)
407def to_codes(kind_id, obj):
408 return g_codes_classes[kind_id](obj)
411def to_codes_guess(s):
412 try:
413 return g_codes_classes_ndot[s.count('.')](s)
414 except KeyError:
415 raise CodesError('Cannot guess codes type: %s' % s)
418# derived list class to enable detection of already preprocessed codes patterns
419class codes_patterns_list(list):
420 pass
423def codes_patterns_for_kind(kind_id, codes):
424 if isinstance(codes, codes_patterns_list):
425 return codes
427 if isinstance(codes, list):
428 lcodes = codes_patterns_list()
429 for sc in codes:
430 lcodes.extend(codes_patterns_for_kind(kind_id, sc))
432 return lcodes
434 codes = to_codes(kind_id, codes)
436 lcodes = codes_patterns_list()
437 lcodes.append(codes)
438 if kind_id == STATION:
439 lcodes.append(codes.replace(location='[*]'))
441 return lcodes
444g_content_kind_ids = (
445 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT,
446 WAVEFORM_PROMISE) = range(len(g_content_kinds))
449g_tmin, g_tmax = util.get_working_system_time_range()[:2]
452try:
453 g_tmin_queries = max(g_tmin, util.str_to_time_fillup('1900-01-01'))
454except Exception:
455 g_tmin_queries = g_tmin
458def to_kind(kind_id):
459 return g_content_kinds[kind_id]
462def to_kinds(kind_ids):
463 return [g_content_kinds[kind_id] for kind_id in kind_ids]
466def to_kind_id(kind):
467 return g_content_kinds.index(kind)
470def to_kind_ids(kinds):
471 return [g_content_kinds.index(kind) for kind in kinds]
474g_kind_mask_all = 0xff
477def to_kind_mask(kinds):
478 if kinds:
479 return sum(1 << kind_id for kind_id in to_kind_ids(kinds))
480 else:
481 return g_kind_mask_all
484def str_or_none(x):
485 if x is None:
486 return None
487 else:
488 return str(x)
491def float_or_none(x):
492 if x is None:
493 return None
494 else:
495 return float(x)
498def int_or_none(x):
499 if x is None:
500 return None
501 else:
502 return int(x)
505def int_or_g_tmin(x):
506 if x is None:
507 return g_tmin
508 else:
509 return int(x)
512def int_or_g_tmax(x):
513 if x is None:
514 return g_tmax
515 else:
516 return int(x)
519def tmin_or_none(x):
520 if x == g_tmin:
521 return None
522 else:
523 return x
526def tmax_or_none(x):
527 if x == g_tmax:
528 return None
529 else:
530 return x
533def time_or_none_to_str(x):
534 if x is None:
535 return '...'.ljust(17)
536 else:
537 return util.time_to_str(x)
540def codes_to_str_abbreviated(codes, indent=' '):
541 codes = [str(x) for x in codes]
543 if len(codes) > 20:
544 scodes = '\n' + util.ewrap(codes[:10], indent=indent) \
545 + '\n%s[%i more]\n' % (indent, len(codes) - 20) \
546 + util.ewrap(codes[-10:], indent=' ')
547 else:
548 scodes = '\n' + util.ewrap(codes, indent=indent) \
549 if codes else '<none>'
551 return scodes
554g_offset_time_unit_inv = 1000000000
555g_offset_time_unit = 1.0 / g_offset_time_unit_inv
558def tsplit(t):
559 if t is None:
560 return None, 0
562 t = util.to_time_float(t)
563 if type(t) is float:
564 t = round(t, 5)
565 else:
566 t = round(t, 9)
568 seconds = num.floor(t)
569 offset = t - seconds
570 return int(seconds), int(round(offset * g_offset_time_unit_inv))
573def tjoin(seconds, offset):
574 if seconds is None:
575 return None
577 return util.to_time_float(seconds) \
578 + util.to_time_float(offset*g_offset_time_unit)
581tscale_min = 1
582tscale_max = 365 * 24 * 3600 # last edge is one above
583tscale_logbase = 20
585tscale_edges = [tscale_min]
586while True:
587 tscale_edges.append(tscale_edges[-1]*tscale_logbase)
588 if tscale_edges[-1] >= tscale_max:
589 break
592tscale_edges = num.array(tscale_edges)
595def tscale_to_kscale(tscale):
597 # 0 <= x < tscale_edges[1]: 0
598 # tscale_edges[1] <= x < tscale_edges[2]: 1
599 # ...
600 # tscale_edges[len(tscale_edges)-1] <= x: len(tscale_edges)
602 return int(num.searchsorted(tscale_edges, tscale))
605@squirrel_content
606class Station(Location):
607 '''
608 A seismic station.
609 '''
611 codes = CodesNSL.T()
613 tmin = Timestamp.T(optional=True)
614 tmax = Timestamp.T(optional=True)
616 description = Unicode.T(optional=True)
618 def __init__(self, **kwargs):
619 kwargs['codes'] = CodesNSL(kwargs['codes'])
620 Location.__init__(self, **kwargs)
622 @property
623 def time_span(self):
624 return (self.tmin, self.tmax)
626 def get_pyrocko_station(self):
627 from pyrocko import model
628 return model.Station(*self._get_pyrocko_station_args())
630 def _get_pyrocko_station_args(self):
631 return (
632 self.codes.network,
633 self.codes.station,
634 self.codes.location,
635 self.lat,
636 self.lon,
637 self.elevation,
638 self.depth,
639 self.north_shift,
640 self.east_shift)
643@squirrel_content
644class ChannelBase(Location):
645 codes = CodesNSLCE.T()
647 tmin = Timestamp.T(optional=True)
648 tmax = Timestamp.T(optional=True)
650 deltat = Float.T(optional=True)
652 def __init__(self, **kwargs):
653 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
654 Location.__init__(self, **kwargs)
656 @property
657 def time_span(self):
658 return (self.tmin, self.tmax)
660 def _get_sensor_codes(self):
661 return self.codes.replace(
662 channel=self.codes.channel[:-1] + '?')
664 def _get_sensor_args(self):
665 def getattr_rep(k):
666 if k == 'codes':
667 return self._get_sensor_codes()
668 else:
669 return getattr(self, k)
671 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
673 def _get_channel_args(self, component):
674 def getattr_rep(k):
675 if k == 'codes':
676 return self.codes.replace(
677 channel=self.codes.channel[:-1] + component)
678 else:
679 return getattr(self, k)
681 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
683 def _get_pyrocko_station_args(self):
684 return (
685 self.codes.network,
686 self.codes.station,
687 self.codes.location,
688 self.lat,
689 self.lon,
690 self.elevation,
691 self.depth,
692 self.north_shift,
693 self.east_shift)
696class Channel(ChannelBase):
697 '''
698 A channel of a seismic station.
699 '''
701 dip = Float.T(optional=True)
702 azimuth = Float.T(optional=True)
704 def get_pyrocko_channel(self):
705 from pyrocko import model
706 return model.Channel(*self._get_pyrocko_channel_args())
708 def _get_pyrocko_channel_args(self):
709 return (
710 self.codes.channel,
711 self.azimuth,
712 self.dip)
714 @property
715 def orientation_enz(self):
716 if None in (self.azimuth, self.dip):
717 return None
719 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
720 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
721 d = math.sin(self.dip*d2r)
722 return mkvec(e, n, -d)
725def cut_intervals(channels):
726 channels = list(channels)
727 times = set()
728 for channel in channels:
729 if channel.tmin is not None:
730 times.add(channel.tmin)
731 if channel.tmax is not None:
732 times.add(channel.tmax)
734 times = sorted(times)
736 if any(channel.tmin is None for channel in channels):
737 times[0:0] = [None]
739 if any(channel.tmax is None for channel in channels):
740 times.append(None)
742 if len(times) <= 2:
743 return channels
745 channels_out = []
746 for channel in channels:
747 for i in range(len(times)-1):
748 tmin = times[i]
749 tmax = times[i+1]
750 if ((channel.tmin is None or (
751 tmin is not None and channel.tmin <= tmin))
752 and (channel.tmax is None or (
753 tmax is not None and tmax <= channel.tmax))):
755 channel_out = clone(channel)
756 channel_out.tmin = tmin
757 channel_out.tmax = tmax
758 channels_out.append(channel_out)
760 return channels_out
763class Sensor(ChannelBase):
764 '''
765 Representation of a channel group.
766 '''
768 channels = List.T(Channel.T())
770 @classmethod
771 def from_channels(cls, channels):
772 groups = defaultdict(list)
773 for channel in channels:
774 groups[channel._get_sensor_codes()].append(channel)
776 channels_cut = []
777 for group in groups.values():
778 channels_cut.extend(cut_intervals(group))
780 groups = defaultdict(list)
781 for channel in channels_cut:
782 groups[channel._get_sensor_args()].append(channel)
784 return [
785 cls(channels=channels,
786 **dict(zip(ChannelBase.T.propnames, args)))
787 for args, channels in groups.items()]
789 def channel_vectors(self):
790 return num.vstack(
791 [channel.orientation_enz for channel in self.channels])
793 def projected_channels(self, system, component_names):
794 return [
795 Channel(
796 azimuth=math.atan2(e, n) * r2d,
797 dip=-math.asin(z) * r2d,
798 **dict(zip(
799 ChannelBase.T.propnames,
800 self._get_channel_args(comp))))
801 for comp, (e, n, z) in zip(component_names, system)]
803 def matrix_to(self, system, epsilon=1e-7):
804 m = num.dot(system, self.channel_vectors().T)
805 m[num.abs(m) < epsilon] = 0.0
806 return m
808 def projection_to(self, system, component_names):
809 return (
810 self.matrix_to(system),
811 self.channels,
812 self.projected_channels(system, component_names))
814 def projection_to_enz(self):
815 return self.projection_to(num.identity(3), 'ENZ')
817 def projection_to_trz(self, source, azimuth=None):
818 if azimuth is not None:
819 assert source is None
820 else:
821 azimuth = source.azibazi_to(self)[1] + 180.
823 return self.projection_to(num.array([
824 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
825 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
826 [0., 0., 1.]], dtype=float), 'TRZ')
828 def project_to_enz(self, traces):
829 from pyrocko import trace
831 matrix, in_channels, out_channels = self.projection_to_enz()
833 return trace.project(traces, matrix, in_channels, out_channels)
835 def project_to_trz(self, source, traces, azimuth=None):
836 from pyrocko import trace
838 matrix, in_channels, out_channels = self.projection_to_trz(
839 source, azimuth=azimuth)
841 return trace.project(traces, matrix, in_channels, out_channels)
844observational_quantities = [
845 'acceleration', 'velocity', 'displacement', 'pressure',
846 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
847 'temperature']
850technical_quantities = [
851 'voltage', 'counts']
854class QuantityType(StringChoice):
855 '''
856 Choice of observational or technical quantity.
858 SI units are used for all quantities, where applicable.
859 '''
860 choices = observational_quantities + technical_quantities
863class ResponseStage(Object):
864 '''
865 Representation of a response stage.
867 Components of a seismic recording system are represented as a sequence of
868 response stages, e.g. sensor, pre-amplifier, digitizer, digital
869 downsampling.
870 '''
871 input_quantity = QuantityType.T(optional=True)
872 input_sample_rate = Float.T(optional=True)
873 output_quantity = QuantityType.T(optional=True)
874 output_sample_rate = Float.T(optional=True)
875 elements = List.T(FrequencyResponse.T())
876 log = List.T(Tuple.T(3, String.T()))
878 @property
879 def stage_type(self):
880 if self.input_quantity in observational_quantities \
881 and self.output_quantity in observational_quantities:
882 return 'conversion'
884 if self.input_quantity in observational_quantities \
885 and self.output_quantity == 'voltage':
886 return 'sensor'
888 elif self.input_quantity == 'voltage' \
889 and self.output_quantity == 'voltage':
890 return 'electronics'
892 elif self.input_quantity == 'voltage' \
893 and self.output_quantity == 'counts':
894 return 'digitizer'
896 elif self.decimation_factor is not None \
897 and (self.input_quantity is None or self.input_quantity == 'counts') \
898 and (self.output_quantity is None or self.output_quantity == 'counts'): # noqa
899 return 'decimation'
901 elif self.input_quantity in observational_quantities \
902 and self.output_quantity == 'counts':
903 return 'combination'
905 else:
906 return 'unknown'
908 @property
909 def decimation_factor(self):
910 irate = self.input_sample_rate
911 orate = self.output_sample_rate
912 if irate is not None and orate is not None \
913 and irate > orate and irate / orate > 1.0:
915 return irate / orate
916 else:
917 return None
919 @property
920 def summary_quantities(self):
921 return '%s -> %s' % (
922 self.input_quantity or '?',
923 self.output_quantity or '?')
925 @property
926 def summary_rates(self):
927 irate = self.input_sample_rate
928 orate = self.output_sample_rate
929 factor = self.decimation_factor
931 if irate and orate is None:
932 return '%g Hz' % irate
934 elif orate and irate is None:
935 return '%g Hz' % orate
937 elif irate and orate and irate == orate:
938 return '%g Hz' % irate
940 elif any(x for x in (irate, orate, factor)):
941 return '%s -> %s Hz (%s)' % (
942 '%g' % irate if irate else '?',
943 '%g' % orate if orate else '?',
944 ':%g' % factor if factor else '?')
945 else:
946 return ''
948 @property
949 def summary_elements(self):
950 xs = [x.summary for x in self.elements]
951 return '%s' % ('*'.join(x for x in xs if x != 'one') or 'one')
953 @property
954 def summary_log(self):
955 return ''.join(sorted(set(x[0][0].upper() for x in self.log)))
957 @property
958 def summary_entries(self):
959 return (
960 self.__class__.__name__,
961 self.stage_type,
962 self.summary_log,
963 self.summary_quantities,
964 self.summary_rates,
965 self.summary_elements)
967 @property
968 def summary(self):
969 return util.fmt_summary(self.summary_entries, (10, 15, 3, 30, 30, 0))
971 def get_effective(self):
972 return MultiplyResponse(responses=list(self.elements))
975D = 'displacement'
976V = 'velocity'
977A = 'acceleration'
979g_converters = {
980 (V, D): IntegrationResponse(1),
981 (A, D): IntegrationResponse(2),
982 (D, V): DifferentiationResponse(1),
983 (A, V): IntegrationResponse(1),
984 (D, A): DifferentiationResponse(2),
985 (V, A): DifferentiationResponse(1)}
988def response_converters(input_quantity, output_quantity):
989 if input_quantity is None or input_quantity == output_quantity:
990 return []
992 if output_quantity is None:
993 raise ConversionError('Unspecified target quantity.')
995 try:
996 return [g_converters[input_quantity, output_quantity]]
998 except KeyError:
999 raise ConversionError('No rule to convert from "%s" to "%s".' % (
1000 input_quantity, output_quantity))
1003@squirrel_content
1004class Response(Object):
1005 '''
1006 The instrument response of a seismic station channel.
1007 '''
1009 codes = CodesNSLCE.T()
1010 tmin = Timestamp.T(optional=True)
1011 tmax = Timestamp.T(optional=True)
1013 stages = List.T(ResponseStage.T())
1014 checkpoints = List.T(FrequencyResponseCheckpoint.T())
1016 deltat = Float.T(optional=True)
1017 log = List.T(Tuple.T(3, String.T()))
1019 def __init__(self, **kwargs):
1020 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1021 Object.__init__(self, **kwargs)
1023 @property
1024 def time_span(self):
1025 return (self.tmin, self.tmax)
1027 @property
1028 def nstages(self):
1029 return len(self.stages)
1031 @property
1032 def input_quantity(self):
1033 return self.stages[0].input_quantity if self.stages else None
1035 @property
1036 def output_quantity(self):
1037 return self.stages[-1].output_quantity if self.stages else None
1039 @property
1040 def output_sample_rate(self):
1041 return self.stages[-1].output_sample_rate if self.stages else None
1043 @property
1044 def summary_stages(self):
1045 def grouped(xs):
1046 xs = list(xs)
1047 g = []
1048 for i in range(len(xs)):
1049 g.append(xs[i])
1050 if i+1 < len(xs) and xs[i+1] != xs[i]:
1051 yield g
1052 g = []
1054 if g:
1055 yield g
1057 return '+'.join(
1058 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1059 for g in grouped(stage.stage_type for stage in self.stages))
1061 @property
1062 def summary_quantities(self):
1063 orate = self.output_sample_rate
1064 return '%s -> %s%s' % (
1065 self.input_quantity or '?',
1066 self.output_quantity or '?',
1067 ' @ %g Hz' % orate if orate else '')
1069 @property
1070 def summary_log(self):
1071 return ''.join(sorted(set(x[0][0].upper() for x in self.log)))
1073 @property
1074 def summary_entries(self):
1075 return (
1076 self.__class__.__name__,
1077 str(self.codes),
1078 self.str_time_span,
1079 self.summary_log,
1080 self.summary_quantities,
1081 self.summary_stages)
1083 @property
1084 def summary(self):
1085 return util.fmt_summary(self.summary_entries, (10, 20, 55, 3, 35, 0))
1087 def get_effective(self, input_quantity=None):
1088 try:
1089 elements = response_converters(input_quantity, self.input_quantity)
1090 except ConversionError as e:
1091 raise ConversionError(str(e) + ' (%s)' % self.summary)
1093 elements.extend(
1094 stage.get_effective() for stage in self.stages)
1096 if input_quantity is None \
1097 or input_quantity == self.input_quantity:
1098 checkpoints = self.checkpoints
1099 else:
1100 checkpoints = []
1102 return MultiplyResponse(
1103 responses=simplify_responses(elements),
1104 checkpoints=checkpoints)
1107@squirrel_content
1108class Event(Object):
1109 '''
1110 A seismic event.
1111 '''
1113 name = String.T(optional=True)
1114 time = Timestamp.T()
1115 duration = Float.T(optional=True)
1117 lat = Float.T()
1118 lon = Float.T()
1119 depth = Float.T(optional=True)
1121 magnitude = Float.T(optional=True)
1123 def get_pyrocko_event(self):
1124 from pyrocko import model
1125 model.Event(
1126 name=self.name,
1127 time=self.time,
1128 lat=self.lat,
1129 lon=self.lon,
1130 depth=self.depth,
1131 magnitude=self.magnitude,
1132 duration=self.duration)
1134 @property
1135 def time_span(self):
1136 return (self.time, self.time)
1139def ehash(s):
1140 return hashlib.sha1(s.encode('utf8')).hexdigest()
1143def random_name(n=8):
1144 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1147@squirrel_content
1148class WaveformPromise(Object):
1149 '''
1150 Information about a waveform potentially downloadable from a remote site.
1152 In the Squirrel framework, waveform promises are used to permit download of
1153 selected waveforms from a remote site. They are typically generated by
1154 calls to
1155 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1156 Waveform promises are inserted and indexed in the database similar to
1157 normal waveforms. When processing a waveform query, e.g. from
1158 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1159 waveform is available for the queried time span, a matching promise can be
1160 resolved, i.e. an attempt is made to download the waveform from the remote
1161 site. The promise is removed after the download attempt (except when a
1162 network error occurs). This prevents Squirrel from making unnecessary
1163 queries for waveforms missing at the remote site.
1164 '''
1166 codes = CodesNSLCE.T()
1167 tmin = Timestamp.T()
1168 tmax = Timestamp.T()
1170 deltat = Float.T(optional=True)
1172 source_hash = String.T()
1174 def __init__(self, **kwargs):
1175 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1176 Object.__init__(self, **kwargs)
1178 @property
1179 def time_span(self):
1180 return (self.tmin, self.tmax)
1183class InvalidWaveform(Exception):
1184 pass
1187class WaveformOrder(Object):
1188 '''
1189 Waveform request information.
1190 '''
1192 source_id = String.T()
1193 codes = CodesNSLCE.T()
1194 deltat = Float.T()
1195 tmin = Timestamp.T()
1196 tmax = Timestamp.T()
1197 gaps = List.T(Tuple.T(2, Timestamp.T()))
1199 @property
1200 def client(self):
1201 return self.source_id.split(':')[1]
1203 def describe(self, site='?'):
1204 return '%s:%s %s [%s - %s]' % (
1205 self.client, site, str(self.codes),
1206 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1208 def validate(self, tr):
1209 if tr.ydata.size == 0:
1210 raise InvalidWaveform(
1211 'waveform with zero data samples')
1213 if tr.deltat != self.deltat:
1214 raise InvalidWaveform(
1215 'incorrect sampling interval - waveform: %g s, '
1216 'meta-data: %g s' % (
1217 tr.deltat, self.deltat))
1219 if not num.all(num.isfinite(tr.ydata)):
1220 raise InvalidWaveform('waveform has NaN values')
1223def order_summary(orders):
1224 codes_list = sorted(set(order.codes for order in orders))
1225 if len(codes_list) > 3:
1226 return '%i order%s: %s - %s' % (
1227 len(orders),
1228 util.plural_s(orders),
1229 str(codes_list[0]),
1230 str(codes_list[1]))
1232 else:
1233 return '%i order%s: %s' % (
1234 len(orders),
1235 util.plural_s(orders),
1236 ', '.join(str(codes) for codes in codes_list))
1239class Nut(Object):
1240 '''
1241 Index entry referencing an elementary piece of content.
1243 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1244 meta-information about individual pieces of waveforms, stations, channels,
1245 etc. together with the information where it was found or generated.
1246 '''
1248 file_path = String.T(optional=True)
1249 file_format = String.T(optional=True)
1250 file_mtime = Timestamp.T(optional=True)
1251 file_size = Int.T(optional=True)
1253 file_segment = Int.T(optional=True)
1254 file_element = Int.T(optional=True)
1256 kind_id = Int.T()
1257 codes = Codes.T()
1259 tmin_seconds = Int.T(default=0)
1260 tmin_offset = Int.T(default=0, optional=True)
1261 tmax_seconds = Int.T(default=0)
1262 tmax_offset = Int.T(default=0, optional=True)
1264 deltat = Float.T(default=0.0)
1266 content = Any.T(optional=True)
1267 raw_content = Dict.T(String.T(), Any.T())
1269 content_in_db = False
1271 def __init__(
1272 self,
1273 file_path=None,
1274 file_format=None,
1275 file_mtime=None,
1276 file_size=None,
1277 file_segment=None,
1278 file_element=None,
1279 kind_id=0,
1280 codes=CodesX(''),
1281 tmin_seconds=None,
1282 tmin_offset=0,
1283 tmax_seconds=None,
1284 tmax_offset=0,
1285 deltat=None,
1286 content=None,
1287 raw_content=None,
1288 tmin=None,
1289 tmax=None,
1290 values_nocheck=None):
1292 if values_nocheck is not None:
1293 (self.file_path, self.file_format, self.file_mtime,
1294 self.file_size,
1295 self.file_segment, self.file_element,
1296 self.kind_id, codes_safe_str,
1297 self.tmin_seconds, self.tmin_offset,
1298 self.tmax_seconds, self.tmax_offset,
1299 self.deltat) = values_nocheck
1301 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1302 self.deltat = self.deltat or None
1303 self.raw_content = {}
1304 self.content = None
1305 else:
1306 if tmin is not None:
1307 tmin_seconds, tmin_offset = tsplit(tmin)
1309 if tmax is not None:
1310 tmax_seconds, tmax_offset = tsplit(tmax)
1312 self.kind_id = int(kind_id)
1313 self.codes = codes
1314 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1315 self.tmin_offset = int(tmin_offset)
1316 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1317 self.tmax_offset = int(tmax_offset)
1318 self.deltat = float_or_none(deltat)
1319 self.file_path = str_or_none(file_path)
1320 self.file_segment = int_or_none(file_segment)
1321 self.file_element = int_or_none(file_element)
1322 self.file_format = str_or_none(file_format)
1323 self.file_mtime = float_or_none(file_mtime)
1324 self.file_size = int_or_none(file_size)
1325 self.content = content
1326 if raw_content is None:
1327 self.raw_content = {}
1328 else:
1329 self.raw_content = raw_content
1331 Object.__init__(self, init_props=False)
1333 def __eq__(self, other):
1334 return (isinstance(other, Nut) and
1335 self.equality_values == other.equality_values)
1337 def hash(self):
1338 return ehash(','.join(str(x) for x in self.key))
1340 def __ne__(self, other):
1341 return not (self == other)
1343 def get_io_backend(self):
1344 from . import io
1345 return io.get_backend(self.file_format)
1347 def file_modified(self):
1348 return self.get_io_backend().get_stats(self.file_path) \
1349 != (self.file_mtime, self.file_size)
1351 @property
1352 def dkey(self):
1353 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1355 @property
1356 def key(self):
1357 return (
1358 self.file_path,
1359 self.file_segment,
1360 self.file_element,
1361 self.file_mtime)
1363 @property
1364 def equality_values(self):
1365 return (
1366 self.file_segment, self.file_element,
1367 self.kind_id, self.codes,
1368 self.tmin_seconds, self.tmin_offset,
1369 self.tmax_seconds, self.tmax_offset, self.deltat)
1371 def diff(self, other):
1372 names = [
1373 'file_segment', 'file_element', 'kind_id', 'codes',
1374 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1375 'deltat']
1377 d = []
1378 for name, a, b in zip(
1379 names, self.equality_values, other.equality_values):
1381 if a != b:
1382 d.append((name, a, b))
1384 return d
1386 @property
1387 def tmin(self):
1388 return tjoin(self.tmin_seconds, self.tmin_offset)
1390 @tmin.setter
1391 def tmin(self, t):
1392 self.tmin_seconds, self.tmin_offset = tsplit(t)
1394 @property
1395 def tmax(self):
1396 return tjoin(self.tmax_seconds, self.tmax_offset)
1398 @tmax.setter
1399 def tmax(self, t):
1400 self.tmax_seconds, self.tmax_offset = tsplit(t)
1402 @property
1403 def kscale(self):
1404 if self.tmin_seconds is None or self.tmax_seconds is None:
1405 return 0
1406 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1408 @property
1409 def waveform_kwargs(self):
1410 network, station, location, channel, extra = self.codes
1412 return dict(
1413 network=network,
1414 station=station,
1415 location=location,
1416 channel=channel,
1417 extra=extra,
1418 tmin=self.tmin,
1419 tmax=self.tmax,
1420 deltat=self.deltat)
1422 @property
1423 def waveform_promise_kwargs(self):
1424 return dict(
1425 codes=self.codes,
1426 tmin=self.tmin,
1427 tmax=self.tmax,
1428 deltat=self.deltat)
1430 @property
1431 def station_kwargs(self):
1432 network, station, location = self.codes
1433 return dict(
1434 codes=self.codes,
1435 tmin=tmin_or_none(self.tmin),
1436 tmax=tmax_or_none(self.tmax))
1438 @property
1439 def channel_kwargs(self):
1440 network, station, location, channel, extra = self.codes
1441 return dict(
1442 codes=self.codes,
1443 tmin=tmin_or_none(self.tmin),
1444 tmax=tmax_or_none(self.tmax),
1445 deltat=self.deltat)
1447 @property
1448 def response_kwargs(self):
1449 return dict(
1450 codes=self.codes,
1451 tmin=tmin_or_none(self.tmin),
1452 tmax=tmax_or_none(self.tmax),
1453 deltat=self.deltat)
1455 @property
1456 def event_kwargs(self):
1457 return dict(
1458 name=self.codes,
1459 time=self.tmin,
1460 duration=(self.tmax - self.tmin) or None)
1462 @property
1463 def trace_kwargs(self):
1464 network, station, location, channel, extra = self.codes
1466 return dict(
1467 network=network,
1468 station=station,
1469 location=location,
1470 channel=channel,
1471 extra=extra,
1472 tmin=self.tmin,
1473 tmax=self.tmax-self.deltat,
1474 deltat=self.deltat)
1476 @property
1477 def dummy_trace(self):
1478 return DummyTrace(self)
1480 @property
1481 def summary_entries(self):
1482 if self.tmin == self.tmax:
1483 ts = util.time_to_str(self.tmin)
1484 else:
1485 ts = '%s - %s' % (
1486 util.time_to_str(self.tmin),
1487 util.time_to_str(self.tmax))
1489 return (
1490 self.__class__.__name__,
1491 to_kind(self.kind_id),
1492 str(self.codes),
1493 ts)
1495 @property
1496 def summary(self):
1497 return util.fmt_summary(self.summary_entries, (10, 16, 20, 0))
1500def make_waveform_nut(**kwargs):
1501 return Nut(kind_id=WAVEFORM, **kwargs)
1504def make_waveform_promise_nut(**kwargs):
1505 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1508def make_station_nut(**kwargs):
1509 return Nut(kind_id=STATION, **kwargs)
1512def make_channel_nut(**kwargs):
1513 return Nut(kind_id=CHANNEL, **kwargs)
1516def make_response_nut(**kwargs):
1517 return Nut(kind_id=RESPONSE, **kwargs)
1520def make_event_nut(**kwargs):
1521 return Nut(kind_id=EVENT, **kwargs)
1524def group_channels(nuts):
1525 by_ansl = {}
1526 for nut in nuts:
1527 if nut.kind_id != CHANNEL:
1528 continue
1530 ansl = nut.codes[:4]
1532 if ansl not in by_ansl:
1533 by_ansl[ansl] = {}
1535 group = by_ansl[ansl]
1537 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1539 if k not in group:
1540 group[k] = set()
1542 group.add(nut.codes[4])
1544 return by_ansl
1547class DummyTrace(object):
1549 def __init__(self, nut):
1550 self.nut = nut
1551 self.codes = nut.codes
1552 self.meta = {}
1554 @property
1555 def tmin(self):
1556 return self.nut.tmin
1558 @property
1559 def tmax(self):
1560 return self.nut.tmax
1562 @property
1563 def deltat(self):
1564 return self.nut.deltat
1566 @property
1567 def nslc_id(self):
1568 return self.codes.nslc
1570 @property
1571 def network(self):
1572 return self.codes.network
1574 @property
1575 def station(self):
1576 return self.codes.station
1578 @property
1579 def location(self):
1580 return self.codes.location
1582 @property
1583 def channel(self):
1584 return self.codes.channel
1586 @property
1587 def extra(self):
1588 return self.codes.extra
1590 def overlaps(self, tmin, tmax):
1591 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1594def duration_to_str(t):
1595 if t > 24*3600:
1596 return '%gd' % (t / (24.*3600.))
1597 elif t > 3600:
1598 return '%gh' % (t / 3600.)
1599 elif t > 60:
1600 return '%gm' % (t / 60.)
1601 else:
1602 return '%gs' % t
1605class Coverage(Object):
1606 '''
1607 Information about times covered by a waveform or other time series data.
1608 '''
1609 kind_id = Int.T(
1610 help='Content type.')
1611 pattern = Codes.T(
1612 help='The codes pattern in the request, which caused this entry to '
1613 'match.')
1614 codes = Codes.T(
1615 help='NSLCE or NSL codes identifier of the time series.')
1616 deltat = Float.T(
1617 help='Sampling interval [s]',
1618 optional=True)
1619 tmin = Timestamp.T(
1620 help='Global start time of time series.',
1621 optional=True)
1622 tmax = Timestamp.T(
1623 help='Global end time of time series.',
1624 optional=True)
1625 changes = List.T(
1626 Tuple.T(2, Any.T()),
1627 help='List of change points, with entries of the form '
1628 '``(time, count)``, where a ``count`` of zero indicates start of '
1629 'a gap, a value of 1 start of normal data coverage and a higher '
1630 'value duplicate or redundant data coverage.')
1632 @classmethod
1633 def from_values(cls, args):
1634 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1635 return cls(
1636 kind_id=kind_id,
1637 pattern=pattern,
1638 codes=codes,
1639 deltat=deltat,
1640 tmin=tmin,
1641 tmax=tmax,
1642 changes=changes)
1644 @property
1645 def summary_entries(self):
1646 ts = '%s - %s' % (
1647 util.time_to_str(self.tmin),
1648 util.time_to_str(self.tmax))
1650 srate = self.sample_rate
1652 total = self.total
1654 return (
1655 self.__class__.__name__,
1656 to_kind(self.kind_id),
1657 str(self.codes),
1658 ts,
1659 '%10.3g' % srate if srate else '',
1660 '%i' % len(self.changes),
1661 '%s' % duration_to_str(total) if total else 'none')
1663 @property
1664 def summary(self):
1665 return util.fmt_summary(
1666 self.summary_entries,
1667 (10, 16, 20, 55, 10, 4, 0))
1669 @property
1670 def sample_rate(self):
1671 if self.deltat is None:
1672 return None
1673 elif self.deltat == 0.0:
1674 return 0.0
1675 else:
1676 return 1.0 / self.deltat
1678 @property
1679 def labels(self):
1680 srate = self.sample_rate
1681 return (
1682 ('%s' % str(self.codes)),
1683 '%.3g' % srate if srate else '')
1685 @property
1686 def total(self):
1687 total_t = None
1688 for tmin, tmax, _ in self.iter_spans():
1689 total_t = (total_t or 0.0) + (tmax - tmin)
1691 return total_t
1693 def iter_spans(self):
1694 last = None
1695 for (t, count) in self.changes:
1696 if last is not None:
1697 last_t, last_count = last
1698 if last_count > 0:
1699 yield last_t, t, last_count
1701 last = (t, count)
1703 def iter_uncovered_by(self, other):
1704 a = self
1705 b = other
1706 ia = ib = -1
1707 ca = cb = 0
1708 last = None
1709 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1710 if ib + 1 == len(b.changes):
1711 ia += 1
1712 t, ca = a.changes[ia]
1713 elif ia + 1 == len(a.changes):
1714 ib += 1
1715 t, cb = b.changes[ib]
1716 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1717 ia += 1
1718 t, ca = a.changes[ia]
1719 else:
1720 ib += 1
1721 t, cb = b.changes[ib]
1723 if last is not None:
1724 tl, cal, cbl = last
1725 if tl < t and cal > 0 and cbl == 0:
1726 yield tl, t, ia, ib
1728 last = (t, ca, cb)
1730 def iter_uncovered_by_combined(self, other):
1731 ib_last = None
1732 group = None
1733 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1734 if ib_last is None or ib != ib_last:
1735 if group:
1736 yield (group[0][0], group[-1][1])
1738 group = []
1740 group.append((tmin, tmax))
1741 ib_last = ib
1743 if group:
1744 yield (group[0][0], group[-1][1])
1747__all__ = [
1748 'to_codes',
1749 'to_codes_guess',
1750 'to_codes_simple',
1751 'to_kind',
1752 'to_kinds',
1753 'to_kind_id',
1754 'to_kind_ids',
1755 'CodesError',
1756 'Codes',
1757 'CodesNSLCE',
1758 'CodesNSL',
1759 'CodesX',
1760 'Station',
1761 'Channel',
1762 'Sensor',
1763 'Response',
1764 'Nut',
1765 'Coverage',
1766 'WaveformPromise',
1767]