1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import sys
7import time
8import logging
9import datetime
10import calendar
11import math
12import copy
14import numpy as num
16from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
17 Unicode, Int, Float, List, Object, Timestamp,
18 ValidationError, TBase, re_tz, Any, Tuple)
19from pyrocko.guts import load_xml # noqa
20from pyrocko.util import hpfloat, time_to_str, get_time_float
22import pyrocko.model
23from pyrocko import util, response
25guts_prefix = 'sx'
27guts_xmlns = 'http://www.fdsn.org/xml/station/1'
29logger = logging.getLogger('pyrocko.io.stationxml')
31conversion = {
32 ('M', 'M'): None,
33 ('M/S', 'M'): response.IntegrationResponse(1),
34 ('M/S**2', 'M'): response.IntegrationResponse(2),
35 ('M', 'M/S'): response.DifferentiationResponse(1),
36 ('M/S', 'M/S'): None,
37 ('M/S**2', 'M/S'): response.IntegrationResponse(1),
38 ('M', 'M/S**2'): response.DifferentiationResponse(2),
39 ('M/S', 'M/S**2'): response.DifferentiationResponse(1),
40 ('M/S**2', 'M/S**2'): None}
43unit_to_quantity = {
44 'M/S': 'velocity',
45 'M': 'displacement',
46 'M/S**2': 'acceleration',
47 'V': 'voltage',
48 'COUNTS': 'counts',
49 'COUNT': 'counts',
50 'PA': 'pressure',
51 'RAD': 'rotation-displacement',
52 'R': 'rotation-displacement',
53 'RAD/S': 'rotation-velocity',
54 'R/S': 'rotation-velocity',
55 'RAD/S**2': 'rotation-acceleration',
56 'R/S**2': 'rotation-acceleration'}
59def to_quantity(unit, context, delivery):
61 if unit is None:
62 return None
64 name = unit.name.upper()
65 if name in unit_to_quantity:
66 return unit_to_quantity[name]
67 else:
68 delivery.log.append((
69 'warning',
70 'Units not supported by Squirrel framework: %s' % (
71 unit.name.upper() + (
72 ' (%s)' % unit.description
73 if unit.description else '')),
74 context))
76 return 'unsupported_quantity(%s)' % unit
79class StationXMLError(Exception):
80 pass
83class Inconsistencies(StationXMLError):
84 pass
87class NoResponseInformation(StationXMLError):
88 pass
91class MultipleResponseInformation(StationXMLError):
92 pass
95class InconsistentResponseInformation(StationXMLError):
96 pass
99class InconsistentChannelLocations(StationXMLError):
100 pass
103class InvalidRecord(StationXMLError):
104 def __init__(self, line):
105 StationXMLError.__init__(self)
106 self._line = line
108 def __str__(self):
109 return 'Invalid record: "%s"' % self._line
112_exceptions = dict(
113 Inconsistencies=Inconsistencies,
114 NoResponseInformation=NoResponseInformation,
115 MultipleResponseInformation=MultipleResponseInformation,
116 InconsistentResponseInformation=InconsistentResponseInformation,
117 InconsistentChannelLocations=InconsistentChannelLocations,
118 InvalidRecord=InvalidRecord,
119 ValueError=ValueError)
122_logs = dict(
123 info=logger.info,
124 warning=logger.warning,
125 error=logger.error)
128class DeliveryError(StationXMLError):
129 pass
132class Delivery(Object):
134 def __init__(self, payload=None, log=None, errors=None, error=None):
135 if payload is None:
136 payload = []
138 if log is None:
139 log = []
141 if errors is None:
142 errors = []
144 if error is not None:
145 errors.append(error)
147 Object.__init__(self, payload=payload, log=log, errors=errors)
149 payload = List.T(Any.T())
150 log = List.T(Tuple.T(3, String.T()))
151 errors = List.T(Tuple.T(3, String.T()))
153 def extend(self, other):
154 self.payload.extend(other.payload)
155 self.log.extend(other.log)
156 self.errors.extend(other.errors)
158 def extend_without_payload(self, other):
159 self.log.extend(other.log)
160 self.errors.extend(other.errors)
161 return other.payload
163 def emit_log(self):
164 for name, message, context in self.log:
165 message = '%s: %s' % (context, message)
166 _logs[name](message)
168 def expect(self, quiet=False):
169 if not quiet:
170 self.emit_log()
172 if self.errors:
173 name, message, context = self.errors[0]
174 if context:
175 message += ' (%s)' % context
177 if len(self.errors) > 1:
178 message += ' Additional errors pending.'
180 raise _exceptions[name](message)
182 return self.payload
184 def expect_one(self, quiet=False):
185 payload = self.expect(quiet=quiet)
186 if len(payload) != 1:
187 raise DeliveryError(
188 'Expected 1 element but got %i.' % len(payload))
190 return payload[0]
193def wrap(s, width=80, indent=4):
194 words = s.split()
195 lines = []
196 t = []
197 n = 0
198 for w in words:
199 if n + len(w) >= width:
200 lines.append(' '.join(t))
201 n = indent
202 t = [' '*(indent-1)]
204 t.append(w)
205 n += len(w) + 1
207 lines.append(' '.join(t))
208 return '\n'.join(lines)
211def same(x, eps=0.0):
212 if any(type(x[0]) != type(r) for r in x):
213 return False
215 if isinstance(x[0], float):
216 return all(abs(r-x[0]) <= eps for r in x)
217 else:
218 return all(r == x[0] for r in x)
221def same_sample_rate(a, b, eps=1.0e-6):
222 return abs(a - b) < min(a, b)*eps
225def evaluate1(resp, f):
226 return resp.evaluate(num.array([f], dtype=float))[0]
229def check_resp(resp, value, frequency, limit_db, prelude, context):
231 try:
232 value_resp = num.abs(evaluate1(resp, frequency))
233 except response.InvalidResponseError as e:
234 return Delivery(log=[(
235 'warning',
236 'Could not check response: %s' % str(e),
237 context)])
239 if value_resp == 0.0:
240 return Delivery(log=[(
241 'warning',
242 '%s\n'
243 ' computed response is zero' % prelude,
244 context)])
246 diff_db = 20.0 * num.log10(value_resp/value)
248 if num.abs(diff_db) > limit_db:
249 return Delivery(log=[(
250 'warning',
251 '%s\n'
252 ' reported value: %g\n'
253 ' computed value: %g\n'
254 ' at frequency [Hz]: %g\n'
255 ' factor: %g\n'
256 ' difference [dB]: %g\n'
257 ' limit [dB]: %g' % (
258 prelude,
259 value,
260 value_resp,
261 frequency,
262 value_resp/value,
263 diff_db,
264 limit_db),
265 context)])
267 return Delivery()
270def tts(t):
271 if t is None:
272 return '<none>'
273 else:
274 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
277def le_open_left(a, b):
278 return a is None or (b is not None and a <= b)
281def le_open_right(a, b):
282 return b is None or (a is not None and a <= b)
285def eq_open(a, b):
286 return (a is None and b is None) \
287 or (a is not None and b is not None and a == b)
290def contains(a, b):
291 return le_open_left(a.start_date, b.start_date) \
292 and le_open_right(b.end_date, a.end_date)
295def find_containing(candidates, node):
296 for candidate in candidates:
297 if contains(candidate, node):
298 return candidate
300 return None
303this_year = time.gmtime()[0]
306class DummyAwareOptionalTimestamp(Object):
307 '''
308 Optional timestamp with support for some common placeholder values.
310 Some StationXML files contain arbitrary placeholder values for open end
311 intervals, like "2100-01-01". Depending on the time range supported by the
312 system, these dates are translated into ``None`` to prevent crashes with
313 this type.
314 '''
315 dummy_for = (hpfloat, float)
316 dummy_for_description = 'time_float'
318 class __T(TBase):
320 def regularize_extra(self, val):
321 time_float = get_time_float()
323 if isinstance(val, datetime.datetime):
324 tt = val.utctimetuple()
325 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
327 elif isinstance(val, datetime.date):
328 tt = val.timetuple()
329 val = time_float(calendar.timegm(tt))
331 elif isinstance(val, str):
332 val = val.strip()
334 tz_offset = 0
336 m = re_tz.search(val)
337 if m:
338 sh = m.group(2)
339 sm = m.group(4)
340 tz_offset = (int(sh)*3600 if sh else 0) \
341 + (int(sm)*60 if sm else 0)
343 val = re_tz.sub('', val)
345 if len(val) > 10 and val[10] == 'T':
346 val = val.replace('T', ' ', 1)
348 try:
349 val = util.str_to_time(val) - tz_offset
351 except util.TimeStrError:
352 year = int(val[:4])
353 if sys.maxsize > 2**32: # if we're on 64bit
354 if year > this_year + 100:
355 return None # StationXML contained a dummy date
357 if year < 1903: # for macOS, 1900-01-01 dummy dates
358 return None
360 else: # 32bit end of time is in 2038
361 if this_year < 2037 and year > 2037 or year < 1903:
362 return None # StationXML contained a dummy date
364 raise
366 elif isinstance(val, (int, float)):
367 val = time_float(val)
369 else:
370 raise ValidationError(
371 '%s: cannot convert "%s" to type %s' % (
372 self.xname(), val, time_float))
374 return val
376 def to_save(self, val):
377 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
378 .rstrip('0').rstrip('.')
380 def to_save_xml(self, val):
381 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
382 .rstrip('0').rstrip('.') + 'Z'
385class Nominal(StringChoice):
386 choices = [
387 'NOMINAL',
388 'CALCULATED']
391class Email(UnicodePattern):
392 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
395class RestrictedStatus(StringChoice):
396 choices = [
397 'open',
398 'closed',
399 'partial']
402class Type(StringChoice):
403 choices = [
404 'TRIGGERED',
405 'CONTINUOUS',
406 'HEALTH',
407 'GEOPHYSICAL',
408 'WEATHER',
409 'FLAG',
410 'SYNTHESIZED',
411 'INPUT',
412 'EXPERIMENTAL',
413 'MAINTENANCE',
414 'BEAM']
416 class __T(StringChoice.T):
417 def validate_extra(self, val):
418 if val not in self.choices:
419 logger.warning(
420 'channel type: "%s" is not a valid choice out of %s' %
421 (val, repr(self.choices)))
424class PzTransferFunction(StringChoice):
425 choices = [
426 'LAPLACE (RADIANS/SECOND)',
427 'LAPLACE (HERTZ)',
428 'DIGITAL (Z-TRANSFORM)']
431class Symmetry(StringChoice):
432 choices = [
433 'NONE',
434 'EVEN',
435 'ODD']
438class CfTransferFunction(StringChoice):
440 class __T(StringChoice.T):
441 def validate(self, val, regularize=False, depth=-1):
442 if regularize:
443 try:
444 val = str(val)
445 except ValueError:
446 raise ValidationError(
447 '%s: cannot convert to string %s' % (self.xname,
448 repr(val)))
450 val = self.dummy_cls.replacements.get(val, val)
452 self.validate_extra(val)
453 return val
455 choices = [
456 'ANALOG (RADIANS/SECOND)',
457 'ANALOG (HERTZ)',
458 'DIGITAL']
460 replacements = {
461 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
462 'ANALOG (HZ)': 'ANALOG (HERTZ)',
463 }
466class Approximation(StringChoice):
467 choices = [
468 'MACLAURIN']
471class PhoneNumber(StringPattern):
472 pattern = '[0-9]+-[0-9]+'
475class Site(Object):
476 '''
477 Description of a site location using name and optional geopolitical
478 boundaries (country, city, etc.).
479 '''
481 name = Unicode.T(default='', xmltagname='Name')
482 description = Unicode.T(optional=True, xmltagname='Description')
483 town = Unicode.T(optional=True, xmltagname='Town')
484 county = Unicode.T(optional=True, xmltagname='County')
485 region = Unicode.T(optional=True, xmltagname='Region')
486 country = Unicode.T(optional=True, xmltagname='Country')
489class ExternalReference(Object):
490 '''
491 This type contains a URI and description for external data that users may
492 want to reference in StationXML.
493 '''
495 uri = String.T(xmltagname='URI')
496 description = Unicode.T(xmltagname='Description')
499class Units(Object):
500 '''
501 A type to document units. Corresponds to SEED blockette 34.
502 '''
504 def __init__(self, name=None, **kwargs):
505 Object.__init__(self, name=name, **kwargs)
507 name = String.T(xmltagname='Name')
508 description = Unicode.T(optional=True, xmltagname='Description')
511class Counter(Int):
512 pass
515class SampleRateRatio(Object):
516 '''
517 Sample rate expressed as number of samples in a number of seconds.
518 '''
520 number_samples = Int.T(xmltagname='NumberSamples')
521 number_seconds = Int.T(xmltagname='NumberSeconds')
524class Gain(Object):
525 '''
526 Complex type for sensitivity and frequency ranges. This complex type can be
527 used to represent both overall sensitivities and individual stage gains.
528 The FrequencyRangeGroup is an optional construct that defines a pass band
529 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
530 is valid within the number of decibels specified in FrequencyDBVariation.
531 '''
533 def __init__(self, value=None, **kwargs):
534 Object.__init__(self, value=value, **kwargs)
536 value = Float.T(optional=True, xmltagname='Value')
537 frequency = Float.T(optional=True, xmltagname='Frequency')
539 def summary(self):
540 return 'gain(%g @ %g)' % (self.value, self.frequency)
543class NumeratorCoefficient(Object):
544 i = Int.T(optional=True, xmlstyle='attribute')
545 value = Float.T(xmlstyle='content')
548class FloatNoUnit(Object):
549 def __init__(self, value=None, **kwargs):
550 Object.__init__(self, value=value, **kwargs)
552 plus_error = Float.T(optional=True, xmlstyle='attribute')
553 minus_error = Float.T(optional=True, xmlstyle='attribute')
554 value = Float.T(xmlstyle='content')
557class FloatWithUnit(FloatNoUnit):
558 unit = String.T(optional=True, xmlstyle='attribute')
561class Equipment(Object):
562 resource_id = String.T(optional=True, xmlstyle='attribute')
563 type = String.T(optional=True, xmltagname='Type')
564 description = Unicode.T(optional=True, xmltagname='Description')
565 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
566 vendor = Unicode.T(optional=True, xmltagname='Vendor')
567 model = Unicode.T(optional=True, xmltagname='Model')
568 serial_number = String.T(optional=True, xmltagname='SerialNumber')
569 installation_date = DummyAwareOptionalTimestamp.T(
570 optional=True,
571 xmltagname='InstallationDate')
572 removal_date = DummyAwareOptionalTimestamp.T(
573 optional=True,
574 xmltagname='RemovalDate')
575 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
578class PhoneNumber(Object):
579 description = Unicode.T(optional=True, xmlstyle='attribute')
580 country_code = Int.T(optional=True, xmltagname='CountryCode')
581 area_code = Int.T(xmltagname='AreaCode')
582 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
585class BaseFilter(Object):
586 '''
587 The BaseFilter is derived by all filters.
588 '''
590 resource_id = String.T(optional=True, xmlstyle='attribute')
591 name = String.T(optional=True, xmlstyle='attribute')
592 description = Unicode.T(optional=True, xmltagname='Description')
593 input_units = Units.T(optional=True, xmltagname='InputUnits')
594 output_units = Units.T(optional=True, xmltagname='OutputUnits')
597class Sensitivity(Gain):
598 '''
599 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
600 construct that defines a pass band in Hertz (FrequencyStart and
601 FrequencyEnd) in which the SensitivityValue is valid within the number of
602 decibels specified in FrequencyDBVariation.
603 '''
605 input_units = Units.T(optional=True, xmltagname='InputUnits')
606 output_units = Units.T(optional=True, xmltagname='OutputUnits')
607 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
608 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
609 frequency_db_variation = Float.T(optional=True,
610 xmltagname='FrequencyDBVariation')
612 def get_pyrocko_response(self):
613 return Delivery(
614 [response.PoleZeroResponse(constant=self.value)])
617class Coefficient(FloatNoUnit):
618 number = Counter.T(optional=True, xmlstyle='attribute')
621class PoleZero(Object):
622 '''
623 Complex numbers used as poles or zeros in channel response.
624 '''
626 number = Int.T(optional=True, xmlstyle='attribute')
627 real = FloatNoUnit.T(xmltagname='Real')
628 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
630 def value(self):
631 return self.real.value + 1J * self.imaginary.value
634class ClockDrift(FloatWithUnit):
635 unit = String.T(default='SECONDS/SAMPLE', optional=True,
636 xmlstyle='attribute') # fixed
639class Second(FloatWithUnit):
640 '''
641 A time value in seconds.
642 '''
644 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
645 # fixed unit
648class Voltage(FloatWithUnit):
649 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
650 # fixed unit
653class Angle(FloatWithUnit):
654 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
655 # fixed unit
658class Azimuth(FloatWithUnit):
659 '''
660 Instrument azimuth, degrees clockwise from North.
661 '''
663 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
664 # fixed unit
667class Dip(FloatWithUnit):
668 '''
669 Instrument dip in degrees down from horizontal. Together azimuth and dip
670 describe the direction of the sensitive axis of the instrument.
671 '''
673 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
674 # fixed unit
677class Distance(FloatWithUnit):
678 '''
679 Extension of FloatWithUnit for distances, elevations, and depths.
680 '''
682 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
683 # NOT fixed unit!
686class Frequency(FloatWithUnit):
687 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
688 # fixed unit
691class SampleRate(FloatWithUnit):
692 '''
693 Sample rate in samples per second.
694 '''
696 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
697 # fixed unit
700class Person(Object):
701 '''
702 Representation of a person's contact information. A person can belong to
703 multiple agencies and have multiple email addresses and phone numbers.
704 '''
706 name_list = List.T(Unicode.T(xmltagname='Name'))
707 agency_list = List.T(Unicode.T(xmltagname='Agency'))
708 email_list = List.T(Email.T(xmltagname='Email'))
709 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
712class FIR(BaseFilter):
713 '''
714 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
715 also commonly documented using the Coefficients element.
716 '''
718 symmetry = Symmetry.T(xmltagname='Symmetry')
719 numerator_coefficient_list = List.T(
720 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
722 def summary(self):
723 return 'fir(%i%s)' % (
724 self.get_ncoefficiencs(),
725 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
727 def get_effective_coefficients(self):
728 b = num.array(
729 [v.value for v in self.numerator_coefficient_list],
730 dtype=float)
732 if self.symmetry == 'ODD':
733 b = num.concatenate((b, b[-2::-1]))
734 elif self.symmetry == 'EVEN':
735 b = num.concatenate((b, b[::-1]))
737 return b
739 def get_effective_symmetry(self):
740 if self.symmetry == 'NONE':
741 b = self.get_effective_coefficients()
742 if num.all(b - b[::-1] == 0):
743 return ['EVEN', 'ODD'][b.size % 2]
745 return self.symmetry
747 def get_ncoefficiencs(self):
748 nf = len(self.numerator_coefficient_list)
749 if self.symmetry == 'ODD':
750 nc = nf * 2 + 1
751 elif self.symmetry == 'EVEN':
752 nc = nf * 2
753 else:
754 nc = nf
756 return nc
758 def estimate_delay(self, deltat):
759 nc = self.get_ncoefficiencs()
760 if nc > 0:
761 return deltat * (nc - 1) / 2.0
762 else:
763 return 0.0
765 def get_pyrocko_response(
766 self, context, deltat, delay_responses, normalization_frequency):
768 context += self.summary()
770 if not self.numerator_coefficient_list:
771 return Delivery([])
773 b = self.get_effective_coefficients()
775 log = []
776 drop_phase = self.get_effective_symmetry() != 'NONE'
778 if not deltat:
779 log.append((
780 'error',
781 'Digital response requires knowledge about sampling '
782 'interval. Response will be unusable.',
783 context))
785 resp = response.DigitalFilterResponse(
786 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
788 if normalization_frequency is not None and deltat is not None:
789 normalization_frequency = 0.0
790 normalization = num.abs(evaluate1(resp, normalization_frequency))
792 if num.abs(normalization - 1.0) > 1e-2:
793 log.append((
794 'warning',
795 'FIR filter coefficients are not normalized. Normalizing '
796 'them. Factor: %g' % normalization, context))
798 resp = response.DigitalFilterResponse(
799 (b/normalization).tolist(), [1.0], deltat,
800 drop_phase=drop_phase)
802 resps = [resp]
804 if not drop_phase:
805 resps.extend(delay_responses)
807 return Delivery(resps, log=log)
810class Coefficients(BaseFilter):
811 '''
812 Response: coefficients for FIR filter. Laplace transforms or IIR filters
813 can be expressed using type as well but the PolesAndZeros should be used
814 instead. Corresponds to SEED blockette 54.
815 '''
817 cf_transfer_function_type = CfTransferFunction.T(
818 xmltagname='CfTransferFunctionType')
819 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
820 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
822 def summary(self):
823 return 'coef_%s(%i,%i%s)' % (
824 'ABC?'[
825 CfTransferFunction.choices.index(
826 self.cf_transfer_function_type)],
827 len(self.numerator_list),
828 len(self.denominator_list),
829 ',sym' if self.is_symmetric_fir else '')
831 def estimate_delay(self, deltat):
832 nc = len(self.numerator_list)
833 if nc > 0:
834 return deltat * (len(self.numerator_list) - 1) / 2.0
835 else:
836 return 0.0
838 def is_symmetric_fir(self):
839 if len(self.denominator_list) != 0:
840 return False
841 b = [v.value for v in self.numerator_list]
842 return b == b[::-1]
844 def get_pyrocko_response(
845 self, context, deltat, delay_responses, normalization_frequency):
847 context += self.summary()
849 factor = 1.0
850 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
851 factor = 2.0*math.pi
853 if not self.numerator_list and not self.denominator_list:
854 return Delivery(payload=[])
856 b = num.array(
857 [v.value*factor for v in self.numerator_list], dtype=float)
859 a = num.array(
860 [1.0] + [v.value*factor for v in self.denominator_list],
861 dtype=float)
863 log = []
864 resps = []
865 if self.cf_transfer_function_type in [
866 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
867 resps.append(response.AnalogFilterResponse(b, a))
869 elif self.cf_transfer_function_type == 'DIGITAL':
870 if not deltat:
871 log.append((
872 'error',
873 'Digital response requires knowledge about sampling '
874 'interval. Response will be unusable.',
875 context))
877 drop_phase = self.is_symmetric_fir()
878 resp = response.DigitalFilterResponse(
879 b, a, deltat or 0.0, drop_phase=drop_phase)
881 if normalization_frequency is not None and deltat is not None:
882 normalization = num.abs(evaluate1(
883 resp, normalization_frequency))
885 if num.abs(normalization - 1.0) > 1e-2:
886 log.append((
887 'warning',
888 'FIR filter coefficients are not normalized. '
889 'Normalizing them. Factor: %g' % normalization,
890 context))
892 resp = response.DigitalFilterResponse(
893 (b/normalization).tolist(), [1.0], deltat,
894 drop_phase=drop_phase)
896 resps.append(resp)
898 if not drop_phase:
899 resps.extend(delay_responses)
901 else:
902 return Delivery(error=(
903 'ValueError',
904 'Unknown transfer function type: %s' % (
905 self.cf_transfer_function_type)))
907 return Delivery(payload=resps, log=log)
910class Latitude(FloatWithUnit):
911 '''
912 Type for latitude coordinate.
913 '''
915 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
916 # fixed unit
917 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
920class Longitude(FloatWithUnit):
921 '''
922 Type for longitude coordinate.
923 '''
925 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
926 # fixed unit
927 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
930class PolesZeros(BaseFilter):
931 '''
932 Response: complex poles and zeros. Corresponds to SEED blockette 53.
933 '''
935 pz_transfer_function_type = PzTransferFunction.T(
936 xmltagname='PzTransferFunctionType')
937 normalization_factor = Float.T(default=1.0,
938 xmltagname='NormalizationFactor')
939 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
940 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
941 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
943 def summary(self):
944 return 'pz_%s(%i,%i)' % (
945 'ABC?'[
946 PzTransferFunction.choices.index(
947 self.pz_transfer_function_type)],
948 len(self.pole_list),
949 len(self.zero_list))
951 def get_pyrocko_response(self, context='', deltat=None):
953 context += self.summary()
955 factor = 1.0
956 cfactor = 1.0
957 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
958 factor = 2. * math.pi
959 cfactor = (2. * math.pi)**(
960 len(self.pole_list) - len(self.zero_list))
962 log = []
963 if self.normalization_factor is None \
964 or self.normalization_factor == 0.0:
966 log.append((
967 'warning',
968 'No pole-zero normalization factor given. '
969 'Assuming a value of 1.0',
970 context))
972 nfactor = 1.0
973 else:
974 nfactor = self.normalization_factor
976 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
977 if not is_digital:
978 resp = response.PoleZeroResponse(
979 constant=nfactor*cfactor,
980 zeros=[z.value()*factor for z in self.zero_list],
981 poles=[p.value()*factor for p in self.pole_list])
982 else:
983 if not deltat:
984 log.append((
985 'error',
986 'Digital response requires knowledge about sampling '
987 'interval. Response will be unusable.',
988 context))
990 resp = response.DigitalPoleZeroResponse(
991 constant=nfactor*cfactor,
992 zeros=[z.value()*factor for z in self.zero_list],
993 poles=[p.value()*factor for p in self.pole_list],
994 deltat=deltat or 0.0)
996 if not self.normalization_frequency.value:
997 log.append((
998 'warning',
999 'Cannot check pole-zero normalization factor: '
1000 'No normalization frequency given.',
1001 context))
1003 else:
1004 if is_digital and not deltat:
1005 log.append((
1006 'warning',
1007 'Cannot check computed vs reported normalization '
1008 'factor without knowing the sampling interval.',
1009 context))
1010 else:
1011 computed_normalization_factor = nfactor / abs(evaluate1(
1012 resp, self.normalization_frequency.value))
1014 db = 20.0 * num.log10(
1015 computed_normalization_factor / nfactor)
1017 if abs(db) > 0.17:
1018 log.append((
1019 'warning',
1020 'Computed and reported normalization factors differ '
1021 'by %g dB: computed: %g, reported: %g' % (
1022 db,
1023 computed_normalization_factor,
1024 nfactor),
1025 context))
1027 return Delivery([resp], log)
1030class ResponseListElement(Object):
1031 frequency = Frequency.T(xmltagname='Frequency')
1032 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1033 phase = Angle.T(xmltagname='Phase')
1036class Polynomial(BaseFilter):
1037 '''
1038 Response: expressed as a polynomial (allows non-linear sensors to be
1039 described). Corresponds to SEED blockette 62. Can be used to describe a
1040 stage of acquisition or a complete system.
1041 '''
1043 approximation_type = Approximation.T(default='MACLAURIN',
1044 xmltagname='ApproximationType')
1045 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1046 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1047 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1048 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1049 maximum_error = Float.T(xmltagname='MaximumError')
1050 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1052 def summary(self):
1053 return 'poly(%i)' % len(self.coefficient_list)
1056class Decimation(Object):
1057 '''
1058 Corresponds to SEED blockette 57.
1059 '''
1061 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1062 factor = Int.T(xmltagname='Factor')
1063 offset = Int.T(xmltagname='Offset')
1064 delay = FloatWithUnit.T(xmltagname='Delay')
1065 correction = FloatWithUnit.T(xmltagname='Correction')
1067 def summary(self):
1068 return 'deci(%i, %g -> %g, %g)' % (
1069 self.factor,
1070 self.input_sample_rate.value,
1071 self.input_sample_rate.value / self.factor,
1072 self.delay.value)
1074 def get_pyrocko_response(self):
1075 if self.delay and self.delay.value != 0.0:
1076 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1077 else:
1078 return Delivery([])
1081class Operator(Object):
1082 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1083 contact_list = List.T(Person.T(xmltagname='Contact'))
1084 web_site = String.T(optional=True, xmltagname='WebSite')
1087class Comment(Object):
1088 '''
1089 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1090 and 59.
1091 '''
1093 id = Counter.T(optional=True, xmlstyle='attribute')
1094 value = Unicode.T(xmltagname='Value')
1095 begin_effective_time = DummyAwareOptionalTimestamp.T(
1096 optional=True,
1097 xmltagname='BeginEffectiveTime')
1098 end_effective_time = DummyAwareOptionalTimestamp.T(
1099 optional=True,
1100 xmltagname='EndEffectiveTime')
1101 author_list = List.T(Person.T(xmltagname='Author'))
1104class ResponseList(BaseFilter):
1105 '''
1106 Response: list of frequency, amplitude and phase values. Corresponds to
1107 SEED blockette 55.
1108 '''
1110 response_list_element_list = List.T(
1111 ResponseListElement.T(xmltagname='ResponseListElement'))
1113 def summary(self):
1114 return 'list(%i)' % len(self.response_list_element_list)
1117class Log(Object):
1118 '''
1119 Container for log entries.
1120 '''
1122 entry_list = List.T(Comment.T(xmltagname='Entry'))
1125class ResponseStage(Object):
1126 '''
1127 This complex type represents channel response and covers SEED blockettes 53
1128 to 56.
1129 '''
1131 number = Counter.T(xmlstyle='attribute')
1132 resource_id = String.T(optional=True, xmlstyle='attribute')
1133 poles_zeros_list = List.T(
1134 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1135 coefficients_list = List.T(
1136 Coefficients.T(optional=True, xmltagname='Coefficients'))
1137 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1138 fir = FIR.T(optional=True, xmltagname='FIR')
1139 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1140 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1141 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1143 def summary(self):
1144 elements = []
1146 for stuff in [
1147 self.poles_zeros_list, self.coefficients_list,
1148 self.response_list, self.fir, self.polynomial,
1149 self.decimation, self.stage_gain]:
1151 if stuff:
1152 if isinstance(stuff, Object):
1153 elements.append(stuff.summary())
1154 else:
1155 elements.extend(obj.summary() for obj in stuff)
1157 return '%i: %s %s -> %s' % (
1158 self.number,
1159 ', '.join(elements),
1160 self.input_units.name.upper() if self.input_units else '?',
1161 self.output_units.name.upper() if self.output_units else '?')
1163 def get_squirrel_response_stage(self, context):
1164 from pyrocko.squirrel.model import ResponseStage
1166 delivery = Delivery()
1167 delivery_pr = self.get_pyrocko_response(context)
1168 log = delivery_pr.log
1169 delivery_pr.log = []
1170 elements = delivery.extend_without_payload(delivery_pr)
1172 delivery.payload = [ResponseStage(
1173 input_quantity=to_quantity(self.input_units, context, delivery),
1174 output_quantity=to_quantity(self.output_units, context, delivery),
1175 input_sample_rate=self.input_sample_rate,
1176 output_sample_rate=self.output_sample_rate,
1177 elements=elements,
1178 log=log)]
1180 return delivery
1182 def get_pyrocko_response(self, context, gain_only=False):
1184 context = context + ', stage %i' % self.number
1186 responses = []
1187 log = []
1188 if self.stage_gain:
1189 normalization_frequency = self.stage_gain.frequency or 0.0
1190 else:
1191 normalization_frequency = 0.0
1193 if not gain_only:
1194 deltat = None
1195 delay_responses = []
1196 if self.decimation:
1197 rate = self.decimation.input_sample_rate.value
1198 if rate > 0.0:
1199 deltat = 1.0 / rate
1200 delivery = self.decimation.get_pyrocko_response()
1201 if delivery.errors:
1202 return delivery
1204 delay_responses = delivery.payload
1205 log.extend(delivery.log)
1207 for pzs in self.poles_zeros_list:
1208 delivery = pzs.get_pyrocko_response(context, deltat)
1209 if delivery.errors:
1210 return delivery
1212 pz_resps = delivery.payload
1213 log.extend(delivery.log)
1214 responses.extend(pz_resps)
1216 # emulate incorrect? evalresp behaviour
1217 if pzs.normalization_frequency != normalization_frequency \
1218 and normalization_frequency != 0.0:
1220 try:
1221 trial = response.MultiplyResponse(pz_resps)
1222 anorm = num.abs(evaluate1(
1223 trial, pzs.normalization_frequency.value))
1224 asens = num.abs(
1225 evaluate1(trial, normalization_frequency))
1227 factor = anorm/asens
1229 if abs(factor - 1.0) > 0.01:
1230 log.append((
1231 'warning',
1232 'PZ normalization frequency (%g) is different '
1233 'from stage gain frequency (%s) -> Emulating '
1234 'possibly incorrect evalresp behaviour. '
1235 'Correction factor: %g' % (
1236 pzs.normalization_frequency.value,
1237 normalization_frequency,
1238 factor),
1239 context))
1241 responses.append(
1242 response.PoleZeroResponse(constant=factor))
1243 except response.InvalidResponseError as e:
1244 log.append((
1245 'warning',
1246 'Could not check response: %s' % str(e),
1247 context))
1249 if len(self.poles_zeros_list) > 1:
1250 log.append((
1251 'warning',
1252 'Multiple poles and zeros records in single response '
1253 'stage.',
1254 context))
1256 for cfs in self.coefficients_list + (
1257 [self.fir] if self.fir else []):
1259 delivery = cfs.get_pyrocko_response(
1260 context, deltat, delay_responses,
1261 normalization_frequency)
1263 if delivery.errors:
1264 return delivery
1266 responses.extend(delivery.payload)
1267 log.extend(delivery.log)
1269 if len(self.coefficients_list) > 1:
1270 log.append((
1271 'warning',
1272 'Multiple filter coefficients lists in single response '
1273 'stage.',
1274 context))
1276 if self.response_list:
1277 log.append((
1278 'warning',
1279 'Unhandled response element of type: ResponseList',
1280 context))
1282 if self.polynomial:
1283 log.append((
1284 'warning',
1285 'Unhandled response element of type: Polynomial',
1286 context))
1288 if self.stage_gain:
1289 responses.append(
1290 response.PoleZeroResponse(constant=self.stage_gain.value))
1292 return Delivery(responses, log)
1294 @property
1295 def input_units(self):
1296 for e in (self.poles_zeros_list + self.coefficients_list +
1297 [self.response_list, self.fir, self.polynomial]):
1298 if e is not None:
1299 return e.input_units
1301 return None
1303 @property
1304 def output_units(self):
1305 for e in (self.poles_zeros_list + self.coefficients_list +
1306 [self.response_list, self.fir, self.polynomial]):
1307 if e is not None:
1308 return e.output_units
1310 return None
1312 @property
1313 def input_sample_rate(self):
1314 if self.decimation:
1315 return self.decimation.input_sample_rate.value
1317 return None
1319 @property
1320 def output_sample_rate(self):
1321 if self.decimation:
1322 return self.decimation.input_sample_rate.value \
1323 / self.decimation.factor
1325 return None
1328class Response(Object):
1329 resource_id = String.T(optional=True, xmlstyle='attribute')
1330 instrument_sensitivity = Sensitivity.T(optional=True,
1331 xmltagname='InstrumentSensitivity')
1332 instrument_polynomial = Polynomial.T(optional=True,
1333 xmltagname='InstrumentPolynomial')
1334 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1336 @property
1337 def output_sample_rate(self):
1338 if self.stage_list:
1339 return self.stage_list[-1].output_sample_rate
1340 else:
1341 return None
1343 def check_sample_rates(self, channel):
1345 if self.stage_list:
1346 sample_rate = None
1348 for stage in self.stage_list:
1349 if stage.decimation:
1350 input_sample_rate = \
1351 stage.decimation.input_sample_rate.value
1353 if sample_rate is not None and not same_sample_rate(
1354 sample_rate, input_sample_rate):
1356 logger.warning(
1357 'Response stage %i has unexpected input sample '
1358 'rate: %g Hz (expected: %g Hz)' % (
1359 stage.number,
1360 input_sample_rate,
1361 sample_rate))
1363 sample_rate = input_sample_rate / stage.decimation.factor
1365 if sample_rate is not None and channel.sample_rate \
1366 and not same_sample_rate(
1367 sample_rate, channel.sample_rate.value):
1369 logger.warning(
1370 'Channel sample rate (%g Hz) does not match sample rate '
1371 'deducted from response stages information (%g Hz).' % (
1372 channel.sample_rate.value,
1373 sample_rate))
1375 def check_units(self):
1377 if self.instrument_sensitivity \
1378 and self.instrument_sensitivity.input_units:
1380 units = self.instrument_sensitivity.input_units.name.upper()
1382 if self.stage_list:
1383 for stage in self.stage_list:
1384 if units and stage.input_units \
1385 and stage.input_units.name.upper() != units:
1387 logger.warning(
1388 'Input units of stage %i (%s) do not match %s (%s).'
1389 % (
1390 stage.number,
1391 units,
1392 'output units of stage %i'
1393 if stage.number == 0
1394 else 'sensitivity input units',
1395 units))
1397 if stage.output_units:
1398 units = stage.output_units.name.upper()
1399 else:
1400 units = None
1402 sout_units = self.instrument_sensitivity.output_units
1403 if self.instrument_sensitivity and sout_units:
1404 if units is not None and units != sout_units.name.upper():
1405 logger.warning(
1406 'Output units of stage %i (%s) do not match %s (%s).'
1407 % (
1408 stage.number,
1409 units,
1410 'sensitivity output units',
1411 sout_units.name.upper()))
1413 def _sensitivity_checkpoints(self, responses, context):
1414 delivery = Delivery()
1416 if self.instrument_sensitivity:
1417 sval = self.instrument_sensitivity.value
1418 sfreq = self.instrument_sensitivity.frequency
1419 if sval is None:
1420 delivery.log.append((
1421 'warning',
1422 'No sensitivity value given.',
1423 context))
1425 elif sval is None:
1426 delivery.log.append((
1427 'warning',
1428 'Reported sensitivity value is zero.',
1429 context))
1431 elif sfreq is None:
1432 delivery.log.append((
1433 'warning',
1434 'Sensitivity frequency not given.',
1435 context))
1437 else:
1438 trial = response.MultiplyResponse(responses)
1440 delivery.extend(
1441 check_resp(
1442 trial, sval, sfreq, 0.1,
1443 'Instrument sensitivity value inconsistent with '
1444 'sensitivity computed from complete response.',
1445 context))
1447 delivery.payload.append(response.FrequencyResponseCheckpoint(
1448 frequency=sfreq, value=sval))
1450 return delivery
1452 def get_squirrel_response(self, context, **kwargs):
1453 from pyrocko.squirrel.model import Response
1455 if self.stage_list:
1456 delivery = Delivery()
1457 for istage, stage in enumerate(self.stage_list):
1458 delivery.extend(stage.get_squirrel_response_stage(context))
1460 checkpoints = []
1461 if not delivery.errors:
1462 all_responses = []
1463 for sq_stage in delivery.payload:
1464 all_responses.extend(sq_stage.elements)
1466 checkpoints.extend(delivery.extend_without_payload(
1467 self._sensitivity_checkpoints(all_responses, context)))
1469 sq_stages = delivery.payload
1470 if sq_stages:
1471 if sq_stages[0].input_quantity is None \
1472 and self.instrument_sensitivity is not None:
1474 sq_stages[0].input_quantity = to_quantity(
1475 self.instrument_sensitivity.input_units,
1476 context, delivery)
1477 sq_stages[-1].output_quantity = to_quantity(
1478 self.instrument_sensitivity.output_units,
1479 context, delivery)
1481 sq_stages = delivery.expect()
1483 return Response(
1484 stages=sq_stages,
1485 log=delivery.log,
1486 checkpoints=checkpoints,
1487 **kwargs)
1489 elif self.instrument_sensitivity:
1490 raise NoResponseInformation(
1491 "Only instrument sensitivity given (won't use it). (%s)."
1492 % context)
1493 else:
1494 raise NoResponseInformation(
1495 'Empty instrument response (%s).'
1496 % context)
1498 def get_pyrocko_response(
1499 self, context, fake_input_units=None, stages=(0, 1)):
1501 delivery = Delivery()
1502 if self.stage_list:
1503 for istage, stage in enumerate(self.stage_list):
1504 delivery.extend(stage.get_pyrocko_response(
1505 context, gain_only=not (
1506 stages is None or stages[0] <= istage < stages[1])))
1508 elif self.instrument_sensitivity:
1509 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1511 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1512 checkpoints = delivery.extend_without_payload(delivery_cp)
1513 if delivery.errors:
1514 return delivery
1516 if fake_input_units is not None:
1517 if not self.instrument_sensitivity or \
1518 self.instrument_sensitivity.input_units is None:
1520 delivery.errors.append((
1521 'NoResponseInformation',
1522 'No input units given, so cannot convert to requested '
1523 'units: %s' % fake_input_units.upper(),
1524 context))
1526 return delivery
1528 input_units = self.instrument_sensitivity.input_units.name.upper()
1530 conresp = None
1531 try:
1532 conresp = conversion[
1533 fake_input_units.upper(), input_units]
1535 except KeyError:
1536 delivery.errors.append((
1537 'NoResponseInformation',
1538 'Cannot convert between units: %s, %s'
1539 % (fake_input_units, input_units),
1540 context))
1542 if conresp is not None:
1543 delivery.payload.append(conresp)
1544 for checkpoint in checkpoints:
1545 checkpoint.value *= num.abs(evaluate1(
1546 conresp, checkpoint.frequency))
1548 delivery.payload = [
1549 response.MultiplyResponse(
1550 delivery.payload,
1551 checkpoints=checkpoints)]
1553 return delivery
1555 @classmethod
1556 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1557 normalization_frequency=1.0):
1558 '''
1559 Convert Pyrocko pole-zero response to StationXML response.
1561 :param presponse: Pyrocko pole-zero response
1562 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1563 :param input_unit: Input unit to be reported in the StationXML
1564 response.
1565 :type input_unit: str
1566 :param output_unit: Output unit to be reported in the StationXML
1567 response.
1568 :type output_unit: str
1569 :param normalization_frequency: Frequency where the normalization
1570 factor for the StationXML response should be computed.
1571 :type normalization_frequency: float
1572 '''
1574 norm_factor = 1.0/float(abs(
1575 evaluate1(presponse, normalization_frequency)
1576 / presponse.constant))
1578 pzs = PolesZeros(
1579 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1580 normalization_factor=norm_factor,
1581 normalization_frequency=Frequency(normalization_frequency),
1582 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1583 imaginary=FloatNoUnit(z.imag))
1584 for z in presponse.zeros],
1585 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1586 imaginary=FloatNoUnit(z.imag))
1587 for z in presponse.poles])
1589 pzs.validate()
1591 stage = ResponseStage(
1592 number=1,
1593 poles_zeros_list=[pzs],
1594 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1596 resp = Response(
1597 instrument_sensitivity=Sensitivity(
1598 value=stage.stage_gain.value,
1599 frequency=normalization_frequency,
1600 input_units=Units(input_unit),
1601 output_units=Units(output_unit)),
1603 stage_list=[stage])
1605 return resp
1608class BaseNode(Object):
1609 '''
1610 A base node type for derivation from: Network, Station and Channel types.
1611 '''
1613 code = String.T(xmlstyle='attribute')
1614 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1615 xmlstyle='attribute')
1616 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1617 xmlstyle='attribute')
1618 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1619 alternate_code = String.T(optional=True, xmlstyle='attribute')
1620 historical_code = String.T(optional=True, xmlstyle='attribute')
1621 description = Unicode.T(optional=True, xmltagname='Description')
1622 comment_list = List.T(Comment.T(xmltagname='Comment'))
1624 def spans(self, *args):
1625 if len(args) == 0:
1626 return True
1627 elif len(args) == 1:
1628 return ((self.start_date is None or
1629 self.start_date <= args[0]) and
1630 (self.end_date is None or
1631 args[0] <= self.end_date))
1633 elif len(args) == 2:
1634 return ((self.start_date is None or
1635 args[1] >= self.start_date) and
1636 (self.end_date is None or
1637 self.end_date >= args[0]))
1640def overlaps(a, b):
1641 return (
1642 a.start_date is None or b.end_date is None
1643 or a.start_date < b.end_date
1644 ) and (
1645 b.start_date is None or a.end_date is None
1646 or b.start_date < a.end_date
1647 )
1650def check_overlaps(node_type_name, codes, nodes):
1651 errors = []
1652 for ia, a in enumerate(nodes):
1653 for b in nodes[ia+1:]:
1654 if overlaps(a, b):
1655 errors.append(
1656 '%s epochs overlap for %s:\n'
1657 ' %s - %s\n %s - %s' % (
1658 node_type_name,
1659 '.'.join(codes),
1660 tts(a.start_date), tts(a.end_date),
1661 tts(b.start_date), tts(b.end_date)))
1663 return errors
1666class Channel(BaseNode):
1667 '''
1668 Equivalent to SEED blockette 52 and parent element for the related the
1669 response blockettes.
1670 '''
1672 location_code = String.T(xmlstyle='attribute')
1673 external_reference_list = List.T(
1674 ExternalReference.T(xmltagname='ExternalReference'))
1675 latitude = Latitude.T(xmltagname='Latitude')
1676 longitude = Longitude.T(xmltagname='Longitude')
1677 elevation = Distance.T(xmltagname='Elevation')
1678 depth = Distance.T(xmltagname='Depth')
1679 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1680 dip = Dip.T(optional=True, xmltagname='Dip')
1681 type_list = List.T(Type.T(xmltagname='Type'))
1682 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1683 sample_rate_ratio = SampleRateRatio.T(optional=True,
1684 xmltagname='SampleRateRatio')
1685 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1686 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1687 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1688 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1689 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1690 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1691 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1692 response = Response.T(optional=True, xmltagname='Response')
1694 @property
1695 def position_values(self):
1696 lat = self.latitude.value
1697 lon = self.longitude.value
1698 elevation = value_or_none(self.elevation)
1699 depth = value_or_none(self.depth)
1700 return lat, lon, elevation, depth
1703class Station(BaseNode):
1704 '''
1705 This type represents a Station epoch. It is common to only have a single
1706 station epoch with the station's creation and termination dates as the
1707 epoch start and end dates.
1708 '''
1710 latitude = Latitude.T(xmltagname='Latitude')
1711 longitude = Longitude.T(xmltagname='Longitude')
1712 elevation = Distance.T(xmltagname='Elevation')
1713 site = Site.T(optional=True, xmltagname='Site')
1714 vault = Unicode.T(optional=True, xmltagname='Vault')
1715 geology = Unicode.T(optional=True, xmltagname='Geology')
1716 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1717 operator_list = List.T(Operator.T(xmltagname='Operator'))
1718 creation_date = DummyAwareOptionalTimestamp.T(
1719 optional=True, xmltagname='CreationDate')
1720 termination_date = DummyAwareOptionalTimestamp.T(
1721 optional=True, xmltagname='TerminationDate')
1722 total_number_channels = Counter.T(
1723 optional=True, xmltagname='TotalNumberChannels')
1724 selected_number_channels = Counter.T(
1725 optional=True, xmltagname='SelectedNumberChannels')
1726 external_reference_list = List.T(
1727 ExternalReference.T(xmltagname='ExternalReference'))
1728 channel_list = List.T(Channel.T(xmltagname='Channel'))
1730 @property
1731 def position_values(self):
1732 lat = self.latitude.value
1733 lon = self.longitude.value
1734 elevation = value_or_none(self.elevation)
1735 return lat, lon, elevation
1738class Network(BaseNode):
1739 '''
1740 This type represents the Network layer, all station metadata is contained
1741 within this element. The official name of the network or other descriptive
1742 information can be included in the Description element. The Network can
1743 contain 0 or more Stations.
1744 '''
1746 total_number_stations = Counter.T(optional=True,
1747 xmltagname='TotalNumberStations')
1748 selected_number_stations = Counter.T(optional=True,
1749 xmltagname='SelectedNumberStations')
1750 station_list = List.T(Station.T(xmltagname='Station'))
1752 @property
1753 def station_code_list(self):
1754 return sorted(set(s.code for s in self.station_list))
1756 @property
1757 def sl_code_list(self):
1758 sls = set()
1759 for station in self.station_list:
1760 for channel in station.channel_list:
1761 sls.add((station.code, channel.location_code))
1763 return sorted(sls)
1765 def summary(self, width=80, indent=4):
1766 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1767 lines = ['%s (%i):' % (self.code, len(sls))]
1768 if sls:
1769 ssls = ['.'.join(x for x in c if x) for c in sls]
1770 w = max(len(x) for x in ssls)
1771 n = (width - indent) / (w+1)
1772 while ssls:
1773 lines.append(
1774 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1776 ssls[:n] = []
1778 return '\n'.join(lines)
1781def value_or_none(x):
1782 if x is not None:
1783 return x.value
1784 else:
1785 return None
1788def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1790 pos = lat, lon, elevation, depth = \
1791 channels[0].position_values
1793 if not all(pos == x.position_values for x in channels):
1794 info = '\n'.join(
1795 ' %s: %s' % (x.code, x.position_values) for
1796 x in channels)
1798 mess = 'encountered inconsistencies in channel ' \
1799 'lat/lon/elevation/depth ' \
1800 'for %s.%s.%s: \n%s' % (nsl + (info,))
1802 if inconsistencies == 'raise':
1803 raise InconsistentChannelLocations(mess)
1805 elif inconsistencies == 'warn':
1806 logger.warning(mess)
1807 logger.warning(' -> using mean values')
1809 apos = num.array([x.position_values for x in channels], dtype=float)
1810 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1811 / num.sum(num.isfinite(apos), axis=0)
1813 groups = {}
1814 for channel in channels:
1815 if channel.code not in groups:
1816 groups[channel.code] = []
1818 groups[channel.code].append(channel)
1820 pchannels = []
1821 for code in sorted(groups.keys()):
1822 data = [
1823 (channel.code, value_or_none(channel.azimuth),
1824 value_or_none(channel.dip))
1825 for channel in groups[code]]
1827 azimuth, dip = util.consistency_merge(
1828 data,
1829 message='channel orientation values differ:',
1830 error=inconsistencies)
1832 pchannels.append(
1833 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1835 return pyrocko.model.Station(
1836 *nsl,
1837 lat=mlat,
1838 lon=mlon,
1839 elevation=mele,
1840 depth=mdep,
1841 channels=pchannels)
1844class FDSNStationXML(Object):
1845 '''
1846 Top-level type for Station XML. Required field are Source (network ID of
1847 the institution sending the message) and one or more Network containers or
1848 one or more Station containers.
1849 '''
1851 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1852 source = String.T(xmltagname='Source')
1853 sender = String.T(optional=True, xmltagname='Sender')
1854 module = String.T(optional=True, xmltagname='Module')
1855 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1856 created = Timestamp.T(optional=True, xmltagname='Created')
1857 network_list = List.T(Network.T(xmltagname='Network'))
1859 xmltagname = 'FDSNStationXML'
1860 guessable_xmlns = [guts_xmlns]
1862 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1863 time=None, timespan=None,
1864 inconsistencies='warn'):
1866 assert inconsistencies in ('raise', 'warn')
1868 if nslcs is not None:
1869 nslcs = set(nslcs)
1871 if nsls is not None:
1872 nsls = set(nsls)
1874 tt = ()
1875 if time is not None:
1876 tt = (time,)
1877 elif timespan is not None:
1878 tt = timespan
1880 pstations = []
1881 for network in self.network_list:
1882 if not network.spans(*tt):
1883 continue
1885 for station in network.station_list:
1886 if not station.spans(*tt):
1887 continue
1889 if station.channel_list:
1890 loc_to_channels = {}
1891 for channel in station.channel_list:
1892 if not channel.spans(*tt):
1893 continue
1895 loc = channel.location_code.strip()
1896 if loc not in loc_to_channels:
1897 loc_to_channels[loc] = []
1899 loc_to_channels[loc].append(channel)
1901 for loc in sorted(loc_to_channels.keys()):
1902 channels = loc_to_channels[loc]
1903 if nslcs is not None:
1904 channels = [channel for channel in channels
1905 if (network.code, station.code, loc,
1906 channel.code) in nslcs]
1908 if not channels:
1909 continue
1911 nsl = network.code, station.code, loc
1912 if nsls is not None and nsl not in nsls:
1913 continue
1915 pstations.append(
1916 pyrocko_station_from_channels(
1917 nsl,
1918 channels,
1919 inconsistencies=inconsistencies))
1920 else:
1921 pstations.append(pyrocko.model.Station(
1922 network.code, station.code, '*',
1923 lat=station.latitude.value,
1924 lon=station.longitude.value,
1925 elevation=value_or_none(station.elevation),
1926 name=station.description or ''))
1928 return pstations
1930 @classmethod
1931 def from_pyrocko_stations(
1932 cls, pyrocko_stations, add_flat_responses_from=None):
1934 '''
1935 Generate :py:class:`FDSNStationXML` from list of
1936 :py:class;`pyrocko.model.Station` instances.
1938 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1939 instances.
1940 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1941 '''
1942 from collections import defaultdict
1943 network_dict = defaultdict(list)
1945 if add_flat_responses_from:
1946 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1947 extra = dict(
1948 response=Response(
1949 instrument_sensitivity=Sensitivity(
1950 value=1.0,
1951 frequency=1.0,
1952 input_units=Units(name=add_flat_responses_from))))
1953 else:
1954 extra = {}
1956 have_offsets = set()
1957 for s in pyrocko_stations:
1959 if s.north_shift != 0.0 or s.east_shift != 0.0:
1960 have_offsets.add(s.nsl())
1962 network, station, location = s.nsl()
1963 channel_list = []
1964 for c in s.channels:
1965 channel_list.append(
1966 Channel(
1967 location_code=location,
1968 code=c.name,
1969 latitude=Latitude(value=s.effective_lat),
1970 longitude=Longitude(value=s.effective_lon),
1971 elevation=Distance(value=s.elevation),
1972 depth=Distance(value=s.depth),
1973 azimuth=Azimuth(value=c.azimuth),
1974 dip=Dip(value=c.dip),
1975 **extra
1976 )
1977 )
1979 network_dict[network].append(
1980 Station(
1981 code=station,
1982 latitude=Latitude(value=s.effective_lat),
1983 longitude=Longitude(value=s.effective_lon),
1984 elevation=Distance(value=s.elevation),
1985 channel_list=channel_list)
1986 )
1988 if have_offsets:
1989 logger.warning(
1990 'StationXML does not support Cartesian offsets in '
1991 'coordinates. Storing effective lat/lon for stations: %s' %
1992 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1994 timestamp = util.to_time_float(time.time())
1995 network_list = []
1996 for k, station_list in network_dict.items():
1998 network_list.append(
1999 Network(
2000 code=k, station_list=station_list,
2001 total_number_stations=len(station_list)))
2003 sxml = FDSNStationXML(
2004 source='from pyrocko stations list', created=timestamp,
2005 network_list=network_list)
2007 sxml.validate()
2008 return sxml
2010 def iter_network_stations(
2011 self, net=None, sta=None, time=None, timespan=None):
2013 tt = ()
2014 if time is not None:
2015 tt = (time,)
2016 elif timespan is not None:
2017 tt = timespan
2019 for network in self.network_list:
2020 if not network.spans(*tt) or (
2021 net is not None and network.code != net):
2022 continue
2024 for station in network.station_list:
2025 if not station.spans(*tt) or (
2026 sta is not None and station.code != sta):
2027 continue
2029 yield (network, station)
2031 def iter_network_station_channels(
2032 self, net=None, sta=None, loc=None, cha=None,
2033 time=None, timespan=None):
2035 if loc is not None:
2036 loc = loc.strip()
2038 tt = ()
2039 if time is not None:
2040 tt = (time,)
2041 elif timespan is not None:
2042 tt = timespan
2044 for network in self.network_list:
2045 if not network.spans(*tt) or (
2046 net is not None and network.code != net):
2047 continue
2049 for station in network.station_list:
2050 if not station.spans(*tt) or (
2051 sta is not None and station.code != sta):
2052 continue
2054 if station.channel_list:
2055 for channel in station.channel_list:
2056 if (not channel.spans(*tt) or
2057 (cha is not None and channel.code != cha) or
2058 (loc is not None and
2059 channel.location_code.strip() != loc)):
2060 continue
2062 yield (network, station, channel)
2064 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2065 time=None, timespan=None):
2067 groups = {}
2068 for network, station, channel in self.iter_network_station_channels(
2069 net, sta, loc, cha, time=time, timespan=timespan):
2071 net = network.code
2072 sta = station.code
2073 cha = channel.code
2074 loc = channel.location_code.strip()
2075 if len(cha) == 3:
2076 bic = cha[:2] # band and intrument code according to SEED
2077 elif len(cha) == 1:
2078 bic = ''
2079 else:
2080 bic = cha
2082 if channel.response and \
2083 channel.response.instrument_sensitivity and \
2084 channel.response.instrument_sensitivity.input_units:
2086 unit = channel.response.instrument_sensitivity\
2087 .input_units.name.upper()
2088 else:
2089 unit = None
2091 bic = (bic, unit)
2093 k = net, sta, loc
2094 if k not in groups:
2095 groups[k] = {}
2097 if bic not in groups[k]:
2098 groups[k][bic] = []
2100 groups[k][bic].append(channel)
2102 for nsl, bic_to_channels in groups.items():
2103 bad_bics = []
2104 for bic, channels in bic_to_channels.items():
2105 sample_rates = []
2106 for channel in channels:
2107 sample_rates.append(channel.sample_rate.value)
2109 if not same(sample_rates):
2110 scs = ','.join(channel.code for channel in channels)
2111 srs = ', '.join('%e' % x for x in sample_rates)
2112 err = 'ignoring channels with inconsistent sampling ' + \
2113 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2115 logger.warning(err)
2116 bad_bics.append(bic)
2118 for bic in bad_bics:
2119 del bic_to_channels[bic]
2121 return groups
2123 def choose_channels(
2124 self,
2125 target_sample_rate=None,
2126 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2127 priority_units=['M/S', 'M/S**2'],
2128 priority_instrument_code=['H', 'L'],
2129 time=None,
2130 timespan=None):
2132 nslcs = {}
2133 for nsl, bic_to_channels in self.get_channel_groups(
2134 time=time, timespan=timespan).items():
2136 useful_bics = []
2137 for bic, channels in bic_to_channels.items():
2138 rate = channels[0].sample_rate.value
2140 if target_sample_rate is not None and \
2141 rate < target_sample_rate*0.99999:
2142 continue
2144 if len(bic[0]) == 2:
2145 if bic[0][0] not in priority_band_code:
2146 continue
2148 if bic[0][1] not in priority_instrument_code:
2149 continue
2151 unit = bic[1]
2153 prio_unit = len(priority_units)
2154 try:
2155 prio_unit = priority_units.index(unit)
2156 except ValueError:
2157 pass
2159 prio_inst = len(priority_instrument_code)
2160 prio_band = len(priority_band_code)
2161 if len(channels[0].code) == 3:
2162 try:
2163 prio_inst = priority_instrument_code.index(
2164 channels[0].code[1])
2165 except ValueError:
2166 pass
2168 try:
2169 prio_band = priority_band_code.index(
2170 channels[0].code[0])
2171 except ValueError:
2172 pass
2174 if target_sample_rate is None:
2175 rate = -rate
2177 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2178 prio_inst, bic))
2180 useful_bics.sort()
2182 for _, _, rate, _, _, bic in useful_bics:
2183 channels = sorted(
2184 bic_to_channels[bic],
2185 key=lambda channel: channel.code)
2187 if channels:
2188 for channel in channels:
2189 nslcs[nsl + (channel.code,)] = channel
2191 break
2193 return nslcs
2195 def get_pyrocko_response(
2196 self, nslc,
2197 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2199 net, sta, loc, cha = nslc
2200 resps = []
2201 for _, _, channel in self.iter_network_station_channels(
2202 net, sta, loc, cha, time=time, timespan=timespan):
2203 resp = channel.response
2204 if resp:
2205 resp.check_sample_rates(channel)
2206 resp.check_units()
2207 resps.append(resp.get_pyrocko_response(
2208 '.'.join(nslc),
2209 fake_input_units=fake_input_units,
2210 stages=stages).expect_one())
2212 if not resps:
2213 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2214 elif len(resps) > 1:
2215 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2217 return resps[0]
2219 @property
2220 def n_code_list(self):
2221 return sorted(set(x.code for x in self.network_list))
2223 @property
2224 def ns_code_list(self):
2225 nss = set()
2226 for network in self.network_list:
2227 for station in network.station_list:
2228 nss.add((network.code, station.code))
2230 return sorted(nss)
2232 @property
2233 def nsl_code_list(self):
2234 nsls = set()
2235 for network in self.network_list:
2236 for station in network.station_list:
2237 for channel in station.channel_list:
2238 nsls.add(
2239 (network.code, station.code, channel.location_code))
2241 return sorted(nsls)
2243 @property
2244 def nslc_code_list(self):
2245 nslcs = set()
2246 for network in self.network_list:
2247 for station in network.station_list:
2248 for channel in station.channel_list:
2249 nslcs.add(
2250 (network.code, station.code, channel.location_code,
2251 channel.code))
2253 return sorted(nslcs)
2255 def summary(self):
2256 lst = [
2257 'number of n codes: %i' % len(self.n_code_list),
2258 'number of ns codes: %i' % len(self.ns_code_list),
2259 'number of nsl codes: %i' % len(self.nsl_code_list),
2260 'number of nslc codes: %i' % len(self.nslc_code_list)
2261 ]
2262 return '\n'.join(lst)
2264 def summary_stages(self):
2265 data = []
2266 for network, station, channel in self.iter_network_station_channels():
2267 nslc = (network.code, station.code, channel.location_code,
2268 channel.code)
2270 stages = []
2271 in_units = '?'
2272 out_units = '?'
2273 if channel.response:
2274 sens = channel.response.instrument_sensitivity
2275 if sens:
2276 in_units = sens.input_units.name.upper()
2277 out_units = sens.output_units.name.upper()
2279 for stage in channel.response.stage_list:
2280 stages.append(stage.summary())
2282 data.append(
2283 (nslc, tts(channel.start_date), tts(channel.end_date),
2284 in_units, out_units, stages))
2286 data.sort()
2288 lst = []
2289 for nslc, stmin, stmax, in_units, out_units, stages in data:
2290 lst.append(' %s: %s - %s, %s -> %s' % (
2291 '.'.join(nslc), stmin, stmax, in_units, out_units))
2292 for stage in stages:
2293 lst.append(' %s' % stage)
2295 return '\n'.join(lst)
2297 def _check_overlaps(self):
2298 by_nslc = {}
2299 for network in self.network_list:
2300 for station in network.station_list:
2301 for channel in station.channel_list:
2302 nslc = (network.code, station.code, channel.location_code,
2303 channel.code)
2304 if nslc not in by_nslc:
2305 by_nslc[nslc] = []
2307 by_nslc[nslc].append(channel)
2309 errors = []
2310 for nslc, channels in by_nslc.items():
2311 errors.extend(check_overlaps('Channel', nslc, channels))
2313 return errors
2315 def check(self):
2316 errors = []
2317 for meth in [self._check_overlaps]:
2318 errors.extend(meth())
2320 if errors:
2321 raise Inconsistencies(
2322 'Inconsistencies found in StationXML:\n '
2323 + '\n '.join(errors))
2326def load_channel_table(stream):
2328 networks = {}
2329 stations = {}
2331 for line in stream:
2332 line = str(line.decode('ascii'))
2333 if line.startswith('#'):
2334 continue
2336 t = line.rstrip().split('|')
2338 if len(t) != 17:
2339 logger.warning('Invalid channel record: %s' % line)
2340 continue
2342 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2343 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2345 try:
2346 scale = float(scale)
2347 except ValueError:
2348 scale = None
2350 try:
2351 scale_freq = float(scale_freq)
2352 except ValueError:
2353 scale_freq = None
2355 try:
2356 depth = float(dep)
2357 except ValueError:
2358 depth = 0.0
2360 try:
2361 azi = float(azi)
2362 dip = float(dip)
2363 except ValueError:
2364 azi = None
2365 dip = None
2367 try:
2368 if net not in networks:
2369 network = Network(code=net)
2370 else:
2371 network = networks[net]
2373 if (net, sta) not in stations:
2374 station = Station(
2375 code=sta, latitude=lat, longitude=lon, elevation=ele)
2377 station.regularize()
2378 else:
2379 station = stations[net, sta]
2381 if scale:
2382 resp = Response(
2383 instrument_sensitivity=Sensitivity(
2384 value=scale,
2385 frequency=scale_freq,
2386 input_units=scale_units))
2387 else:
2388 resp = None
2390 channel = Channel(
2391 code=cha,
2392 location_code=loc.strip(),
2393 latitude=lat,
2394 longitude=lon,
2395 elevation=ele,
2396 depth=depth,
2397 azimuth=azi,
2398 dip=dip,
2399 sensor=Equipment(description=sens),
2400 response=resp,
2401 sample_rate=sample_rate,
2402 start_date=start_date,
2403 end_date=end_date or None)
2405 channel.regularize()
2407 except ValidationError:
2408 raise InvalidRecord(line)
2410 if net not in networks:
2411 networks[net] = network
2413 if (net, sta) not in stations:
2414 stations[net, sta] = station
2415 network.station_list.append(station)
2417 station.channel_list.append(channel)
2419 return FDSNStationXML(
2420 source='created from table input',
2421 created=time.time(),
2422 network_list=sorted(networks.values(), key=lambda x: x.code))
2425def primitive_merge(sxs):
2426 networks = []
2427 for sx in sxs:
2428 networks.extend(sx.network_list)
2430 return FDSNStationXML(
2431 source='merged from different sources',
2432 created=time.time(),
2433 network_list=copy.deepcopy(
2434 sorted(networks, key=lambda x: x.code)))
2437def split_channels(sx):
2438 for nslc in sx.nslc_code_list:
2439 network_list = sx.network_list
2440 network_list_filtered = [
2441 network for network in network_list
2442 if network.code == nslc[0]]
2444 for network in network_list_filtered:
2445 sx.network_list = [network]
2446 station_list = network.station_list
2447 station_list_filtered = [
2448 station for station in station_list
2449 if station.code == nslc[1]]
2451 for station in station_list_filtered:
2452 network.station_list = [station]
2453 channel_list = station.channel_list
2454 station.channel_list = [
2455 channel for channel in channel_list
2456 if (channel.location_code, channel.code) == nslc[2:4]]
2458 yield nslc, copy.deepcopy(sx)
2460 station.channel_list = channel_list
2462 network.station_list = station_list
2464 sx.network_list = network_list
2467if __name__ == '__main__':
2468 from optparse import OptionParser
2470 util.setup_logging('pyrocko.io.stationxml', 'warning')
2472 usage = \
2473 'python -m pyrocko.io.stationxml check|stats|stages ' \
2474 '<filename> [options]'
2476 description = '''Torture StationXML file.'''
2478 parser = OptionParser(
2479 usage=usage,
2480 description=description,
2481 formatter=util.BetterHelpFormatter())
2483 (options, args) = parser.parse_args(sys.argv[1:])
2485 if len(args) != 2:
2486 parser.print_help()
2487 sys.exit(1)
2489 action, path = args
2491 sx = load_xml(filename=path)
2492 if action == 'check':
2493 try:
2494 sx.check()
2495 except Inconsistencies as e:
2496 logger.error(e)
2497 sys.exit(1)
2499 elif action == 'stats':
2500 print(sx.summary())
2502 elif action == 'stages':
2503 print(sx.summary_stages())
2505 else:
2506 parser.print_help()
2507 sys.exit('unknown action: %s' % action)