1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import sys
8import time
9import logging
10import datetime
11import calendar
12import math
13import copy
15import numpy as num
17from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
18 Unicode, Int, Float, List, Object, Timestamp,
19 ValidationError, TBase, re_tz, Any, Tuple)
20from pyrocko.guts import load_xml # noqa
21from pyrocko.util import hpfloat, time_to_str, get_time_float
23import pyrocko.model
24from pyrocko import util, response
26try:
27 newstr = unicode
28except NameError:
29 newstr = str
31guts_prefix = 'sx'
33guts_xmlns = 'http://www.fdsn.org/xml/station/1'
35logger = logging.getLogger('pyrocko.io.stationxml')
37conversion = {
38 ('M', 'M'): None,
39 ('M/S', 'M'): response.IntegrationResponse(1),
40 ('M/S**2', 'M'): response.IntegrationResponse(2),
41 ('M', 'M/S'): response.DifferentiationResponse(1),
42 ('M/S', 'M/S'): None,
43 ('M/S**2', 'M/S'): response.IntegrationResponse(1),
44 ('M', 'M/S**2'): response.DifferentiationResponse(2),
45 ('M/S', 'M/S**2'): response.DifferentiationResponse(1),
46 ('M/S**2', 'M/S**2'): None}
49unit_to_quantity = {
50 'M/S': 'velocity',
51 'M': 'displacement',
52 'M/S**2': 'acceleration',
53 'V': 'voltage',
54 'COUNTS': 'counts',
55 'COUNT': 'counts',
56 'PA': 'pressure',
57 'RAD': 'rotation-displacement',
58 'R': 'rotation-displacement',
59 'RAD/S': 'rotation-velocity',
60 'R/S': 'rotation-velocity',
61 'RAD/S**2': 'rotation-acceleration',
62 'R/S**2': 'rotation-acceleration'}
65def to_quantity(unit, context, delivery):
67 if unit is None:
68 return None
70 name = unit.name.upper()
71 if name in unit_to_quantity:
72 return unit_to_quantity[name]
73 else:
74 delivery.log.append((
75 'warning',
76 'Units not supported by Squirrel framework: %s' % (
77 unit.name.upper() + (
78 ' (%s)' % unit.description
79 if unit.description else '')),
80 context))
82 return 'unsupported_quantity(%s)' % unit
85class StationXMLError(Exception):
86 pass
89class Inconsistencies(StationXMLError):
90 pass
93class NoResponseInformation(StationXMLError):
94 pass
97class MultipleResponseInformation(StationXMLError):
98 pass
101class InconsistentResponseInformation(StationXMLError):
102 pass
105class InconsistentChannelLocations(StationXMLError):
106 pass
109class InvalidRecord(StationXMLError):
110 def __init__(self, line):
111 StationXMLError.__init__(self)
112 self._line = line
114 def __str__(self):
115 return 'Invalid record: "%s"' % self._line
118_exceptions = dict(
119 Inconsistencies=Inconsistencies,
120 NoResponseInformation=NoResponseInformation,
121 MultipleResponseInformation=MultipleResponseInformation,
122 InconsistentResponseInformation=InconsistentResponseInformation,
123 InconsistentChannelLocations=InconsistentChannelLocations,
124 InvalidRecord=InvalidRecord,
125 ValueError=ValueError)
128_logs = dict(
129 info=logger.info,
130 warning=logger.warning,
131 error=logger.error)
134class DeliveryError(StationXMLError):
135 pass
138class Delivery(Object):
140 def __init__(self, payload=None, log=None, errors=None, error=None):
141 if payload is None:
142 payload = []
144 if log is None:
145 log = []
147 if errors is None:
148 errors = []
150 if error is not None:
151 errors.append(error)
153 Object.__init__(self, payload=payload, log=log, errors=errors)
155 payload = List.T(Any.T())
156 log = List.T(Tuple.T(3, String.T()))
157 errors = List.T(Tuple.T(3, String.T()))
159 def extend(self, other):
160 self.payload.extend(other.payload)
161 self.log.extend(other.log)
162 self.errors.extend(other.errors)
164 def extend_without_payload(self, other):
165 self.log.extend(other.log)
166 self.errors.extend(other.errors)
167 return other.payload
169 def emit_log(self):
170 for name, message, context in self.log:
171 message = '%s: %s' % (context, message)
172 _logs[name](message)
174 def expect(self, quiet=False):
175 if not quiet:
176 self.emit_log()
178 if self.errors:
179 name, message, context = self.errors[0]
180 if context:
181 message += ' (%s)' % context
183 if len(self.errors) > 1:
184 message += ' Additional errors pending.'
186 raise _exceptions[name](message)
188 return self.payload
190 def expect_one(self, quiet=False):
191 payload = self.expect(quiet=quiet)
192 if len(payload) != 1:
193 raise DeliveryError(
194 'Expected 1 element but got %i.' % len(payload))
196 return payload[0]
199def wrap(s, width=80, indent=4):
200 words = s.split()
201 lines = []
202 t = []
203 n = 0
204 for w in words:
205 if n + len(w) >= width:
206 lines.append(' '.join(t))
207 n = indent
208 t = [' '*(indent-1)]
210 t.append(w)
211 n += len(w) + 1
213 lines.append(' '.join(t))
214 return '\n'.join(lines)
217def same(x, eps=0.0):
218 if any(type(x[0]) != type(r) for r in x):
219 return False
221 if isinstance(x[0], float):
222 return all(abs(r-x[0]) <= eps for r in x)
223 else:
224 return all(r == x[0] for r in x)
227def same_sample_rate(a, b, eps=1.0e-6):
228 return abs(a - b) < min(a, b)*eps
231def evaluate1(resp, f):
232 return resp.evaluate(num.array([f], dtype=float))[0]
235def check_resp(resp, value, frequency, limit_db, prelude, context):
237 try:
238 value_resp = num.abs(evaluate1(resp, frequency))
239 except response.InvalidResponseError as e:
240 return Delivery(log=[(
241 'warning',
242 'Could not check response: %s' % str(e),
243 context)])
245 if value_resp == 0.0:
246 return Delivery(log=[(
247 'warning',
248 '%s\n'
249 ' computed response is zero' % prelude,
250 context)])
252 diff_db = 20.0 * num.log10(value_resp/value)
254 if num.abs(diff_db) > limit_db:
255 return Delivery(log=[(
256 'warning',
257 '%s\n'
258 ' reported value: %g\n'
259 ' computed value: %g\n'
260 ' at frequency [Hz]: %g\n'
261 ' factor: %g\n'
262 ' difference [dB]: %g\n'
263 ' limit [dB]: %g' % (
264 prelude,
265 value,
266 value_resp,
267 frequency,
268 value_resp/value,
269 diff_db,
270 limit_db),
271 context)])
273 return Delivery()
276def tts(t):
277 if t is None:
278 return '?'
279 else:
280 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
283def le_open_left(a, b):
284 return a is None or (b is not None and a <= b)
287def le_open_right(a, b):
288 return b is None or (a is not None and a <= b)
291def eq_open(a, b):
292 return (a is None and b is None) \
293 or (a is not None and b is not None and a == b)
296def contains(a, b):
297 return le_open_left(a.start_date, b.start_date) \
298 and le_open_right(b.end_date, a.end_date)
301def find_containing(candidates, node):
302 for candidate in candidates:
303 if contains(candidate, node):
304 return candidate
306 return None
309this_year = time.gmtime()[0]
312class DummyAwareOptionalTimestamp(Object):
313 '''
314 Optional timestamp with support for some common placeholder values.
316 Some StationXML files contain arbitrary placeholder values for open end
317 intervals, like "2100-01-01". Depending on the time range supported by the
318 system, these dates are translated into ``None`` to prevent crashes with
319 this type.
320 '''
321 dummy_for = (hpfloat, float)
322 dummy_for_description = 'time_float'
324 class __T(TBase):
326 def regularize_extra(self, val):
327 time_float = get_time_float()
329 if isinstance(val, datetime.datetime):
330 tt = val.utctimetuple()
331 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
333 elif isinstance(val, datetime.date):
334 tt = val.timetuple()
335 val = time_float(calendar.timegm(tt))
337 elif isinstance(val, (str, newstr)):
338 val = val.strip()
340 tz_offset = 0
342 m = re_tz.search(val)
343 if m:
344 sh = m.group(2)
345 sm = m.group(4)
346 tz_offset = (int(sh)*3600 if sh else 0) \
347 + (int(sm)*60 if sm else 0)
349 val = re_tz.sub('', val)
351 if len(val) > 10 and val[10] == 'T':
352 val = val.replace('T', ' ', 1)
354 try:
355 val = util.str_to_time(val) - tz_offset
357 except util.TimeStrError:
358 year = int(val[:4])
359 if sys.maxsize > 2**32: # if we're on 64bit
360 if year > this_year + 100:
361 return None # StationXML contained a dummy date
363 if year < 1903: # for macOS, 1900-01-01 dummy dates
364 return None
366 else: # 32bit end of time is in 2038
367 if this_year < 2037 and year > 2037 or year < 1903:
368 return None # StationXML contained a dummy date
370 raise
372 elif isinstance(val, (int, float)):
373 val = time_float(val)
375 else:
376 raise ValidationError(
377 '%s: cannot convert "%s" to type %s' % (
378 self.xname(), val, time_float))
380 return val
382 def to_save(self, val):
383 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
384 .rstrip('0').rstrip('.')
386 def to_save_xml(self, val):
387 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
388 .rstrip('0').rstrip('.') + 'Z'
391class Nominal(StringChoice):
392 choices = [
393 'NOMINAL',
394 'CALCULATED']
397class Email(UnicodePattern):
398 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
401class RestrictedStatus(StringChoice):
402 choices = [
403 'open',
404 'closed',
405 'partial']
408class Type(StringChoice):
409 choices = [
410 'TRIGGERED',
411 'CONTINUOUS',
412 'HEALTH',
413 'GEOPHYSICAL',
414 'WEATHER',
415 'FLAG',
416 'SYNTHESIZED',
417 'INPUT',
418 'EXPERIMENTAL',
419 'MAINTENANCE',
420 'BEAM']
422 class __T(StringChoice.T):
423 def validate_extra(self, val):
424 if val not in self.choices:
425 logger.warning(
426 'channel type: "%s" is not a valid choice out of %s' %
427 (val, repr(self.choices)))
430class PzTransferFunction(StringChoice):
431 choices = [
432 'LAPLACE (RADIANS/SECOND)',
433 'LAPLACE (HERTZ)',
434 'DIGITAL (Z-TRANSFORM)']
437class Symmetry(StringChoice):
438 choices = [
439 'NONE',
440 'EVEN',
441 'ODD']
444class CfTransferFunction(StringChoice):
446 class __T(StringChoice.T):
447 def validate(self, val, regularize=False, depth=-1):
448 if regularize:
449 try:
450 val = str(val)
451 except ValueError:
452 raise ValidationError(
453 '%s: cannot convert to string %s' % (self.xname,
454 repr(val)))
456 val = self.dummy_cls.replacements.get(val, val)
458 self.validate_extra(val)
459 return val
461 choices = [
462 'ANALOG (RADIANS/SECOND)',
463 'ANALOG (HERTZ)',
464 'DIGITAL']
466 replacements = {
467 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
468 'ANALOG (HZ)': 'ANALOG (HERTZ)',
469 }
472class Approximation(StringChoice):
473 choices = [
474 'MACLAURIN']
477class PhoneNumber(StringPattern):
478 pattern = '[0-9]+-[0-9]+'
481class Site(Object):
482 '''
483 Description of a site location using name and optional geopolitical
484 boundaries (country, city, etc.).
485 '''
487 name = Unicode.T(default='', xmltagname='Name')
488 description = Unicode.T(optional=True, xmltagname='Description')
489 town = Unicode.T(optional=True, xmltagname='Town')
490 county = Unicode.T(optional=True, xmltagname='County')
491 region = Unicode.T(optional=True, xmltagname='Region')
492 country = Unicode.T(optional=True, xmltagname='Country')
495class ExternalReference(Object):
496 '''
497 This type contains a URI and description for external data that users may
498 want to reference in StationXML.
499 '''
501 uri = String.T(xmltagname='URI')
502 description = Unicode.T(xmltagname='Description')
505class Units(Object):
506 '''
507 A type to document units. Corresponds to SEED blockette 34.
508 '''
510 def __init__(self, name=None, **kwargs):
511 Object.__init__(self, name=name, **kwargs)
513 name = String.T(xmltagname='Name')
514 description = Unicode.T(optional=True, xmltagname='Description')
517class Counter(Int):
518 pass
521class SampleRateRatio(Object):
522 '''
523 Sample rate expressed as number of samples in a number of seconds.
524 '''
526 number_samples = Int.T(xmltagname='NumberSamples')
527 number_seconds = Int.T(xmltagname='NumberSeconds')
530class Gain(Object):
531 '''
532 Complex type for sensitivity and frequency ranges. This complex type can be
533 used to represent both overall sensitivities and individual stage gains.
534 The FrequencyRangeGroup is an optional construct that defines a pass band
535 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
536 is valid within the number of decibels specified in FrequencyDBVariation.
537 '''
539 def __init__(self, value=None, **kwargs):
540 Object.__init__(self, value=value, **kwargs)
542 value = Float.T(optional=True, xmltagname='Value')
543 frequency = Float.T(optional=True, xmltagname='Frequency')
545 def summary(self):
546 return 'gain(%g @ %g)' % (self.value, self.frequency)
549class NumeratorCoefficient(Object):
550 i = Int.T(optional=True, xmlstyle='attribute')
551 value = Float.T(xmlstyle='content')
554class FloatNoUnit(Object):
555 def __init__(self, value=None, **kwargs):
556 Object.__init__(self, value=value, **kwargs)
558 plus_error = Float.T(optional=True, xmlstyle='attribute')
559 minus_error = Float.T(optional=True, xmlstyle='attribute')
560 value = Float.T(xmlstyle='content')
563class FloatWithUnit(FloatNoUnit):
564 unit = String.T(optional=True, xmlstyle='attribute')
567class Equipment(Object):
568 resource_id = String.T(optional=True, xmlstyle='attribute')
569 type = String.T(optional=True, xmltagname='Type')
570 description = Unicode.T(optional=True, xmltagname='Description')
571 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
572 vendor = Unicode.T(optional=True, xmltagname='Vendor')
573 model = Unicode.T(optional=True, xmltagname='Model')
574 serial_number = String.T(optional=True, xmltagname='SerialNumber')
575 installation_date = DummyAwareOptionalTimestamp.T(
576 optional=True,
577 xmltagname='InstallationDate')
578 removal_date = DummyAwareOptionalTimestamp.T(
579 optional=True,
580 xmltagname='RemovalDate')
581 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
584class PhoneNumber(Object):
585 description = Unicode.T(optional=True, xmlstyle='attribute')
586 country_code = Int.T(optional=True, xmltagname='CountryCode')
587 area_code = Int.T(xmltagname='AreaCode')
588 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
591class BaseFilter(Object):
592 '''
593 The BaseFilter is derived by all filters.
594 '''
596 resource_id = String.T(optional=True, xmlstyle='attribute')
597 name = String.T(optional=True, xmlstyle='attribute')
598 description = Unicode.T(optional=True, xmltagname='Description')
599 input_units = Units.T(optional=True, xmltagname='InputUnits')
600 output_units = Units.T(optional=True, xmltagname='OutputUnits')
603class Sensitivity(Gain):
604 '''
605 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
606 construct that defines a pass band in Hertz (FrequencyStart and
607 FrequencyEnd) in which the SensitivityValue is valid within the number of
608 decibels specified in FrequencyDBVariation.
609 '''
611 input_units = Units.T(optional=True, xmltagname='InputUnits')
612 output_units = Units.T(optional=True, xmltagname='OutputUnits')
613 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
614 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
615 frequency_db_variation = Float.T(optional=True,
616 xmltagname='FrequencyDBVariation')
618 def get_pyrocko_response(self):
619 return Delivery(
620 [response.PoleZeroResponse(constant=self.value)])
623class Coefficient(FloatNoUnit):
624 number = Counter.T(optional=True, xmlstyle='attribute')
627class PoleZero(Object):
628 '''
629 Complex numbers used as poles or zeros in channel response.
630 '''
632 number = Int.T(optional=True, xmlstyle='attribute')
633 real = FloatNoUnit.T(xmltagname='Real')
634 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
636 def value(self):
637 return self.real.value + 1J * self.imaginary.value
640class ClockDrift(FloatWithUnit):
641 unit = String.T(default='SECONDS/SAMPLE', optional=True,
642 xmlstyle='attribute') # fixed
645class Second(FloatWithUnit):
646 '''
647 A time value in seconds.
648 '''
650 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
651 # fixed unit
654class Voltage(FloatWithUnit):
655 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
656 # fixed unit
659class Angle(FloatWithUnit):
660 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
661 # fixed unit
664class Azimuth(FloatWithUnit):
665 '''
666 Instrument azimuth, degrees clockwise from North.
667 '''
669 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
670 # fixed unit
673class Dip(FloatWithUnit):
674 '''
675 Instrument dip in degrees down from horizontal. Together azimuth and dip
676 describe the direction of the sensitive axis of the instrument.
677 '''
679 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
680 # fixed unit
683class Distance(FloatWithUnit):
684 '''
685 Extension of FloatWithUnit for distances, elevations, and depths.
686 '''
688 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
689 # NOT fixed unit!
692class Frequency(FloatWithUnit):
693 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
694 # fixed unit
697class SampleRate(FloatWithUnit):
698 '''
699 Sample rate in samples per second.
700 '''
702 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
703 # fixed unit
706class Person(Object):
707 '''
708 Representation of a person's contact information. A person can belong to
709 multiple agencies and have multiple email addresses and phone numbers.
710 '''
712 name_list = List.T(Unicode.T(xmltagname='Name'))
713 agency_list = List.T(Unicode.T(xmltagname='Agency'))
714 email_list = List.T(Email.T(xmltagname='Email'))
715 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
718class FIR(BaseFilter):
719 '''
720 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
721 also commonly documented using the Coefficients element.
722 '''
724 symmetry = Symmetry.T(xmltagname='Symmetry')
725 numerator_coefficient_list = List.T(
726 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
728 def summary(self):
729 return 'fir(%i%s)' % (
730 self.get_ncoefficiencs(),
731 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
733 def get_effective_coefficients(self):
734 b = num.array(
735 [v.value for v in self.numerator_coefficient_list],
736 dtype=float)
738 if self.symmetry == 'ODD':
739 b = num.concatenate((b, b[-2::-1]))
740 elif self.symmetry == 'EVEN':
741 b = num.concatenate((b, b[::-1]))
743 return b
745 def get_effective_symmetry(self):
746 if self.symmetry == 'NONE':
747 b = self.get_effective_coefficients()
748 if num.all(b - b[::-1] == 0):
749 return ['EVEN', 'ODD'][b.size % 2]
751 return self.symmetry
753 def get_ncoefficiencs(self):
754 nf = len(self.numerator_coefficient_list)
755 if self.symmetry == 'ODD':
756 nc = nf * 2 + 1
757 elif self.symmetry == 'EVEN':
758 nc = nf * 2
759 else:
760 nc = nf
762 return nc
764 def estimate_delay(self, deltat):
765 nc = self.get_ncoefficiencs()
766 if nc > 0:
767 return deltat * (nc - 1) / 2.0
768 else:
769 return 0.0
771 def get_pyrocko_response(
772 self, context, deltat, delay_responses, normalization_frequency):
774 context += self.summary()
776 if not self.numerator_coefficient_list:
777 return Delivery([])
779 b = self.get_effective_coefficients()
781 log = []
782 drop_phase = self.get_effective_symmetry() != 'NONE'
784 if not deltat:
785 log.append((
786 'error',
787 'Digital response requires knowledge about sampling '
788 'interval. Response will be unusable.',
789 context))
791 resp = response.DigitalFilterResponse(
792 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
794 if normalization_frequency is not None and deltat is not None:
795 normalization_frequency = 0.0
796 normalization = num.abs(evaluate1(resp, normalization_frequency))
798 if num.abs(normalization - 1.0) > 1e-2:
799 log.append((
800 'warning',
801 'FIR filter coefficients are not normalized. Normalizing '
802 'them. Factor: %g' % normalization, context))
804 resp = response.DigitalFilterResponse(
805 (b/normalization).tolist(), [1.0], deltat,
806 drop_phase=drop_phase)
808 resps = [resp]
810 if not drop_phase:
811 resps.extend(delay_responses)
813 return Delivery(resps, log=log)
816class Coefficients(BaseFilter):
817 '''
818 Response: coefficients for FIR filter. Laplace transforms or IIR filters
819 can be expressed using type as well but the PolesAndZeros should be used
820 instead. Corresponds to SEED blockette 54.
821 '''
823 cf_transfer_function_type = CfTransferFunction.T(
824 xmltagname='CfTransferFunctionType')
825 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
826 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
828 def summary(self):
829 return 'coef_%s(%i,%i%s)' % (
830 'ABC?'[
831 CfTransferFunction.choices.index(
832 self.cf_transfer_function_type)],
833 len(self.numerator_list),
834 len(self.denominator_list),
835 ',sym' if self.is_symmetric_fir else '')
837 def estimate_delay(self, deltat):
838 nc = len(self.numerator_list)
839 if nc > 0:
840 return deltat * (len(self.numerator_list) - 1) / 2.0
841 else:
842 return 0.0
844 def is_symmetric_fir(self):
845 if len(self.denominator_list) != 0:
846 return False
847 b = [v.value for v in self.numerator_list]
848 return b == b[::-1]
850 def get_pyrocko_response(
851 self, context, deltat, delay_responses, normalization_frequency):
853 context += self.summary()
855 factor = 1.0
856 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
857 factor = 2.0*math.pi
859 if not self.numerator_list and not self.denominator_list:
860 return Delivery(payload=[])
862 b = num.array(
863 [v.value*factor for v in self.numerator_list], dtype=float)
865 a = num.array(
866 [1.0] + [v.value*factor for v in self.denominator_list],
867 dtype=float)
869 log = []
870 resps = []
871 if self.cf_transfer_function_type in [
872 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
873 resps.append(response.AnalogFilterResponse(b, a))
875 elif self.cf_transfer_function_type == 'DIGITAL':
876 if not deltat:
877 log.append((
878 'error',
879 'Digital response requires knowledge about sampling '
880 'interval. Response will be unusable.',
881 context))
883 drop_phase = self.is_symmetric_fir()
884 resp = response.DigitalFilterResponse(
885 b, a, deltat or 0.0, drop_phase=drop_phase)
887 if normalization_frequency is not None and deltat is not None:
888 normalization = num.abs(evaluate1(
889 resp, normalization_frequency))
891 if num.abs(normalization - 1.0) > 1e-2:
892 log.append((
893 'warning',
894 'FIR filter coefficients are not normalized. '
895 'Normalizing them. Factor: %g' % normalization,
896 context))
898 resp = response.DigitalFilterResponse(
899 (b/normalization).tolist(), [1.0], deltat,
900 drop_phase=drop_phase)
902 resps.append(resp)
904 if not drop_phase:
905 resps.extend(delay_responses)
907 else:
908 return Delivery(error=(
909 'ValueError',
910 'Unknown transfer function type: %s' % (
911 self.cf_transfer_function_type)))
913 return Delivery(payload=resps, log=log)
916class Latitude(FloatWithUnit):
917 '''
918 Type for latitude coordinate.
919 '''
921 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
922 # fixed unit
923 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
926class Longitude(FloatWithUnit):
927 '''
928 Type for longitude coordinate.
929 '''
931 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
932 # fixed unit
933 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
936class PolesZeros(BaseFilter):
937 '''
938 Response: complex poles and zeros. Corresponds to SEED blockette 53.
939 '''
941 pz_transfer_function_type = PzTransferFunction.T(
942 xmltagname='PzTransferFunctionType')
943 normalization_factor = Float.T(default=1.0,
944 xmltagname='NormalizationFactor')
945 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
946 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
947 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
949 def summary(self):
950 return 'pz_%s(%i,%i)' % (
951 'ABC?'[
952 PzTransferFunction.choices.index(
953 self.pz_transfer_function_type)],
954 len(self.pole_list),
955 len(self.zero_list))
957 def get_pyrocko_response(self, context='', deltat=None):
959 context += self.summary()
961 factor = 1.0
962 cfactor = 1.0
963 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
964 factor = 2. * math.pi
965 cfactor = (2. * math.pi)**(
966 len(self.pole_list) - len(self.zero_list))
968 log = []
969 if self.normalization_factor is None \
970 or self.normalization_factor == 0.0:
972 log.append((
973 'warning',
974 'No pole-zero normalization factor given. '
975 'Assuming a value of 1.0',
976 context))
978 nfactor = 1.0
979 else:
980 nfactor = self.normalization_factor
982 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
983 if not is_digital:
984 resp = response.PoleZeroResponse(
985 constant=nfactor*cfactor,
986 zeros=[z.value()*factor for z in self.zero_list],
987 poles=[p.value()*factor for p in self.pole_list])
988 else:
989 if not deltat:
990 log.append((
991 'error',
992 'Digital response requires knowledge about sampling '
993 'interval. Response will be unusable.',
994 context))
996 resp = response.DigitalPoleZeroResponse(
997 constant=nfactor*cfactor,
998 zeros=[z.value()*factor for z in self.zero_list],
999 poles=[p.value()*factor for p in self.pole_list],
1000 deltat=deltat or 0.0)
1002 if not self.normalization_frequency.value:
1003 log.append((
1004 'warning',
1005 'Cannot check pole-zero normalization factor: '
1006 'No normalization frequency given.',
1007 context))
1009 else:
1010 if is_digital and not deltat:
1011 log.append((
1012 'warning',
1013 'Cannot check computed vs reported normalization '
1014 'factor without knowing the sampling interval.',
1015 context))
1016 else:
1017 computed_normalization_factor = nfactor / abs(evaluate1(
1018 resp, self.normalization_frequency.value))
1020 db = 20.0 * num.log10(
1021 computed_normalization_factor / nfactor)
1023 if abs(db) > 0.17:
1024 log.append((
1025 'warning',
1026 'Computed and reported normalization factors differ '
1027 'by %g dB: computed: %g, reported: %g' % (
1028 db,
1029 computed_normalization_factor,
1030 nfactor),
1031 context))
1033 return Delivery([resp], log)
1036class ResponseListElement(Object):
1037 frequency = Frequency.T(xmltagname='Frequency')
1038 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1039 phase = Angle.T(xmltagname='Phase')
1042class Polynomial(BaseFilter):
1043 '''
1044 Response: expressed as a polynomial (allows non-linear sensors to be
1045 described). Corresponds to SEED blockette 62. Can be used to describe a
1046 stage of acquisition or a complete system.
1047 '''
1049 approximation_type = Approximation.T(default='MACLAURIN',
1050 xmltagname='ApproximationType')
1051 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1052 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1053 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1054 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1055 maximum_error = Float.T(xmltagname='MaximumError')
1056 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1058 def summary(self):
1059 return 'poly(%i)' % len(self.coefficient_list)
1062class Decimation(Object):
1063 '''
1064 Corresponds to SEED blockette 57.
1065 '''
1067 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1068 factor = Int.T(xmltagname='Factor')
1069 offset = Int.T(xmltagname='Offset')
1070 delay = FloatWithUnit.T(xmltagname='Delay')
1071 correction = FloatWithUnit.T(xmltagname='Correction')
1073 def summary(self):
1074 return 'deci(%i, %g -> %g, %g)' % (
1075 self.factor,
1076 self.input_sample_rate.value,
1077 self.input_sample_rate.value / self.factor,
1078 self.delay.value)
1080 def get_pyrocko_response(self):
1081 if self.delay and self.delay.value != 0.0:
1082 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1083 else:
1084 return Delivery([])
1087class Operator(Object):
1088 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1089 contact_list = List.T(Person.T(xmltagname='Contact'))
1090 web_site = String.T(optional=True, xmltagname='WebSite')
1093class Comment(Object):
1094 '''
1095 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1096 and 59.
1097 '''
1099 id = Counter.T(optional=True, xmlstyle='attribute')
1100 value = Unicode.T(xmltagname='Value')
1101 begin_effective_time = DummyAwareOptionalTimestamp.T(
1102 optional=True,
1103 xmltagname='BeginEffectiveTime')
1104 end_effective_time = DummyAwareOptionalTimestamp.T(
1105 optional=True,
1106 xmltagname='EndEffectiveTime')
1107 author_list = List.T(Person.T(xmltagname='Author'))
1110class ResponseList(BaseFilter):
1111 '''
1112 Response: list of frequency, amplitude and phase values. Corresponds to
1113 SEED blockette 55.
1114 '''
1116 response_list_element_list = List.T(
1117 ResponseListElement.T(xmltagname='ResponseListElement'))
1119 def summary(self):
1120 return 'list(%i)' % len(self.response_list_element_list)
1123class Log(Object):
1124 '''
1125 Container for log entries.
1126 '''
1128 entry_list = List.T(Comment.T(xmltagname='Entry'))
1131class ResponseStage(Object):
1132 '''
1133 This complex type represents channel response and covers SEED blockettes 53
1134 to 56.
1135 '''
1137 number = Counter.T(xmlstyle='attribute')
1138 resource_id = String.T(optional=True, xmlstyle='attribute')
1139 poles_zeros_list = List.T(
1140 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1141 coefficients_list = List.T(
1142 Coefficients.T(optional=True, xmltagname='Coefficients'))
1143 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1144 fir = FIR.T(optional=True, xmltagname='FIR')
1145 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1146 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1147 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1149 def summary(self):
1150 elements = []
1152 for stuff in [
1153 self.poles_zeros_list, self.coefficients_list,
1154 self.response_list, self.fir, self.polynomial,
1155 self.decimation, self.stage_gain]:
1157 if stuff:
1158 if isinstance(stuff, Object):
1159 elements.append(stuff.summary())
1160 else:
1161 elements.extend(obj.summary() for obj in stuff)
1163 return '%i: %s %s -> %s' % (
1164 self.number,
1165 ', '.join(elements),
1166 self.input_units.name.upper() if self.input_units else '?',
1167 self.output_units.name.upper() if self.output_units else '?')
1169 def get_squirrel_response_stage(self, context):
1170 from pyrocko.squirrel.model import ResponseStage
1172 delivery = Delivery()
1173 delivery_pr = self.get_pyrocko_response(context)
1174 log = delivery_pr.log
1175 delivery_pr.log = []
1176 elements = delivery.extend_without_payload(delivery_pr)
1178 delivery.payload = [ResponseStage(
1179 input_quantity=to_quantity(self.input_units, context, delivery),
1180 output_quantity=to_quantity(self.output_units, context, delivery),
1181 input_sample_rate=self.input_sample_rate,
1182 output_sample_rate=self.output_sample_rate,
1183 elements=elements,
1184 log=log)]
1186 return delivery
1188 def get_pyrocko_response(self, context, gain_only=False):
1190 context = context + ', stage %i' % self.number
1192 responses = []
1193 log = []
1194 if self.stage_gain:
1195 normalization_frequency = self.stage_gain.frequency or 0.0
1196 else:
1197 normalization_frequency = 0.0
1199 if not gain_only:
1200 deltat = None
1201 delay_responses = []
1202 if self.decimation:
1203 rate = self.decimation.input_sample_rate.value
1204 if rate > 0.0:
1205 deltat = 1.0 / rate
1206 delivery = self.decimation.get_pyrocko_response()
1207 if delivery.errors:
1208 return delivery
1210 delay_responses = delivery.payload
1211 log.extend(delivery.log)
1213 for pzs in self.poles_zeros_list:
1214 delivery = pzs.get_pyrocko_response(context, deltat)
1215 if delivery.errors:
1216 return delivery
1218 pz_resps = delivery.payload
1219 log.extend(delivery.log)
1220 responses.extend(pz_resps)
1222 # emulate incorrect? evalresp behaviour
1223 if pzs.normalization_frequency != normalization_frequency \
1224 and normalization_frequency != 0.0:
1226 try:
1227 trial = response.MultiplyResponse(pz_resps)
1228 anorm = num.abs(evaluate1(
1229 trial, pzs.normalization_frequency.value))
1230 asens = num.abs(
1231 evaluate1(trial, normalization_frequency))
1233 factor = anorm/asens
1235 if abs(factor - 1.0) > 0.01:
1236 log.append((
1237 'warning',
1238 'PZ normalization frequency (%g) is different '
1239 'from stage gain frequency (%s) -> Emulating '
1240 'possibly incorrect evalresp behaviour. '
1241 'Correction factor: %g' % (
1242 pzs.normalization_frequency.value,
1243 normalization_frequency,
1244 factor),
1245 context))
1247 responses.append(
1248 response.PoleZeroResponse(constant=factor))
1249 except response.InvalidResponseError as e:
1250 log.append((
1251 'warning',
1252 'Could not check response: %s' % str(e),
1253 context))
1255 if len(self.poles_zeros_list) > 1:
1256 log.append((
1257 'warning',
1258 'Multiple poles and zeros records in single response '
1259 'stage.',
1260 context))
1262 for cfs in self.coefficients_list + (
1263 [self.fir] if self.fir else []):
1265 delivery = cfs.get_pyrocko_response(
1266 context, deltat, delay_responses,
1267 normalization_frequency)
1269 if delivery.errors:
1270 return delivery
1272 responses.extend(delivery.payload)
1273 log.extend(delivery.log)
1275 if len(self.coefficients_list) > 1:
1276 log.append((
1277 'warning',
1278 'Multiple filter coefficients lists in single response '
1279 'stage.',
1280 context))
1282 if self.response_list:
1283 log.append((
1284 'warning',
1285 'Unhandled response element of type: ResponseList',
1286 context))
1288 if self.polynomial:
1289 log.append((
1290 'warning',
1291 'Unhandled response element of type: Polynomial',
1292 context))
1294 if self.stage_gain:
1295 responses.append(
1296 response.PoleZeroResponse(constant=self.stage_gain.value))
1298 return Delivery(responses, log)
1300 @property
1301 def input_units(self):
1302 for e in (self.poles_zeros_list + self.coefficients_list +
1303 [self.response_list, self.fir, self.polynomial]):
1304 if e is not None:
1305 return e.input_units
1307 return None
1309 @property
1310 def output_units(self):
1311 for e in (self.poles_zeros_list + self.coefficients_list +
1312 [self.response_list, self.fir, self.polynomial]):
1313 if e is not None:
1314 return e.output_units
1316 return None
1318 @property
1319 def input_sample_rate(self):
1320 if self.decimation:
1321 return self.decimation.input_sample_rate.value
1323 return None
1325 @property
1326 def output_sample_rate(self):
1327 if self.decimation:
1328 return self.decimation.input_sample_rate.value \
1329 / self.decimation.factor
1331 return None
1334class Response(Object):
1335 resource_id = String.T(optional=True, xmlstyle='attribute')
1336 instrument_sensitivity = Sensitivity.T(optional=True,
1337 xmltagname='InstrumentSensitivity')
1338 instrument_polynomial = Polynomial.T(optional=True,
1339 xmltagname='InstrumentPolynomial')
1340 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1342 def check_sample_rates(self, channel):
1344 if self.stage_list:
1345 sample_rate = None
1347 for stage in self.stage_list:
1348 if stage.decimation:
1349 input_sample_rate = \
1350 stage.decimation.input_sample_rate.value
1352 if sample_rate is not None and not same_sample_rate(
1353 sample_rate, input_sample_rate):
1355 logger.warning(
1356 'Response stage %i has unexpected input sample '
1357 'rate: %g Hz (expected: %g Hz)' % (
1358 stage.number,
1359 input_sample_rate,
1360 sample_rate))
1362 sample_rate = input_sample_rate / stage.decimation.factor
1364 if sample_rate is not None and channel.sample_rate \
1365 and not same_sample_rate(
1366 sample_rate, channel.sample_rate.value):
1368 logger.warning(
1369 'Channel sample rate (%g Hz) does not match sample rate '
1370 'deducted from response stages information (%g Hz).' % (
1371 channel.sample_rate.value,
1372 sample_rate))
1374 def check_units(self):
1376 if self.instrument_sensitivity \
1377 and self.instrument_sensitivity.input_units:
1379 units = self.instrument_sensitivity.input_units.name.upper()
1381 if self.stage_list:
1382 for stage in self.stage_list:
1383 if units and stage.input_units \
1384 and stage.input_units.name.upper() != units:
1386 logger.warning(
1387 'Input units of stage %i (%s) do not match %s (%s).'
1388 % (
1389 stage.number,
1390 units,
1391 'output units of stage %i'
1392 if stage.number == 0
1393 else 'sensitivity input units',
1394 units))
1396 if stage.output_units:
1397 units = stage.output_units.name.upper()
1398 else:
1399 units = None
1401 sout_units = self.instrument_sensitivity.output_units
1402 if self.instrument_sensitivity and sout_units:
1403 if units is not None and units != sout_units.name.upper():
1404 logger.warning(
1405 'Output units of stage %i (%s) do not match %s (%s).'
1406 % (
1407 stage.number,
1408 units,
1409 'sensitivity output units',
1410 sout_units.name.upper()))
1412 def _sensitivity_checkpoints(self, responses, context):
1413 delivery = Delivery()
1415 if self.instrument_sensitivity:
1416 sval = self.instrument_sensitivity.value
1417 sfreq = self.instrument_sensitivity.frequency
1418 if sval is None:
1419 delivery.log.append((
1420 'warning',
1421 'No sensitivity value given.',
1422 context))
1424 elif sval is None:
1425 delivery.log.append((
1426 'warning',
1427 'Reported sensitivity value is zero.',
1428 context))
1430 elif sfreq is None:
1431 delivery.log.append((
1432 'warning',
1433 'Sensitivity frequency not given.',
1434 context))
1436 else:
1437 trial = response.MultiplyResponse(responses)
1439 delivery.extend(
1440 check_resp(
1441 trial, sval, sfreq, 0.1,
1442 'Instrument sensitivity value inconsistent with '
1443 'sensitivity computed from complete response',
1444 context))
1446 delivery.payload.append(response.FrequencyResponseCheckpoint(
1447 frequency=sfreq, value=sval))
1449 return delivery
1451 def get_squirrel_response(self, context, **kwargs):
1452 from pyrocko.squirrel.model import Response
1454 if self.stage_list:
1455 delivery = Delivery()
1456 for istage, stage in enumerate(self.stage_list):
1457 delivery.extend(stage.get_squirrel_response_stage(context))
1459 checkpoints = []
1460 if not delivery.errors:
1461 all_responses = []
1462 for sq_stage in delivery.payload:
1463 all_responses.extend(sq_stage.elements)
1465 checkpoints.extend(delivery.extend_without_payload(
1466 self._sensitivity_checkpoints(all_responses, context)))
1468 sq_stages = delivery.payload
1469 if sq_stages:
1470 if sq_stages[0].input_quantity is None \
1471 and self.instrument_sensitivity is not None:
1473 sq_stages[0].input_quantity = to_quantity(
1474 self.instrument_sensitivity.input_units,
1475 context, delivery)
1476 sq_stages[-1].output_quantity = to_quantity(
1477 self.instrument_sensitivity.output_units,
1478 context, delivery)
1480 sq_stages = delivery.expect()
1482 return Response(
1483 stages=sq_stages,
1484 log=delivery.log,
1485 checkpoints=checkpoints,
1486 **kwargs)
1488 elif self.instrument_sensitivity:
1489 raise NoResponseInformation(
1490 'Only instrument sensitivity given (won\'t use it). (%s).'
1491 % context)
1492 else:
1493 raise NoResponseInformation(
1494 'Empty instrument response (%s).'
1495 % context)
1497 def get_pyrocko_response(
1498 self, context, fake_input_units=None, stages=(0, 1)):
1500 delivery = Delivery()
1501 if self.stage_list:
1502 for istage, stage in enumerate(self.stage_list):
1503 delivery.extend(stage.get_pyrocko_response(
1504 context, gain_only=not (
1505 stages is None or stages[0] <= istage < stages[1])))
1507 elif self.instrument_sensitivity:
1508 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1510 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1511 checkpoints = delivery.extend_without_payload(delivery_cp)
1512 if delivery.errors:
1513 return delivery
1515 if fake_input_units is not None:
1516 if not self.instrument_sensitivity or \
1517 self.instrument_sensitivity.input_units is None:
1519 delivery.errors.append((
1520 'NoResponseInformation',
1521 'No input units given, so cannot convert to requested '
1522 'units: %s' % fake_input_units.upper(),
1523 context))
1525 return delivery
1527 input_units = self.instrument_sensitivity.input_units.name.upper()
1529 conresp = None
1530 try:
1531 conresp = conversion[
1532 fake_input_units.upper(), input_units]
1534 except KeyError:
1535 delivery.errors.append((
1536 'NoResponseInformation',
1537 'Cannot convert between units: %s, %s'
1538 % (fake_input_units, input_units),
1539 context))
1541 if conresp is not None:
1542 delivery.payload.append(conresp)
1543 for checkpoint in checkpoints:
1544 checkpoint.value *= num.abs(evaluate1(
1545 conresp, checkpoint.frequency))
1547 delivery.payload = [
1548 response.MultiplyResponse(
1549 delivery.payload,
1550 checkpoints=checkpoints)]
1552 return delivery
1554 @classmethod
1555 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1556 normalization_frequency=1.0):
1557 '''
1558 Convert Pyrocko pole-zero response to StationXML response.
1560 :param presponse: Pyrocko pole-zero response
1561 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1562 :param input_unit: Input unit to be reported in the StationXML
1563 response.
1564 :type input_unit: str
1565 :param output_unit: Output unit to be reported in the StationXML
1566 response.
1567 :type output_unit: str
1568 :param normalization_frequency: Frequency where the normalization
1569 factor for the StationXML response should be computed.
1570 :type normalization_frequency: float
1571 '''
1573 norm_factor = 1.0/float(abs(
1574 evaluate1(presponse, normalization_frequency)
1575 / presponse.constant))
1577 pzs = PolesZeros(
1578 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1579 normalization_factor=norm_factor,
1580 normalization_frequency=Frequency(normalization_frequency),
1581 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1582 imaginary=FloatNoUnit(z.imag))
1583 for z in presponse.zeros],
1584 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1585 imaginary=FloatNoUnit(z.imag))
1586 for z in presponse.poles])
1588 pzs.validate()
1590 stage = ResponseStage(
1591 number=1,
1592 poles_zeros_list=[pzs],
1593 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1595 resp = Response(
1596 instrument_sensitivity=Sensitivity(
1597 value=stage.stage_gain.value,
1598 frequency=normalization_frequency,
1599 input_units=Units(input_unit),
1600 output_units=Units(output_unit)),
1602 stage_list=[stage])
1604 return resp
1607class BaseNode(Object):
1608 '''
1609 A base node type for derivation from: Network, Station and Channel types.
1610 '''
1612 code = String.T(xmlstyle='attribute')
1613 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1614 xmlstyle='attribute')
1615 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1616 xmlstyle='attribute')
1617 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1618 alternate_code = String.T(optional=True, xmlstyle='attribute')
1619 historical_code = String.T(optional=True, xmlstyle='attribute')
1620 description = Unicode.T(optional=True, xmltagname='Description')
1621 comment_list = List.T(Comment.T(xmltagname='Comment'))
1623 def spans(self, *args):
1624 if len(args) == 0:
1625 return True
1626 elif len(args) == 1:
1627 return ((self.start_date is None or
1628 self.start_date <= args[0]) and
1629 (self.end_date is None or
1630 args[0] <= self.end_date))
1632 elif len(args) == 2:
1633 return ((self.start_date is None or
1634 args[1] >= self.start_date) and
1635 (self.end_date is None or
1636 self.end_date >= args[0]))
1639def overlaps(a, b):
1640 return (
1641 a.start_date is None or b.end_date is None
1642 or a.start_date < b.end_date
1643 ) and (
1644 b.start_date is None or a.end_date is None
1645 or b.start_date < a.end_date
1646 )
1649def check_overlaps(node_type_name, codes, nodes):
1650 errors = []
1651 for ia, a in enumerate(nodes):
1652 for b in nodes[ia+1:]:
1653 if overlaps(a, b):
1654 errors.append(
1655 '%s epochs overlap for %s:\n'
1656 ' %s - %s\n %s - %s' % (
1657 node_type_name,
1658 '.'.join(codes),
1659 tts(a.start_date), tts(a.end_date),
1660 tts(b.start_date), tts(b.end_date)))
1662 return errors
1665class Channel(BaseNode):
1666 '''
1667 Equivalent to SEED blockette 52 and parent element for the related the
1668 response blockettes.
1669 '''
1671 location_code = String.T(xmlstyle='attribute')
1672 external_reference_list = List.T(
1673 ExternalReference.T(xmltagname='ExternalReference'))
1674 latitude = Latitude.T(xmltagname='Latitude')
1675 longitude = Longitude.T(xmltagname='Longitude')
1676 elevation = Distance.T(xmltagname='Elevation')
1677 depth = Distance.T(xmltagname='Depth')
1678 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1679 dip = Dip.T(optional=True, xmltagname='Dip')
1680 type_list = List.T(Type.T(xmltagname='Type'))
1681 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1682 sample_rate_ratio = SampleRateRatio.T(optional=True,
1683 xmltagname='SampleRateRatio')
1684 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1685 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1686 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1687 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1688 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1689 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1690 equipment = Equipment.T(optional=True, xmltagname='Equipment')
1691 response = Response.T(optional=True, xmltagname='Response')
1693 @property
1694 def position_values(self):
1695 lat = self.latitude.value
1696 lon = self.longitude.value
1697 elevation = value_or_none(self.elevation)
1698 depth = value_or_none(self.depth)
1699 return lat, lon, elevation, depth
1702class Station(BaseNode):
1703 '''
1704 This type represents a Station epoch. It is common to only have a single
1705 station epoch with the station's creation and termination dates as the
1706 epoch start and end dates.
1707 '''
1709 latitude = Latitude.T(xmltagname='Latitude')
1710 longitude = Longitude.T(xmltagname='Longitude')
1711 elevation = Distance.T(xmltagname='Elevation')
1712 site = Site.T(optional=True, xmltagname='Site')
1713 vault = Unicode.T(optional=True, xmltagname='Vault')
1714 geology = Unicode.T(optional=True, xmltagname='Geology')
1715 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1716 operator_list = List.T(Operator.T(xmltagname='Operator'))
1717 creation_date = DummyAwareOptionalTimestamp.T(
1718 optional=True, xmltagname='CreationDate')
1719 termination_date = DummyAwareOptionalTimestamp.T(
1720 optional=True, xmltagname='TerminationDate')
1721 total_number_channels = Counter.T(
1722 optional=True, xmltagname='TotalNumberChannels')
1723 selected_number_channels = Counter.T(
1724 optional=True, xmltagname='SelectedNumberChannels')
1725 external_reference_list = List.T(
1726 ExternalReference.T(xmltagname='ExternalReference'))
1727 channel_list = List.T(Channel.T(xmltagname='Channel'))
1729 @property
1730 def position_values(self):
1731 lat = self.latitude.value
1732 lon = self.longitude.value
1733 elevation = value_or_none(self.elevation)
1734 return lat, lon, elevation
1737class Network(BaseNode):
1738 '''
1739 This type represents the Network layer, all station metadata is contained
1740 within this element. The official name of the network or other descriptive
1741 information can be included in the Description element. The Network can
1742 contain 0 or more Stations.
1743 '''
1745 total_number_stations = Counter.T(optional=True,
1746 xmltagname='TotalNumberStations')
1747 selected_number_stations = Counter.T(optional=True,
1748 xmltagname='SelectedNumberStations')
1749 station_list = List.T(Station.T(xmltagname='Station'))
1751 @property
1752 def station_code_list(self):
1753 return sorted(set(s.code for s in self.station_list))
1755 @property
1756 def sl_code_list(self):
1757 sls = set()
1758 for station in self.station_list:
1759 for channel in station.channel_list:
1760 sls.add((station.code, channel.location_code))
1762 return sorted(sls)
1764 def summary(self, width=80, indent=4):
1765 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1766 lines = ['%s (%i):' % (self.code, len(sls))]
1767 if sls:
1768 ssls = ['.'.join(x for x in c if x) for c in sls]
1769 w = max(len(x) for x in ssls)
1770 n = (width - indent) / (w+1)
1771 while ssls:
1772 lines.append(
1773 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1775 ssls[:n] = []
1777 return '\n'.join(lines)
1780def value_or_none(x):
1781 if x is not None:
1782 return x.value
1783 else:
1784 return None
1787def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1789 pos = lat, lon, elevation, depth = \
1790 channels[0].position_values
1792 if not all(pos == x.position_values for x in channels):
1793 info = '\n'.join(
1794 ' %s: %s' % (x.code, x.position_values) for
1795 x in channels)
1797 mess = 'encountered inconsistencies in channel ' \
1798 'lat/lon/elevation/depth ' \
1799 'for %s.%s.%s: \n%s' % (nsl + (info,))
1801 if inconsistencies == 'raise':
1802 raise InconsistentChannelLocations(mess)
1804 elif inconsistencies == 'warn':
1805 logger.warning(mess)
1806 logger.warning(' -> using mean values')
1808 apos = num.array([x.position_values for x in channels], dtype=float)
1809 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1810 / num.sum(num.isfinite(apos), axis=0)
1812 groups = {}
1813 for channel in channels:
1814 if channel.code not in groups:
1815 groups[channel.code] = []
1817 groups[channel.code].append(channel)
1819 pchannels = []
1820 for code in sorted(groups.keys()):
1821 data = [
1822 (channel.code, value_or_none(channel.azimuth),
1823 value_or_none(channel.dip))
1824 for channel in groups[code]]
1826 azimuth, dip = util.consistency_merge(
1827 data,
1828 message='channel orientation values differ:',
1829 error=inconsistencies)
1831 pchannels.append(
1832 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1834 return pyrocko.model.Station(
1835 *nsl,
1836 lat=mlat,
1837 lon=mlon,
1838 elevation=mele,
1839 depth=mdep,
1840 channels=pchannels)
1843class FDSNStationXML(Object):
1844 '''
1845 Top-level type for Station XML. Required field are Source (network ID of
1846 the institution sending the message) and one or more Network containers or
1847 one or more Station containers.
1848 '''
1850 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1851 source = String.T(xmltagname='Source')
1852 sender = String.T(optional=True, xmltagname='Sender')
1853 module = String.T(optional=True, xmltagname='Module')
1854 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1855 created = Timestamp.T(optional=True, xmltagname='Created')
1856 network_list = List.T(Network.T(xmltagname='Network'))
1858 xmltagname = 'FDSNStationXML'
1859 guessable_xmlns = [guts_xmlns]
1861 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1862 time=None, timespan=None,
1863 inconsistencies='warn'):
1865 assert inconsistencies in ('raise', 'warn')
1867 if nslcs is not None:
1868 nslcs = set(nslcs)
1870 if nsls is not None:
1871 nsls = set(nsls)
1873 tt = ()
1874 if time is not None:
1875 tt = (time,)
1876 elif timespan is not None:
1877 tt = timespan
1879 pstations = []
1880 for network in self.network_list:
1881 if not network.spans(*tt):
1882 continue
1884 for station in network.station_list:
1885 if not station.spans(*tt):
1886 continue
1888 if station.channel_list:
1889 loc_to_channels = {}
1890 for channel in station.channel_list:
1891 if not channel.spans(*tt):
1892 continue
1894 loc = channel.location_code.strip()
1895 if loc not in loc_to_channels:
1896 loc_to_channels[loc] = []
1898 loc_to_channels[loc].append(channel)
1900 for loc in sorted(loc_to_channels.keys()):
1901 channels = loc_to_channels[loc]
1902 if nslcs is not None:
1903 channels = [channel for channel in channels
1904 if (network.code, station.code, loc,
1905 channel.code) in nslcs]
1907 if not channels:
1908 continue
1910 nsl = network.code, station.code, loc
1911 if nsls is not None and nsl not in nsls:
1912 continue
1914 pstations.append(
1915 pyrocko_station_from_channels(
1916 nsl,
1917 channels,
1918 inconsistencies=inconsistencies))
1919 else:
1920 pstations.append(pyrocko.model.Station(
1921 network.code, station.code, '*',
1922 lat=station.latitude.value,
1923 lon=station.longitude.value,
1924 elevation=value_or_none(station.elevation),
1925 name=station.description or ''))
1927 return pstations
1929 @classmethod
1930 def from_pyrocko_stations(
1931 cls, pyrocko_stations, add_flat_responses_from=None):
1933 '''
1934 Generate :py:class:`FDSNStationXML` from list of
1935 :py:class;`pyrocko.model.Station` instances.
1937 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1938 instances.
1939 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1940 '''
1941 from collections import defaultdict
1942 network_dict = defaultdict(list)
1944 if add_flat_responses_from:
1945 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1946 extra = dict(
1947 response=Response(
1948 instrument_sensitivity=Sensitivity(
1949 value=1.0,
1950 frequency=1.0,
1951 input_units=Units(name=add_flat_responses_from))))
1952 else:
1953 extra = {}
1955 have_offsets = set()
1956 for s in pyrocko_stations:
1958 if s.north_shift != 0.0 or s.east_shift != 0.0:
1959 have_offsets.add(s.nsl())
1961 network, station, location = s.nsl()
1962 channel_list = []
1963 for c in s.channels:
1964 channel_list.append(
1965 Channel(
1966 location_code=location,
1967 code=c.name,
1968 latitude=Latitude(value=s.effective_lat),
1969 longitude=Longitude(value=s.effective_lon),
1970 elevation=Distance(value=s.elevation),
1971 depth=Distance(value=s.depth),
1972 azimuth=Azimuth(value=c.azimuth),
1973 dip=Dip(value=c.dip),
1974 **extra
1975 )
1976 )
1978 network_dict[network].append(
1979 Station(
1980 code=station,
1981 latitude=Latitude(value=s.effective_lat),
1982 longitude=Longitude(value=s.effective_lon),
1983 elevation=Distance(value=s.elevation),
1984 channel_list=channel_list)
1985 )
1987 if have_offsets:
1988 logger.warning(
1989 'StationXML does not support Cartesian offsets in '
1990 'coordinates. Storing effective lat/lon for stations: %s' %
1991 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1993 timestamp = util.to_time_float(time.time())
1994 network_list = []
1995 for k, station_list in network_dict.items():
1997 network_list.append(
1998 Network(
1999 code=k, station_list=station_list,
2000 total_number_stations=len(station_list)))
2002 sxml = FDSNStationXML(
2003 source='from pyrocko stations list', created=timestamp,
2004 network_list=network_list)
2006 sxml.validate()
2007 return sxml
2009 def iter_network_stations(
2010 self, net=None, sta=None, time=None, timespan=None):
2012 tt = ()
2013 if time is not None:
2014 tt = (time,)
2015 elif timespan is not None:
2016 tt = timespan
2018 for network in self.network_list:
2019 if not network.spans(*tt) or (
2020 net is not None and network.code != net):
2021 continue
2023 for station in network.station_list:
2024 if not station.spans(*tt) or (
2025 sta is not None and station.code != sta):
2026 continue
2028 yield (network, station)
2030 def iter_network_station_channels(
2031 self, net=None, sta=None, loc=None, cha=None,
2032 time=None, timespan=None):
2034 if loc is not None:
2035 loc = loc.strip()
2037 tt = ()
2038 if time is not None:
2039 tt = (time,)
2040 elif timespan is not None:
2041 tt = timespan
2043 for network in self.network_list:
2044 if not network.spans(*tt) or (
2045 net is not None and network.code != net):
2046 continue
2048 for station in network.station_list:
2049 if not station.spans(*tt) or (
2050 sta is not None and station.code != sta):
2051 continue
2053 if station.channel_list:
2054 for channel in station.channel_list:
2055 if (not channel.spans(*tt) or
2056 (cha is not None and channel.code != cha) or
2057 (loc is not None and
2058 channel.location_code.strip() != loc)):
2059 continue
2061 yield (network, station, channel)
2063 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2064 time=None, timespan=None):
2066 groups = {}
2067 for network, station, channel in self.iter_network_station_channels(
2068 net, sta, loc, cha, time=time, timespan=timespan):
2070 net = network.code
2071 sta = station.code
2072 cha = channel.code
2073 loc = channel.location_code.strip()
2074 if len(cha) == 3:
2075 bic = cha[:2] # band and intrument code according to SEED
2076 elif len(cha) == 1:
2077 bic = ''
2078 else:
2079 bic = cha
2081 if channel.response and \
2082 channel.response.instrument_sensitivity and \
2083 channel.response.instrument_sensitivity.input_units:
2085 unit = channel.response.instrument_sensitivity\
2086 .input_units.name.upper()
2087 else:
2088 unit = None
2090 bic = (bic, unit)
2092 k = net, sta, loc
2093 if k not in groups:
2094 groups[k] = {}
2096 if bic not in groups[k]:
2097 groups[k][bic] = []
2099 groups[k][bic].append(channel)
2101 for nsl, bic_to_channels in groups.items():
2102 bad_bics = []
2103 for bic, channels in bic_to_channels.items():
2104 sample_rates = []
2105 for channel in channels:
2106 sample_rates.append(channel.sample_rate.value)
2108 if not same(sample_rates):
2109 scs = ','.join(channel.code for channel in channels)
2110 srs = ', '.join('%e' % x for x in sample_rates)
2111 err = 'ignoring channels with inconsistent sampling ' + \
2112 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2114 logger.warning(err)
2115 bad_bics.append(bic)
2117 for bic in bad_bics:
2118 del bic_to_channels[bic]
2120 return groups
2122 def choose_channels(
2123 self,
2124 target_sample_rate=None,
2125 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2126 priority_units=['M/S', 'M/S**2'],
2127 priority_instrument_code=['H', 'L'],
2128 time=None,
2129 timespan=None):
2131 nslcs = {}
2132 for nsl, bic_to_channels in self.get_channel_groups(
2133 time=time, timespan=timespan).items():
2135 useful_bics = []
2136 for bic, channels in bic_to_channels.items():
2137 rate = channels[0].sample_rate.value
2139 if target_sample_rate is not None and \
2140 rate < target_sample_rate*0.99999:
2141 continue
2143 if len(bic[0]) == 2:
2144 if bic[0][0] not in priority_band_code:
2145 continue
2147 if bic[0][1] not in priority_instrument_code:
2148 continue
2150 unit = bic[1]
2152 prio_unit = len(priority_units)
2153 try:
2154 prio_unit = priority_units.index(unit)
2155 except ValueError:
2156 pass
2158 prio_inst = len(priority_instrument_code)
2159 prio_band = len(priority_band_code)
2160 if len(channels[0].code) == 3:
2161 try:
2162 prio_inst = priority_instrument_code.index(
2163 channels[0].code[1])
2164 except ValueError:
2165 pass
2167 try:
2168 prio_band = priority_band_code.index(
2169 channels[0].code[0])
2170 except ValueError:
2171 pass
2173 if target_sample_rate is None:
2174 rate = -rate
2176 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2177 prio_inst, bic))
2179 useful_bics.sort()
2181 for _, _, rate, _, _, bic in useful_bics:
2182 channels = sorted(
2183 bic_to_channels[bic],
2184 key=lambda channel: channel.code)
2186 if channels:
2187 for channel in channels:
2188 nslcs[nsl + (channel.code,)] = channel
2190 break
2192 return nslcs
2194 def get_pyrocko_response(
2195 self, nslc,
2196 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2198 net, sta, loc, cha = nslc
2199 resps = []
2200 for _, _, channel in self.iter_network_station_channels(
2201 net, sta, loc, cha, time=time, timespan=timespan):
2202 resp = channel.response
2203 if resp:
2204 resp.check_sample_rates(channel)
2205 resp.check_units()
2206 resps.append(resp.get_pyrocko_response(
2207 '.'.join(nslc),
2208 fake_input_units=fake_input_units,
2209 stages=stages).expect_one())
2211 if not resps:
2212 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2213 elif len(resps) > 1:
2214 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2216 return resps[0]
2218 @property
2219 def n_code_list(self):
2220 return sorted(set(x.code for x in self.network_list))
2222 @property
2223 def ns_code_list(self):
2224 nss = set()
2225 for network in self.network_list:
2226 for station in network.station_list:
2227 nss.add((network.code, station.code))
2229 return sorted(nss)
2231 @property
2232 def nsl_code_list(self):
2233 nsls = set()
2234 for network in self.network_list:
2235 for station in network.station_list:
2236 for channel in station.channel_list:
2237 nsls.add(
2238 (network.code, station.code, channel.location_code))
2240 return sorted(nsls)
2242 @property
2243 def nslc_code_list(self):
2244 nslcs = set()
2245 for network in self.network_list:
2246 for station in network.station_list:
2247 for channel in station.channel_list:
2248 nslcs.add(
2249 (network.code, station.code, channel.location_code,
2250 channel.code))
2252 return sorted(nslcs)
2254 def summary(self):
2255 lst = [
2256 'number of n codes: %i' % len(self.n_code_list),
2257 'number of ns codes: %i' % len(self.ns_code_list),
2258 'number of nsl codes: %i' % len(self.nsl_code_list),
2259 'number of nslc codes: %i' % len(self.nslc_code_list)
2260 ]
2261 return '\n'.join(lst)
2263 def summary_stages(self):
2264 data = []
2265 for network, station, channel in self.iter_network_station_channels():
2266 nslc = (network.code, station.code, channel.location_code,
2267 channel.code)
2269 stages = []
2270 in_units = '?'
2271 out_units = '?'
2272 if channel.response:
2273 sens = channel.response.instrument_sensitivity
2274 if sens:
2275 in_units = sens.input_units.name.upper()
2276 out_units = sens.output_units.name.upper()
2278 for stage in channel.response.stage_list:
2279 stages.append(stage.summary())
2281 data.append(
2282 (nslc, tts(channel.start_date), tts(channel.end_date),
2283 in_units, out_units, stages))
2285 data.sort()
2287 lst = []
2288 for nslc, stmin, stmax, in_units, out_units, stages in data:
2289 lst.append(' %s: %s - %s, %s -> %s' % (
2290 '.'.join(nslc), stmin, stmax, in_units, out_units))
2291 for stage in stages:
2292 lst.append(' %s' % stage)
2294 return '\n'.join(lst)
2296 def _check_overlaps(self):
2297 by_nslc = {}
2298 for network in self.network_list:
2299 for station in network.station_list:
2300 for channel in station.channel_list:
2301 nslc = (network.code, station.code, channel.location_code,
2302 channel.code)
2303 if nslc not in by_nslc:
2304 by_nslc[nslc] = []
2306 by_nslc[nslc].append(channel)
2308 errors = []
2309 for nslc, channels in by_nslc.items():
2310 errors.extend(check_overlaps('Channel', nslc, channels))
2312 return errors
2314 def check(self):
2315 errors = []
2316 for meth in [self._check_overlaps]:
2317 errors.extend(meth())
2319 if errors:
2320 raise Inconsistencies(
2321 'Inconsistencies found in StationXML:\n '
2322 + '\n '.join(errors))
2325def load_channel_table(stream):
2327 networks = {}
2328 stations = {}
2330 for line in stream:
2331 line = str(line.decode('ascii'))
2332 if line.startswith('#'):
2333 continue
2335 t = line.rstrip().split('|')
2337 if len(t) != 17:
2338 logger.warning('Invalid channel record: %s' % line)
2339 continue
2341 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2342 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2344 try:
2345 scale = float(scale)
2346 except ValueError:
2347 scale = None
2349 try:
2350 scale_freq = float(scale_freq)
2351 except ValueError:
2352 scale_freq = None
2354 try:
2355 depth = float(dep)
2356 except ValueError:
2357 depth = 0.0
2359 try:
2360 azi = float(azi)
2361 dip = float(dip)
2362 except ValueError:
2363 azi = None
2364 dip = None
2366 try:
2367 if net not in networks:
2368 network = Network(code=net)
2369 else:
2370 network = networks[net]
2372 if (net, sta) not in stations:
2373 station = Station(
2374 code=sta, latitude=lat, longitude=lon, elevation=ele)
2376 station.regularize()
2377 else:
2378 station = stations[net, sta]
2380 if scale:
2381 resp = Response(
2382 instrument_sensitivity=Sensitivity(
2383 value=scale,
2384 frequency=scale_freq,
2385 input_units=scale_units))
2386 else:
2387 resp = None
2389 channel = Channel(
2390 code=cha,
2391 location_code=loc.strip(),
2392 latitude=lat,
2393 longitude=lon,
2394 elevation=ele,
2395 depth=depth,
2396 azimuth=azi,
2397 dip=dip,
2398 sensor=Equipment(description=sens),
2399 response=resp,
2400 sample_rate=sample_rate,
2401 start_date=start_date,
2402 end_date=end_date or None)
2404 channel.regularize()
2406 except ValidationError:
2407 raise InvalidRecord(line)
2409 if net not in networks:
2410 networks[net] = network
2412 if (net, sta) not in stations:
2413 stations[net, sta] = station
2414 network.station_list.append(station)
2416 station.channel_list.append(channel)
2418 return FDSNStationXML(
2419 source='created from table input',
2420 created=time.time(),
2421 network_list=sorted(networks.values(), key=lambda x: x.code))
2424def primitive_merge(sxs):
2425 networks = []
2426 for sx in sxs:
2427 networks.extend(sx.network_list)
2429 return FDSNStationXML(
2430 source='merged from different sources',
2431 created=time.time(),
2432 network_list=copy.deepcopy(
2433 sorted(networks, key=lambda x: x.code)))
2436def split_channels(sx):
2437 for nslc in sx.nslc_code_list:
2438 network_list = sx.network_list
2439 network_list_filtered = [
2440 network for network in network_list
2441 if network.code == nslc[0]]
2443 for network in network_list_filtered:
2444 sx.network_list = [network]
2445 station_list = network.station_list
2446 station_list_filtered = [
2447 station for station in station_list
2448 if station.code == nslc[1]]
2450 for station in station_list_filtered:
2451 network.station_list = [station]
2452 channel_list = station.channel_list
2453 station.channel_list = [
2454 channel for channel in channel_list
2455 if (channel.location_code, channel.code) == nslc[2:4]]
2457 yield nslc, copy.deepcopy(sx)
2459 station.channel_list = channel_list
2461 network.station_list = station_list
2463 sx.network_list = network_list
2466if __name__ == '__main__':
2467 from optparse import OptionParser
2469 util.setup_logging('pyrocko.io.stationxml', 'warning')
2471 usage = \
2472 'python -m pyrocko.io.stationxml check|stats|stages ' \
2473 '<filename> [options]'
2475 description = '''Torture StationXML file.'''
2477 parser = OptionParser(
2478 usage=usage,
2479 description=description,
2480 formatter=util.BetterHelpFormatter())
2482 (options, args) = parser.parse_args(sys.argv[1:])
2484 if len(args) != 2:
2485 parser.print_help()
2486 sys.exit(1)
2488 action, path = args
2490 sx = load_xml(filename=path)
2491 if action == 'check':
2492 try:
2493 sx.check()
2494 except Inconsistencies as e:
2495 logger.error(e)
2496 sys.exit(1)
2498 elif action == 'stats':
2499 print(sx.summary())
2501 elif action == 'stages':
2502 print(sx.summary_stages())
2504 else:
2505 parser.print_help()
2506 sys.exit('unknown action: %s' % action)