Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/stationxml.py: 70%
1177 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-12 12:00 +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
18from collections import defaultdict
20import numpy as num
22from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
23 Unicode, Int, Float, List, Object, Timestamp,
24 ValidationError, TBase, re_tz, Any, Tuple)
25from pyrocko.guts import load_xml # noqa
26from pyrocko.util import hpfloat, time_to_str, get_time_float
28import pyrocko.model
29from pyrocko import util, response
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': 'rotation-displacement',
58 'R': 'rotation-displacement',
59 'RAD/S': 'rotation-velocity',
60 'R/S': 'rotation-velocity',
61 'RAD/S**2': 'rotation-acceleration',
62 'R/S**2': 'rotation-acceleration'}
65def to_quantity(unit, context, delivery):
67 if unit is None:
68 return None
70 name = unit.name.upper()
71 if name in unit_to_quantity:
72 return unit_to_quantity[name]
73 else:
74 delivery.log.append((
75 'warning',
76 'Units not supported by Squirrel framework: %s' % (
77 unit.name.upper() + (
78 ' (%s)' % unit.description
79 if unit.description else '')),
80 context))
82 return 'unsupported_quantity(%s)' % unit
85class StationXMLError(Exception):
86 pass
89class Inconsistencies(StationXMLError):
90 pass
93class NoResponseInformation(StationXMLError):
94 pass
97class MultipleResponseInformation(StationXMLError):
98 pass
101class InconsistentResponseInformation(StationXMLError):
102 pass
105class InconsistentChannelLocations(StationXMLError):
106 pass
109class InvalidRecord(StationXMLError):
110 def __init__(self, line):
111 StationXMLError.__init__(self)
112 self._line = line
114 def __str__(self):
115 return 'Invalid record: "%s"' % self._line
118_exceptions = dict(
119 Inconsistencies=Inconsistencies,
120 NoResponseInformation=NoResponseInformation,
121 MultipleResponseInformation=MultipleResponseInformation,
122 InconsistentResponseInformation=InconsistentResponseInformation,
123 InconsistentChannelLocations=InconsistentChannelLocations,
124 InvalidRecord=InvalidRecord,
125 ValueError=ValueError)
128_logs = dict(
129 info=logger.info,
130 warning=logger.warning,
131 error=logger.error)
134class DeliveryError(StationXMLError):
135 pass
138class Delivery(Object):
140 def __init__(self, payload=None, log=None, errors=None, error=None):
141 if payload is None:
142 payload = []
144 if log is None:
145 log = []
147 if errors is None:
148 errors = []
150 if error is not None:
151 errors.append(error)
153 Object.__init__(self, payload=payload, log=log, errors=errors)
155 payload = List.T(Any.T())
156 log = List.T(Tuple.T(3, String.T()))
157 errors = List.T(Tuple.T(3, String.T()))
159 def extend(self, other):
160 self.payload.extend(other.payload)
161 self.log.extend(other.log)
162 self.errors.extend(other.errors)
164 def extend_without_payload(self, other):
165 self.log.extend(other.log)
166 self.errors.extend(other.errors)
167 return other.payload
169 def emit_log(self):
170 for name, message, context in self.log:
171 message = '%s: %s' % (context, message)
172 _logs[name](message)
174 def expect(self, quiet=False):
175 if not quiet:
176 self.emit_log()
178 if self.errors:
179 name, message, context = self.errors[0]
180 if context:
181 message += ' (%s)' % context
183 if len(self.errors) > 1:
184 message += ' Additional errors pending.'
186 raise _exceptions[name](message)
188 return self.payload
190 def expect_one(self, quiet=False):
191 payload = self.expect(quiet=quiet)
192 if len(payload) != 1:
193 raise DeliveryError(
194 'Expected 1 element but got %i.' % len(payload))
196 return payload[0]
199def wrap(s, width=80, indent=4):
200 words = s.split()
201 lines = []
202 t = []
203 n = 0
204 for w in words:
205 if n + len(w) >= width:
206 lines.append(' '.join(t))
207 n = indent
208 t = [' '*(indent-1)]
210 t.append(w)
211 n += len(w) + 1
213 lines.append(' '.join(t))
214 return '\n'.join(lines)
217def same(x, eps=0.0):
218 if any(type(x[0]) != type(r) for r in x):
219 return False
221 if isinstance(x[0], float):
222 return all(abs(r-x[0]) <= eps for r in x)
223 else:
224 return all(r == x[0] for r in x)
227def same_sample_rate(a, b, eps=1.0e-6):
228 return abs(a - b) < min(a, b)*eps
231def evaluate1(resp, f):
232 return resp.evaluate(num.array([f], dtype=float))[0]
235def check_resp(resp, value, frequency, limit_db, prelude, context):
237 try:
238 value_resp = num.abs(evaluate1(resp, frequency))
239 except response.InvalidResponseError as e:
240 return Delivery(log=[(
241 'warning',
242 'Could not check response: %s' % str(e),
243 context)])
245 if value_resp == 0.0:
246 return Delivery(log=[(
247 'warning',
248 '%s\n'
249 ' computed response is zero' % prelude,
250 context)])
252 diff_db = 20.0 * num.log10(value_resp/value)
254 if num.abs(diff_db) > limit_db:
255 return Delivery(log=[(
256 'warning',
257 '%s\n'
258 ' reported value: %g\n'
259 ' computed value: %g\n'
260 ' at frequency [Hz]: %g\n'
261 ' factor: %g\n'
262 ' difference [dB]: %g\n'
263 ' limit [dB]: %g' % (
264 prelude,
265 value,
266 value_resp,
267 frequency,
268 value_resp/value,
269 diff_db,
270 limit_db),
271 context)])
273 return Delivery()
276def tts(t):
277 if t is None:
278 return '<none>'
279 else:
280 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
283def le_open_left(a, b):
284 return a is None or (b is not None and a <= b)
287def le_open_right(a, b):
288 return b is None or (a is not None and a <= b)
291def eq_open(a, b):
292 return (a is None and b is None) \
293 or (a is not None and b is not None and a == b)
296def contains(a, b):
297 return le_open_left(a.start_date, b.start_date) \
298 and le_open_right(b.end_date, a.end_date)
301def find_containing(candidates, node):
302 for candidate in candidates:
303 if contains(candidate, node):
304 return candidate
306 return None
309this_year = time.gmtime()[0]
312class DummyAwareOptionalTimestamp(Object):
313 '''
314 Optional timestamp with support for some common placeholder values.
316 Some StationXML files contain arbitrary placeholder values for open end
317 intervals, like "2100-01-01". Depending on the time range supported by the
318 system, these dates are translated into ``None`` to prevent crashes with
319 this type.
320 '''
321 dummy_for = (hpfloat, float)
322 dummy_for_description = 'pyrocko.util.get_time_float'
324 class __T(TBase):
326 def regularize_extra(self, val):
327 time_float = get_time_float()
329 if isinstance(val, datetime.datetime):
330 tt = val.utctimetuple()
331 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
333 elif isinstance(val, datetime.date):
334 tt = val.timetuple()
335 val = time_float(calendar.timegm(tt))
337 elif isinstance(val, str):
338 val = val.strip()
340 tz_offset = 0
342 m = re_tz.search(val)
343 if m:
344 sh = m.group(2)
345 sm = m.group(4)
346 tz_offset = (int(sh)*3600 if sh else 0) \
347 + (int(sm)*60 if sm else 0)
349 val = re_tz.sub('', val)
351 if len(val) > 10 and val[10] == 'T':
352 val = val.replace('T', ' ', 1)
354 try:
355 val = util.str_to_time(val) - tz_offset
357 except util.TimeStrError:
358 year = int(val[:4])
359 if sys.maxsize > 2**32: # if we're on 64bit
360 if year > this_year + 100:
361 return None # StationXML contained a dummy date
363 if year < 1903: # for macOS, 1900-01-01 dummy dates
364 return None
366 else: # 32bit end of time is in 2038
367 if this_year < 2037 and year > 2037 or year < 1903:
368 return None # StationXML contained a dummy date
370 raise
372 elif isinstance(val, (int, float)):
373 val = time_float(val)
375 else:
376 raise ValidationError(
377 '%s: cannot convert "%s" to type %s' % (
378 self.xname(), val, time_float))
380 return val
382 def to_save(self, val):
383 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
384 .rstrip('0').rstrip('.')
386 def to_save_xml(self, val):
387 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
388 .rstrip('0').rstrip('.') + 'Z'
391class Nominal(StringChoice):
392 choices = [
393 'NOMINAL',
394 'CALCULATED']
397class Email(UnicodePattern):
398 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
401class RestrictedStatus(StringChoice):
402 choices = [
403 'open',
404 'closed',
405 'partial']
408class Type(StringChoice):
409 choices = [
410 'TRIGGERED',
411 'CONTINUOUS',
412 'HEALTH',
413 'GEOPHYSICAL',
414 'WEATHER',
415 'FLAG',
416 'SYNTHESIZED',
417 'INPUT',
418 'EXPERIMENTAL',
419 'MAINTENANCE',
420 'BEAM']
422 class __T(StringChoice.T):
423 def validate_extra(self, val):
424 if val not in self.choices:
425 logger.warning(
426 'channel type: "%s" is not a valid choice out of %s' %
427 (val, repr(self.choices)))
430class PzTransferFunction(StringChoice):
431 choices = [
432 'LAPLACE (RADIANS/SECOND)',
433 'LAPLACE (HERTZ)',
434 'DIGITAL (Z-TRANSFORM)']
437class Symmetry(StringChoice):
438 choices = [
439 'NONE',
440 'EVEN',
441 'ODD']
444class CfTransferFunction(StringChoice):
446 class __T(StringChoice.T):
447 def validate(self, val, regularize=False, depth=-1):
448 if regularize:
449 try:
450 val = str(val)
451 except ValueError:
452 raise ValidationError(
453 '%s: cannot convert to string %s' % (self.xname,
454 repr(val)))
456 val = self._dummy_cls.replacements.get(val, val)
458 self.validate_extra(val)
459 return val
461 choices = [
462 'ANALOG (RADIANS/SECOND)',
463 'ANALOG (HERTZ)',
464 'DIGITAL']
466 replacements = {
467 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
468 'ANALOG (HZ)': 'ANALOG (HERTZ)',
469 }
472class Approximation(StringChoice):
473 choices = [
474 'MACLAURIN']
477class PhoneNumber(StringPattern):
478 pattern = '[0-9]+-[0-9]+'
481class Site(Object):
482 '''
483 Description of a site location using name and optional geopolitical
484 boundaries (country, city, etc.).
485 '''
487 name = Unicode.T(default='', xmltagname='Name')
488 description = Unicode.T(optional=True, xmltagname='Description')
489 town = Unicode.T(optional=True, xmltagname='Town')
490 county = Unicode.T(optional=True, xmltagname='County')
491 region = Unicode.T(optional=True, xmltagname='Region')
492 country = Unicode.T(optional=True, xmltagname='Country')
495class ExternalReference(Object):
496 '''
497 This type contains a URI and description for external data that users may
498 want to reference in StationXML.
499 '''
501 uri = String.T(xmltagname='URI')
502 description = Unicode.T(xmltagname='Description')
505class Units(Object):
506 '''
507 A type to document units. Corresponds to SEED blockette 34.
508 '''
510 def __init__(self, name=None, **kwargs):
511 Object.__init__(self, name=name, **kwargs)
513 name = String.T(xmltagname='Name')
514 description = Unicode.T(optional=True, xmltagname='Description')
517class Counter(Int):
518 pass
521class SampleRateRatio(Object):
522 '''
523 Sample rate expressed as number of samples in a number of seconds.
524 '''
526 number_samples = Int.T(xmltagname='NumberSamples')
527 number_seconds = Int.T(xmltagname='NumberSeconds')
530class Gain(Object):
531 '''
532 Complex type for sensitivity and frequency ranges. This complex type can be
533 used to represent both overall sensitivities and individual stage gains.
534 The FrequencyRangeGroup is an optional construct that defines a pass band
535 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
536 is valid within the number of decibels specified in FrequencyDBVariation.
537 '''
539 def __init__(self, value=None, **kwargs):
540 Object.__init__(self, value=value, **kwargs)
542 value = Float.T(optional=True, xmltagname='Value')
543 frequency = Float.T(optional=True, xmltagname='Frequency')
545 def summary(self):
546 return 'gain(%g @ %g)' % (self.value, self.frequency)
549class NumeratorCoefficient(Object):
550 i = Int.T(optional=True, xmlstyle='attribute')
551 value = Float.T(xmlstyle='content')
554class FloatNoUnit(Object):
555 def __init__(self, value=None, **kwargs):
556 Object.__init__(self, value=value, **kwargs)
558 plus_error = Float.T(optional=True, xmlstyle='attribute')
559 minus_error = Float.T(optional=True, xmlstyle='attribute')
560 value = Float.T(xmlstyle='content')
563class FloatWithUnit(FloatNoUnit):
564 unit = String.T(optional=True, xmlstyle='attribute')
567class Equipment(Object):
568 resource_id = String.T(optional=True, xmlstyle='attribute')
569 type = String.T(optional=True, xmltagname='Type')
570 description = Unicode.T(optional=True, xmltagname='Description')
571 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
572 vendor = Unicode.T(optional=True, xmltagname='Vendor')
573 model = Unicode.T(optional=True, xmltagname='Model')
574 serial_number = String.T(optional=True, xmltagname='SerialNumber')
575 installation_date = DummyAwareOptionalTimestamp.T(
576 optional=True,
577 xmltagname='InstallationDate')
578 removal_date = DummyAwareOptionalTimestamp.T(
579 optional=True,
580 xmltagname='RemovalDate')
581 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
584class PhoneNumber(Object):
585 description = Unicode.T(optional=True, xmlstyle='attribute')
586 country_code = Int.T(optional=True, xmltagname='CountryCode')
587 area_code = Int.T(xmltagname='AreaCode')
588 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
591class BaseFilter(Object):
592 '''
593 The BaseFilter is derived by all filters.
594 '''
596 resource_id = String.T(optional=True, xmlstyle='attribute')
597 name = String.T(optional=True, xmlstyle='attribute')
598 description = Unicode.T(optional=True, xmltagname='Description')
599 input_units = Units.T(optional=True, xmltagname='InputUnits')
600 output_units = Units.T(optional=True, xmltagname='OutputUnits')
603class Sensitivity(Gain):
604 '''
605 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
606 construct that defines a pass band in Hertz (FrequencyStart and
607 FrequencyEnd) in which the SensitivityValue is valid within the number of
608 decibels specified in FrequencyDBVariation.
609 '''
611 input_units = Units.T(optional=True, xmltagname='InputUnits')
612 output_units = Units.T(optional=True, xmltagname='OutputUnits')
613 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
614 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
615 frequency_db_variation = Float.T(optional=True,
616 xmltagname='FrequencyDBVariation')
618 def get_pyrocko_response(self):
619 return Delivery(
620 [response.PoleZeroResponse(constant=self.value)])
623class Coefficient(FloatNoUnit):
624 number = Counter.T(optional=True, xmlstyle='attribute')
627class PoleZero(Object):
628 '''
629 Complex numbers used as poles or zeros in channel response.
630 '''
632 number = Int.T(optional=True, xmlstyle='attribute')
633 real = FloatNoUnit.T(xmltagname='Real')
634 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
636 def value(self):
637 return self.real.value + 1J * self.imaginary.value
640class ClockDrift(FloatWithUnit):
641 unit = String.T(default='SECONDS/SAMPLE', optional=True,
642 xmlstyle='attribute') # fixed
645class Second(FloatWithUnit):
646 '''
647 A time value in seconds.
648 '''
650 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
651 # fixed unit
654class Voltage(FloatWithUnit):
655 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
656 # fixed unit
659class Angle(FloatWithUnit):
660 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
661 # fixed unit
664class Azimuth(FloatWithUnit):
665 '''
666 Instrument azimuth, degrees clockwise from North.
667 '''
669 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
670 # fixed unit
673class Dip(FloatWithUnit):
674 '''
675 Instrument dip in degrees down from horizontal. Together azimuth and dip
676 describe the direction of the sensitive axis of the instrument.
677 '''
679 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
680 # fixed unit
683class Distance(FloatWithUnit):
684 '''
685 Extension of FloatWithUnit for distances, elevations, and depths.
686 '''
688 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
689 # NOT fixed unit!
692class Frequency(FloatWithUnit):
693 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
694 # fixed unit
697class SampleRate(FloatWithUnit):
698 '''
699 Sample rate in samples per second.
700 '''
702 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
703 # fixed unit
706class Person(Object):
707 '''
708 Representation of a person's contact information. A person can belong to
709 multiple agencies and have multiple email addresses and phone numbers.
710 '''
712 name_list = List.T(Unicode.T(xmltagname='Name'))
713 agency_list = List.T(Unicode.T(xmltagname='Agency'))
714 email_list = List.T(Email.T(xmltagname='Email'))
715 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
718class FIR(BaseFilter):
719 '''
720 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
721 also commonly documented using the Coefficients element.
722 '''
724 symmetry = Symmetry.T(xmltagname='Symmetry')
725 numerator_coefficient_list = List.T(
726 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
728 def summary(self):
729 return 'fir(%i%s)' % (
730 self.get_ncoefficiencs(),
731 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
733 def get_effective_coefficients(self):
734 b = num.array(
735 [v.value for v in self.numerator_coefficient_list],
736 dtype=float)
738 if self.symmetry == 'ODD':
739 b = num.concatenate((b, b[-2::-1]))
740 elif self.symmetry == 'EVEN':
741 b = num.concatenate((b, b[::-1]))
743 return b
745 def get_effective_symmetry(self):
746 if self.symmetry == 'NONE':
747 b = self.get_effective_coefficients()
748 if num.all(b - b[::-1] == 0):
749 return ['EVEN', 'ODD'][b.size % 2]
751 return self.symmetry
753 def get_ncoefficiencs(self):
754 nf = len(self.numerator_coefficient_list)
755 if self.symmetry == 'ODD':
756 nc = nf * 2 + 1
757 elif self.symmetry == 'EVEN':
758 nc = nf * 2
759 else:
760 nc = nf
762 return nc
764 def estimate_delay(self, deltat):
765 nc = self.get_ncoefficiencs()
766 if nc > 0:
767 return deltat * (nc - 1) / 2.0
768 else:
769 return 0.0
771 def get_pyrocko_response(
772 self, context, deltat, delay_responses, normalization_frequency):
774 context += self.summary()
776 if not self.numerator_coefficient_list:
777 return Delivery([])
779 b = self.get_effective_coefficients()
781 log = []
782 drop_phase = self.get_effective_symmetry() != 'NONE'
784 if not deltat:
785 log.append((
786 'error',
787 'Digital response requires knowledge about sampling '
788 'interval. Response will be unusable.',
789 context))
791 resp = response.DigitalFilterResponse(
792 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
794 if normalization_frequency is not None and deltat is not None:
795 normalization_frequency = 0.0
796 normalization = num.abs(evaluate1(resp, normalization_frequency))
798 if num.abs(normalization - 1.0) > 1e-2:
799 log.append((
800 'warning',
801 'FIR filter coefficients are not normalized. Normalizing '
802 'them. Factor: %g' % normalization, context))
804 resp = response.DigitalFilterResponse(
805 (b/normalization).tolist(), [1.0], deltat,
806 drop_phase=drop_phase)
808 resps = [resp]
810 if not drop_phase:
811 resps.extend(delay_responses)
813 return Delivery(resps, log=log)
816class Coefficients(BaseFilter):
817 '''
818 Response: coefficients for FIR filter. Laplace transforms or IIR filters
819 can be expressed using type as well but the PolesAndZeros should be used
820 instead. Corresponds to SEED blockette 54.
821 '''
823 cf_transfer_function_type = CfTransferFunction.T(
824 xmltagname='CfTransferFunctionType')
825 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
826 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
828 def summary(self):
829 return 'coef_%s(%i,%i%s)' % (
830 'ABC?'[
831 CfTransferFunction.choices.index(
832 self.cf_transfer_function_type)],
833 len(self.numerator_list),
834 len(self.denominator_list),
835 ',sym' if self.is_symmetric_fir else '')
837 def estimate_delay(self, deltat):
838 nc = len(self.numerator_list)
839 if nc > 0:
840 return deltat * (len(self.numerator_list) - 1) / 2.0
841 else:
842 return 0.0
844 def is_symmetric_fir(self):
845 if len(self.denominator_list) != 0:
846 return False
847 b = [v.value for v in self.numerator_list]
848 return b == b[::-1]
850 def get_pyrocko_response(
851 self, context, deltat, delay_responses, normalization_frequency):
853 context += self.summary()
855 factor = 1.0
856 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
857 factor = 2.0*math.pi
859 if not self.numerator_list and not self.denominator_list:
860 return Delivery(payload=[])
862 b = num.array(
863 [v.value*factor for v in self.numerator_list], dtype=float)
865 a = num.array(
866 [1.0] + [v.value*factor for v in self.denominator_list],
867 dtype=float)
869 log = []
870 resps = []
871 if self.cf_transfer_function_type in [
872 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
873 resps.append(response.AnalogFilterResponse(b, a))
875 elif self.cf_transfer_function_type == 'DIGITAL':
876 if not deltat:
877 log.append((
878 'error',
879 'Digital response requires knowledge about sampling '
880 'interval. Response will be unusable.',
881 context))
883 drop_phase = self.is_symmetric_fir()
884 resp = response.DigitalFilterResponse(
885 b, a, deltat or 0.0, drop_phase=drop_phase)
887 if normalization_frequency is not None and deltat is not None:
888 normalization = num.abs(evaluate1(
889 resp, normalization_frequency))
891 if num.abs(normalization - 1.0) > 1e-2:
892 log.append((
893 'warning',
894 'FIR filter coefficients are not normalized. '
895 'Normalizing them. Factor: %g' % normalization,
896 context))
898 resp = response.DigitalFilterResponse(
899 (b/normalization).tolist(), [1.0], deltat,
900 drop_phase=drop_phase)
902 resps.append(resp)
904 if not drop_phase:
905 resps.extend(delay_responses)
907 else:
908 return Delivery(error=(
909 'ValueError',
910 'Unknown transfer function type: %s' % (
911 self.cf_transfer_function_type)))
913 return Delivery(payload=resps, log=log)
916class Latitude(FloatWithUnit):
917 '''
918 Type for latitude coordinate.
919 '''
921 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
922 # fixed unit
923 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
926class Longitude(FloatWithUnit):
927 '''
928 Type for longitude coordinate.
929 '''
931 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
932 # fixed unit
933 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
936class PolesZeros(BaseFilter):
937 '''
938 Response: complex poles and zeros. Corresponds to SEED blockette 53.
939 '''
941 pz_transfer_function_type = PzTransferFunction.T(
942 xmltagname='PzTransferFunctionType')
943 normalization_factor = Float.T(
944 default=1.0,
945 xmltagname='NormalizationFactor')
946 normalization_frequency = Frequency.T(
947 optional=True, # but required by standard
948 xmltagname='NormalizationFrequency')
949 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
950 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
952 def summary(self):
953 return 'pz_%s(%i,%i)' % (
954 'ABC?'[
955 PzTransferFunction.choices.index(
956 self.pz_transfer_function_type)],
957 len(self.pole_list),
958 len(self.zero_list))
960 def get_pyrocko_response(self, context='', deltat=None):
962 context += self.summary()
964 factor = 1.0
965 cfactor = 1.0
966 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
967 factor = 2. * math.pi
968 cfactor = (2. * math.pi)**(
969 len(self.pole_list) - len(self.zero_list))
971 log = []
972 if self.normalization_factor is None \
973 or self.normalization_factor == 0.0:
975 log.append((
976 'warning',
977 'No pole-zero normalization factor given. '
978 'Assuming a value of 1.0',
979 context))
981 nfactor = 1.0
982 else:
983 nfactor = self.normalization_factor
985 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
986 if not is_digital:
987 resp = response.PoleZeroResponse(
988 constant=nfactor*cfactor,
989 zeros=[z.value()*factor for z in self.zero_list],
990 poles=[p.value()*factor for p in self.pole_list])
991 else:
992 if not deltat:
993 log.append((
994 'error',
995 'Digital response requires knowledge about sampling '
996 'interval. Response will be unusable.',
997 context))
999 resp = response.DigitalPoleZeroResponse(
1000 constant=nfactor*cfactor,
1001 zeros=[z.value()*factor for z in self.zero_list],
1002 poles=[p.value()*factor for p in self.pole_list],
1003 deltat=deltat or 0.0)
1005 if not self.normalization_frequency:
1006 log.append((
1007 'warning',
1008 'Cannot check pole-zero normalization factor: '
1009 'No normalization frequency given.',
1010 context))
1012 else:
1013 if is_digital and not deltat:
1014 log.append((
1015 'warning',
1016 'Cannot check computed vs reported normalization '
1017 'factor without knowing the sampling interval.',
1018 context))
1019 else:
1020 computed_normalization_factor = nfactor / abs(evaluate1(
1021 resp, self.normalization_frequency.value))
1023 db = 20.0 * num.log10(
1024 computed_normalization_factor / nfactor)
1026 if abs(db) > 0.17:
1027 log.append((
1028 'warning',
1029 'Computed and reported normalization factors differ '
1030 'by %g dB: computed: %g, reported: %g' % (
1031 db,
1032 computed_normalization_factor,
1033 nfactor),
1034 context))
1036 return Delivery([resp], log)
1039class ResponseListElement(Object):
1040 frequency = Frequency.T(xmltagname='Frequency')
1041 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1042 phase = Angle.T(xmltagname='Phase')
1045class Polynomial(BaseFilter):
1046 '''
1047 Response: expressed as a polynomial (allows non-linear sensors to be
1048 described). Corresponds to SEED blockette 62. Can be used to describe a
1049 stage of acquisition or a complete system.
1050 '''
1052 approximation_type = Approximation.T(default='MACLAURIN',
1053 xmltagname='ApproximationType')
1054 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1055 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1056 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1057 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1058 maximum_error = Float.T(xmltagname='MaximumError')
1059 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1061 def summary(self):
1062 return 'poly(%i)' % len(self.coefficient_list)
1065class Decimation(Object):
1066 '''
1067 Corresponds to SEED blockette 57.
1068 '''
1070 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1071 factor = Int.T(xmltagname='Factor')
1072 offset = Int.T(xmltagname='Offset')
1073 delay = FloatWithUnit.T(xmltagname='Delay')
1074 correction = FloatWithUnit.T(xmltagname='Correction')
1076 def summary(self):
1077 return 'deci(%i, %g -> %g, %g)' % (
1078 self.factor,
1079 self.input_sample_rate.value,
1080 self.input_sample_rate.value / self.factor,
1081 self.delay.value)
1083 def get_pyrocko_response(self):
1084 if self.delay and self.delay.value != 0.0:
1085 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1086 else:
1087 return Delivery([])
1090class Operator(Object):
1091 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1092 contact_list = List.T(Person.T(xmltagname='Contact'))
1093 web_site = String.T(optional=True, xmltagname='WebSite')
1096class Comment(Object):
1097 '''
1098 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1099 and 59.
1100 '''
1102 id = Counter.T(optional=True, xmlstyle='attribute')
1103 value = Unicode.T(xmltagname='Value')
1104 begin_effective_time = DummyAwareOptionalTimestamp.T(
1105 optional=True,
1106 xmltagname='BeginEffectiveTime')
1107 end_effective_time = DummyAwareOptionalTimestamp.T(
1108 optional=True,
1109 xmltagname='EndEffectiveTime')
1110 author_list = List.T(Person.T(xmltagname='Author'))
1113class ResponseList(BaseFilter):
1114 '''
1115 Response: list of frequency, amplitude and phase values. Corresponds to
1116 SEED blockette 55.
1117 '''
1119 response_list_element_list = List.T(
1120 ResponseListElement.T(xmltagname='ResponseListElement'))
1122 def summary(self):
1123 return 'list(%i)' % len(self.response_list_element_list)
1126class Log(Object):
1127 '''
1128 Container for log entries.
1129 '''
1131 entry_list = List.T(Comment.T(xmltagname='Entry'))
1134class ResponseStage(Object):
1135 '''
1136 This complex type represents channel response and covers SEED blockettes 53
1137 to 56.
1138 '''
1140 number = Counter.T(xmlstyle='attribute')
1141 resource_id = String.T(optional=True, xmlstyle='attribute')
1142 poles_zeros_list = List.T(
1143 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1144 coefficients_list = List.T(
1145 Coefficients.T(optional=True, xmltagname='Coefficients'))
1146 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1147 fir = FIR.T(optional=True, xmltagname='FIR')
1148 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1149 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1150 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1152 def summary(self):
1153 elements = []
1155 for stuff in [
1156 self.poles_zeros_list, self.coefficients_list,
1157 self.response_list, self.fir, self.polynomial,
1158 self.decimation, self.stage_gain]:
1160 if stuff:
1161 if isinstance(stuff, Object):
1162 elements.append(stuff.summary())
1163 else:
1164 elements.extend(obj.summary() for obj in stuff)
1166 return '%i: %s %s -> %s' % (
1167 self.number,
1168 ', '.join(elements),
1169 self.input_units.name.upper() if self.input_units else '?',
1170 self.output_units.name.upper() if self.output_units else '?')
1172 def get_squirrel_response_stage(self, context):
1173 from pyrocko.squirrel.model import ResponseStage
1175 delivery = Delivery()
1176 delivery_pr = self.get_pyrocko_response(context)
1177 log = delivery_pr.log
1178 delivery_pr.log = []
1179 elements = delivery.extend_without_payload(delivery_pr)
1181 delivery.payload = [ResponseStage(
1182 input_quantity=to_quantity(self.input_units, context, delivery),
1183 output_quantity=to_quantity(self.output_units, context, delivery),
1184 input_sample_rate=self.input_sample_rate,
1185 output_sample_rate=self.output_sample_rate,
1186 elements=elements,
1187 log=log)]
1189 return delivery
1191 def get_pyrocko_response(self, context, gain_only=False):
1193 context = context + ', stage %i' % self.number
1195 responses = []
1196 log = []
1197 if self.stage_gain:
1198 normalization_frequency = self.stage_gain.frequency or 0.0
1199 else:
1200 normalization_frequency = 0.0
1202 if not gain_only:
1203 deltat = None
1204 delay_responses = []
1205 if self.decimation:
1206 rate = self.decimation.input_sample_rate.value
1207 if rate > 0.0:
1208 deltat = 1.0 / rate
1209 delivery = self.decimation.get_pyrocko_response()
1210 if delivery.errors:
1211 return delivery
1213 delay_responses = delivery.payload
1214 log.extend(delivery.log)
1216 for pzs in self.poles_zeros_list:
1217 delivery = pzs.get_pyrocko_response(context, deltat)
1218 if delivery.errors:
1219 return delivery
1221 pz_resps = delivery.payload
1222 log.extend(delivery.log)
1223 responses.extend(pz_resps)
1225 # emulate incorrect? evalresp behaviour
1226 if pzs.normalization_frequency is not None and \
1227 pzs.normalization_frequency.value \
1228 != normalization_frequency \
1229 and normalization_frequency != 0.0:
1231 try:
1232 trial = response.MultiplyResponse(pz_resps)
1233 anorm = num.abs(evaluate1(
1234 trial, pzs.normalization_frequency.value))
1235 asens = num.abs(
1236 evaluate1(trial, normalization_frequency))
1238 factor = anorm/asens
1240 if abs(factor - 1.0) > 0.01:
1241 log.append((
1242 'warning',
1243 'PZ normalization frequency (%g) is different '
1244 'from stage gain frequency (%s) -> Emulating '
1245 'possibly incorrect evalresp behaviour. '
1246 'Correction factor: %g' % (
1247 pzs.normalization_frequency.value,
1248 normalization_frequency,
1249 factor),
1250 context))
1252 responses.append(
1253 response.PoleZeroResponse(constant=factor))
1254 except response.InvalidResponseError as e:
1255 log.append((
1256 'warning',
1257 'Could not check response: %s' % str(e),
1258 context))
1260 if len(self.poles_zeros_list) > 1:
1261 log.append((
1262 'warning',
1263 'Multiple poles and zeros records in single response '
1264 'stage.',
1265 context))
1267 for cfs in self.coefficients_list + (
1268 [self.fir] if self.fir else []):
1270 delivery = cfs.get_pyrocko_response(
1271 context, deltat, delay_responses,
1272 normalization_frequency)
1274 if delivery.errors:
1275 return delivery
1277 responses.extend(delivery.payload)
1278 log.extend(delivery.log)
1280 if len(self.coefficients_list) > 1:
1281 log.append((
1282 'warning',
1283 'Multiple filter coefficients lists in single response '
1284 'stage.',
1285 context))
1287 if self.response_list:
1288 log.append((
1289 'warning',
1290 'Unhandled response element of type: ResponseList',
1291 context))
1293 if self.polynomial:
1294 log.append((
1295 'warning',
1296 'Unhandled response element of type: Polynomial',
1297 context))
1299 if self.stage_gain:
1300 responses.append(
1301 response.PoleZeroResponse(constant=self.stage_gain.value))
1303 return Delivery(responses, log)
1305 @property
1306 def input_units(self):
1307 for e in (self.poles_zeros_list + self.coefficients_list +
1308 [self.response_list, self.fir, self.polynomial]):
1309 if e is not None:
1310 return e.input_units
1312 return None
1314 @property
1315 def output_units(self):
1316 for e in (self.poles_zeros_list + self.coefficients_list +
1317 [self.response_list, self.fir, self.polynomial]):
1318 if e is not None:
1319 return e.output_units
1321 return None
1323 @property
1324 def input_sample_rate(self):
1325 if self.decimation:
1326 return self.decimation.input_sample_rate.value
1328 return None
1330 @property
1331 def output_sample_rate(self):
1332 if self.decimation:
1333 return self.decimation.input_sample_rate.value \
1334 / self.decimation.factor
1336 return None
1339class Response(Object):
1340 resource_id = String.T(optional=True, xmlstyle='attribute')
1341 instrument_sensitivity = Sensitivity.T(optional=True,
1342 xmltagname='InstrumentSensitivity')
1343 instrument_polynomial = Polynomial.T(optional=True,
1344 xmltagname='InstrumentPolynomial')
1345 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1347 @property
1348 def output_sample_rate(self):
1349 if self.stage_list:
1350 return self.stage_list[-1].output_sample_rate
1351 else:
1352 return None
1354 def check_sample_rates(self, channel):
1356 if self.stage_list:
1357 sample_rate = None
1359 for stage in self.stage_list:
1360 if stage.decimation:
1361 input_sample_rate = \
1362 stage.decimation.input_sample_rate.value
1364 if sample_rate is not None and not same_sample_rate(
1365 sample_rate, input_sample_rate):
1367 logger.warning(
1368 'Response stage %i has unexpected input sample '
1369 'rate: %g Hz (expected: %g Hz)' % (
1370 stage.number,
1371 input_sample_rate,
1372 sample_rate))
1374 sample_rate = input_sample_rate / stage.decimation.factor
1376 if sample_rate is not None and channel.sample_rate \
1377 and not same_sample_rate(
1378 sample_rate, channel.sample_rate.value):
1380 logger.warning(
1381 'Channel sample rate (%g Hz) does not match sample rate '
1382 'deducted from response stages information (%g Hz).' % (
1383 channel.sample_rate.value,
1384 sample_rate))
1386 def check_units(self):
1388 if self.instrument_sensitivity \
1389 and self.instrument_sensitivity.input_units:
1391 units = self.instrument_sensitivity.input_units.name.upper()
1393 if self.stage_list:
1394 for stage in self.stage_list:
1395 if units and stage.input_units \
1396 and stage.input_units.name.upper() != units:
1398 logger.warning(
1399 'Input units of stage %i (%s) do not match %s (%s).'
1400 % (
1401 stage.number,
1402 units,
1403 'output units of stage %i'
1404 if stage.number == 0
1405 else 'sensitivity input units',
1406 units))
1408 if stage.output_units:
1409 units = stage.output_units.name.upper()
1410 else:
1411 units = None
1413 sout_units = self.instrument_sensitivity.output_units
1414 if self.instrument_sensitivity and sout_units:
1415 if units is not None and units != sout_units.name.upper():
1416 logger.warning(
1417 'Output units of stage %i (%s) do not match %s (%s).'
1418 % (
1419 stage.number,
1420 units,
1421 'sensitivity output units',
1422 sout_units.name.upper()))
1424 def _sensitivity_checkpoints(self, responses, context):
1425 delivery = Delivery()
1427 if self.instrument_sensitivity:
1428 sval = self.instrument_sensitivity.value
1429 sfreq = self.instrument_sensitivity.frequency
1430 if sval is None:
1431 delivery.log.append((
1432 'warning',
1433 'No sensitivity value given.',
1434 context))
1436 elif sval is None:
1437 delivery.log.append((
1438 'warning',
1439 'Reported sensitivity value is zero.',
1440 context))
1442 elif sfreq is None:
1443 delivery.log.append((
1444 'warning',
1445 'Sensitivity frequency not given.',
1446 context))
1448 else:
1449 trial = response.MultiplyResponse(responses)
1451 delivery.extend(
1452 check_resp(
1453 trial, sval, sfreq, 0.1,
1454 'Instrument sensitivity value inconsistent with '
1455 'sensitivity computed from complete response.',
1456 context))
1458 delivery.payload.append(response.FrequencyResponseCheckpoint(
1459 frequency=sfreq, value=sval))
1461 return delivery
1463 def get_squirrel_response(self, context, **kwargs):
1464 from pyrocko.squirrel.model import Response
1466 if self.stage_list:
1467 delivery = Delivery()
1468 for istage, stage in enumerate(self.stage_list):
1469 delivery.extend(stage.get_squirrel_response_stage(context))
1471 checkpoints = []
1472 if not delivery.errors:
1473 all_responses = []
1474 for sq_stage in delivery.payload:
1475 all_responses.extend(sq_stage.elements)
1477 checkpoints.extend(delivery.extend_without_payload(
1478 self._sensitivity_checkpoints(all_responses, context)))
1480 sq_stages = delivery.payload
1481 if sq_stages:
1482 if sq_stages[0].input_quantity is None \
1483 and self.instrument_sensitivity is not None:
1485 sq_stages[0].input_quantity = to_quantity(
1486 self.instrument_sensitivity.input_units,
1487 context, delivery)
1488 sq_stages[-1].output_quantity = to_quantity(
1489 self.instrument_sensitivity.output_units,
1490 context, delivery)
1492 sq_stages = delivery.expect()
1494 return Response(
1495 stages=sq_stages,
1496 log=delivery.log,
1497 checkpoints=checkpoints,
1498 **kwargs)
1500 elif self.instrument_sensitivity:
1501 raise NoResponseInformation(
1502 "Only instrument sensitivity given (won't use it). (%s)."
1503 % context)
1504 else:
1505 raise NoResponseInformation(
1506 'Empty instrument response (%s).'
1507 % context)
1509 def get_pyrocko_response(
1510 self, context, fake_input_units=None, stages=(0, 1)):
1512 delivery = Delivery()
1513 if self.stage_list:
1514 for istage, stage in enumerate(self.stage_list):
1515 delivery.extend(stage.get_pyrocko_response(
1516 context, gain_only=not (
1517 stages is None or stages[0] <= istage < stages[1])))
1519 elif self.instrument_sensitivity:
1520 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1522 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1523 checkpoints = delivery.extend_without_payload(delivery_cp)
1524 if delivery.errors:
1525 return delivery
1527 if fake_input_units is not None:
1528 if not self.instrument_sensitivity or \
1529 self.instrument_sensitivity.input_units is None:
1531 delivery.errors.append((
1532 'NoResponseInformation',
1533 'No input units given, so cannot convert to requested '
1534 'units: %s' % fake_input_units.upper(),
1535 context))
1537 return delivery
1539 input_units = self.instrument_sensitivity.input_units.name.upper()
1541 conresp = None
1542 try:
1543 conresp = conversion[
1544 fake_input_units.upper(), input_units]
1546 except KeyError:
1547 delivery.errors.append((
1548 'NoResponseInformation',
1549 'Cannot convert between units: %s, %s'
1550 % (fake_input_units, input_units),
1551 context))
1553 if conresp is not None:
1554 delivery.payload.append(conresp)
1555 for checkpoint in checkpoints:
1556 checkpoint.value *= num.abs(evaluate1(
1557 conresp, checkpoint.frequency))
1559 delivery.payload = [
1560 response.MultiplyResponse(
1561 delivery.payload,
1562 checkpoints=checkpoints)]
1564 return delivery
1566 @classmethod
1567 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1568 normalization_frequency=1.0):
1569 '''
1570 Convert Pyrocko pole-zero response to StationXML response.
1572 :param presponse: Pyrocko pole-zero response
1573 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1574 :param input_unit: Input unit to be reported in the StationXML
1575 response.
1576 :type input_unit: str
1577 :param output_unit: Output unit to be reported in the StationXML
1578 response.
1579 :type output_unit: str
1580 :param normalization_frequency: Frequency where the normalization
1581 factor for the StationXML response should be computed.
1582 :type normalization_frequency: float
1583 '''
1585 norm_factor = 1.0/float(abs(
1586 evaluate1(presponse, normalization_frequency)
1587 / presponse.constant))
1589 pzs = PolesZeros(
1590 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1591 normalization_factor=norm_factor,
1592 normalization_frequency=Frequency(normalization_frequency),
1593 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1594 imaginary=FloatNoUnit(z.imag))
1595 for z in presponse.zeros],
1596 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1597 imaginary=FloatNoUnit(z.imag))
1598 for z in presponse.poles])
1600 pzs.validate()
1602 stage = ResponseStage(
1603 number=1,
1604 poles_zeros_list=[pzs],
1605 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1607 resp = Response(
1608 instrument_sensitivity=Sensitivity(
1609 value=stage.stage_gain.value,
1610 frequency=normalization_frequency,
1611 input_units=Units(input_unit),
1612 output_units=Units(output_unit)),
1614 stage_list=[stage])
1616 return resp
1619class BaseNode(Object):
1620 '''
1621 A base node type for derivation from: Network, Station and Channel types.
1622 '''
1624 code = String.T(xmlstyle='attribute')
1625 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1626 xmlstyle='attribute')
1627 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1628 xmlstyle='attribute')
1629 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1630 alternate_code = String.T(optional=True, xmlstyle='attribute')
1631 historical_code = String.T(optional=True, xmlstyle='attribute')
1632 description = Unicode.T(optional=True, xmltagname='Description')
1633 comment_list = List.T(Comment.T(xmltagname='Comment'))
1635 def spans(self, *args):
1636 if len(args) == 0:
1637 return True
1638 elif len(args) == 1:
1639 return ((self.start_date is None or
1640 self.start_date <= args[0]) and
1641 (self.end_date is None or
1642 args[0] <= self.end_date))
1644 elif len(args) == 2:
1645 return ((self.start_date is None or
1646 args[1] >= self.start_date) and
1647 (self.end_date is None or
1648 self.end_date >= args[0]))
1651def overlaps(a, b):
1652 return (
1653 a.start_date is None or b.end_date is None
1654 or a.start_date < b.end_date
1655 ) and (
1656 b.start_date is None or a.end_date is None
1657 or b.start_date < a.end_date
1658 )
1661def check_overlaps(node_type_name, codes, nodes):
1662 errors = []
1663 for ia, a in enumerate(nodes):
1664 for b in nodes[ia+1:]:
1665 if overlaps(a, b):
1666 errors.append(
1667 '%s epochs overlap for %s:\n'
1668 ' %s - %s\n %s - %s' % (
1669 node_type_name,
1670 '.'.join(codes),
1671 tts(a.start_date), tts(a.end_date),
1672 tts(b.start_date), tts(b.end_date)))
1674 return errors
1677class Channel(BaseNode):
1678 '''
1679 Equivalent to SEED blockette 52 and parent element for the related the
1680 response blockettes.
1681 '''
1683 location_code = String.T(xmlstyle='attribute')
1684 external_reference_list = List.T(
1685 ExternalReference.T(xmltagname='ExternalReference'))
1686 latitude = Latitude.T(xmltagname='Latitude')
1687 longitude = Longitude.T(xmltagname='Longitude')
1688 elevation = Distance.T(xmltagname='Elevation')
1689 depth = Distance.T(xmltagname='Depth')
1690 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1691 dip = Dip.T(optional=True, xmltagname='Dip')
1692 type_list = List.T(Type.T(xmltagname='Type'))
1693 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1694 sample_rate_ratio = SampleRateRatio.T(optional=True,
1695 xmltagname='SampleRateRatio')
1696 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1697 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1698 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1699 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1700 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1701 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1702 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1703 response = Response.T(optional=True, xmltagname='Response')
1705 @property
1706 def position_values(self):
1707 lat = self.latitude.value
1708 lon = self.longitude.value
1709 elevation = value_or_none(self.elevation)
1710 depth = value_or_none(self.depth)
1711 return lat, lon, elevation, depth
1714class Station(BaseNode):
1715 '''
1716 This type represents a Station epoch. It is common to only have a single
1717 station epoch with the station's creation and termination dates as the
1718 epoch start and end dates.
1719 '''
1721 latitude = Latitude.T(xmltagname='Latitude')
1722 longitude = Longitude.T(xmltagname='Longitude')
1723 elevation = Distance.T(xmltagname='Elevation')
1724 site = Site.T(default=Site.D(name=''), xmltagname='Site')
1725 vault = Unicode.T(optional=True, xmltagname='Vault')
1726 geology = Unicode.T(optional=True, xmltagname='Geology')
1727 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1728 operator_list = List.T(Operator.T(xmltagname='Operator'))
1729 creation_date = DummyAwareOptionalTimestamp.T(
1730 optional=True, xmltagname='CreationDate')
1731 termination_date = DummyAwareOptionalTimestamp.T(
1732 optional=True, xmltagname='TerminationDate')
1733 total_number_channels = Counter.T(
1734 optional=True, xmltagname='TotalNumberChannels')
1735 selected_number_channels = Counter.T(
1736 optional=True, xmltagname='SelectedNumberChannels')
1737 external_reference_list = List.T(
1738 ExternalReference.T(xmltagname='ExternalReference'))
1739 channel_list = List.T(Channel.T(xmltagname='Channel'))
1741 @property
1742 def position_values(self):
1743 lat = self.latitude.value
1744 lon = self.longitude.value
1745 elevation = value_or_none(self.elevation)
1746 return lat, lon, elevation
1749class Network(BaseNode):
1750 '''
1751 This type represents the Network layer, all station metadata is contained
1752 within this element. The official name of the network or other descriptive
1753 information can be included in the Description element. The Network can
1754 contain 0 or more Stations.
1755 '''
1757 total_number_stations = Counter.T(optional=True,
1758 xmltagname='TotalNumberStations')
1759 selected_number_stations = Counter.T(optional=True,
1760 xmltagname='SelectedNumberStations')
1761 station_list = List.T(Station.T(xmltagname='Station'))
1763 @property
1764 def station_code_list(self):
1765 return sorted(set(s.code for s in self.station_list))
1767 @property
1768 def sl_code_list(self):
1769 sls = set()
1770 for station in self.station_list:
1771 for channel in station.channel_list:
1772 sls.add((station.code, channel.location_code))
1774 return sorted(sls)
1776 def summary(self, width=80, indent=4):
1777 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1778 lines = ['%s (%i):' % (self.code, len(sls))]
1779 if sls:
1780 ssls = ['.'.join(x for x in c if x) for c in sls]
1781 w = max(len(x) for x in ssls)
1782 n = (width - indent) / (w+1)
1783 while ssls:
1784 lines.append(
1785 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1787 ssls[:n] = []
1789 return '\n'.join(lines)
1792def value_or_none(x):
1793 if x is not None:
1794 return x.value
1795 else:
1796 return None
1799def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1801 pos = lat, lon, elevation, depth = \
1802 channels[0].position_values
1804 if not all(pos == x.position_values for x in channels):
1805 info = '\n'.join(
1806 ' %s: %s' % (x.code, x.position_values) for
1807 x in channels)
1809 mess = 'encountered inconsistencies in channel ' \
1810 'lat/lon/elevation/depth ' \
1811 'for %s.%s.%s: \n%s' % (nsl + (info,))
1813 if inconsistencies == 'raise':
1814 raise InconsistentChannelLocations(mess)
1816 elif inconsistencies == 'warn':
1817 logger.warning(mess)
1818 logger.warning(' -> using mean values')
1820 apos = num.array([x.position_values for x in channels], dtype=float)
1821 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1822 / num.sum(num.isfinite(apos), axis=0)
1824 groups = {}
1825 for channel in channels:
1826 if channel.code not in groups:
1827 groups[channel.code] = []
1829 groups[channel.code].append(channel)
1831 pchannels = []
1832 for code in sorted(groups.keys()):
1833 data = [
1834 (channel.code, value_or_none(channel.azimuth),
1835 value_or_none(channel.dip))
1836 for channel in groups[code]]
1838 azimuth, dip = util.consistency_merge(
1839 data,
1840 message='channel orientation values differ:',
1841 error=inconsistencies)
1843 pchannels.append(
1844 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1846 return pyrocko.model.Station(
1847 *nsl,
1848 lat=mlat,
1849 lon=mlon,
1850 elevation=mele,
1851 depth=mdep,
1852 channels=pchannels)
1855class FDSNStationXML(Object):
1856 '''
1857 Top-level type for Station XML. Required field are Source (network ID of
1858 the institution sending the message) and one or more Network containers or
1859 one or more Station containers.
1860 '''
1862 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1863 source = String.T(xmltagname='Source')
1864 sender = String.T(optional=True, xmltagname='Sender')
1865 module = String.T(optional=True, xmltagname='Module')
1866 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1867 created = Timestamp.T(optional=True, xmltagname='Created')
1868 network_list = List.T(Network.T(xmltagname='Network'))
1870 xmltagname = 'FDSNStationXML'
1871 guessable_xmlns = [guts_xmlns]
1873 def __init__(self, *args, **kwargs):
1874 Object.__init__(self, *args, **kwargs)
1875 if self.created is None:
1876 self.created = time.time()
1878 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1879 time=None, timespan=None,
1880 inconsistencies='warn'):
1882 assert inconsistencies in ('raise', 'warn')
1884 if nslcs is not None:
1885 nslcs = set(nslcs)
1887 if nsls is not None:
1888 nsls = set(nsls)
1890 tt = ()
1891 if time is not None:
1892 tt = (time,)
1893 elif timespan is not None:
1894 tt = timespan
1896 ns_have = set()
1897 pstations = []
1898 nsltt_to_channels = defaultdict(list)
1899 for network in self.network_list:
1900 if not network.spans(*tt):
1901 continue
1903 for station in network.station_list:
1904 if not station.spans(*tt):
1905 continue
1907 if station.channel_list:
1908 loc_to_channels = {}
1909 for channel in station.channel_list:
1910 if not channel.spans(*tt):
1911 continue
1913 loc = channel.location_code.strip()
1914 if loc not in loc_to_channels:
1915 loc_to_channels[loc] = []
1917 loc_to_channels[loc].append(channel)
1919 for loc in sorted(loc_to_channels.keys()):
1920 channels = loc_to_channels[loc]
1921 if nslcs is not None:
1922 channels = [channel for channel in channels
1923 if (network.code, station.code, loc,
1924 channel.code) in nslcs]
1926 if not channels:
1927 continue
1929 nsl = network.code, station.code, loc
1930 if nsls is not None and nsl not in nsls:
1931 continue
1933 for channel in channels:
1934 k = nsl, channel.start_date, channel.end_date
1935 nsltt_to_channels[k].append(channel)
1937 else:
1938 ns = (network.code, station.code)
1939 if ns in ns_have:
1940 message = 'Duplicate station ' \
1941 '(multiple epochs match): %s.%s ' % ns
1943 if inconsistencies == 'raise':
1944 raise Inconsistencies(message)
1945 else:
1946 logger.warning(message)
1948 ns_have.add(ns)
1950 pstations.append(pyrocko.model.Station(
1951 network.code, station.code, '*',
1952 lat=station.latitude.value,
1953 lon=station.longitude.value,
1954 elevation=value_or_none(station.elevation),
1955 name=station.description or ''))
1957 nsl_have = set()
1958 for (nsl, _, _), channels in nsltt_to_channels.items():
1959 if nsl in nsl_have:
1960 message = 'Duplicate station ' \
1961 '(multiple epochs match): %s.%s.%s' % nsl
1963 if inconsistencies == 'raise':
1964 raise Inconsistencies(message)
1965 else:
1966 logger.warning(message)
1968 nsl_have.add(nsl)
1970 pstations.append(
1971 pyrocko_station_from_channels(
1972 nsl,
1973 channels,
1974 inconsistencies=inconsistencies))
1976 return pstations
1978 @classmethod
1979 def from_pyrocko_stations(
1980 cls, pyrocko_stations, add_flat_responses_from=None):
1982 '''
1983 Generate :py:class:`FDSNStationXML` from list of
1984 :py:class;`pyrocko.model.Station` instances.
1986 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1987 instances.
1988 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1989 '''
1990 network_dict = defaultdict(list)
1992 if add_flat_responses_from:
1993 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1994 extra = dict(
1995 response=Response(
1996 instrument_sensitivity=Sensitivity(
1997 value=1.0,
1998 frequency=1.0,
1999 input_units=Units(name=add_flat_responses_from))))
2000 else:
2001 extra = {}
2003 have_offsets = set()
2004 for s in pyrocko_stations:
2006 if s.north_shift != 0.0 or s.east_shift != 0.0:
2007 have_offsets.add(s.nsl())
2009 network, station, location = s.nsl()
2010 channel_list = []
2011 for c in s.channels:
2012 channel_list.append(
2013 Channel(
2014 location_code=location,
2015 code=c.name,
2016 latitude=Latitude(value=s.effective_lat),
2017 longitude=Longitude(value=s.effective_lon),
2018 elevation=Distance(value=s.elevation),
2019 depth=Distance(value=s.depth),
2020 azimuth=Azimuth(value=c.azimuth),
2021 dip=Dip(value=c.dip),
2022 **extra
2023 )
2024 )
2026 network_dict[network].append(
2027 Station(
2028 code=station,
2029 latitude=Latitude(value=s.effective_lat),
2030 longitude=Longitude(value=s.effective_lon),
2031 elevation=Distance(value=s.elevation),
2032 channel_list=channel_list)
2033 )
2035 if have_offsets:
2036 logger.warning(
2037 'StationXML does not support Cartesian offsets in '
2038 'coordinates. Storing effective lat/lon for stations: %s' %
2039 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2041 timestamp = util.to_time_float(time.time())
2042 network_list = []
2043 for k, station_list in network_dict.items():
2045 network_list.append(
2046 Network(
2047 code=k, station_list=station_list,
2048 total_number_stations=len(station_list)))
2050 sxml = FDSNStationXML(
2051 source='from pyrocko stations list', created=timestamp,
2052 network_list=network_list)
2054 sxml.validate()
2055 return sxml
2057 def iter_network_stations(
2058 self, net=None, sta=None, time=None, timespan=None):
2060 tt = ()
2061 if time is not None:
2062 tt = (time,)
2063 elif timespan is not None:
2064 tt = timespan
2066 for network in self.network_list:
2067 if not network.spans(*tt) or (
2068 net is not None and network.code != net):
2069 continue
2071 for station in network.station_list:
2072 if not station.spans(*tt) or (
2073 sta is not None and station.code != sta):
2074 continue
2076 yield (network, station)
2078 def iter_network_station_channels(
2079 self, net=None, sta=None, loc=None, cha=None,
2080 time=None, timespan=None):
2082 if loc is not None:
2083 loc = loc.strip()
2085 tt = ()
2086 if time is not None:
2087 tt = (time,)
2088 elif timespan is not None:
2089 tt = timespan
2091 for network in self.network_list:
2092 if not network.spans(*tt) or (
2093 net is not None and network.code != net):
2094 continue
2096 for station in network.station_list:
2097 if not station.spans(*tt) or (
2098 sta is not None and station.code != sta):
2099 continue
2101 if station.channel_list:
2102 for channel in station.channel_list:
2103 if (not channel.spans(*tt) or
2104 (cha is not None and channel.code != cha) or
2105 (loc is not None and
2106 channel.location_code.strip() != loc)):
2107 continue
2109 yield (network, station, channel)
2111 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2112 time=None, timespan=None):
2114 groups = {}
2115 for network, station, channel in self.iter_network_station_channels(
2116 net, sta, loc, cha, time=time, timespan=timespan):
2118 net = network.code
2119 sta = station.code
2120 cha = channel.code
2121 loc = channel.location_code.strip()
2122 if len(cha) == 3:
2123 bic = cha[:2] # band and intrument code according to SEED
2124 elif len(cha) == 1:
2125 bic = ''
2126 else:
2127 bic = cha
2129 if channel.response and \
2130 channel.response.instrument_sensitivity and \
2131 channel.response.instrument_sensitivity.input_units:
2133 unit = channel.response.instrument_sensitivity\
2134 .input_units.name.upper()
2135 else:
2136 unit = None
2138 bic = (bic, unit)
2140 k = net, sta, loc
2141 if k not in groups:
2142 groups[k] = {}
2144 if bic not in groups[k]:
2145 groups[k][bic] = []
2147 groups[k][bic].append(channel)
2149 for nsl, bic_to_channels in groups.items():
2150 bad_bics = []
2151 for bic, channels in bic_to_channels.items():
2152 sample_rates = []
2153 for channel in channels:
2154 sample_rates.append(channel.sample_rate.value)
2156 if not same(sample_rates):
2157 scs = ','.join(channel.code for channel in channels)
2158 srs = ', '.join('%e' % x for x in sample_rates)
2159 err = 'ignoring channels with inconsistent sampling ' + \
2160 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2162 logger.warning(err)
2163 bad_bics.append(bic)
2165 for bic in bad_bics:
2166 del bic_to_channels[bic]
2168 return groups
2170 def choose_channels(
2171 self,
2172 target_sample_rate=None,
2173 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2174 priority_units=['M/S', 'M/S**2'],
2175 priority_instrument_code=['H', 'L'],
2176 time=None,
2177 timespan=None):
2179 nslcs = {}
2180 for nsl, bic_to_channels in self.get_channel_groups(
2181 time=time, timespan=timespan).items():
2183 useful_bics = []
2184 for bic, channels in bic_to_channels.items():
2185 rate = channels[0].sample_rate.value
2187 if target_sample_rate is not None and \
2188 rate < target_sample_rate*0.99999:
2189 continue
2191 if len(bic[0]) == 2:
2192 if bic[0][0] not in priority_band_code:
2193 continue
2195 if bic[0][1] not in priority_instrument_code:
2196 continue
2198 unit = bic[1]
2200 prio_unit = len(priority_units)
2201 try:
2202 prio_unit = priority_units.index(unit)
2203 except ValueError:
2204 pass
2206 prio_inst = len(priority_instrument_code)
2207 prio_band = len(priority_band_code)
2208 if len(channels[0].code) == 3:
2209 try:
2210 prio_inst = priority_instrument_code.index(
2211 channels[0].code[1])
2212 except ValueError:
2213 pass
2215 try:
2216 prio_band = priority_band_code.index(
2217 channels[0].code[0])
2218 except ValueError:
2219 pass
2221 if target_sample_rate is None:
2222 rate = -rate
2224 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2225 prio_inst, bic))
2227 useful_bics.sort()
2229 for _, _, rate, _, _, bic in useful_bics:
2230 channels = sorted(
2231 bic_to_channels[bic],
2232 key=lambda channel: channel.code)
2234 if channels:
2235 for channel in channels:
2236 nslcs[nsl + (channel.code,)] = channel
2238 break
2240 return nslcs
2242 def get_pyrocko_response(
2243 self, nslc,
2244 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2246 net, sta, loc, cha = nslc
2247 resps = []
2248 for _, _, channel in self.iter_network_station_channels(
2249 net, sta, loc, cha, time=time, timespan=timespan):
2250 resp = channel.response
2251 if resp:
2252 resp.check_sample_rates(channel)
2253 resp.check_units()
2254 resps.append(resp.get_pyrocko_response(
2255 '.'.join(nslc),
2256 fake_input_units=fake_input_units,
2257 stages=stages).expect_one())
2259 if not resps:
2260 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2261 elif len(resps) > 1:
2262 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2264 return resps[0]
2266 @property
2267 def n_code_list(self):
2268 return sorted(set(x.code for x in self.network_list))
2270 @property
2271 def ns_code_list(self):
2272 nss = set()
2273 for network in self.network_list:
2274 for station in network.station_list:
2275 nss.add((network.code, station.code))
2277 return sorted(nss)
2279 @property
2280 def nsl_code_list(self):
2281 nsls = set()
2282 for network in self.network_list:
2283 for station in network.station_list:
2284 for channel in station.channel_list:
2285 nsls.add(
2286 (network.code, station.code, channel.location_code))
2288 return sorted(nsls)
2290 @property
2291 def nslc_code_list(self):
2292 nslcs = set()
2293 for network in self.network_list:
2294 for station in network.station_list:
2295 for channel in station.channel_list:
2296 nslcs.add(
2297 (network.code, station.code, channel.location_code,
2298 channel.code))
2300 return sorted(nslcs)
2302 def summary(self):
2303 lst = [
2304 'number of n codes: %i' % len(self.n_code_list),
2305 'number of ns codes: %i' % len(self.ns_code_list),
2306 'number of nsl codes: %i' % len(self.nsl_code_list),
2307 'number of nslc codes: %i' % len(self.nslc_code_list)
2308 ]
2309 return '\n'.join(lst)
2311 def summary_stages(self):
2312 data = []
2313 for network, station, channel in self.iter_network_station_channels():
2314 nslc = (network.code, station.code, channel.location_code,
2315 channel.code)
2317 stages = []
2318 in_units = '?'
2319 out_units = '?'
2320 if channel.response:
2321 sens = channel.response.instrument_sensitivity
2322 if sens:
2323 in_units = sens.input_units.name.upper()
2324 out_units = sens.output_units.name.upper()
2326 for stage in channel.response.stage_list:
2327 stages.append(stage.summary())
2329 data.append(
2330 (nslc, tts(channel.start_date), tts(channel.end_date),
2331 in_units, out_units, stages))
2333 data.sort()
2335 lst = []
2336 for nslc, stmin, stmax, in_units, out_units, stages in data:
2337 lst.append(' %s: %s - %s, %s -> %s' % (
2338 '.'.join(nslc), stmin, stmax, in_units, out_units))
2339 for stage in stages:
2340 lst.append(' %s' % stage)
2342 return '\n'.join(lst)
2344 def _check_overlaps(self):
2345 by_nslc = {}
2346 for network in self.network_list:
2347 for station in network.station_list:
2348 for channel in station.channel_list:
2349 nslc = (network.code, station.code, channel.location_code,
2350 channel.code)
2351 if nslc not in by_nslc:
2352 by_nslc[nslc] = []
2354 by_nslc[nslc].append(channel)
2356 errors = []
2357 for nslc, channels in by_nslc.items():
2358 errors.extend(check_overlaps('Channel', nslc, channels))
2360 return errors
2362 def check(self):
2363 errors = []
2364 for meth in [self._check_overlaps]:
2365 errors.extend(meth())
2367 if errors:
2368 raise Inconsistencies(
2369 'Inconsistencies found in StationXML:\n '
2370 + '\n '.join(errors))
2373def load_channel_table(stream):
2375 networks = {}
2376 stations = {}
2378 for line in stream:
2379 line = str(line.decode('ascii'))
2380 if line.startswith('#'):
2381 continue
2383 t = line.rstrip().split('|')
2385 if len(t) != 17:
2386 logger.warning('Invalid channel record: %s' % line)
2387 continue
2389 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2390 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2392 try:
2393 scale = float(scale)
2394 except ValueError:
2395 scale = None
2397 try:
2398 scale_freq = float(scale_freq)
2399 except ValueError:
2400 scale_freq = None
2402 try:
2403 depth = float(dep)
2404 except ValueError:
2405 depth = 0.0
2407 try:
2408 azi = float(azi)
2409 dip = float(dip)
2410 except ValueError:
2411 azi = None
2412 dip = None
2414 try:
2415 if net not in networks:
2416 network = Network(code=net)
2417 else:
2418 network = networks[net]
2420 if (net, sta) not in stations:
2421 station = Station(
2422 code=sta, latitude=lat, longitude=lon, elevation=ele)
2424 station.regularize()
2425 else:
2426 station = stations[net, sta]
2428 if scale:
2429 resp = Response(
2430 instrument_sensitivity=Sensitivity(
2431 value=scale,
2432 frequency=scale_freq,
2433 input_units=scale_units))
2434 else:
2435 resp = None
2437 channel = Channel(
2438 code=cha,
2439 location_code=loc.strip(),
2440 latitude=lat,
2441 longitude=lon,
2442 elevation=ele,
2443 depth=depth,
2444 azimuth=azi,
2445 dip=dip,
2446 sensor=Equipment(description=sens),
2447 response=resp,
2448 sample_rate=sample_rate,
2449 start_date=start_date,
2450 end_date=end_date or None)
2452 channel.regularize()
2454 except ValidationError:
2455 raise InvalidRecord(line)
2457 if net not in networks:
2458 networks[net] = network
2460 if (net, sta) not in stations:
2461 stations[net, sta] = station
2462 network.station_list.append(station)
2464 station.channel_list.append(channel)
2466 return FDSNStationXML(
2467 source='created from table input',
2468 created=time.time(),
2469 network_list=sorted(networks.values(), key=lambda x: x.code))
2472def primitive_merge(sxs):
2473 networks = []
2474 for sx in sxs:
2475 networks.extend(sx.network_list)
2477 return FDSNStationXML(
2478 source='merged from different sources',
2479 created=time.time(),
2480 network_list=copy.deepcopy(
2481 sorted(networks, key=lambda x: x.code)))
2484def split_channels(sx):
2485 for nslc in sx.nslc_code_list:
2486 network_list = sx.network_list
2487 network_list_filtered = [
2488 network for network in network_list
2489 if network.code == nslc[0]]
2491 for network in network_list_filtered:
2492 sx.network_list = [network]
2493 station_list = network.station_list
2494 station_list_filtered = [
2495 station for station in station_list
2496 if station.code == nslc[1]]
2498 for station in station_list_filtered:
2499 network.station_list = [station]
2500 channel_list = station.channel_list
2501 station.channel_list = [
2502 channel for channel in channel_list
2503 if (channel.location_code, channel.code) == nslc[2:4]]
2505 yield nslc, copy.deepcopy(sx)
2507 station.channel_list = channel_list
2509 network.station_list = station_list
2511 sx.network_list = network_list
2514if __name__ == '__main__':
2515 from optparse import OptionParser
2517 util.setup_logging('pyrocko.io.stationxml', 'warning')
2519 usage = \
2520 'python -m pyrocko.io.stationxml check|stats|stages ' \
2521 '<filename> [options]'
2523 description = '''Torture StationXML file.'''
2525 parser = OptionParser(
2526 usage=usage,
2527 description=description,
2528 formatter=util.BetterHelpFormatter())
2530 (options, args) = parser.parse_args(sys.argv[1:])
2532 if len(args) != 2:
2533 parser.print_help()
2534 sys.exit(1)
2536 action, path = args
2538 sx = load_xml(filename=path)
2539 if action == 'check':
2540 try:
2541 sx.check()
2542 except Inconsistencies as e:
2543 logger.error(e)
2544 sys.exit(1)
2546 elif action == 'stats':
2547 print(sx.summary())
2549 elif action == 'stages':
2550 print(sx.summary_stages())
2552 else:
2553 parser.print_help()
2554 sys.exit('unknown action: %s' % action)