Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/stationxml.py: 70%
1161 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-06 06:59 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7`FDSN StationXML <https://www.fdsn.org/xml/station/>`_ input, output and data
8model.
9'''
11import sys
12import time
13import logging
14import datetime
15import calendar
16import math
17import copy
19import numpy as num
21from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
22 Unicode, Int, Float, List, Object, Timestamp,
23 ValidationError, TBase, re_tz, Any, Tuple)
24from pyrocko.guts import load_xml # noqa
25from pyrocko.util import hpfloat, time_to_str, get_time_float
27import pyrocko.model
28from pyrocko import util, response
30guts_prefix = 'sx'
32guts_xmlns = 'http://www.fdsn.org/xml/station/1'
34logger = logging.getLogger('pyrocko.io.stationxml')
36conversion = {
37 ('M', 'M'): None,
38 ('M/S', 'M'): response.IntegrationResponse(1),
39 ('M/S**2', 'M'): response.IntegrationResponse(2),
40 ('M', 'M/S'): response.DifferentiationResponse(1),
41 ('M/S', 'M/S'): None,
42 ('M/S**2', 'M/S'): response.IntegrationResponse(1),
43 ('M', 'M/S**2'): response.DifferentiationResponse(2),
44 ('M/S', 'M/S**2'): response.DifferentiationResponse(1),
45 ('M/S**2', 'M/S**2'): None}
48unit_to_quantity = {
49 'M/S': 'velocity',
50 'M': 'displacement',
51 'M/S**2': 'acceleration',
52 'V': 'voltage',
53 'COUNTS': 'counts',
54 'COUNT': 'counts',
55 'PA': 'pressure',
56 'RAD': 'rotation-displacement',
57 'R': 'rotation-displacement',
58 'RAD/S': 'rotation-velocity',
59 'R/S': 'rotation-velocity',
60 'RAD/S**2': 'rotation-acceleration',
61 'R/S**2': 'rotation-acceleration'}
64def to_quantity(unit, context, delivery):
66 if unit is None:
67 return None
69 name = unit.name.upper()
70 if name in unit_to_quantity:
71 return unit_to_quantity[name]
72 else:
73 delivery.log.append((
74 'warning',
75 'Units not supported by Squirrel framework: %s' % (
76 unit.name.upper() + (
77 ' (%s)' % unit.description
78 if unit.description else '')),
79 context))
81 return 'unsupported_quantity(%s)' % unit
84class StationXMLError(Exception):
85 pass
88class Inconsistencies(StationXMLError):
89 pass
92class NoResponseInformation(StationXMLError):
93 pass
96class MultipleResponseInformation(StationXMLError):
97 pass
100class InconsistentResponseInformation(StationXMLError):
101 pass
104class InconsistentChannelLocations(StationXMLError):
105 pass
108class InvalidRecord(StationXMLError):
109 def __init__(self, line):
110 StationXMLError.__init__(self)
111 self._line = line
113 def __str__(self):
114 return 'Invalid record: "%s"' % self._line
117_exceptions = dict(
118 Inconsistencies=Inconsistencies,
119 NoResponseInformation=NoResponseInformation,
120 MultipleResponseInformation=MultipleResponseInformation,
121 InconsistentResponseInformation=InconsistentResponseInformation,
122 InconsistentChannelLocations=InconsistentChannelLocations,
123 InvalidRecord=InvalidRecord,
124 ValueError=ValueError)
127_logs = dict(
128 info=logger.info,
129 warning=logger.warning,
130 error=logger.error)
133class DeliveryError(StationXMLError):
134 pass
137class Delivery(Object):
139 def __init__(self, payload=None, log=None, errors=None, error=None):
140 if payload is None:
141 payload = []
143 if log is None:
144 log = []
146 if errors is None:
147 errors = []
149 if error is not None:
150 errors.append(error)
152 Object.__init__(self, payload=payload, log=log, errors=errors)
154 payload = List.T(Any.T())
155 log = List.T(Tuple.T(3, String.T()))
156 errors = List.T(Tuple.T(3, String.T()))
158 def extend(self, other):
159 self.payload.extend(other.payload)
160 self.log.extend(other.log)
161 self.errors.extend(other.errors)
163 def extend_without_payload(self, other):
164 self.log.extend(other.log)
165 self.errors.extend(other.errors)
166 return other.payload
168 def emit_log(self):
169 for name, message, context in self.log:
170 message = '%s: %s' % (context, message)
171 _logs[name](message)
173 def expect(self, quiet=False):
174 if not quiet:
175 self.emit_log()
177 if self.errors:
178 name, message, context = self.errors[0]
179 if context:
180 message += ' (%s)' % context
182 if len(self.errors) > 1:
183 message += ' Additional errors pending.'
185 raise _exceptions[name](message)
187 return self.payload
189 def expect_one(self, quiet=False):
190 payload = self.expect(quiet=quiet)
191 if len(payload) != 1:
192 raise DeliveryError(
193 'Expected 1 element but got %i.' % len(payload))
195 return payload[0]
198def wrap(s, width=80, indent=4):
199 words = s.split()
200 lines = []
201 t = []
202 n = 0
203 for w in words:
204 if n + len(w) >= width:
205 lines.append(' '.join(t))
206 n = indent
207 t = [' '*(indent-1)]
209 t.append(w)
210 n += len(w) + 1
212 lines.append(' '.join(t))
213 return '\n'.join(lines)
216def same(x, eps=0.0):
217 if any(type(x[0]) != type(r) for r in x):
218 return False
220 if isinstance(x[0], float):
221 return all(abs(r-x[0]) <= eps for r in x)
222 else:
223 return all(r == x[0] for r in x)
226def same_sample_rate(a, b, eps=1.0e-6):
227 return abs(a - b) < min(a, b)*eps
230def evaluate1(resp, f):
231 return resp.evaluate(num.array([f], dtype=float))[0]
234def check_resp(resp, value, frequency, limit_db, prelude, context):
236 try:
237 value_resp = num.abs(evaluate1(resp, frequency))
238 except response.InvalidResponseError as e:
239 return Delivery(log=[(
240 'warning',
241 'Could not check response: %s' % str(e),
242 context)])
244 if value_resp == 0.0:
245 return Delivery(log=[(
246 'warning',
247 '%s\n'
248 ' computed response is zero' % prelude,
249 context)])
251 diff_db = 20.0 * num.log10(value_resp/value)
253 if num.abs(diff_db) > limit_db:
254 return Delivery(log=[(
255 'warning',
256 '%s\n'
257 ' reported value: %g\n'
258 ' computed value: %g\n'
259 ' at frequency [Hz]: %g\n'
260 ' factor: %g\n'
261 ' difference [dB]: %g\n'
262 ' limit [dB]: %g' % (
263 prelude,
264 value,
265 value_resp,
266 frequency,
267 value_resp/value,
268 diff_db,
269 limit_db),
270 context)])
272 return Delivery()
275def tts(t):
276 if t is None:
277 return '<none>'
278 else:
279 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
282def le_open_left(a, b):
283 return a is None or (b is not None and a <= b)
286def le_open_right(a, b):
287 return b is None or (a is not None and a <= b)
290def eq_open(a, b):
291 return (a is None and b is None) \
292 or (a is not None and b is not None and a == b)
295def contains(a, b):
296 return le_open_left(a.start_date, b.start_date) \
297 and le_open_right(b.end_date, a.end_date)
300def find_containing(candidates, node):
301 for candidate in candidates:
302 if contains(candidate, node):
303 return candidate
305 return None
308this_year = time.gmtime()[0]
311class DummyAwareOptionalTimestamp(Object):
312 '''
313 Optional timestamp with support for some common placeholder values.
315 Some StationXML files contain arbitrary placeholder values for open end
316 intervals, like "2100-01-01". Depending on the time range supported by the
317 system, these dates are translated into ``None`` to prevent crashes with
318 this type.
319 '''
320 dummy_for = (hpfloat, float)
321 dummy_for_description = 'pyrocko.util.get_time_float'
323 class __T(TBase):
325 def regularize_extra(self, val):
326 time_float = get_time_float()
328 if isinstance(val, datetime.datetime):
329 tt = val.utctimetuple()
330 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
332 elif isinstance(val, datetime.date):
333 tt = val.timetuple()
334 val = time_float(calendar.timegm(tt))
336 elif isinstance(val, str):
337 val = val.strip()
339 tz_offset = 0
341 m = re_tz.search(val)
342 if m:
343 sh = m.group(2)
344 sm = m.group(4)
345 tz_offset = (int(sh)*3600 if sh else 0) \
346 + (int(sm)*60 if sm else 0)
348 val = re_tz.sub('', val)
350 if len(val) > 10 and val[10] == 'T':
351 val = val.replace('T', ' ', 1)
353 try:
354 val = util.str_to_time(val) - tz_offset
356 except util.TimeStrError:
357 year = int(val[:4])
358 if sys.maxsize > 2**32: # if we're on 64bit
359 if year > this_year + 100:
360 return None # StationXML contained a dummy date
362 if year < 1903: # for macOS, 1900-01-01 dummy dates
363 return None
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']
407class Type(StringChoice):
408 choices = [
409 'TRIGGERED',
410 'CONTINUOUS',
411 'HEALTH',
412 'GEOPHYSICAL',
413 'WEATHER',
414 'FLAG',
415 'SYNTHESIZED',
416 'INPUT',
417 'EXPERIMENTAL',
418 'MAINTENANCE',
419 'BEAM']
421 class __T(StringChoice.T):
422 def validate_extra(self, val):
423 if val not in self.choices:
424 logger.warning(
425 'channel type: "%s" is not a valid choice out of %s' %
426 (val, repr(self.choices)))
429class PzTransferFunction(StringChoice):
430 choices = [
431 'LAPLACE (RADIANS/SECOND)',
432 'LAPLACE (HERTZ)',
433 'DIGITAL (Z-TRANSFORM)']
436class Symmetry(StringChoice):
437 choices = [
438 'NONE',
439 'EVEN',
440 'ODD']
443class CfTransferFunction(StringChoice):
445 class __T(StringChoice.T):
446 def validate(self, val, regularize=False, depth=-1):
447 if regularize:
448 try:
449 val = str(val)
450 except ValueError:
451 raise ValidationError(
452 '%s: cannot convert to string %s' % (self.xname,
453 repr(val)))
455 val = self._dummy_cls.replacements.get(val, val)
457 self.validate_extra(val)
458 return val
460 choices = [
461 'ANALOG (RADIANS/SECOND)',
462 'ANALOG (HERTZ)',
463 'DIGITAL']
465 replacements = {
466 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
467 'ANALOG (HZ)': 'ANALOG (HERTZ)',
468 }
471class Approximation(StringChoice):
472 choices = [
473 'MACLAURIN']
476class PhoneNumber(StringPattern):
477 pattern = '[0-9]+-[0-9]+'
480class Site(Object):
481 '''
482 Description of a site location using name and optional geopolitical
483 boundaries (country, city, etc.).
484 '''
486 name = Unicode.T(default='', xmltagname='Name')
487 description = Unicode.T(optional=True, xmltagname='Description')
488 town = Unicode.T(optional=True, xmltagname='Town')
489 county = Unicode.T(optional=True, xmltagname='County')
490 region = Unicode.T(optional=True, xmltagname='Region')
491 country = Unicode.T(optional=True, xmltagname='Country')
494class ExternalReference(Object):
495 '''
496 This type contains a URI and description for external data that users may
497 want to reference in StationXML.
498 '''
500 uri = String.T(xmltagname='URI')
501 description = Unicode.T(xmltagname='Description')
504class Units(Object):
505 '''
506 A type to document units. Corresponds to SEED blockette 34.
507 '''
509 def __init__(self, name=None, **kwargs):
510 Object.__init__(self, name=name, **kwargs)
512 name = String.T(xmltagname='Name')
513 description = Unicode.T(optional=True, xmltagname='Description')
516class Counter(Int):
517 pass
520class SampleRateRatio(Object):
521 '''
522 Sample rate expressed as number of samples in a number of seconds.
523 '''
525 number_samples = Int.T(xmltagname='NumberSamples')
526 number_seconds = Int.T(xmltagname='NumberSeconds')
529class Gain(Object):
530 '''
531 Complex type for sensitivity and frequency ranges. This complex type can be
532 used to represent both overall sensitivities and individual stage gains.
533 The FrequencyRangeGroup is an optional construct that defines a pass band
534 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
535 is valid within the number of decibels specified in FrequencyDBVariation.
536 '''
538 def __init__(self, value=None, **kwargs):
539 Object.__init__(self, value=value, **kwargs)
541 value = Float.T(optional=True, xmltagname='Value')
542 frequency = Float.T(optional=True, xmltagname='Frequency')
544 def summary(self):
545 return 'gain(%g @ %g)' % (self.value, self.frequency)
548class NumeratorCoefficient(Object):
549 i = Int.T(optional=True, xmlstyle='attribute')
550 value = Float.T(xmlstyle='content')
553class FloatNoUnit(Object):
554 def __init__(self, value=None, **kwargs):
555 Object.__init__(self, value=value, **kwargs)
557 plus_error = Float.T(optional=True, xmlstyle='attribute')
558 minus_error = Float.T(optional=True, xmlstyle='attribute')
559 value = Float.T(xmlstyle='content')
562class FloatWithUnit(FloatNoUnit):
563 unit = String.T(optional=True, xmlstyle='attribute')
566class Equipment(Object):
567 resource_id = String.T(optional=True, xmlstyle='attribute')
568 type = String.T(optional=True, xmltagname='Type')
569 description = Unicode.T(optional=True, xmltagname='Description')
570 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
571 vendor = Unicode.T(optional=True, xmltagname='Vendor')
572 model = Unicode.T(optional=True, xmltagname='Model')
573 serial_number = String.T(optional=True, xmltagname='SerialNumber')
574 installation_date = DummyAwareOptionalTimestamp.T(
575 optional=True,
576 xmltagname='InstallationDate')
577 removal_date = DummyAwareOptionalTimestamp.T(
578 optional=True,
579 xmltagname='RemovalDate')
580 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
583class PhoneNumber(Object):
584 description = Unicode.T(optional=True, xmlstyle='attribute')
585 country_code = Int.T(optional=True, xmltagname='CountryCode')
586 area_code = Int.T(xmltagname='AreaCode')
587 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
590class BaseFilter(Object):
591 '''
592 The BaseFilter is derived by all filters.
593 '''
595 resource_id = String.T(optional=True, xmlstyle='attribute')
596 name = String.T(optional=True, xmlstyle='attribute')
597 description = Unicode.T(optional=True, xmltagname='Description')
598 input_units = Units.T(optional=True, xmltagname='InputUnits')
599 output_units = Units.T(optional=True, xmltagname='OutputUnits')
602class Sensitivity(Gain):
603 '''
604 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
605 construct that defines a pass band in Hertz (FrequencyStart and
606 FrequencyEnd) in which the SensitivityValue is valid within the number of
607 decibels specified in FrequencyDBVariation.
608 '''
610 input_units = Units.T(optional=True, xmltagname='InputUnits')
611 output_units = Units.T(optional=True, xmltagname='OutputUnits')
612 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
613 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
614 frequency_db_variation = Float.T(optional=True,
615 xmltagname='FrequencyDBVariation')
617 def get_pyrocko_response(self):
618 return Delivery(
619 [response.PoleZeroResponse(constant=self.value)])
622class Coefficient(FloatNoUnit):
623 number = Counter.T(optional=True, xmlstyle='attribute')
626class PoleZero(Object):
627 '''
628 Complex numbers used as poles or zeros in channel response.
629 '''
631 number = Int.T(optional=True, xmlstyle='attribute')
632 real = FloatNoUnit.T(xmltagname='Real')
633 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
635 def value(self):
636 return self.real.value + 1J * self.imaginary.value
639class ClockDrift(FloatWithUnit):
640 unit = String.T(default='SECONDS/SAMPLE', optional=True,
641 xmlstyle='attribute') # fixed
644class Second(FloatWithUnit):
645 '''
646 A time value in seconds.
647 '''
649 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
650 # fixed unit
653class Voltage(FloatWithUnit):
654 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
655 # fixed unit
658class Angle(FloatWithUnit):
659 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
660 # fixed unit
663class Azimuth(FloatWithUnit):
664 '''
665 Instrument azimuth, degrees clockwise from North.
666 '''
668 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
669 # fixed unit
672class Dip(FloatWithUnit):
673 '''
674 Instrument dip in degrees down from horizontal. Together azimuth and dip
675 describe the direction of the sensitive axis of the instrument.
676 '''
678 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
679 # fixed unit
682class Distance(FloatWithUnit):
683 '''
684 Extension of FloatWithUnit for distances, elevations, and depths.
685 '''
687 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
688 # NOT fixed unit!
691class Frequency(FloatWithUnit):
692 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
693 # fixed unit
696class SampleRate(FloatWithUnit):
697 '''
698 Sample rate in samples per second.
699 '''
701 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
702 # fixed unit
705class Person(Object):
706 '''
707 Representation of a person's contact information. A person can belong to
708 multiple agencies and have multiple email addresses and phone numbers.
709 '''
711 name_list = List.T(Unicode.T(xmltagname='Name'))
712 agency_list = List.T(Unicode.T(xmltagname='Agency'))
713 email_list = List.T(Email.T(xmltagname='Email'))
714 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
717class FIR(BaseFilter):
718 '''
719 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
720 also commonly documented using the Coefficients element.
721 '''
723 symmetry = Symmetry.T(xmltagname='Symmetry')
724 numerator_coefficient_list = List.T(
725 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
727 def summary(self):
728 return 'fir(%i%s)' % (
729 self.get_ncoefficiencs(),
730 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
732 def get_effective_coefficients(self):
733 b = num.array(
734 [v.value for v in self.numerator_coefficient_list],
735 dtype=float)
737 if self.symmetry == 'ODD':
738 b = num.concatenate((b, b[-2::-1]))
739 elif self.symmetry == 'EVEN':
740 b = num.concatenate((b, b[::-1]))
742 return b
744 def get_effective_symmetry(self):
745 if self.symmetry == 'NONE':
746 b = self.get_effective_coefficients()
747 if num.all(b - b[::-1] == 0):
748 return ['EVEN', 'ODD'][b.size % 2]
750 return self.symmetry
752 def get_ncoefficiencs(self):
753 nf = len(self.numerator_coefficient_list)
754 if self.symmetry == 'ODD':
755 nc = nf * 2 + 1
756 elif self.symmetry == 'EVEN':
757 nc = nf * 2
758 else:
759 nc = nf
761 return nc
763 def estimate_delay(self, deltat):
764 nc = self.get_ncoefficiencs()
765 if nc > 0:
766 return deltat * (nc - 1) / 2.0
767 else:
768 return 0.0
770 def get_pyrocko_response(
771 self, context, deltat, delay_responses, normalization_frequency):
773 context += self.summary()
775 if not self.numerator_coefficient_list:
776 return Delivery([])
778 b = self.get_effective_coefficients()
780 log = []
781 drop_phase = self.get_effective_symmetry() != 'NONE'
783 if not deltat:
784 log.append((
785 'error',
786 'Digital response requires knowledge about sampling '
787 'interval. Response will be unusable.',
788 context))
790 resp = response.DigitalFilterResponse(
791 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
793 if normalization_frequency is not None and deltat is not None:
794 normalization_frequency = 0.0
795 normalization = num.abs(evaluate1(resp, normalization_frequency))
797 if num.abs(normalization - 1.0) > 1e-2:
798 log.append((
799 'warning',
800 'FIR filter coefficients are not normalized. Normalizing '
801 'them. Factor: %g' % normalization, context))
803 resp = response.DigitalFilterResponse(
804 (b/normalization).tolist(), [1.0], deltat,
805 drop_phase=drop_phase)
807 resps = [resp]
809 if not drop_phase:
810 resps.extend(delay_responses)
812 return Delivery(resps, log=log)
815class Coefficients(BaseFilter):
816 '''
817 Response: coefficients for FIR filter. Laplace transforms or IIR filters
818 can be expressed using type as well but the PolesAndZeros should be used
819 instead. Corresponds to SEED blockette 54.
820 '''
822 cf_transfer_function_type = CfTransferFunction.T(
823 xmltagname='CfTransferFunctionType')
824 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
825 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
827 def summary(self):
828 return 'coef_%s(%i,%i%s)' % (
829 'ABC?'[
830 CfTransferFunction.choices.index(
831 self.cf_transfer_function_type)],
832 len(self.numerator_list),
833 len(self.denominator_list),
834 ',sym' if self.is_symmetric_fir else '')
836 def estimate_delay(self, deltat):
837 nc = len(self.numerator_list)
838 if nc > 0:
839 return deltat * (len(self.numerator_list) - 1) / 2.0
840 else:
841 return 0.0
843 def is_symmetric_fir(self):
844 if len(self.denominator_list) != 0:
845 return False
846 b = [v.value for v in self.numerator_list]
847 return b == b[::-1]
849 def get_pyrocko_response(
850 self, context, deltat, delay_responses, normalization_frequency):
852 context += self.summary()
854 factor = 1.0
855 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
856 factor = 2.0*math.pi
858 if not self.numerator_list and not self.denominator_list:
859 return Delivery(payload=[])
861 b = num.array(
862 [v.value*factor for v in self.numerator_list], dtype=float)
864 a = num.array(
865 [1.0] + [v.value*factor for v in self.denominator_list],
866 dtype=float)
868 log = []
869 resps = []
870 if self.cf_transfer_function_type in [
871 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
872 resps.append(response.AnalogFilterResponse(b, a))
874 elif self.cf_transfer_function_type == 'DIGITAL':
875 if not deltat:
876 log.append((
877 'error',
878 'Digital response requires knowledge about sampling '
879 'interval. Response will be unusable.',
880 context))
882 drop_phase = self.is_symmetric_fir()
883 resp = response.DigitalFilterResponse(
884 b, a, deltat or 0.0, drop_phase=drop_phase)
886 if normalization_frequency is not None and deltat is not None:
887 normalization = num.abs(evaluate1(
888 resp, normalization_frequency))
890 if num.abs(normalization - 1.0) > 1e-2:
891 log.append((
892 'warning',
893 'FIR filter coefficients are not normalized. '
894 'Normalizing them. Factor: %g' % normalization,
895 context))
897 resp = response.DigitalFilterResponse(
898 (b/normalization).tolist(), [1.0], deltat,
899 drop_phase=drop_phase)
901 resps.append(resp)
903 if not drop_phase:
904 resps.extend(delay_responses)
906 else:
907 return Delivery(error=(
908 'ValueError',
909 'Unknown transfer function type: %s' % (
910 self.cf_transfer_function_type)))
912 return Delivery(payload=resps, log=log)
915class Latitude(FloatWithUnit):
916 '''
917 Type for latitude coordinate.
918 '''
920 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
921 # fixed unit
922 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
925class Longitude(FloatWithUnit):
926 '''
927 Type for longitude coordinate.
928 '''
930 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
931 # fixed unit
932 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
935class PolesZeros(BaseFilter):
936 '''
937 Response: complex poles and zeros. Corresponds to SEED blockette 53.
938 '''
940 pz_transfer_function_type = PzTransferFunction.T(
941 xmltagname='PzTransferFunctionType')
942 normalization_factor = Float.T(default=1.0,
943 xmltagname='NormalizationFactor')
944 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
945 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
946 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
948 def summary(self):
949 return 'pz_%s(%i,%i)' % (
950 'ABC?'[
951 PzTransferFunction.choices.index(
952 self.pz_transfer_function_type)],
953 len(self.pole_list),
954 len(self.zero_list))
956 def get_pyrocko_response(self, context='', deltat=None):
958 context += self.summary()
960 factor = 1.0
961 cfactor = 1.0
962 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
963 factor = 2. * math.pi
964 cfactor = (2. * math.pi)**(
965 len(self.pole_list) - len(self.zero_list))
967 log = []
968 if self.normalization_factor is None \
969 or self.normalization_factor == 0.0:
971 log.append((
972 'warning',
973 'No pole-zero normalization factor given. '
974 'Assuming a value of 1.0',
975 context))
977 nfactor = 1.0
978 else:
979 nfactor = self.normalization_factor
981 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
982 if not is_digital:
983 resp = response.PoleZeroResponse(
984 constant=nfactor*cfactor,
985 zeros=[z.value()*factor for z in self.zero_list],
986 poles=[p.value()*factor for p in self.pole_list])
987 else:
988 if not deltat:
989 log.append((
990 'error',
991 'Digital response requires knowledge about sampling '
992 'interval. Response will be unusable.',
993 context))
995 resp = response.DigitalPoleZeroResponse(
996 constant=nfactor*cfactor,
997 zeros=[z.value()*factor for z in self.zero_list],
998 poles=[p.value()*factor for p in self.pole_list],
999 deltat=deltat or 0.0)
1001 if not self.normalization_frequency.value:
1002 log.append((
1003 'warning',
1004 'Cannot check pole-zero normalization factor: '
1005 'No normalization frequency given.',
1006 context))
1008 else:
1009 if is_digital and not deltat:
1010 log.append((
1011 'warning',
1012 'Cannot check computed vs reported normalization '
1013 'factor without knowing the sampling interval.',
1014 context))
1015 else:
1016 computed_normalization_factor = nfactor / abs(evaluate1(
1017 resp, self.normalization_frequency.value))
1019 db = 20.0 * num.log10(
1020 computed_normalization_factor / nfactor)
1022 if abs(db) > 0.17:
1023 log.append((
1024 'warning',
1025 'Computed and reported normalization factors differ '
1026 'by %g dB: computed: %g, reported: %g' % (
1027 db,
1028 computed_normalization_factor,
1029 nfactor),
1030 context))
1032 return Delivery([resp], log)
1035class ResponseListElement(Object):
1036 frequency = Frequency.T(xmltagname='Frequency')
1037 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1038 phase = Angle.T(xmltagname='Phase')
1041class Polynomial(BaseFilter):
1042 '''
1043 Response: expressed as a polynomial (allows non-linear sensors to be
1044 described). Corresponds to SEED blockette 62. Can be used to describe a
1045 stage of acquisition or a complete system.
1046 '''
1048 approximation_type = Approximation.T(default='MACLAURIN',
1049 xmltagname='ApproximationType')
1050 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1051 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1052 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1053 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1054 maximum_error = Float.T(xmltagname='MaximumError')
1055 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1057 def summary(self):
1058 return 'poly(%i)' % len(self.coefficient_list)
1061class Decimation(Object):
1062 '''
1063 Corresponds to SEED blockette 57.
1064 '''
1066 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1067 factor = Int.T(xmltagname='Factor')
1068 offset = Int.T(xmltagname='Offset')
1069 delay = FloatWithUnit.T(xmltagname='Delay')
1070 correction = FloatWithUnit.T(xmltagname='Correction')
1072 def summary(self):
1073 return 'deci(%i, %g -> %g, %g)' % (
1074 self.factor,
1075 self.input_sample_rate.value,
1076 self.input_sample_rate.value / self.factor,
1077 self.delay.value)
1079 def get_pyrocko_response(self):
1080 if self.delay and self.delay.value != 0.0:
1081 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1082 else:
1083 return Delivery([])
1086class Operator(Object):
1087 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1088 contact_list = List.T(Person.T(xmltagname='Contact'))
1089 web_site = String.T(optional=True, xmltagname='WebSite')
1092class Comment(Object):
1093 '''
1094 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1095 and 59.
1096 '''
1098 id = Counter.T(optional=True, xmlstyle='attribute')
1099 value = Unicode.T(xmltagname='Value')
1100 begin_effective_time = DummyAwareOptionalTimestamp.T(
1101 optional=True,
1102 xmltagname='BeginEffectiveTime')
1103 end_effective_time = DummyAwareOptionalTimestamp.T(
1104 optional=True,
1105 xmltagname='EndEffectiveTime')
1106 author_list = List.T(Person.T(xmltagname='Author'))
1109class ResponseList(BaseFilter):
1110 '''
1111 Response: list of frequency, amplitude and phase values. Corresponds to
1112 SEED blockette 55.
1113 '''
1115 response_list_element_list = List.T(
1116 ResponseListElement.T(xmltagname='ResponseListElement'))
1118 def summary(self):
1119 return 'list(%i)' % len(self.response_list_element_list)
1122class Log(Object):
1123 '''
1124 Container for log entries.
1125 '''
1127 entry_list = List.T(Comment.T(xmltagname='Entry'))
1130class ResponseStage(Object):
1131 '''
1132 This complex type represents channel response and covers SEED blockettes 53
1133 to 56.
1134 '''
1136 number = Counter.T(xmlstyle='attribute')
1137 resource_id = String.T(optional=True, xmlstyle='attribute')
1138 poles_zeros_list = List.T(
1139 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1140 coefficients_list = List.T(
1141 Coefficients.T(optional=True, xmltagname='Coefficients'))
1142 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1143 fir = FIR.T(optional=True, xmltagname='FIR')
1144 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1145 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1146 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1148 def summary(self):
1149 elements = []
1151 for stuff in [
1152 self.poles_zeros_list, self.coefficients_list,
1153 self.response_list, self.fir, self.polynomial,
1154 self.decimation, self.stage_gain]:
1156 if stuff:
1157 if isinstance(stuff, Object):
1158 elements.append(stuff.summary())
1159 else:
1160 elements.extend(obj.summary() for obj in stuff)
1162 return '%i: %s %s -> %s' % (
1163 self.number,
1164 ', '.join(elements),
1165 self.input_units.name.upper() if self.input_units else '?',
1166 self.output_units.name.upper() if self.output_units else '?')
1168 def get_squirrel_response_stage(self, context):
1169 from pyrocko.squirrel.model import ResponseStage
1171 delivery = Delivery()
1172 delivery_pr = self.get_pyrocko_response(context)
1173 log = delivery_pr.log
1174 delivery_pr.log = []
1175 elements = delivery.extend_without_payload(delivery_pr)
1177 delivery.payload = [ResponseStage(
1178 input_quantity=to_quantity(self.input_units, context, delivery),
1179 output_quantity=to_quantity(self.output_units, context, delivery),
1180 input_sample_rate=self.input_sample_rate,
1181 output_sample_rate=self.output_sample_rate,
1182 elements=elements,
1183 log=log)]
1185 return delivery
1187 def get_pyrocko_response(self, context, gain_only=False):
1189 context = context + ', stage %i' % self.number
1191 responses = []
1192 log = []
1193 if self.stage_gain:
1194 normalization_frequency = self.stage_gain.frequency or 0.0
1195 else:
1196 normalization_frequency = 0.0
1198 if not gain_only:
1199 deltat = None
1200 delay_responses = []
1201 if self.decimation:
1202 rate = self.decimation.input_sample_rate.value
1203 if rate > 0.0:
1204 deltat = 1.0 / rate
1205 delivery = self.decimation.get_pyrocko_response()
1206 if delivery.errors:
1207 return delivery
1209 delay_responses = delivery.payload
1210 log.extend(delivery.log)
1212 for pzs in self.poles_zeros_list:
1213 delivery = pzs.get_pyrocko_response(context, deltat)
1214 if delivery.errors:
1215 return delivery
1217 pz_resps = delivery.payload
1218 log.extend(delivery.log)
1219 responses.extend(pz_resps)
1221 # emulate incorrect? evalresp behaviour
1222 if pzs.normalization_frequency != normalization_frequency \
1223 and normalization_frequency != 0.0:
1225 try:
1226 trial = response.MultiplyResponse(pz_resps)
1227 anorm = num.abs(evaluate1(
1228 trial, pzs.normalization_frequency.value))
1229 asens = num.abs(
1230 evaluate1(trial, normalization_frequency))
1232 factor = anorm/asens
1234 if abs(factor - 1.0) > 0.01:
1235 log.append((
1236 'warning',
1237 'PZ normalization frequency (%g) is different '
1238 'from stage gain frequency (%s) -> Emulating '
1239 'possibly incorrect evalresp behaviour. '
1240 'Correction factor: %g' % (
1241 pzs.normalization_frequency.value,
1242 normalization_frequency,
1243 factor),
1244 context))
1246 responses.append(
1247 response.PoleZeroResponse(constant=factor))
1248 except response.InvalidResponseError as e:
1249 log.append((
1250 'warning',
1251 'Could not check response: %s' % str(e),
1252 context))
1254 if len(self.poles_zeros_list) > 1:
1255 log.append((
1256 'warning',
1257 'Multiple poles and zeros records in single response '
1258 'stage.',
1259 context))
1261 for cfs in self.coefficients_list + (
1262 [self.fir] if self.fir else []):
1264 delivery = cfs.get_pyrocko_response(
1265 context, deltat, delay_responses,
1266 normalization_frequency)
1268 if delivery.errors:
1269 return delivery
1271 responses.extend(delivery.payload)
1272 log.extend(delivery.log)
1274 if len(self.coefficients_list) > 1:
1275 log.append((
1276 'warning',
1277 'Multiple filter coefficients lists in single response '
1278 'stage.',
1279 context))
1281 if self.response_list:
1282 log.append((
1283 'warning',
1284 'Unhandled response element of type: ResponseList',
1285 context))
1287 if self.polynomial:
1288 log.append((
1289 'warning',
1290 'Unhandled response element of type: Polynomial',
1291 context))
1293 if self.stage_gain:
1294 responses.append(
1295 response.PoleZeroResponse(constant=self.stage_gain.value))
1297 return Delivery(responses, log)
1299 @property
1300 def input_units(self):
1301 for e in (self.poles_zeros_list + self.coefficients_list +
1302 [self.response_list, self.fir, self.polynomial]):
1303 if e is not None:
1304 return e.input_units
1306 return None
1308 @property
1309 def output_units(self):
1310 for e in (self.poles_zeros_list + self.coefficients_list +
1311 [self.response_list, self.fir, self.polynomial]):
1312 if e is not None:
1313 return e.output_units
1315 return None
1317 @property
1318 def input_sample_rate(self):
1319 if self.decimation:
1320 return self.decimation.input_sample_rate.value
1322 return None
1324 @property
1325 def output_sample_rate(self):
1326 if self.decimation:
1327 return self.decimation.input_sample_rate.value \
1328 / self.decimation.factor
1330 return None
1333class Response(Object):
1334 resource_id = String.T(optional=True, xmlstyle='attribute')
1335 instrument_sensitivity = Sensitivity.T(optional=True,
1336 xmltagname='InstrumentSensitivity')
1337 instrument_polynomial = Polynomial.T(optional=True,
1338 xmltagname='InstrumentPolynomial')
1339 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1341 @property
1342 def output_sample_rate(self):
1343 if self.stage_list:
1344 return self.stage_list[-1].output_sample_rate
1345 else:
1346 return None
1348 def check_sample_rates(self, channel):
1350 if self.stage_list:
1351 sample_rate = None
1353 for stage in self.stage_list:
1354 if stage.decimation:
1355 input_sample_rate = \
1356 stage.decimation.input_sample_rate.value
1358 if sample_rate is not None and not same_sample_rate(
1359 sample_rate, input_sample_rate):
1361 logger.warning(
1362 'Response stage %i has unexpected input sample '
1363 'rate: %g Hz (expected: %g Hz)' % (
1364 stage.number,
1365 input_sample_rate,
1366 sample_rate))
1368 sample_rate = input_sample_rate / stage.decimation.factor
1370 if sample_rate is not None and channel.sample_rate \
1371 and not same_sample_rate(
1372 sample_rate, channel.sample_rate.value):
1374 logger.warning(
1375 'Channel sample rate (%g Hz) does not match sample rate '
1376 'deducted from response stages information (%g Hz).' % (
1377 channel.sample_rate.value,
1378 sample_rate))
1380 def check_units(self):
1382 if self.instrument_sensitivity \
1383 and self.instrument_sensitivity.input_units:
1385 units = self.instrument_sensitivity.input_units.name.upper()
1387 if self.stage_list:
1388 for stage in self.stage_list:
1389 if units and stage.input_units \
1390 and stage.input_units.name.upper() != units:
1392 logger.warning(
1393 'Input units of stage %i (%s) do not match %s (%s).'
1394 % (
1395 stage.number,
1396 units,
1397 'output units of stage %i'
1398 if stage.number == 0
1399 else 'sensitivity input units',
1400 units))
1402 if stage.output_units:
1403 units = stage.output_units.name.upper()
1404 else:
1405 units = None
1407 sout_units = self.instrument_sensitivity.output_units
1408 if self.instrument_sensitivity and sout_units:
1409 if units is not None and units != sout_units.name.upper():
1410 logger.warning(
1411 'Output units of stage %i (%s) do not match %s (%s).'
1412 % (
1413 stage.number,
1414 units,
1415 'sensitivity output units',
1416 sout_units.name.upper()))
1418 def _sensitivity_checkpoints(self, responses, context):
1419 delivery = Delivery()
1421 if self.instrument_sensitivity:
1422 sval = self.instrument_sensitivity.value
1423 sfreq = self.instrument_sensitivity.frequency
1424 if sval is None:
1425 delivery.log.append((
1426 'warning',
1427 'No sensitivity value given.',
1428 context))
1430 elif sval is None:
1431 delivery.log.append((
1432 'warning',
1433 'Reported sensitivity value is zero.',
1434 context))
1436 elif sfreq is None:
1437 delivery.log.append((
1438 'warning',
1439 'Sensitivity frequency not given.',
1440 context))
1442 else:
1443 trial = response.MultiplyResponse(responses)
1445 delivery.extend(
1446 check_resp(
1447 trial, sval, sfreq, 0.1,
1448 'Instrument sensitivity value inconsistent with '
1449 'sensitivity computed from complete response.',
1450 context))
1452 delivery.payload.append(response.FrequencyResponseCheckpoint(
1453 frequency=sfreq, value=sval))
1455 return delivery
1457 def get_squirrel_response(self, context, **kwargs):
1458 from pyrocko.squirrel.model import Response
1460 if self.stage_list:
1461 delivery = Delivery()
1462 for istage, stage in enumerate(self.stage_list):
1463 delivery.extend(stage.get_squirrel_response_stage(context))
1465 checkpoints = []
1466 if not delivery.errors:
1467 all_responses = []
1468 for sq_stage in delivery.payload:
1469 all_responses.extend(sq_stage.elements)
1471 checkpoints.extend(delivery.extend_without_payload(
1472 self._sensitivity_checkpoints(all_responses, context)))
1474 sq_stages = delivery.payload
1475 if sq_stages:
1476 if sq_stages[0].input_quantity is None \
1477 and self.instrument_sensitivity is not None:
1479 sq_stages[0].input_quantity = to_quantity(
1480 self.instrument_sensitivity.input_units,
1481 context, delivery)
1482 sq_stages[-1].output_quantity = to_quantity(
1483 self.instrument_sensitivity.output_units,
1484 context, delivery)
1486 sq_stages = delivery.expect()
1488 return Response(
1489 stages=sq_stages,
1490 log=delivery.log,
1491 checkpoints=checkpoints,
1492 **kwargs)
1494 elif self.instrument_sensitivity:
1495 raise NoResponseInformation(
1496 "Only instrument sensitivity given (won't use it). (%s)."
1497 % context)
1498 else:
1499 raise NoResponseInformation(
1500 'Empty instrument response (%s).'
1501 % context)
1503 def get_pyrocko_response(
1504 self, context, fake_input_units=None, stages=(0, 1)):
1506 delivery = Delivery()
1507 if self.stage_list:
1508 for istage, stage in enumerate(self.stage_list):
1509 delivery.extend(stage.get_pyrocko_response(
1510 context, gain_only=not (
1511 stages is None or stages[0] <= istage < stages[1])))
1513 elif self.instrument_sensitivity:
1514 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1516 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1517 checkpoints = delivery.extend_without_payload(delivery_cp)
1518 if delivery.errors:
1519 return delivery
1521 if fake_input_units is not None:
1522 if not self.instrument_sensitivity or \
1523 self.instrument_sensitivity.input_units is None:
1525 delivery.errors.append((
1526 'NoResponseInformation',
1527 'No input units given, so cannot convert to requested '
1528 'units: %s' % fake_input_units.upper(),
1529 context))
1531 return delivery
1533 input_units = self.instrument_sensitivity.input_units.name.upper()
1535 conresp = None
1536 try:
1537 conresp = conversion[
1538 fake_input_units.upper(), input_units]
1540 except KeyError:
1541 delivery.errors.append((
1542 'NoResponseInformation',
1543 'Cannot convert between units: %s, %s'
1544 % (fake_input_units, input_units),
1545 context))
1547 if conresp is not None:
1548 delivery.payload.append(conresp)
1549 for checkpoint in checkpoints:
1550 checkpoint.value *= num.abs(evaluate1(
1551 conresp, checkpoint.frequency))
1553 delivery.payload = [
1554 response.MultiplyResponse(
1555 delivery.payload,
1556 checkpoints=checkpoints)]
1558 return delivery
1560 @classmethod
1561 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1562 normalization_frequency=1.0):
1563 '''
1564 Convert Pyrocko pole-zero response to StationXML response.
1566 :param presponse: Pyrocko pole-zero response
1567 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1568 :param input_unit: Input unit to be reported in the StationXML
1569 response.
1570 :type input_unit: str
1571 :param output_unit: Output unit to be reported in the StationXML
1572 response.
1573 :type output_unit: str
1574 :param normalization_frequency: Frequency where the normalization
1575 factor for the StationXML response should be computed.
1576 :type normalization_frequency: float
1577 '''
1579 norm_factor = 1.0/float(abs(
1580 evaluate1(presponse, normalization_frequency)
1581 / presponse.constant))
1583 pzs = PolesZeros(
1584 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1585 normalization_factor=norm_factor,
1586 normalization_frequency=Frequency(normalization_frequency),
1587 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1588 imaginary=FloatNoUnit(z.imag))
1589 for z in presponse.zeros],
1590 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1591 imaginary=FloatNoUnit(z.imag))
1592 for z in presponse.poles])
1594 pzs.validate()
1596 stage = ResponseStage(
1597 number=1,
1598 poles_zeros_list=[pzs],
1599 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1601 resp = Response(
1602 instrument_sensitivity=Sensitivity(
1603 value=stage.stage_gain.value,
1604 frequency=normalization_frequency,
1605 input_units=Units(input_unit),
1606 output_units=Units(output_unit)),
1608 stage_list=[stage])
1610 return resp
1613class BaseNode(Object):
1614 '''
1615 A base node type for derivation from: Network, Station and Channel types.
1616 '''
1618 code = String.T(xmlstyle='attribute')
1619 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1620 xmlstyle='attribute')
1621 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1622 xmlstyle='attribute')
1623 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1624 alternate_code = String.T(optional=True, xmlstyle='attribute')
1625 historical_code = String.T(optional=True, xmlstyle='attribute')
1626 description = Unicode.T(optional=True, xmltagname='Description')
1627 comment_list = List.T(Comment.T(xmltagname='Comment'))
1629 def spans(self, *args):
1630 if len(args) == 0:
1631 return True
1632 elif len(args) == 1:
1633 return ((self.start_date is None or
1634 self.start_date <= args[0]) and
1635 (self.end_date is None or
1636 args[0] <= self.end_date))
1638 elif len(args) == 2:
1639 return ((self.start_date is None or
1640 args[1] >= self.start_date) and
1641 (self.end_date is None or
1642 self.end_date >= args[0]))
1645def overlaps(a, b):
1646 return (
1647 a.start_date is None or b.end_date is None
1648 or a.start_date < b.end_date
1649 ) and (
1650 b.start_date is None or a.end_date is None
1651 or b.start_date < a.end_date
1652 )
1655def check_overlaps(node_type_name, codes, nodes):
1656 errors = []
1657 for ia, a in enumerate(nodes):
1658 for b in nodes[ia+1:]:
1659 if overlaps(a, b):
1660 errors.append(
1661 '%s epochs overlap for %s:\n'
1662 ' %s - %s\n %s - %s' % (
1663 node_type_name,
1664 '.'.join(codes),
1665 tts(a.start_date), tts(a.end_date),
1666 tts(b.start_date), tts(b.end_date)))
1668 return errors
1671class Channel(BaseNode):
1672 '''
1673 Equivalent to SEED blockette 52 and parent element for the related the
1674 response blockettes.
1675 '''
1677 location_code = String.T(xmlstyle='attribute')
1678 external_reference_list = List.T(
1679 ExternalReference.T(xmltagname='ExternalReference'))
1680 latitude = Latitude.T(xmltagname='Latitude')
1681 longitude = Longitude.T(xmltagname='Longitude')
1682 elevation = Distance.T(xmltagname='Elevation')
1683 depth = Distance.T(xmltagname='Depth')
1684 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1685 dip = Dip.T(optional=True, xmltagname='Dip')
1686 type_list = List.T(Type.T(xmltagname='Type'))
1687 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1688 sample_rate_ratio = SampleRateRatio.T(optional=True,
1689 xmltagname='SampleRateRatio')
1690 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1691 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1692 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1693 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1694 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1695 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1696 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1697 response = Response.T(optional=True, xmltagname='Response')
1699 @property
1700 def position_values(self):
1701 lat = self.latitude.value
1702 lon = self.longitude.value
1703 elevation = value_or_none(self.elevation)
1704 depth = value_or_none(self.depth)
1705 return lat, lon, elevation, depth
1708class Station(BaseNode):
1709 '''
1710 This type represents a Station epoch. It is common to only have a single
1711 station epoch with the station's creation and termination dates as the
1712 epoch start and end dates.
1713 '''
1715 latitude = Latitude.T(xmltagname='Latitude')
1716 longitude = Longitude.T(xmltagname='Longitude')
1717 elevation = Distance.T(xmltagname='Elevation')
1718 site = Site.T(default=Site.D(name=''), xmltagname='Site')
1719 vault = Unicode.T(optional=True, xmltagname='Vault')
1720 geology = Unicode.T(optional=True, xmltagname='Geology')
1721 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1722 operator_list = List.T(Operator.T(xmltagname='Operator'))
1723 creation_date = DummyAwareOptionalTimestamp.T(
1724 optional=True, xmltagname='CreationDate')
1725 termination_date = DummyAwareOptionalTimestamp.T(
1726 optional=True, xmltagname='TerminationDate')
1727 total_number_channels = Counter.T(
1728 optional=True, xmltagname='TotalNumberChannels')
1729 selected_number_channels = Counter.T(
1730 optional=True, xmltagname='SelectedNumberChannels')
1731 external_reference_list = List.T(
1732 ExternalReference.T(xmltagname='ExternalReference'))
1733 channel_list = List.T(Channel.T(xmltagname='Channel'))
1735 @property
1736 def position_values(self):
1737 lat = self.latitude.value
1738 lon = self.longitude.value
1739 elevation = value_or_none(self.elevation)
1740 return lat, lon, elevation
1743class Network(BaseNode):
1744 '''
1745 This type represents the Network layer, all station metadata is contained
1746 within this element. The official name of the network or other descriptive
1747 information can be included in the Description element. The Network can
1748 contain 0 or more Stations.
1749 '''
1751 total_number_stations = Counter.T(optional=True,
1752 xmltagname='TotalNumberStations')
1753 selected_number_stations = Counter.T(optional=True,
1754 xmltagname='SelectedNumberStations')
1755 station_list = List.T(Station.T(xmltagname='Station'))
1757 @property
1758 def station_code_list(self):
1759 return sorted(set(s.code for s in self.station_list))
1761 @property
1762 def sl_code_list(self):
1763 sls = set()
1764 for station in self.station_list:
1765 for channel in station.channel_list:
1766 sls.add((station.code, channel.location_code))
1768 return sorted(sls)
1770 def summary(self, width=80, indent=4):
1771 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1772 lines = ['%s (%i):' % (self.code, len(sls))]
1773 if sls:
1774 ssls = ['.'.join(x for x in c if x) for c in sls]
1775 w = max(len(x) for x in ssls)
1776 n = (width - indent) / (w+1)
1777 while ssls:
1778 lines.append(
1779 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1781 ssls[:n] = []
1783 return '\n'.join(lines)
1786def value_or_none(x):
1787 if x is not None:
1788 return x.value
1789 else:
1790 return None
1793def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1795 pos = lat, lon, elevation, depth = \
1796 channels[0].position_values
1798 if not all(pos == x.position_values for x in channels):
1799 info = '\n'.join(
1800 ' %s: %s' % (x.code, x.position_values) for
1801 x in channels)
1803 mess = 'encountered inconsistencies in channel ' \
1804 'lat/lon/elevation/depth ' \
1805 'for %s.%s.%s: \n%s' % (nsl + (info,))
1807 if inconsistencies == 'raise':
1808 raise InconsistentChannelLocations(mess)
1810 elif inconsistencies == 'warn':
1811 logger.warning(mess)
1812 logger.warning(' -> using mean values')
1814 apos = num.array([x.position_values for x in channels], dtype=float)
1815 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1816 / num.sum(num.isfinite(apos), axis=0)
1818 groups = {}
1819 for channel in channels:
1820 if channel.code not in groups:
1821 groups[channel.code] = []
1823 groups[channel.code].append(channel)
1825 pchannels = []
1826 for code in sorted(groups.keys()):
1827 data = [
1828 (channel.code, value_or_none(channel.azimuth),
1829 value_or_none(channel.dip))
1830 for channel in groups[code]]
1832 azimuth, dip = util.consistency_merge(
1833 data,
1834 message='channel orientation values differ:',
1835 error=inconsistencies)
1837 pchannels.append(
1838 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1840 return pyrocko.model.Station(
1841 *nsl,
1842 lat=mlat,
1843 lon=mlon,
1844 elevation=mele,
1845 depth=mdep,
1846 channels=pchannels)
1849class FDSNStationXML(Object):
1850 '''
1851 Top-level type for Station XML. Required field are Source (network ID of
1852 the institution sending the message) and one or more Network containers or
1853 one or more Station containers.
1854 '''
1856 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1857 source = String.T(xmltagname='Source')
1858 sender = String.T(optional=True, xmltagname='Sender')
1859 module = String.T(optional=True, xmltagname='Module')
1860 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1861 created = Timestamp.T(optional=True, xmltagname='Created')
1862 network_list = List.T(Network.T(xmltagname='Network'))
1864 xmltagname = 'FDSNStationXML'
1865 guessable_xmlns = [guts_xmlns]
1867 def __init__(self, *args, **kwargs):
1868 Object.__init__(self, *args, **kwargs)
1869 if self.created is None:
1870 self.created = time.time()
1872 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1873 time=None, timespan=None,
1874 inconsistencies='warn'):
1876 assert inconsistencies in ('raise', 'warn')
1878 if nslcs is not None:
1879 nslcs = set(nslcs)
1881 if nsls is not None:
1882 nsls = set(nsls)
1884 tt = ()
1885 if time is not None:
1886 tt = (time,)
1887 elif timespan is not None:
1888 tt = timespan
1890 pstations = []
1891 for network in self.network_list:
1892 if not network.spans(*tt):
1893 continue
1895 for station in network.station_list:
1896 if not station.spans(*tt):
1897 continue
1899 if station.channel_list:
1900 loc_to_channels = {}
1901 for channel in station.channel_list:
1902 if not channel.spans(*tt):
1903 continue
1905 loc = channel.location_code.strip()
1906 if loc not in loc_to_channels:
1907 loc_to_channels[loc] = []
1909 loc_to_channels[loc].append(channel)
1911 for loc in sorted(loc_to_channels.keys()):
1912 channels = loc_to_channels[loc]
1913 if nslcs is not None:
1914 channels = [channel for channel in channels
1915 if (network.code, station.code, loc,
1916 channel.code) in nslcs]
1918 if not channels:
1919 continue
1921 nsl = network.code, station.code, loc
1922 if nsls is not None and nsl not in nsls:
1923 continue
1925 pstations.append(
1926 pyrocko_station_from_channels(
1927 nsl,
1928 channels,
1929 inconsistencies=inconsistencies))
1930 else:
1931 pstations.append(pyrocko.model.Station(
1932 network.code, station.code, '*',
1933 lat=station.latitude.value,
1934 lon=station.longitude.value,
1935 elevation=value_or_none(station.elevation),
1936 name=station.description or ''))
1938 return pstations
1940 @classmethod
1941 def from_pyrocko_stations(
1942 cls, pyrocko_stations, add_flat_responses_from=None):
1944 '''
1945 Generate :py:class:`FDSNStationXML` from list of
1946 :py:class;`pyrocko.model.Station` instances.
1948 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1949 instances.
1950 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1951 '''
1952 from collections import defaultdict
1953 network_dict = defaultdict(list)
1955 if add_flat_responses_from:
1956 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1957 extra = dict(
1958 response=Response(
1959 instrument_sensitivity=Sensitivity(
1960 value=1.0,
1961 frequency=1.0,
1962 input_units=Units(name=add_flat_responses_from))))
1963 else:
1964 extra = {}
1966 have_offsets = set()
1967 for s in pyrocko_stations:
1969 if s.north_shift != 0.0 or s.east_shift != 0.0:
1970 have_offsets.add(s.nsl())
1972 network, station, location = s.nsl()
1973 channel_list = []
1974 for c in s.channels:
1975 channel_list.append(
1976 Channel(
1977 location_code=location,
1978 code=c.name,
1979 latitude=Latitude(value=s.effective_lat),
1980 longitude=Longitude(value=s.effective_lon),
1981 elevation=Distance(value=s.elevation),
1982 depth=Distance(value=s.depth),
1983 azimuth=Azimuth(value=c.azimuth),
1984 dip=Dip(value=c.dip),
1985 **extra
1986 )
1987 )
1989 network_dict[network].append(
1990 Station(
1991 code=station,
1992 latitude=Latitude(value=s.effective_lat),
1993 longitude=Longitude(value=s.effective_lon),
1994 elevation=Distance(value=s.elevation),
1995 channel_list=channel_list)
1996 )
1998 if have_offsets:
1999 logger.warning(
2000 'StationXML does not support Cartesian offsets in '
2001 'coordinates. Storing effective lat/lon for stations: %s' %
2002 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2004 timestamp = util.to_time_float(time.time())
2005 network_list = []
2006 for k, station_list in network_dict.items():
2008 network_list.append(
2009 Network(
2010 code=k, station_list=station_list,
2011 total_number_stations=len(station_list)))
2013 sxml = FDSNStationXML(
2014 source='from pyrocko stations list', created=timestamp,
2015 network_list=network_list)
2017 sxml.validate()
2018 return sxml
2020 def iter_network_stations(
2021 self, net=None, sta=None, time=None, timespan=None):
2023 tt = ()
2024 if time is not None:
2025 tt = (time,)
2026 elif timespan is not None:
2027 tt = timespan
2029 for network in self.network_list:
2030 if not network.spans(*tt) or (
2031 net is not None and network.code != net):
2032 continue
2034 for station in network.station_list:
2035 if not station.spans(*tt) or (
2036 sta is not None and station.code != sta):
2037 continue
2039 yield (network, station)
2041 def iter_network_station_channels(
2042 self, net=None, sta=None, loc=None, cha=None,
2043 time=None, timespan=None):
2045 if loc is not None:
2046 loc = loc.strip()
2048 tt = ()
2049 if time is not None:
2050 tt = (time,)
2051 elif timespan is not None:
2052 tt = timespan
2054 for network in self.network_list:
2055 if not network.spans(*tt) or (
2056 net is not None and network.code != net):
2057 continue
2059 for station in network.station_list:
2060 if not station.spans(*tt) or (
2061 sta is not None and station.code != sta):
2062 continue
2064 if station.channel_list:
2065 for channel in station.channel_list:
2066 if (not channel.spans(*tt) or
2067 (cha is not None and channel.code != cha) or
2068 (loc is not None and
2069 channel.location_code.strip() != loc)):
2070 continue
2072 yield (network, station, channel)
2074 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2075 time=None, timespan=None):
2077 groups = {}
2078 for network, station, channel in self.iter_network_station_channels(
2079 net, sta, loc, cha, time=time, timespan=timespan):
2081 net = network.code
2082 sta = station.code
2083 cha = channel.code
2084 loc = channel.location_code.strip()
2085 if len(cha) == 3:
2086 bic = cha[:2] # band and intrument code according to SEED
2087 elif len(cha) == 1:
2088 bic = ''
2089 else:
2090 bic = cha
2092 if channel.response and \
2093 channel.response.instrument_sensitivity and \
2094 channel.response.instrument_sensitivity.input_units:
2096 unit = channel.response.instrument_sensitivity\
2097 .input_units.name.upper()
2098 else:
2099 unit = None
2101 bic = (bic, unit)
2103 k = net, sta, loc
2104 if k not in groups:
2105 groups[k] = {}
2107 if bic not in groups[k]:
2108 groups[k][bic] = []
2110 groups[k][bic].append(channel)
2112 for nsl, bic_to_channels in groups.items():
2113 bad_bics = []
2114 for bic, channels in bic_to_channels.items():
2115 sample_rates = []
2116 for channel in channels:
2117 sample_rates.append(channel.sample_rate.value)
2119 if not same(sample_rates):
2120 scs = ','.join(channel.code for channel in channels)
2121 srs = ', '.join('%e' % x for x in sample_rates)
2122 err = 'ignoring channels with inconsistent sampling ' + \
2123 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2125 logger.warning(err)
2126 bad_bics.append(bic)
2128 for bic in bad_bics:
2129 del bic_to_channels[bic]
2131 return groups
2133 def choose_channels(
2134 self,
2135 target_sample_rate=None,
2136 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2137 priority_units=['M/S', 'M/S**2'],
2138 priority_instrument_code=['H', 'L'],
2139 time=None,
2140 timespan=None):
2142 nslcs = {}
2143 for nsl, bic_to_channels in self.get_channel_groups(
2144 time=time, timespan=timespan).items():
2146 useful_bics = []
2147 for bic, channels in bic_to_channels.items():
2148 rate = channels[0].sample_rate.value
2150 if target_sample_rate is not None and \
2151 rate < target_sample_rate*0.99999:
2152 continue
2154 if len(bic[0]) == 2:
2155 if bic[0][0] not in priority_band_code:
2156 continue
2158 if bic[0][1] not in priority_instrument_code:
2159 continue
2161 unit = bic[1]
2163 prio_unit = len(priority_units)
2164 try:
2165 prio_unit = priority_units.index(unit)
2166 except ValueError:
2167 pass
2169 prio_inst = len(priority_instrument_code)
2170 prio_band = len(priority_band_code)
2171 if len(channels[0].code) == 3:
2172 try:
2173 prio_inst = priority_instrument_code.index(
2174 channels[0].code[1])
2175 except ValueError:
2176 pass
2178 try:
2179 prio_band = priority_band_code.index(
2180 channels[0].code[0])
2181 except ValueError:
2182 pass
2184 if target_sample_rate is None:
2185 rate = -rate
2187 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2188 prio_inst, bic))
2190 useful_bics.sort()
2192 for _, _, rate, _, _, bic in useful_bics:
2193 channels = sorted(
2194 bic_to_channels[bic],
2195 key=lambda channel: channel.code)
2197 if channels:
2198 for channel in channels:
2199 nslcs[nsl + (channel.code,)] = channel
2201 break
2203 return nslcs
2205 def get_pyrocko_response(
2206 self, nslc,
2207 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2209 net, sta, loc, cha = nslc
2210 resps = []
2211 for _, _, channel in self.iter_network_station_channels(
2212 net, sta, loc, cha, time=time, timespan=timespan):
2213 resp = channel.response
2214 if resp:
2215 resp.check_sample_rates(channel)
2216 resp.check_units()
2217 resps.append(resp.get_pyrocko_response(
2218 '.'.join(nslc),
2219 fake_input_units=fake_input_units,
2220 stages=stages).expect_one())
2222 if not resps:
2223 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2224 elif len(resps) > 1:
2225 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2227 return resps[0]
2229 @property
2230 def n_code_list(self):
2231 return sorted(set(x.code for x in self.network_list))
2233 @property
2234 def ns_code_list(self):
2235 nss = set()
2236 for network in self.network_list:
2237 for station in network.station_list:
2238 nss.add((network.code, station.code))
2240 return sorted(nss)
2242 @property
2243 def nsl_code_list(self):
2244 nsls = set()
2245 for network in self.network_list:
2246 for station in network.station_list:
2247 for channel in station.channel_list:
2248 nsls.add(
2249 (network.code, station.code, channel.location_code))
2251 return sorted(nsls)
2253 @property
2254 def nslc_code_list(self):
2255 nslcs = set()
2256 for network in self.network_list:
2257 for station in network.station_list:
2258 for channel in station.channel_list:
2259 nslcs.add(
2260 (network.code, station.code, channel.location_code,
2261 channel.code))
2263 return sorted(nslcs)
2265 def summary(self):
2266 lst = [
2267 'number of n codes: %i' % len(self.n_code_list),
2268 'number of ns codes: %i' % len(self.ns_code_list),
2269 'number of nsl codes: %i' % len(self.nsl_code_list),
2270 'number of nslc codes: %i' % len(self.nslc_code_list)
2271 ]
2272 return '\n'.join(lst)
2274 def summary_stages(self):
2275 data = []
2276 for network, station, channel in self.iter_network_station_channels():
2277 nslc = (network.code, station.code, channel.location_code,
2278 channel.code)
2280 stages = []
2281 in_units = '?'
2282 out_units = '?'
2283 if channel.response:
2284 sens = channel.response.instrument_sensitivity
2285 if sens:
2286 in_units = sens.input_units.name.upper()
2287 out_units = sens.output_units.name.upper()
2289 for stage in channel.response.stage_list:
2290 stages.append(stage.summary())
2292 data.append(
2293 (nslc, tts(channel.start_date), tts(channel.end_date),
2294 in_units, out_units, stages))
2296 data.sort()
2298 lst = []
2299 for nslc, stmin, stmax, in_units, out_units, stages in data:
2300 lst.append(' %s: %s - %s, %s -> %s' % (
2301 '.'.join(nslc), stmin, stmax, in_units, out_units))
2302 for stage in stages:
2303 lst.append(' %s' % stage)
2305 return '\n'.join(lst)
2307 def _check_overlaps(self):
2308 by_nslc = {}
2309 for network in self.network_list:
2310 for station in network.station_list:
2311 for channel in station.channel_list:
2312 nslc = (network.code, station.code, channel.location_code,
2313 channel.code)
2314 if nslc not in by_nslc:
2315 by_nslc[nslc] = []
2317 by_nslc[nslc].append(channel)
2319 errors = []
2320 for nslc, channels in by_nslc.items():
2321 errors.extend(check_overlaps('Channel', nslc, channels))
2323 return errors
2325 def check(self):
2326 errors = []
2327 for meth in [self._check_overlaps]:
2328 errors.extend(meth())
2330 if errors:
2331 raise Inconsistencies(
2332 'Inconsistencies found in StationXML:\n '
2333 + '\n '.join(errors))
2336def load_channel_table(stream):
2338 networks = {}
2339 stations = {}
2341 for line in stream:
2342 line = str(line.decode('ascii'))
2343 if line.startswith('#'):
2344 continue
2346 t = line.rstrip().split('|')
2348 if len(t) != 17:
2349 logger.warning('Invalid channel record: %s' % line)
2350 continue
2352 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2353 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2355 try:
2356 scale = float(scale)
2357 except ValueError:
2358 scale = None
2360 try:
2361 scale_freq = float(scale_freq)
2362 except ValueError:
2363 scale_freq = None
2365 try:
2366 depth = float(dep)
2367 except ValueError:
2368 depth = 0.0
2370 try:
2371 azi = float(azi)
2372 dip = float(dip)
2373 except ValueError:
2374 azi = None
2375 dip = None
2377 try:
2378 if net not in networks:
2379 network = Network(code=net)
2380 else:
2381 network = networks[net]
2383 if (net, sta) not in stations:
2384 station = Station(
2385 code=sta, latitude=lat, longitude=lon, elevation=ele)
2387 station.regularize()
2388 else:
2389 station = stations[net, sta]
2391 if scale:
2392 resp = Response(
2393 instrument_sensitivity=Sensitivity(
2394 value=scale,
2395 frequency=scale_freq,
2396 input_units=scale_units))
2397 else:
2398 resp = None
2400 channel = Channel(
2401 code=cha,
2402 location_code=loc.strip(),
2403 latitude=lat,
2404 longitude=lon,
2405 elevation=ele,
2406 depth=depth,
2407 azimuth=azi,
2408 dip=dip,
2409 sensor=Equipment(description=sens),
2410 response=resp,
2411 sample_rate=sample_rate,
2412 start_date=start_date,
2413 end_date=end_date or None)
2415 channel.regularize()
2417 except ValidationError:
2418 raise InvalidRecord(line)
2420 if net not in networks:
2421 networks[net] = network
2423 if (net, sta) not in stations:
2424 stations[net, sta] = station
2425 network.station_list.append(station)
2427 station.channel_list.append(channel)
2429 return FDSNStationXML(
2430 source='created from table input',
2431 created=time.time(),
2432 network_list=sorted(networks.values(), key=lambda x: x.code))
2435def primitive_merge(sxs):
2436 networks = []
2437 for sx in sxs:
2438 networks.extend(sx.network_list)
2440 return FDSNStationXML(
2441 source='merged from different sources',
2442 created=time.time(),
2443 network_list=copy.deepcopy(
2444 sorted(networks, key=lambda x: x.code)))
2447def split_channels(sx):
2448 for nslc in sx.nslc_code_list:
2449 network_list = sx.network_list
2450 network_list_filtered = [
2451 network for network in network_list
2452 if network.code == nslc[0]]
2454 for network in network_list_filtered:
2455 sx.network_list = [network]
2456 station_list = network.station_list
2457 station_list_filtered = [
2458 station for station in station_list
2459 if station.code == nslc[1]]
2461 for station in station_list_filtered:
2462 network.station_list = [station]
2463 channel_list = station.channel_list
2464 station.channel_list = [
2465 channel for channel in channel_list
2466 if (channel.location_code, channel.code) == nslc[2:4]]
2468 yield nslc, copy.deepcopy(sx)
2470 station.channel_list = channel_list
2472 network.station_list = station_list
2474 sx.network_list = network_list
2477if __name__ == '__main__':
2478 from optparse import OptionParser
2480 util.setup_logging('pyrocko.io.stationxml', 'warning')
2482 usage = \
2483 'python -m pyrocko.io.stationxml check|stats|stages ' \
2484 '<filename> [options]'
2486 description = '''Torture StationXML file.'''
2488 parser = OptionParser(
2489 usage=usage,
2490 description=description,
2491 formatter=util.BetterHelpFormatter())
2493 (options, args) = parser.parse_args(sys.argv[1:])
2495 if len(args) != 2:
2496 parser.print_help()
2497 sys.exit(1)
2499 action, path = args
2501 sx = load_xml(filename=path)
2502 if action == 'check':
2503 try:
2504 sx.check()
2505 except Inconsistencies as e:
2506 logger.error(e)
2507 sys.exit(1)
2509 elif action == 'stats':
2510 print(sx.summary())
2512 elif action == 'stages':
2513 print(sx.summary_stages())
2515 else:
2516 parser.print_help()
2517 sys.exit('unknown action: %s' % action)