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[:4]
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 def as_tuple(self):
186 return tuple(self)
188 def replace(self, **kwargs):
189 x = CodesNSLCEBase._replace(self, **kwargs)
190 return g_codes_pool.setdefault(x, x)
193def normalize_nsl(*args, **kwargs):
194 if args and kwargs:
195 raise ValueError('Either *args or **kwargs accepted, not both.')
197 if len(args) == 1:
198 if isinstance(args[0], str):
199 args = tuple(args[0].split('.'))
200 elif isinstance(args[0], tuple):
201 args = args[0]
202 else:
203 raise ValueError('Invalid argument type: %s' % type(args[0]))
205 nargs = len(args)
206 if nargs == 3:
207 t = args
209 elif nargs == 0:
210 d = dict(
211 network='',
212 station='',
213 location='')
215 d.update(kwargs)
216 t = tuple(kwargs.get(k, '') for k in (
217 'network', 'station', 'location'))
219 else:
220 raise CodesError(
221 'Does not match NSL codes pattern: %s' % '.'.join(args))
223 if '.'.join(t).count('.') != 2:
224 raise CodesError(
225 'Codes may not contain a ".": "%s", "%s", "%s"' % t)
227 return t
230CodesNSLBase = namedtuple(
231 'CodesNSLBase', [
232 'network', 'station', 'location'])
235class CodesNSL(CodesNSLBase, Codes):
236 '''
237 Codes denominating a seismic station (NSL).
239 NET.STA.LOC is accepted, slightly different from SEED/StationXML, where
240 LOC is part of the channel. By setting location='*' is possible to get
241 compatible behaviour in most cases.
242 '''
244 __slots__ = ()
245 __hash__ = CodesNSLBase.__hash__
247 as_dict = CodesNSLBase._asdict
249 def __new__(cls, *args, safe_str=None, **kwargs):
250 nargs = len(args)
251 if nargs == 1 and isinstance(args[0], CodesNSL):
252 return args[0]
253 elif nargs == 1 and isinstance(args[0], CodesNSLCE):
254 t = args[0].nsl
255 elif nargs == 1 and isinstance(args[0], CodesX):
256 t = ('*', '*', '*')
257 elif safe_str is not None:
258 t = safe_str.split('.')
259 else:
260 t = normalize_nsl(*args, **kwargs)
262 x = CodesNSLBase.__new__(cls, *t)
263 return g_codes_pool.setdefault(x, x)
265 def __init__(self, *args, **kwargs):
266 Codes.__init__(self)
268 def __str__(self):
269 return '.'.join(self)
271 def __eq__(self, other):
272 if not isinstance(other, CodesNSL):
273 other = CodesNSL(other)
275 return CodesNSLBase.__eq__(self, other)
277 def matches(self, pattern):
278 if not isinstance(pattern, CodesNSL):
279 pattern = CodesNSL(pattern)
281 return match_codes(pattern, self)
283 @property
284 def safe_str(self):
285 return '.'.join(self)
287 @property
288 def ns(self):
289 return self[:2]
291 @property
292 def nsl(self):
293 return self[:3]
295 def as_tuple(self):
296 return tuple(self)
298 def replace(self, **kwargs):
299 x = CodesNSLBase._replace(self, **kwargs)
300 return g_codes_pool.setdefault(x, x)
303CodesXBase = namedtuple(
304 'CodesXBase', [
305 'name'])
308class CodesX(CodesXBase, Codes):
309 '''
310 General purpose codes for anything other than channels or stations.
311 '''
313 __slots__ = ()
314 __hash__ = CodesXBase.__hash__
315 __eq__ = CodesXBase.__eq__
317 as_dict = CodesXBase._asdict
319 def __new__(cls, name='', safe_str=None):
320 if isinstance(name, CodesX):
321 return name
322 elif isinstance(name, (CodesNSLCE, CodesNSL)):
323 name = '*'
324 elif safe_str is not None:
325 name = safe_str
326 else:
327 if '.' in name:
328 raise CodesError('Code may not contain a ".": %s' % name)
330 x = CodesXBase.__new__(cls, name)
331 return g_codes_pool.setdefault(x, x)
333 def __init__(self, *args, **kwargs):
334 Codes.__init__(self)
336 def __str__(self):
337 return '.'.join(self)
339 @property
340 def safe_str(self):
341 return '.'.join(self)
343 @property
344 def ns(self):
345 return self[:2]
347 def as_tuple(self):
348 return tuple(self)
350 def replace(self, **kwargs):
351 x = CodesXBase._replace(self, **kwargs)
352 return g_codes_pool.setdefault(x, x)
355g_codes_patterns = {}
358def match_codes(pattern, codes):
359 spattern = pattern.safe_str
360 scodes = codes.safe_str
361 if spattern not in g_codes_patterns:
362 rpattern = re.compile(fnmatch.translate(spattern), re.I)
363 g_codes_patterns[spattern] = rpattern
365 rpattern = g_codes_patterns[spattern]
366 return bool(rpattern.match(scodes))
369g_content_kinds = [
370 'undefined',
371 'waveform',
372 'station',
373 'channel',
374 'response',
375 'event',
376 'waveform_promise']
379g_codes_classes = [
380 CodesX,
381 CodesNSLCE,
382 CodesNSL,
383 CodesNSLCE,
384 CodesNSLCE,
385 CodesX,
386 CodesNSLCE]
388g_codes_classes_ndot = {
389 0: CodesX,
390 2: CodesNSL,
391 3: CodesNSLCE,
392 4: CodesNSLCE}
395def to_codes_simple(kind_id, codes_safe_str):
396 return g_codes_classes[kind_id](safe_str=codes_safe_str)
399def to_codes(kind_id, obj):
400 return g_codes_classes[kind_id](obj)
403def to_codes_guess(s):
404 try:
405 return g_codes_classes_ndot[s.count('.')](s)
406 except KeyError:
407 raise CodesError('Cannot guess codes type: %s' % s)
410# derived list class to enable detection of already preprocessed codes patterns
411class codes_patterns_list(list):
412 pass
415def codes_patterns_for_kind(kind_id, codes):
416 if isinstance(codes, codes_patterns_list):
417 return codes
419 if isinstance(codes, list):
420 lcodes = codes_patterns_list()
421 for sc in codes:
422 lcodes.extend(codes_patterns_for_kind(kind_id, sc))
424 return lcodes
426 codes = to_codes(kind_id, codes)
428 lcodes = codes_patterns_list()
429 lcodes.append(codes)
430 if kind_id == STATION:
431 lcodes.append(codes.replace(location='[*]'))
433 return lcodes
436g_content_kind_ids = (
437 UNDEFINED, WAVEFORM, STATION, CHANNEL, RESPONSE, EVENT,
438 WAVEFORM_PROMISE) = range(len(g_content_kinds))
441g_tmin, g_tmax = util.get_working_system_time_range()[:2]
444try:
445 g_tmin_queries = max(g_tmin, util.str_to_time_fillup('1900-01-01'))
446except Exception:
447 g_tmin_queries = g_tmin
450def to_kind(kind_id):
451 return g_content_kinds[kind_id]
454def to_kinds(kind_ids):
455 return [g_content_kinds[kind_id] for kind_id in kind_ids]
458def to_kind_id(kind):
459 return g_content_kinds.index(kind)
462def to_kind_ids(kinds):
463 return [g_content_kinds.index(kind) for kind in kinds]
466g_kind_mask_all = 0xff
469def to_kind_mask(kinds):
470 if kinds:
471 return sum(1 << kind_id for kind_id in to_kind_ids(kinds))
472 else:
473 return g_kind_mask_all
476def str_or_none(x):
477 if x is None:
478 return None
479 else:
480 return str(x)
483def float_or_none(x):
484 if x is None:
485 return None
486 else:
487 return float(x)
490def int_or_none(x):
491 if x is None:
492 return None
493 else:
494 return int(x)
497def int_or_g_tmin(x):
498 if x is None:
499 return g_tmin
500 else:
501 return int(x)
504def int_or_g_tmax(x):
505 if x is None:
506 return g_tmax
507 else:
508 return int(x)
511def tmin_or_none(x):
512 if x == g_tmin:
513 return None
514 else:
515 return x
518def tmax_or_none(x):
519 if x == g_tmax:
520 return None
521 else:
522 return x
525def time_or_none_to_str(x):
526 if x is None:
527 return '...'.ljust(17)
528 else:
529 return util.time_to_str(x)
532def codes_to_str_abbreviated(codes, indent=' '):
533 codes = [str(x) for x in codes]
535 if len(codes) > 20:
536 scodes = '\n' + util.ewrap(codes[:10], indent=indent) \
537 + '\n%s[%i more]\n' % (indent, len(codes) - 20) \
538 + util.ewrap(codes[-10:], indent=' ')
539 else:
540 scodes = '\n' + util.ewrap(codes, indent=indent) \
541 if codes else '<none>'
543 return scodes
546g_offset_time_unit_inv = 1000000000
547g_offset_time_unit = 1.0 / g_offset_time_unit_inv
550def tsplit(t):
551 if t is None:
552 return None, 0
554 t = util.to_time_float(t)
555 if type(t) is float:
556 t = round(t, 5)
557 else:
558 t = round(t, 9)
560 seconds = num.floor(t)
561 offset = t - seconds
562 return int(seconds), int(round(offset * g_offset_time_unit_inv))
565def tjoin(seconds, offset):
566 if seconds is None:
567 return None
569 return util.to_time_float(seconds) \
570 + util.to_time_float(offset*g_offset_time_unit)
573tscale_min = 1
574tscale_max = 365 * 24 * 3600 # last edge is one above
575tscale_logbase = 20
577tscale_edges = [tscale_min]
578while True:
579 tscale_edges.append(tscale_edges[-1]*tscale_logbase)
580 if tscale_edges[-1] >= tscale_max:
581 break
584tscale_edges = num.array(tscale_edges)
587def tscale_to_kscale(tscale):
589 # 0 <= x < tscale_edges[1]: 0
590 # tscale_edges[1] <= x < tscale_edges[2]: 1
591 # ...
592 # tscale_edges[len(tscale_edges)-1] <= x: len(tscale_edges)
594 return int(num.searchsorted(tscale_edges, tscale))
597@squirrel_content
598class Station(Location):
599 '''
600 A seismic station.
601 '''
603 codes = CodesNSL.T()
605 tmin = Timestamp.T(optional=True)
606 tmax = Timestamp.T(optional=True)
608 description = Unicode.T(optional=True)
610 def __init__(self, **kwargs):
611 kwargs['codes'] = CodesNSL(kwargs['codes'])
612 Location.__init__(self, **kwargs)
614 @property
615 def time_span(self):
616 return (self.tmin, self.tmax)
618 def get_pyrocko_station(self):
619 from pyrocko import model
620 return model.Station(*self._get_pyrocko_station_args())
622 def _get_pyrocko_station_args(self):
623 return (
624 self.codes.network,
625 self.codes.station,
626 self.codes.location,
627 self.lat,
628 self.lon,
629 self.elevation,
630 self.depth,
631 self.north_shift,
632 self.east_shift)
635@squirrel_content
636class ChannelBase(Location):
637 codes = CodesNSLCE.T()
639 tmin = Timestamp.T(optional=True)
640 tmax = Timestamp.T(optional=True)
642 deltat = Float.T(optional=True)
644 def __init__(self, **kwargs):
645 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
646 Location.__init__(self, **kwargs)
648 @property
649 def time_span(self):
650 return (self.tmin, self.tmax)
652 def _get_sensor_codes(self):
653 return self.codes.replace(
654 channel=self.codes.channel[:-1] + '?')
656 def _get_sensor_args(self):
657 def getattr_rep(k):
658 if k == 'codes':
659 return self._get_sensor_codes()
660 else:
661 return getattr(self, k)
663 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
665 def _get_channel_args(self, component):
666 def getattr_rep(k):
667 if k == 'codes':
668 return self.codes.replace(
669 channel=self.codes.channel[:-1] + component)
670 else:
671 return getattr(self, k)
673 return tuple(getattr_rep(k) for k in ChannelBase.T.propnames)
675 def _get_pyrocko_station_args(self):
676 return (
677 self.codes.network,
678 self.codes.station,
679 self.codes.location,
680 self.lat,
681 self.lon,
682 self.elevation,
683 self.depth,
684 self.north_shift,
685 self.east_shift)
688class Channel(ChannelBase):
689 '''
690 A channel of a seismic station.
691 '''
693 dip = Float.T(optional=True)
694 azimuth = Float.T(optional=True)
696 def get_pyrocko_channel(self):
697 from pyrocko import model
698 return model.Channel(*self._get_pyrocko_channel_args())
700 def _get_pyrocko_channel_args(self):
701 return (
702 self.codes.channel,
703 self.azimuth,
704 self.dip)
706 @property
707 def orientation_enz(self):
708 if None in (self.azimuth, self.dip):
709 return None
711 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
712 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
713 d = math.sin(self.dip*d2r)
714 return mkvec(e, n, -d)
717def cut_intervals(channels):
718 channels = list(channels)
719 times = set()
720 for channel in channels:
721 if channel.tmin is not None:
722 times.add(channel.tmin)
723 if channel.tmax is not None:
724 times.add(channel.tmax)
726 times = sorted(times)
728 if any(channel.tmin is None for channel in channels):
729 times[0:0] = [None]
731 if any(channel.tmax is None for channel in channels):
732 times.append(None)
734 if len(times) <= 2:
735 return channels
737 channels_out = []
738 for channel in channels:
739 for i in range(len(times)-1):
740 tmin = times[i]
741 tmax = times[i+1]
742 if ((channel.tmin is None or (
743 tmin is not None and channel.tmin <= tmin))
744 and (channel.tmax is None or (
745 tmax is not None and tmax <= channel.tmax))):
747 channel_out = clone(channel)
748 channel_out.tmin = tmin
749 channel_out.tmax = tmax
750 channels_out.append(channel_out)
752 return channels_out
755class Sensor(ChannelBase):
756 '''
757 Representation of a channel group.
758 '''
760 channels = List.T(Channel.T())
762 @classmethod
763 def from_channels(cls, channels):
764 groups = defaultdict(list)
765 for channel in channels:
766 groups[channel._get_sensor_codes()].append(channel)
768 channels_cut = []
769 for group in groups.values():
770 channels_cut.extend(cut_intervals(group))
772 groups = defaultdict(list)
773 for channel in channels_cut:
774 groups[channel._get_sensor_args()].append(channel)
776 return [
777 cls(channels=channels,
778 **dict(zip(ChannelBase.T.propnames, args)))
779 for args, channels in groups.items()]
781 def channel_vectors(self):
782 return num.vstack(
783 [channel.orientation_enz for channel in self.channels])
785 def projected_channels(self, system, component_names):
786 return [
787 Channel(
788 azimuth=math.atan2(e, n) * r2d,
789 dip=-math.asin(z) * r2d,
790 **dict(zip(
791 ChannelBase.T.propnames,
792 self._get_channel_args(comp))))
793 for comp, (e, n, z) in zip(component_names, system)]
795 def matrix_to(self, system, epsilon=1e-7):
796 m = num.dot(system, self.channel_vectors().T)
797 m[num.abs(m) < epsilon] = 0.0
798 return m
800 def projection_to(self, system, component_names):
801 return (
802 self.matrix_to(system),
803 self.channels,
804 self.projected_channels(system, component_names))
806 def projection_to_enz(self):
807 return self.projection_to(num.identity(3), 'ENZ')
809 def projection_to_trz(self, source, azimuth=None):
810 if azimuth is not None:
811 assert source is None
812 else:
813 azimuth = source.azibazi_to(self)[1] + 180.
815 return self.projection_to(num.array([
816 [math.cos(azimuth*d2r), -math.sin(azimuth*d2r), 0.],
817 [math.sin(azimuth*d2r), math.cos(azimuth*d2r), 0.],
818 [0., 0., 1.]], dtype=float), 'TRZ')
820 def project_to_enz(self, traces):
821 from pyrocko import trace
823 matrix, in_channels, out_channels = self.projection_to_enz()
825 return trace.project(traces, matrix, in_channels, out_channels)
827 def project_to_trz(self, source, traces, azimuth=None):
828 from pyrocko import trace
830 matrix, in_channels, out_channels = self.projection_to_trz(
831 source, azimuth=azimuth)
833 return trace.project(traces, matrix, in_channels, out_channels)
836observational_quantities = [
837 'acceleration', 'velocity', 'displacement', 'pressure',
838 'rotation-displacement', 'rotation-velocity', 'rotation-acceleration',
839 'temperature']
842technical_quantities = [
843 'voltage', 'counts']
846class QuantityType(StringChoice):
847 '''
848 Choice of observational or technical quantity.
850 SI units are used for all quantities, where applicable.
851 '''
852 choices = observational_quantities + technical_quantities
855class ResponseStage(Object):
856 '''
857 Representation of a response stage.
859 Components of a seismic recording system are represented as a sequence of
860 response stages, e.g. sensor, pre-amplifier, digitizer, digital
861 downsampling.
862 '''
863 input_quantity = QuantityType.T(optional=True)
864 input_sample_rate = Float.T(optional=True)
865 output_quantity = QuantityType.T(optional=True)
866 output_sample_rate = Float.T(optional=True)
867 elements = List.T(FrequencyResponse.T())
868 log = List.T(Tuple.T(3, String.T()))
870 @property
871 def stage_type(self):
872 if self.input_quantity in observational_quantities \
873 and self.output_quantity in observational_quantities:
874 return 'conversion'
876 if self.input_quantity in observational_quantities \
877 and self.output_quantity == 'voltage':
878 return 'sensor'
880 elif self.input_quantity == 'voltage' \
881 and self.output_quantity == 'voltage':
882 return 'electronics'
884 elif self.input_quantity == 'voltage' \
885 and self.output_quantity == 'counts':
886 return 'digitizer'
888 elif self.input_quantity == 'counts' \
889 and self.output_quantity == 'counts' \
890 and self.input_sample_rate != self.output_sample_rate:
891 return 'decimation'
893 elif self.input_quantity in observational_quantities \
894 and self.output_quantity == 'counts':
895 return 'combination'
897 else:
898 return 'unknown'
900 @property
901 def summary(self):
902 irate = self.input_sample_rate
903 orate = self.output_sample_rate
904 factor = None
905 if irate and orate:
906 factor = irate / orate
907 return 'ResponseStage, ' + (
908 '%s%s => %s%s%s' % (
909 self.input_quantity or '?',
910 ' @ %g Hz' % irate if irate else '',
911 self.output_quantity or '?',
912 ' @ %g Hz' % orate if orate else '',
913 ' :%g' % factor if factor else '')
914 )
916 def get_effective(self):
917 return MultiplyResponse(responses=list(self.elements))
920D = 'displacement'
921V = 'velocity'
922A = 'acceleration'
924g_converters = {
925 (V, D): IntegrationResponse(1),
926 (A, D): IntegrationResponse(2),
927 (D, V): DifferentiationResponse(1),
928 (A, V): IntegrationResponse(1),
929 (D, A): DifferentiationResponse(2),
930 (V, A): DifferentiationResponse(1)}
933def response_converters(input_quantity, output_quantity):
934 if input_quantity is None or input_quantity == output_quantity:
935 return []
937 if output_quantity is None:
938 raise ConversionError('Unspecified target quantity.')
940 try:
941 return [g_converters[input_quantity, output_quantity]]
943 except KeyError:
944 raise ConversionError('No rule to convert from "%s" to "%s".' % (
945 input_quantity, output_quantity))
948@squirrel_content
949class Response(Object):
950 '''
951 The instrument response of a seismic station channel.
952 '''
954 codes = CodesNSLCE.T()
955 tmin = Timestamp.T(optional=True)
956 tmax = Timestamp.T(optional=True)
958 stages = List.T(ResponseStage.T())
959 checkpoints = List.T(FrequencyResponseCheckpoint.T())
961 deltat = Float.T(optional=True)
962 log = List.T(Tuple.T(3, String.T()))
964 def __init__(self, **kwargs):
965 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
966 Object.__init__(self, **kwargs)
968 @property
969 def time_span(self):
970 return (self.tmin, self.tmax)
972 @property
973 def nstages(self):
974 return len(self.stages)
976 @property
977 def input_quantity(self):
978 return self.stages[0].input_quantity if self.stages else None
980 @property
981 def output_quantity(self):
982 return self.stages[-1].input_quantity if self.stages else None
984 @property
985 def output_sample_rate(self):
986 return self.stages[-1].output_sample_rate if self.stages else None
988 @property
989 def stages_summary(self):
990 def grouped(xs):
991 xs = list(xs)
992 g = []
993 for i in range(len(xs)):
994 g.append(xs[i])
995 if i+1 < len(xs) and xs[i+1] != xs[i]:
996 yield g
997 g = []
999 if g:
1000 yield g
1002 return '+'.join(
1003 '%s%s' % (g[0], '(%i)' % len(g) if len(g) > 1 else '')
1004 for g in grouped(stage.stage_type for stage in self.stages))
1006 @property
1007 def summary(self):
1008 orate = self.output_sample_rate
1009 return '%s %-16s %s' % (
1010 self.__class__.__name__, self.str_codes, self.str_time_span) \
1011 + ', ' + ', '.join((
1012 '%s => %s' % (
1013 self.input_quantity or '?', self.output_quantity or '?')
1014 + (' @ %g Hz' % orate if orate else ''),
1015 self.stages_summary,
1016 ))
1018 def get_effective(self, input_quantity=None):
1019 try:
1020 elements = response_converters(input_quantity, self.input_quantity)
1021 except ConversionError as e:
1022 raise ConversionError(str(e) + ' (%s)' % self.summary)
1024 elements.extend(
1025 stage.get_effective() for stage in self.stages)
1027 return MultiplyResponse(responses=simplify_responses(elements))
1030@squirrel_content
1031class Event(Object):
1032 '''
1033 A seismic event.
1034 '''
1036 name = String.T(optional=True)
1037 time = Timestamp.T()
1038 duration = Float.T(optional=True)
1040 lat = Float.T()
1041 lon = Float.T()
1042 depth = Float.T(optional=True)
1044 magnitude = Float.T(optional=True)
1046 def get_pyrocko_event(self):
1047 from pyrocko import model
1048 model.Event(
1049 name=self.name,
1050 time=self.time,
1051 lat=self.lat,
1052 lon=self.lon,
1053 depth=self.depth,
1054 magnitude=self.magnitude,
1055 duration=self.duration)
1057 @property
1058 def time_span(self):
1059 return (self.time, self.time)
1062def ehash(s):
1063 return hashlib.sha1(s.encode('utf8')).hexdigest()
1066def random_name(n=8):
1067 return urlsafe_b64encode(urandom(n)).rstrip(b'=').decode('ascii')
1070@squirrel_content
1071class WaveformPromise(Object):
1072 '''
1073 Information about a waveform potentially downloadable from a remote site.
1075 In the Squirrel framework, waveform promises are used to permit download of
1076 selected waveforms from a remote site. They are typically generated by
1077 calls to
1078 :py:meth:`~pyrocko.squirrel.base.Squirrel.update_waveform_promises`.
1079 Waveform promises are inserted and indexed in the database similar to
1080 normal waveforms. When processing a waveform query, e.g. from
1081 :py:meth:`~pyrocko.squirrel.base.Squirrel.get_waveform`, and no local
1082 waveform is available for the queried time span, a matching promise can be
1083 resolved, i.e. an attempt is made to download the waveform from the remote
1084 site. The promise is removed after the download attempt (except when a
1085 network error occurs). This prevents Squirrel from making unnecessary
1086 queries for waveforms missing at the remote site.
1087 '''
1089 codes = CodesNSLCE.T()
1090 tmin = Timestamp.T()
1091 tmax = Timestamp.T()
1093 deltat = Float.T(optional=True)
1095 source_hash = String.T()
1097 def __init__(self, **kwargs):
1098 kwargs['codes'] = CodesNSLCE(kwargs['codes'])
1099 Object.__init__(self, **kwargs)
1101 @property
1102 def time_span(self):
1103 return (self.tmin, self.tmax)
1106class InvalidWaveform(Exception):
1107 pass
1110class WaveformOrder(Object):
1111 '''
1112 Waveform request information.
1113 '''
1115 source_id = String.T()
1116 codes = CodesNSLCE.T()
1117 deltat = Float.T()
1118 tmin = Timestamp.T()
1119 tmax = Timestamp.T()
1120 gaps = List.T(Tuple.T(2, Timestamp.T()))
1122 @property
1123 def client(self):
1124 return self.source_id.split(':')[1]
1126 def describe(self, site='?'):
1127 return '%s:%s %s [%s - %s]' % (
1128 self.client, site, str(self.codes),
1129 util.time_to_str(self.tmin), util.time_to_str(self.tmax))
1131 def validate(self, tr):
1132 if tr.ydata.size == 0:
1133 raise InvalidWaveform(
1134 'waveform with zero data samples')
1136 if tr.deltat != self.deltat:
1137 raise InvalidWaveform(
1138 'incorrect sampling interval - waveform: %g s, '
1139 'meta-data: %g s' % (
1140 tr.deltat, self.deltat))
1142 if not num.all(num.isfinite(tr.ydata)):
1143 raise InvalidWaveform('waveform has NaN values')
1146def order_summary(orders):
1147 codes_list = sorted(set(order.codes for order in orders))
1148 if len(codes_list) > 3:
1149 return '%i order%s: %s - %s' % (
1150 len(orders),
1151 util.plural_s(orders),
1152 str(codes_list[0]),
1153 str(codes_list[1]))
1155 else:
1156 return '%i order%s: %s' % (
1157 len(orders),
1158 util.plural_s(orders),
1159 ', '.join(str(codes) for codes in codes_list))
1162class Nut(Object):
1163 '''
1164 Index entry referencing an elementary piece of content.
1166 So-called *nuts* are used in Pyrocko's Squirrel framework to hold common
1167 meta-information about individual pieces of waveforms, stations, channels,
1168 etc. together with the information where it was found or generated.
1169 '''
1171 file_path = String.T(optional=True)
1172 file_format = String.T(optional=True)
1173 file_mtime = Timestamp.T(optional=True)
1174 file_size = Int.T(optional=True)
1176 file_segment = Int.T(optional=True)
1177 file_element = Int.T(optional=True)
1179 kind_id = Int.T()
1180 codes = Codes.T()
1182 tmin_seconds = Int.T(default=0)
1183 tmin_offset = Int.T(default=0, optional=True)
1184 tmax_seconds = Int.T(default=0)
1185 tmax_offset = Int.T(default=0, optional=True)
1187 deltat = Float.T(default=0.0)
1189 content = Any.T(optional=True)
1190 raw_content = Dict.T(String.T(), Any.T())
1192 content_in_db = False
1194 def __init__(
1195 self,
1196 file_path=None,
1197 file_format=None,
1198 file_mtime=None,
1199 file_size=None,
1200 file_segment=None,
1201 file_element=None,
1202 kind_id=0,
1203 codes=CodesX(''),
1204 tmin_seconds=None,
1205 tmin_offset=0,
1206 tmax_seconds=None,
1207 tmax_offset=0,
1208 deltat=None,
1209 content=None,
1210 raw_content=None,
1211 tmin=None,
1212 tmax=None,
1213 values_nocheck=None):
1215 if values_nocheck is not None:
1216 (self.file_path, self.file_format, self.file_mtime,
1217 self.file_size,
1218 self.file_segment, self.file_element,
1219 self.kind_id, codes_safe_str,
1220 self.tmin_seconds, self.tmin_offset,
1221 self.tmax_seconds, self.tmax_offset,
1222 self.deltat) = values_nocheck
1224 self.codes = to_codes_simple(self.kind_id, codes_safe_str)
1225 self.deltat = self.deltat or None
1226 self.raw_content = {}
1227 self.content = None
1228 else:
1229 if tmin is not None:
1230 tmin_seconds, tmin_offset = tsplit(tmin)
1232 if tmax is not None:
1233 tmax_seconds, tmax_offset = tsplit(tmax)
1235 self.kind_id = int(kind_id)
1236 self.codes = codes
1237 self.tmin_seconds = int_or_g_tmin(tmin_seconds)
1238 self.tmin_offset = int(tmin_offset)
1239 self.tmax_seconds = int_or_g_tmax(tmax_seconds)
1240 self.tmax_offset = int(tmax_offset)
1241 self.deltat = float_or_none(deltat)
1242 self.file_path = str_or_none(file_path)
1243 self.file_segment = int_or_none(file_segment)
1244 self.file_element = int_or_none(file_element)
1245 self.file_format = str_or_none(file_format)
1246 self.file_mtime = float_or_none(file_mtime)
1247 self.file_size = int_or_none(file_size)
1248 self.content = content
1249 if raw_content is None:
1250 self.raw_content = {}
1251 else:
1252 self.raw_content = raw_content
1254 Object.__init__(self, init_props=False)
1256 def __eq__(self, other):
1257 return (isinstance(other, Nut) and
1258 self.equality_values == other.equality_values)
1260 def hash(self):
1261 return ehash(','.join(str(x) for x in self.key))
1263 def __ne__(self, other):
1264 return not (self == other)
1266 def get_io_backend(self):
1267 from . import io
1268 return io.get_backend(self.file_format)
1270 def file_modified(self):
1271 return self.get_io_backend().get_stats(self.file_path) \
1272 != (self.file_mtime, self.file_size)
1274 @property
1275 def dkey(self):
1276 return (self.kind_id, self.tmin_seconds, self.tmin_offset, self.codes)
1278 @property
1279 def key(self):
1280 return (
1281 self.file_path,
1282 self.file_segment,
1283 self.file_element,
1284 self.file_mtime)
1286 @property
1287 def equality_values(self):
1288 return (
1289 self.file_segment, self.file_element,
1290 self.kind_id, self.codes,
1291 self.tmin_seconds, self.tmin_offset,
1292 self.tmax_seconds, self.tmax_offset, self.deltat)
1294 def diff(self, other):
1295 names = [
1296 'file_segment', 'file_element', 'kind_id', 'codes',
1297 'tmin_seconds', 'tmin_offset', 'tmax_seconds', 'tmax_offset',
1298 'deltat']
1300 d = []
1301 for name, a, b in zip(
1302 names, self.equality_values, other.equality_values):
1304 if a != b:
1305 d.append((name, a, b))
1307 return d
1309 @property
1310 def tmin(self):
1311 return tjoin(self.tmin_seconds, self.tmin_offset)
1313 @tmin.setter
1314 def tmin(self, t):
1315 self.tmin_seconds, self.tmin_offset = tsplit(t)
1317 @property
1318 def tmax(self):
1319 return tjoin(self.tmax_seconds, self.tmax_offset)
1321 @tmax.setter
1322 def tmax(self, t):
1323 self.tmax_seconds, self.tmax_offset = tsplit(t)
1325 @property
1326 def kscale(self):
1327 if self.tmin_seconds is None or self.tmax_seconds is None:
1328 return 0
1329 return tscale_to_kscale(self.tmax_seconds - self.tmin_seconds)
1331 @property
1332 def waveform_kwargs(self):
1333 network, station, location, channel, extra = self.codes
1335 return dict(
1336 network=network,
1337 station=station,
1338 location=location,
1339 channel=channel,
1340 extra=extra,
1341 tmin=self.tmin,
1342 tmax=self.tmax,
1343 deltat=self.deltat)
1345 @property
1346 def waveform_promise_kwargs(self):
1347 return dict(
1348 codes=self.codes,
1349 tmin=self.tmin,
1350 tmax=self.tmax,
1351 deltat=self.deltat)
1353 @property
1354 def station_kwargs(self):
1355 network, station, location = self.codes
1356 return dict(
1357 codes=self.codes,
1358 tmin=tmin_or_none(self.tmin),
1359 tmax=tmax_or_none(self.tmax))
1361 @property
1362 def channel_kwargs(self):
1363 network, station, location, channel, extra = self.codes
1364 return dict(
1365 codes=self.codes,
1366 tmin=tmin_or_none(self.tmin),
1367 tmax=tmax_or_none(self.tmax),
1368 deltat=self.deltat)
1370 @property
1371 def response_kwargs(self):
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 event_kwargs(self):
1380 return dict(
1381 name=self.codes,
1382 time=self.tmin,
1383 duration=(self.tmax - self.tmin) or None)
1385 @property
1386 def trace_kwargs(self):
1387 network, station, location, channel, extra = self.codes
1389 return dict(
1390 network=network,
1391 station=station,
1392 location=location,
1393 channel=channel,
1394 extra=extra,
1395 tmin=self.tmin,
1396 tmax=self.tmax-self.deltat,
1397 deltat=self.deltat)
1399 @property
1400 def dummy_trace(self):
1401 return DummyTrace(self)
1403 @property
1404 def summary(self):
1405 if self.tmin == self.tmax:
1406 ts = util.time_to_str(self.tmin)
1407 else:
1408 ts = '%s - %s' % (
1409 util.time_to_str(self.tmin),
1410 util.time_to_str(self.tmax))
1412 return ' '.join((
1413 ('%s,' % to_kind(self.kind_id)).ljust(9),
1414 ('%s,' % str(self.codes)).ljust(18),
1415 ts))
1418def make_waveform_nut(**kwargs):
1419 return Nut(kind_id=WAVEFORM, **kwargs)
1422def make_waveform_promise_nut(**kwargs):
1423 return Nut(kind_id=WAVEFORM_PROMISE, **kwargs)
1426def make_station_nut(**kwargs):
1427 return Nut(kind_id=STATION, **kwargs)
1430def make_channel_nut(**kwargs):
1431 return Nut(kind_id=CHANNEL, **kwargs)
1434def make_response_nut(**kwargs):
1435 return Nut(kind_id=RESPONSE, **kwargs)
1438def make_event_nut(**kwargs):
1439 return Nut(kind_id=EVENT, **kwargs)
1442def group_channels(nuts):
1443 by_ansl = {}
1444 for nut in nuts:
1445 if nut.kind_id != CHANNEL:
1446 continue
1448 ansl = nut.codes[:4]
1450 if ansl not in by_ansl:
1451 by_ansl[ansl] = {}
1453 group = by_ansl[ansl]
1455 k = nut.codes[4][:-1], nut.deltat, nut.tmin, nut.tmax
1457 if k not in group:
1458 group[k] = set()
1460 group.add(nut.codes[4])
1462 return by_ansl
1465class DummyTrace(object):
1467 def __init__(self, nut):
1468 self.nut = nut
1469 self.codes = nut.codes
1470 self.meta = {}
1472 @property
1473 def tmin(self):
1474 return self.nut.tmin
1476 @property
1477 def tmax(self):
1478 return self.nut.tmax
1480 @property
1481 def deltat(self):
1482 return self.nut.deltat
1484 @property
1485 def nslc_id(self):
1486 return self.codes.nslc
1488 @property
1489 def network(self):
1490 return self.codes.network
1492 @property
1493 def station(self):
1494 return self.codes.station
1496 @property
1497 def location(self):
1498 return self.codes.location
1500 @property
1501 def channel(self):
1502 return self.codes.channel
1504 @property
1505 def extra(self):
1506 return self.codes.extra
1508 def overlaps(self, tmin, tmax):
1509 return not (tmax < self.nut.tmin or self.nut.tmax < tmin)
1512def duration_to_str(t):
1513 if t > 24*3600:
1514 return '%gd' % (t / (24.*3600.))
1515 elif t > 3600:
1516 return '%gh' % (t / 3600.)
1517 elif t > 60:
1518 return '%gm' % (t / 60.)
1519 else:
1520 return '%gs' % t
1523class Coverage(Object):
1524 '''
1525 Information about times covered by a waveform or other time series data.
1526 '''
1527 kind_id = Int.T(
1528 help='Content type.')
1529 pattern = Codes.T(
1530 help='The codes pattern in the request, which caused this entry to '
1531 'match.')
1532 codes = Codes.T(
1533 help='NSLCE or NSL codes identifier of the time series.')
1534 deltat = Float.T(
1535 help='Sampling interval [s]',
1536 optional=True)
1537 tmin = Timestamp.T(
1538 help='Global start time of time series.',
1539 optional=True)
1540 tmax = Timestamp.T(
1541 help='Global end time of time series.',
1542 optional=True)
1543 changes = List.T(
1544 Tuple.T(2, Any.T()),
1545 help='List of change points, with entries of the form '
1546 '``(time, count)``, where a ``count`` of zero indicates start of '
1547 'a gap, a value of 1 start of normal data coverage and a higher '
1548 'value duplicate or redundant data coverage.')
1550 @classmethod
1551 def from_values(cls, args):
1552 pattern, codes, deltat, tmin, tmax, changes, kind_id = args
1553 return cls(
1554 kind_id=kind_id,
1555 pattern=pattern,
1556 codes=codes,
1557 deltat=deltat,
1558 tmin=tmin,
1559 tmax=tmax,
1560 changes=changes)
1562 @property
1563 def summary(self):
1564 ts = '%s - %s,' % (
1565 util.time_to_str(self.tmin),
1566 util.time_to_str(self.tmax))
1568 srate = self.sample_rate
1570 total = self.total
1572 return ' '.join((
1573 ('%s,' % to_kind(self.kind_id)).ljust(9),
1574 ('%s,' % str(self.codes)).ljust(18),
1575 ts,
1576 '%10.3g,' % srate if srate else '',
1577 '%4i' % len(self.changes),
1578 '%s' % duration_to_str(total) if total else 'none'))
1580 @property
1581 def sample_rate(self):
1582 if self.deltat is None:
1583 return None
1584 elif self.deltat == 0.0:
1585 return 0.0
1586 else:
1587 return 1.0 / self.deltat
1589 @property
1590 def labels(self):
1591 srate = self.sample_rate
1592 return (
1593 ('%s' % str(self.codes)),
1594 '%.3g' % srate if srate else '')
1596 @property
1597 def total(self):
1598 total_t = None
1599 for tmin, tmax, _ in self.iter_spans():
1600 total_t = (total_t or 0.0) + (tmax - tmin)
1602 return total_t
1604 def iter_spans(self):
1605 last = None
1606 for (t, count) in self.changes:
1607 if last is not None:
1608 last_t, last_count = last
1609 if last_count > 0:
1610 yield last_t, t, last_count
1612 last = (t, count)
1614 def iter_uncovered_by(self, other):
1615 a = self
1616 b = other
1617 ia = ib = -1
1618 ca = cb = 0
1619 last = None
1620 while not (ib + 1 == len(b.changes) and ia + 1 == len(a.changes)):
1621 if ib + 1 == len(b.changes):
1622 ia += 1
1623 t, ca = a.changes[ia]
1624 elif ia + 1 == len(a.changes):
1625 ib += 1
1626 t, cb = b.changes[ib]
1627 elif a.changes[ia+1][0] < b.changes[ib+1][0]:
1628 ia += 1
1629 t, ca = a.changes[ia]
1630 else:
1631 ib += 1
1632 t, cb = b.changes[ib]
1634 if last is not None:
1635 tl, cal, cbl = last
1636 if tl < t and cal > 0 and cbl == 0:
1637 yield tl, t, ia, ib
1639 last = (t, ca, cb)
1641 def iter_uncovered_by_combined(self, other):
1642 ib_last = None
1643 group = None
1644 for tmin, tmax, _, ib in self.iter_uncovered_by(other):
1645 if ib_last is None or ib != ib_last:
1646 if group:
1647 yield (group[0][0], group[-1][1])
1649 group = []
1651 group.append((tmin, tmax))
1652 ib_last = ib
1654 if group:
1655 yield (group[0][0], group[-1][1])
1658__all__ = [
1659 'to_codes',
1660 'to_codes_guess',
1661 'to_codes_simple',
1662 'to_kind',
1663 'to_kinds',
1664 'to_kind_id',
1665 'to_kind_ids',
1666 'CodesError',
1667 'Codes',
1668 'CodesNSLCE',
1669 'CodesNSL',
1670 'CodesX',
1671 'Station',
1672 'Channel',
1673 'Sensor',
1674 'Response',
1675 'Nut',
1676 'Coverage',
1677 'WaveformPromise',
1678]