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