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, Duration, 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)
1022 self._effective_responses_cache = {}
1024 @property
1025 def time_span(self):
1026 return (self.tmin, self.tmax)
1028 @property
1029 def nstages(self):
1030 return len(self.stages)
1032 @property
1033 def input_quantity(self):
1034 return self.stages[0].input_quantity if self.stages else None
1036 @property
1037 def output_quantity(self):
1038 return self.stages[-1].output_quantity if self.stages else None
1040 @property
1041 def output_sample_rate(self):
1042 return self.stages[-1].output_sample_rate if self.stages else None
1044 @property
1045 def summary_stages(self):
1046 def grouped(xs):
1047 xs = list(xs)
1048 g = []
1049 for i in range(len(xs)):
1050 g.append(xs[i])
1051 if i+1 < len(xs) and xs[i+1] != xs[i]:
1052 yield g
1053 g = []
1055 if g:
1056 yield g
1058 return '+'.join(
1059 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1060 for g in grouped(stage.stage_type for stage in self.stages))
1062 @property
1063 def summary_quantities(self):
1064 orate = self.output_sample_rate
1065 return '%s -> %s%s' % (
1066 self.input_quantity or '?',
1067 self.output_quantity or '?',
1068 ' @ %g Hz' % orate if orate else '')
1070 @property
1071 def summary_log(self):
1072 return ''.join(sorted(set(x[0][0].upper() for x in self.log)))
1074 @property
1075 def summary_entries(self):
1076 return (
1077 self.__class__.__name__,
1078 str(self.codes),
1079 self.str_time_span,
1080 self.summary_log,
1081 self.summary_quantities,
1082 self.summary_stages)
1084 @property
1085 def summary(self):
1086 return util.fmt_summary(self.summary_entries, (10, 20, 55, 3, 35, 0))
1088 def get_effective(self, input_quantity=None, stages=(None, None)):
1089 cache_key = (input_quantity, stages)
1090 if cache_key in self._effective_responses_cache:
1091 return self._effective_responses_cache[cache_key]
1093 try:
1094 elements = response_converters(input_quantity, self.input_quantity)
1095 except ConversionError as e:
1096 raise ConversionError(str(e) + ' (%s)' % self.summary)
1098 elements.extend(
1099 stage.get_effective() for stage in self.stages[slice(*stages)])
1101 if input_quantity is None \
1102 or input_quantity == self.input_quantity:
1103 checkpoints = self.checkpoints
1104 else:
1105 checkpoints = []
1107 resp = MultiplyResponse(
1108 responses=simplify_responses(elements),
1109 checkpoints=checkpoints)
1111 self._effective_responses_cache[cache_key] = resp
1112 return resp
1115@squirrel_content
1116class Event(Object):
1117 '''
1118 A seismic event.
1119 '''
1121 name = String.T(optional=True)
1122 time = Timestamp.T()
1123 duration = Float.T(optional=True)
1125 lat = Float.T()
1126 lon = Float.T()
1127 depth = Float.T(optional=True)
1129 magnitude = Float.T(optional=True)
1131 def get_pyrocko_event(self):
1132 from pyrocko import model
1133 model.Event(
1134 name=self.name,
1135 time=self.time,
1136 lat=self.lat,
1137 lon=self.lon,
1138 depth=self.depth,
1139 magnitude=self.magnitude,
1140 duration=self.duration)
1142 @property
1143 def time_span(self):
1144 return (self.time, self.time)
1147def ehash(s):
1148 return hashlib.sha1(s.encode('utf8')).hexdigest()
1151def random_name(n=8):
1152 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1155@squirrel_content
1156class WaveformPromise(Object):
1157 '''
1158 Information about a waveform potentially downloadable from a remote site.
1160 In the Squirrel framework, waveform promises are used to permit download of
1161 selected waveforms from a remote site. They are typically generated by
1162 calls to
1163 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1164 Waveform promises are inserted and indexed in the database similar to
1165 normal waveforms. When processing a waveform query, e.g. from
1166 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1167 waveform is available for the queried time span, a matching promise can be
1168 resolved, i.e. an attempt is made to download the waveform from the remote
1169 site. The promise is removed after the download attempt (except when a
1170 network error occurs). This prevents Squirrel from making unnecessary
1171 queries for waveforms missing at the remote site.
1172 '''
1174 codes = CodesNSLCE.T()
1175 tmin = Timestamp.T()
1176 tmax = Timestamp.T()
1178 deltat = Float.T(optional=True)
1180 source_hash = String.T()
1182 def __init__(self, **kwargs):
1183 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1184 Object.__init__(self, **kwargs)
1186 @property
1187 def time_span(self):
1188 return (self.tmin, self.tmax)
1191class InvalidWaveform(Exception):
1192 pass
1195class WaveformOrder(Object):
1196 '''
1197 Waveform request information.
1198 '''
1200 source_id = String.T()
1201 codes = CodesNSLCE.T()
1202 deltat = Float.T()
1203 tmin = Timestamp.T()
1204 tmax = Timestamp.T()
1205 gaps = List.T(Tuple.T(2, Timestamp.T()))
1206 time_created = Timestamp.T()
1207 anxious = Duration.T(default=600.)
1209 @property
1210 def client(self):
1211 return self.source_id.split(':')[1]
1213 def describe(self, site='?'):
1214 return '%s:%s %s [%s - %s]' % (
1215 self.client, site, str(self.codes),
1216 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1218 def validate(self, tr):
1219 if tr.ydata.size == 0:
1220 raise InvalidWaveform(
1221 'waveform with zero data samples')
1223 if tr.deltat != self.deltat:
1224 raise InvalidWaveform(
1225 'incorrect sampling interval - waveform: %g s, '
1226 'meta-data: %g s' % (
1227 tr.deltat, self.deltat))
1229 if not num.all(num.isfinite(tr.ydata)):
1230 raise InvalidWaveform('waveform has NaN values')
1232 def estimate_nsamples(self):
1233 return int(round((self.tmax - self.tmin) / self.deltat))+1
1235 def is_near_real_time(self):
1236 return self.tmax > self.time_created - self.anxious
1239def order_summary(orders):
1240 codes_list = sorted(set(order.codes for order in orders))
1241 if len(codes_list) > 3:
1242 return '%i order%s: %s - %s' % (
1243 len(orders),
1244 util.plural_s(orders),
1245 str(codes_list[0]),
1246 str(codes_list[1]))
1248 else:
1249 return '%i order%s: %s' % (
1250 len(orders),
1251 util.plural_s(orders),
1252 ', '.join(str(codes) for codes in codes_list))
1255class Nut(Object):
1256 '''
1257 Index entry referencing an elementary piece of content.
1259 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1260 meta-information about individual pieces of waveforms, stations, channels,
1261 etc. together with the information where it was found or generated.
1262 '''
1264 file_path = String.T(optional=True)
1265 file_format = String.T(optional=True)
1266 file_mtime = Timestamp.T(optional=True)
1267 file_size = Int.T(optional=True)
1269 file_segment = Int.T(optional=True)
1270 file_element = Int.T(optional=True)
1272 kind_id = Int.T()
1273 codes = Codes.T()
1275 tmin_seconds = Int.T(default=0)
1276 tmin_offset = Int.T(default=0, optional=True)
1277 tmax_seconds = Int.T(default=0)
1278 tmax_offset = Int.T(default=0, optional=True)
1280 deltat = Float.T(default=0.0)
1282 content = Any.T(optional=True)
1283 raw_content = Dict.T(String.T(), Any.T())
1285 content_in_db = False
1287 def __init__(
1288 self,
1289 file_path=None,
1290 file_format=None,
1291 file_mtime=None,
1292 file_size=None,
1293 file_segment=None,
1294 file_element=None,
1295 kind_id=0,
1296 codes=CodesX(''),
1297 tmin_seconds=None,
1298 tmin_offset=0,
1299 tmax_seconds=None,
1300 tmax_offset=0,
1301 deltat=None,
1302 content=None,
1303 raw_content=None,
1304 tmin=None,
1305 tmax=None,
1306 values_nocheck=None):
1308 if values_nocheck is not None:
1309 (self.file_path, self.file_format, self.file_mtime,
1310 self.file_size,
1311 self.file_segment, self.file_element,
1312 self.kind_id, codes_safe_str,
1313 self.tmin_seconds, self.tmin_offset,
1314 self.tmax_seconds, self.tmax_offset,
1315 self.deltat) = values_nocheck
1317 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1318 self.deltat = self.deltat or None
1319 self.raw_content = {}
1320 self.content = None
1321 else:
1322 if tmin is not None:
1323 tmin_seconds, tmin_offset = tsplit(tmin)
1325 if tmax is not None:
1326 tmax_seconds, tmax_offset = tsplit(tmax)
1328 self.kind_id = int(kind_id)
1329 self.codes = codes
1330 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1331 self.tmin_offset = int(tmin_offset)
1332 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1333 self.tmax_offset = int(tmax_offset)
1334 self.deltat = float_or_none(deltat)
1335 self.file_path = str_or_none(file_path)
1336 self.file_segment = int_or_none(file_segment)
1337 self.file_element = int_or_none(file_element)
1338 self.file_format = str_or_none(file_format)
1339 self.file_mtime = float_or_none(file_mtime)
1340 self.file_size = int_or_none(file_size)
1341 self.content = content
1342 if raw_content is None:
1343 self.raw_content = {}
1344 else:
1345 self.raw_content = raw_content
1347 Object.__init__(self, init_props=False)
1349 def __eq__(self, other):
1350 return (isinstance(other, Nut) and
1351 self.equality_values == other.equality_values)
1353 def hash(self):
1354 return ehash(','.join(str(x) for x in self.key))
1356 def __ne__(self, other):
1357 return not (self == other)
1359 def get_io_backend(self):
1360 from . import io
1361 return io.get_backend(self.file_format)
1363 def file_modified(self):
1364 return self.get_io_backend().get_stats(self.file_path) \
1365 != (self.file_mtime, self.file_size)
1367 @property
1368 def dkey(self):
1369 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1371 @property
1372 def key(self):
1373 return (
1374 self.file_path,
1375 self.file_segment,
1376 self.file_element,
1377 self.file_mtime)
1379 @property
1380 def equality_values(self):
1381 return (
1382 self.file_segment, self.file_element,
1383 self.kind_id, self.codes,
1384 self.tmin_seconds, self.tmin_offset,
1385 self.tmax_seconds, self.tmax_offset, self.deltat)
1387 def diff(self, other):
1388 names = [
1389 'file_segment', 'file_element', 'kind_id', 'codes',
1390 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1391 'deltat']
1393 d = []
1394 for name, a, b in zip(
1395 names, self.equality_values, other.equality_values):
1397 if a != b:
1398 d.append((name, a, b))
1400 return d
1402 @property
1403 def tmin(self):
1404 return tjoin(self.tmin_seconds, self.tmin_offset)
1406 @tmin.setter
1407 def tmin(self, t):
1408 self.tmin_seconds, self.tmin_offset = tsplit(t)
1410 @property
1411 def tmax(self):
1412 return tjoin(self.tmax_seconds, self.tmax_offset)
1414 @tmax.setter
1415 def tmax(self, t):
1416 self.tmax_seconds, self.tmax_offset = tsplit(t)
1418 @property
1419 def kscale(self):
1420 if self.tmin_seconds is None or self.tmax_seconds is None:
1421 return 0
1422 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1424 @property
1425 def waveform_kwargs(self):
1426 network, station, location, channel, extra = self.codes
1428 return dict(
1429 network=network,
1430 station=station,
1431 location=location,
1432 channel=channel,
1433 extra=extra,
1434 tmin=self.tmin,
1435 tmax=self.tmax,
1436 deltat=self.deltat)
1438 @property
1439 def waveform_promise_kwargs(self):
1440 return dict(
1441 codes=self.codes,
1442 tmin=self.tmin,
1443 tmax=self.tmax,
1444 deltat=self.deltat)
1446 @property
1447 def station_kwargs(self):
1448 network, station, location = self.codes
1449 return dict(
1450 codes=self.codes,
1451 tmin=tmin_or_none(self.tmin),
1452 tmax=tmax_or_none(self.tmax))
1454 @property
1455 def channel_kwargs(self):
1456 network, station, location, channel, extra = self.codes
1457 return dict(
1458 codes=self.codes,
1459 tmin=tmin_or_none(self.tmin),
1460 tmax=tmax_or_none(self.tmax),
1461 deltat=self.deltat)
1463 @property
1464 def response_kwargs(self):
1465 return dict(
1466 codes=self.codes,
1467 tmin=tmin_or_none(self.tmin),
1468 tmax=tmax_or_none(self.tmax),
1469 deltat=self.deltat)
1471 @property
1472 def event_kwargs(self):
1473 return dict(
1474 name=self.codes,
1475 time=self.tmin,
1476 duration=(self.tmax - self.tmin) or None)
1478 @property
1479 def trace_kwargs(self):
1480 network, station, location, channel, extra = self.codes
1482 return dict(
1483 network=network,
1484 station=station,
1485 location=location,
1486 channel=channel,
1487 extra=extra,
1488 tmin=self.tmin,
1489 tmax=self.tmax-self.deltat,
1490 deltat=self.deltat)
1492 @property
1493 def dummy_trace(self):
1494 return DummyTrace(self)
1496 @property
1497 def summary_entries(self):
1498 if self.tmin == self.tmax:
1499 ts = util.time_to_str(self.tmin)
1500 else:
1501 ts = '%s - %s' % (
1502 util.time_to_str(self.tmin),
1503 util.time_to_str(self.tmax))
1505 return (
1506 self.__class__.__name__,
1507 to_kind(self.kind_id),
1508 str(self.codes),
1509 ts)
1511 @property
1512 def summary(self):
1513 return util.fmt_summary(self.summary_entries, (10, 16, 20, 0))
1516def make_waveform_nut(**kwargs):
1517 return Nut(kind_id=WAVEFORM, **kwargs)
1520def make_waveform_promise_nut(**kwargs):
1521 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1524def make_station_nut(**kwargs):
1525 return Nut(kind_id=STATION, **kwargs)
1528def make_channel_nut(**kwargs):
1529 return Nut(kind_id=CHANNEL, **kwargs)
1532def make_response_nut(**kwargs):
1533 return Nut(kind_id=RESPONSE, **kwargs)
1536def make_event_nut(**kwargs):
1537 return Nut(kind_id=EVENT, **kwargs)
1540def group_channels(nuts):
1541 by_ansl = {}
1542 for nut in nuts:
1543 if nut.kind_id != CHANNEL:
1544 continue
1546 ansl = nut.codes[:4]
1548 if ansl not in by_ansl:
1549 by_ansl[ansl] = {}
1551 group = by_ansl[ansl]
1553 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1555 if k not in group:
1556 group[k] = set()
1558 group.add(nut.codes[4])
1560 return by_ansl
1563class DummyTrace(object):
1565 def __init__(self, nut):
1566 self.nut = nut
1567 self.codes = nut.codes
1568 self.meta = {}
1570 @property
1571 def tmin(self):
1572 return self.nut.tmin
1574 @property
1575 def tmax(self):
1576 return self.nut.tmax
1578 @property
1579 def deltat(self):
1580 return self.nut.deltat
1582 @property
1583 def nslc_id(self):
1584 return self.codes.nslc
1586 @property
1587 def network(self):
1588 return self.codes.network
1590 @property
1591 def station(self):
1592 return self.codes.station
1594 @property
1595 def location(self):
1596 return self.codes.location
1598 @property
1599 def channel(self):
1600 return self.codes.channel
1602 @property
1603 def extra(self):
1604 return self.codes.extra
1606 def overlaps(self, tmin, tmax):
1607 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1610def duration_to_str(t):
1611 if t > 24*3600:
1612 return '%gd' % (t / (24.*3600.))
1613 elif t > 3600:
1614 return '%gh' % (t / 3600.)
1615 elif t > 60:
1616 return '%gm' % (t / 60.)
1617 else:
1618 return '%gs' % t
1621class Coverage(Object):
1622 '''
1623 Information about times covered by a waveform or other time series data.
1624 '''
1625 kind_id = Int.T(
1626 help='Content type.')
1627 pattern = Codes.T(
1628 help='The codes pattern in the request, which caused this entry to '
1629 'match.')
1630 codes = Codes.T(
1631 help='NSLCE or NSL codes identifier of the time series.')
1632 deltat = Float.T(
1633 help='Sampling interval [s]',
1634 optional=True)
1635 tmin = Timestamp.T(
1636 help='Global start time of time series.',
1637 optional=True)
1638 tmax = Timestamp.T(
1639 help='Global end time of time series.',
1640 optional=True)
1641 changes = List.T(
1642 Tuple.T(2, Any.T()),
1643 help='List of change points, with entries of the form '
1644 '``(time, count)``, where a ``count`` of zero indicates start of '
1645 'a gap, a value of 1 start of normal data coverage and a higher '
1646 'value duplicate or redundant data coverage.')
1648 @classmethod
1649 def from_values(cls, args):
1650 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1651 return cls(
1652 kind_id=kind_id,
1653 pattern=pattern,
1654 codes=codes,
1655 deltat=deltat,
1656 tmin=tmin,
1657 tmax=tmax,
1658 changes=changes)
1660 @property
1661 def summary_entries(self):
1662 ts = '%s - %s' % (
1663 util.time_to_str(self.tmin),
1664 util.time_to_str(self.tmax))
1666 srate = self.sample_rate
1668 total = self.total
1670 return (
1671 self.__class__.__name__,
1672 to_kind(self.kind_id),
1673 str(self.codes),
1674 ts,
1675 '%10.3g' % srate if srate else '',
1676 '%i' % len(self.changes),
1677 '%s' % duration_to_str(total) if total else 'none')
1679 @property
1680 def summary(self):
1681 return util.fmt_summary(
1682 self.summary_entries,
1683 (10, 16, 20, 55, 10, 4, 0))
1685 @property
1686 def sample_rate(self):
1687 if self.deltat is None:
1688 return None
1689 elif self.deltat == 0.0:
1690 return 0.0
1691 else:
1692 return 1.0 / self.deltat
1694 @property
1695 def labels(self):
1696 srate = self.sample_rate
1697 return (
1698 ('%s' % str(self.codes)),
1699 '%.4g' % srate if srate else '')
1701 @property
1702 def total(self):
1703 total_t = None
1704 for tmin, tmax, _ in self.iter_spans():
1705 total_t = (total_t or 0.0) + (tmax - tmin)
1707 return total_t
1709 def iter_spans(self):
1710 last = None
1711 for (t, count) in self.changes:
1712 if last is not None:
1713 last_t, last_count = last
1714 if last_count > 0:
1715 yield last_t, t, last_count
1717 last = (t, count)
1719 def iter_uncovered_by(self, other):
1720 a = self
1721 b = other
1722 ia = ib = -1
1723 ca = cb = 0
1724 last = None
1725 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1726 if ib + 1 == len(b.changes):
1727 ia += 1
1728 t, ca = a.changes[ia]
1729 elif ia + 1 == len(a.changes):
1730 ib += 1
1731 t, cb = b.changes[ib]
1732 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1733 ia += 1
1734 t, ca = a.changes[ia]
1735 else:
1736 ib += 1
1737 t, cb = b.changes[ib]
1739 if last is not None:
1740 tl, cal, cbl = last
1741 if tl < t and cal > 0 and cbl == 0:
1742 yield tl, t, ia, ib
1744 last = (t, ca, cb)
1746 def iter_uncovered_by_combined(self, other):
1747 ib_last = None
1748 group = None
1749 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1750 if ib_last is None or ib != ib_last:
1751 if group:
1752 yield (group[0][0], group[-1][1])
1754 group = []
1756 group.append((tmin, tmax))
1757 ib_last = ib
1759 if group:
1760 yield (group[0][0], group[-1][1])
1763__all__ = [
1764 'to_codes',
1765 'to_codes_guess',
1766 'to_codes_simple',
1767 'to_kind',
1768 'to_kinds',
1769 'to_kind_id',
1770 'to_kind_ids',
1771 'CodesError',
1772 'Codes',
1773 'CodesNSLCE',
1774 'CodesNSL',
1775 'CodesX',
1776 'Station',
1777 'Channel',
1778 'Sensor',
1779 'Response',
1780 'Nut',
1781 'Coverage',
1782 'WaveformPromise',
1783]