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.input_quantity == 'counts' \
897 and self.output_quantity == 'counts' \
898 and self.input_sample_rate != self.output_sample_rate:
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 summary(self):
910 irate = self.input_sample_rate
911 orate = self.output_sample_rate
912 factor = None
913 if irate and orate:
914 factor = irate / orate
915 return 'ResponseStage, ' + (
916 '%s%s => %s%s%s' % (
917 self.input_quantity or '?',
918 ' @ %g Hz' % irate if irate else '',
919 self.output_quantity or '?',
920 ' @ %g Hz' % orate if orate else '',
921 ' :%g' % factor if factor else '')
922 )
924 def get_effective(self):
925 return MultiplyResponse(responses=list(self.elements))
928D = 'displacement'
929V = 'velocity'
930A = 'acceleration'
932g_converters = {
933 (V, D): IntegrationResponse(1),
934 (A, D): IntegrationResponse(2),
935 (D, V): DifferentiationResponse(1),
936 (A, V): IntegrationResponse(1),
937 (D, A): DifferentiationResponse(2),
938 (V, A): DifferentiationResponse(1)}
941def response_converters(input_quantity, output_quantity):
942 if input_quantity is None or input_quantity == output_quantity:
943 return []
945 if output_quantity is None:
946 raise ConversionError('Unspecified target quantity.')
948 try:
949 return [g_converters[input_quantity, output_quantity]]
951 except KeyError:
952 raise ConversionError('No rule to convert from "%s" to "%s".' % (
953 input_quantity, output_quantity))
956@squirrel_content
957class Response(Object):
958 '''
959 The instrument response of a seismic station channel.
960 '''
962 codes = CodesNSLCE.T()
963 tmin = Timestamp.T(optional=True)
964 tmax = Timestamp.T(optional=True)
966 stages = List.T(ResponseStage.T())
967 checkpoints = List.T(FrequencyResponseCheckpoint.T())
969 deltat = Float.T(optional=True)
970 log = List.T(Tuple.T(3, String.T()))
972 def __init__(self, **kwargs):
973 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
974 Object.__init__(self, **kwargs)
976 @property
977 def time_span(self):
978 return (self.tmin, self.tmax)
980 @property
981 def nstages(self):
982 return len(self.stages)
984 @property
985 def input_quantity(self):
986 return self.stages[0].input_quantity if self.stages else None
988 @property
989 def output_quantity(self):
990 return self.stages[-1].input_quantity if self.stages else None
992 @property
993 def output_sample_rate(self):
994 return self.stages[-1].output_sample_rate if self.stages else None
996 @property
997 def stages_summary(self):
998 def grouped(xs):
999 xs = list(xs)
1000 g = []
1001 for i in range(len(xs)):
1002 g.append(xs[i])
1003 if i+1 < len(xs) and xs[i+1] != xs[i]:
1004 yield g
1005 g = []
1007 if g:
1008 yield g
1010 return '+'.join(
1011 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1012 for g in grouped(stage.stage_type for stage in self.stages))
1014 @property
1015 def summary(self):
1016 orate = self.output_sample_rate
1017 return '%s %-16s %s' % (
1018 self.__class__.__name__, self.str_codes, self.str_time_span) \
1019 + ', ' + ', '.join((
1020 '%s => %s' % (
1021 self.input_quantity or '?', self.output_quantity or '?')
1022 + (' @ %g Hz' % orate if orate else ''),
1023 self.stages_summary,
1024 ))
1026 def get_effective(self, input_quantity=None):
1027 try:
1028 elements = response_converters(input_quantity, self.input_quantity)
1029 except ConversionError as e:
1030 raise ConversionError(str(e) + ' (%s)' % self.summary)
1032 elements.extend(
1033 stage.get_effective() for stage in self.stages)
1035 return MultiplyResponse(responses=simplify_responses(elements))
1038@squirrel_content
1039class Event(Object):
1040 '''
1041 A seismic event.
1042 '''
1044 name = String.T(optional=True)
1045 time = Timestamp.T()
1046 duration = Float.T(optional=True)
1048 lat = Float.T()
1049 lon = Float.T()
1050 depth = Float.T(optional=True)
1052 magnitude = Float.T(optional=True)
1054 def get_pyrocko_event(self):
1055 from pyrocko import model
1056 model.Event(
1057 name=self.name,
1058 time=self.time,
1059 lat=self.lat,
1060 lon=self.lon,
1061 depth=self.depth,
1062 magnitude=self.magnitude,
1063 duration=self.duration)
1065 @property
1066 def time_span(self):
1067 return (self.time, self.time)
1070def ehash(s):
1071 return hashlib.sha1(s.encode('utf8')).hexdigest()
1074def random_name(n=8):
1075 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1078@squirrel_content
1079class WaveformPromise(Object):
1080 '''
1081 Information about a waveform potentially downloadable from a remote site.
1083 In the Squirrel framework, waveform promises are used to permit download of
1084 selected waveforms from a remote site. They are typically generated by
1085 calls to
1086 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1087 Waveform promises are inserted and indexed in the database similar to
1088 normal waveforms. When processing a waveform query, e.g. from
1089 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1090 waveform is available for the queried time span, a matching promise can be
1091 resolved, i.e. an attempt is made to download the waveform from the remote
1092 site. The promise is removed after the download attempt (except when a
1093 network error occurs). This prevents Squirrel from making unnecessary
1094 queries for waveforms missing at the remote site.
1095 '''
1097 codes = CodesNSLCE.T()
1098 tmin = Timestamp.T()
1099 tmax = Timestamp.T()
1101 deltat = Float.T(optional=True)
1103 source_hash = String.T()
1105 def __init__(self, **kwargs):
1106 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1107 Object.__init__(self, **kwargs)
1109 @property
1110 def time_span(self):
1111 return (self.tmin, self.tmax)
1114class InvalidWaveform(Exception):
1115 pass
1118class WaveformOrder(Object):
1119 '''
1120 Waveform request information.
1121 '''
1123 source_id = String.T()
1124 codes = CodesNSLCE.T()
1125 deltat = Float.T()
1126 tmin = Timestamp.T()
1127 tmax = Timestamp.T()
1128 gaps = List.T(Tuple.T(2, Timestamp.T()))
1130 @property
1131 def client(self):
1132 return self.source_id.split(':')[1]
1134 def describe(self, site='?'):
1135 return '%s:%s %s [%s - %s]' % (
1136 self.client, site, str(self.codes),
1137 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1139 def validate(self, tr):
1140 if tr.ydata.size == 0:
1141 raise InvalidWaveform(
1142 'waveform with zero data samples')
1144 if tr.deltat != self.deltat:
1145 raise InvalidWaveform(
1146 'incorrect sampling interval - waveform: %g s, '
1147 'meta-data: %g s' % (
1148 tr.deltat, self.deltat))
1150 if not num.all(num.isfinite(tr.ydata)):
1151 raise InvalidWaveform('waveform has NaN values')
1154def order_summary(orders):
1155 codes_list = sorted(set(order.codes for order in orders))
1156 if len(codes_list) > 3:
1157 return '%i order%s: %s - %s' % (
1158 len(orders),
1159 util.plural_s(orders),
1160 str(codes_list[0]),
1161 str(codes_list[1]))
1163 else:
1164 return '%i order%s: %s' % (
1165 len(orders),
1166 util.plural_s(orders),
1167 ', '.join(str(codes) for codes in codes_list))
1170class Nut(Object):
1171 '''
1172 Index entry referencing an elementary piece of content.
1174 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1175 meta-information about individual pieces of waveforms, stations, channels,
1176 etc. together with the information where it was found or generated.
1177 '''
1179 file_path = String.T(optional=True)
1180 file_format = String.T(optional=True)
1181 file_mtime = Timestamp.T(optional=True)
1182 file_size = Int.T(optional=True)
1184 file_segment = Int.T(optional=True)
1185 file_element = Int.T(optional=True)
1187 kind_id = Int.T()
1188 codes = Codes.T()
1190 tmin_seconds = Int.T(default=0)
1191 tmin_offset = Int.T(default=0, optional=True)
1192 tmax_seconds = Int.T(default=0)
1193 tmax_offset = Int.T(default=0, optional=True)
1195 deltat = Float.T(default=0.0)
1197 content = Any.T(optional=True)
1198 raw_content = Dict.T(String.T(), Any.T())
1200 content_in_db = False
1202 def __init__(
1203 self,
1204 file_path=None,
1205 file_format=None,
1206 file_mtime=None,
1207 file_size=None,
1208 file_segment=None,
1209 file_element=None,
1210 kind_id=0,
1211 codes=CodesX(''),
1212 tmin_seconds=None,
1213 tmin_offset=0,
1214 tmax_seconds=None,
1215 tmax_offset=0,
1216 deltat=None,
1217 content=None,
1218 raw_content=None,
1219 tmin=None,
1220 tmax=None,
1221 values_nocheck=None):
1223 if values_nocheck is not None:
1224 (self.file_path, self.file_format, self.file_mtime,
1225 self.file_size,
1226 self.file_segment, self.file_element,
1227 self.kind_id, codes_safe_str,
1228 self.tmin_seconds, self.tmin_offset,
1229 self.tmax_seconds, self.tmax_offset,
1230 self.deltat) = values_nocheck
1232 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1233 self.deltat = self.deltat or None
1234 self.raw_content = {}
1235 self.content = None
1236 else:
1237 if tmin is not None:
1238 tmin_seconds, tmin_offset = tsplit(tmin)
1240 if tmax is not None:
1241 tmax_seconds, tmax_offset = tsplit(tmax)
1243 self.kind_id = int(kind_id)
1244 self.codes = codes
1245 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1246 self.tmin_offset = int(tmin_offset)
1247 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1248 self.tmax_offset = int(tmax_offset)
1249 self.deltat = float_or_none(deltat)
1250 self.file_path = str_or_none(file_path)
1251 self.file_segment = int_or_none(file_segment)
1252 self.file_element = int_or_none(file_element)
1253 self.file_format = str_or_none(file_format)
1254 self.file_mtime = float_or_none(file_mtime)
1255 self.file_size = int_or_none(file_size)
1256 self.content = content
1257 if raw_content is None:
1258 self.raw_content = {}
1259 else:
1260 self.raw_content = raw_content
1262 Object.__init__(self, init_props=False)
1264 def __eq__(self, other):
1265 return (isinstance(other, Nut) and
1266 self.equality_values == other.equality_values)
1268 def hash(self):
1269 return ehash(','.join(str(x) for x in self.key))
1271 def __ne__(self, other):
1272 return not (self == other)
1274 def get_io_backend(self):
1275 from . import io
1276 return io.get_backend(self.file_format)
1278 def file_modified(self):
1279 return self.get_io_backend().get_stats(self.file_path) \
1280 != (self.file_mtime, self.file_size)
1282 @property
1283 def dkey(self):
1284 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1286 @property
1287 def key(self):
1288 return (
1289 self.file_path,
1290 self.file_segment,
1291 self.file_element,
1292 self.file_mtime)
1294 @property
1295 def equality_values(self):
1296 return (
1297 self.file_segment, self.file_element,
1298 self.kind_id, self.codes,
1299 self.tmin_seconds, self.tmin_offset,
1300 self.tmax_seconds, self.tmax_offset, self.deltat)
1302 def diff(self, other):
1303 names = [
1304 'file_segment', 'file_element', 'kind_id', 'codes',
1305 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1306 'deltat']
1308 d = []
1309 for name, a, b in zip(
1310 names, self.equality_values, other.equality_values):
1312 if a != b:
1313 d.append((name, a, b))
1315 return d
1317 @property
1318 def tmin(self):
1319 return tjoin(self.tmin_seconds, self.tmin_offset)
1321 @tmin.setter
1322 def tmin(self, t):
1323 self.tmin_seconds, self.tmin_offset = tsplit(t)
1325 @property
1326 def tmax(self):
1327 return tjoin(self.tmax_seconds, self.tmax_offset)
1329 @tmax.setter
1330 def tmax(self, t):
1331 self.tmax_seconds, self.tmax_offset = tsplit(t)
1333 @property
1334 def kscale(self):
1335 if self.tmin_seconds is None or self.tmax_seconds is None:
1336 return 0
1337 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1339 @property
1340 def waveform_kwargs(self):
1341 network, station, location, channel, extra = self.codes
1343 return dict(
1344 network=network,
1345 station=station,
1346 location=location,
1347 channel=channel,
1348 extra=extra,
1349 tmin=self.tmin,
1350 tmax=self.tmax,
1351 deltat=self.deltat)
1353 @property
1354 def waveform_promise_kwargs(self):
1355 return dict(
1356 codes=self.codes,
1357 tmin=self.tmin,
1358 tmax=self.tmax,
1359 deltat=self.deltat)
1361 @property
1362 def station_kwargs(self):
1363 network, station, location = self.codes
1364 return dict(
1365 codes=self.codes,
1366 tmin=tmin_or_none(self.tmin),
1367 tmax=tmax_or_none(self.tmax))
1369 @property
1370 def channel_kwargs(self):
1371 network, station, location, channel, extra = self.codes
1372 return dict(
1373 codes=self.codes,
1374 tmin=tmin_or_none(self.tmin),
1375 tmax=tmax_or_none(self.tmax),
1376 deltat=self.deltat)
1378 @property
1379 def response_kwargs(self):
1380 return dict(
1381 codes=self.codes,
1382 tmin=tmin_or_none(self.tmin),
1383 tmax=tmax_or_none(self.tmax),
1384 deltat=self.deltat)
1386 @property
1387 def event_kwargs(self):
1388 return dict(
1389 name=self.codes,
1390 time=self.tmin,
1391 duration=(self.tmax - self.tmin) or None)
1393 @property
1394 def trace_kwargs(self):
1395 network, station, location, channel, extra = self.codes
1397 return dict(
1398 network=network,
1399 station=station,
1400 location=location,
1401 channel=channel,
1402 extra=extra,
1403 tmin=self.tmin,
1404 tmax=self.tmax-self.deltat,
1405 deltat=self.deltat)
1407 @property
1408 def dummy_trace(self):
1409 return DummyTrace(self)
1411 @property
1412 def summary(self):
1413 if self.tmin == self.tmax:
1414 ts = util.time_to_str(self.tmin)
1415 else:
1416 ts = '%s - %s' % (
1417 util.time_to_str(self.tmin),
1418 util.time_to_str(self.tmax))
1420 return ' '.join((
1421 ('%s,' % to_kind(self.kind_id)).ljust(9),
1422 ('%s,' % str(self.codes)).ljust(18),
1423 ts))
1426def make_waveform_nut(**kwargs):
1427 return Nut(kind_id=WAVEFORM, **kwargs)
1430def make_waveform_promise_nut(**kwargs):
1431 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1434def make_station_nut(**kwargs):
1435 return Nut(kind_id=STATION, **kwargs)
1438def make_channel_nut(**kwargs):
1439 return Nut(kind_id=CHANNEL, **kwargs)
1442def make_response_nut(**kwargs):
1443 return Nut(kind_id=RESPONSE, **kwargs)
1446def make_event_nut(**kwargs):
1447 return Nut(kind_id=EVENT, **kwargs)
1450def group_channels(nuts):
1451 by_ansl = {}
1452 for nut in nuts:
1453 if nut.kind_id != CHANNEL:
1454 continue
1456 ansl = nut.codes[:4]
1458 if ansl not in by_ansl:
1459 by_ansl[ansl] = {}
1461 group = by_ansl[ansl]
1463 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1465 if k not in group:
1466 group[k] = set()
1468 group.add(nut.codes[4])
1470 return by_ansl
1473class DummyTrace(object):
1475 def __init__(self, nut):
1476 self.nut = nut
1477 self.codes = nut.codes
1478 self.meta = {}
1480 @property
1481 def tmin(self):
1482 return self.nut.tmin
1484 @property
1485 def tmax(self):
1486 return self.nut.tmax
1488 @property
1489 def deltat(self):
1490 return self.nut.deltat
1492 @property
1493 def nslc_id(self):
1494 return self.codes.nslc
1496 @property
1497 def network(self):
1498 return self.codes.network
1500 @property
1501 def station(self):
1502 return self.codes.station
1504 @property
1505 def location(self):
1506 return self.codes.location
1508 @property
1509 def channel(self):
1510 return self.codes.channel
1512 @property
1513 def extra(self):
1514 return self.codes.extra
1516 def overlaps(self, tmin, tmax):
1517 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1520def duration_to_str(t):
1521 if t > 24*3600:
1522 return '%gd' % (t / (24.*3600.))
1523 elif t > 3600:
1524 return '%gh' % (t / 3600.)
1525 elif t > 60:
1526 return '%gm' % (t / 60.)
1527 else:
1528 return '%gs' % t
1531class Coverage(Object):
1532 '''
1533 Information about times covered by a waveform or other time series data.
1534 '''
1535 kind_id = Int.T(
1536 help='Content type.')
1537 pattern = Codes.T(
1538 help='The codes pattern in the request, which caused this entry to '
1539 'match.')
1540 codes = Codes.T(
1541 help='NSLCE or NSL codes identifier of the time series.')
1542 deltat = Float.T(
1543 help='Sampling interval [s]',
1544 optional=True)
1545 tmin = Timestamp.T(
1546 help='Global start time of time series.',
1547 optional=True)
1548 tmax = Timestamp.T(
1549 help='Global end time of time series.',
1550 optional=True)
1551 changes = List.T(
1552 Tuple.T(2, Any.T()),
1553 help='List of change points, with entries of the form '
1554 '``(time, count)``, where a ``count`` of zero indicates start of '
1555 'a gap, a value of 1 start of normal data coverage and a higher '
1556 'value duplicate or redundant data coverage.')
1558 @classmethod
1559 def from_values(cls, args):
1560 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1561 return cls(
1562 kind_id=kind_id,
1563 pattern=pattern,
1564 codes=codes,
1565 deltat=deltat,
1566 tmin=tmin,
1567 tmax=tmax,
1568 changes=changes)
1570 @property
1571 def summary(self):
1572 ts = '%s - %s,' % (
1573 util.time_to_str(self.tmin),
1574 util.time_to_str(self.tmax))
1576 srate = self.sample_rate
1578 total = self.total
1580 return ' '.join((
1581 ('%s,' % to_kind(self.kind_id)).ljust(9),
1582 ('%s,' % str(self.codes)).ljust(18),
1583 ts,
1584 '%10.3g,' % srate if srate else '',
1585 '%4i' % len(self.changes),
1586 '%s' % duration_to_str(total) if total else 'none'))
1588 @property
1589 def sample_rate(self):
1590 if self.deltat is None:
1591 return None
1592 elif self.deltat == 0.0:
1593 return 0.0
1594 else:
1595 return 1.0 / self.deltat
1597 @property
1598 def labels(self):
1599 srate = self.sample_rate
1600 return (
1601 ('%s' % str(self.codes)),
1602 '%.3g' % srate if srate else '')
1604 @property
1605 def total(self):
1606 total_t = None
1607 for tmin, tmax, _ in self.iter_spans():
1608 total_t = (total_t or 0.0) + (tmax - tmin)
1610 return total_t
1612 def iter_spans(self):
1613 last = None
1614 for (t, count) in self.changes:
1615 if last is not None:
1616 last_t, last_count = last
1617 if last_count > 0:
1618 yield last_t, t, last_count
1620 last = (t, count)
1622 def iter_uncovered_by(self, other):
1623 a = self
1624 b = other
1625 ia = ib = -1
1626 ca = cb = 0
1627 last = None
1628 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1629 if ib + 1 == len(b.changes):
1630 ia += 1
1631 t, ca = a.changes[ia]
1632 elif ia + 1 == len(a.changes):
1633 ib += 1
1634 t, cb = b.changes[ib]
1635 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1636 ia += 1
1637 t, ca = a.changes[ia]
1638 else:
1639 ib += 1
1640 t, cb = b.changes[ib]
1642 if last is not None:
1643 tl, cal, cbl = last
1644 if tl < t and cal > 0 and cbl == 0:
1645 yield tl, t, ia, ib
1647 last = (t, ca, cb)
1649 def iter_uncovered_by_combined(self, other):
1650 ib_last = None
1651 group = None
1652 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1653 if ib_last is None or ib != ib_last:
1654 if group:
1655 yield (group[0][0], group[-1][1])
1657 group = []
1659 group.append((tmin, tmax))
1660 ib_last = ib
1662 if group:
1663 yield (group[0][0], group[-1][1])
1666__all__ = [
1667 'to_codes',
1668 'to_codes_guess',
1669 'to_codes_simple',
1670 'to_kind',
1671 'to_kinds',
1672 'to_kind_id',
1673 'to_kind_ids',
1674 'CodesError',
1675 'Codes',
1676 'CodesNSLCE',
1677 'CodesNSL',
1678 'CodesX',
1679 'Station',
1680 'Channel',
1681 'Sensor',
1682 'Response',
1683 'Nut',
1684 'Coverage',
1685 'WaveformPromise',
1686]