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)
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, stages=(None, 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[slice(*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()))
1198 time_created = Timestamp.T()
1199 anxious = Duration.T(default=600.)
1201 @property
1202 def client(self):
1203 return self.source_id.split(':')[1]
1205 def describe(self, site='?'):
1206 return '%s:%s %s [%s - %s]' % (
1207 self.client, site, str(self.codes),
1208 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1210 def validate(self, tr):
1211 if tr.ydata.size == 0:
1212 raise InvalidWaveform(
1213 'waveform with zero data samples')
1215 if tr.deltat != self.deltat:
1216 raise InvalidWaveform(
1217 'incorrect sampling interval - waveform: %g s, '
1218 'meta-data: %g s' % (
1219 tr.deltat, self.deltat))
1221 if not num.all(num.isfinite(tr.ydata)):
1222 raise InvalidWaveform('waveform has NaN values')
1224 def estimate_nsamples(self):
1225 return int(round((self.tmax - self.tmin) / self.deltat))+1
1227 def is_near_real_time(self):
1228 return self.tmax > self.time_created - self.anxious
1231def order_summary(orders):
1232 codes_list = sorted(set(order.codes for order in orders))
1233 if len(codes_list) > 3:
1234 return '%i order%s: %s - %s' % (
1235 len(orders),
1236 util.plural_s(orders),
1237 str(codes_list[0]),
1238 str(codes_list[1]))
1240 else:
1241 return '%i order%s: %s' % (
1242 len(orders),
1243 util.plural_s(orders),
1244 ', '.join(str(codes) for codes in codes_list))
1247class Nut(Object):
1248 '''
1249 Index entry referencing an elementary piece of content.
1251 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1252 meta-information about individual pieces of waveforms, stations, channels,
1253 etc. together with the information where it was found or generated.
1254 '''
1256 file_path = String.T(optional=True)
1257 file_format = String.T(optional=True)
1258 file_mtime = Timestamp.T(optional=True)
1259 file_size = Int.T(optional=True)
1261 file_segment = Int.T(optional=True)
1262 file_element = Int.T(optional=True)
1264 kind_id = Int.T()
1265 codes = Codes.T()
1267 tmin_seconds = Int.T(default=0)
1268 tmin_offset = Int.T(default=0, optional=True)
1269 tmax_seconds = Int.T(default=0)
1270 tmax_offset = Int.T(default=0, optional=True)
1272 deltat = Float.T(default=0.0)
1274 content = Any.T(optional=True)
1275 raw_content = Dict.T(String.T(), Any.T())
1277 content_in_db = False
1279 def __init__(
1280 self,
1281 file_path=None,
1282 file_format=None,
1283 file_mtime=None,
1284 file_size=None,
1285 file_segment=None,
1286 file_element=None,
1287 kind_id=0,
1288 codes=CodesX(''),
1289 tmin_seconds=None,
1290 tmin_offset=0,
1291 tmax_seconds=None,
1292 tmax_offset=0,
1293 deltat=None,
1294 content=None,
1295 raw_content=None,
1296 tmin=None,
1297 tmax=None,
1298 values_nocheck=None):
1300 if values_nocheck is not None:
1301 (self.file_path, self.file_format, self.file_mtime,
1302 self.file_size,
1303 self.file_segment, self.file_element,
1304 self.kind_id, codes_safe_str,
1305 self.tmin_seconds, self.tmin_offset,
1306 self.tmax_seconds, self.tmax_offset,
1307 self.deltat) = values_nocheck
1309 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1310 self.deltat = self.deltat or None
1311 self.raw_content = {}
1312 self.content = None
1313 else:
1314 if tmin is not None:
1315 tmin_seconds, tmin_offset = tsplit(tmin)
1317 if tmax is not None:
1318 tmax_seconds, tmax_offset = tsplit(tmax)
1320 self.kind_id = int(kind_id)
1321 self.codes = codes
1322 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1323 self.tmin_offset = int(tmin_offset)
1324 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1325 self.tmax_offset = int(tmax_offset)
1326 self.deltat = float_or_none(deltat)
1327 self.file_path = str_or_none(file_path)
1328 self.file_segment = int_or_none(file_segment)
1329 self.file_element = int_or_none(file_element)
1330 self.file_format = str_or_none(file_format)
1331 self.file_mtime = float_or_none(file_mtime)
1332 self.file_size = int_or_none(file_size)
1333 self.content = content
1334 if raw_content is None:
1335 self.raw_content = {}
1336 else:
1337 self.raw_content = raw_content
1339 Object.__init__(self, init_props=False)
1341 def __eq__(self, other):
1342 return (isinstance(other, Nut) and
1343 self.equality_values == other.equality_values)
1345 def hash(self):
1346 return ehash(','.join(str(x) for x in self.key))
1348 def __ne__(self, other):
1349 return not (self == other)
1351 def get_io_backend(self):
1352 from . import io
1353 return io.get_backend(self.file_format)
1355 def file_modified(self):
1356 return self.get_io_backend().get_stats(self.file_path) \
1357 != (self.file_mtime, self.file_size)
1359 @property
1360 def dkey(self):
1361 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1363 @property
1364 def key(self):
1365 return (
1366 self.file_path,
1367 self.file_segment,
1368 self.file_element,
1369 self.file_mtime)
1371 @property
1372 def equality_values(self):
1373 return (
1374 self.file_segment, self.file_element,
1375 self.kind_id, self.codes,
1376 self.tmin_seconds, self.tmin_offset,
1377 self.tmax_seconds, self.tmax_offset, self.deltat)
1379 def diff(self, other):
1380 names = [
1381 'file_segment', 'file_element', 'kind_id', 'codes',
1382 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1383 'deltat']
1385 d = []
1386 for name, a, b in zip(
1387 names, self.equality_values, other.equality_values):
1389 if a != b:
1390 d.append((name, a, b))
1392 return d
1394 @property
1395 def tmin(self):
1396 return tjoin(self.tmin_seconds, self.tmin_offset)
1398 @tmin.setter
1399 def tmin(self, t):
1400 self.tmin_seconds, self.tmin_offset = tsplit(t)
1402 @property
1403 def tmax(self):
1404 return tjoin(self.tmax_seconds, self.tmax_offset)
1406 @tmax.setter
1407 def tmax(self, t):
1408 self.tmax_seconds, self.tmax_offset = tsplit(t)
1410 @property
1411 def kscale(self):
1412 if self.tmin_seconds is None or self.tmax_seconds is None:
1413 return 0
1414 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1416 @property
1417 def waveform_kwargs(self):
1418 network, station, location, channel, extra = self.codes
1420 return dict(
1421 network=network,
1422 station=station,
1423 location=location,
1424 channel=channel,
1425 extra=extra,
1426 tmin=self.tmin,
1427 tmax=self.tmax,
1428 deltat=self.deltat)
1430 @property
1431 def waveform_promise_kwargs(self):
1432 return dict(
1433 codes=self.codes,
1434 tmin=self.tmin,
1435 tmax=self.tmax,
1436 deltat=self.deltat)
1438 @property
1439 def station_kwargs(self):
1440 network, station, location = self.codes
1441 return dict(
1442 codes=self.codes,
1443 tmin=tmin_or_none(self.tmin),
1444 tmax=tmax_or_none(self.tmax))
1446 @property
1447 def channel_kwargs(self):
1448 network, station, location, channel, extra = self.codes
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 response_kwargs(self):
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 event_kwargs(self):
1465 return dict(
1466 name=self.codes,
1467 time=self.tmin,
1468 duration=(self.tmax - self.tmin) or None)
1470 @property
1471 def trace_kwargs(self):
1472 network, station, location, channel, extra = self.codes
1474 return dict(
1475 network=network,
1476 station=station,
1477 location=location,
1478 channel=channel,
1479 extra=extra,
1480 tmin=self.tmin,
1481 tmax=self.tmax-self.deltat,
1482 deltat=self.deltat)
1484 @property
1485 def dummy_trace(self):
1486 return DummyTrace(self)
1488 @property
1489 def summary_entries(self):
1490 if self.tmin == self.tmax:
1491 ts = util.time_to_str(self.tmin)
1492 else:
1493 ts = '%s - %s' % (
1494 util.time_to_str(self.tmin),
1495 util.time_to_str(self.tmax))
1497 return (
1498 self.__class__.__name__,
1499 to_kind(self.kind_id),
1500 str(self.codes),
1501 ts)
1503 @property
1504 def summary(self):
1505 return util.fmt_summary(self.summary_entries, (10, 16, 20, 0))
1508def make_waveform_nut(**kwargs):
1509 return Nut(kind_id=WAVEFORM, **kwargs)
1512def make_waveform_promise_nut(**kwargs):
1513 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1516def make_station_nut(**kwargs):
1517 return Nut(kind_id=STATION, **kwargs)
1520def make_channel_nut(**kwargs):
1521 return Nut(kind_id=CHANNEL, **kwargs)
1524def make_response_nut(**kwargs):
1525 return Nut(kind_id=RESPONSE, **kwargs)
1528def make_event_nut(**kwargs):
1529 return Nut(kind_id=EVENT, **kwargs)
1532def group_channels(nuts):
1533 by_ansl = {}
1534 for nut in nuts:
1535 if nut.kind_id != CHANNEL:
1536 continue
1538 ansl = nut.codes[:4]
1540 if ansl not in by_ansl:
1541 by_ansl[ansl] = {}
1543 group = by_ansl[ansl]
1545 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1547 if k not in group:
1548 group[k] = set()
1550 group.add(nut.codes[4])
1552 return by_ansl
1555class DummyTrace(object):
1557 def __init__(self, nut):
1558 self.nut = nut
1559 self.codes = nut.codes
1560 self.meta = {}
1562 @property
1563 def tmin(self):
1564 return self.nut.tmin
1566 @property
1567 def tmax(self):
1568 return self.nut.tmax
1570 @property
1571 def deltat(self):
1572 return self.nut.deltat
1574 @property
1575 def nslc_id(self):
1576 return self.codes.nslc
1578 @property
1579 def network(self):
1580 return self.codes.network
1582 @property
1583 def station(self):
1584 return self.codes.station
1586 @property
1587 def location(self):
1588 return self.codes.location
1590 @property
1591 def channel(self):
1592 return self.codes.channel
1594 @property
1595 def extra(self):
1596 return self.codes.extra
1598 def overlaps(self, tmin, tmax):
1599 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1602def duration_to_str(t):
1603 if t > 24*3600:
1604 return '%gd' % (t / (24.*3600.))
1605 elif t > 3600:
1606 return '%gh' % (t / 3600.)
1607 elif t > 60:
1608 return '%gm' % (t / 60.)
1609 else:
1610 return '%gs' % t
1613class Coverage(Object):
1614 '''
1615 Information about times covered by a waveform or other time series data.
1616 '''
1617 kind_id = Int.T(
1618 help='Content type.')
1619 pattern = Codes.T(
1620 help='The codes pattern in the request, which caused this entry to '
1621 'match.')
1622 codes = Codes.T(
1623 help='NSLCE or NSL codes identifier of the time series.')
1624 deltat = Float.T(
1625 help='Sampling interval [s]',
1626 optional=True)
1627 tmin = Timestamp.T(
1628 help='Global start time of time series.',
1629 optional=True)
1630 tmax = Timestamp.T(
1631 help='Global end time of time series.',
1632 optional=True)
1633 changes = List.T(
1634 Tuple.T(2, Any.T()),
1635 help='List of change points, with entries of the form '
1636 '``(time, count)``, where a ``count`` of zero indicates start of '
1637 'a gap, a value of 1 start of normal data coverage and a higher '
1638 'value duplicate or redundant data coverage.')
1640 @classmethod
1641 def from_values(cls, args):
1642 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1643 return cls(
1644 kind_id=kind_id,
1645 pattern=pattern,
1646 codes=codes,
1647 deltat=deltat,
1648 tmin=tmin,
1649 tmax=tmax,
1650 changes=changes)
1652 @property
1653 def summary_entries(self):
1654 ts = '%s - %s' % (
1655 util.time_to_str(self.tmin),
1656 util.time_to_str(self.tmax))
1658 srate = self.sample_rate
1660 total = self.total
1662 return (
1663 self.__class__.__name__,
1664 to_kind(self.kind_id),
1665 str(self.codes),
1666 ts,
1667 '%10.3g' % srate if srate else '',
1668 '%i' % len(self.changes),
1669 '%s' % duration_to_str(total) if total else 'none')
1671 @property
1672 def summary(self):
1673 return util.fmt_summary(
1674 self.summary_entries,
1675 (10, 16, 20, 55, 10, 4, 0))
1677 @property
1678 def sample_rate(self):
1679 if self.deltat is None:
1680 return None
1681 elif self.deltat == 0.0:
1682 return 0.0
1683 else:
1684 return 1.0 / self.deltat
1686 @property
1687 def labels(self):
1688 srate = self.sample_rate
1689 return (
1690 ('%s' % str(self.codes)),
1691 '%.4g' % srate if srate else '')
1693 @property
1694 def total(self):
1695 total_t = None
1696 for tmin, tmax, _ in self.iter_spans():
1697 total_t = (total_t or 0.0) + (tmax - tmin)
1699 return total_t
1701 def iter_spans(self):
1702 last = None
1703 for (t, count) in self.changes:
1704 if last is not None:
1705 last_t, last_count = last
1706 if last_count > 0:
1707 yield last_t, t, last_count
1709 last = (t, count)
1711 def iter_uncovered_by(self, other):
1712 a = self
1713 b = other
1714 ia = ib = -1
1715 ca = cb = 0
1716 last = None
1717 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1718 if ib + 1 == len(b.changes):
1719 ia += 1
1720 t, ca = a.changes[ia]
1721 elif ia + 1 == len(a.changes):
1722 ib += 1
1723 t, cb = b.changes[ib]
1724 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1725 ia += 1
1726 t, ca = a.changes[ia]
1727 else:
1728 ib += 1
1729 t, cb = b.changes[ib]
1731 if last is not None:
1732 tl, cal, cbl = last
1733 if tl < t and cal > 0 and cbl == 0:
1734 yield tl, t, ia, ib
1736 last = (t, ca, cb)
1738 def iter_uncovered_by_combined(self, other):
1739 ib_last = None
1740 group = None
1741 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1742 if ib_last is None or ib != ib_last:
1743 if group:
1744 yield (group[0][0], group[-1][1])
1746 group = []
1748 group.append((tmin, tmax))
1749 ib_last = ib
1751 if group:
1752 yield (group[0][0], group[-1][1])
1755__all__ = [
1756 'to_codes',
1757 'to_codes_guess',
1758 'to_codes_simple',
1759 'to_kind',
1760 'to_kinds',
1761 'to_kind_id',
1762 'to_kind_ids',
1763 'CodesError',
1764 'Codes',
1765 'CodesNSLCE',
1766 'CodesNSL',
1767 'CodesX',
1768 'Station',
1769 'Channel',
1770 'Sensor',
1771 'Response',
1772 'Nut',
1773 'Coverage',
1774 'WaveformPromise',
1775]