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 if val == 'None':
353 return None # StationXML contained "None"
354 else:
355 year = int(val[:4])
357 if sys.maxsize > 2**32: # if we're on 64bit
358 if year > this_year + 100:
359 return None # StationXML contained a dummy date
361 if year < 1903: # for macOS, 1900-01-01 dummy dates
362 return None
364 return None # time-stamp format incomplete
365 else: # 32bit end of time is in 2038
366 if this_year < 2037 and year > 2037 or year < 1903:
367 return None # StationXML contained a dummy date
369 raise
371 elif isinstance(val, (int, float)):
372 val = time_float(val)
374 else:
375 raise ValidationError(
376 '%s: cannot convert "%s" to type %s' % (
377 self.xname(), val, time_float))
379 return val
381 def to_save(self, val):
382 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
383 .rstrip('0').rstrip('.')
385 def to_save_xml(self, val):
386 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
387 .rstrip('0').rstrip('.') + 'Z'
390class Nominal(StringChoice):
391 choices = [
392 'NOMINAL',
393 'CALCULATED']
396class Email(UnicodePattern):
397 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
400class RestrictedStatus(StringChoice):
401 choices = [
402 'open',
403 'closed',
404 'partial',
405 'public']
408class Type(StringChoice):
409 choices = [
410 'TRIGGERED',
411 'CONTINUOUS',
412 'HEALTH',
413 'GEOPHYSICAL',
414 'WEATHER',
415 'FLAG',
416 'SYNTHESIZED',
417 'INPUT',
418 'EXPERIMENTAL',
419 'MAINTENANCE',
420 'BEAM']
422 class __T(StringChoice.T):
423 def validate_extra(self, val):
424 if val not in self.choices:
425 logger.warning(
426 'channel type: "%s" is not a valid choice out of %s' %
427 (val, repr(self.choices)))
430class PzTransferFunction(StringChoice):
431 choices = [
432 'LAPLACE (RADIANS/SECOND)',
433 'LAPLACE (HERTZ)',
434 'DIGITAL (Z-TRANSFORM)']
437class Symmetry(StringChoice):
438 choices = [
439 'NONE',
440 'EVEN',
441 'ODD']
444class CfTransferFunction(StringChoice):
446 class __T(StringChoice.T):
447 def validate(self, val, regularize=False, depth=-1):
448 if regularize:
449 try:
450 val = str(val)
451 except ValueError:
452 raise ValidationError(
453 '%s: cannot convert to string %s' % (self.xname,
454 repr(val)))
456 val = self.dummy_cls.replacements.get(val, val)
458 self.validate_extra(val)
459 return val
461 choices = [
462 'ANALOG (RADIANS/SECOND)',
463 'ANALOG (HERTZ)',
464 'DIGITAL']
466 replacements = {
467 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
468 'ANALOG (HZ)': 'ANALOG (HERTZ)',
469 }
472class Approximation(StringChoice):
473 choices = [
474 'MACLAURIN']
477class PhoneNumber(StringPattern):
478 pattern = '[0-9]+-[0-9]+'
481class Site(Object):
482 '''
483 Description of a site location using name and optional geopolitical
484 boundaries (country, city, etc.).
485 '''
487 name = Unicode.T(default='', xmltagname='Name')
488 description = Unicode.T(optional=True, xmltagname='Description')
489 town = Unicode.T(optional=True, xmltagname='Town')
490 county = Unicode.T(optional=True, xmltagname='County')
491 region = Unicode.T(optional=True, xmltagname='Region')
492 country = Unicode.T(optional=True, xmltagname='Country')
495class ExternalReference(Object):
496 '''
497 This type contains a URI and description for external data that users may
498 want to reference in StationXML.
499 '''
501 uri = String.T(xmltagname='URI')
502 description = Unicode.T(xmltagname='Description')
505class Units(Object):
506 '''
507 A type to document units. Corresponds to SEED blockette 34.
508 '''
510 def __init__(self, name=None, **kwargs):
511 Object.__init__(self, name=name, **kwargs)
513 name = String.T(xmltagname='Name')
514 description = Unicode.T(optional=True, xmltagname='Description')
517class Counter(Int):
519 class __T(TBase):
521 def regularize_extra(self, val):
522 if isinstance(val, str) and len(val) == 0:
523 return None # empty Total Number of Stations empty
524 else:
525 return int(val)
528class SampleRateRatio(Object):
529 '''
530 Sample rate expressed as number of samples in a number of seconds.
531 '''
533 number_samples = Int.T(xmltagname='NumberSamples')
534 number_seconds = Int.T(xmltagname='NumberSeconds')
537class Gain(Object):
538 '''
539 Complex type for sensitivity and frequency ranges. This complex type can be
540 used to represent both overall sensitivities and individual stage gains.
541 The FrequencyRangeGroup is an optional construct that defines a pass band
542 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
543 is valid within the number of decibels specified in FrequencyDBVariation.
544 '''
546 def __init__(self, value=None, **kwargs):
547 Object.__init__(self, value=value, **kwargs)
549 value = Float.T(optional=True, xmltagname='Value')
550 frequency = Float.T(optional=True, xmltagname='Frequency')
552 def summary(self):
553 return 'gain(%g @ %g)' % (self.value, self.frequency)
556class NumeratorCoefficient(Object):
557 i = Int.T(optional=True, xmlstyle='attribute')
558 value = Float.T(xmlstyle='content')
561class FloatNoUnit(Object):
562 def __init__(self, value=None, **kwargs):
563 Object.__init__(self, value=value, **kwargs)
565 plus_error = Float.T(optional=True, xmlstyle='attribute')
566 minus_error = Float.T(optional=True, xmlstyle='attribute')
567 value = Float.T(xmlstyle='content')
570class FloatWithUnit(FloatNoUnit):
571 unit = String.T(optional=True, xmlstyle='attribute')
574class Equipment(Object):
575 resource_id = String.T(optional=True, xmlstyle='attribute')
576 type = String.T(optional=True, xmltagname='Type')
577 description = Unicode.T(optional=True, xmltagname='Description')
578 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
579 vendor = Unicode.T(optional=True, xmltagname='Vendor')
580 model = Unicode.T(optional=True, xmltagname='Model')
581 serial_number = String.T(optional=True, xmltagname='SerialNumber')
582 installation_date = DummyAwareOptionalTimestamp.T(
583 optional=True,
584 xmltagname='InstallationDate')
585 removal_date = DummyAwareOptionalTimestamp.T(
586 optional=True,
587 xmltagname='RemovalDate')
588 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
591class PhoneNumber(Object):
592 description = Unicode.T(optional=True, xmlstyle='attribute')
593 country_code = Int.T(optional=True, xmltagname='CountryCode')
594 area_code = Int.T(xmltagname='AreaCode')
595 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
598class BaseFilter(Object):
599 '''
600 The BaseFilter is derived by all filters.
601 '''
603 resource_id = String.T(optional=True, xmlstyle='attribute')
604 name = String.T(optional=True, xmlstyle='attribute')
605 description = Unicode.T(optional=True, xmltagname='Description')
606 input_units = Units.T(optional=True, xmltagname='InputUnits')
607 output_units = Units.T(optional=True, xmltagname='OutputUnits')
610class Sensitivity(Gain):
611 '''
612 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
613 construct that defines a pass band in Hertz (FrequencyStart and
614 FrequencyEnd) in which the SensitivityValue is valid within the number of
615 decibels specified in FrequencyDBVariation.
616 '''
618 input_units = Units.T(optional=True, xmltagname='InputUnits')
619 output_units = Units.T(optional=True, xmltagname='OutputUnits')
620 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
621 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
622 frequency_db_variation = Float.T(optional=True,
623 xmltagname='FrequencyDBVariation')
625 def get_pyrocko_response(self):
626 return Delivery(
627 [response.PoleZeroResponse(constant=self.value)])
630class Coefficient(FloatNoUnit):
631 number = Counter.T(optional=True, xmlstyle='attribute')
634class PoleZero(Object):
635 '''
636 Complex numbers used as poles or zeros in channel response.
637 '''
639 number = Int.T(optional=True, xmlstyle='attribute')
640 real = FloatNoUnit.T(xmltagname='Real')
641 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
643 def value(self):
644 return self.real.value + 1J * self.imaginary.value
647class ClockDrift(FloatWithUnit):
648 unit = String.T(default='SECONDS/SAMPLE', optional=True,
649 xmlstyle='attribute') # fixed
652class Second(FloatWithUnit):
653 '''
654 A time value in seconds.
655 '''
657 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
658 # fixed unit
661class Voltage(FloatWithUnit):
662 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
663 # fixed unit
666class Angle(FloatWithUnit):
667 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
668 # fixed unit
671class Azimuth(FloatWithUnit):
672 '''
673 Instrument azimuth, degrees clockwise from North.
674 '''
676 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
677 # fixed unit
680class Dip(FloatWithUnit):
681 '''
682 Instrument dip in degrees down from horizontal. Together azimuth and dip
683 describe the direction of the sensitive axis of the instrument.
684 '''
686 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
687 # fixed unit
690class Distance(FloatWithUnit):
691 '''
692 Extension of FloatWithUnit for distances, elevations, and depths.
693 '''
695 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
696 # NOT fixed unit!
699class Frequency(FloatWithUnit):
700 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
701 # fixed unit
704class SampleRate(FloatWithUnit):
705 '''
706 Sample rate in samples per second.
707 '''
709 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
710 # fixed unit
713class Person(Object):
714 '''
715 Representation of a person's contact information. A person can belong to
716 multiple agencies and have multiple email addresses and phone numbers.
717 '''
719 name_list = List.T(Unicode.T(xmltagname='Name'))
720 agency_list = List.T(Unicode.T(xmltagname='Agency'))
721 email_list = List.T(Email.T(xmltagname='Email'))
722 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
725class FIR(BaseFilter):
726 '''
727 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
728 also commonly documented using the Coefficients element.
729 '''
731 symmetry = Symmetry.T(xmltagname='Symmetry')
732 numerator_coefficient_list = List.T(
733 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
735 def summary(self):
736 return 'fir(%i%s)' % (
737 self.get_ncoefficiencs(),
738 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
740 def get_effective_coefficients(self):
741 b = num.array(
742 [v.value for v in self.numerator_coefficient_list],
743 dtype=float)
745 if self.symmetry == 'ODD':
746 b = num.concatenate((b, b[-2::-1]))
747 elif self.symmetry == 'EVEN':
748 b = num.concatenate((b, b[::-1]))
750 return b
752 def get_effective_symmetry(self):
753 if self.symmetry == 'NONE':
754 b = self.get_effective_coefficients()
755 if num.all(b - b[::-1] == 0):
756 return ['EVEN', 'ODD'][b.size % 2]
758 return self.symmetry
760 def get_ncoefficiencs(self):
761 nf = len(self.numerator_coefficient_list)
762 if self.symmetry == 'ODD':
763 nc = nf * 2 + 1
764 elif self.symmetry == 'EVEN':
765 nc = nf * 2
766 else:
767 nc = nf
769 return nc
771 def estimate_delay(self, deltat):
772 nc = self.get_ncoefficiencs()
773 if nc > 0:
774 return deltat * (nc - 1) / 2.0
775 else:
776 return 0.0
778 def get_pyrocko_response(
779 self, context, deltat, delay_responses, normalization_frequency):
781 context += self.summary()
783 if not self.numerator_coefficient_list:
784 return Delivery([])
786 b = self.get_effective_coefficients()
788 log = []
789 drop_phase = self.get_effective_symmetry() != 'NONE'
791 if not deltat:
792 log.append((
793 'error',
794 'Digital response requires knowledge about sampling '
795 'interval. Response will be unusable.',
796 context))
798 resp = response.DigitalFilterResponse(
799 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
801 if normalization_frequency is not None and deltat is not None:
802 normalization_frequency = 0.0
803 normalization = num.abs(evaluate1(resp, normalization_frequency))
805 if num.abs(normalization - 1.0) > 1e-2:
806 log.append((
807 'warning',
808 'FIR filter coefficients are not normalized. Normalizing '
809 'them. Factor: %g' % normalization, context))
811 resp = response.DigitalFilterResponse(
812 (b/normalization).tolist(), [1.0], deltat,
813 drop_phase=drop_phase)
815 resps = [resp]
817 if not drop_phase:
818 resps.extend(delay_responses)
820 return Delivery(resps, log=log)
823class Coefficients(BaseFilter):
824 '''
825 Response: coefficients for FIR filter. Laplace transforms or IIR filters
826 can be expressed using type as well but the PolesAndZeros should be used
827 instead. Corresponds to SEED blockette 54.
828 '''
830 cf_transfer_function_type = CfTransferFunction.T(
831 xmltagname='CfTransferFunctionType')
832 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
833 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
835 def summary(self):
836 return 'coef_%s(%i,%i%s)' % (
837 'ABC?'[
838 CfTransferFunction.choices.index(
839 self.cf_transfer_function_type)],
840 len(self.numerator_list),
841 len(self.denominator_list),
842 ',sym' if self.is_symmetric_fir else '')
844 def estimate_delay(self, deltat):
845 nc = len(self.numerator_list)
846 if nc > 0:
847 return deltat * (len(self.numerator_list) - 1) / 2.0
848 else:
849 return 0.0
851 def is_symmetric_fir(self):
852 if len(self.denominator_list) != 0:
853 return False
854 b = [v.value for v in self.numerator_list]
855 return b == b[::-1]
857 def get_pyrocko_response(
858 self, context, deltat, delay_responses, normalization_frequency):
860 context += self.summary()
862 factor = 1.0
863 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
864 factor = 2.0*math.pi
866 if not self.numerator_list and not self.denominator_list:
867 return Delivery(payload=[])
869 b = num.array(
870 [v.value*factor for v in self.numerator_list], dtype=float)
872 a = num.array(
873 [1.0] + [v.value*factor for v in self.denominator_list],
874 dtype=float)
876 log = []
877 resps = []
878 if self.cf_transfer_function_type in [
879 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
880 resps.append(response.AnalogFilterResponse(b, a))
882 elif self.cf_transfer_function_type == 'DIGITAL':
883 if not deltat:
884 log.append((
885 'error',
886 'Digital response requires knowledge about sampling '
887 'interval. Response will be unusable.',
888 context))
890 drop_phase = self.is_symmetric_fir()
891 resp = response.DigitalFilterResponse(
892 b, a, deltat or 0.0, drop_phase=drop_phase)
894 if normalization_frequency is not None and deltat is not None:
895 normalization = num.abs(evaluate1(
896 resp, normalization_frequency))
898 if num.abs(normalization - 1.0) > 1e-2:
899 log.append((
900 'warning',
901 'FIR filter coefficients are not normalized. '
902 'Normalizing them. Factor: %g' % normalization,
903 context))
905 resp = response.DigitalFilterResponse(
906 (b/normalization).tolist(), [1.0], deltat,
907 drop_phase=drop_phase)
909 resps.append(resp)
911 if not drop_phase:
912 resps.extend(delay_responses)
914 else:
915 return Delivery(error=(
916 'ValueError',
917 'Unknown transfer function type: %s' % (
918 self.cf_transfer_function_type)))
920 return Delivery(payload=resps, log=log)
923class Latitude(FloatWithUnit):
924 '''
925 Type for latitude coordinate.
926 '''
928 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
929 # fixed unit
930 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
933class Longitude(FloatWithUnit):
934 '''
935 Type for longitude coordinate.
936 '''
938 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
939 # fixed unit
940 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
943class PolesZeros(BaseFilter):
944 '''
945 Response: complex poles and zeros. Corresponds to SEED blockette 53.
946 '''
948 pz_transfer_function_type = PzTransferFunction.T(
949 xmltagname='PzTransferFunctionType')
950 normalization_factor = Float.T(default=1.0,
951 xmltagname='NormalizationFactor')
952 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
953 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
954 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
956 def summary(self):
957 return 'pz_%s(%i,%i)' % (
958 'ABC?'[
959 PzTransferFunction.choices.index(
960 self.pz_transfer_function_type)],
961 len(self.pole_list),
962 len(self.zero_list))
964 def get_pyrocko_response(self, context='', deltat=None):
966 context += self.summary()
968 factor = 1.0
969 cfactor = 1.0
970 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
971 factor = 2. * math.pi
972 cfactor = (2. * math.pi)**(
973 len(self.pole_list) - len(self.zero_list))
975 log = []
976 if self.normalization_factor is None \
977 or self.normalization_factor == 0.0:
979 log.append((
980 'warning',
981 'No pole-zero normalization factor given. '
982 'Assuming a value of 1.0',
983 context))
985 nfactor = 1.0
986 else:
987 nfactor = self.normalization_factor
989 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
990 if not is_digital:
991 resp = response.PoleZeroResponse(
992 constant=nfactor*cfactor,
993 zeros=[z.value()*factor for z in self.zero_list],
994 poles=[p.value()*factor for p in self.pole_list])
995 else:
996 if not deltat:
997 log.append((
998 'error',
999 'Digital response requires knowledge about sampling '
1000 'interval. Response will be unusable.',
1001 context))
1003 resp = response.DigitalPoleZeroResponse(
1004 constant=nfactor*cfactor,
1005 zeros=[z.value()*factor for z in self.zero_list],
1006 poles=[p.value()*factor for p in self.pole_list],
1007 deltat=deltat or 0.0)
1009 if not self.normalization_frequency.value:
1010 log.append((
1011 'warning',
1012 'Cannot check pole-zero normalization factor: '
1013 'No normalization frequency given.',
1014 context))
1016 else:
1017 if is_digital and not deltat:
1018 log.append((
1019 'warning',
1020 'Cannot check computed vs reported normalization '
1021 'factor without knowing the sampling interval.',
1022 context))
1023 else:
1024 computed_normalization_factor = nfactor / abs(evaluate1(
1025 resp, self.normalization_frequency.value))
1027 db = 20.0 * num.log10(
1028 computed_normalization_factor / nfactor)
1030 if abs(db) > 0.17:
1031 log.append((
1032 'warning',
1033 'Computed and reported normalization factors differ '
1034 'by %g dB: computed: %g, reported: %g' % (
1035 db,
1036 computed_normalization_factor,
1037 nfactor),
1038 context))
1040 return Delivery([resp], log)
1043class ResponseListElement(Object):
1044 frequency = Frequency.T(xmltagname='Frequency')
1045 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1046 phase = Angle.T(xmltagname='Phase')
1049class Polynomial(BaseFilter):
1050 '''
1051 Response: expressed as a polynomial (allows non-linear sensors to be
1052 described). Corresponds to SEED blockette 62. Can be used to describe a
1053 stage of acquisition or a complete system.
1054 '''
1056 approximation_type = Approximation.T(default='MACLAURIN',
1057 xmltagname='ApproximationType')
1058 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1059 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1060 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1061 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1062 maximum_error = Float.T(xmltagname='MaximumError')
1063 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1065 def summary(self):
1066 return 'poly(%i)' % len(self.coefficient_list)
1069class Decimation(Object):
1070 '''
1071 Corresponds to SEED blockette 57.
1072 '''
1074 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1075 factor = Int.T(xmltagname='Factor')
1076 offset = Int.T(xmltagname='Offset')
1077 delay = FloatWithUnit.T(xmltagname='Delay')
1078 correction = FloatWithUnit.T(xmltagname='Correction')
1080 def summary(self):
1081 return 'deci(%i, %g -> %g, %g)' % (
1082 self.factor,
1083 self.input_sample_rate.value,
1084 self.input_sample_rate.value / self.factor,
1085 self.delay.value)
1087 def get_pyrocko_response(self):
1088 if self.delay and self.delay.value != 0.0:
1089 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1090 else:
1091 return Delivery([])
1094class Operator(Object):
1095 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1096 contact_list = List.T(Person.T(xmltagname='Contact'))
1097 web_site = String.T(optional=True, xmltagname='WebSite')
1100class Comment(Object):
1101 '''
1102 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1103 and 59.
1104 '''
1106 id = Counter.T(optional=True, xmlstyle='attribute')
1107 value = Unicode.T(xmltagname='Value')
1108 begin_effective_time = DummyAwareOptionalTimestamp.T(
1109 optional=True,
1110 xmltagname='BeginEffectiveTime')
1111 end_effective_time = DummyAwareOptionalTimestamp.T(
1112 optional=True,
1113 xmltagname='EndEffectiveTime')
1114 author_list = List.T(Person.T(xmltagname='Author'))
1117class ResponseList(BaseFilter):
1118 '''
1119 Response: list of frequency, amplitude and phase values. Corresponds to
1120 SEED blockette 55.
1121 '''
1123 response_list_element_list = List.T(
1124 ResponseListElement.T(xmltagname='ResponseListElement'))
1126 def summary(self):
1127 return 'list(%i)' % len(self.response_list_element_list)
1130class Log(Object):
1131 '''
1132 Container for log entries.
1133 '''
1135 entry_list = List.T(Comment.T(xmltagname='Entry'))
1138class ResponseStage(Object):
1139 '''
1140 This complex type represents channel response and covers SEED blockettes 53
1141 to 56.
1142 '''
1144 number = Counter.T(xmlstyle='attribute')
1145 resource_id = String.T(optional=True, xmlstyle='attribute')
1146 poles_zeros_list = List.T(
1147 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1148 coefficients_list = List.T(
1149 Coefficients.T(optional=True, xmltagname='Coefficients'))
1150 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1151 fir = FIR.T(optional=True, xmltagname='FIR')
1152 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1153 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1154 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1156 def summary(self):
1157 elements = []
1159 for stuff in [
1160 self.poles_zeros_list, self.coefficients_list,
1161 self.response_list, self.fir, self.polynomial,
1162 self.decimation, self.stage_gain]:
1164 if stuff:
1165 if isinstance(stuff, Object):
1166 elements.append(stuff.summary())
1167 else:
1168 elements.extend(obj.summary() for obj in stuff)
1170 return '%i: %s %s -> %s' % (
1171 self.number,
1172 ', '.join(elements),
1173 self.input_units.name.upper() if self.input_units else '?',
1174 self.output_units.name.upper() if self.output_units else '?')
1176 def get_squirrel_response_stage(self, context):
1177 from pyrocko.squirrel.model import ResponseStage
1179 delivery = Delivery()
1180 delivery_pr = self.get_pyrocko_response(context)
1181 log = delivery_pr.log
1182 delivery_pr.log = []
1183 elements = delivery.extend_without_payload(delivery_pr)
1185 delivery.payload = [ResponseStage(
1186 input_quantity=to_quantity(self.input_units, context, delivery),
1187 output_quantity=to_quantity(self.output_units, context, delivery),
1188 input_sample_rate=self.input_sample_rate,
1189 output_sample_rate=self.output_sample_rate,
1190 elements=elements,
1191 log=log)]
1193 return delivery
1195 def get_pyrocko_response(self, context, gain_only=False):
1197 context = context + ', stage %i' % self.number
1199 responses = []
1200 log = []
1201 if self.stage_gain:
1202 normalization_frequency = self.stage_gain.frequency or 0.0
1203 else:
1204 normalization_frequency = 0.0
1206 if not gain_only:
1207 deltat = None
1208 delay_responses = []
1209 if self.decimation:
1210 rate = self.decimation.input_sample_rate.value
1211 if rate > 0.0:
1212 deltat = 1.0 / rate
1213 delivery = self.decimation.get_pyrocko_response()
1214 if delivery.errors:
1215 return delivery
1217 delay_responses = delivery.payload
1218 log.extend(delivery.log)
1220 for pzs in self.poles_zeros_list:
1221 delivery = pzs.get_pyrocko_response(context, deltat)
1222 if delivery.errors:
1223 return delivery
1225 pz_resps = delivery.payload
1226 log.extend(delivery.log)
1227 responses.extend(pz_resps)
1229 # emulate incorrect? evalresp behaviour
1230 if pzs.normalization_frequency != normalization_frequency \
1231 and normalization_frequency != 0.0:
1233 try:
1234 trial = response.MultiplyResponse(pz_resps)
1235 anorm = num.abs(evaluate1(
1236 trial, pzs.normalization_frequency.value))
1237 asens = num.abs(
1238 evaluate1(trial, normalization_frequency))
1240 factor = anorm/asens
1242 if abs(factor - 1.0) > 0.01:
1243 log.append((
1244 'warning',
1245 'PZ normalization frequency (%g) is different '
1246 'from stage gain frequency (%s) -> Emulating '
1247 'possibly incorrect evalresp behaviour. '
1248 'Correction factor: %g' % (
1249 pzs.normalization_frequency.value,
1250 normalization_frequency,
1251 factor),
1252 context))
1254 responses.append(
1255 response.PoleZeroResponse(constant=factor))
1256 except response.InvalidResponseError as e:
1257 log.append((
1258 'warning',
1259 'Could not check response: %s' % str(e),
1260 context))
1262 if len(self.poles_zeros_list) > 1:
1263 log.append((
1264 'warning',
1265 'Multiple poles and zeros records in single response '
1266 'stage.',
1267 context))
1269 for cfs in self.coefficients_list + (
1270 [self.fir] if self.fir else []):
1272 delivery = cfs.get_pyrocko_response(
1273 context, deltat, delay_responses,
1274 normalization_frequency)
1276 if delivery.errors:
1277 return delivery
1279 responses.extend(delivery.payload)
1280 log.extend(delivery.log)
1282 if len(self.coefficients_list) > 1:
1283 log.append((
1284 'warning',
1285 'Multiple filter coefficients lists in single response '
1286 'stage.',
1287 context))
1289 if self.response_list:
1290 log.append((
1291 'warning',
1292 'Unhandled response element of type: ResponseList',
1293 context))
1295 if self.polynomial:
1296 log.append((
1297 'warning',
1298 'Unhandled response element of type: Polynomial',
1299 context))
1301 if self.stage_gain:
1302 responses.append(
1303 response.PoleZeroResponse(constant=self.stage_gain.value))
1305 return Delivery(responses, log)
1307 @property
1308 def input_units(self):
1309 for e in (self.poles_zeros_list + self.coefficients_list +
1310 [self.response_list, self.fir, self.polynomial]):
1311 if e is not None:
1312 return e.input_units
1314 return None
1316 @property
1317 def output_units(self):
1318 for e in (self.poles_zeros_list + self.coefficients_list +
1319 [self.response_list, self.fir, self.polynomial]):
1320 if e is not None:
1321 return e.output_units
1323 return None
1325 @property
1326 def input_sample_rate(self):
1327 if self.decimation:
1328 return self.decimation.input_sample_rate.value
1330 return None
1332 @property
1333 def output_sample_rate(self):
1334 if self.decimation:
1335 return self.decimation.input_sample_rate.value \
1336 / self.decimation.factor
1338 return None
1341class Response(Object):
1342 resource_id = String.T(optional=True, xmlstyle='attribute')
1343 instrument_sensitivity = Sensitivity.T(optional=True,
1344 xmltagname='InstrumentSensitivity')
1345 instrument_polynomial = Polynomial.T(optional=True,
1346 xmltagname='InstrumentPolynomial')
1347 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1349 @property
1350 def output_sample_rate(self):
1351 if self.stage_list:
1352 return self.stage_list[-1].output_sample_rate
1353 else:
1354 return None
1356 def check_sample_rates(self, channel):
1358 if self.stage_list:
1359 sample_rate = None
1361 for stage in self.stage_list:
1362 if stage.decimation:
1363 input_sample_rate = \
1364 stage.decimation.input_sample_rate.value
1366 if sample_rate is not None and not same_sample_rate(
1367 sample_rate, input_sample_rate):
1369 logger.warning(
1370 'Response stage %i has unexpected input sample '
1371 'rate: %g Hz (expected: %g Hz)' % (
1372 stage.number,
1373 input_sample_rate,
1374 sample_rate))
1376 sample_rate = input_sample_rate / stage.decimation.factor
1378 if sample_rate is not None and channel.sample_rate \
1379 and not same_sample_rate(
1380 sample_rate, channel.sample_rate.value):
1382 logger.warning(
1383 'Channel sample rate (%g Hz) does not match sample rate '
1384 'deducted from response stages information (%g Hz).' % (
1385 channel.sample_rate.value,
1386 sample_rate))
1388 def check_units(self):
1390 if self.instrument_sensitivity \
1391 and self.instrument_sensitivity.input_units:
1393 units = self.instrument_sensitivity.input_units.name.upper()
1395 if self.stage_list:
1396 for stage in self.stage_list:
1397 if units and stage.input_units \
1398 and stage.input_units.name.upper() != units:
1400 logger.warning(
1401 'Input units of stage %i (%s) do not match %s (%s).'
1402 % (
1403 stage.number,
1404 units,
1405 'output units of stage %i'
1406 if stage.number == 0
1407 else 'sensitivity input units',
1408 units))
1410 if stage.output_units:
1411 units = stage.output_units.name.upper()
1412 else:
1413 units = None
1415 sout_units = self.instrument_sensitivity.output_units
1416 if self.instrument_sensitivity and sout_units:
1417 if units is not None and units != sout_units.name.upper():
1418 logger.warning(
1419 'Output units of stage %i (%s) do not match %s (%s).'
1420 % (
1421 stage.number,
1422 units,
1423 'sensitivity output units',
1424 sout_units.name.upper()))
1426 def _sensitivity_checkpoints(self, responses, context):
1427 delivery = Delivery()
1429 if self.instrument_sensitivity:
1430 sval = self.instrument_sensitivity.value
1431 sfreq = self.instrument_sensitivity.frequency
1432 if sval is None:
1433 delivery.log.append((
1434 'warning',
1435 'No sensitivity value given.',
1436 context))
1438 elif sval is None:
1439 delivery.log.append((
1440 'warning',
1441 'Reported sensitivity value is zero.',
1442 context))
1444 elif sfreq is None:
1445 delivery.log.append((
1446 'warning',
1447 'Sensitivity frequency not given.',
1448 context))
1450 else:
1451 trial = response.MultiplyResponse(responses)
1453 delivery.extend(
1454 check_resp(
1455 trial, sval, sfreq, 0.1,
1456 'Instrument sensitivity value inconsistent with '
1457 'sensitivity computed from complete response',
1458 context))
1460 delivery.payload.append(response.FrequencyResponseCheckpoint(
1461 frequency=sfreq, value=sval))
1463 return delivery
1465 def get_squirrel_response(self, context, **kwargs):
1466 from pyrocko.squirrel.model import Response
1468 if self.stage_list:
1469 delivery = Delivery()
1470 for istage, stage in enumerate(self.stage_list):
1471 delivery.extend(stage.get_squirrel_response_stage(context))
1473 checkpoints = []
1474 if not delivery.errors:
1475 all_responses = []
1476 for sq_stage in delivery.payload:
1477 all_responses.extend(sq_stage.elements)
1479 checkpoints.extend(delivery.extend_without_payload(
1480 self._sensitivity_checkpoints(all_responses, context)))
1482 sq_stages = delivery.payload
1483 if sq_stages:
1484 if sq_stages[0].input_quantity is None \
1485 and self.instrument_sensitivity is not None:
1487 sq_stages[0].input_quantity = to_quantity(
1488 self.instrument_sensitivity.input_units,
1489 context, delivery)
1490 sq_stages[-1].output_quantity = to_quantity(
1491 self.instrument_sensitivity.output_units,
1492 context, delivery)
1494 sq_stages = delivery.expect()
1496 return Response(
1497 stages=sq_stages,
1498 log=delivery.log,
1499 checkpoints=checkpoints,
1500 **kwargs)
1502 elif self.instrument_sensitivity:
1503 raise NoResponseInformation(
1504 "Only instrument sensitivity given (won't use it). (%s)."
1505 % context)
1506 else:
1507 raise NoResponseInformation(
1508 'Empty instrument response (%s).'
1509 % context)
1511 def get_pyrocko_response(
1512 self, context, fake_input_units=None, stages=(0, 1)):
1514 delivery = Delivery()
1515 if self.stage_list:
1516 for istage, stage in enumerate(self.stage_list):
1517 delivery.extend(stage.get_pyrocko_response(
1518 context, gain_only=not (
1519 stages is None or stages[0] <= istage < stages[1])))
1521 elif self.instrument_sensitivity:
1522 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1524 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1525 checkpoints = delivery.extend_without_payload(delivery_cp)
1526 if delivery.errors:
1527 return delivery
1529 if fake_input_units is not None:
1530 if not self.instrument_sensitivity or \
1531 self.instrument_sensitivity.input_units is None:
1533 delivery.errors.append((
1534 'NoResponseInformation',
1535 'No input units given, so cannot convert to requested '
1536 'units: %s' % fake_input_units.upper(),
1537 context))
1539 return delivery
1541 input_units = self.instrument_sensitivity.input_units.name.upper()
1543 conresp = None
1544 try:
1545 conresp = conversion[
1546 fake_input_units.upper(), input_units]
1548 except KeyError:
1549 delivery.errors.append((
1550 'NoResponseInformation',
1551 'Cannot convert between units: %s, %s'
1552 % (fake_input_units, input_units),
1553 context))
1555 if conresp is not None:
1556 delivery.payload.append(conresp)
1557 for checkpoint in checkpoints:
1558 checkpoint.value *= num.abs(evaluate1(
1559 conresp, checkpoint.frequency))
1561 delivery.payload = [
1562 response.MultiplyResponse(
1563 delivery.payload,
1564 checkpoints=checkpoints)]
1566 return delivery
1568 @classmethod
1569 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1570 normalization_frequency=1.0):
1571 '''
1572 Convert Pyrocko pole-zero response to StationXML response.
1574 :param presponse: Pyrocko pole-zero response
1575 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1576 :param input_unit: Input unit to be reported in the StationXML
1577 response.
1578 :type input_unit: str
1579 :param output_unit: Output unit to be reported in the StationXML
1580 response.
1581 :type output_unit: str
1582 :param normalization_frequency: Frequency where the normalization
1583 factor for the StationXML response should be computed.
1584 :type normalization_frequency: float
1585 '''
1587 norm_factor = 1.0/float(abs(
1588 evaluate1(presponse, normalization_frequency)
1589 / presponse.constant))
1591 pzs = PolesZeros(
1592 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1593 normalization_factor=norm_factor,
1594 normalization_frequency=Frequency(normalization_frequency),
1595 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1596 imaginary=FloatNoUnit(z.imag))
1597 for z in presponse.zeros],
1598 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1599 imaginary=FloatNoUnit(z.imag))
1600 for z in presponse.poles])
1602 pzs.validate()
1604 stage = ResponseStage(
1605 number=1,
1606 poles_zeros_list=[pzs],
1607 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1609 resp = Response(
1610 instrument_sensitivity=Sensitivity(
1611 value=stage.stage_gain.value,
1612 frequency=normalization_frequency,
1613 input_units=Units(input_unit),
1614 output_units=Units(output_unit)),
1616 stage_list=[stage])
1618 return resp
1621class BaseNode(Object):
1622 '''
1623 A base node type for derivation from: Network, Station and Channel types.
1624 '''
1626 code = String.T(xmlstyle='attribute')
1627 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1628 xmlstyle='attribute')
1629 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1630 xmlstyle='attribute')
1631 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1632 alternate_code = String.T(optional=True, xmlstyle='attribute')
1633 historical_code = String.T(optional=True, xmlstyle='attribute')
1634 description = Unicode.T(optional=True, xmltagname='Description')
1635 comment_list = List.T(Comment.T(xmltagname='Comment'))
1637 def spans(self, *args):
1638 if len(args) == 0:
1639 return True
1640 elif len(args) == 1:
1641 return ((self.start_date is None or
1642 self.start_date <= args[0]) and
1643 (self.end_date is None or
1644 args[0] <= self.end_date))
1646 elif len(args) == 2:
1647 return ((self.start_date is None or
1648 args[1] >= self.start_date) and
1649 (self.end_date is None or
1650 self.end_date >= args[0]))
1653def overlaps(a, b):
1654 return (
1655 a.start_date is None or b.end_date is None
1656 or a.start_date < b.end_date
1657 ) and (
1658 b.start_date is None or a.end_date is None
1659 or b.start_date < a.end_date
1660 )
1663def check_overlaps(node_type_name, codes, nodes):
1664 errors = []
1665 for ia, a in enumerate(nodes):
1666 for b in nodes[ia+1:]:
1667 if overlaps(a, b):
1668 errors.append(
1669 '%s epochs overlap for %s:\n'
1670 ' %s - %s\n %s - %s' % (
1671 node_type_name,
1672 '.'.join(codes),
1673 tts(a.start_date), tts(a.end_date),
1674 tts(b.start_date), tts(b.end_date)))
1676 return errors
1679class Channel(BaseNode):
1680 '''
1681 Equivalent to SEED blockette 52 and parent element for the related the
1682 response blockettes.
1683 '''
1685 location_code = String.T(xmlstyle='attribute')
1686 external_reference_list = List.T(
1687 ExternalReference.T(xmltagname='ExternalReference'))
1688 latitude = Latitude.T(xmltagname='Latitude')
1689 longitude = Longitude.T(xmltagname='Longitude')
1690 elevation = Distance.T(xmltagname='Elevation')
1691 depth = Distance.T(xmltagname='Depth')
1692 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1693 dip = Dip.T(optional=True, xmltagname='Dip')
1694 type_list = List.T(Type.T(xmltagname='Type'))
1695 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1696 sample_rate_ratio = SampleRateRatio.T(optional=True,
1697 xmltagname='SampleRateRatio')
1698 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1699 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1700 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1701 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1702 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1703 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1704 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1705 response = Response.T(optional=True, xmltagname='Response')
1707 @property
1708 def position_values(self):
1709 lat = self.latitude.value
1710 lon = self.longitude.value
1711 elevation = value_or_none(self.elevation)
1712 depth = value_or_none(self.depth)
1713 return lat, lon, elevation, depth
1716class Station(BaseNode):
1717 '''
1718 This type represents a Station epoch. It is common to only have a single
1719 station epoch with the station's creation and termination dates as the
1720 epoch start and end dates.
1721 '''
1723 latitude = Latitude.T(xmltagname='Latitude')
1724 longitude = Longitude.T(xmltagname='Longitude')
1725 elevation = Distance.T(xmltagname='Elevation')
1726 site = Site.T(optional=True, xmltagname='Site')
1727 vault = Unicode.T(optional=True, xmltagname='Vault')
1728 geology = Unicode.T(optional=True, xmltagname='Geology')
1729 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1730 operator_list = List.T(Operator.T(xmltagname='Operator'))
1731 creation_date = DummyAwareOptionalTimestamp.T(
1732 optional=True, xmltagname='CreationDate')
1733 termination_date = DummyAwareOptionalTimestamp.T(
1734 optional=True, xmltagname='TerminationDate')
1735 total_number_channels = Counter.T(
1736 optional=True, xmltagname='TotalNumberChannels')
1737 selected_number_channels = Counter.T(
1738 optional=True, xmltagname='SelectedNumberChannels')
1739 external_reference_list = List.T(
1740 ExternalReference.T(xmltagname='ExternalReference'))
1741 channel_list = List.T(Channel.T(xmltagname='Channel'))
1743 @property
1744 def position_values(self):
1745 lat = self.latitude.value
1746 lon = self.longitude.value
1747 elevation = value_or_none(self.elevation)
1748 return lat, lon, elevation
1751class Network(BaseNode):
1752 '''
1753 This type represents the Network layer, all station metadata is contained
1754 within this element. The official name of the network or other descriptive
1755 information can be included in the Description element. The Network can
1756 contain 0 or more Stations.
1757 '''
1759 total_number_stations = Counter.T(optional=True,
1760 xmltagname='TotalNumberStations')
1761 selected_number_stations = Counter.T(optional=True,
1762 xmltagname='SelectedNumberStations')
1763 station_list = List.T(Station.T(xmltagname='Station'))
1765 @property
1766 def station_code_list(self):
1767 return sorted(set(s.code for s in self.station_list))
1769 @property
1770 def sl_code_list(self):
1771 sls = set()
1772 for station in self.station_list:
1773 for channel in station.channel_list:
1774 sls.add((station.code, channel.location_code))
1776 return sorted(sls)
1778 def summary(self, width=80, indent=4):
1779 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1780 lines = ['%s (%i):' % (self.code, len(sls))]
1781 if sls:
1782 ssls = ['.'.join(x for x in c if x) for c in sls]
1783 w = max(len(x) for x in ssls)
1784 n = (width - indent) / (w+1)
1785 while ssls:
1786 lines.append(
1787 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1789 ssls[:n] = []
1791 return '\n'.join(lines)
1794def value_or_none(x):
1795 if x is not None:
1796 return x.value
1797 else:
1798 return None
1801def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1803 pos = lat, lon, elevation, depth = \
1804 channels[0].position_values
1806 if not all(pos == x.position_values for x in channels):
1807 info = '\n'.join(
1808 ' %s: %s' % (x.code, x.position_values) for
1809 x in channels)
1811 mess = 'encountered inconsistencies in channel ' \
1812 'lat/lon/elevation/depth ' \
1813 'for %s.%s.%s: \n%s' % (nsl + (info,))
1815 if inconsistencies == 'raise':
1816 raise InconsistentChannelLocations(mess)
1818 elif inconsistencies == 'warn':
1819 logger.warning(mess)
1820 logger.warning(' -> using mean values')
1822 apos = num.array([x.position_values for x in channels], dtype=float)
1823 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1824 / num.sum(num.isfinite(apos), axis=0)
1826 groups = {}
1827 for channel in channels:
1828 if channel.code not in groups:
1829 groups[channel.code] = []
1831 groups[channel.code].append(channel)
1833 pchannels = []
1834 for code in sorted(groups.keys()):
1835 data = [
1836 (channel.code, value_or_none(channel.azimuth),
1837 value_or_none(channel.dip))
1838 for channel in groups[code]]
1840 azimuth, dip = util.consistency_merge(
1841 data,
1842 message='channel orientation values differ:',
1843 error=inconsistencies)
1845 pchannels.append(
1846 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1848 return pyrocko.model.Station(
1849 *nsl,
1850 lat=mlat,
1851 lon=mlon,
1852 elevation=mele,
1853 depth=mdep,
1854 channels=pchannels)
1857class FDSNStationXML(Object):
1858 '''
1859 Top-level type for Station XML. Required field are Source (network ID of
1860 the institution sending the message) and one or more Network containers or
1861 one or more Station containers.
1862 '''
1864 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1865 source = String.T(xmltagname='Source')
1866 sender = String.T(optional=True, xmltagname='Sender')
1867 module = String.T(optional=True, xmltagname='Module')
1868 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1869 created = Timestamp.T(optional=True, xmltagname='Created')
1870 network_list = List.T(Network.T(xmltagname='Network'))
1872 xmltagname = 'FDSNStationXML'
1873 guessable_xmlns = [guts_xmlns]
1875 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1876 time=None, timespan=None,
1877 inconsistencies='warn'):
1879 assert inconsistencies in ('raise', 'warn')
1881 if nslcs is not None:
1882 nslcs = set(nslcs)
1884 if nsls is not None:
1885 nsls = set(nsls)
1887 tt = ()
1888 if time is not None:
1889 tt = (time,)
1890 elif timespan is not None:
1891 tt = timespan
1893 pstations = []
1894 for network in self.network_list:
1895 if not network.spans(*tt):
1896 continue
1898 for station in network.station_list:
1899 if not station.spans(*tt):
1900 continue
1902 if station.channel_list:
1903 loc_to_channels = {}
1904 for channel in station.channel_list:
1905 if not channel.spans(*tt):
1906 continue
1908 loc = channel.location_code.strip()
1909 if loc not in loc_to_channels:
1910 loc_to_channels[loc] = []
1912 loc_to_channels[loc].append(channel)
1914 for loc in sorted(loc_to_channels.keys()):
1915 channels = loc_to_channels[loc]
1916 if nslcs is not None:
1917 channels = [channel for channel in channels
1918 if (network.code, station.code, loc,
1919 channel.code) in nslcs]
1921 if not channels:
1922 continue
1924 nsl = network.code, station.code, loc
1925 if nsls is not None and nsl not in nsls:
1926 continue
1928 pstations.append(
1929 pyrocko_station_from_channels(
1930 nsl,
1931 channels,
1932 inconsistencies=inconsistencies))
1933 else:
1934 pstations.append(pyrocko.model.Station(
1935 network.code, station.code, '*',
1936 lat=station.latitude.value,
1937 lon=station.longitude.value,
1938 elevation=value_or_none(station.elevation),
1939 name=station.description or ''))
1941 return pstations
1943 @classmethod
1944 def from_pyrocko_stations(
1945 cls, pyrocko_stations, add_flat_responses_from=None):
1947 '''
1948 Generate :py:class:`FDSNStationXML` from list of
1949 :py:class;`pyrocko.model.Station` instances.
1951 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1952 instances.
1953 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1954 '''
1955 from collections import defaultdict
1956 network_dict = defaultdict(list)
1958 if add_flat_responses_from:
1959 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1960 extra = dict(
1961 response=Response(
1962 instrument_sensitivity=Sensitivity(
1963 value=1.0,
1964 frequency=1.0,
1965 input_units=Units(name=add_flat_responses_from))))
1966 else:
1967 extra = {}
1969 have_offsets = set()
1970 for s in pyrocko_stations:
1972 if s.north_shift != 0.0 or s.east_shift != 0.0:
1973 have_offsets.add(s.nsl())
1975 network, station, location = s.nsl()
1976 channel_list = []
1977 for c in s.channels:
1978 channel_list.append(
1979 Channel(
1980 location_code=location,
1981 code=c.name,
1982 latitude=Latitude(value=s.effective_lat),
1983 longitude=Longitude(value=s.effective_lon),
1984 elevation=Distance(value=s.elevation),
1985 depth=Distance(value=s.depth),
1986 azimuth=Azimuth(value=c.azimuth),
1987 dip=Dip(value=c.dip),
1988 **extra
1989 )
1990 )
1992 network_dict[network].append(
1993 Station(
1994 code=station,
1995 latitude=Latitude(value=s.effective_lat),
1996 longitude=Longitude(value=s.effective_lon),
1997 elevation=Distance(value=s.elevation),
1998 channel_list=channel_list)
1999 )
2001 if have_offsets:
2002 logger.warning(
2003 'StationXML does not support Cartesian offsets in '
2004 'coordinates. Storing effective lat/lon for stations: %s' %
2005 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2007 timestamp = util.to_time_float(time.time())
2008 network_list = []
2009 for k, station_list in network_dict.items():
2011 network_list.append(
2012 Network(
2013 code=k, station_list=station_list,
2014 total_number_stations=len(station_list)))
2016 sxml = FDSNStationXML(
2017 source='from pyrocko stations list', created=timestamp,
2018 network_list=network_list)
2020 sxml.validate()
2021 return sxml
2023 def iter_network_stations(
2024 self, net=None, sta=None, time=None, timespan=None):
2026 tt = ()
2027 if time is not None:
2028 tt = (time,)
2029 elif timespan is not None:
2030 tt = timespan
2032 for network in self.network_list:
2033 if not network.spans(*tt) or (
2034 net is not None and network.code != net):
2035 continue
2037 for station in network.station_list:
2038 if not station.spans(*tt) or (
2039 sta is not None and station.code != sta):
2040 continue
2042 yield (network, station)
2044 def iter_network_station_channels(
2045 self, net=None, sta=None, loc=None, cha=None,
2046 time=None, timespan=None):
2048 if loc is not None:
2049 loc = loc.strip()
2051 tt = ()
2052 if time is not None:
2053 tt = (time,)
2054 elif timespan is not None:
2055 tt = timespan
2057 for network in self.network_list:
2058 if not network.spans(*tt) or (
2059 net is not None and network.code != net):
2060 continue
2062 for station in network.station_list:
2063 if not station.spans(*tt) or (
2064 sta is not None and station.code != sta):
2065 continue
2067 if station.channel_list:
2068 for channel in station.channel_list:
2069 if (not channel.spans(*tt) or
2070 (cha is not None and channel.code != cha) or
2071 (loc is not None and
2072 channel.location_code.strip() != loc)):
2073 continue
2075 yield (network, station, channel)
2077 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2078 time=None, timespan=None):
2080 groups = {}
2081 for network, station, channel in self.iter_network_station_channels(
2082 net, sta, loc, cha, time=time, timespan=timespan):
2084 net = network.code
2085 sta = station.code
2086 cha = channel.code
2087 loc = channel.location_code.strip()
2088 if len(cha) == 3:
2089 bic = cha[:2] # band and intrument code according to SEED
2090 elif len(cha) == 1:
2091 bic = ''
2092 else:
2093 bic = cha
2095 if channel.response and \
2096 channel.response.instrument_sensitivity and \
2097 channel.response.instrument_sensitivity.input_units:
2099 unit = channel.response.instrument_sensitivity\
2100 .input_units.name.upper()
2101 else:
2102 unit = None
2104 bic = (bic, unit)
2106 k = net, sta, loc
2107 if k not in groups:
2108 groups[k] = {}
2110 if bic not in groups[k]:
2111 groups[k][bic] = []
2113 groups[k][bic].append(channel)
2115 for nsl, bic_to_channels in groups.items():
2116 bad_bics = []
2117 for bic, channels in bic_to_channels.items():
2118 sample_rates = []
2119 for channel in channels:
2120 sample_rates.append(channel.sample_rate.value)
2122 if not same(sample_rates):
2123 scs = ','.join(channel.code for channel in channels)
2124 srs = ', '.join('%e' % x for x in sample_rates)
2125 err = 'ignoring channels with inconsistent sampling ' + \
2126 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2128 logger.warning(err)
2129 bad_bics.append(bic)
2131 for bic in bad_bics:
2132 del bic_to_channels[bic]
2134 return groups
2136 def choose_channels(
2137 self,
2138 target_sample_rate=None,
2139 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2140 priority_units=['M/S', 'M/S**2'],
2141 priority_instrument_code=['H', 'L'],
2142 time=None,
2143 timespan=None):
2145 nslcs = {}
2146 for nsl, bic_to_channels in self.get_channel_groups(
2147 time=time, timespan=timespan).items():
2149 useful_bics = []
2150 for bic, channels in bic_to_channels.items():
2151 rate = channels[0].sample_rate.value
2153 if target_sample_rate is not None and \
2154 rate < target_sample_rate*0.99999:
2155 continue
2157 if len(bic[0]) == 2:
2158 if bic[0][0] not in priority_band_code:
2159 continue
2161 if bic[0][1] not in priority_instrument_code:
2162 continue
2164 unit = bic[1]
2166 prio_unit = len(priority_units)
2167 try:
2168 prio_unit = priority_units.index(unit)
2169 except ValueError:
2170 pass
2172 prio_inst = len(priority_instrument_code)
2173 prio_band = len(priority_band_code)
2174 if len(channels[0].code) == 3:
2175 try:
2176 prio_inst = priority_instrument_code.index(
2177 channels[0].code[1])
2178 except ValueError:
2179 pass
2181 try:
2182 prio_band = priority_band_code.index(
2183 channels[0].code[0])
2184 except ValueError:
2185 pass
2187 if target_sample_rate is None:
2188 rate = -rate
2190 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2191 prio_inst, bic))
2193 useful_bics.sort()
2195 for _, _, rate, _, _, bic in useful_bics:
2196 channels = sorted(
2197 bic_to_channels[bic],
2198 key=lambda channel: channel.code)
2200 if channels:
2201 for channel in channels:
2202 nslcs[nsl + (channel.code,)] = channel
2204 break
2206 return nslcs
2208 def get_pyrocko_response(
2209 self, nslc,
2210 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2212 net, sta, loc, cha = nslc
2213 resps = []
2214 for _, _, channel in self.iter_network_station_channels(
2215 net, sta, loc, cha, time=time, timespan=timespan):
2216 resp = channel.response
2217 if resp:
2218 resp.check_sample_rates(channel)
2219 resp.check_units()
2220 resps.append(resp.get_pyrocko_response(
2221 '.'.join(nslc),
2222 fake_input_units=fake_input_units,
2223 stages=stages).expect_one())
2225 if not resps:
2226 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2227 elif len(resps) > 1:
2228 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2230 return resps[0]
2232 @property
2233 def n_code_list(self):
2234 return sorted(set(x.code for x in self.network_list))
2236 @property
2237 def ns_code_list(self):
2238 nss = set()
2239 for network in self.network_list:
2240 for station in network.station_list:
2241 nss.add((network.code, station.code))
2243 return sorted(nss)
2245 @property
2246 def nsl_code_list(self):
2247 nsls = set()
2248 for network in self.network_list:
2249 for station in network.station_list:
2250 for channel in station.channel_list:
2251 nsls.add(
2252 (network.code, station.code, channel.location_code))
2254 return sorted(nsls)
2256 @property
2257 def nslc_code_list(self):
2258 nslcs = set()
2259 for network in self.network_list:
2260 for station in network.station_list:
2261 for channel in station.channel_list:
2262 nslcs.add(
2263 (network.code, station.code, channel.location_code,
2264 channel.code))
2266 return sorted(nslcs)
2268 def summary(self):
2269 lst = [
2270 'number of n codes: %i' % len(self.n_code_list),
2271 'number of ns codes: %i' % len(self.ns_code_list),
2272 'number of nsl codes: %i' % len(self.nsl_code_list),
2273 'number of nslc codes: %i' % len(self.nslc_code_list)
2274 ]
2275 return '\n'.join(lst)
2277 def summary_stages(self):
2278 data = []
2279 for network, station, channel in self.iter_network_station_channels():
2280 nslc = (network.code, station.code, channel.location_code,
2281 channel.code)
2283 stages = []
2284 in_units = '?'
2285 out_units = '?'
2286 if channel.response:
2287 sens = channel.response.instrument_sensitivity
2288 if sens:
2289 in_units = sens.input_units.name.upper()
2290 out_units = sens.output_units.name.upper()
2292 for stage in channel.response.stage_list:
2293 stages.append(stage.summary())
2295 data.append(
2296 (nslc, tts(channel.start_date), tts(channel.end_date),
2297 in_units, out_units, stages))
2299 data.sort()
2301 lst = []
2302 for nslc, stmin, stmax, in_units, out_units, stages in data:
2303 lst.append(' %s: %s - %s, %s -> %s' % (
2304 '.'.join(nslc), stmin, stmax, in_units, out_units))
2305 for stage in stages:
2306 lst.append(' %s' % stage)
2308 return '\n'.join(lst)
2310 def _check_overlaps(self):
2311 by_nslc = {}
2312 for network in self.network_list:
2313 for station in network.station_list:
2314 for channel in station.channel_list:
2315 nslc = (network.code, station.code, channel.location_code,
2316 channel.code)
2317 if nslc not in by_nslc:
2318 by_nslc[nslc] = []
2320 by_nslc[nslc].append(channel)
2322 errors = []
2323 for nslc, channels in by_nslc.items():
2324 errors.extend(check_overlaps('Channel', nslc, channels))
2326 return errors
2328 def check(self):
2329 errors = []
2330 for meth in [self._check_overlaps]:
2331 errors.extend(meth())
2333 if errors:
2334 raise Inconsistencies(
2335 'Inconsistencies found in StationXML:\n '
2336 + '\n '.join(errors))
2339def load_channel_table(stream):
2341 networks = {}
2342 stations = {}
2344 for line in stream:
2345 line = str(line.decode('ascii'))
2346 if line.startswith('#'):
2347 continue
2349 t = line.rstrip().split('|')
2351 if len(t) != 17:
2352 logger.warning('Invalid channel record: %s' % line)
2353 continue
2355 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2356 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2358 try:
2359 scale = float(scale)
2360 except ValueError:
2361 scale = None
2363 try:
2364 scale_freq = float(scale_freq)
2365 except ValueError:
2366 scale_freq = None
2368 try:
2369 depth = float(dep)
2370 except ValueError:
2371 depth = 0.0
2373 try:
2374 azi = float(azi)
2375 dip = float(dip)
2376 except ValueError:
2377 azi = None
2378 dip = None
2380 try:
2381 if net not in networks:
2382 network = Network(code=net)
2383 else:
2384 network = networks[net]
2386 if (net, sta) not in stations:
2387 station = Station(
2388 code=sta, latitude=lat, longitude=lon, elevation=ele)
2390 station.regularize()
2391 else:
2392 station = stations[net, sta]
2394 if scale:
2395 resp = Response(
2396 instrument_sensitivity=Sensitivity(
2397 value=scale,
2398 frequency=scale_freq,
2399 input_units=scale_units))
2400 else:
2401 resp = None
2403 channel = Channel(
2404 code=cha,
2405 location_code=loc.strip(),
2406 latitude=lat,
2407 longitude=lon,
2408 elevation=ele,
2409 depth=depth,
2410 azimuth=azi,
2411 dip=dip,
2412 sensor=Equipment(description=sens),
2413 response=resp,
2414 sample_rate=sample_rate,
2415 start_date=start_date,
2416 end_date=end_date or None)
2418 channel.regularize()
2420 except ValidationError:
2421 raise InvalidRecord(line)
2423 if net not in networks:
2424 networks[net] = network
2426 if (net, sta) not in stations:
2427 stations[net, sta] = station
2428 network.station_list.append(station)
2430 station.channel_list.append(channel)
2432 return FDSNStationXML(
2433 source='created from table input',
2434 created=time.time(),
2435 network_list=sorted(networks.values(), key=lambda x: x.code))
2438def primitive_merge(sxs):
2439 networks = []
2440 for sx in sxs:
2441 networks.extend(sx.network_list)
2443 return FDSNStationXML(
2444 source='merged from different sources',
2445 created=time.time(),
2446 network_list=copy.deepcopy(
2447 sorted(networks, key=lambda x: x.code)))
2450def split_channels(sx):
2451 for nslc in sx.nslc_code_list:
2452 network_list = sx.network_list
2453 network_list_filtered = [
2454 network for network in network_list
2455 if network.code == nslc[0]]
2457 for network in network_list_filtered:
2458 sx.network_list = [network]
2459 station_list = network.station_list
2460 station_list_filtered = [
2461 station for station in station_list
2462 if station.code == nslc[1]]
2464 for station in station_list_filtered:
2465 network.station_list = [station]
2466 channel_list = station.channel_list
2467 station.channel_list = [
2468 channel for channel in channel_list
2469 if (channel.location_code, channel.code) == nslc[2:4]]
2471 yield nslc, copy.deepcopy(sx)
2473 station.channel_list = channel_list
2475 network.station_list = station_list
2477 sx.network_list = network_list
2480if __name__ == '__main__':
2481 from optparse import OptionParser
2483 util.setup_logging('pyrocko.io.stationxml', 'warning')
2485 usage = \
2486 'python -m pyrocko.io.stationxml check|stats|stages ' \
2487 '<filename> [options]'
2489 description = '''Torture StationXML file.'''
2491 parser = OptionParser(
2492 usage=usage,
2493 description=description,
2494 formatter=util.BetterHelpFormatter())
2496 (options, args) = parser.parse_args(sys.argv[1:])
2498 if len(args) != 2:
2499 parser.print_help()
2500 sys.exit(1)
2502 action, path = args
2504 sx = load_xml(filename=path)
2505 if action == 'check':
2506 try:
2507 sx.check()
2508 except Inconsistencies as e:
2509 logger.error(e)
2510 sys.exit(1)
2512 elif action == 'stats':
2513 print(sx.summary())
2515 elif action == 'stages':
2516 print(sx.summary_stages())
2518 else:
2519 parser.print_help()
2520 sys.exit('unknown action: %s' % action)