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