Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/stationxml.py: 70%
1187 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-03-07 11:54 +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,
47 ('RAD', 'RAD'): None,
48 ('RAD/S', 'RAD'): response.IntegrationResponse(1),
49 ('RAD/S**2', 'RAD'): response.IntegrationResponse(2),
50 ('RAD', 'RAD/S'): response.DifferentiationResponse(1),
51 ('RAD/S', 'RAD/S'): None,
52 ('RAD/S**2', 'RAD/S'): response.IntegrationResponse(1),
53 ('RAD', 'RAD/S**2'): response.DifferentiationResponse(2),
54 ('RAD/S', 'RAD/S**2'): response.DifferentiationResponse(1),
55 ('RAD/S**2', 'RAD/S**2'): None}
58units_to_quantity = {
59 'M/S': 'velocity',
60 'M': 'displacement',
61 'M/S**2': 'acceleration',
62 'V': 'voltage',
63 'COUNT': 'counts',
64 'PA': 'pressure',
65 'RAD': 'rotation_displacement',
66 'RAD/S': 'rotation_velocity',
67 'RAD/S**2': 'rotation_acceleration'}
70quantity_to_units = dict((v, k) for (k, v) in units_to_quantity.items())
73units_fixes = {
74 'R': 'RAD',
75 'R/S': 'RAD/S',
76 'R/S**2': 'RAD/S**2',
77 'COUNTS': 'COUNT'}
80def sanitize_units(s):
81 s = s.upper()
82 return units_fixes.get(s, s)
85def to_quantity(units, context, delivery):
87 if units is None:
88 return None
90 name = sanitize_units(units.name)
91 if name in units_to_quantity:
92 return units_to_quantity[name]
93 else:
94 delivery.log.append((
95 'warning',
96 'Units not supported by Squirrel framework: %s' % (
97 name + (
98 ' (%s)' % units.description if units.description else '')),
99 context))
101 return 'unsupported_quantity(%s)' % units
104class StationXMLError(Exception):
105 pass
108class Inconsistencies(StationXMLError):
109 pass
112class NoResponseInformation(StationXMLError):
113 pass
116class MultipleResponseInformation(StationXMLError):
117 pass
120class InconsistentResponseInformation(StationXMLError):
121 pass
124class InconsistentChannelLocations(StationXMLError):
125 pass
128class InvalidRecord(StationXMLError):
129 def __init__(self, line):
130 StationXMLError.__init__(self)
131 self._line = line
133 def __str__(self):
134 return 'Invalid record: "%s"' % self._line
137_exceptions = dict(
138 Inconsistencies=Inconsistencies,
139 NoResponseInformation=NoResponseInformation,
140 MultipleResponseInformation=MultipleResponseInformation,
141 InconsistentResponseInformation=InconsistentResponseInformation,
142 InconsistentChannelLocations=InconsistentChannelLocations,
143 InvalidRecord=InvalidRecord,
144 ValueError=ValueError)
147_logs = dict(
148 info=logger.info,
149 warning=logger.warning,
150 error=logger.error)
153class DeliveryError(StationXMLError):
154 pass
157class Delivery(Object):
159 def __init__(self, payload=None, log=None, errors=None, error=None):
160 if payload is None:
161 payload = []
163 if log is None:
164 log = []
166 if errors is None:
167 errors = []
169 if error is not None:
170 errors.append(error)
172 Object.__init__(self, payload=payload, log=log, errors=errors)
174 payload = List.T(Any.T())
175 log = List.T(Tuple.T(3, String.T()))
176 errors = List.T(Tuple.T(3, String.T()))
178 def extend(self, other):
179 self.payload.extend(other.payload)
180 self.log.extend(other.log)
181 self.errors.extend(other.errors)
183 def extend_without_payload(self, other):
184 self.log.extend(other.log)
185 self.errors.extend(other.errors)
186 return other.payload
188 def emit_log(self):
189 for name, message, context in self.log:
190 message = '%s: %s' % (context, message)
191 _logs[name](message)
193 def expect(self, quiet=False):
194 if not quiet:
195 self.emit_log()
197 if self.errors:
198 name, message, context = self.errors[0]
199 if context:
200 message += ' (%s)' % context
202 if len(self.errors) > 1:
203 message += ' Additional errors pending.'
205 raise _exceptions[name](message)
207 return self.payload
209 def expect_one(self, quiet=False):
210 payload = self.expect(quiet=quiet)
211 if len(payload) != 1:
212 raise DeliveryError(
213 'Expected 1 element but got %i.' % len(payload))
215 return payload[0]
218def wrap(s, width=80, indent=4):
219 words = s.split()
220 lines = []
221 t = []
222 n = 0
223 for w in words:
224 if n + len(w) >= width:
225 lines.append(' '.join(t))
226 n = indent
227 t = [' '*(indent-1)]
229 t.append(w)
230 n += len(w) + 1
232 lines.append(' '.join(t))
233 return '\n'.join(lines)
236def same(x, eps=0.0):
237 if any(type(x[0]) is not type(r) for r in x):
238 return False
240 if isinstance(x[0], float):
241 return all(abs(r-x[0]) <= eps for r in x)
242 else:
243 return all(r == x[0] for r in x)
246def same_sample_rate(a, b, eps=1.0e-6):
247 return abs(a - b) < min(a, b)*eps
250def evaluate1(resp, f):
251 return resp.evaluate(num.array([f], dtype=float))[0]
254def check_resp(resp, value, frequency, limit_db, prelude, context):
256 try:
257 value_resp = num.abs(evaluate1(resp, frequency))
258 except response.InvalidResponseError as e:
259 return Delivery(log=[(
260 'warning',
261 'Could not check response: %s' % str(e),
262 context)])
264 if value_resp == 0.0:
265 return Delivery(log=[(
266 'warning',
267 '%s\n'
268 ' computed response is zero' % prelude,
269 context)])
271 diff_db = 20.0 * num.log10(value_resp/value)
273 if num.abs(diff_db) > limit_db:
274 return Delivery(log=[(
275 'warning',
276 '%s\n'
277 ' reported value: %g\n'
278 ' computed value: %g\n'
279 ' at frequency [Hz]: %g\n'
280 ' factor: %g\n'
281 ' difference [dB]: %g\n'
282 ' limit [dB]: %g' % (
283 prelude,
284 value,
285 value_resp,
286 frequency,
287 value_resp/value,
288 diff_db,
289 limit_db),
290 context)])
292 return Delivery()
295def tts(t):
296 if t is None:
297 return '<none>'
298 else:
299 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
302def le_open_left(a, b):
303 return a is None or (b is not None and a <= b)
306def le_open_right(a, b):
307 return b is None or (a is not None and a <= b)
310def eq_open(a, b):
311 return (a is None and b is None) \
312 or (a is not None and b is not None and a == b)
315def contains(a, b):
316 return le_open_left(a.start_date, b.start_date) \
317 and le_open_right(b.end_date, a.end_date)
320def find_containing(candidates, node):
321 for candidate in candidates:
322 if contains(candidate, node):
323 return candidate
325 return None
328this_year = time.gmtime()[0]
331class DummyAwareOptionalTimestamp(Object):
332 '''
333 Optional timestamp with support for some common placeholder values.
335 Some StationXML files contain arbitrary placeholder values for open end
336 intervals, like "2100-01-01". Depending on the time range supported by the
337 system, these dates are translated into ``None`` to prevent crashes with
338 this type.
339 '''
340 dummy_for = (hpfloat, float)
341 dummy_for_description = 'pyrocko.util.get_time_float'
343 class __T(TBase):
345 def regularize_extra(self, val):
346 time_float = get_time_float()
348 if isinstance(val, datetime.datetime):
349 tt = val.utctimetuple()
350 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
352 elif isinstance(val, datetime.date):
353 tt = val.timetuple()
354 val = time_float(calendar.timegm(tt))
356 elif isinstance(val, str):
357 val = val.strip()
359 tz_offset = 0
361 m = re_tz.search(val)
362 if m:
363 sh = m.group(2)
364 sm = m.group(4)
365 tz_offset = (int(sh)*3600 if sh else 0) \
366 + (int(sm)*60 if sm else 0)
368 val = re_tz.sub('', val)
370 if len(val) > 10 and val[10] == 'T':
371 val = val.replace('T', ' ', 1)
373 try:
374 val = util.str_to_time(val) - tz_offset
376 except util.TimeStrError:
377 year = int(val[:4])
378 if sys.maxsize > 2**32: # if we're on 64bit
379 if year > this_year + 100:
380 return None # StationXML contained a dummy date
382 if year < 1903: # for macOS, 1900-01-01 dummy dates
383 return None
385 else: # 32bit end of time is in 2038
386 if this_year < 2037 and year > 2037 or year < 1903:
387 return None # StationXML contained a dummy date
389 raise
391 elif isinstance(val, (int, float)):
392 val = time_float(val)
394 else:
395 raise ValidationError(
396 '%s: cannot convert "%s" to type %s' % (
397 self.xname(), val, time_float))
399 return val
401 def to_save(self, val):
402 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
403 .rstrip('0').rstrip('.')
405 def to_save_xml(self, val):
406 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
407 .rstrip('0').rstrip('.') + 'Z'
410class Nominal(StringChoice):
411 choices = [
412 'NOMINAL',
413 'CALCULATED']
416class Email(UnicodePattern):
417 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
420class RestrictedStatus(StringChoice):
421 choices = [
422 'open',
423 'closed',
424 'partial']
427class Type(StringChoice):
428 choices = [
429 'TRIGGERED',
430 'CONTINUOUS',
431 'HEALTH',
432 'GEOPHYSICAL',
433 'WEATHER',
434 'FLAG',
435 'SYNTHESIZED',
436 'INPUT',
437 'EXPERIMENTAL',
438 'MAINTENANCE',
439 'BEAM']
441 class __T(StringChoice.T):
442 def validate_extra(self, val):
443 if val not in self.choices:
444 logger.warning(
445 'channel type: "%s" is not a valid choice out of %s' %
446 (val, repr(self.choices)))
449class PzTransferFunction(StringChoice):
450 choices = [
451 'LAPLACE (RADIANS/SECOND)',
452 'LAPLACE (HERTZ)',
453 'DIGITAL (Z-TRANSFORM)']
456class Symmetry(StringChoice):
457 choices = [
458 'NONE',
459 'EVEN',
460 'ODD']
463class CfTransferFunction(StringChoice):
465 class __T(StringChoice.T):
466 def validate(self, val, regularize=False, depth=-1):
467 if regularize:
468 try:
469 val = str(val)
470 except ValueError:
471 raise ValidationError(
472 '%s: cannot convert to string %s' % (self.xname,
473 repr(val)))
475 val = self._dummy_cls.replacements.get(val, val)
477 self.validate_extra(val)
478 return val
480 choices = [
481 'ANALOG (RADIANS/SECOND)',
482 'ANALOG (HERTZ)',
483 'DIGITAL']
485 replacements = {
486 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
487 'ANALOG (HZ)': 'ANALOG (HERTZ)',
488 }
491class Approximation(StringChoice):
492 choices = [
493 'MACLAURIN']
496class PhoneNumber(StringPattern):
497 pattern = '[0-9]+-[0-9]+'
500class Site(Object):
501 '''
502 Description of a site location using name and optional geopolitical
503 boundaries (country, city, etc.).
504 '''
506 name = Unicode.T(default='', xmltagname='Name')
507 description = Unicode.T(optional=True, xmltagname='Description')
508 town = Unicode.T(optional=True, xmltagname='Town')
509 county = Unicode.T(optional=True, xmltagname='County')
510 region = Unicode.T(optional=True, xmltagname='Region')
511 country = Unicode.T(optional=True, xmltagname='Country')
514class ExternalReference(Object):
515 '''
516 This type contains a URI and description for external data that users may
517 want to reference in StationXML.
518 '''
520 uri = String.T(xmltagname='URI')
521 description = Unicode.T(xmltagname='Description')
524class Units(Object):
525 '''
526 A type to document units. Corresponds to SEED blockette 34.
527 '''
529 def __init__(self, name=None, **kwargs):
530 Object.__init__(self, name=name, **kwargs)
532 name = String.T(xmltagname='Name')
533 description = Unicode.T(optional=True, xmltagname='Description')
536class Counter(Int):
537 pass
540class SampleRateRatio(Object):
541 '''
542 Sample rate expressed as number of samples in a number of seconds.
543 '''
545 number_samples = Int.T(xmltagname='NumberSamples')
546 number_seconds = Int.T(xmltagname='NumberSeconds')
549class Gain(Object):
550 '''
551 Complex type for sensitivity and frequency ranges. This complex type can be
552 used to represent both overall sensitivities and individual stage gains.
553 The FrequencyRangeGroup is an optional construct that defines a pass band
554 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
555 is valid within the number of decibels specified in FrequencyDBVariation.
556 '''
558 def __init__(self, value=None, **kwargs):
559 Object.__init__(self, value=value, **kwargs)
561 value = Float.T(optional=True, xmltagname='Value')
562 frequency = Float.T(optional=True, xmltagname='Frequency')
564 def summary(self):
565 return 'gain(%g @ %g)' % (self.value, self.frequency)
568class NumeratorCoefficient(Object):
569 i = Int.T(optional=True, xmlstyle='attribute')
570 value = Float.T(xmlstyle='content')
573class FloatNoUnit(Object):
574 def __init__(self, value=None, **kwargs):
575 Object.__init__(self, value=value, **kwargs)
577 plus_error = Float.T(optional=True, xmlstyle='attribute')
578 minus_error = Float.T(optional=True, xmlstyle='attribute')
579 value = Float.T(xmlstyle='content')
582class FloatWithUnit(FloatNoUnit):
583 unit = String.T(optional=True, xmlstyle='attribute')
586class Equipment(Object):
587 resource_id = String.T(optional=True, xmlstyle='attribute')
588 type = String.T(optional=True, xmltagname='Type')
589 description = Unicode.T(optional=True, xmltagname='Description')
590 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
591 vendor = Unicode.T(optional=True, xmltagname='Vendor')
592 model = Unicode.T(optional=True, xmltagname='Model')
593 serial_number = String.T(optional=True, xmltagname='SerialNumber')
594 installation_date = DummyAwareOptionalTimestamp.T(
595 optional=True,
596 xmltagname='InstallationDate')
597 removal_date = DummyAwareOptionalTimestamp.T(
598 optional=True,
599 xmltagname='RemovalDate')
600 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
603class PhoneNumber(Object):
604 description = Unicode.T(optional=True, xmlstyle='attribute')
605 country_code = Int.T(optional=True, xmltagname='CountryCode')
606 area_code = Int.T(xmltagname='AreaCode')
607 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
610class BaseFilter(Object):
611 '''
612 The BaseFilter is derived by all filters.
613 '''
615 resource_id = String.T(optional=True, xmlstyle='attribute')
616 name = String.T(optional=True, xmlstyle='attribute')
617 description = Unicode.T(optional=True, xmltagname='Description')
618 input_units = Units.T(optional=True, xmltagname='InputUnits')
619 output_units = Units.T(optional=True, xmltagname='OutputUnits')
622class Sensitivity(Gain):
623 '''
624 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
625 construct that defines a pass band in Hertz (FrequencyStart and
626 FrequencyEnd) in which the SensitivityValue is valid within the number of
627 decibels specified in FrequencyDBVariation.
628 '''
630 input_units = Units.T(optional=True, xmltagname='InputUnits')
631 output_units = Units.T(optional=True, xmltagname='OutputUnits')
632 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
633 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
634 frequency_db_variation = Float.T(optional=True,
635 xmltagname='FrequencyDBVariation')
637 def get_pyrocko_response(self):
638 return Delivery(
639 [response.PoleZeroResponse(constant=self.value)])
642class Coefficient(FloatNoUnit):
643 number = Counter.T(optional=True, xmlstyle='attribute')
646class PoleZero(Object):
647 '''
648 Complex numbers used as poles or zeros in channel response.
649 '''
651 number = Int.T(optional=True, xmlstyle='attribute')
652 real = FloatNoUnit.T(xmltagname='Real')
653 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
655 def value(self):
656 return self.real.value + 1J * self.imaginary.value
659class ClockDrift(FloatWithUnit):
660 unit = String.T(default='SECONDS/SAMPLE', optional=True,
661 xmlstyle='attribute') # fixed
664class Second(FloatWithUnit):
665 '''
666 A time value in seconds.
667 '''
669 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
670 # fixed unit
673class Voltage(FloatWithUnit):
674 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
675 # fixed unit
678class Angle(FloatWithUnit):
679 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
680 # fixed unit
683class Azimuth(FloatWithUnit):
684 '''
685 Instrument azimuth, degrees clockwise from North.
686 '''
688 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
689 # fixed unit
692class Dip(FloatWithUnit):
693 '''
694 Instrument dip in degrees down from horizontal. Together azimuth and dip
695 describe the direction of the sensitive axis of the instrument.
696 '''
698 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
699 # fixed unit
702class Distance(FloatWithUnit):
703 '''
704 Extension of FloatWithUnit for distances, elevations, and depths.
705 '''
707 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
708 # NOT fixed unit!
711class Frequency(FloatWithUnit):
712 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
713 # fixed unit
716class SampleRate(FloatWithUnit):
717 '''
718 Sample rate in samples per second.
719 '''
721 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
722 # fixed unit
725class Person(Object):
726 '''
727 Representation of a person's contact information. A person can belong to
728 multiple agencies and have multiple email addresses and phone numbers.
729 '''
731 name_list = List.T(Unicode.T(xmltagname='Name'))
732 agency_list = List.T(Unicode.T(xmltagname='Agency'))
733 email_list = List.T(Email.T(xmltagname='Email'))
734 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
737class FIR(BaseFilter):
738 '''
739 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
740 also commonly documented using the Coefficients element.
741 '''
743 symmetry = Symmetry.T(xmltagname='Symmetry')
744 numerator_coefficient_list = List.T(
745 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
747 def summary(self):
748 return 'fir(%i%s)' % (
749 self.get_ncoefficiencs(),
750 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
752 def get_effective_coefficients(self):
753 b = num.array(
754 [v.value for v in self.numerator_coefficient_list],
755 dtype=float)
757 if self.symmetry == 'ODD':
758 b = num.concatenate((b, b[-2::-1]))
759 elif self.symmetry == 'EVEN':
760 b = num.concatenate((b, b[::-1]))
762 return b
764 def get_effective_symmetry(self):
765 if self.symmetry == 'NONE':
766 b = self.get_effective_coefficients()
767 if num.all(b - b[::-1] == 0):
768 return ['EVEN', 'ODD'][b.size % 2]
770 return self.symmetry
772 def get_ncoefficiencs(self):
773 nf = len(self.numerator_coefficient_list)
774 if self.symmetry == 'ODD':
775 nc = nf * 2 + 1
776 elif self.symmetry == 'EVEN':
777 nc = nf * 2
778 else:
779 nc = nf
781 return nc
783 def estimate_delay(self, deltat):
784 nc = self.get_ncoefficiencs()
785 if nc > 0:
786 return deltat * (nc - 1) / 2.0
787 else:
788 return 0.0
790 def get_pyrocko_response(
791 self, context, deltat, delay_responses, normalization_frequency):
793 context += self.summary()
795 if not self.numerator_coefficient_list:
796 return Delivery([])
798 b = self.get_effective_coefficients()
800 log = []
801 drop_phase = self.get_effective_symmetry() != 'NONE'
803 if not deltat:
804 log.append((
805 'error',
806 'Digital response requires knowledge about sampling '
807 'interval. Response will be unusable.',
808 context))
810 resp = response.DigitalFilterResponse(
811 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
813 if normalization_frequency is not None and deltat is not None:
814 normalization_frequency = 0.0
815 normalization = num.abs(evaluate1(resp, normalization_frequency))
817 if num.abs(normalization - 1.0) > 1e-2:
818 log.append((
819 'warning',
820 'FIR filter coefficients are not normalized. Normalizing '
821 'them. Factor: %g' % normalization, context))
823 resp = response.DigitalFilterResponse(
824 (b/normalization).tolist(), [1.0], deltat,
825 drop_phase=drop_phase)
827 resps = [resp]
829 if not drop_phase:
830 resps.extend(delay_responses)
832 return Delivery(resps, log=log)
835class Coefficients(BaseFilter):
836 '''
837 Response: coefficients for FIR filter. Laplace transforms or IIR filters
838 can be expressed using type as well but the PolesAndZeros should be used
839 instead. Corresponds to SEED blockette 54.
840 '''
842 cf_transfer_function_type = CfTransferFunction.T(
843 xmltagname='CfTransferFunctionType')
844 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
845 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
847 def summary(self):
848 return 'coef_%s(%i,%i%s)' % (
849 'ABC?'[
850 CfTransferFunction.choices.index(
851 self.cf_transfer_function_type)],
852 len(self.numerator_list),
853 len(self.denominator_list),
854 ',sym' if self.is_symmetric_fir else '')
856 def estimate_delay(self, deltat):
857 nc = len(self.numerator_list)
858 if nc > 0:
859 return deltat * (len(self.numerator_list) - 1) / 2.0
860 else:
861 return 0.0
863 def is_symmetric_fir(self):
864 if len(self.denominator_list) != 0:
865 return False
866 b = [v.value for v in self.numerator_list]
867 return b == b[::-1]
869 def get_pyrocko_response(
870 self, context, deltat, delay_responses, normalization_frequency):
872 context += self.summary()
874 factor = 1.0
875 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
876 factor = 2.0*math.pi
878 if not self.numerator_list and not self.denominator_list:
879 return Delivery(payload=[])
881 b = num.array(
882 [v.value*factor for v in self.numerator_list], dtype=float)
884 a = num.array(
885 [1.0] + [v.value*factor for v in self.denominator_list],
886 dtype=float)
888 log = []
889 resps = []
890 if self.cf_transfer_function_type in [
891 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
892 resps.append(response.AnalogFilterResponse(b, a))
894 elif self.cf_transfer_function_type == 'DIGITAL':
895 if not deltat:
896 log.append((
897 'error',
898 'Digital response requires knowledge about sampling '
899 'interval. Response will be unusable.',
900 context))
902 drop_phase = self.is_symmetric_fir()
903 resp = response.DigitalFilterResponse(
904 b, a, deltat or 0.0, drop_phase=drop_phase)
906 if normalization_frequency is not None and deltat is not None:
907 normalization = num.abs(evaluate1(
908 resp, normalization_frequency))
910 if num.abs(normalization - 1.0) > 1e-2:
911 log.append((
912 'warning',
913 'FIR filter coefficients are not normalized. '
914 'Normalizing them. Factor: %g' % normalization,
915 context))
917 resp = response.DigitalFilterResponse(
918 (b/normalization).tolist(), [1.0], deltat,
919 drop_phase=drop_phase)
921 resps.append(resp)
923 if not drop_phase:
924 resps.extend(delay_responses)
926 else:
927 return Delivery(error=(
928 'ValueError',
929 'Unknown transfer function type: %s' % (
930 self.cf_transfer_function_type)))
932 return Delivery(payload=resps, log=log)
935class Latitude(FloatWithUnit):
936 '''
937 Type for latitude coordinate.
938 '''
940 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
941 # fixed unit
942 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
945class Longitude(FloatWithUnit):
946 '''
947 Type for longitude coordinate.
948 '''
950 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
951 # fixed unit
952 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
955class PolesZeros(BaseFilter):
956 '''
957 Response: complex poles and zeros. Corresponds to SEED blockette 53.
958 '''
960 pz_transfer_function_type = PzTransferFunction.T(
961 xmltagname='PzTransferFunctionType')
962 normalization_factor = Float.T(
963 default=1.0,
964 xmltagname='NormalizationFactor')
965 normalization_frequency = Frequency.T(
966 optional=True, # but required by standard
967 xmltagname='NormalizationFrequency')
968 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
969 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
971 def summary(self):
972 return 'pz_%s(%i,%i)' % (
973 'ABC?'[
974 PzTransferFunction.choices.index(
975 self.pz_transfer_function_type)],
976 len(self.pole_list),
977 len(self.zero_list))
979 def get_pyrocko_response(self, context='', deltat=None):
981 context += self.summary()
983 factor = 1.0
984 cfactor = 1.0
985 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
986 factor = 2. * math.pi
987 cfactor = (2. * math.pi)**(
988 len(self.pole_list) - len(self.zero_list))
990 log = []
991 if self.normalization_factor is None \
992 or self.normalization_factor == 0.0:
994 log.append((
995 'warning',
996 'No pole-zero normalization factor given. '
997 'Assuming a value of 1.0',
998 context))
1000 nfactor = 1.0
1001 else:
1002 nfactor = self.normalization_factor
1004 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
1005 if not is_digital:
1006 resp = response.PoleZeroResponse(
1007 constant=nfactor*cfactor,
1008 zeros=[z.value()*factor for z in self.zero_list],
1009 poles=[p.value()*factor for p in self.pole_list])
1010 else:
1011 if not deltat:
1012 log.append((
1013 'error',
1014 'Digital response requires knowledge about sampling '
1015 'interval. Response will be unusable.',
1016 context))
1018 resp = response.DigitalPoleZeroResponse(
1019 constant=nfactor*cfactor,
1020 zeros=[z.value()*factor for z in self.zero_list],
1021 poles=[p.value()*factor for p in self.pole_list],
1022 deltat=deltat or 0.0)
1024 if not self.normalization_frequency:
1025 log.append((
1026 'warning',
1027 'Cannot check pole-zero normalization factor: '
1028 'No normalization frequency given.',
1029 context))
1031 else:
1032 if is_digital and not deltat:
1033 log.append((
1034 'warning',
1035 'Cannot check computed vs reported normalization '
1036 'factor without knowing the sampling interval.',
1037 context))
1038 else:
1039 computed_normalization_factor = nfactor / abs(evaluate1(
1040 resp, self.normalization_frequency.value))
1042 db = 20.0 * num.log10(
1043 computed_normalization_factor / nfactor)
1045 if abs(db) > 0.17:
1046 log.append((
1047 'warning',
1048 'Computed and reported normalization factors differ '
1049 'by %g dB: computed: %g, reported: %g' % (
1050 db,
1051 computed_normalization_factor,
1052 nfactor),
1053 context))
1055 return Delivery([resp], log)
1058class ResponseListElement(Object):
1059 frequency = Frequency.T(xmltagname='Frequency')
1060 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1061 phase = Angle.T(xmltagname='Phase')
1064class Polynomial(BaseFilter):
1065 '''
1066 Response: expressed as a polynomial (allows non-linear sensors to be
1067 described). Corresponds to SEED blockette 62. Can be used to describe a
1068 stage of acquisition or a complete system.
1069 '''
1071 approximation_type = Approximation.T(default='MACLAURIN',
1072 xmltagname='ApproximationType')
1073 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1074 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1075 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1076 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1077 maximum_error = Float.T(xmltagname='MaximumError')
1078 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1080 def summary(self):
1081 return 'poly(%i)' % len(self.coefficient_list)
1084class Decimation(Object):
1085 '''
1086 Corresponds to SEED blockette 57.
1087 '''
1089 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1090 factor = Int.T(xmltagname='Factor')
1091 offset = Int.T(xmltagname='Offset')
1092 delay = FloatWithUnit.T(xmltagname='Delay')
1093 correction = FloatWithUnit.T(xmltagname='Correction')
1095 def summary(self):
1096 return 'deci(%i, %g -> %g, %g)' % (
1097 self.factor,
1098 self.input_sample_rate.value,
1099 self.input_sample_rate.value / self.factor,
1100 self.delay.value)
1102 def get_pyrocko_response(self):
1103 if self.delay and self.delay.value != 0.0:
1104 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1105 else:
1106 return Delivery([])
1109class Operator(Object):
1110 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1111 contact_list = List.T(Person.T(xmltagname='Contact'))
1112 web_site = String.T(optional=True, xmltagname='WebSite')
1115class Comment(Object):
1116 '''
1117 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1118 and 59.
1119 '''
1121 id = Counter.T(optional=True, xmlstyle='attribute')
1122 value = Unicode.T(xmltagname='Value')
1123 begin_effective_time = DummyAwareOptionalTimestamp.T(
1124 optional=True,
1125 xmltagname='BeginEffectiveTime')
1126 end_effective_time = DummyAwareOptionalTimestamp.T(
1127 optional=True,
1128 xmltagname='EndEffectiveTime')
1129 author_list = List.T(Person.T(xmltagname='Author'))
1132class ResponseList(BaseFilter):
1133 '''
1134 Response: list of frequency, amplitude and phase values. Corresponds to
1135 SEED blockette 55.
1136 '''
1138 response_list_element_list = List.T(
1139 ResponseListElement.T(xmltagname='ResponseListElement'))
1141 def summary(self):
1142 return 'list(%i)' % len(self.response_list_element_list)
1145class Log(Object):
1146 '''
1147 Container for log entries.
1148 '''
1150 entry_list = List.T(Comment.T(xmltagname='Entry'))
1153class ResponseStage(Object):
1154 '''
1155 This complex type represents channel response and covers SEED blockettes 53
1156 to 56.
1157 '''
1159 number = Counter.T(xmlstyle='attribute')
1160 resource_id = String.T(optional=True, xmlstyle='attribute')
1161 poles_zeros_list = List.T(
1162 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1163 coefficients_list = List.T(
1164 Coefficients.T(optional=True, xmltagname='Coefficients'))
1165 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1166 fir = FIR.T(optional=True, xmltagname='FIR')
1167 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1168 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1169 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1171 def summary(self):
1172 elements = []
1174 for stuff in [
1175 self.poles_zeros_list, self.coefficients_list,
1176 self.response_list, self.fir, self.polynomial,
1177 self.decimation, self.stage_gain]:
1179 if stuff:
1180 if isinstance(stuff, Object):
1181 elements.append(stuff.summary())
1182 else:
1183 elements.extend(obj.summary() for obj in stuff)
1185 return '%i: %s %s -> %s' % (
1186 self.number,
1187 ', '.join(elements),
1188 sanitize_units(self.input_units.name)
1189 if self.input_units else '?',
1190 sanitize_units(self.output_units.name)
1191 if self.output_units else '?')
1193 def get_squirrel_response_stage(self, context):
1194 from pyrocko.squirrel.model import ResponseStage
1196 delivery = Delivery()
1197 delivery_pr = self.get_pyrocko_response(context)
1198 log = delivery_pr.log
1199 delivery_pr.log = []
1200 elements = delivery.extend_without_payload(delivery_pr)
1202 delivery.payload = [ResponseStage(
1203 input_quantity=to_quantity(self.input_units, context, delivery),
1204 output_quantity=to_quantity(self.output_units, context, delivery),
1205 input_sample_rate=self.input_sample_rate,
1206 output_sample_rate=self.output_sample_rate,
1207 elements=elements,
1208 log=log)]
1210 return delivery
1212 def get_pyrocko_response(self, context, gain_only=False):
1214 context = context + ', stage %i' % self.number
1216 responses = []
1217 log = []
1218 if self.stage_gain:
1219 normalization_frequency = self.stage_gain.frequency or 0.0
1220 else:
1221 normalization_frequency = 0.0
1223 if not gain_only:
1224 deltat = None
1225 delay_responses = []
1226 if self.decimation:
1227 rate = self.decimation.input_sample_rate.value
1228 if rate > 0.0:
1229 deltat = 1.0 / rate
1230 delivery = self.decimation.get_pyrocko_response()
1231 if delivery.errors:
1232 return delivery
1234 delay_responses = delivery.payload
1235 log.extend(delivery.log)
1237 for pzs in self.poles_zeros_list:
1238 delivery = pzs.get_pyrocko_response(context, deltat)
1239 if delivery.errors:
1240 return delivery
1242 pz_resps = delivery.payload
1243 log.extend(delivery.log)
1244 responses.extend(pz_resps)
1246 # emulate incorrect? evalresp behaviour
1247 if pzs.normalization_frequency is not None and \
1248 pzs.normalization_frequency.value \
1249 != normalization_frequency \
1250 and normalization_frequency != 0.0:
1252 try:
1253 trial = response.MultiplyResponse(pz_resps)
1254 anorm = num.abs(evaluate1(
1255 trial, pzs.normalization_frequency.value))
1256 asens = num.abs(
1257 evaluate1(trial, normalization_frequency))
1259 factor = anorm/asens
1261 if abs(factor - 1.0) > 0.01:
1262 log.append((
1263 'warning',
1264 'PZ normalization frequency (%g) is different '
1265 'from stage gain frequency (%s) -> Emulating '
1266 'possibly incorrect evalresp behaviour. '
1267 'Correction factor: %g' % (
1268 pzs.normalization_frequency.value,
1269 normalization_frequency,
1270 factor),
1271 context))
1273 responses.append(
1274 response.PoleZeroResponse(constant=factor))
1275 except response.InvalidResponseError as e:
1276 log.append((
1277 'warning',
1278 'Could not check response: %s' % str(e),
1279 context))
1281 if len(self.poles_zeros_list) > 1:
1282 log.append((
1283 'warning',
1284 'Multiple poles and zeros records in single response '
1285 'stage.',
1286 context))
1288 for cfs in self.coefficients_list + (
1289 [self.fir] if self.fir else []):
1291 delivery = cfs.get_pyrocko_response(
1292 context, deltat, delay_responses,
1293 normalization_frequency)
1295 if delivery.errors:
1296 return delivery
1298 responses.extend(delivery.payload)
1299 log.extend(delivery.log)
1301 if len(self.coefficients_list) > 1:
1302 log.append((
1303 'warning',
1304 'Multiple filter coefficients lists in single response '
1305 'stage.',
1306 context))
1308 if self.response_list:
1309 log.append((
1310 'warning',
1311 'Unhandled response element of type: ResponseList',
1312 context))
1314 if self.polynomial:
1315 log.append((
1316 'warning',
1317 'Unhandled response element of type: Polynomial',
1318 context))
1320 if self.stage_gain:
1321 responses.append(
1322 response.PoleZeroResponse(constant=self.stage_gain.value))
1324 return Delivery(responses, log)
1326 @property
1327 def input_units(self):
1328 for e in (self.poles_zeros_list + self.coefficients_list +
1329 [self.response_list, self.fir, self.polynomial]):
1330 if e is not None:
1331 return e.input_units
1333 return None
1335 @property
1336 def output_units(self):
1337 for e in (self.poles_zeros_list + self.coefficients_list +
1338 [self.response_list, self.fir, self.polynomial]):
1339 if e is not None:
1340 return e.output_units
1342 return None
1344 @property
1345 def input_sample_rate(self):
1346 if self.decimation:
1347 return self.decimation.input_sample_rate.value
1349 return None
1351 @property
1352 def output_sample_rate(self):
1353 if self.decimation:
1354 return self.decimation.input_sample_rate.value \
1355 / self.decimation.factor
1357 return None
1360class Response(Object):
1361 resource_id = String.T(optional=True, xmlstyle='attribute')
1362 instrument_sensitivity = Sensitivity.T(optional=True,
1363 xmltagname='InstrumentSensitivity')
1364 instrument_polynomial = Polynomial.T(optional=True,
1365 xmltagname='InstrumentPolynomial')
1366 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1368 @property
1369 def output_sample_rate(self):
1370 if self.stage_list:
1371 return self.stage_list[-1].output_sample_rate
1372 else:
1373 return None
1375 def check_sample_rates(self, channel):
1377 if self.stage_list:
1378 sample_rate = None
1380 for stage in self.stage_list:
1381 if stage.decimation:
1382 input_sample_rate = \
1383 stage.decimation.input_sample_rate.value
1385 if sample_rate is not None and not same_sample_rate(
1386 sample_rate, input_sample_rate):
1388 logger.warning(
1389 'Response stage %i has unexpected input sample '
1390 'rate: %g Hz (expected: %g Hz)' % (
1391 stage.number,
1392 input_sample_rate,
1393 sample_rate))
1395 sample_rate = input_sample_rate / stage.decimation.factor
1397 if sample_rate is not None and channel.sample_rate \
1398 and not same_sample_rate(
1399 sample_rate, channel.sample_rate.value):
1401 logger.warning(
1402 'Channel sample rate (%g Hz) does not match sample rate '
1403 'deducted from response stages information (%g Hz).' % (
1404 channel.sample_rate.value,
1405 sample_rate))
1407 def check_units(self):
1409 if self.instrument_sensitivity \
1410 and self.instrument_sensitivity.input_units:
1412 units = sanitize_units(
1413 self.instrument_sensitivity.input_units.name)
1415 if self.stage_list:
1416 for stage in self.stage_list:
1417 if units and stage.input_units \
1418 and sanitize_units(stage.input_units.name) != units:
1420 logger.warning(
1421 'Input units of stage %i (%s) do not match %s (%s).'
1422 % (
1423 stage.number,
1424 units,
1425 'output units of stage %i'
1426 if stage.number == 0
1427 else 'sensitivity input units',
1428 units))
1430 if stage.output_units:
1431 units = sanitize_units(stage.output_units.name)
1432 else:
1433 units = None
1435 sout_units = self.instrument_sensitivity.output_units
1436 if self.instrument_sensitivity and sout_units:
1437 if units is not None and units != sanitize_units(
1438 sout_units.name):
1439 logger.warning(
1440 'Output units of stage %i (%s) do not match %s (%s).'
1441 % (
1442 stage.number,
1443 units,
1444 'sensitivity output units',
1445 sanitize_units(sout_units.name)))
1447 def _sensitivity_checkpoints(self, responses, context):
1448 delivery = Delivery()
1450 if self.instrument_sensitivity:
1451 sval = self.instrument_sensitivity.value
1452 sfreq = self.instrument_sensitivity.frequency
1453 if sval is None:
1454 delivery.log.append((
1455 'warning',
1456 'No sensitivity value given.',
1457 context))
1459 elif sval is None:
1460 delivery.log.append((
1461 'warning',
1462 'Reported sensitivity value is zero.',
1463 context))
1465 elif sfreq is None:
1466 delivery.log.append((
1467 'warning',
1468 'Sensitivity frequency not given.',
1469 context))
1471 else:
1472 trial = response.MultiplyResponse(responses)
1474 delivery.extend(
1475 check_resp(
1476 trial, sval, sfreq, 0.1,
1477 'Instrument sensitivity value inconsistent with '
1478 'sensitivity computed from complete response.',
1479 context))
1481 delivery.payload.append(response.FrequencyResponseCheckpoint(
1482 frequency=sfreq, value=sval))
1484 return delivery
1486 def get_squirrel_response(self, context, **kwargs):
1487 from pyrocko.squirrel.model import Response
1489 if self.stage_list:
1490 delivery = Delivery()
1491 for istage, stage in enumerate(self.stage_list):
1492 delivery.extend(stage.get_squirrel_response_stage(context))
1494 checkpoints = []
1495 if not delivery.errors:
1496 all_responses = []
1497 for sq_stage in delivery.payload:
1498 all_responses.extend(sq_stage.elements)
1500 checkpoints.extend(delivery.extend_without_payload(
1501 self._sensitivity_checkpoints(all_responses, context)))
1503 sq_stages = delivery.payload
1504 if sq_stages:
1505 if sq_stages[0].input_quantity is None \
1506 and self.instrument_sensitivity is not None:
1508 sq_stages[0].input_quantity = to_quantity(
1509 self.instrument_sensitivity.input_units,
1510 context, delivery)
1511 sq_stages[-1].output_quantity = to_quantity(
1512 self.instrument_sensitivity.output_units,
1513 context, delivery)
1515 sq_stages = delivery.expect()
1517 return Response(
1518 stages=sq_stages,
1519 log=delivery.log,
1520 checkpoints=checkpoints,
1521 **kwargs)
1523 elif self.instrument_sensitivity:
1524 raise NoResponseInformation(
1525 "Only instrument sensitivity given (won't use it). (%s)."
1526 % context)
1527 else:
1528 raise NoResponseInformation(
1529 'Empty instrument response (%s).'
1530 % context)
1532 def get_pyrocko_response(
1533 self, context, fake_input_units=None, stages=(0, 1)):
1535 if fake_input_units is not None:
1536 fake_input_units = sanitize_units(fake_input_units)
1538 delivery = Delivery()
1539 if self.stage_list:
1540 for istage, stage in enumerate(self.stage_list):
1541 delivery.extend(stage.get_pyrocko_response(
1542 context, gain_only=not (
1543 stages is None or stages[0] <= istage < stages[1])))
1545 elif self.instrument_sensitivity:
1546 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1548 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1549 checkpoints = delivery.extend_without_payload(delivery_cp)
1550 if delivery.errors:
1551 return delivery
1553 if fake_input_units is not None:
1554 if not self.instrument_sensitivity or \
1555 self.instrument_sensitivity.input_units is None:
1557 delivery.errors.append((
1558 'NoResponseInformation',
1559 'No input units given, so cannot convert to requested '
1560 'units: %s' % fake_input_units,
1561 context))
1563 return delivery
1565 input_units = sanitize_units(
1566 self.instrument_sensitivity.input_units.name)
1568 conresp = None
1569 try:
1570 conresp = conversion[
1571 fake_input_units, input_units]
1573 except KeyError:
1574 delivery.errors.append((
1575 'NoResponseInformation',
1576 'Cannot convert between units: %s, %s'
1577 % (fake_input_units, input_units),
1578 context))
1580 if conresp is not None:
1581 delivery.payload.append(conresp)
1582 for checkpoint in checkpoints:
1583 checkpoint.value *= num.abs(evaluate1(
1584 conresp, checkpoint.frequency))
1586 delivery.payload = [
1587 response.MultiplyResponse(
1588 delivery.payload,
1589 checkpoints=checkpoints)]
1591 return delivery
1593 @classmethod
1594 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1595 normalization_frequency=1.0):
1596 '''
1597 Convert Pyrocko pole-zero response to StationXML response.
1599 :param presponse: Pyrocko pole-zero response
1600 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1601 :param input_unit: Input unit to be reported in the StationXML
1602 response.
1603 :type input_unit: str
1604 :param output_unit: Output unit to be reported in the StationXML
1605 response.
1606 :type output_unit: str
1607 :param normalization_frequency: Frequency where the normalization
1608 factor for the StationXML response should be computed.
1609 :type normalization_frequency: float
1610 '''
1612 norm_factor = 1.0/float(abs(
1613 evaluate1(presponse, normalization_frequency)
1614 / presponse.constant))
1616 pzs = PolesZeros(
1617 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1618 normalization_factor=norm_factor,
1619 normalization_frequency=Frequency(normalization_frequency),
1620 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1621 imaginary=FloatNoUnit(z.imag))
1622 for z in presponse.zeros],
1623 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1624 imaginary=FloatNoUnit(z.imag))
1625 for z in presponse.poles])
1627 pzs.validate()
1629 stage = ResponseStage(
1630 number=1,
1631 poles_zeros_list=[pzs],
1632 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1634 resp = Response(
1635 instrument_sensitivity=Sensitivity(
1636 value=stage.stage_gain.value,
1637 frequency=normalization_frequency,
1638 input_units=Units(input_unit),
1639 output_units=Units(output_unit)),
1641 stage_list=[stage])
1643 return resp
1646class BaseNode(Object):
1647 '''
1648 A base node type for derivation from: Network, Station and Channel types.
1649 '''
1651 code = String.T(xmlstyle='attribute')
1652 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1653 xmlstyle='attribute')
1654 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1655 xmlstyle='attribute')
1656 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1657 alternate_code = String.T(optional=True, xmlstyle='attribute')
1658 historical_code = String.T(optional=True, xmlstyle='attribute')
1659 description = Unicode.T(optional=True, xmltagname='Description')
1660 comment_list = List.T(Comment.T(xmltagname='Comment'))
1662 def spans(self, *args):
1663 if len(args) == 0:
1664 return True
1665 elif len(args) == 1:
1666 return ((self.start_date is None or
1667 self.start_date <= args[0]) and
1668 (self.end_date is None or
1669 args[0] <= self.end_date))
1671 elif len(args) == 2:
1672 return ((self.start_date is None or
1673 args[1] >= self.start_date) and
1674 (self.end_date is None or
1675 self.end_date >= args[0]))
1678def overlaps(a, b):
1679 return (
1680 a.start_date is None or b.end_date is None
1681 or a.start_date < b.end_date
1682 ) and (
1683 b.start_date is None or a.end_date is None
1684 or b.start_date < a.end_date
1685 )
1688def check_overlaps(node_type_name, codes, nodes):
1689 errors = []
1690 for ia, a in enumerate(nodes):
1691 for b in nodes[ia+1:]:
1692 if overlaps(a, b):
1693 errors.append(
1694 '%s epochs overlap for %s:\n'
1695 ' %s - %s\n %s - %s' % (
1696 node_type_name,
1697 '.'.join(codes),
1698 tts(a.start_date), tts(a.end_date),
1699 tts(b.start_date), tts(b.end_date)))
1701 return errors
1704class Channel(BaseNode):
1705 '''
1706 Equivalent to SEED blockette 52 and parent element for the related the
1707 response blockettes.
1708 '''
1710 location_code = String.T(xmlstyle='attribute')
1711 external_reference_list = List.T(
1712 ExternalReference.T(xmltagname='ExternalReference'))
1713 latitude = Latitude.T(xmltagname='Latitude')
1714 longitude = Longitude.T(xmltagname='Longitude')
1715 elevation = Distance.T(xmltagname='Elevation')
1716 depth = Distance.T(xmltagname='Depth')
1717 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1718 dip = Dip.T(optional=True, xmltagname='Dip')
1719 type_list = List.T(Type.T(xmltagname='Type'))
1720 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1721 sample_rate_ratio = SampleRateRatio.T(optional=True,
1722 xmltagname='SampleRateRatio')
1723 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1724 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1725 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1726 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1727 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1728 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1729 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1730 response = Response.T(optional=True, xmltagname='Response')
1732 @property
1733 def position_values(self):
1734 lat = self.latitude.value
1735 lon = self.longitude.value
1736 elevation = value_or_none(self.elevation)
1737 depth = value_or_none(self.depth)
1738 return lat, lon, elevation, depth
1741class Station(BaseNode):
1742 '''
1743 This type represents a Station epoch. It is common to only have a single
1744 station epoch with the station's creation and termination dates as the
1745 epoch start and end dates.
1746 '''
1748 latitude = Latitude.T(xmltagname='Latitude')
1749 longitude = Longitude.T(xmltagname='Longitude')
1750 elevation = Distance.T(xmltagname='Elevation')
1751 site = Site.T(default=Site.D(name=''), xmltagname='Site')
1752 vault = Unicode.T(optional=True, xmltagname='Vault')
1753 geology = Unicode.T(optional=True, xmltagname='Geology')
1754 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1755 operator_list = List.T(Operator.T(xmltagname='Operator'))
1756 creation_date = DummyAwareOptionalTimestamp.T(
1757 optional=True, xmltagname='CreationDate')
1758 termination_date = DummyAwareOptionalTimestamp.T(
1759 optional=True, xmltagname='TerminationDate')
1760 total_number_channels = Counter.T(
1761 optional=True, xmltagname='TotalNumberChannels')
1762 selected_number_channels = Counter.T(
1763 optional=True, xmltagname='SelectedNumberChannels')
1764 external_reference_list = List.T(
1765 ExternalReference.T(xmltagname='ExternalReference'))
1766 channel_list = List.T(Channel.T(xmltagname='Channel'))
1768 @property
1769 def position_values(self):
1770 lat = self.latitude.value
1771 lon = self.longitude.value
1772 elevation = value_or_none(self.elevation)
1773 return lat, lon, elevation
1776class Network(BaseNode):
1777 '''
1778 This type represents the Network layer, all station metadata is contained
1779 within this element. The official name of the network or other descriptive
1780 information can be included in the Description element. The Network can
1781 contain 0 or more Stations.
1782 '''
1784 total_number_stations = Counter.T(optional=True,
1785 xmltagname='TotalNumberStations')
1786 selected_number_stations = Counter.T(optional=True,
1787 xmltagname='SelectedNumberStations')
1788 station_list = List.T(Station.T(xmltagname='Station'))
1790 @property
1791 def station_code_list(self):
1792 return sorted(set(s.code for s in self.station_list))
1794 @property
1795 def sl_code_list(self):
1796 sls = set()
1797 for station in self.station_list:
1798 for channel in station.channel_list:
1799 sls.add((station.code, channel.location_code))
1801 return sorted(sls)
1803 def summary(self, width=80, indent=4):
1804 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1805 lines = ['%s (%i):' % (self.code, len(sls))]
1806 if sls:
1807 ssls = ['.'.join(x for x in c if x) for c in sls]
1808 w = max(len(x) for x in ssls)
1809 n = (width - indent) / (w+1)
1810 while ssls:
1811 lines.append(
1812 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1814 ssls[:n] = []
1816 return '\n'.join(lines)
1819def value_or_none(x):
1820 if x is not None:
1821 return x.value
1822 else:
1823 return None
1826def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1828 pos = lat, lon, elevation, depth = \
1829 channels[0].position_values
1831 if not all(pos == x.position_values for x in channels):
1832 info = '\n'.join(
1833 ' %s: %s' % (x.code, x.position_values) for
1834 x in channels)
1836 mess = 'encountered inconsistencies in channel ' \
1837 'lat/lon/elevation/depth ' \
1838 'for %s.%s.%s: \n%s' % (nsl + (info,))
1840 if inconsistencies == 'raise':
1841 raise InconsistentChannelLocations(mess)
1843 elif inconsistencies == 'warn':
1844 logger.warning(mess)
1845 logger.warning(' -> using mean values')
1847 apos = num.array([x.position_values for x in channels], dtype=float)
1848 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1849 / num.sum(num.isfinite(apos), axis=0)
1851 groups = {}
1852 for channel in channels:
1853 if channel.code not in groups:
1854 groups[channel.code] = []
1856 groups[channel.code].append(channel)
1858 pchannels = []
1859 for code in sorted(groups.keys()):
1860 data = [
1861 (channel.code, value_or_none(channel.azimuth),
1862 value_or_none(channel.dip))
1863 for channel in groups[code]]
1865 azimuth, dip = util.consistency_merge(
1866 data,
1867 message='channel orientation values differ:',
1868 error=inconsistencies)
1870 pchannels.append(
1871 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1873 return pyrocko.model.Station(
1874 *nsl,
1875 lat=mlat,
1876 lon=mlon,
1877 elevation=mele,
1878 depth=mdep,
1879 channels=pchannels)
1882class FDSNStationXML(Object):
1883 '''
1884 Top-level type for Station XML. Required field are Source (network ID of
1885 the institution sending the message) and one or more Network containers or
1886 one or more Station containers.
1887 '''
1889 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1890 source = String.T(xmltagname='Source')
1891 sender = String.T(optional=True, xmltagname='Sender')
1892 module = String.T(optional=True, xmltagname='Module')
1893 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1894 created = Timestamp.T(optional=True, xmltagname='Created')
1895 network_list = List.T(Network.T(xmltagname='Network'))
1897 xmltagname = 'FDSNStationXML'
1898 guessable_xmlns = [guts_xmlns]
1900 def __init__(self, *args, **kwargs):
1901 Object.__init__(self, *args, **kwargs)
1902 if self.created is None:
1903 self.created = time.time()
1905 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1906 time=None, timespan=None,
1907 inconsistencies='warn',
1908 sloppy=False):
1910 assert inconsistencies in ('raise', 'warn')
1912 if nslcs is not None:
1913 nslcs = set(nslcs)
1915 if nsls is not None:
1916 nsls = set(nsls)
1918 tt = ()
1919 if time is not None:
1920 tt = (time,)
1921 elif timespan is not None:
1922 tt = timespan
1924 ns_have = set()
1925 pstations = []
1926 sensor_to_channels = defaultdict(list)
1927 for network in self.network_list:
1928 if not network.spans(*tt):
1929 continue
1931 for station in network.station_list:
1932 if not station.spans(*tt):
1933 continue
1935 if station.channel_list:
1936 loc_to_channels = {}
1937 for channel in station.channel_list:
1938 if not channel.spans(*tt):
1939 continue
1941 loc = channel.location_code.strip()
1942 if loc not in loc_to_channels:
1943 loc_to_channels[loc] = []
1945 loc_to_channels[loc].append(channel)
1947 for loc in sorted(loc_to_channels.keys()):
1948 channels = loc_to_channels[loc]
1949 if nslcs is not None:
1950 channels = [channel for channel in channels
1951 if (network.code, station.code, loc,
1952 channel.code) in nslcs]
1954 if not channels:
1955 continue
1957 nsl = network.code, station.code, loc
1958 if nsls is not None and nsl not in nsls:
1959 continue
1961 for channel in channels:
1962 k = (nsl, channel.code[:-1])
1963 if not sloppy:
1964 k += (channel.start_date, channel.end_date)
1966 sensor_to_channels[k].append(channel)
1968 else:
1969 ns = (network.code, station.code)
1970 if ns in ns_have:
1971 message = 'Duplicate station ' \
1972 '(multiple epochs match): %s.%s ' % ns
1974 if inconsistencies == 'raise':
1975 raise Inconsistencies(message)
1976 else:
1977 logger.warning(message)
1979 ns_have.add(ns)
1981 pstations.append(pyrocko.model.Station(
1982 network.code, station.code, '*',
1983 lat=station.latitude.value,
1984 lon=station.longitude.value,
1985 elevation=value_or_none(station.elevation),
1986 name=station.description or ''))
1988 sensor_have = set()
1989 for k, channels in sensor_to_channels.items():
1990 (nsl, bi) = k[:2]
1991 if (nsl, bi) in sensor_have:
1992 message = 'Duplicate station ' \
1993 '(multiple epochs match): %s.%s.%s' % nsl
1995 if inconsistencies == 'raise':
1996 raise Inconsistencies(message)
1997 else:
1998 logger.warning(message)
2000 sensor_have.add((nsl, bi))
2002 pstations.append(
2003 pyrocko_station_from_channels(
2004 nsl,
2005 channels,
2006 inconsistencies=inconsistencies))
2008 return pstations
2010 @classmethod
2011 def from_pyrocko_stations(
2012 cls, pyrocko_stations, add_flat_responses_from=None):
2014 '''
2015 Generate :py:class:`FDSNStationXML` from list of
2016 :py:class;`pyrocko.model.Station` instances.
2018 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
2019 instances.
2020 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
2021 '''
2022 network_dict = defaultdict(list)
2024 if add_flat_responses_from:
2025 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
2026 extra = dict(
2027 response=Response(
2028 instrument_sensitivity=Sensitivity(
2029 value=1.0,
2030 frequency=1.0,
2031 input_units=Units(name=add_flat_responses_from))))
2032 else:
2033 extra = {}
2035 have_offsets = set()
2036 for s in pyrocko_stations:
2038 if s.north_shift != 0.0 or s.east_shift != 0.0:
2039 have_offsets.add(s.nsl())
2041 network, station, location = s.nsl()
2042 channel_list = []
2043 for c in s.channels:
2044 channel_list.append(
2045 Channel(
2046 location_code=location,
2047 code=c.name,
2048 latitude=Latitude(value=s.effective_lat),
2049 longitude=Longitude(value=s.effective_lon),
2050 elevation=Distance(value=s.elevation),
2051 depth=Distance(value=s.depth),
2052 azimuth=Azimuth(value=c.azimuth),
2053 dip=Dip(value=c.dip),
2054 **extra
2055 )
2056 )
2058 network_dict[network].append(
2059 Station(
2060 code=station,
2061 latitude=Latitude(value=s.effective_lat),
2062 longitude=Longitude(value=s.effective_lon),
2063 elevation=Distance(value=s.elevation),
2064 channel_list=channel_list)
2065 )
2067 if have_offsets:
2068 logger.warning(
2069 'StationXML does not support Cartesian offsets in '
2070 'coordinates. Storing effective lat/lon for stations: %s' %
2071 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2073 timestamp = util.to_time_float(time.time())
2074 network_list = []
2075 for k, station_list in network_dict.items():
2077 network_list.append(
2078 Network(
2079 code=k, station_list=station_list,
2080 total_number_stations=len(station_list)))
2082 sxml = FDSNStationXML(
2083 source='from pyrocko stations list', created=timestamp,
2084 network_list=network_list)
2086 sxml.validate()
2087 return sxml
2089 def iter_network_stations(
2090 self, net=None, sta=None, time=None, timespan=None):
2092 tt = ()
2093 if time is not None:
2094 tt = (time,)
2095 elif timespan is not None:
2096 tt = timespan
2098 for network in self.network_list:
2099 if not network.spans(*tt) or (
2100 net is not None and network.code != net):
2101 continue
2103 for station in network.station_list:
2104 if not station.spans(*tt) or (
2105 sta is not None and station.code != sta):
2106 continue
2108 yield (network, station)
2110 def iter_network_station_channels(
2111 self, net=None, sta=None, loc=None, cha=None,
2112 time=None, timespan=None):
2114 if loc is not None:
2115 loc = loc.strip()
2117 tt = ()
2118 if time is not None:
2119 tt = (time,)
2120 elif timespan is not None:
2121 tt = timespan
2123 for network in self.network_list:
2124 if not network.spans(*tt) or (
2125 net is not None and network.code != net):
2126 continue
2128 for station in network.station_list:
2129 if not station.spans(*tt) or (
2130 sta is not None and station.code != sta):
2131 continue
2133 if station.channel_list:
2134 for channel in station.channel_list:
2135 if (not channel.spans(*tt) or
2136 (cha is not None and channel.code != cha) or
2137 (loc is not None and
2138 channel.location_code.strip() != loc)):
2139 continue
2141 yield (network, station, channel)
2143 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2144 time=None, timespan=None):
2146 groups = {}
2147 for network, station, channel in self.iter_network_station_channels(
2148 net, sta, loc, cha, time=time, timespan=timespan):
2150 net = network.code
2151 sta = station.code
2152 cha = channel.code
2153 loc = channel.location_code.strip()
2154 if len(cha) == 3:
2155 bic = cha[:2] # band and intrument code according to SEED
2156 elif len(cha) == 1:
2157 bic = ''
2158 else:
2159 bic = cha
2161 if channel.response and \
2162 channel.response.instrument_sensitivity and \
2163 channel.response.instrument_sensitivity.input_units:
2165 unit = sanitize_units(
2166 channel.response.instrument_sensitivity.input_units.name)
2167 else:
2168 unit = None
2170 bic = (bic, unit)
2172 k = net, sta, loc
2173 if k not in groups:
2174 groups[k] = {}
2176 if bic not in groups[k]:
2177 groups[k][bic] = []
2179 groups[k][bic].append(channel)
2181 for nsl, bic_to_channels in groups.items():
2182 bad_bics = []
2183 for bic, channels in bic_to_channels.items():
2184 sample_rates = []
2185 for channel in channels:
2186 sample_rates.append(channel.sample_rate.value)
2188 if not same(sample_rates):
2189 scs = ','.join(channel.code for channel in channels)
2190 srs = ', '.join('%e' % x for x in sample_rates)
2191 err = 'ignoring channels with inconsistent sampling ' + \
2192 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2194 logger.warning(err)
2195 bad_bics.append(bic)
2197 for bic in bad_bics:
2198 del bic_to_channels[bic]
2200 return groups
2202 def choose_channels(
2203 self,
2204 target_sample_rate=None,
2205 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2206 priority_units=['M/S', 'M/S**2'],
2207 priority_instrument_code=['H', 'L'],
2208 time=None,
2209 timespan=None):
2211 nslcs = {}
2212 for nsl, bic_to_channels in self.get_channel_groups(
2213 time=time, timespan=timespan).items():
2215 useful_bics = []
2216 for bic, channels in bic_to_channels.items():
2217 rate = channels[0].sample_rate.value
2219 if target_sample_rate is not None and \
2220 rate < target_sample_rate*0.99999:
2221 continue
2223 if len(bic[0]) == 2:
2224 if bic[0][0] not in priority_band_code:
2225 continue
2227 if bic[0][1] not in priority_instrument_code:
2228 continue
2230 unit = bic[1]
2232 prio_unit = len(priority_units)
2233 try:
2234 prio_unit = priority_units.index(unit)
2235 except ValueError:
2236 pass
2238 prio_inst = len(priority_instrument_code)
2239 prio_band = len(priority_band_code)
2240 if len(channels[0].code) == 3:
2241 try:
2242 prio_inst = priority_instrument_code.index(
2243 channels[0].code[1])
2244 except ValueError:
2245 pass
2247 try:
2248 prio_band = priority_band_code.index(
2249 channels[0].code[0])
2250 except ValueError:
2251 pass
2253 if target_sample_rate is None:
2254 rate = -rate
2256 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2257 prio_inst, bic))
2259 useful_bics.sort()
2261 for _, _, rate, _, _, bic in useful_bics:
2262 channels = sorted(
2263 bic_to_channels[bic],
2264 key=lambda channel: channel.code)
2266 if channels:
2267 for channel in channels:
2268 nslcs[nsl + (channel.code,)] = channel
2270 break
2272 return nslcs
2274 def get_pyrocko_response(
2275 self, nslc,
2276 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2278 net, sta, loc, cha = nslc
2279 resps = []
2280 for _, _, channel in self.iter_network_station_channels(
2281 net, sta, loc, cha, time=time, timespan=timespan):
2282 resp = channel.response
2283 if resp:
2284 resp.check_sample_rates(channel)
2285 resp.check_units()
2286 resps.append(resp.get_pyrocko_response(
2287 '.'.join(nslc),
2288 fake_input_units=fake_input_units,
2289 stages=stages).expect_one())
2291 if not resps:
2292 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2293 elif len(resps) > 1:
2294 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2296 return resps[0]
2298 @property
2299 def n_code_list(self):
2300 return sorted(set(x.code for x in self.network_list))
2302 @property
2303 def ns_code_list(self):
2304 nss = set()
2305 for network in self.network_list:
2306 for station in network.station_list:
2307 nss.add((network.code, station.code))
2309 return sorted(nss)
2311 @property
2312 def nsl_code_list(self):
2313 nsls = set()
2314 for network in self.network_list:
2315 for station in network.station_list:
2316 for channel in station.channel_list:
2317 nsls.add(
2318 (network.code, station.code, channel.location_code))
2320 return sorted(nsls)
2322 @property
2323 def nslc_code_list(self):
2324 nslcs = set()
2325 for network in self.network_list:
2326 for station in network.station_list:
2327 for channel in station.channel_list:
2328 nslcs.add(
2329 (network.code, station.code, channel.location_code,
2330 channel.code))
2332 return sorted(nslcs)
2334 def summary(self):
2335 lst = [
2336 'number of n codes: %i' % len(self.n_code_list),
2337 'number of ns codes: %i' % len(self.ns_code_list),
2338 'number of nsl codes: %i' % len(self.nsl_code_list),
2339 'number of nslc codes: %i' % len(self.nslc_code_list)
2340 ]
2341 return '\n'.join(lst)
2343 def summary_stages(self):
2344 data = []
2345 for network, station, channel in self.iter_network_station_channels():
2346 nslc = (network.code, station.code, channel.location_code,
2347 channel.code)
2349 stages = []
2350 in_units = '?'
2351 out_units = '?'
2352 if channel.response:
2353 sens = channel.response.instrument_sensitivity
2354 if sens:
2355 in_units = sanitize_units(sens.input_units.name)
2356 out_units = sanitize_units(sens.output_units.name)
2358 for stage in channel.response.stage_list:
2359 stages.append(stage.summary())
2361 data.append(
2362 (nslc, tts(channel.start_date), tts(channel.end_date),
2363 in_units, out_units, stages))
2365 data.sort()
2367 lst = []
2368 for nslc, stmin, stmax, in_units, out_units, stages in data:
2369 lst.append(' %s: %s - %s, %s -> %s' % (
2370 '.'.join(nslc), stmin, stmax, in_units, out_units))
2371 for stage in stages:
2372 lst.append(' %s' % stage)
2374 return '\n'.join(lst)
2376 def _check_overlaps(self):
2377 by_nslc = {}
2378 for network in self.network_list:
2379 for station in network.station_list:
2380 for channel in station.channel_list:
2381 nslc = (network.code, station.code, channel.location_code,
2382 channel.code)
2383 if nslc not in by_nslc:
2384 by_nslc[nslc] = []
2386 by_nslc[nslc].append(channel)
2388 errors = []
2389 for nslc, channels in by_nslc.items():
2390 errors.extend(check_overlaps('Channel', nslc, channels))
2392 return errors
2394 def check(self):
2395 errors = []
2396 for meth in [self._check_overlaps]:
2397 errors.extend(meth())
2399 if errors:
2400 raise Inconsistencies(
2401 'Inconsistencies found in StationXML:\n '
2402 + '\n '.join(errors))
2405def load_channel_table(stream):
2407 networks = {}
2408 stations = {}
2410 for line in stream:
2411 line = str(line.decode('ascii'))
2412 if line.startswith('#'):
2413 continue
2415 t = line.rstrip().split('|')
2417 if len(t) != 17:
2418 logger.warning('Invalid channel record: %s' % line)
2419 continue
2421 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2422 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2424 try:
2425 scale = float(scale)
2426 except ValueError:
2427 scale = None
2429 try:
2430 scale_freq = float(scale_freq)
2431 except ValueError:
2432 scale_freq = None
2434 try:
2435 depth = float(dep)
2436 except ValueError:
2437 depth = 0.0
2439 try:
2440 azi = float(azi)
2441 dip = float(dip)
2442 except ValueError:
2443 azi = None
2444 dip = None
2446 try:
2447 if net not in networks:
2448 network = Network(code=net)
2449 else:
2450 network = networks[net]
2452 if (net, sta) not in stations:
2453 station = Station(
2454 code=sta, latitude=lat, longitude=lon, elevation=ele)
2456 station.regularize()
2457 else:
2458 station = stations[net, sta]
2460 if scale:
2461 resp = Response(
2462 instrument_sensitivity=Sensitivity(
2463 value=scale,
2464 frequency=scale_freq,
2465 input_units=scale_units))
2466 else:
2467 resp = None
2469 channel = Channel(
2470 code=cha,
2471 location_code=loc.strip(),
2472 latitude=lat,
2473 longitude=lon,
2474 elevation=ele,
2475 depth=depth,
2476 azimuth=azi,
2477 dip=dip,
2478 sensor=Equipment(description=sens),
2479 response=resp,
2480 sample_rate=sample_rate,
2481 start_date=start_date,
2482 end_date=end_date or None)
2484 channel.regularize()
2486 except ValidationError:
2487 raise InvalidRecord(line)
2489 if net not in networks:
2490 networks[net] = network
2492 if (net, sta) not in stations:
2493 stations[net, sta] = station
2494 network.station_list.append(station)
2496 station.channel_list.append(channel)
2498 return FDSNStationXML(
2499 source='created from table input',
2500 created=time.time(),
2501 network_list=sorted(networks.values(), key=lambda x: x.code))
2504def primitive_merge(sxs):
2505 networks = []
2506 for sx in sxs:
2507 networks.extend(sx.network_list)
2509 return FDSNStationXML(
2510 source='merged from different sources',
2511 created=time.time(),
2512 network_list=copy.deepcopy(
2513 sorted(networks, key=lambda x: x.code)))
2516def split_channels(sx):
2517 for nslc in sx.nslc_code_list:
2518 network_list = sx.network_list
2519 network_list_filtered = [
2520 network for network in network_list
2521 if network.code == nslc[0]]
2523 for network in network_list_filtered:
2524 sx.network_list = [network]
2525 station_list = network.station_list
2526 station_list_filtered = [
2527 station for station in station_list
2528 if station.code == nslc[1]]
2530 for station in station_list_filtered:
2531 network.station_list = [station]
2532 channel_list = station.channel_list
2533 station.channel_list = [
2534 channel for channel in channel_list
2535 if (channel.location_code, channel.code) == nslc[2:4]]
2537 yield nslc, copy.deepcopy(sx)
2539 station.channel_list = channel_list
2541 network.station_list = station_list
2543 sx.network_list = network_list
2546if __name__ == '__main__':
2547 from optparse import OptionParser
2549 util.setup_logging('pyrocko.io.stationxml', 'warning')
2551 usage = \
2552 'python -m pyrocko.io.stationxml check|stats|stages ' \
2553 '<filename> [options]'
2555 description = '''Torture StationXML file.'''
2557 parser = OptionParser(
2558 usage=usage,
2559 description=description,
2560 formatter=util.BetterHelpFormatter())
2562 (options, args) = parser.parse_args(sys.argv[1:])
2564 if len(args) != 2:
2565 parser.print_help()
2566 sys.exit(1)
2568 action, path = args
2570 sx = load_xml(filename=path)
2571 if action == 'check':
2572 try:
2573 sx.check()
2574 except Inconsistencies as e:
2575 logger.error(e)
2576 sys.exit(1)
2578 elif action == 'stats':
2579 print(sx.summary())
2581 elif action == 'stages':
2582 print(sx.summary_stages())
2584 else:
2585 parser.print_help()
2586 sys.exit('unknown action: %s' % action)