1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import sys
8import time
9import logging
10import datetime
11import calendar
12import math
13import copy
15import numpy as num
17from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
18 Unicode, Int, Float, List, Object, Timestamp,
19 ValidationError, TBase, re_tz, Any, Tuple)
20from pyrocko.guts import load_xml # noqa
21from pyrocko.util import hpfloat, time_to_str, get_time_float
23import pyrocko.model
24from pyrocko import util, response
26try:
27 newstr = unicode
28except NameError:
29 newstr = str
31guts_prefix = 'sx'
33guts_xmlns = 'http://www.fdsn.org/xml/station/1'
35logger = logging.getLogger('pyrocko.io.stationxml')
37conversion = {
38 ('M', 'M'): None,
39 ('M/S', 'M'): response.IntegrationResponse(1),
40 ('M/S**2', 'M'): response.IntegrationResponse(2),
41 ('M', 'M/S'): response.DifferentiationResponse(1),
42 ('M/S', 'M/S'): None,
43 ('M/S**2', 'M/S'): response.IntegrationResponse(1),
44 ('M', 'M/S**2'): response.DifferentiationResponse(2),
45 ('M/S', 'M/S**2'): response.DifferentiationResponse(1),
46 ('M/S**2', 'M/S**2'): None}
49unit_to_quantity = {
50 'M/S': 'velocity',
51 'M': 'displacement',
52 'M/S**2': 'acceleration',
53 'V': 'voltage',
54 'COUNTS': 'counts',
55 'COUNT': 'counts',
56 'PA': 'pressure',
57 'RAD/S': 'rotation'}
60def to_quantity(unit, context, delivery):
62 if unit is None:
63 return None
65 name = unit.name.upper()
66 if name in unit_to_quantity:
67 return unit_to_quantity[name]
68 else:
69 delivery.log.append((
70 'warning',
71 'Units not supported by Squirrel framework: %s' % (
72 unit.name.upper() + (
73 ' (%s)' % unit.description
74 if unit.description else '')),
75 context))
77 return 'unsupported_quantity(%s)' % unit
80class StationXMLError(Exception):
81 pass
84class Inconsistencies(StationXMLError):
85 pass
88class NoResponseInformation(StationXMLError):
89 pass
92class MultipleResponseInformation(StationXMLError):
93 pass
96class InconsistentResponseInformation(StationXMLError):
97 pass
100class InconsistentChannelLocations(StationXMLError):
101 pass
104class InvalidRecord(StationXMLError):
105 def __init__(self, line):
106 StationXMLError.__init__(self)
107 self._line = line
109 def __str__(self):
110 return 'Invalid record: "%s"' % self._line
113_exceptions = dict(
114 Inconsistencies=Inconsistencies,
115 NoResponseInformation=NoResponseInformation,
116 MultipleResponseInformation=MultipleResponseInformation,
117 InconsistentResponseInformation=InconsistentResponseInformation,
118 InconsistentChannelLocations=InconsistentChannelLocations,
119 InvalidRecord=InvalidRecord,
120 ValueError=ValueError)
123_logs = dict(
124 info=logger.info,
125 warning=logger.warning,
126 error=logger.error)
129class DeliveryError(StationXMLError):
130 pass
133class Delivery(Object):
135 def __init__(self, payload=None, log=None, errors=None, error=None):
136 if payload is None:
137 payload = []
139 if log is None:
140 log = []
142 if errors is None:
143 errors = []
145 if error is not None:
146 errors.append(error)
148 Object.__init__(self, payload=payload, log=log, errors=errors)
150 payload = List.T(Any.T())
151 log = List.T(Tuple.T(3, String.T()))
152 errors = List.T(Tuple.T(3, String.T()))
154 def extend(self, other):
155 self.payload.extend(other.payload)
156 self.log.extend(other.log)
157 self.errors.extend(other.errors)
159 def extend_without_payload(self, other):
160 self.log.extend(other.log)
161 self.errors.extend(other.errors)
162 return other.payload
164 def emit_log(self):
165 for name, message, context in self.log:
166 message = '%s: %s' % (context, message)
167 _logs[name](message)
169 def expect(self, quiet=False):
170 if not quiet:
171 self.emit_log()
173 if self.errors:
174 name, message, context = self.errors[0]
175 if context:
176 message += ' (%s)' % context
178 if len(self.errors) > 1:
179 message += ' Additional errors pending.'
181 raise _exceptions[name](message)
183 return self.payload
185 def expect_one(self, quiet=False):
186 payload = self.expect(quiet=quiet)
187 if len(payload) != 1:
188 raise DeliveryError(
189 'Expected 1 element but got %i.' % len(payload))
191 return payload[0]
194def wrap(s, width=80, indent=4):
195 words = s.split()
196 lines = []
197 t = []
198 n = 0
199 for w in words:
200 if n + len(w) >= width:
201 lines.append(' '.join(t))
202 n = indent
203 t = [' '*(indent-1)]
205 t.append(w)
206 n += len(w) + 1
208 lines.append(' '.join(t))
209 return '\n'.join(lines)
212def same(x, eps=0.0):
213 if any(type(x[0]) != type(r) for r in x):
214 return False
216 if isinstance(x[0], float):
217 return all(abs(r-x[0]) <= eps for r in x)
218 else:
219 return all(r == x[0] for r in x)
222def same_sample_rate(a, b, eps=1.0e-6):
223 return abs(a - b) < min(a, b)*eps
226def evaluate1(resp, f):
227 return resp.evaluate(num.array([f], dtype=float))[0]
230def check_resp(resp, value, frequency, limit_db, prelude, context):
232 try:
233 value_resp = num.abs(evaluate1(resp, frequency))
234 except response.InvalidResponseError as e:
235 return Delivery(log=[(
236 'warning',
237 'Could not check response: %s' % str(e),
238 context)])
240 if value_resp == 0.0:
241 return Delivery(log=[(
242 'warning',
243 '%s\n'
244 ' computed response is zero' % prelude,
245 context)])
247 diff_db = 20.0 * num.log10(value_resp/value)
249 if num.abs(diff_db) > limit_db:
250 return Delivery(log=[(
251 'warning',
252 '%s\n'
253 ' reported value: %g\n'
254 ' computed value: %g\n'
255 ' at frequency [Hz]: %g\n'
256 ' factor: %g\n'
257 ' difference [dB]: %g\n'
258 ' limit [dB]: %g' % (
259 prelude,
260 value,
261 value_resp,
262 frequency,
263 value_resp/value,
264 diff_db,
265 limit_db),
266 context)])
268 return Delivery()
271def tts(t):
272 if t is None:
273 return '?'
274 else:
275 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
278def le_open_left(a, b):
279 return a is None or (b is not None and a <= b)
282def le_open_right(a, b):
283 return b is None or (a is not None and a <= b)
286def eq_open(a, b):
287 return (a is None and b is None) \
288 or (a is not None and b is not None and a == b)
291def contains(a, b):
292 return le_open_left(a.start_date, b.start_date) \
293 and le_open_right(b.end_date, a.end_date)
296def find_containing(candidates, node):
297 for candidate in candidates:
298 if contains(candidate, node):
299 return candidate
301 return None
304this_year = time.gmtime()[0]
307class DummyAwareOptionalTimestamp(Object):
308 '''
309 Optional timestamp with support for some common placeholder values.
311 Some StationXML files contain arbitrary placeholder values for open end
312 intervals, like "2100-01-01". Depending on the time range supported by the
313 system, these dates are translated into ``None`` to prevent crashes with
314 this type.
315 '''
316 dummy_for = (hpfloat, float)
317 dummy_for_description = 'time_float'
319 class __T(TBase):
321 def regularize_extra(self, val):
322 time_float = get_time_float()
324 if isinstance(val, datetime.datetime):
325 tt = val.utctimetuple()
326 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
328 elif isinstance(val, datetime.date):
329 tt = val.timetuple()
330 val = time_float(calendar.timegm(tt))
332 elif isinstance(val, (str, newstr)):
333 val = val.strip()
335 tz_offset = 0
337 m = re_tz.search(val)
338 if m:
339 sh = m.group(2)
340 sm = m.group(4)
341 tz_offset = (int(sh)*3600 if sh else 0) \
342 + (int(sm)*60 if sm else 0)
344 val = re_tz.sub('', val)
346 if len(val) > 10 and val[10] == 'T':
347 val = val.replace('T', ' ', 1)
349 try:
350 val = util.str_to_time(val) - tz_offset
352 except util.TimeStrError:
353 year = int(val[:4])
354 if sys.maxsize > 2**32: # if we're on 64bit
355 if year > this_year + 100:
356 return None # StationXML contained a dummy date
358 if year < 1903: # for macOS, 1900-01-01 dummy dates
359 return None
361 else: # 32bit end of time is in 2038
362 if this_year < 2037 and year > 2037 or year < 1903:
363 return None # StationXML contained a dummy date
365 raise
367 elif isinstance(val, (int, float)):
368 val = time_float(val)
370 else:
371 raise ValidationError(
372 '%s: cannot convert "%s" to type %s' % (
373 self.xname(), val, time_float))
375 return val
377 def to_save(self, val):
378 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
379 .rstrip('0').rstrip('.')
381 def to_save_xml(self, val):
382 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
383 .rstrip('0').rstrip('.') + 'Z'
386class Nominal(StringChoice):
387 choices = [
388 'NOMINAL',
389 'CALCULATED']
392class Email(UnicodePattern):
393 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
396class RestrictedStatus(StringChoice):
397 choices = [
398 'open',
399 'closed',
400 'partial']
403class Type(StringChoice):
404 choices = [
405 'TRIGGERED',
406 'CONTINUOUS',
407 'HEALTH',
408 'GEOPHYSICAL',
409 'WEATHER',
410 'FLAG',
411 'SYNTHESIZED',
412 'INPUT',
413 'EXPERIMENTAL',
414 'MAINTENANCE',
415 'BEAM']
417 class __T(StringChoice.T):
418 def validate_extra(self, val):
419 if val not in self.choices:
420 logger.warning(
421 'channel type: "%s" is not a valid choice out of %s' %
422 (val, repr(self.choices)))
425class PzTransferFunction(StringChoice):
426 choices = [
427 'LAPLACE (RADIANS/SECOND)',
428 'LAPLACE (HERTZ)',
429 'DIGITAL (Z-TRANSFORM)']
432class Symmetry(StringChoice):
433 choices = [
434 'NONE',
435 'EVEN',
436 'ODD']
439class CfTransferFunction(StringChoice):
441 class __T(StringChoice.T):
442 def validate(self, val, regularize=False, depth=-1):
443 if regularize:
444 try:
445 val = str(val)
446 except ValueError:
447 raise ValidationError(
448 '%s: cannot convert to string %s' % (self.xname,
449 repr(val)))
451 val = self.dummy_cls.replacements.get(val, val)
453 self.validate_extra(val)
454 return val
456 choices = [
457 'ANALOG (RADIANS/SECOND)',
458 'ANALOG (HERTZ)',
459 'DIGITAL']
461 replacements = {
462 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
463 'ANALOG (HZ)': 'ANALOG (HERTZ)',
464 }
467class Approximation(StringChoice):
468 choices = [
469 'MACLAURIN']
472class PhoneNumber(StringPattern):
473 pattern = '[0-9]+-[0-9]+'
476class Site(Object):
477 '''
478 Description of a site location using name and optional geopolitical
479 boundaries (country, city, etc.).
480 '''
482 name = Unicode.T(default='', xmltagname='Name')
483 description = Unicode.T(optional=True, xmltagname='Description')
484 town = Unicode.T(optional=True, xmltagname='Town')
485 county = Unicode.T(optional=True, xmltagname='County')
486 region = Unicode.T(optional=True, xmltagname='Region')
487 country = Unicode.T(optional=True, xmltagname='Country')
490class ExternalReference(Object):
491 '''
492 This type contains a URI and description for external data that users may
493 want to reference in StationXML.
494 '''
496 uri = String.T(xmltagname='URI')
497 description = Unicode.T(xmltagname='Description')
500class Units(Object):
501 '''
502 A type to document units. Corresponds to SEED blockette 34.
503 '''
505 def __init__(self, name=None, **kwargs):
506 Object.__init__(self, name=name, **kwargs)
508 name = String.T(xmltagname='Name')
509 description = Unicode.T(optional=True, xmltagname='Description')
512class Counter(Int):
513 pass
516class SampleRateRatio(Object):
517 '''
518 Sample rate expressed as number of samples in a number of seconds.
519 '''
521 number_samples = Int.T(xmltagname='NumberSamples')
522 number_seconds = Int.T(xmltagname='NumberSeconds')
525class Gain(Object):
526 '''
527 Complex type for sensitivity and frequency ranges. This complex type can be
528 used to represent both overall sensitivities and individual stage gains.
529 The FrequencyRangeGroup is an optional construct that defines a pass band
530 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
531 is valid within the number of decibels specified in FrequencyDBVariation.
532 '''
534 def __init__(self, value=None, **kwargs):
535 Object.__init__(self, value=value, **kwargs)
537 value = Float.T(optional=True, xmltagname='Value')
538 frequency = Float.T(optional=True, xmltagname='Frequency')
540 def summary(self):
541 return 'gain(%g @ %g)' % (self.value, self.frequency)
544class NumeratorCoefficient(Object):
545 i = Int.T(optional=True, xmlstyle='attribute')
546 value = Float.T(xmlstyle='content')
549class FloatNoUnit(Object):
550 def __init__(self, value=None, **kwargs):
551 Object.__init__(self, value=value, **kwargs)
553 plus_error = Float.T(optional=True, xmlstyle='attribute')
554 minus_error = Float.T(optional=True, xmlstyle='attribute')
555 value = Float.T(xmlstyle='content')
558class FloatWithUnit(FloatNoUnit):
559 unit = String.T(optional=True, xmlstyle='attribute')
562class Equipment(Object):
563 resource_id = String.T(optional=True, xmlstyle='attribute')
564 type = String.T(optional=True, xmltagname='Type')
565 description = Unicode.T(optional=True, xmltagname='Description')
566 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
567 vendor = Unicode.T(optional=True, xmltagname='Vendor')
568 model = Unicode.T(optional=True, xmltagname='Model')
569 serial_number = String.T(optional=True, xmltagname='SerialNumber')
570 installation_date = DummyAwareOptionalTimestamp.T(
571 optional=True,
572 xmltagname='InstallationDate')
573 removal_date = DummyAwareOptionalTimestamp.T(
574 optional=True,
575 xmltagname='RemovalDate')
576 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
579class PhoneNumber(Object):
580 description = Unicode.T(optional=True, xmlstyle='attribute')
581 country_code = Int.T(optional=True, xmltagname='CountryCode')
582 area_code = Int.T(xmltagname='AreaCode')
583 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
586class BaseFilter(Object):
587 '''
588 The BaseFilter is derived by all filters.
589 '''
591 resource_id = String.T(optional=True, xmlstyle='attribute')
592 name = String.T(optional=True, xmlstyle='attribute')
593 description = Unicode.T(optional=True, xmltagname='Description')
594 input_units = Units.T(optional=True, xmltagname='InputUnits')
595 output_units = Units.T(optional=True, xmltagname='OutputUnits')
598class Sensitivity(Gain):
599 '''
600 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
601 construct that defines a pass band in Hertz (FrequencyStart and
602 FrequencyEnd) in which the SensitivityValue is valid within the number of
603 decibels specified in FrequencyDBVariation.
604 '''
606 input_units = Units.T(optional=True, xmltagname='InputUnits')
607 output_units = Units.T(optional=True, xmltagname='OutputUnits')
608 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
609 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
610 frequency_db_variation = Float.T(optional=True,
611 xmltagname='FrequencyDBVariation')
613 def get_pyrocko_response(self):
614 return Delivery(
615 [response.PoleZeroResponse(constant=self.value)])
618class Coefficient(FloatNoUnit):
619 number = Counter.T(optional=True, xmlstyle='attribute')
622class PoleZero(Object):
623 '''
624 Complex numbers used as poles or zeros in channel response.
625 '''
627 number = Int.T(optional=True, xmlstyle='attribute')
628 real = FloatNoUnit.T(xmltagname='Real')
629 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
631 def value(self):
632 return self.real.value + 1J * self.imaginary.value
635class ClockDrift(FloatWithUnit):
636 unit = String.T(default='SECONDS/SAMPLE', optional=True,
637 xmlstyle='attribute') # fixed
640class Second(FloatWithUnit):
641 '''
642 A time value in seconds.
643 '''
645 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
646 # fixed unit
649class Voltage(FloatWithUnit):
650 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
651 # fixed unit
654class Angle(FloatWithUnit):
655 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
656 # fixed unit
659class Azimuth(FloatWithUnit):
660 '''
661 Instrument azimuth, degrees clockwise from North.
662 '''
664 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
665 # fixed unit
668class Dip(FloatWithUnit):
669 '''
670 Instrument dip in degrees down from horizontal. Together azimuth and dip
671 describe the direction of the sensitive axis of the instrument.
672 '''
674 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
675 # fixed unit
678class Distance(FloatWithUnit):
679 '''
680 Extension of FloatWithUnit for distances, elevations, and depths.
681 '''
683 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
684 # NOT fixed unit!
687class Frequency(FloatWithUnit):
688 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
689 # fixed unit
692class SampleRate(FloatWithUnit):
693 '''
694 Sample rate in samples per second.
695 '''
697 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
698 # fixed unit
701class Person(Object):
702 '''
703 Representation of a person's contact information. A person can belong to
704 multiple agencies and have multiple email addresses and phone numbers.
705 '''
707 name_list = List.T(Unicode.T(xmltagname='Name'))
708 agency_list = List.T(Unicode.T(xmltagname='Agency'))
709 email_list = List.T(Email.T(xmltagname='Email'))
710 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
713class FIR(BaseFilter):
714 '''
715 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
716 also commonly documented using the Coefficients element.
717 '''
719 symmetry = Symmetry.T(xmltagname='Symmetry')
720 numerator_coefficient_list = List.T(
721 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
723 def summary(self):
724 return 'fir(%i%s)' % (
725 self.get_ncoefficiencs(),
726 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
728 def get_effective_coefficients(self):
729 b = num.array(
730 [v.value for v in self.numerator_coefficient_list],
731 dtype=float)
733 if self.symmetry == 'ODD':
734 b = num.concatenate((b, b[-2::-1]))
735 elif self.symmetry == 'EVEN':
736 b = num.concatenate((b, b[::-1]))
738 return b
740 def get_effective_symmetry(self):
741 if self.symmetry == 'NONE':
742 b = self.get_effective_coefficients()
743 if num.all(b - b[::-1] == 0):
744 return ['EVEN', 'ODD'][b.size % 2]
746 return self.symmetry
748 def get_ncoefficiencs(self):
749 nf = len(self.numerator_coefficient_list)
750 if self.symmetry == 'ODD':
751 nc = nf * 2 + 1
752 elif self.symmetry == 'EVEN':
753 nc = nf * 2
754 else:
755 nc = nf
757 return nc
759 def estimate_delay(self, deltat):
760 nc = self.get_ncoefficiencs()
761 if nc > 0:
762 return deltat * (nc - 1) / 2.0
763 else:
764 return 0.0
766 def get_pyrocko_response(
767 self, context, deltat, delay_responses, normalization_frequency):
769 context += self.summary()
771 if not self.numerator_coefficient_list:
772 return Delivery([])
774 b = self.get_effective_coefficients()
776 log = []
777 drop_phase = self.get_effective_symmetry() != 'NONE'
779 if not deltat:
780 log.append((
781 'error',
782 'Digital response requires knowledge about sampling '
783 'interval. Response will be unusable.',
784 context))
786 resp = response.DigitalFilterResponse(
787 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
789 if normalization_frequency is not None and deltat is not None:
790 normalization_frequency = 0.0
791 normalization = num.abs(evaluate1(resp, normalization_frequency))
793 if num.abs(normalization - 1.0) > 1e-2:
794 log.append((
795 'warning',
796 'FIR filter coefficients are not normalized. Normalizing '
797 'them. Factor: %g' % normalization, context))
799 resp = response.DigitalFilterResponse(
800 (b/normalization).tolist(), [1.0], deltat,
801 drop_phase=drop_phase)
803 resps = [resp]
805 if not drop_phase:
806 resps.extend(delay_responses)
808 return Delivery(resps, log=log)
811class Coefficients(BaseFilter):
812 '''
813 Response: coefficients for FIR filter. Laplace transforms or IIR filters
814 can be expressed using type as well but the PolesAndZeros should be used
815 instead. Corresponds to SEED blockette 54.
816 '''
818 cf_transfer_function_type = CfTransferFunction.T(
819 xmltagname='CfTransferFunctionType')
820 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
821 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
823 def summary(self):
824 return 'coef_%s(%i,%i%s)' % (
825 'ABC?'[
826 CfTransferFunction.choices.index(
827 self.cf_transfer_function_type)],
828 len(self.numerator_list),
829 len(self.denominator_list),
830 ',sym' if self.is_symmetric_fir else '')
832 def estimate_delay(self, deltat):
833 nc = len(self.numerator_list)
834 if nc > 0:
835 return deltat * (len(self.numerator_list) - 1) / 2.0
836 else:
837 return 0.0
839 def is_symmetric_fir(self):
840 if len(self.denominator_list) != 0:
841 return False
842 b = [v.value for v in self.numerator_list]
843 return b == b[::-1]
845 def get_pyrocko_response(
846 self, context, deltat, delay_responses, normalization_frequency):
848 context += self.summary()
850 factor = 1.0
851 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
852 factor = 2.0*math.pi
854 if not self.numerator_list and not self.denominator_list:
855 return Delivery(payload=[])
857 b = num.array(
858 [v.value*factor for v in self.numerator_list], dtype=float)
860 a = num.array(
861 [1.0] + [v.value*factor for v in self.denominator_list],
862 dtype=float)
864 log = []
865 resps = []
866 if self.cf_transfer_function_type in [
867 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
868 resps.append(response.AnalogFilterResponse(b, a))
870 elif self.cf_transfer_function_type == 'DIGITAL':
871 if not deltat:
872 log.append((
873 'error',
874 'Digital response requires knowledge about sampling '
875 'interval. Response will be unusable.',
876 context))
878 drop_phase = self.is_symmetric_fir()
879 resp = response.DigitalFilterResponse(
880 b, a, deltat or 0.0, drop_phase=drop_phase)
882 if normalization_frequency is not None and deltat is not None:
883 normalization = num.abs(evaluate1(
884 resp, normalization_frequency))
886 if num.abs(normalization - 1.0) > 1e-2:
887 log.append((
888 'warning',
889 'FIR filter coefficients are not normalized. '
890 'Normalizing them. Factor: %g' % normalization,
891 context))
893 resp = response.DigitalFilterResponse(
894 (b/normalization).tolist(), [1.0], deltat,
895 drop_phase=drop_phase)
897 resps.append(resp)
899 if not drop_phase:
900 resps.extend(delay_responses)
902 else:
903 return Delivery(error=(
904 'ValueError',
905 'Unknown transfer function type: %s' % (
906 self.cf_transfer_function_type)))
908 return Delivery(payload=resps, log=log)
911class Latitude(FloatWithUnit):
912 '''
913 Type for latitude coordinate.
914 '''
916 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
917 # fixed unit
918 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
921class Longitude(FloatWithUnit):
922 '''
923 Type for longitude coordinate.
924 '''
926 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
927 # fixed unit
928 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
931class PolesZeros(BaseFilter):
932 '''
933 Response: complex poles and zeros. Corresponds to SEED blockette 53.
934 '''
936 pz_transfer_function_type = PzTransferFunction.T(
937 xmltagname='PzTransferFunctionType')
938 normalization_factor = Float.T(default=1.0,
939 xmltagname='NormalizationFactor')
940 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
941 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
942 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
944 def summary(self):
945 return 'pz_%s(%i,%i)' % (
946 'ABC?'[
947 PzTransferFunction.choices.index(
948 self.pz_transfer_function_type)],
949 len(self.pole_list),
950 len(self.zero_list))
952 def get_pyrocko_response(self, context='', deltat=None):
954 context += self.summary()
956 factor = 1.0
957 cfactor = 1.0
958 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
959 factor = 2. * math.pi
960 cfactor = (2. * math.pi)**(
961 len(self.pole_list) - len(self.zero_list))
963 log = []
964 if self.normalization_factor is None \
965 or self.normalization_factor == 0.0:
967 log.append((
968 'warning',
969 'No pole-zero normalization factor given. '
970 'Assuming a value of 1.0',
971 context))
973 nfactor = 1.0
974 else:
975 nfactor = self.normalization_factor
977 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
978 if not is_digital:
979 resp = response.PoleZeroResponse(
980 constant=nfactor*cfactor,
981 zeros=[z.value()*factor for z in self.zero_list],
982 poles=[p.value()*factor for p in self.pole_list])
983 else:
984 if not deltat:
985 log.append((
986 'error',
987 'Digital response requires knowledge about sampling '
988 'interval. Response will be unusable.',
989 context))
991 resp = response.DigitalPoleZeroResponse(
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 deltat=deltat or 0.0)
997 if not self.normalization_frequency.value:
998 log.append((
999 'warning',
1000 'Cannot check pole-zero normalization factor: '
1001 'No normalization frequency given.',
1002 context))
1004 else:
1005 if is_digital and not deltat:
1006 log.append((
1007 'warning',
1008 'Cannot check computed vs reported normalization '
1009 'factor without knowing the sampling interval.',
1010 context))
1011 else:
1012 computed_normalization_factor = nfactor / abs(evaluate1(
1013 resp, self.normalization_frequency.value))
1015 db = 20.0 * num.log10(
1016 computed_normalization_factor / nfactor)
1018 if abs(db) > 0.17:
1019 log.append((
1020 'warning',
1021 'Computed and reported normalization factors differ '
1022 'by %g dB: computed: %g, reported: %g' % (
1023 db,
1024 computed_normalization_factor,
1025 nfactor),
1026 context))
1028 return Delivery([resp], log)
1031class ResponseListElement(Object):
1032 frequency = Frequency.T(xmltagname='Frequency')
1033 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1034 phase = Angle.T(xmltagname='Phase')
1037class Polynomial(BaseFilter):
1038 '''
1039 Response: expressed as a polynomial (allows non-linear sensors to be
1040 described). Corresponds to SEED blockette 62. Can be used to describe a
1041 stage of acquisition or a complete system.
1042 '''
1044 approximation_type = Approximation.T(default='MACLAURIN',
1045 xmltagname='ApproximationType')
1046 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1047 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1048 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1049 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1050 maximum_error = Float.T(xmltagname='MaximumError')
1051 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1053 def summary(self):
1054 return 'poly(%i)' % len(self.coefficient_list)
1057class Decimation(Object):
1058 '''
1059 Corresponds to SEED blockette 57.
1060 '''
1062 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1063 factor = Int.T(xmltagname='Factor')
1064 offset = Int.T(xmltagname='Offset')
1065 delay = FloatWithUnit.T(xmltagname='Delay')
1066 correction = FloatWithUnit.T(xmltagname='Correction')
1068 def summary(self):
1069 return 'deci(%i, %g -> %g, %g)' % (
1070 self.factor,
1071 self.input_sample_rate.value,
1072 self.input_sample_rate.value / self.factor,
1073 self.delay.value)
1075 def get_pyrocko_response(self):
1076 if self.delay and self.delay.value != 0.0:
1077 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1078 else:
1079 return Delivery([])
1082class Operator(Object):
1083 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1084 contact_list = List.T(Person.T(xmltagname='Contact'))
1085 web_site = String.T(optional=True, xmltagname='WebSite')
1088class Comment(Object):
1089 '''
1090 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1091 and 59.
1092 '''
1094 id = Counter.T(optional=True, xmlstyle='attribute')
1095 value = Unicode.T(xmltagname='Value')
1096 begin_effective_time = DummyAwareOptionalTimestamp.T(
1097 optional=True,
1098 xmltagname='BeginEffectiveTime')
1099 end_effective_time = DummyAwareOptionalTimestamp.T(
1100 optional=True,
1101 xmltagname='EndEffectiveTime')
1102 author_list = List.T(Person.T(xmltagname='Author'))
1105class ResponseList(BaseFilter):
1106 '''
1107 Response: list of frequency, amplitude and phase values. Corresponds to
1108 SEED blockette 55.
1109 '''
1111 response_list_element_list = List.T(
1112 ResponseListElement.T(xmltagname='ResponseListElement'))
1114 def summary(self):
1115 return 'list(%i)' % len(self.response_list_element_list)
1118class Log(Object):
1119 '''
1120 Container for log entries.
1121 '''
1123 entry_list = List.T(Comment.T(xmltagname='Entry'))
1126class ResponseStage(Object):
1127 '''
1128 This complex type represents channel response and covers SEED blockettes 53
1129 to 56.
1130 '''
1132 number = Counter.T(xmlstyle='attribute')
1133 resource_id = String.T(optional=True, xmlstyle='attribute')
1134 poles_zeros_list = List.T(
1135 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1136 coefficients_list = List.T(
1137 Coefficients.T(optional=True, xmltagname='Coefficients'))
1138 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1139 fir = FIR.T(optional=True, xmltagname='FIR')
1140 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1141 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1142 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1144 def summary(self):
1145 elements = []
1147 for stuff in [
1148 self.poles_zeros_list, self.coefficients_list,
1149 self.response_list, self.fir, self.polynomial,
1150 self.decimation, self.stage_gain]:
1152 if stuff:
1153 if isinstance(stuff, Object):
1154 elements.append(stuff.summary())
1155 else:
1156 elements.extend(obj.summary() for obj in stuff)
1158 return '%i: %s %s -> %s' % (
1159 self.number,
1160 ', '.join(elements),
1161 self.input_units.name.upper() if self.input_units else '?',
1162 self.output_units.name.upper() if self.output_units else '?')
1164 def get_squirrel_response_stage(self, context):
1165 from pyrocko.squirrel.model import ResponseStage
1167 delivery = Delivery()
1168 delivery_pr = self.get_pyrocko_response(context)
1169 log = delivery_pr.log
1170 delivery_pr.log = []
1171 elements = delivery.extend_without_payload(delivery_pr)
1173 delivery.payload = [ResponseStage(
1174 input_quantity=to_quantity(self.input_units, context, delivery),
1175 output_quantity=to_quantity(self.output_units, context, delivery),
1176 input_sample_rate=self.input_sample_rate,
1177 output_sample_rate=self.output_sample_rate,
1178 elements=elements,
1179 log=log)]
1181 return delivery
1183 def get_pyrocko_response(self, context, gain_only=False):
1185 context = context + ', stage %i' % self.number
1187 responses = []
1188 log = []
1189 if self.stage_gain:
1190 normalization_frequency = self.stage_gain.frequency or 0.0
1191 else:
1192 normalization_frequency = 0.0
1194 if not gain_only:
1195 deltat = None
1196 delay_responses = []
1197 if self.decimation:
1198 rate = self.decimation.input_sample_rate.value
1199 if rate > 0.0:
1200 deltat = 1.0 / rate
1201 delivery = self.decimation.get_pyrocko_response()
1202 if delivery.errors:
1203 return delivery
1205 delay_responses = delivery.payload
1206 log.extend(delivery.log)
1208 for pzs in self.poles_zeros_list:
1209 delivery = pzs.get_pyrocko_response(context, deltat)
1210 if delivery.errors:
1211 return delivery
1213 pz_resps = delivery.payload
1214 log.extend(delivery.log)
1215 responses.extend(pz_resps)
1217 # emulate incorrect? evalresp behaviour
1218 if pzs.normalization_frequency != normalization_frequency \
1219 and normalization_frequency != 0.0:
1221 try:
1222 trial = response.MultiplyResponse(pz_resps)
1223 anorm = num.abs(evaluate1(
1224 trial, pzs.normalization_frequency.value))
1225 asens = num.abs(
1226 evaluate1(trial, normalization_frequency))
1228 factor = anorm/asens
1230 if abs(factor - 1.0) > 0.01:
1231 log.append((
1232 'warning',
1233 'PZ normalization frequency (%g) is different '
1234 'from stage gain frequency (%s) -> Emulating '
1235 'possibly incorrect evalresp behaviour. '
1236 'Correction factor: %g' % (
1237 pzs.normalization_frequency.value,
1238 normalization_frequency,
1239 factor),
1240 context))
1242 responses.append(
1243 response.PoleZeroResponse(constant=factor))
1244 except response.InvalidResponseError as e:
1245 log.append((
1246 'warning',
1247 'Could not check response: %s' % str(e),
1248 context))
1250 if len(self.poles_zeros_list) > 1:
1251 log.append((
1252 'warning',
1253 'Multiple poles and zeros records in single response '
1254 'stage.',
1255 context))
1257 for cfs in self.coefficients_list + (
1258 [self.fir] if self.fir else []):
1260 delivery = cfs.get_pyrocko_response(
1261 context, deltat, delay_responses,
1262 normalization_frequency)
1264 if delivery.errors:
1265 return delivery
1267 responses.extend(delivery.payload)
1268 log.extend(delivery.log)
1270 if len(self.coefficients_list) > 1:
1271 log.append((
1272 'warning',
1273 'Multiple filter coefficients lists in single response '
1274 'stage.',
1275 context))
1277 if self.response_list:
1278 log.append((
1279 'warning',
1280 'Unhandled response element of type: ResponseList',
1281 context))
1283 if self.polynomial:
1284 log.append((
1285 'warning',
1286 'Unhandled response element of type: Polynomial',
1287 context))
1289 if self.stage_gain:
1290 responses.append(
1291 response.PoleZeroResponse(constant=self.stage_gain.value))
1293 return Delivery(responses, log)
1295 @property
1296 def input_units(self):
1297 for e in (self.poles_zeros_list + self.coefficients_list +
1298 [self.response_list, self.fir, self.polynomial]):
1299 if e is not None:
1300 return e.input_units
1302 return None
1304 @property
1305 def output_units(self):
1306 for e in (self.poles_zeros_list + self.coefficients_list +
1307 [self.response_list, self.fir, self.polynomial]):
1308 if e is not None:
1309 return e.output_units
1311 return None
1313 @property
1314 def input_sample_rate(self):
1315 if self.decimation:
1316 return self.decimation.input_sample_rate.value
1318 return None
1320 @property
1321 def output_sample_rate(self):
1322 if self.decimation:
1323 return self.decimation.input_sample_rate.value \
1324 / self.decimation.factor
1326 return None
1329class Response(Object):
1330 resource_id = String.T(optional=True, xmlstyle='attribute')
1331 instrument_sensitivity = Sensitivity.T(optional=True,
1332 xmltagname='InstrumentSensitivity')
1333 instrument_polynomial = Polynomial.T(optional=True,
1334 xmltagname='InstrumentPolynomial')
1335 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1337 def check_sample_rates(self, channel):
1339 if self.stage_list:
1340 sample_rate = None
1342 for stage in self.stage_list:
1343 if stage.decimation:
1344 input_sample_rate = \
1345 stage.decimation.input_sample_rate.value
1347 if sample_rate is not None and not same_sample_rate(
1348 sample_rate, input_sample_rate):
1350 logger.warning(
1351 'Response stage %i has unexpected input sample '
1352 'rate: %g Hz (expected: %g Hz)' % (
1353 stage.number,
1354 input_sample_rate,
1355 sample_rate))
1357 sample_rate = input_sample_rate / stage.decimation.factor
1359 if sample_rate is not None and channel.sample_rate \
1360 and not same_sample_rate(
1361 sample_rate, channel.sample_rate.value):
1363 logger.warning(
1364 'Channel sample rate (%g Hz) does not match sample rate '
1365 'deducted from response stages information (%g Hz).' % (
1366 channel.sample_rate.value,
1367 sample_rate))
1369 def check_units(self):
1371 if self.instrument_sensitivity \
1372 and self.instrument_sensitivity.input_units:
1374 units = self.instrument_sensitivity.input_units.name.upper()
1376 if self.stage_list:
1377 for stage in self.stage_list:
1378 if units and stage.input_units \
1379 and stage.input_units.name.upper() != units:
1381 logger.warning(
1382 'Input units of stage %i (%s) do not match %s (%s).'
1383 % (
1384 stage.number,
1385 units,
1386 'output units of stage %i'
1387 if stage.number == 0
1388 else 'sensitivity input units',
1389 units))
1391 if stage.output_units:
1392 units = stage.output_units.name.upper()
1393 else:
1394 units = None
1396 sout_units = self.instrument_sensitivity.output_units
1397 if self.instrument_sensitivity and sout_units:
1398 if units is not None and units != sout_units.name.upper():
1399 logger.warning(
1400 'Output units of stage %i (%s) do not match %s (%s).'
1401 % (
1402 stage.number,
1403 units,
1404 'sensitivity output units',
1405 sout_units.name.upper()))
1407 def _sensitivity_checkpoints(self, responses, context):
1408 delivery = Delivery()
1410 if self.instrument_sensitivity:
1411 sval = self.instrument_sensitivity.value
1412 sfreq = self.instrument_sensitivity.frequency
1413 if sval is None:
1414 delivery.log.append((
1415 'warning',
1416 'No sensitivity value given.',
1417 context))
1419 elif sval is None:
1420 delivery.log.append((
1421 'warning',
1422 'Reported sensitivity value is zero.',
1423 context))
1425 elif sfreq is None:
1426 delivery.log.append((
1427 'warning',
1428 'Sensitivity frequency not given.',
1429 context))
1431 else:
1432 trial = response.MultiplyResponse(responses)
1434 delivery.extend(
1435 check_resp(
1436 trial, sval, sfreq, 0.1,
1437 'Instrument sensitivity value inconsistent with '
1438 'sensitivity computed from complete response',
1439 context))
1441 delivery.payload.append(response.FrequencyResponseCheckpoint(
1442 frequency=sfreq, value=sval))
1444 return delivery
1446 def get_squirrel_response(self, context, **kwargs):
1447 from pyrocko.squirrel.model import Response
1449 if self.stage_list:
1450 delivery = Delivery()
1451 for istage, stage in enumerate(self.stage_list):
1452 delivery.extend(stage.get_squirrel_response_stage(context))
1454 checkpoints = []
1455 if not delivery.errors:
1456 all_responses = []
1457 for sq_stage in delivery.payload:
1458 all_responses.extend(sq_stage.elements)
1460 checkpoints.extend(delivery.extend_without_payload(
1461 self._sensitivity_checkpoints(all_responses, context)))
1463 sq_stages = delivery.payload
1464 if sq_stages:
1465 if sq_stages[0].input_quantity is None \
1466 and self.instrument_sensitivity is not None:
1468 sq_stages[0].input_quantity = to_quantity(
1469 self.instrument_sensitivity.input_units,
1470 context, delivery)
1471 sq_stages[-1].output_quantity = to_quantity(
1472 self.instrument_sensitivity.output_units,
1473 context, delivery)
1475 sq_stages = delivery.expect()
1477 return Response(
1478 stages=sq_stages,
1479 log=delivery.log,
1480 checkpoints=checkpoints,
1481 **kwargs)
1483 elif self.instrument_sensitivity:
1484 raise NoResponseInformation(
1485 'Only instrument sensitivity given (won\'t use it). (%s).'
1486 % context)
1487 else:
1488 raise NoResponseInformation(
1489 'Empty instrument response (%s).'
1490 % context)
1492 def get_pyrocko_response(
1493 self, context, fake_input_units=None, stages=(0, 1)):
1495 delivery = Delivery()
1496 if self.stage_list:
1497 for istage, stage in enumerate(self.stage_list):
1498 delivery.extend(stage.get_pyrocko_response(
1499 context, gain_only=not (
1500 stages is None or stages[0] <= istage < stages[1])))
1502 elif self.instrument_sensitivity:
1503 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1505 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1506 checkpoints = delivery.extend_without_payload(delivery_cp)
1507 if delivery.errors:
1508 return delivery
1510 if fake_input_units is not None:
1511 if not self.instrument_sensitivity or \
1512 self.instrument_sensitivity.input_units is None:
1514 delivery.errors.append((
1515 'NoResponseInformation',
1516 'No input units given, so cannot convert to requested '
1517 'units: %s' % fake_input_units.upper(),
1518 context))
1520 return delivery
1522 input_units = self.instrument_sensitivity.input_units.name.upper()
1524 conresp = None
1525 try:
1526 conresp = conversion[
1527 fake_input_units.upper(), input_units]
1529 except KeyError:
1530 delivery.errors.append((
1531 'NoResponseInformation',
1532 'Cannot convert between units: %s, %s'
1533 % (fake_input_units, input_units),
1534 context))
1536 if conresp is not None:
1537 delivery.payload.append(conresp)
1538 for checkpoint in checkpoints:
1539 checkpoint.value *= num.abs(evaluate1(
1540 conresp, checkpoint.frequency))
1542 delivery.payload = [
1543 response.MultiplyResponse(
1544 delivery.payload,
1545 checkpoints=checkpoints)]
1547 return delivery
1549 @classmethod
1550 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1551 normalization_frequency=1.0):
1552 '''
1553 Convert Pyrocko pole-zero response to StationXML response.
1555 :param presponse: Pyrocko pole-zero response
1556 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1557 :param input_unit: Input unit to be reported in the StationXML
1558 response.
1559 :type input_unit: str
1560 :param output_unit: Output unit to be reported in the StationXML
1561 response.
1562 :type output_unit: str
1563 :param normalization_frequency: Frequency where the normalization
1564 factor for the StationXML response should be computed.
1565 :type normalization_frequency: float
1566 '''
1568 norm_factor = 1.0/float(abs(
1569 evaluate1(presponse, normalization_frequency)
1570 / presponse.constant))
1572 pzs = PolesZeros(
1573 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1574 normalization_factor=norm_factor,
1575 normalization_frequency=Frequency(normalization_frequency),
1576 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1577 imaginary=FloatNoUnit(z.imag))
1578 for z in presponse.zeros],
1579 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1580 imaginary=FloatNoUnit(z.imag))
1581 for z in presponse.poles])
1583 pzs.validate()
1585 stage = ResponseStage(
1586 number=1,
1587 poles_zeros_list=[pzs],
1588 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1590 resp = Response(
1591 instrument_sensitivity=Sensitivity(
1592 value=stage.stage_gain.value,
1593 frequency=normalization_frequency,
1594 input_units=Units(input_unit),
1595 output_units=Units(output_unit)),
1597 stage_list=[stage])
1599 return resp
1602class BaseNode(Object):
1603 '''
1604 A base node type for derivation from: Network, Station and Channel types.
1605 '''
1607 code = String.T(xmlstyle='attribute')
1608 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1609 xmlstyle='attribute')
1610 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1611 xmlstyle='attribute')
1612 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1613 alternate_code = String.T(optional=True, xmlstyle='attribute')
1614 historical_code = String.T(optional=True, xmlstyle='attribute')
1615 description = Unicode.T(optional=True, xmltagname='Description')
1616 comment_list = List.T(Comment.T(xmltagname='Comment'))
1618 def spans(self, *args):
1619 if len(args) == 0:
1620 return True
1621 elif len(args) == 1:
1622 return ((self.start_date is None or
1623 self.start_date <= args[0]) and
1624 (self.end_date is None or
1625 args[0] <= self.end_date))
1627 elif len(args) == 2:
1628 return ((self.start_date is None or
1629 args[1] >= self.start_date) and
1630 (self.end_date is None or
1631 self.end_date >= args[0]))
1634def overlaps(a, b):
1635 return (
1636 a.start_date is None or b.end_date is None
1637 or a.start_date < b.end_date
1638 ) and (
1639 b.start_date is None or a.end_date is None
1640 or b.start_date < a.end_date
1641 )
1644def check_overlaps(node_type_name, codes, nodes):
1645 errors = []
1646 for ia, a in enumerate(nodes):
1647 for b in nodes[ia+1:]:
1648 if overlaps(a, b):
1649 errors.append(
1650 '%s epochs overlap for %s:\n'
1651 ' %s - %s\n %s - %s' % (
1652 node_type_name,
1653 '.'.join(codes),
1654 tts(a.start_date), tts(a.end_date),
1655 tts(b.start_date), tts(b.end_date)))
1657 return errors
1660class Channel(BaseNode):
1661 '''
1662 Equivalent to SEED blockette 52 and parent element for the related the
1663 response blockettes.
1664 '''
1666 location_code = String.T(xmlstyle='attribute')
1667 external_reference_list = List.T(
1668 ExternalReference.T(xmltagname='ExternalReference'))
1669 latitude = Latitude.T(xmltagname='Latitude')
1670 longitude = Longitude.T(xmltagname='Longitude')
1671 elevation = Distance.T(xmltagname='Elevation')
1672 depth = Distance.T(xmltagname='Depth')
1673 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1674 dip = Dip.T(optional=True, xmltagname='Dip')
1675 type_list = List.T(Type.T(xmltagname='Type'))
1676 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1677 sample_rate_ratio = SampleRateRatio.T(optional=True,
1678 xmltagname='SampleRateRatio')
1679 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1680 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1681 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1682 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1683 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1684 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1685 equipment = Equipment.T(optional=True, xmltagname='Equipment')
1686 response = Response.T(optional=True, xmltagname='Response')
1688 @property
1689 def position_values(self):
1690 lat = self.latitude.value
1691 lon = self.longitude.value
1692 elevation = value_or_none(self.elevation)
1693 depth = value_or_none(self.depth)
1694 return lat, lon, elevation, depth
1697class Station(BaseNode):
1698 '''
1699 This type represents a Station epoch. It is common to only have a single
1700 station epoch with the station's creation and termination dates as the
1701 epoch start and end dates.
1702 '''
1704 latitude = Latitude.T(xmltagname='Latitude')
1705 longitude = Longitude.T(xmltagname='Longitude')
1706 elevation = Distance.T(xmltagname='Elevation')
1707 site = Site.T(optional=True, xmltagname='Site')
1708 vault = Unicode.T(optional=True, xmltagname='Vault')
1709 geology = Unicode.T(optional=True, xmltagname='Geology')
1710 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1711 operator_list = List.T(Operator.T(xmltagname='Operator'))
1712 creation_date = DummyAwareOptionalTimestamp.T(
1713 optional=True, xmltagname='CreationDate')
1714 termination_date = DummyAwareOptionalTimestamp.T(
1715 optional=True, xmltagname='TerminationDate')
1716 total_number_channels = Counter.T(
1717 optional=True, xmltagname='TotalNumberChannels')
1718 selected_number_channels = Counter.T(
1719 optional=True, xmltagname='SelectedNumberChannels')
1720 external_reference_list = List.T(
1721 ExternalReference.T(xmltagname='ExternalReference'))
1722 channel_list = List.T(Channel.T(xmltagname='Channel'))
1724 @property
1725 def position_values(self):
1726 lat = self.latitude.value
1727 lon = self.longitude.value
1728 elevation = value_or_none(self.elevation)
1729 return lat, lon, elevation
1732class Network(BaseNode):
1733 '''
1734 This type represents the Network layer, all station metadata is contained
1735 within this element. The official name of the network or other descriptive
1736 information can be included in the Description element. The Network can
1737 contain 0 or more Stations.
1738 '''
1740 total_number_stations = Counter.T(optional=True,
1741 xmltagname='TotalNumberStations')
1742 selected_number_stations = Counter.T(optional=True,
1743 xmltagname='SelectedNumberStations')
1744 station_list = List.T(Station.T(xmltagname='Station'))
1746 @property
1747 def station_code_list(self):
1748 return sorted(set(s.code for s in self.station_list))
1750 @property
1751 def sl_code_list(self):
1752 sls = set()
1753 for station in self.station_list:
1754 for channel in station.channel_list:
1755 sls.add((station.code, channel.location_code))
1757 return sorted(sls)
1759 def summary(self, width=80, indent=4):
1760 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1761 lines = ['%s (%i):' % (self.code, len(sls))]
1762 if sls:
1763 ssls = ['.'.join(x for x in c if x) for c in sls]
1764 w = max(len(x) for x in ssls)
1765 n = (width - indent) / (w+1)
1766 while ssls:
1767 lines.append(
1768 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1770 ssls[:n] = []
1772 return '\n'.join(lines)
1775def value_or_none(x):
1776 if x is not None:
1777 return x.value
1778 else:
1779 return None
1782def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1784 pos = lat, lon, elevation, depth = \
1785 channels[0].position_values
1787 if not all(pos == x.position_values for x in channels):
1788 info = '\n'.join(
1789 ' %s: %s' % (x.code, x.position_values) for
1790 x in channels)
1792 mess = 'encountered inconsistencies in channel ' \
1793 'lat/lon/elevation/depth ' \
1794 'for %s.%s.%s: \n%s' % (nsl + (info,))
1796 if inconsistencies == 'raise':
1797 raise InconsistentChannelLocations(mess)
1799 elif inconsistencies == 'warn':
1800 logger.warning(mess)
1801 logger.warning(' -> using mean values')
1803 apos = num.array([x.position_values for x in channels], dtype=float)
1804 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1805 / num.sum(num.isfinite(apos), axis=0)
1807 groups = {}
1808 for channel in channels:
1809 if channel.code not in groups:
1810 groups[channel.code] = []
1812 groups[channel.code].append(channel)
1814 pchannels = []
1815 for code in sorted(groups.keys()):
1816 data = [
1817 (channel.code, value_or_none(channel.azimuth),
1818 value_or_none(channel.dip))
1819 for channel in groups[code]]
1821 azimuth, dip = util.consistency_merge(
1822 data,
1823 message='channel orientation values differ:',
1824 error=inconsistencies)
1826 pchannels.append(
1827 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1829 return pyrocko.model.Station(
1830 *nsl,
1831 lat=mlat,
1832 lon=mlon,
1833 elevation=mele,
1834 depth=mdep,
1835 channels=pchannels)
1838class FDSNStationXML(Object):
1839 '''
1840 Top-level type for Station XML. Required field are Source (network ID of
1841 the institution sending the message) and one or more Network containers or
1842 one or more Station containers.
1843 '''
1845 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1846 source = String.T(xmltagname='Source')
1847 sender = String.T(optional=True, xmltagname='Sender')
1848 module = String.T(optional=True, xmltagname='Module')
1849 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1850 created = Timestamp.T(optional=True, xmltagname='Created')
1851 network_list = List.T(Network.T(xmltagname='Network'))
1853 xmltagname = 'FDSNStationXML'
1854 guessable_xmlns = [guts_xmlns]
1856 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1857 time=None, timespan=None,
1858 inconsistencies='warn'):
1860 assert inconsistencies in ('raise', 'warn')
1862 if nslcs is not None:
1863 nslcs = set(nslcs)
1865 if nsls is not None:
1866 nsls = set(nsls)
1868 tt = ()
1869 if time is not None:
1870 tt = (time,)
1871 elif timespan is not None:
1872 tt = timespan
1874 pstations = []
1875 for network in self.network_list:
1876 if not network.spans(*tt):
1877 continue
1879 for station in network.station_list:
1880 if not station.spans(*tt):
1881 continue
1883 if station.channel_list:
1884 loc_to_channels = {}
1885 for channel in station.channel_list:
1886 if not channel.spans(*tt):
1887 continue
1889 loc = channel.location_code.strip()
1890 if loc not in loc_to_channels:
1891 loc_to_channels[loc] = []
1893 loc_to_channels[loc].append(channel)
1895 for loc in sorted(loc_to_channels.keys()):
1896 channels = loc_to_channels[loc]
1897 if nslcs is not None:
1898 channels = [channel for channel in channels
1899 if (network.code, station.code, loc,
1900 channel.code) in nslcs]
1902 if not channels:
1903 continue
1905 nsl = network.code, station.code, loc
1906 if nsls is not None and nsl not in nsls:
1907 continue
1909 pstations.append(
1910 pyrocko_station_from_channels(
1911 nsl,
1912 channels,
1913 inconsistencies=inconsistencies))
1914 else:
1915 pstations.append(pyrocko.model.Station(
1916 network.code, station.code, '*',
1917 lat=station.latitude.value,
1918 lon=station.longitude.value,
1919 elevation=value_or_none(station.elevation),
1920 name=station.description or ''))
1922 return pstations
1924 @classmethod
1925 def from_pyrocko_stations(
1926 cls, pyrocko_stations, add_flat_responses_from=None):
1928 '''
1929 Generate :py:class:`FDSNStationXML` from list of
1930 :py:class;`pyrocko.model.Station` instances.
1932 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1933 instances.
1934 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1935 '''
1936 from collections import defaultdict
1937 network_dict = defaultdict(list)
1939 if add_flat_responses_from:
1940 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1941 extra = dict(
1942 response=Response(
1943 instrument_sensitivity=Sensitivity(
1944 value=1.0,
1945 frequency=1.0,
1946 input_units=Units(name=add_flat_responses_from))))
1947 else:
1948 extra = {}
1950 have_offsets = set()
1951 for s in pyrocko_stations:
1953 if s.north_shift != 0.0 or s.east_shift != 0.0:
1954 have_offsets.add(s.nsl())
1956 network, station, location = s.nsl()
1957 channel_list = []
1958 for c in s.channels:
1959 channel_list.append(
1960 Channel(
1961 location_code=location,
1962 code=c.name,
1963 latitude=Latitude(value=s.effective_lat),
1964 longitude=Longitude(value=s.effective_lon),
1965 elevation=Distance(value=s.elevation),
1966 depth=Distance(value=s.depth),
1967 azimuth=Azimuth(value=c.azimuth),
1968 dip=Dip(value=c.dip),
1969 **extra
1970 )
1971 )
1973 network_dict[network].append(
1974 Station(
1975 code=station,
1976 latitude=Latitude(value=s.effective_lat),
1977 longitude=Longitude(value=s.effective_lon),
1978 elevation=Distance(value=s.elevation),
1979 channel_list=channel_list)
1980 )
1982 if have_offsets:
1983 logger.warning(
1984 'StationXML does not support Cartesian offsets in '
1985 'coordinates. Storing effective lat/lon for stations: %s' %
1986 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1988 timestamp = util.to_time_float(time.time())
1989 network_list = []
1990 for k, station_list in network_dict.items():
1992 network_list.append(
1993 Network(
1994 code=k, station_list=station_list,
1995 total_number_stations=len(station_list)))
1997 sxml = FDSNStationXML(
1998 source='from pyrocko stations list', created=timestamp,
1999 network_list=network_list)
2001 sxml.validate()
2002 return sxml
2004 def iter_network_stations(
2005 self, net=None, sta=None, time=None, timespan=None):
2007 tt = ()
2008 if time is not None:
2009 tt = (time,)
2010 elif timespan is not None:
2011 tt = timespan
2013 for network in self.network_list:
2014 if not network.spans(*tt) or (
2015 net is not None and network.code != net):
2016 continue
2018 for station in network.station_list:
2019 if not station.spans(*tt) or (
2020 sta is not None and station.code != sta):
2021 continue
2023 yield (network, station)
2025 def iter_network_station_channels(
2026 self, net=None, sta=None, loc=None, cha=None,
2027 time=None, timespan=None):
2029 if loc is not None:
2030 loc = loc.strip()
2032 tt = ()
2033 if time is not None:
2034 tt = (time,)
2035 elif timespan is not None:
2036 tt = timespan
2038 for network in self.network_list:
2039 if not network.spans(*tt) or (
2040 net is not None and network.code != net):
2041 continue
2043 for station in network.station_list:
2044 if not station.spans(*tt) or (
2045 sta is not None and station.code != sta):
2046 continue
2048 if station.channel_list:
2049 for channel in station.channel_list:
2050 if (not channel.spans(*tt) or
2051 (cha is not None and channel.code != cha) or
2052 (loc is not None and
2053 channel.location_code.strip() != loc)):
2054 continue
2056 yield (network, station, channel)
2058 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2059 time=None, timespan=None):
2061 groups = {}
2062 for network, station, channel in self.iter_network_station_channels(
2063 net, sta, loc, cha, time=time, timespan=timespan):
2065 net = network.code
2066 sta = station.code
2067 cha = channel.code
2068 loc = channel.location_code.strip()
2069 if len(cha) == 3:
2070 bic = cha[:2] # band and intrument code according to SEED
2071 elif len(cha) == 1:
2072 bic = ''
2073 else:
2074 bic = cha
2076 if channel.response and \
2077 channel.response.instrument_sensitivity and \
2078 channel.response.instrument_sensitivity.input_units:
2080 unit = channel.response.instrument_sensitivity\
2081 .input_units.name.upper()
2082 else:
2083 unit = None
2085 bic = (bic, unit)
2087 k = net, sta, loc
2088 if k not in groups:
2089 groups[k] = {}
2091 if bic not in groups[k]:
2092 groups[k][bic] = []
2094 groups[k][bic].append(channel)
2096 for nsl, bic_to_channels in groups.items():
2097 bad_bics = []
2098 for bic, channels in bic_to_channels.items():
2099 sample_rates = []
2100 for channel in channels:
2101 sample_rates.append(channel.sample_rate.value)
2103 if not same(sample_rates):
2104 scs = ','.join(channel.code for channel in channels)
2105 srs = ', '.join('%e' % x for x in sample_rates)
2106 err = 'ignoring channels with inconsistent sampling ' + \
2107 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2109 logger.warning(err)
2110 bad_bics.append(bic)
2112 for bic in bad_bics:
2113 del bic_to_channels[bic]
2115 return groups
2117 def choose_channels(
2118 self,
2119 target_sample_rate=None,
2120 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2121 priority_units=['M/S', 'M/S**2'],
2122 priority_instrument_code=['H', 'L'],
2123 time=None,
2124 timespan=None):
2126 nslcs = {}
2127 for nsl, bic_to_channels in self.get_channel_groups(
2128 time=time, timespan=timespan).items():
2130 useful_bics = []
2131 for bic, channels in bic_to_channels.items():
2132 rate = channels[0].sample_rate.value
2134 if target_sample_rate is not None and \
2135 rate < target_sample_rate*0.99999:
2136 continue
2138 if len(bic[0]) == 2:
2139 if bic[0][0] not in priority_band_code:
2140 continue
2142 if bic[0][1] not in priority_instrument_code:
2143 continue
2145 unit = bic[1]
2147 prio_unit = len(priority_units)
2148 try:
2149 prio_unit = priority_units.index(unit)
2150 except ValueError:
2151 pass
2153 prio_inst = len(priority_instrument_code)
2154 prio_band = len(priority_band_code)
2155 if len(channels[0].code) == 3:
2156 try:
2157 prio_inst = priority_instrument_code.index(
2158 channels[0].code[1])
2159 except ValueError:
2160 pass
2162 try:
2163 prio_band = priority_band_code.index(
2164 channels[0].code[0])
2165 except ValueError:
2166 pass
2168 if target_sample_rate is None:
2169 rate = -rate
2171 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2172 prio_inst, bic))
2174 useful_bics.sort()
2176 for _, _, rate, _, _, bic in useful_bics:
2177 channels = sorted(
2178 bic_to_channels[bic],
2179 key=lambda channel: channel.code)
2181 if channels:
2182 for channel in channels:
2183 nslcs[nsl + (channel.code,)] = channel
2185 break
2187 return nslcs
2189 def get_pyrocko_response(
2190 self, nslc,
2191 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2193 net, sta, loc, cha = nslc
2194 resps = []
2195 for _, _, channel in self.iter_network_station_channels(
2196 net, sta, loc, cha, time=time, timespan=timespan):
2197 resp = channel.response
2198 if resp:
2199 resp.check_sample_rates(channel)
2200 resp.check_units()
2201 resps.append(resp.get_pyrocko_response(
2202 '.'.join(nslc),
2203 fake_input_units=fake_input_units,
2204 stages=stages).expect_one())
2206 if not resps:
2207 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2208 elif len(resps) > 1:
2209 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2211 return resps[0]
2213 @property
2214 def n_code_list(self):
2215 return sorted(set(x.code for x in self.network_list))
2217 @property
2218 def ns_code_list(self):
2219 nss = set()
2220 for network in self.network_list:
2221 for station in network.station_list:
2222 nss.add((network.code, station.code))
2224 return sorted(nss)
2226 @property
2227 def nsl_code_list(self):
2228 nsls = set()
2229 for network in self.network_list:
2230 for station in network.station_list:
2231 for channel in station.channel_list:
2232 nsls.add(
2233 (network.code, station.code, channel.location_code))
2235 return sorted(nsls)
2237 @property
2238 def nslc_code_list(self):
2239 nslcs = set()
2240 for network in self.network_list:
2241 for station in network.station_list:
2242 for channel in station.channel_list:
2243 nslcs.add(
2244 (network.code, station.code, channel.location_code,
2245 channel.code))
2247 return sorted(nslcs)
2249 def summary(self):
2250 lst = [
2251 'number of n codes: %i' % len(self.n_code_list),
2252 'number of ns codes: %i' % len(self.ns_code_list),
2253 'number of nsl codes: %i' % len(self.nsl_code_list),
2254 'number of nslc codes: %i' % len(self.nslc_code_list)
2255 ]
2256 return '\n'.join(lst)
2258 def summary_stages(self):
2259 data = []
2260 for network, station, channel in self.iter_network_station_channels():
2261 nslc = (network.code, station.code, channel.location_code,
2262 channel.code)
2264 stages = []
2265 in_units = '?'
2266 out_units = '?'
2267 if channel.response:
2268 sens = channel.response.instrument_sensitivity
2269 if sens:
2270 in_units = sens.input_units.name.upper()
2271 out_units = sens.output_units.name.upper()
2273 for stage in channel.response.stage_list:
2274 stages.append(stage.summary())
2276 data.append(
2277 (nslc, tts(channel.start_date), tts(channel.end_date),
2278 in_units, out_units, stages))
2280 data.sort()
2282 lst = []
2283 for nslc, stmin, stmax, in_units, out_units, stages in data:
2284 lst.append(' %s: %s - %s, %s -> %s' % (
2285 '.'.join(nslc), stmin, stmax, in_units, out_units))
2286 for stage in stages:
2287 lst.append(' %s' % stage)
2289 return '\n'.join(lst)
2291 def _check_overlaps(self):
2292 by_nslc = {}
2293 for network in self.network_list:
2294 for station in network.station_list:
2295 for channel in station.channel_list:
2296 nslc = (network.code, station.code, channel.location_code,
2297 channel.code)
2298 if nslc not in by_nslc:
2299 by_nslc[nslc] = []
2301 by_nslc[nslc].append(channel)
2303 errors = []
2304 for nslc, channels in by_nslc.items():
2305 errors.extend(check_overlaps('Channel', nslc, channels))
2307 return errors
2309 def check(self):
2310 errors = []
2311 for meth in [self._check_overlaps]:
2312 errors.extend(meth())
2314 if errors:
2315 raise Inconsistencies(
2316 'Inconsistencies found in StationXML:\n '
2317 + '\n '.join(errors))
2320def load_channel_table(stream):
2322 networks = {}
2323 stations = {}
2325 for line in stream:
2326 line = str(line.decode('ascii'))
2327 if line.startswith('#'):
2328 continue
2330 t = line.rstrip().split('|')
2332 if len(t) != 17:
2333 logger.warning('Invalid channel record: %s' % line)
2334 continue
2336 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2337 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2339 try:
2340 scale = float(scale)
2341 except ValueError:
2342 scale = None
2344 try:
2345 scale_freq = float(scale_freq)
2346 except ValueError:
2347 scale_freq = None
2349 try:
2350 depth = float(dep)
2351 except ValueError:
2352 depth = 0.0
2354 try:
2355 azi = float(azi)
2356 dip = float(dip)
2357 except ValueError:
2358 azi = None
2359 dip = None
2361 try:
2362 if net not in networks:
2363 network = Network(code=net)
2364 else:
2365 network = networks[net]
2367 if (net, sta) not in stations:
2368 station = Station(
2369 code=sta, latitude=lat, longitude=lon, elevation=ele)
2371 station.regularize()
2372 else:
2373 station = stations[net, sta]
2375 if scale:
2376 resp = Response(
2377 instrument_sensitivity=Sensitivity(
2378 value=scale,
2379 frequency=scale_freq,
2380 input_units=scale_units))
2381 else:
2382 resp = None
2384 channel = Channel(
2385 code=cha,
2386 location_code=loc.strip(),
2387 latitude=lat,
2388 longitude=lon,
2389 elevation=ele,
2390 depth=depth,
2391 azimuth=azi,
2392 dip=dip,
2393 sensor=Equipment(description=sens),
2394 response=resp,
2395 sample_rate=sample_rate,
2396 start_date=start_date,
2397 end_date=end_date or None)
2399 channel.regularize()
2401 except ValidationError:
2402 raise InvalidRecord(line)
2404 if net not in networks:
2405 networks[net] = network
2407 if (net, sta) not in stations:
2408 stations[net, sta] = station
2409 network.station_list.append(station)
2411 station.channel_list.append(channel)
2413 return FDSNStationXML(
2414 source='created from table input',
2415 created=time.time(),
2416 network_list=sorted(networks.values(), key=lambda x: x.code))
2419def primitive_merge(sxs):
2420 networks = []
2421 for sx in sxs:
2422 networks.extend(sx.network_list)
2424 return FDSNStationXML(
2425 source='merged from different sources',
2426 created=time.time(),
2427 network_list=copy.deepcopy(
2428 sorted(networks, key=lambda x: x.code)))
2431def split_channels(sx):
2432 for nslc in sx.nslc_code_list:
2433 network_list = sx.network_list
2434 network_list_filtered = [
2435 network for network in network_list
2436 if network.code == nslc[0]]
2438 for network in network_list_filtered:
2439 sx.network_list = [network]
2440 station_list = network.station_list
2441 station_list_filtered = [
2442 station for station in station_list
2443 if station.code == nslc[1]]
2445 for station in station_list_filtered:
2446 network.station_list = [station]
2447 channel_list = station.channel_list
2448 station.channel_list = [
2449 channel for channel in channel_list
2450 if (channel.location_code, channel.code) == nslc[2:4]]
2452 yield nslc, copy.deepcopy(sx)
2454 station.channel_list = channel_list
2456 network.station_list = station_list
2458 sx.network_list = network_list
2461if __name__ == '__main__':
2462 from optparse import OptionParser
2464 util.setup_logging('pyrocko.io.stationxml', 'warning')
2466 usage = \
2467 'python -m pyrocko.io.stationxml check|stats|stages ' \
2468 '<filename> [options]'
2470 description = '''Torture StationXML file.'''
2472 parser = OptionParser(
2473 usage=usage,
2474 description=description,
2475 formatter=util.BetterHelpFormatter())
2477 (options, args) = parser.parse_args(sys.argv[1:])
2479 if len(args) != 2:
2480 parser.print_help()
2481 sys.exit(1)
2483 action, path = args
2485 sx = load_xml(filename=path)
2486 if action == 'check':
2487 try:
2488 sx.check()
2489 except Inconsistencies as e:
2490 logger.error(e)
2491 sys.exit(1)
2493 elif action == 'stats':
2494 print(sx.summary())
2496 elif action == 'stages':
2497 print(sx.summary_stages())
2499 else:
2500 parser.print_help()
2501 sys.exit('unknown action: %s' % action)