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-23 12:04 +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 sensor_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.code[:-1], \
1935 channel.start_date, channel.end_date
1936 sensor_to_channels[k].append(channel)
1938 else:
1939 ns = (network.code, station.code)
1940 if ns in ns_have:
1941 message = 'Duplicate station ' \
1942 '(multiple epochs match): %s.%s ' % ns
1944 if inconsistencies == 'raise':
1945 raise Inconsistencies(message)
1946 else:
1947 logger.warning(message)
1949 ns_have.add(ns)
1951 pstations.append(pyrocko.model.Station(
1952 network.code, station.code, '*',
1953 lat=station.latitude.value,
1954 lon=station.longitude.value,
1955 elevation=value_or_none(station.elevation),
1956 name=station.description or ''))
1958 sensor_have = set()
1959 for (nsl, bi, _, _), channels in sensor_to_channels.items():
1960 if (nsl, bi) in sensor_have:
1961 message = 'Duplicate station ' \
1962 '(multiple epochs match): %s.%s.%s' % nsl
1964 if inconsistencies == 'raise':
1965 raise Inconsistencies(message)
1966 else:
1967 logger.warning(message)
1969 sensor_have.add((nsl, bi))
1971 pstations.append(
1972 pyrocko_station_from_channels(
1973 nsl,
1974 channels,
1975 inconsistencies=inconsistencies))
1977 return pstations
1979 @classmethod
1980 def from_pyrocko_stations(
1981 cls, pyrocko_stations, add_flat_responses_from=None):
1983 '''
1984 Generate :py:class:`FDSNStationXML` from list of
1985 :py:class;`pyrocko.model.Station` instances.
1987 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1988 instances.
1989 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1990 '''
1991 network_dict = defaultdict(list)
1993 if add_flat_responses_from:
1994 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1995 extra = dict(
1996 response=Response(
1997 instrument_sensitivity=Sensitivity(
1998 value=1.0,
1999 frequency=1.0,
2000 input_units=Units(name=add_flat_responses_from))))
2001 else:
2002 extra = {}
2004 have_offsets = set()
2005 for s in pyrocko_stations:
2007 if s.north_shift != 0.0 or s.east_shift != 0.0:
2008 have_offsets.add(s.nsl())
2010 network, station, location = s.nsl()
2011 channel_list = []
2012 for c in s.channels:
2013 channel_list.append(
2014 Channel(
2015 location_code=location,
2016 code=c.name,
2017 latitude=Latitude(value=s.effective_lat),
2018 longitude=Longitude(value=s.effective_lon),
2019 elevation=Distance(value=s.elevation),
2020 depth=Distance(value=s.depth),
2021 azimuth=Azimuth(value=c.azimuth),
2022 dip=Dip(value=c.dip),
2023 **extra
2024 )
2025 )
2027 network_dict[network].append(
2028 Station(
2029 code=station,
2030 latitude=Latitude(value=s.effective_lat),
2031 longitude=Longitude(value=s.effective_lon),
2032 elevation=Distance(value=s.elevation),
2033 channel_list=channel_list)
2034 )
2036 if have_offsets:
2037 logger.warning(
2038 'StationXML does not support Cartesian offsets in '
2039 'coordinates. Storing effective lat/lon for stations: %s' %
2040 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2042 timestamp = util.to_time_float(time.time())
2043 network_list = []
2044 for k, station_list in network_dict.items():
2046 network_list.append(
2047 Network(
2048 code=k, station_list=station_list,
2049 total_number_stations=len(station_list)))
2051 sxml = FDSNStationXML(
2052 source='from pyrocko stations list', created=timestamp,
2053 network_list=network_list)
2055 sxml.validate()
2056 return sxml
2058 def iter_network_stations(
2059 self, net=None, sta=None, time=None, timespan=None):
2061 tt = ()
2062 if time is not None:
2063 tt = (time,)
2064 elif timespan is not None:
2065 tt = timespan
2067 for network in self.network_list:
2068 if not network.spans(*tt) or (
2069 net is not None and network.code != net):
2070 continue
2072 for station in network.station_list:
2073 if not station.spans(*tt) or (
2074 sta is not None and station.code != sta):
2075 continue
2077 yield (network, station)
2079 def iter_network_station_channels(
2080 self, net=None, sta=None, loc=None, cha=None,
2081 time=None, timespan=None):
2083 if loc is not None:
2084 loc = loc.strip()
2086 tt = ()
2087 if time is not None:
2088 tt = (time,)
2089 elif timespan is not None:
2090 tt = timespan
2092 for network in self.network_list:
2093 if not network.spans(*tt) or (
2094 net is not None and network.code != net):
2095 continue
2097 for station in network.station_list:
2098 if not station.spans(*tt) or (
2099 sta is not None and station.code != sta):
2100 continue
2102 if station.channel_list:
2103 for channel in station.channel_list:
2104 if (not channel.spans(*tt) or
2105 (cha is not None and channel.code != cha) or
2106 (loc is not None and
2107 channel.location_code.strip() != loc)):
2108 continue
2110 yield (network, station, channel)
2112 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2113 time=None, timespan=None):
2115 groups = {}
2116 for network, station, channel in self.iter_network_station_channels(
2117 net, sta, loc, cha, time=time, timespan=timespan):
2119 net = network.code
2120 sta = station.code
2121 cha = channel.code
2122 loc = channel.location_code.strip()
2123 if len(cha) == 3:
2124 bic = cha[:2] # band and intrument code according to SEED
2125 elif len(cha) == 1:
2126 bic = ''
2127 else:
2128 bic = cha
2130 if channel.response and \
2131 channel.response.instrument_sensitivity and \
2132 channel.response.instrument_sensitivity.input_units:
2134 unit = channel.response.instrument_sensitivity\
2135 .input_units.name.upper()
2136 else:
2137 unit = None
2139 bic = (bic, unit)
2141 k = net, sta, loc
2142 if k not in groups:
2143 groups[k] = {}
2145 if bic not in groups[k]:
2146 groups[k][bic] = []
2148 groups[k][bic].append(channel)
2150 for nsl, bic_to_channels in groups.items():
2151 bad_bics = []
2152 for bic, channels in bic_to_channels.items():
2153 sample_rates = []
2154 for channel in channels:
2155 sample_rates.append(channel.sample_rate.value)
2157 if not same(sample_rates):
2158 scs = ','.join(channel.code for channel in channels)
2159 srs = ', '.join('%e' % x for x in sample_rates)
2160 err = 'ignoring channels with inconsistent sampling ' + \
2161 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2163 logger.warning(err)
2164 bad_bics.append(bic)
2166 for bic in bad_bics:
2167 del bic_to_channels[bic]
2169 return groups
2171 def choose_channels(
2172 self,
2173 target_sample_rate=None,
2174 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2175 priority_units=['M/S', 'M/S**2'],
2176 priority_instrument_code=['H', 'L'],
2177 time=None,
2178 timespan=None):
2180 nslcs = {}
2181 for nsl, bic_to_channels in self.get_channel_groups(
2182 time=time, timespan=timespan).items():
2184 useful_bics = []
2185 for bic, channels in bic_to_channels.items():
2186 rate = channels[0].sample_rate.value
2188 if target_sample_rate is not None and \
2189 rate < target_sample_rate*0.99999:
2190 continue
2192 if len(bic[0]) == 2:
2193 if bic[0][0] not in priority_band_code:
2194 continue
2196 if bic[0][1] not in priority_instrument_code:
2197 continue
2199 unit = bic[1]
2201 prio_unit = len(priority_units)
2202 try:
2203 prio_unit = priority_units.index(unit)
2204 except ValueError:
2205 pass
2207 prio_inst = len(priority_instrument_code)
2208 prio_band = len(priority_band_code)
2209 if len(channels[0].code) == 3:
2210 try:
2211 prio_inst = priority_instrument_code.index(
2212 channels[0].code[1])
2213 except ValueError:
2214 pass
2216 try:
2217 prio_band = priority_band_code.index(
2218 channels[0].code[0])
2219 except ValueError:
2220 pass
2222 if target_sample_rate is None:
2223 rate = -rate
2225 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2226 prio_inst, bic))
2228 useful_bics.sort()
2230 for _, _, rate, _, _, bic in useful_bics:
2231 channels = sorted(
2232 bic_to_channels[bic],
2233 key=lambda channel: channel.code)
2235 if channels:
2236 for channel in channels:
2237 nslcs[nsl + (channel.code,)] = channel
2239 break
2241 return nslcs
2243 def get_pyrocko_response(
2244 self, nslc,
2245 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2247 net, sta, loc, cha = nslc
2248 resps = []
2249 for _, _, channel in self.iter_network_station_channels(
2250 net, sta, loc, cha, time=time, timespan=timespan):
2251 resp = channel.response
2252 if resp:
2253 resp.check_sample_rates(channel)
2254 resp.check_units()
2255 resps.append(resp.get_pyrocko_response(
2256 '.'.join(nslc),
2257 fake_input_units=fake_input_units,
2258 stages=stages).expect_one())
2260 if not resps:
2261 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2262 elif len(resps) > 1:
2263 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2265 return resps[0]
2267 @property
2268 def n_code_list(self):
2269 return sorted(set(x.code for x in self.network_list))
2271 @property
2272 def ns_code_list(self):
2273 nss = set()
2274 for network in self.network_list:
2275 for station in network.station_list:
2276 nss.add((network.code, station.code))
2278 return sorted(nss)
2280 @property
2281 def nsl_code_list(self):
2282 nsls = set()
2283 for network in self.network_list:
2284 for station in network.station_list:
2285 for channel in station.channel_list:
2286 nsls.add(
2287 (network.code, station.code, channel.location_code))
2289 return sorted(nsls)
2291 @property
2292 def nslc_code_list(self):
2293 nslcs = set()
2294 for network in self.network_list:
2295 for station in network.station_list:
2296 for channel in station.channel_list:
2297 nslcs.add(
2298 (network.code, station.code, channel.location_code,
2299 channel.code))
2301 return sorted(nslcs)
2303 def summary(self):
2304 lst = [
2305 'number of n codes: %i' % len(self.n_code_list),
2306 'number of ns codes: %i' % len(self.ns_code_list),
2307 'number of nsl codes: %i' % len(self.nsl_code_list),
2308 'number of nslc codes: %i' % len(self.nslc_code_list)
2309 ]
2310 return '\n'.join(lst)
2312 def summary_stages(self):
2313 data = []
2314 for network, station, channel in self.iter_network_station_channels():
2315 nslc = (network.code, station.code, channel.location_code,
2316 channel.code)
2318 stages = []
2319 in_units = '?'
2320 out_units = '?'
2321 if channel.response:
2322 sens = channel.response.instrument_sensitivity
2323 if sens:
2324 in_units = sens.input_units.name.upper()
2325 out_units = sens.output_units.name.upper()
2327 for stage in channel.response.stage_list:
2328 stages.append(stage.summary())
2330 data.append(
2331 (nslc, tts(channel.start_date), tts(channel.end_date),
2332 in_units, out_units, stages))
2334 data.sort()
2336 lst = []
2337 for nslc, stmin, stmax, in_units, out_units, stages in data:
2338 lst.append(' %s: %s - %s, %s -> %s' % (
2339 '.'.join(nslc), stmin, stmax, in_units, out_units))
2340 for stage in stages:
2341 lst.append(' %s' % stage)
2343 return '\n'.join(lst)
2345 def _check_overlaps(self):
2346 by_nslc = {}
2347 for network in self.network_list:
2348 for station in network.station_list:
2349 for channel in station.channel_list:
2350 nslc = (network.code, station.code, channel.location_code,
2351 channel.code)
2352 if nslc not in by_nslc:
2353 by_nslc[nslc] = []
2355 by_nslc[nslc].append(channel)
2357 errors = []
2358 for nslc, channels in by_nslc.items():
2359 errors.extend(check_overlaps('Channel', nslc, channels))
2361 return errors
2363 def check(self):
2364 errors = []
2365 for meth in [self._check_overlaps]:
2366 errors.extend(meth())
2368 if errors:
2369 raise Inconsistencies(
2370 'Inconsistencies found in StationXML:\n '
2371 + '\n '.join(errors))
2374def load_channel_table(stream):
2376 networks = {}
2377 stations = {}
2379 for line in stream:
2380 line = str(line.decode('ascii'))
2381 if line.startswith('#'):
2382 continue
2384 t = line.rstrip().split('|')
2386 if len(t) != 17:
2387 logger.warning('Invalid channel record: %s' % line)
2388 continue
2390 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2391 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2393 try:
2394 scale = float(scale)
2395 except ValueError:
2396 scale = None
2398 try:
2399 scale_freq = float(scale_freq)
2400 except ValueError:
2401 scale_freq = None
2403 try:
2404 depth = float(dep)
2405 except ValueError:
2406 depth = 0.0
2408 try:
2409 azi = float(azi)
2410 dip = float(dip)
2411 except ValueError:
2412 azi = None
2413 dip = None
2415 try:
2416 if net not in networks:
2417 network = Network(code=net)
2418 else:
2419 network = networks[net]
2421 if (net, sta) not in stations:
2422 station = Station(
2423 code=sta, latitude=lat, longitude=lon, elevation=ele)
2425 station.regularize()
2426 else:
2427 station = stations[net, sta]
2429 if scale:
2430 resp = Response(
2431 instrument_sensitivity=Sensitivity(
2432 value=scale,
2433 frequency=scale_freq,
2434 input_units=scale_units))
2435 else:
2436 resp = None
2438 channel = Channel(
2439 code=cha,
2440 location_code=loc.strip(),
2441 latitude=lat,
2442 longitude=lon,
2443 elevation=ele,
2444 depth=depth,
2445 azimuth=azi,
2446 dip=dip,
2447 sensor=Equipment(description=sens),
2448 response=resp,
2449 sample_rate=sample_rate,
2450 start_date=start_date,
2451 end_date=end_date or None)
2453 channel.regularize()
2455 except ValidationError:
2456 raise InvalidRecord(line)
2458 if net not in networks:
2459 networks[net] = network
2461 if (net, sta) not in stations:
2462 stations[net, sta] = station
2463 network.station_list.append(station)
2465 station.channel_list.append(channel)
2467 return FDSNStationXML(
2468 source='created from table input',
2469 created=time.time(),
2470 network_list=sorted(networks.values(), key=lambda x: x.code))
2473def primitive_merge(sxs):
2474 networks = []
2475 for sx in sxs:
2476 networks.extend(sx.network_list)
2478 return FDSNStationXML(
2479 source='merged from different sources',
2480 created=time.time(),
2481 network_list=copy.deepcopy(
2482 sorted(networks, key=lambda x: x.code)))
2485def split_channels(sx):
2486 for nslc in sx.nslc_code_list:
2487 network_list = sx.network_list
2488 network_list_filtered = [
2489 network for network in network_list
2490 if network.code == nslc[0]]
2492 for network in network_list_filtered:
2493 sx.network_list = [network]
2494 station_list = network.station_list
2495 station_list_filtered = [
2496 station for station in station_list
2497 if station.code == nslc[1]]
2499 for station in station_list_filtered:
2500 network.station_list = [station]
2501 channel_list = station.channel_list
2502 station.channel_list = [
2503 channel for channel in channel_list
2504 if (channel.location_code, channel.code) == nslc[2:4]]
2506 yield nslc, copy.deepcopy(sx)
2508 station.channel_list = channel_list
2510 network.station_list = station_list
2512 sx.network_list = network_list
2515if __name__ == '__main__':
2516 from optparse import OptionParser
2518 util.setup_logging('pyrocko.io.stationxml', 'warning')
2520 usage = \
2521 'python -m pyrocko.io.stationxml check|stats|stages ' \
2522 '<filename> [options]'
2524 description = '''Torture StationXML file.'''
2526 parser = OptionParser(
2527 usage=usage,
2528 description=description,
2529 formatter=util.BetterHelpFormatter())
2531 (options, args) = parser.parse_args(sys.argv[1:])
2533 if len(args) != 2:
2534 parser.print_help()
2535 sys.exit(1)
2537 action, path = args
2539 sx = load_xml(filename=path)
2540 if action == 'check':
2541 try:
2542 sx.check()
2543 except Inconsistencies as e:
2544 logger.error(e)
2545 sys.exit(1)
2547 elif action == 'stats':
2548 print(sx.summary())
2550 elif action == 'stages':
2551 print(sx.summary_stages())
2553 else:
2554 parser.print_help()
2555 sys.exit('unknown action: %s' % action)