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 '?'
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 def check_sample_rates(self, channel):
1338 if self.stage_list:
1339 sample_rate = None
1341 for stage in self.stage_list:
1342 if stage.decimation:
1343 input_sample_rate = \
1344 stage.decimation.input_sample_rate.value
1346 if sample_rate is not None and not same_sample_rate(
1347 sample_rate, input_sample_rate):
1349 logger.warning(
1350 'Response stage %i has unexpected input sample '
1351 'rate: %g Hz (expected: %g Hz)' % (
1352 stage.number,
1353 input_sample_rate,
1354 sample_rate))
1356 sample_rate = input_sample_rate / stage.decimation.factor
1358 if sample_rate is not None and channel.sample_rate \
1359 and not same_sample_rate(
1360 sample_rate, channel.sample_rate.value):
1362 logger.warning(
1363 'Channel sample rate (%g Hz) does not match sample rate '
1364 'deducted from response stages information (%g Hz).' % (
1365 channel.sample_rate.value,
1366 sample_rate))
1368 def check_units(self):
1370 if self.instrument_sensitivity \
1371 and self.instrument_sensitivity.input_units:
1373 units = self.instrument_sensitivity.input_units.name.upper()
1375 if self.stage_list:
1376 for stage in self.stage_list:
1377 if units and stage.input_units \
1378 and stage.input_units.name.upper() != units:
1380 logger.warning(
1381 'Input units of stage %i (%s) do not match %s (%s).'
1382 % (
1383 stage.number,
1384 units,
1385 'output units of stage %i'
1386 if stage.number == 0
1387 else 'sensitivity input units',
1388 units))
1390 if stage.output_units:
1391 units = stage.output_units.name.upper()
1392 else:
1393 units = None
1395 sout_units = self.instrument_sensitivity.output_units
1396 if self.instrument_sensitivity and sout_units:
1397 if units is not None and units != sout_units.name.upper():
1398 logger.warning(
1399 'Output units of stage %i (%s) do not match %s (%s).'
1400 % (
1401 stage.number,
1402 units,
1403 'sensitivity output units',
1404 sout_units.name.upper()))
1406 def _sensitivity_checkpoints(self, responses, context):
1407 delivery = Delivery()
1409 if self.instrument_sensitivity:
1410 sval = self.instrument_sensitivity.value
1411 sfreq = self.instrument_sensitivity.frequency
1412 if sval is None:
1413 delivery.log.append((
1414 'warning',
1415 'No sensitivity value given.',
1416 context))
1418 elif sval is None:
1419 delivery.log.append((
1420 'warning',
1421 'Reported sensitivity value is zero.',
1422 context))
1424 elif sfreq is None:
1425 delivery.log.append((
1426 'warning',
1427 'Sensitivity frequency not given.',
1428 context))
1430 else:
1431 trial = response.MultiplyResponse(responses)
1433 delivery.extend(
1434 check_resp(
1435 trial, sval, sfreq, 0.1,
1436 'Instrument sensitivity value inconsistent with '
1437 'sensitivity computed from complete response',
1438 context))
1440 delivery.payload.append(response.FrequencyResponseCheckpoint(
1441 frequency=sfreq, value=sval))
1443 return delivery
1445 def get_squirrel_response(self, context, **kwargs):
1446 from pyrocko.squirrel.model import Response
1448 if self.stage_list:
1449 delivery = Delivery()
1450 for istage, stage in enumerate(self.stage_list):
1451 delivery.extend(stage.get_squirrel_response_stage(context))
1453 checkpoints = []
1454 if not delivery.errors:
1455 all_responses = []
1456 for sq_stage in delivery.payload:
1457 all_responses.extend(sq_stage.elements)
1459 checkpoints.extend(delivery.extend_without_payload(
1460 self._sensitivity_checkpoints(all_responses, context)))
1462 sq_stages = delivery.payload
1463 if sq_stages:
1464 if sq_stages[0].input_quantity is None \
1465 and self.instrument_sensitivity is not None:
1467 sq_stages[0].input_quantity = to_quantity(
1468 self.instrument_sensitivity.input_units,
1469 context, delivery)
1470 sq_stages[-1].output_quantity = to_quantity(
1471 self.instrument_sensitivity.output_units,
1472 context, delivery)
1474 sq_stages = delivery.expect()
1476 return Response(
1477 stages=sq_stages,
1478 log=delivery.log,
1479 checkpoints=checkpoints,
1480 **kwargs)
1482 elif self.instrument_sensitivity:
1483 raise NoResponseInformation(
1484 'Only instrument sensitivity given (won\'t use it). (%s).'
1485 % context)
1486 else:
1487 raise NoResponseInformation(
1488 'Empty instrument response (%s).'
1489 % context)
1491 def get_pyrocko_response(
1492 self, context, fake_input_units=None, stages=(0, 1)):
1494 delivery = Delivery()
1495 if self.stage_list:
1496 for istage, stage in enumerate(self.stage_list):
1497 delivery.extend(stage.get_pyrocko_response(
1498 context, gain_only=not (
1499 stages is None or stages[0] <= istage < stages[1])))
1501 elif self.instrument_sensitivity:
1502 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1504 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1505 checkpoints = delivery.extend_without_payload(delivery_cp)
1506 if delivery.errors:
1507 return delivery
1509 if fake_input_units is not None:
1510 if not self.instrument_sensitivity or \
1511 self.instrument_sensitivity.input_units is None:
1513 delivery.errors.append((
1514 'NoResponseInformation',
1515 'No input units given, so cannot convert to requested '
1516 'units: %s' % fake_input_units.upper(),
1517 context))
1519 return delivery
1521 input_units = self.instrument_sensitivity.input_units.name.upper()
1523 conresp = None
1524 try:
1525 conresp = conversion[
1526 fake_input_units.upper(), input_units]
1528 except KeyError:
1529 delivery.errors.append((
1530 'NoResponseInformation',
1531 'Cannot convert between units: %s, %s'
1532 % (fake_input_units, input_units),
1533 context))
1535 if conresp is not None:
1536 delivery.payload.append(conresp)
1537 for checkpoint in checkpoints:
1538 checkpoint.value *= num.abs(evaluate1(
1539 conresp, checkpoint.frequency))
1541 delivery.payload = [
1542 response.MultiplyResponse(
1543 delivery.payload,
1544 checkpoints=checkpoints)]
1546 return delivery
1548 @classmethod
1549 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1550 normalization_frequency=1.0):
1551 '''
1552 Convert Pyrocko pole-zero response to StationXML response.
1554 :param presponse: Pyrocko pole-zero response
1555 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1556 :param input_unit: Input unit to be reported in the StationXML
1557 response.
1558 :type input_unit: str
1559 :param output_unit: Output unit to be reported in the StationXML
1560 response.
1561 :type output_unit: str
1562 :param normalization_frequency: Frequency where the normalization
1563 factor for the StationXML response should be computed.
1564 :type normalization_frequency: float
1565 '''
1567 norm_factor = 1.0/float(abs(
1568 evaluate1(presponse, normalization_frequency)
1569 / presponse.constant))
1571 pzs = PolesZeros(
1572 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1573 normalization_factor=norm_factor,
1574 normalization_frequency=Frequency(normalization_frequency),
1575 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1576 imaginary=FloatNoUnit(z.imag))
1577 for z in presponse.zeros],
1578 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1579 imaginary=FloatNoUnit(z.imag))
1580 for z in presponse.poles])
1582 pzs.validate()
1584 stage = ResponseStage(
1585 number=1,
1586 poles_zeros_list=[pzs],
1587 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1589 resp = Response(
1590 instrument_sensitivity=Sensitivity(
1591 value=stage.stage_gain.value,
1592 frequency=normalization_frequency,
1593 input_units=Units(input_unit),
1594 output_units=Units(output_unit)),
1596 stage_list=[stage])
1598 return resp
1601class BaseNode(Object):
1602 '''
1603 A base node type for derivation from: Network, Station and Channel types.
1604 '''
1606 code = String.T(xmlstyle='attribute')
1607 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1608 xmlstyle='attribute')
1609 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1610 xmlstyle='attribute')
1611 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1612 alternate_code = String.T(optional=True, xmlstyle='attribute')
1613 historical_code = String.T(optional=True, xmlstyle='attribute')
1614 description = Unicode.T(optional=True, xmltagname='Description')
1615 comment_list = List.T(Comment.T(xmltagname='Comment'))
1617 def spans(self, *args):
1618 if len(args) == 0:
1619 return True
1620 elif len(args) == 1:
1621 return ((self.start_date is None or
1622 self.start_date <= args[0]) and
1623 (self.end_date is None or
1624 args[0] <= self.end_date))
1626 elif len(args) == 2:
1627 return ((self.start_date is None or
1628 args[1] >= self.start_date) and
1629 (self.end_date is None or
1630 self.end_date >= args[0]))
1633def overlaps(a, b):
1634 return (
1635 a.start_date is None or b.end_date is None
1636 or a.start_date < b.end_date
1637 ) and (
1638 b.start_date is None or a.end_date is None
1639 or b.start_date < a.end_date
1640 )
1643def check_overlaps(node_type_name, codes, nodes):
1644 errors = []
1645 for ia, a in enumerate(nodes):
1646 for b in nodes[ia+1:]:
1647 if overlaps(a, b):
1648 errors.append(
1649 '%s epochs overlap for %s:\n'
1650 ' %s - %s\n %s - %s' % (
1651 node_type_name,
1652 '.'.join(codes),
1653 tts(a.start_date), tts(a.end_date),
1654 tts(b.start_date), tts(b.end_date)))
1656 return errors
1659class Channel(BaseNode):
1660 '''
1661 Equivalent to SEED blockette 52 and parent element for the related the
1662 response blockettes.
1663 '''
1665 location_code = String.T(xmlstyle='attribute')
1666 external_reference_list = List.T(
1667 ExternalReference.T(xmltagname='ExternalReference'))
1668 latitude = Latitude.T(xmltagname='Latitude')
1669 longitude = Longitude.T(xmltagname='Longitude')
1670 elevation = Distance.T(xmltagname='Elevation')
1671 depth = Distance.T(xmltagname='Depth')
1672 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1673 dip = Dip.T(optional=True, xmltagname='Dip')
1674 type_list = List.T(Type.T(xmltagname='Type'))
1675 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1676 sample_rate_ratio = SampleRateRatio.T(optional=True,
1677 xmltagname='SampleRateRatio')
1678 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1679 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1680 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1681 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1682 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1683 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1684 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1685 response = Response.T(optional=True, xmltagname='Response')
1687 @property
1688 def position_values(self):
1689 lat = self.latitude.value
1690 lon = self.longitude.value
1691 elevation = value_or_none(self.elevation)
1692 depth = value_or_none(self.depth)
1693 return lat, lon, elevation, depth
1696class Station(BaseNode):
1697 '''
1698 This type represents a Station epoch. It is common to only have a single
1699 station epoch with the station's creation and termination dates as the
1700 epoch start and end dates.
1701 '''
1703 latitude = Latitude.T(xmltagname='Latitude')
1704 longitude = Longitude.T(xmltagname='Longitude')
1705 elevation = Distance.T(xmltagname='Elevation')
1706 site = Site.T(optional=True, xmltagname='Site')
1707 vault = Unicode.T(optional=True, xmltagname='Vault')
1708 geology = Unicode.T(optional=True, xmltagname='Geology')
1709 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1710 operator_list = List.T(Operator.T(xmltagname='Operator'))
1711 creation_date = DummyAwareOptionalTimestamp.T(
1712 optional=True, xmltagname='CreationDate')
1713 termination_date = DummyAwareOptionalTimestamp.T(
1714 optional=True, xmltagname='TerminationDate')
1715 total_number_channels = Counter.T(
1716 optional=True, xmltagname='TotalNumberChannels')
1717 selected_number_channels = Counter.T(
1718 optional=True, xmltagname='SelectedNumberChannels')
1719 external_reference_list = List.T(
1720 ExternalReference.T(xmltagname='ExternalReference'))
1721 channel_list = List.T(Channel.T(xmltagname='Channel'))
1723 @property
1724 def position_values(self):
1725 lat = self.latitude.value
1726 lon = self.longitude.value
1727 elevation = value_or_none(self.elevation)
1728 return lat, lon, elevation
1731class Network(BaseNode):
1732 '''
1733 This type represents the Network layer, all station metadata is contained
1734 within this element. The official name of the network or other descriptive
1735 information can be included in the Description element. The Network can
1736 contain 0 or more Stations.
1737 '''
1739 total_number_stations = Counter.T(optional=True,
1740 xmltagname='TotalNumberStations')
1741 selected_number_stations = Counter.T(optional=True,
1742 xmltagname='SelectedNumberStations')
1743 station_list = List.T(Station.T(xmltagname='Station'))
1745 @property
1746 def station_code_list(self):
1747 return sorted(set(s.code for s in self.station_list))
1749 @property
1750 def sl_code_list(self):
1751 sls = set()
1752 for station in self.station_list:
1753 for channel in station.channel_list:
1754 sls.add((station.code, channel.location_code))
1756 return sorted(sls)
1758 def summary(self, width=80, indent=4):
1759 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1760 lines = ['%s (%i):' % (self.code, len(sls))]
1761 if sls:
1762 ssls = ['.'.join(x for x in c if x) for c in sls]
1763 w = max(len(x) for x in ssls)
1764 n = (width - indent) / (w+1)
1765 while ssls:
1766 lines.append(
1767 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1769 ssls[:n] = []
1771 return '\n'.join(lines)
1774def value_or_none(x):
1775 if x is not None:
1776 return x.value
1777 else:
1778 return None
1781def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1783 pos = lat, lon, elevation, depth = \
1784 channels[0].position_values
1786 if not all(pos == x.position_values for x in channels):
1787 info = '\n'.join(
1788 ' %s: %s' % (x.code, x.position_values) for
1789 x in channels)
1791 mess = 'encountered inconsistencies in channel ' \
1792 'lat/lon/elevation/depth ' \
1793 'for %s.%s.%s: \n%s' % (nsl + (info,))
1795 if inconsistencies == 'raise':
1796 raise InconsistentChannelLocations(mess)
1798 elif inconsistencies == 'warn':
1799 logger.warning(mess)
1800 logger.warning(' -> using mean values')
1802 apos = num.array([x.position_values for x in channels], dtype=float)
1803 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1804 / num.sum(num.isfinite(apos), axis=0)
1806 groups = {}
1807 for channel in channels:
1808 if channel.code not in groups:
1809 groups[channel.code] = []
1811 groups[channel.code].append(channel)
1813 pchannels = []
1814 for code in sorted(groups.keys()):
1815 data = [
1816 (channel.code, value_or_none(channel.azimuth),
1817 value_or_none(channel.dip))
1818 for channel in groups[code]]
1820 azimuth, dip = util.consistency_merge(
1821 data,
1822 message='channel orientation values differ:',
1823 error=inconsistencies)
1825 pchannels.append(
1826 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1828 return pyrocko.model.Station(
1829 *nsl,
1830 lat=mlat,
1831 lon=mlon,
1832 elevation=mele,
1833 depth=mdep,
1834 channels=pchannels)
1837class FDSNStationXML(Object):
1838 '''
1839 Top-level type for Station XML. Required field are Source (network ID of
1840 the institution sending the message) and one or more Network containers or
1841 one or more Station containers.
1842 '''
1844 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1845 source = String.T(xmltagname='Source')
1846 sender = String.T(optional=True, xmltagname='Sender')
1847 module = String.T(optional=True, xmltagname='Module')
1848 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1849 created = Timestamp.T(optional=True, xmltagname='Created')
1850 network_list = List.T(Network.T(xmltagname='Network'))
1852 xmltagname = 'FDSNStationXML'
1853 guessable_xmlns = [guts_xmlns]
1855 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1856 time=None, timespan=None,
1857 inconsistencies='warn'):
1859 assert inconsistencies in ('raise', 'warn')
1861 if nslcs is not None:
1862 nslcs = set(nslcs)
1864 if nsls is not None:
1865 nsls = set(nsls)
1867 tt = ()
1868 if time is not None:
1869 tt = (time,)
1870 elif timespan is not None:
1871 tt = timespan
1873 pstations = []
1874 for network in self.network_list:
1875 if not network.spans(*tt):
1876 continue
1878 for station in network.station_list:
1879 if not station.spans(*tt):
1880 continue
1882 if station.channel_list:
1883 loc_to_channels = {}
1884 for channel in station.channel_list:
1885 if not channel.spans(*tt):
1886 continue
1888 loc = channel.location_code.strip()
1889 if loc not in loc_to_channels:
1890 loc_to_channels[loc] = []
1892 loc_to_channels[loc].append(channel)
1894 for loc in sorted(loc_to_channels.keys()):
1895 channels = loc_to_channels[loc]
1896 if nslcs is not None:
1897 channels = [channel for channel in channels
1898 if (network.code, station.code, loc,
1899 channel.code) in nslcs]
1901 if not channels:
1902 continue
1904 nsl = network.code, station.code, loc
1905 if nsls is not None and nsl not in nsls:
1906 continue
1908 pstations.append(
1909 pyrocko_station_from_channels(
1910 nsl,
1911 channels,
1912 inconsistencies=inconsistencies))
1913 else:
1914 pstations.append(pyrocko.model.Station(
1915 network.code, station.code, '*',
1916 lat=station.latitude.value,
1917 lon=station.longitude.value,
1918 elevation=value_or_none(station.elevation),
1919 name=station.description or ''))
1921 return pstations
1923 @classmethod
1924 def from_pyrocko_stations(
1925 cls, pyrocko_stations, add_flat_responses_from=None):
1927 '''
1928 Generate :py:class:`FDSNStationXML` from list of
1929 :py:class;`pyrocko.model.Station` instances.
1931 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1932 instances.
1933 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1934 '''
1935 from collections import defaultdict
1936 network_dict = defaultdict(list)
1938 if add_flat_responses_from:
1939 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1940 extra = dict(
1941 response=Response(
1942 instrument_sensitivity=Sensitivity(
1943 value=1.0,
1944 frequency=1.0,
1945 input_units=Units(name=add_flat_responses_from))))
1946 else:
1947 extra = {}
1949 have_offsets = set()
1950 for s in pyrocko_stations:
1952 if s.north_shift != 0.0 or s.east_shift != 0.0:
1953 have_offsets.add(s.nsl())
1955 network, station, location = s.nsl()
1956 channel_list = []
1957 for c in s.channels:
1958 channel_list.append(
1959 Channel(
1960 location_code=location,
1961 code=c.name,
1962 latitude=Latitude(value=s.effective_lat),
1963 longitude=Longitude(value=s.effective_lon),
1964 elevation=Distance(value=s.elevation),
1965 depth=Distance(value=s.depth),
1966 azimuth=Azimuth(value=c.azimuth),
1967 dip=Dip(value=c.dip),
1968 **extra
1969 )
1970 )
1972 network_dict[network].append(
1973 Station(
1974 code=station,
1975 latitude=Latitude(value=s.effective_lat),
1976 longitude=Longitude(value=s.effective_lon),
1977 elevation=Distance(value=s.elevation),
1978 channel_list=channel_list)
1979 )
1981 if have_offsets:
1982 logger.warning(
1983 'StationXML does not support Cartesian offsets in '
1984 'coordinates. Storing effective lat/lon for stations: %s' %
1985 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1987 timestamp = util.to_time_float(time.time())
1988 network_list = []
1989 for k, station_list in network_dict.items():
1991 network_list.append(
1992 Network(
1993 code=k, station_list=station_list,
1994 total_number_stations=len(station_list)))
1996 sxml = FDSNStationXML(
1997 source='from pyrocko stations list', created=timestamp,
1998 network_list=network_list)
2000 sxml.validate()
2001 return sxml
2003 def iter_network_stations(
2004 self, net=None, sta=None, time=None, timespan=None):
2006 tt = ()
2007 if time is not None:
2008 tt = (time,)
2009 elif timespan is not None:
2010 tt = timespan
2012 for network in self.network_list:
2013 if not network.spans(*tt) or (
2014 net is not None and network.code != net):
2015 continue
2017 for station in network.station_list:
2018 if not station.spans(*tt) or (
2019 sta is not None and station.code != sta):
2020 continue
2022 yield (network, station)
2024 def iter_network_station_channels(
2025 self, net=None, sta=None, loc=None, cha=None,
2026 time=None, timespan=None):
2028 if loc is not None:
2029 loc = loc.strip()
2031 tt = ()
2032 if time is not None:
2033 tt = (time,)
2034 elif timespan is not None:
2035 tt = timespan
2037 for network in self.network_list:
2038 if not network.spans(*tt) or (
2039 net is not None and network.code != net):
2040 continue
2042 for station in network.station_list:
2043 if not station.spans(*tt) or (
2044 sta is not None and station.code != sta):
2045 continue
2047 if station.channel_list:
2048 for channel in station.channel_list:
2049 if (not channel.spans(*tt) or
2050 (cha is not None and channel.code != cha) or
2051 (loc is not None and
2052 channel.location_code.strip() != loc)):
2053 continue
2055 yield (network, station, channel)
2057 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2058 time=None, timespan=None):
2060 groups = {}
2061 for network, station, channel in self.iter_network_station_channels(
2062 net, sta, loc, cha, time=time, timespan=timespan):
2064 net = network.code
2065 sta = station.code
2066 cha = channel.code
2067 loc = channel.location_code.strip()
2068 if len(cha) == 3:
2069 bic = cha[:2] # band and intrument code according to SEED
2070 elif len(cha) == 1:
2071 bic = ''
2072 else:
2073 bic = cha
2075 if channel.response and \
2076 channel.response.instrument_sensitivity and \
2077 channel.response.instrument_sensitivity.input_units:
2079 unit = channel.response.instrument_sensitivity\
2080 .input_units.name.upper()
2081 else:
2082 unit = None
2084 bic = (bic, unit)
2086 k = net, sta, loc
2087 if k not in groups:
2088 groups[k] = {}
2090 if bic not in groups[k]:
2091 groups[k][bic] = []
2093 groups[k][bic].append(channel)
2095 for nsl, bic_to_channels in groups.items():
2096 bad_bics = []
2097 for bic, channels in bic_to_channels.items():
2098 sample_rates = []
2099 for channel in channels:
2100 sample_rates.append(channel.sample_rate.value)
2102 if not same(sample_rates):
2103 scs = ','.join(channel.code for channel in channels)
2104 srs = ', '.join('%e' % x for x in sample_rates)
2105 err = 'ignoring channels with inconsistent sampling ' + \
2106 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2108 logger.warning(err)
2109 bad_bics.append(bic)
2111 for bic in bad_bics:
2112 del bic_to_channels[bic]
2114 return groups
2116 def choose_channels(
2117 self,
2118 target_sample_rate=None,
2119 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2120 priority_units=['M/S', 'M/S**2'],
2121 priority_instrument_code=['H', 'L'],
2122 time=None,
2123 timespan=None):
2125 nslcs = {}
2126 for nsl, bic_to_channels in self.get_channel_groups(
2127 time=time, timespan=timespan).items():
2129 useful_bics = []
2130 for bic, channels in bic_to_channels.items():
2131 rate = channels[0].sample_rate.value
2133 if target_sample_rate is not None and \
2134 rate < target_sample_rate*0.99999:
2135 continue
2137 if len(bic[0]) == 2:
2138 if bic[0][0] not in priority_band_code:
2139 continue
2141 if bic[0][1] not in priority_instrument_code:
2142 continue
2144 unit = bic[1]
2146 prio_unit = len(priority_units)
2147 try:
2148 prio_unit = priority_units.index(unit)
2149 except ValueError:
2150 pass
2152 prio_inst = len(priority_instrument_code)
2153 prio_band = len(priority_band_code)
2154 if len(channels[0].code) == 3:
2155 try:
2156 prio_inst = priority_instrument_code.index(
2157 channels[0].code[1])
2158 except ValueError:
2159 pass
2161 try:
2162 prio_band = priority_band_code.index(
2163 channels[0].code[0])
2164 except ValueError:
2165 pass
2167 if target_sample_rate is None:
2168 rate = -rate
2170 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2171 prio_inst, bic))
2173 useful_bics.sort()
2175 for _, _, rate, _, _, bic in useful_bics:
2176 channels = sorted(
2177 bic_to_channels[bic],
2178 key=lambda channel: channel.code)
2180 if channels:
2181 for channel in channels:
2182 nslcs[nsl + (channel.code,)] = channel
2184 break
2186 return nslcs
2188 def get_pyrocko_response(
2189 self, nslc,
2190 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2192 net, sta, loc, cha = nslc
2193 resps = []
2194 for _, _, channel in self.iter_network_station_channels(
2195 net, sta, loc, cha, time=time, timespan=timespan):
2196 resp = channel.response
2197 if resp:
2198 resp.check_sample_rates(channel)
2199 resp.check_units()
2200 resps.append(resp.get_pyrocko_response(
2201 '.'.join(nslc),
2202 fake_input_units=fake_input_units,
2203 stages=stages).expect_one())
2205 if not resps:
2206 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2207 elif len(resps) > 1:
2208 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2210 return resps[0]
2212 @property
2213 def n_code_list(self):
2214 return sorted(set(x.code for x in self.network_list))
2216 @property
2217 def ns_code_list(self):
2218 nss = set()
2219 for network in self.network_list:
2220 for station in network.station_list:
2221 nss.add((network.code, station.code))
2223 return sorted(nss)
2225 @property
2226 def nsl_code_list(self):
2227 nsls = set()
2228 for network in self.network_list:
2229 for station in network.station_list:
2230 for channel in station.channel_list:
2231 nsls.add(
2232 (network.code, station.code, channel.location_code))
2234 return sorted(nsls)
2236 @property
2237 def nslc_code_list(self):
2238 nslcs = set()
2239 for network in self.network_list:
2240 for station in network.station_list:
2241 for channel in station.channel_list:
2242 nslcs.add(
2243 (network.code, station.code, channel.location_code,
2244 channel.code))
2246 return sorted(nslcs)
2248 def summary(self):
2249 lst = [
2250 'number of n codes: %i' % len(self.n_code_list),
2251 'number of ns codes: %i' % len(self.ns_code_list),
2252 'number of nsl codes: %i' % len(self.nsl_code_list),
2253 'number of nslc codes: %i' % len(self.nslc_code_list)
2254 ]
2255 return '\n'.join(lst)
2257 def summary_stages(self):
2258 data = []
2259 for network, station, channel in self.iter_network_station_channels():
2260 nslc = (network.code, station.code, channel.location_code,
2261 channel.code)
2263 stages = []
2264 in_units = '?'
2265 out_units = '?'
2266 if channel.response:
2267 sens = channel.response.instrument_sensitivity
2268 if sens:
2269 in_units = sens.input_units.name.upper()
2270 out_units = sens.output_units.name.upper()
2272 for stage in channel.response.stage_list:
2273 stages.append(stage.summary())
2275 data.append(
2276 (nslc, tts(channel.start_date), tts(channel.end_date),
2277 in_units, out_units, stages))
2279 data.sort()
2281 lst = []
2282 for nslc, stmin, stmax, in_units, out_units, stages in data:
2283 lst.append(' %s: %s - %s, %s -> %s' % (
2284 '.'.join(nslc), stmin, stmax, in_units, out_units))
2285 for stage in stages:
2286 lst.append(' %s' % stage)
2288 return '\n'.join(lst)
2290 def _check_overlaps(self):
2291 by_nslc = {}
2292 for network in self.network_list:
2293 for station in network.station_list:
2294 for channel in station.channel_list:
2295 nslc = (network.code, station.code, channel.location_code,
2296 channel.code)
2297 if nslc not in by_nslc:
2298 by_nslc[nslc] = []
2300 by_nslc[nslc].append(channel)
2302 errors = []
2303 for nslc, channels in by_nslc.items():
2304 errors.extend(check_overlaps('Channel', nslc, channels))
2306 return errors
2308 def check(self):
2309 errors = []
2310 for meth in [self._check_overlaps]:
2311 errors.extend(meth())
2313 if errors:
2314 raise Inconsistencies(
2315 'Inconsistencies found in StationXML:\n '
2316 + '\n '.join(errors))
2319def load_channel_table(stream):
2321 networks = {}
2322 stations = {}
2324 for line in stream:
2325 line = str(line.decode('ascii'))
2326 if line.startswith('#'):
2327 continue
2329 t = line.rstrip().split('|')
2331 if len(t) != 17:
2332 logger.warning('Invalid channel record: %s' % line)
2333 continue
2335 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2336 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2338 try:
2339 scale = float(scale)
2340 except ValueError:
2341 scale = None
2343 try:
2344 scale_freq = float(scale_freq)
2345 except ValueError:
2346 scale_freq = None
2348 try:
2349 depth = float(dep)
2350 except ValueError:
2351 depth = 0.0
2353 try:
2354 azi = float(azi)
2355 dip = float(dip)
2356 except ValueError:
2357 azi = None
2358 dip = None
2360 try:
2361 if net not in networks:
2362 network = Network(code=net)
2363 else:
2364 network = networks[net]
2366 if (net, sta) not in stations:
2367 station = Station(
2368 code=sta, latitude=lat, longitude=lon, elevation=ele)
2370 station.regularize()
2371 else:
2372 station = stations[net, sta]
2374 if scale:
2375 resp = Response(
2376 instrument_sensitivity=Sensitivity(
2377 value=scale,
2378 frequency=scale_freq,
2379 input_units=scale_units))
2380 else:
2381 resp = None
2383 channel = Channel(
2384 code=cha,
2385 location_code=loc.strip(),
2386 latitude=lat,
2387 longitude=lon,
2388 elevation=ele,
2389 depth=depth,
2390 azimuth=azi,
2391 dip=dip,
2392 sensor=Equipment(description=sens),
2393 response=resp,
2394 sample_rate=sample_rate,
2395 start_date=start_date,
2396 end_date=end_date or None)
2398 channel.regularize()
2400 except ValidationError:
2401 raise InvalidRecord(line)
2403 if net not in networks:
2404 networks[net] = network
2406 if (net, sta) not in stations:
2407 stations[net, sta] = station
2408 network.station_list.append(station)
2410 station.channel_list.append(channel)
2412 return FDSNStationXML(
2413 source='created from table input',
2414 created=time.time(),
2415 network_list=sorted(networks.values(), key=lambda x: x.code))
2418def primitive_merge(sxs):
2419 networks = []
2420 for sx in sxs:
2421 networks.extend(sx.network_list)
2423 return FDSNStationXML(
2424 source='merged from different sources',
2425 created=time.time(),
2426 network_list=copy.deepcopy(
2427 sorted(networks, key=lambda x: x.code)))
2430def split_channels(sx):
2431 for nslc in sx.nslc_code_list:
2432 network_list = sx.network_list
2433 network_list_filtered = [
2434 network for network in network_list
2435 if network.code == nslc[0]]
2437 for network in network_list_filtered:
2438 sx.network_list = [network]
2439 station_list = network.station_list
2440 station_list_filtered = [
2441 station for station in station_list
2442 if station.code == nslc[1]]
2444 for station in station_list_filtered:
2445 network.station_list = [station]
2446 channel_list = station.channel_list
2447 station.channel_list = [
2448 channel for channel in channel_list
2449 if (channel.location_code, channel.code) == nslc[2:4]]
2451 yield nslc, copy.deepcopy(sx)
2453 station.channel_list = channel_list
2455 network.station_list = station_list
2457 sx.network_list = network_list
2460if __name__ == '__main__':
2461 from optparse import OptionParser
2463 util.setup_logging('pyrocko.io.stationxml', 'warning')
2465 usage = \
2466 'python -m pyrocko.io.stationxml check|stats|stages ' \
2467 '<filename> [options]'
2469 description = '''Torture StationXML file.'''
2471 parser = OptionParser(
2472 usage=usage,
2473 description=description,
2474 formatter=util.BetterHelpFormatter())
2476 (options, args) = parser.parse_args(sys.argv[1:])
2478 if len(args) != 2:
2479 parser.print_help()
2480 sys.exit(1)
2482 action, path = args
2484 sx = load_xml(filename=path)
2485 if action == 'check':
2486 try:
2487 sx.check()
2488 except Inconsistencies as e:
2489 logger.error(e)
2490 sys.exit(1)
2492 elif action == 'stats':
2493 print(sx.summary())
2495 elif action == 'stages':
2496 print(sx.summary_stages())
2498 else:
2499 parser.print_help()
2500 sys.exit('unknown action: %s' % action)