Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/io/stationxml.py: 70%
1190 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +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)' % name
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 if value == 0.0:
272 return Delivery(log=[(
273 'warning',
274 '%s\n'
275 ' response value reported in file is zero' % prelude,
276 context)])
278 if value < 0.0:
279 return Delivery(log=[(
280 'warning',
281 '%s\n'
282 ' response value reported in file is negative' % prelude,
283 context)])
285 diff_db = 20.0 * num.log10(value_resp/abs(value))
287 if num.abs(diff_db) > limit_db:
288 return Delivery(log=[(
289 'warning',
290 '%s\n'
291 ' reported value: %g\n'
292 ' computed value: %g\n'
293 ' at frequency [Hz]: %g\n'
294 ' factor: %g\n'
295 ' difference [dB]: %g\n'
296 ' limit [dB]: %g' % (
297 prelude,
298 abs(value),
299 value_resp,
300 frequency,
301 value_resp/abs(value),
302 diff_db,
303 limit_db),
304 context)])
306 return Delivery()
309def tts(t):
310 if t is None:
311 return '<none>'
312 else:
313 return util.tts(t, format='%Y-%m-%d %H:%M:%S')
316def le_open_left(a, b):
317 return a is None or (b is not None and a <= b)
320def le_open_right(a, b):
321 return b is None or (a is not None and a <= b)
324def eq_open(a, b):
325 return (a is None and b is None) \
326 or (a is not None and b is not None and a == b)
329def contains(a, b):
330 return le_open_left(a.start_date, b.start_date) \
331 and le_open_right(b.end_date, a.end_date)
334def find_containing(candidates, node):
335 for candidate in candidates:
336 if contains(candidate, node):
337 return candidate
339 return None
342this_year = time.gmtime()[0]
345class DummyAwareOptionalTimestamp(Object):
346 '''
347 Optional timestamp with support for some common placeholder values.
349 Some StationXML files contain arbitrary placeholder values for open end
350 intervals, like "2100-01-01". Depending on the time range supported by the
351 system, these dates are translated into ``None`` to prevent crashes with
352 this type.
353 '''
354 dummy_for = (hpfloat, float)
355 dummy_for_description = 'pyrocko.util.get_time_float'
357 class __T(TBase):
359 def regularize_extra(self, val):
360 time_float = get_time_float()
362 if isinstance(val, datetime.datetime):
363 tt = val.utctimetuple()
364 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
366 elif isinstance(val, datetime.date):
367 tt = val.timetuple()
368 val = time_float(calendar.timegm(tt))
370 elif isinstance(val, str):
371 val = val.strip()
373 tz_offset = 0
375 m = re_tz.search(val)
376 if m:
377 sh = m.group(2)
378 sm = m.group(4)
379 tz_offset = (int(sh)*3600 if sh else 0) \
380 + (int(sm)*60 if sm else 0)
382 val = re_tz.sub('', val)
384 if len(val) > 10 and val[10] == 'T':
385 val = val.replace('T', ' ', 1)
387 try:
388 val = util.str_to_time(val) - tz_offset
390 except util.TimeStrError:
391 year = int(val[:4])
392 if sys.maxsize > 2**32: # if we're on 64bit
393 if year > this_year + 100:
394 return None # StationXML contained a dummy date
396 if year < 1903: # for macOS, 1900-01-01 dummy dates
397 return None
399 else: # 32bit end of time is in 2038
400 if this_year < 2037 and year > 2037 or year < 1903:
401 return None # StationXML contained a dummy date
403 raise
405 elif isinstance(val, (int, float)):
406 val = time_float(val)
408 else:
409 raise ValidationError(
410 '%s: cannot convert "%s" to type %s' % (
411 self.xname(), val, time_float))
413 return val
415 def to_save(self, val):
416 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
417 .rstrip('0').rstrip('.')
419 def to_save_xml(self, val):
420 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
421 .rstrip('0').rstrip('.') + 'Z'
424class Nominal(StringChoice):
425 choices = [
426 'NOMINAL',
427 'CALCULATED']
430class Email(UnicodePattern):
431 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
434class RestrictedStatus(StringChoice):
435 choices = [
436 'open',
437 'closed',
438 'partial']
441class Type(StringChoice):
442 choices = [
443 'TRIGGERED',
444 'CONTINUOUS',
445 'HEALTH',
446 'GEOPHYSICAL',
447 'WEATHER',
448 'FLAG',
449 'SYNTHESIZED',
450 'INPUT',
451 'EXPERIMENTAL',
452 'MAINTENANCE',
453 'BEAM']
455 class __T(StringChoice.T):
456 def validate_extra(self, val):
457 if val not in self.choices:
458 logger.warning(
459 'channel type: "%s" is not a valid choice out of %s' %
460 (val, repr(self.choices)))
463class PzTransferFunction(StringChoice):
464 choices = [
465 'LAPLACE (RADIANS/SECOND)',
466 'LAPLACE (HERTZ)',
467 'DIGITAL (Z-TRANSFORM)']
470class Symmetry(StringChoice):
471 choices = [
472 'NONE',
473 'EVEN',
474 'ODD']
477class CfTransferFunction(StringChoice):
479 class __T(StringChoice.T):
480 def validate(self, val, regularize=False, depth=-1):
481 if regularize:
482 try:
483 val = str(val)
484 except ValueError:
485 raise ValidationError(
486 '%s: cannot convert to string %s' % (self.xname,
487 repr(val)))
489 val = self._dummy_cls.replacements.get(val, val)
491 self.validate_extra(val)
492 return val
494 choices = [
495 'ANALOG (RADIANS/SECOND)',
496 'ANALOG (HERTZ)',
497 'DIGITAL']
499 replacements = {
500 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
501 'ANALOG (HZ)': 'ANALOG (HERTZ)',
502 }
505class Approximation(StringChoice):
506 choices = [
507 'MACLAURIN']
510class PhoneNumber(StringPattern):
511 pattern = '[0-9]+-[0-9]+'
514class Site(Object):
515 '''
516 Description of a site location using name and optional geopolitical
517 boundaries (country, city, etc.).
518 '''
520 name = Unicode.T(default='', xmltagname='Name')
521 description = Unicode.T(optional=True, xmltagname='Description')
522 town = Unicode.T(optional=True, xmltagname='Town')
523 county = Unicode.T(optional=True, xmltagname='County')
524 region = Unicode.T(optional=True, xmltagname='Region')
525 country = Unicode.T(optional=True, xmltagname='Country')
528class ExternalReference(Object):
529 '''
530 This type contains a URI and description for external data that users may
531 want to reference in StationXML.
532 '''
534 uri = String.T(xmltagname='URI')
535 description = Unicode.T(xmltagname='Description')
538class Units(Object):
539 '''
540 A type to document units. Corresponds to SEED blockette 34.
541 '''
543 def __init__(self, name=None, **kwargs):
544 Object.__init__(self, name=name, **kwargs)
546 name = String.T(xmltagname='Name')
547 description = Unicode.T(optional=True, xmltagname='Description')
550class Counter(Int):
551 pass
554class SampleRateRatio(Object):
555 '''
556 Sample rate expressed as number of samples in a number of seconds.
557 '''
559 number_samples = Int.T(xmltagname='NumberSamples')
560 number_seconds = Int.T(xmltagname='NumberSeconds')
563class Gain(Object):
564 '''
565 Complex type for sensitivity and frequency ranges. This complex type can be
566 used to represent both overall sensitivities and individual stage gains.
567 The FrequencyRangeGroup is an optional construct that defines a pass band
568 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
569 is valid within the number of decibels specified in FrequencyDBVariation.
570 '''
572 def __init__(self, value=None, **kwargs):
573 Object.__init__(self, value=value, **kwargs)
575 value = Float.T(optional=True, xmltagname='Value')
576 frequency = Float.T(optional=True, xmltagname='Frequency')
578 def summary(self):
579 return 'gain(%g @ %g)' % (self.value, self.frequency)
582class NumeratorCoefficient(Object):
583 i = Int.T(optional=True, xmlstyle='attribute')
584 value = Float.T(xmlstyle='content')
587class FloatNoUnit(Object):
588 def __init__(self, value=None, **kwargs):
589 Object.__init__(self, value=value, **kwargs)
591 plus_error = Float.T(optional=True, xmlstyle='attribute')
592 minus_error = Float.T(optional=True, xmlstyle='attribute')
593 value = Float.T(xmlstyle='content')
596class FloatWithUnit(FloatNoUnit):
597 unit = String.T(optional=True, xmlstyle='attribute')
600class Equipment(Object):
601 resource_id = String.T(optional=True, xmlstyle='attribute')
602 type = String.T(optional=True, xmltagname='Type')
603 description = Unicode.T(optional=True, xmltagname='Description')
604 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
605 vendor = Unicode.T(optional=True, xmltagname='Vendor')
606 model = Unicode.T(optional=True, xmltagname='Model')
607 serial_number = String.T(optional=True, xmltagname='SerialNumber')
608 installation_date = DummyAwareOptionalTimestamp.T(
609 optional=True,
610 xmltagname='InstallationDate')
611 removal_date = DummyAwareOptionalTimestamp.T(
612 optional=True,
613 xmltagname='RemovalDate')
614 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
617class PhoneNumber(Object):
618 description = Unicode.T(optional=True, xmlstyle='attribute')
619 country_code = Int.T(optional=True, xmltagname='CountryCode')
620 area_code = Int.T(xmltagname='AreaCode')
621 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
624class BaseFilter(Object):
625 '''
626 The BaseFilter is derived by all filters.
627 '''
629 resource_id = String.T(optional=True, xmlstyle='attribute')
630 name = String.T(optional=True, xmlstyle='attribute')
631 description = Unicode.T(optional=True, xmltagname='Description')
632 input_units = Units.T(optional=True, xmltagname='InputUnits')
633 output_units = Units.T(optional=True, xmltagname='OutputUnits')
636class Sensitivity(Gain):
637 '''
638 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
639 construct that defines a pass band in Hertz (FrequencyStart and
640 FrequencyEnd) in which the SensitivityValue is valid within the number of
641 decibels specified in FrequencyDBVariation.
642 '''
644 input_units = Units.T(optional=True, xmltagname='InputUnits')
645 output_units = Units.T(optional=True, xmltagname='OutputUnits')
646 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
647 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
648 frequency_db_variation = Float.T(optional=True,
649 xmltagname='FrequencyDBVariation')
651 def get_pyrocko_response(self):
652 return Delivery(
653 [response.PoleZeroResponse(constant=self.value)])
656class Coefficient(FloatNoUnit):
657 number = Counter.T(optional=True, xmlstyle='attribute')
660class PoleZero(Object):
661 '''
662 Complex numbers used as poles or zeros in channel response.
663 '''
665 number = Int.T(optional=True, xmlstyle='attribute')
666 real = FloatNoUnit.T(xmltagname='Real')
667 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
669 def value(self):
670 return self.real.value + 1J * self.imaginary.value
673class ClockDrift(FloatWithUnit):
674 unit = String.T(default='SECONDS/SAMPLE', optional=True,
675 xmlstyle='attribute') # fixed
678class Second(FloatWithUnit):
679 '''
680 A time value in seconds.
681 '''
683 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
684 # fixed unit
687class Voltage(FloatWithUnit):
688 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
689 # fixed unit
692class Angle(FloatWithUnit):
693 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
694 # fixed unit
697class Azimuth(FloatWithUnit):
698 '''
699 Instrument azimuth, degrees clockwise from North.
700 '''
702 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
703 # fixed unit
706class Dip(FloatWithUnit):
707 '''
708 Instrument dip in degrees down from horizontal. Together azimuth and dip
709 describe the direction of the sensitive axis of the instrument.
710 '''
712 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
713 # fixed unit
716class Distance(FloatWithUnit):
717 '''
718 Extension of FloatWithUnit for distances, elevations, and depths.
719 '''
721 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
722 # NOT fixed unit!
725class Frequency(FloatWithUnit):
726 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
727 # fixed unit
730class SampleRate(FloatWithUnit):
731 '''
732 Sample rate in samples per second.
733 '''
735 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
736 # fixed unit
739class Person(Object):
740 '''
741 Representation of a person's contact information. A person can belong to
742 multiple agencies and have multiple email addresses and phone numbers.
743 '''
745 name_list = List.T(Unicode.T(xmltagname='Name'))
746 agency_list = List.T(Unicode.T(xmltagname='Agency'))
747 email_list = List.T(Email.T(xmltagname='Email'))
748 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
751class FIR(BaseFilter):
752 '''
753 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
754 also commonly documented using the Coefficients element.
755 '''
757 symmetry = Symmetry.T(xmltagname='Symmetry')
758 numerator_coefficient_list = List.T(
759 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
761 def summary(self):
762 return 'fir(%i%s)' % (
763 self.get_ncoefficiencs(),
764 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
766 def get_effective_coefficients(self):
767 b = num.array(
768 [v.value for v in self.numerator_coefficient_list],
769 dtype=float)
771 if self.symmetry == 'ODD':
772 b = num.concatenate((b, b[-2::-1]))
773 elif self.symmetry == 'EVEN':
774 b = num.concatenate((b, b[::-1]))
776 return b
778 def get_effective_symmetry(self):
779 if self.symmetry == 'NONE':
780 b = self.get_effective_coefficients()
781 if num.all(b - b[::-1] == 0):
782 return ['EVEN', 'ODD'][b.size % 2]
784 return self.symmetry
786 def get_ncoefficiencs(self):
787 nf = len(self.numerator_coefficient_list)
788 if self.symmetry == 'ODD':
789 nc = nf * 2 + 1
790 elif self.symmetry == 'EVEN':
791 nc = nf * 2
792 else:
793 nc = nf
795 return nc
797 def estimate_delay(self, deltat):
798 nc = self.get_ncoefficiencs()
799 if nc > 0:
800 return deltat * (nc - 1) / 2.0
801 else:
802 return 0.0
804 def get_pyrocko_response(
805 self, context, deltat, delay_responses, normalization_frequency):
807 context += self.summary()
809 if not self.numerator_coefficient_list:
810 return Delivery([])
812 b = self.get_effective_coefficients()
814 log = []
815 drop_phase = self.get_effective_symmetry() != 'NONE'
817 if not deltat:
818 log.append((
819 'error',
820 'Digital response requires knowledge about sampling '
821 'interval. Response will be unusable.',
822 context))
824 resp = response.DigitalFilterResponse(
825 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
827 if normalization_frequency is not None and deltat is not None:
828 normalization_frequency = 0.0
829 normalization = num.abs(evaluate1(resp, normalization_frequency))
831 if num.abs(normalization - 1.0) > 1e-2:
832 log.append((
833 'warning',
834 'FIR filter coefficients are not normalized. Normalizing '
835 'them. Factor: %g' % normalization, context))
837 resp = response.DigitalFilterResponse(
838 (b/normalization).tolist(), [1.0], deltat,
839 drop_phase=drop_phase)
841 resps = [resp]
843 if not drop_phase:
844 resps.extend(delay_responses)
846 return Delivery(resps, log=log)
849class Coefficients(BaseFilter):
850 '''
851 Response: coefficients for FIR filter. Laplace transforms or IIR filters
852 can be expressed using type as well but the PolesAndZeros should be used
853 instead. Corresponds to SEED blockette 54.
854 '''
856 cf_transfer_function_type = CfTransferFunction.T(
857 xmltagname='CfTransferFunctionType')
858 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
859 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
861 def summary(self):
862 return 'coef_%s(%i,%i%s)' % (
863 'ABC?'[
864 CfTransferFunction.choices.index(
865 self.cf_transfer_function_type)],
866 len(self.numerator_list),
867 len(self.denominator_list),
868 ',sym' if self.is_symmetric_fir else '')
870 def estimate_delay(self, deltat):
871 nc = len(self.numerator_list)
872 if nc > 0:
873 return deltat * (len(self.numerator_list) - 1) / 2.0
874 else:
875 return 0.0
877 def is_symmetric_fir(self):
878 if len(self.denominator_list) != 0:
879 return False
880 b = [v.value for v in self.numerator_list]
881 return b == b[::-1]
883 def get_pyrocko_response(
884 self, context, deltat, delay_responses, normalization_frequency):
886 context += self.summary()
888 factor = 1.0
889 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
890 factor = 2.0*math.pi
892 if not self.numerator_list and not self.denominator_list:
893 return Delivery(payload=[])
895 b = num.array(
896 [v.value*factor for v in self.numerator_list], dtype=float)
898 a = num.array(
899 [1.0] + [v.value*factor for v in self.denominator_list],
900 dtype=float)
902 log = []
903 resps = []
904 if self.cf_transfer_function_type in [
905 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
906 resps.append(response.AnalogFilterResponse(b, a))
908 elif self.cf_transfer_function_type == 'DIGITAL':
909 if not deltat:
910 log.append((
911 'error',
912 'Digital response requires knowledge about sampling '
913 'interval. Response will be unusable.',
914 context))
916 drop_phase = self.is_symmetric_fir()
917 resp = response.DigitalFilterResponse(
918 b, a, deltat or 0.0, drop_phase=drop_phase)
920 if normalization_frequency is not None and deltat is not None:
921 normalization = num.abs(evaluate1(
922 resp, normalization_frequency))
924 if num.abs(normalization - 1.0) > 1e-2:
925 log.append((
926 'warning',
927 'FIR filter coefficients are not normalized. '
928 'Normalizing them. Factor: %g' % normalization,
929 context))
931 resp = response.DigitalFilterResponse(
932 (b/normalization).tolist(), [1.0], deltat,
933 drop_phase=drop_phase)
935 resps.append(resp)
937 if not drop_phase:
938 resps.extend(delay_responses)
940 else:
941 return Delivery(error=(
942 'ValueError',
943 'Unknown transfer function type: %s' % (
944 self.cf_transfer_function_type)))
946 return Delivery(payload=resps, log=log)
949class Latitude(FloatWithUnit):
950 '''
951 Type for latitude coordinate.
952 '''
954 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
955 # fixed unit
956 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
959class Longitude(FloatWithUnit):
960 '''
961 Type for longitude coordinate.
962 '''
964 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
965 # fixed unit
966 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
969class PolesZeros(BaseFilter):
970 '''
971 Response: complex poles and zeros. Corresponds to SEED blockette 53.
972 '''
974 pz_transfer_function_type = PzTransferFunction.T(
975 xmltagname='PzTransferFunctionType')
976 normalization_factor = Float.T(
977 default=1.0,
978 xmltagname='NormalizationFactor')
979 normalization_frequency = Frequency.T(
980 optional=True, # but required by standard
981 xmltagname='NormalizationFrequency')
982 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
983 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
985 def summary(self):
986 return 'pz_%s(%i,%i)' % (
987 'ABC?'[
988 PzTransferFunction.choices.index(
989 self.pz_transfer_function_type)],
990 len(self.pole_list),
991 len(self.zero_list))
993 def get_pyrocko_response(self, context='', deltat=None):
995 context += self.summary()
997 factor = 1.0
998 cfactor = 1.0
999 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
1000 factor = 2. * math.pi
1001 cfactor = (2. * math.pi)**(
1002 len(self.pole_list) - len(self.zero_list))
1004 log = []
1005 if self.normalization_factor is None \
1006 or self.normalization_factor == 0.0:
1008 log.append((
1009 'warning',
1010 'No pole-zero normalization factor given. '
1011 'Assuming a value of 1.0',
1012 context))
1014 nfactor = 1.0
1015 else:
1016 nfactor = self.normalization_factor
1018 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
1019 if not is_digital:
1020 resp = response.PoleZeroResponse(
1021 constant=nfactor*cfactor,
1022 zeros=[z.value()*factor for z in self.zero_list],
1023 poles=[p.value()*factor for p in self.pole_list])
1024 else:
1025 if not deltat:
1026 log.append((
1027 'error',
1028 'Digital response requires knowledge about sampling '
1029 'interval. Response will be unusable.',
1030 context))
1032 resp = response.DigitalPoleZeroResponse(
1033 constant=nfactor*cfactor,
1034 zeros=[z.value()*factor for z in self.zero_list],
1035 poles=[p.value()*factor for p in self.pole_list],
1036 deltat=deltat or 0.0)
1038 if not self.normalization_frequency:
1039 log.append((
1040 'warning',
1041 'Cannot check pole-zero normalization factor: '
1042 'No normalization frequency given.',
1043 context))
1045 else:
1046 if is_digital and not deltat:
1047 log.append((
1048 'warning',
1049 'Cannot check computed vs reported normalization '
1050 'factor without knowing the sampling interval.',
1051 context))
1052 else:
1053 computed_normalization_factor = nfactor / abs(evaluate1(
1054 resp, self.normalization_frequency.value))
1056 db = 20.0 * num.log10(
1057 computed_normalization_factor / nfactor)
1059 if abs(db) > 0.17:
1060 log.append((
1061 'warning',
1062 'Computed and reported normalization factors differ '
1063 'by %g dB: computed: %g, reported: %g' % (
1064 db,
1065 computed_normalization_factor,
1066 nfactor),
1067 context))
1069 return Delivery([resp], log)
1072class ResponseListElement(Object):
1073 frequency = Frequency.T(xmltagname='Frequency')
1074 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1075 phase = Angle.T(xmltagname='Phase')
1078class Polynomial(BaseFilter):
1079 '''
1080 Response: expressed as a polynomial (allows non-linear sensors to be
1081 described). Corresponds to SEED blockette 62. Can be used to describe a
1082 stage of acquisition or a complete system.
1083 '''
1085 approximation_type = Approximation.T(default='MACLAURIN',
1086 xmltagname='ApproximationType')
1087 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1088 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1089 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1090 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1091 maximum_error = Float.T(xmltagname='MaximumError')
1092 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1094 def summary(self):
1095 return 'poly(%i)' % len(self.coefficient_list)
1098class Decimation(Object):
1099 '''
1100 Corresponds to SEED blockette 57.
1101 '''
1103 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1104 factor = Int.T(xmltagname='Factor')
1105 offset = Int.T(xmltagname='Offset')
1106 delay = FloatWithUnit.T(xmltagname='Delay')
1107 correction = FloatWithUnit.T(xmltagname='Correction')
1109 def summary(self):
1110 return 'deci(%i, %g -> %g, %g)' % (
1111 self.factor,
1112 self.input_sample_rate.value,
1113 self.input_sample_rate.value / self.factor,
1114 self.delay.value)
1116 def get_pyrocko_response(self):
1117 if self.delay and self.delay.value != 0.0:
1118 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1119 else:
1120 return Delivery([])
1123class Operator(Object):
1124 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1125 contact_list = List.T(Person.T(xmltagname='Contact'))
1126 web_site = String.T(optional=True, xmltagname='WebSite')
1129class Comment(Object):
1130 '''
1131 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1132 and 59.
1133 '''
1135 id = Counter.T(optional=True, xmlstyle='attribute')
1136 value = Unicode.T(xmltagname='Value')
1137 begin_effective_time = DummyAwareOptionalTimestamp.T(
1138 optional=True,
1139 xmltagname='BeginEffectiveTime')
1140 end_effective_time = DummyAwareOptionalTimestamp.T(
1141 optional=True,
1142 xmltagname='EndEffectiveTime')
1143 author_list = List.T(Person.T(xmltagname='Author'))
1146class ResponseList(BaseFilter):
1147 '''
1148 Response: list of frequency, amplitude and phase values. Corresponds to
1149 SEED blockette 55.
1150 '''
1152 response_list_element_list = List.T(
1153 ResponseListElement.T(xmltagname='ResponseListElement'))
1155 def summary(self):
1156 return 'list(%i)' % len(self.response_list_element_list)
1159class Log(Object):
1160 '''
1161 Container for log entries.
1162 '''
1164 entry_list = List.T(Comment.T(xmltagname='Entry'))
1167class ResponseStage(Object):
1168 '''
1169 This complex type represents channel response and covers SEED blockettes 53
1170 to 56.
1171 '''
1173 number = Counter.T(xmlstyle='attribute')
1174 resource_id = String.T(optional=True, xmlstyle='attribute')
1175 poles_zeros_list = List.T(
1176 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1177 coefficients_list = List.T(
1178 Coefficients.T(optional=True, xmltagname='Coefficients'))
1179 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1180 fir = FIR.T(optional=True, xmltagname='FIR')
1181 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1182 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1183 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1185 def summary(self):
1186 elements = []
1188 for stuff in [
1189 self.poles_zeros_list, self.coefficients_list,
1190 self.response_list, self.fir, self.polynomial,
1191 self.decimation, self.stage_gain]:
1193 if stuff:
1194 if isinstance(stuff, Object):
1195 elements.append(stuff.summary())
1196 else:
1197 elements.extend(obj.summary() for obj in stuff)
1199 return '%i: %s %s -> %s' % (
1200 self.number,
1201 ', '.join(elements),
1202 sanitize_units(self.input_units.name)
1203 if self.input_units else '?',
1204 sanitize_units(self.output_units.name)
1205 if self.output_units else '?')
1207 def get_squirrel_response_stage(self, context):
1208 from pyrocko.squirrel.model import ResponseStage
1210 delivery = Delivery()
1211 delivery_pr = self.get_pyrocko_response(context)
1212 log = delivery_pr.log
1213 delivery_pr.log = []
1214 elements = delivery.extend_without_payload(delivery_pr)
1216 delivery.payload = [ResponseStage(
1217 input_quantity=to_quantity(self.input_units, context, delivery),
1218 output_quantity=to_quantity(self.output_units, context, delivery),
1219 input_sample_rate=self.input_sample_rate,
1220 output_sample_rate=self.output_sample_rate,
1221 elements=elements,
1222 log=log)]
1224 return delivery
1226 def get_pyrocko_response(self, context, gain_only=False):
1228 context = context + ', stage %i' % self.number
1230 responses = []
1231 log = []
1232 if self.stage_gain:
1233 normalization_frequency = self.stage_gain.frequency or 0.0
1234 else:
1235 normalization_frequency = 0.0
1237 if not gain_only:
1238 deltat = None
1239 delay_responses = []
1240 if self.decimation:
1241 rate = self.decimation.input_sample_rate.value
1242 if rate > 0.0:
1243 deltat = 1.0 / rate
1244 delivery = self.decimation.get_pyrocko_response()
1245 if delivery.errors:
1246 return delivery
1248 delay_responses = delivery.payload
1249 log.extend(delivery.log)
1251 for pzs in self.poles_zeros_list:
1252 delivery = pzs.get_pyrocko_response(context, deltat)
1253 if delivery.errors:
1254 return delivery
1256 pz_resps = delivery.payload
1257 log.extend(delivery.log)
1258 responses.extend(pz_resps)
1260 # emulate incorrect? evalresp behaviour
1261 if pzs.normalization_frequency is not None and \
1262 pzs.normalization_frequency.value \
1263 != normalization_frequency \
1264 and normalization_frequency != 0.0:
1266 try:
1267 trial = response.MultiplyResponse(pz_resps)
1268 anorm = num.abs(evaluate1(
1269 trial, pzs.normalization_frequency.value))
1270 asens = num.abs(
1271 evaluate1(trial, normalization_frequency))
1273 factor = anorm/asens
1275 if abs(factor - 1.0) > 0.01:
1276 log.append((
1277 'warning',
1278 'PZ normalization frequency (%g) is different '
1279 'from stage gain frequency (%s) -> Emulating '
1280 'possibly incorrect evalresp behaviour. '
1281 'Correction factor: %g' % (
1282 pzs.normalization_frequency.value,
1283 normalization_frequency,
1284 factor),
1285 context))
1287 responses.append(
1288 response.PoleZeroResponse(constant=factor))
1289 except response.InvalidResponseError as e:
1290 log.append((
1291 'warning',
1292 'Could not check response: %s' % str(e),
1293 context))
1295 if len(self.poles_zeros_list) > 1:
1296 log.append((
1297 'warning',
1298 'Multiple poles and zeros records in single response '
1299 'stage.',
1300 context))
1302 for cfs in self.coefficients_list + (
1303 [self.fir] if self.fir else []):
1305 delivery = cfs.get_pyrocko_response(
1306 context, deltat, delay_responses,
1307 normalization_frequency)
1309 if delivery.errors:
1310 return delivery
1312 responses.extend(delivery.payload)
1313 log.extend(delivery.log)
1315 if len(self.coefficients_list) > 1:
1316 log.append((
1317 'warning',
1318 'Multiple filter coefficients lists in single response '
1319 'stage.',
1320 context))
1322 if self.response_list:
1323 log.append((
1324 'warning',
1325 'Unhandled response element of type: ResponseList',
1326 context))
1328 if self.polynomial:
1329 log.append((
1330 'warning',
1331 'Unhandled response element of type: Polynomial',
1332 context))
1334 if self.stage_gain:
1335 responses.append(
1336 response.PoleZeroResponse(constant=self.stage_gain.value))
1338 return Delivery(responses, log)
1340 @property
1341 def input_units(self):
1342 for e in (self.poles_zeros_list + self.coefficients_list +
1343 [self.response_list, self.fir, self.polynomial]):
1344 if e is not None:
1345 return e.input_units
1347 return None
1349 @property
1350 def output_units(self):
1351 for e in (self.poles_zeros_list + self.coefficients_list +
1352 [self.response_list, self.fir, self.polynomial]):
1353 if e is not None:
1354 return e.output_units
1356 return None
1358 @property
1359 def input_sample_rate(self):
1360 if self.decimation:
1361 return self.decimation.input_sample_rate.value
1363 return None
1365 @property
1366 def output_sample_rate(self):
1367 if self.decimation:
1368 return self.decimation.input_sample_rate.value \
1369 / self.decimation.factor
1371 return None
1374class Response(Object):
1375 resource_id = String.T(optional=True, xmlstyle='attribute')
1376 instrument_sensitivity = Sensitivity.T(optional=True,
1377 xmltagname='InstrumentSensitivity')
1378 instrument_polynomial = Polynomial.T(optional=True,
1379 xmltagname='InstrumentPolynomial')
1380 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1382 @property
1383 def output_sample_rate(self):
1384 if self.stage_list:
1385 return self.stage_list[-1].output_sample_rate
1386 else:
1387 return None
1389 def check_sample_rates(self, channel):
1391 if self.stage_list:
1392 sample_rate = None
1394 for stage in self.stage_list:
1395 if stage.decimation:
1396 input_sample_rate = \
1397 stage.decimation.input_sample_rate.value
1399 if sample_rate is not None and not same_sample_rate(
1400 sample_rate, input_sample_rate):
1402 logger.warning(
1403 'Response stage %i has unexpected input sample '
1404 'rate: %g Hz (expected: %g Hz)' % (
1405 stage.number,
1406 input_sample_rate,
1407 sample_rate))
1409 sample_rate = input_sample_rate / stage.decimation.factor
1411 if sample_rate is not None and channel.sample_rate \
1412 and not same_sample_rate(
1413 sample_rate, channel.sample_rate.value):
1415 logger.warning(
1416 'Channel sample rate (%g Hz) does not match sample rate '
1417 'deducted from response stages information (%g Hz).' % (
1418 channel.sample_rate.value,
1419 sample_rate))
1421 def check_units(self):
1423 if self.instrument_sensitivity \
1424 and self.instrument_sensitivity.input_units:
1426 units = sanitize_units(
1427 self.instrument_sensitivity.input_units.name)
1429 if self.stage_list:
1430 for stage in self.stage_list:
1431 if units and stage.input_units \
1432 and sanitize_units(stage.input_units.name) != units:
1434 logger.warning(
1435 'Input units of stage %i (%s) do not match %s (%s).'
1436 % (
1437 stage.number,
1438 units,
1439 'output units of stage %i'
1440 if stage.number == 0
1441 else 'sensitivity input units',
1442 units))
1444 if stage.output_units:
1445 units = sanitize_units(stage.output_units.name)
1446 else:
1447 units = None
1449 sout_units = self.instrument_sensitivity.output_units
1450 if self.instrument_sensitivity and sout_units:
1451 if units is not None and units != sanitize_units(
1452 sout_units.name):
1453 logger.warning(
1454 'Output units of stage %i (%s) do not match %s (%s).'
1455 % (
1456 stage.number,
1457 units,
1458 'sensitivity output units',
1459 sanitize_units(sout_units.name)))
1461 def _sensitivity_checkpoints(self, responses, context):
1462 delivery = Delivery()
1464 if self.instrument_sensitivity:
1465 sval = self.instrument_sensitivity.value
1466 sfreq = self.instrument_sensitivity.frequency
1467 if sval is None:
1468 delivery.log.append((
1469 'warning',
1470 'No sensitivity value given.',
1471 context))
1473 elif sval is None:
1474 delivery.log.append((
1475 'warning',
1476 'Reported sensitivity value is zero.',
1477 context))
1479 elif sfreq is None:
1480 delivery.log.append((
1481 'warning',
1482 'Sensitivity frequency not given.',
1483 context))
1485 else:
1486 trial = response.MultiplyResponse(responses)
1488 delivery.extend(
1489 check_resp(
1490 trial, sval, sfreq, 0.1,
1491 'Instrument sensitivity value inconsistent with '
1492 'sensitivity computed from complete response.',
1493 context))
1495 delivery.payload.append(response.FrequencyResponseCheckpoint(
1496 frequency=sfreq, value=sval))
1498 return delivery
1500 def get_squirrel_response(self, context, **kwargs):
1501 from pyrocko.squirrel.model import Response
1503 if self.stage_list:
1504 delivery = Delivery()
1505 for istage, stage in enumerate(self.stage_list):
1506 delivery.extend(stage.get_squirrel_response_stage(context))
1508 checkpoints = []
1509 if not delivery.errors:
1510 all_responses = []
1511 for sq_stage in delivery.payload:
1512 all_responses.extend(sq_stage.elements)
1514 checkpoints.extend(delivery.extend_without_payload(
1515 self._sensitivity_checkpoints(all_responses, context)))
1517 sq_stages = delivery.payload
1518 if sq_stages:
1519 if sq_stages[0].input_quantity is None \
1520 and self.instrument_sensitivity is not None:
1522 sq_stages[0].input_quantity = to_quantity(
1523 self.instrument_sensitivity.input_units,
1524 context, delivery)
1525 sq_stages[-1].output_quantity = to_quantity(
1526 self.instrument_sensitivity.output_units,
1527 context, delivery)
1529 sq_stages = delivery.expect()
1531 return Response(
1532 stages=sq_stages,
1533 log=delivery.log,
1534 checkpoints=checkpoints,
1535 **kwargs)
1537 elif self.instrument_sensitivity:
1538 raise NoResponseInformation(
1539 "Only instrument sensitivity given (won't use it). (%s)."
1540 % context)
1541 else:
1542 raise NoResponseInformation(
1543 'Empty instrument response (%s).'
1544 % context)
1546 def get_pyrocko_response(
1547 self, context, fake_input_units=None, stages=(0, 1)):
1549 if fake_input_units is not None:
1550 fake_input_units = sanitize_units(fake_input_units)
1552 delivery = Delivery()
1553 if self.stage_list:
1554 for istage, stage in enumerate(self.stage_list):
1555 delivery.extend(stage.get_pyrocko_response(
1556 context, gain_only=not (
1557 stages is None or stages[0] <= istage < stages[1])))
1559 elif self.instrument_sensitivity:
1560 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1562 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1563 checkpoints = delivery.extend_without_payload(delivery_cp)
1564 if delivery.errors:
1565 return delivery
1567 if fake_input_units is not None:
1568 if not self.instrument_sensitivity or \
1569 self.instrument_sensitivity.input_units is None:
1571 delivery.errors.append((
1572 'NoResponseInformation',
1573 'No input units given, so cannot convert to requested '
1574 'units: %s' % fake_input_units,
1575 context))
1577 return delivery
1579 input_units = sanitize_units(
1580 self.instrument_sensitivity.input_units.name)
1582 conresp = None
1583 try:
1584 conresp = conversion[
1585 fake_input_units, input_units]
1587 except KeyError:
1588 delivery.errors.append((
1589 'NoResponseInformation',
1590 'Cannot convert between units: %s, %s'
1591 % (fake_input_units, input_units),
1592 context))
1594 if conresp is not None:
1595 delivery.payload.append(conresp)
1596 for checkpoint in checkpoints:
1597 checkpoint.value *= num.abs(evaluate1(
1598 conresp, checkpoint.frequency))
1600 delivery.payload = [
1601 response.MultiplyResponse(
1602 delivery.payload,
1603 checkpoints=checkpoints)]
1605 return delivery
1607 @classmethod
1608 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1609 normalization_frequency=1.0):
1610 '''
1611 Convert Pyrocko pole-zero response to StationXML response.
1613 :param presponse: Pyrocko pole-zero response
1614 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1615 :param input_unit: Input unit to be reported in the StationXML
1616 response.
1617 :type input_unit: str
1618 :param output_unit: Output unit to be reported in the StationXML
1619 response.
1620 :type output_unit: str
1621 :param normalization_frequency: Frequency where the normalization
1622 factor for the StationXML response should be computed.
1623 :type normalization_frequency: float
1624 '''
1626 norm_factor = 1.0/float(abs(
1627 evaluate1(presponse, normalization_frequency)
1628 / presponse.constant))
1630 pzs = PolesZeros(
1631 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1632 normalization_factor=norm_factor,
1633 normalization_frequency=Frequency(normalization_frequency),
1634 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1635 imaginary=FloatNoUnit(z.imag))
1636 for z in presponse.zeros],
1637 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1638 imaginary=FloatNoUnit(z.imag))
1639 for z in presponse.poles])
1641 pzs.validate()
1643 stage = ResponseStage(
1644 number=1,
1645 poles_zeros_list=[pzs],
1646 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1648 resp = Response(
1649 instrument_sensitivity=Sensitivity(
1650 value=stage.stage_gain.value,
1651 frequency=normalization_frequency,
1652 input_units=Units(input_unit),
1653 output_units=Units(output_unit)),
1655 stage_list=[stage])
1657 return resp
1660class BaseNode(Object):
1661 '''
1662 A base node type for derivation from: Network, Station and Channel types.
1663 '''
1665 code = String.T(xmlstyle='attribute')
1666 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1667 xmlstyle='attribute')
1668 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1669 xmlstyle='attribute')
1670 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1671 alternate_code = String.T(optional=True, xmlstyle='attribute')
1672 historical_code = String.T(optional=True, xmlstyle='attribute')
1673 description = Unicode.T(optional=True, xmltagname='Description')
1674 comment_list = List.T(Comment.T(xmltagname='Comment'))
1676 def spans(self, *args):
1677 if len(args) == 0:
1678 return True
1679 elif len(args) == 1:
1680 return ((self.start_date is None or
1681 self.start_date <= args[0]) and
1682 (self.end_date is None or
1683 args[0] <= self.end_date))
1685 elif len(args) == 2:
1686 return ((self.start_date is None or
1687 args[1] >= self.start_date) and
1688 (self.end_date is None or
1689 self.end_date >= args[0]))
1692def overlaps(a, b):
1693 return (
1694 a.start_date is None or b.end_date is None
1695 or a.start_date < b.end_date
1696 ) and (
1697 b.start_date is None or a.end_date is None
1698 or b.start_date < a.end_date
1699 )
1702def check_overlaps(node_type_name, codes, nodes):
1703 errors = []
1704 for ia, a in enumerate(nodes):
1705 for b in nodes[ia+1:]:
1706 if overlaps(a, b):
1707 errors.append(
1708 '%s epochs overlap for %s:\n'
1709 ' %s - %s\n %s - %s' % (
1710 node_type_name,
1711 '.'.join(codes),
1712 tts(a.start_date), tts(a.end_date),
1713 tts(b.start_date), tts(b.end_date)))
1715 return errors
1718class Channel(BaseNode):
1719 '''
1720 Equivalent to SEED blockette 52 and parent element for the related the
1721 response blockettes.
1722 '''
1724 location_code = String.T(xmlstyle='attribute')
1725 external_reference_list = List.T(
1726 ExternalReference.T(xmltagname='ExternalReference'))
1727 latitude = Latitude.T(xmltagname='Latitude')
1728 longitude = Longitude.T(xmltagname='Longitude')
1729 elevation = Distance.T(xmltagname='Elevation')
1730 depth = Distance.T(xmltagname='Depth')
1731 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1732 dip = Dip.T(optional=True, xmltagname='Dip')
1733 type_list = List.T(Type.T(xmltagname='Type'))
1734 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1735 sample_rate_ratio = SampleRateRatio.T(optional=True,
1736 xmltagname='SampleRateRatio')
1737 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1738 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1739 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1740 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1741 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1742 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1743 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1744 response = Response.T(optional=True, xmltagname='Response')
1746 @property
1747 def position_values(self):
1748 lat = self.latitude.value
1749 lon = self.longitude.value
1750 elevation = value_or_none(self.elevation)
1751 depth = value_or_none(self.depth)
1752 return lat, lon, elevation, depth
1755class Station(BaseNode):
1756 '''
1757 This type represents a Station epoch. It is common to only have a single
1758 station epoch with the station's creation and termination dates as the
1759 epoch start and end dates.
1760 '''
1762 latitude = Latitude.T(xmltagname='Latitude')
1763 longitude = Longitude.T(xmltagname='Longitude')
1764 elevation = Distance.T(xmltagname='Elevation')
1765 site = Site.T(default=Site.D(name=''), xmltagname='Site')
1766 vault = Unicode.T(optional=True, xmltagname='Vault')
1767 geology = Unicode.T(optional=True, xmltagname='Geology')
1768 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1769 operator_list = List.T(Operator.T(xmltagname='Operator'))
1770 creation_date = DummyAwareOptionalTimestamp.T(
1771 optional=True, xmltagname='CreationDate')
1772 termination_date = DummyAwareOptionalTimestamp.T(
1773 optional=True, xmltagname='TerminationDate')
1774 total_number_channels = Counter.T(
1775 optional=True, xmltagname='TotalNumberChannels')
1776 selected_number_channels = Counter.T(
1777 optional=True, xmltagname='SelectedNumberChannels')
1778 external_reference_list = List.T(
1779 ExternalReference.T(xmltagname='ExternalReference'))
1780 channel_list = List.T(Channel.T(xmltagname='Channel'))
1782 @property
1783 def position_values(self):
1784 lat = self.latitude.value
1785 lon = self.longitude.value
1786 elevation = value_or_none(self.elevation)
1787 return lat, lon, elevation
1790class Network(BaseNode):
1791 '''
1792 This type represents the Network layer, all station metadata is contained
1793 within this element. The official name of the network or other descriptive
1794 information can be included in the Description element. The Network can
1795 contain 0 or more Stations.
1796 '''
1798 total_number_stations = Counter.T(optional=True,
1799 xmltagname='TotalNumberStations')
1800 selected_number_stations = Counter.T(optional=True,
1801 xmltagname='SelectedNumberStations')
1802 station_list = List.T(Station.T(xmltagname='Station'))
1804 @property
1805 def station_code_list(self):
1806 return sorted(set(s.code for s in self.station_list))
1808 @property
1809 def sl_code_list(self):
1810 sls = set()
1811 for station in self.station_list:
1812 for channel in station.channel_list:
1813 sls.add((station.code, channel.location_code))
1815 return sorted(sls)
1817 def summary(self, width=80, indent=4):
1818 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1819 lines = ['%s (%i):' % (self.code, len(sls))]
1820 if sls:
1821 ssls = ['.'.join(x for x in c if x) for c in sls]
1822 w = max(len(x) for x in ssls)
1823 n = (width - indent) / (w+1)
1824 while ssls:
1825 lines.append(
1826 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1828 ssls[:n] = []
1830 return '\n'.join(lines)
1833def value_or_none(x):
1834 if x is not None:
1835 return x.value
1836 else:
1837 return None
1840def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1842 pos = lat, lon, elevation, depth = \
1843 channels[0].position_values
1845 if not all(pos == x.position_values for x in channels):
1846 info = '\n'.join(
1847 ' %s: %s' % (x.code, x.position_values) for
1848 x in channels)
1850 mess = 'encountered inconsistencies in channel ' \
1851 'lat/lon/elevation/depth ' \
1852 'for %s.%s.%s: \n%s' % (nsl + (info,))
1854 if inconsistencies == 'raise':
1855 raise InconsistentChannelLocations(mess)
1857 elif inconsistencies == 'warn':
1858 logger.warning(mess)
1859 logger.warning(' -> using mean values')
1861 apos = num.array([x.position_values for x in channels], dtype=float)
1862 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1863 / num.sum(num.isfinite(apos), axis=0)
1865 groups = {}
1866 for channel in channels:
1867 if channel.code not in groups:
1868 groups[channel.code] = []
1870 groups[channel.code].append(channel)
1872 pchannels = []
1873 for code in sorted(groups.keys()):
1874 data = [
1875 (channel.code, value_or_none(channel.azimuth),
1876 value_or_none(channel.dip))
1877 for channel in groups[code]]
1879 azimuth, dip = util.consistency_merge(
1880 data,
1881 message='channel orientation values differ:',
1882 error=inconsistencies)
1884 pchannels.append(
1885 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1887 return pyrocko.model.Station(
1888 *nsl,
1889 lat=mlat,
1890 lon=mlon,
1891 elevation=mele,
1892 depth=mdep,
1893 channels=pchannels)
1896class FDSNStationXML(Object):
1897 '''
1898 Top-level type for Station XML. Required field are Source (network ID of
1899 the institution sending the message) and one or more Network containers or
1900 one or more Station containers.
1901 '''
1903 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1904 source = String.T(xmltagname='Source')
1905 sender = String.T(optional=True, xmltagname='Sender')
1906 module = String.T(optional=True, xmltagname='Module')
1907 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1908 created = Timestamp.T(optional=True, xmltagname='Created')
1909 network_list = List.T(Network.T(xmltagname='Network'))
1911 xmltagname = 'FDSNStationXML'
1912 guessable_xmlns = [guts_xmlns]
1914 def __init__(self, *args, **kwargs):
1915 Object.__init__(self, *args, **kwargs)
1916 if self.created is None:
1917 self.created = time.time()
1919 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1920 time=None, timespan=None,
1921 inconsistencies='warn',
1922 sloppy=False):
1924 assert inconsistencies in ('raise', 'warn')
1926 if nslcs is not None:
1927 nslcs = set(nslcs)
1929 if nsls is not None:
1930 nsls = set(nsls)
1932 tt = ()
1933 if time is not None:
1934 tt = (time,)
1935 elif timespan is not None:
1936 tt = timespan
1938 ns_have = set()
1939 pstations = []
1940 sensor_to_channels = defaultdict(list)
1941 for network in self.network_list:
1942 if not network.spans(*tt):
1943 continue
1945 for station in network.station_list:
1946 if not station.spans(*tt):
1947 continue
1949 if station.channel_list:
1950 loc_to_channels = {}
1951 for channel in station.channel_list:
1952 if not channel.spans(*tt):
1953 continue
1955 loc = channel.location_code.strip()
1956 if loc not in loc_to_channels:
1957 loc_to_channels[loc] = []
1959 loc_to_channels[loc].append(channel)
1961 for loc in sorted(loc_to_channels.keys()):
1962 channels = loc_to_channels[loc]
1963 if nslcs is not None:
1964 channels = [channel for channel in channels
1965 if (network.code, station.code, loc,
1966 channel.code) in nslcs]
1968 if not channels:
1969 continue
1971 nsl = network.code, station.code, loc
1972 if nsls is not None and nsl not in nsls:
1973 continue
1975 for channel in channels:
1976 k = (nsl, channel.code[:-1])
1977 if not sloppy:
1978 k += (channel.start_date, channel.end_date)
1980 sensor_to_channels[k].append(channel)
1982 else:
1983 ns = (network.code, station.code)
1984 if ns in ns_have:
1985 message = 'Duplicate station ' \
1986 '(multiple epochs match): %s.%s ' % ns
1988 if inconsistencies == 'raise':
1989 raise Inconsistencies(message)
1990 else:
1991 logger.warning(message)
1993 ns_have.add(ns)
1995 pstations.append(pyrocko.model.Station(
1996 network.code, station.code, '*',
1997 lat=station.latitude.value,
1998 lon=station.longitude.value,
1999 elevation=value_or_none(station.elevation),
2000 name=station.description or ''))
2002 sensor_have = set()
2003 for k, channels in sensor_to_channels.items():
2004 (nsl, bi) = k[:2]
2005 if (nsl, bi) in sensor_have:
2006 message = 'Duplicate station ' \
2007 '(multiple epochs match): %s.%s.%s' % nsl
2009 if inconsistencies == 'raise':
2010 raise Inconsistencies(message)
2011 else:
2012 logger.warning(message)
2014 sensor_have.add((nsl, bi))
2016 pstations.append(
2017 pyrocko_station_from_channels(
2018 nsl,
2019 channels,
2020 inconsistencies=inconsistencies))
2022 return pstations
2024 @classmethod
2025 def from_pyrocko_stations(
2026 cls, pyrocko_stations, add_flat_responses_from=None):
2028 '''
2029 Generate :py:class:`FDSNStationXML` from list of
2030 :py:class;`pyrocko.model.Station` instances.
2032 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
2033 instances.
2034 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
2035 '''
2036 network_dict = defaultdict(list)
2038 if add_flat_responses_from:
2039 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
2040 extra = dict(
2041 response=Response(
2042 instrument_sensitivity=Sensitivity(
2043 value=1.0,
2044 frequency=1.0,
2045 input_units=Units(name=add_flat_responses_from))))
2046 else:
2047 extra = {}
2049 have_offsets = set()
2050 for s in pyrocko_stations:
2052 if s.north_shift != 0.0 or s.east_shift != 0.0:
2053 have_offsets.add(s.nsl())
2055 network, station, location = s.nsl()
2056 channel_list = []
2057 for c in s.channels:
2058 channel_list.append(
2059 Channel(
2060 location_code=location,
2061 code=c.name,
2062 latitude=Latitude(value=s.effective_lat),
2063 longitude=Longitude(value=s.effective_lon),
2064 elevation=Distance(value=s.elevation),
2065 depth=Distance(value=s.depth),
2066 azimuth=Azimuth(value=c.azimuth),
2067 dip=Dip(value=c.dip),
2068 **extra
2069 )
2070 )
2072 network_dict[network].append(
2073 Station(
2074 code=station,
2075 latitude=Latitude(value=s.effective_lat),
2076 longitude=Longitude(value=s.effective_lon),
2077 elevation=Distance(value=s.elevation),
2078 channel_list=channel_list)
2079 )
2081 if have_offsets:
2082 logger.warning(
2083 'StationXML does not support Cartesian offsets in '
2084 'coordinates. Storing effective lat/lon for stations: %s' %
2085 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
2087 timestamp = util.to_time_float(time.time())
2088 network_list = []
2089 for k, station_list in network_dict.items():
2091 network_list.append(
2092 Network(
2093 code=k, station_list=station_list,
2094 total_number_stations=len(station_list)))
2096 sxml = FDSNStationXML(
2097 source='from pyrocko stations list', created=timestamp,
2098 network_list=network_list)
2100 sxml.validate()
2101 return sxml
2103 def iter_network_stations(
2104 self, net=None, sta=None, time=None, timespan=None):
2106 tt = ()
2107 if time is not None:
2108 tt = (time,)
2109 elif timespan is not None:
2110 tt = timespan
2112 for network in self.network_list:
2113 if not network.spans(*tt) or (
2114 net is not None and network.code != net):
2115 continue
2117 for station in network.station_list:
2118 if not station.spans(*tt) or (
2119 sta is not None and station.code != sta):
2120 continue
2122 yield (network, station)
2124 def iter_network_station_channels(
2125 self, net=None, sta=None, loc=None, cha=None,
2126 time=None, timespan=None):
2128 if loc is not None:
2129 loc = loc.strip()
2131 tt = ()
2132 if time is not None:
2133 tt = (time,)
2134 elif timespan is not None:
2135 tt = timespan
2137 for network in self.network_list:
2138 if not network.spans(*tt) or (
2139 net is not None and network.code != net):
2140 continue
2142 for station in network.station_list:
2143 if not station.spans(*tt) or (
2144 sta is not None and station.code != sta):
2145 continue
2147 if station.channel_list:
2148 for channel in station.channel_list:
2149 if (not channel.spans(*tt) or
2150 (cha is not None and channel.code != cha) or
2151 (loc is not None and
2152 channel.location_code.strip() != loc)):
2153 continue
2155 yield (network, station, channel)
2157 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2158 time=None, timespan=None):
2160 groups = {}
2161 for network, station, channel in self.iter_network_station_channels(
2162 net, sta, loc, cha, time=time, timespan=timespan):
2164 net = network.code
2165 sta = station.code
2166 cha = channel.code
2167 loc = channel.location_code.strip()
2168 if len(cha) == 3:
2169 bic = cha[:2] # band and intrument code according to SEED
2170 elif len(cha) == 1:
2171 bic = ''
2172 else:
2173 bic = cha
2175 if channel.response and \
2176 channel.response.instrument_sensitivity and \
2177 channel.response.instrument_sensitivity.input_units:
2179 unit = sanitize_units(
2180 channel.response.instrument_sensitivity.input_units.name)
2181 else:
2182 unit = None
2184 bic = (bic, unit)
2186 k = net, sta, loc
2187 if k not in groups:
2188 groups[k] = {}
2190 if bic not in groups[k]:
2191 groups[k][bic] = []
2193 groups[k][bic].append(channel)
2195 for nsl, bic_to_channels in groups.items():
2196 bad_bics = []
2197 for bic, channels in bic_to_channels.items():
2198 sample_rates = []
2199 for channel in channels:
2200 sample_rates.append(channel.sample_rate.value)
2202 if not same(sample_rates):
2203 scs = ','.join(channel.code for channel in channels)
2204 srs = ', '.join('%e' % x for x in sample_rates)
2205 err = 'ignoring channels with inconsistent sampling ' + \
2206 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2208 logger.warning(err)
2209 bad_bics.append(bic)
2211 for bic in bad_bics:
2212 del bic_to_channels[bic]
2214 return groups
2216 def choose_channels(
2217 self,
2218 target_sample_rate=None,
2219 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2220 priority_units=['M/S', 'M/S**2'],
2221 priority_instrument_code=['H', 'L'],
2222 time=None,
2223 timespan=None):
2225 nslcs = {}
2226 for nsl, bic_to_channels in self.get_channel_groups(
2227 time=time, timespan=timespan).items():
2229 useful_bics = []
2230 for bic, channels in bic_to_channels.items():
2231 rate = channels[0].sample_rate.value
2233 if target_sample_rate is not None and \
2234 rate < target_sample_rate*0.99999:
2235 continue
2237 if len(bic[0]) == 2:
2238 if bic[0][0] not in priority_band_code:
2239 continue
2241 if bic[0][1] not in priority_instrument_code:
2242 continue
2244 unit = bic[1]
2246 prio_unit = len(priority_units)
2247 try:
2248 prio_unit = priority_units.index(unit)
2249 except ValueError:
2250 pass
2252 prio_inst = len(priority_instrument_code)
2253 prio_band = len(priority_band_code)
2254 if len(channels[0].code) == 3:
2255 try:
2256 prio_inst = priority_instrument_code.index(
2257 channels[0].code[1])
2258 except ValueError:
2259 pass
2261 try:
2262 prio_band = priority_band_code.index(
2263 channels[0].code[0])
2264 except ValueError:
2265 pass
2267 if target_sample_rate is None:
2268 rate = -rate
2270 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2271 prio_inst, bic))
2273 useful_bics.sort()
2275 for _, _, rate, _, _, bic in useful_bics:
2276 channels = sorted(
2277 bic_to_channels[bic],
2278 key=lambda channel: channel.code)
2280 if channels:
2281 for channel in channels:
2282 nslcs[nsl + (channel.code,)] = channel
2284 break
2286 return nslcs
2288 def get_pyrocko_response(
2289 self, nslc,
2290 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2292 net, sta, loc, cha = nslc
2293 resps = []
2294 for _, _, channel in self.iter_network_station_channels(
2295 net, sta, loc, cha, time=time, timespan=timespan):
2296 resp = channel.response
2297 if resp:
2298 resp.check_sample_rates(channel)
2299 resp.check_units()
2300 resps.append(resp.get_pyrocko_response(
2301 '.'.join(nslc),
2302 fake_input_units=fake_input_units,
2303 stages=stages).expect_one())
2305 if not resps:
2306 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2307 elif len(resps) > 1:
2308 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2310 return resps[0]
2312 @property
2313 def n_code_list(self):
2314 return sorted(set(x.code for x in self.network_list))
2316 @property
2317 def ns_code_list(self):
2318 nss = set()
2319 for network in self.network_list:
2320 for station in network.station_list:
2321 nss.add((network.code, station.code))
2323 return sorted(nss)
2325 @property
2326 def nsl_code_list(self):
2327 nsls = set()
2328 for network in self.network_list:
2329 for station in network.station_list:
2330 for channel in station.channel_list:
2331 nsls.add(
2332 (network.code, station.code, channel.location_code))
2334 return sorted(nsls)
2336 @property
2337 def nslc_code_list(self):
2338 nslcs = set()
2339 for network in self.network_list:
2340 for station in network.station_list:
2341 for channel in station.channel_list:
2342 nslcs.add(
2343 (network.code, station.code, channel.location_code,
2344 channel.code))
2346 return sorted(nslcs)
2348 def summary(self):
2349 lst = [
2350 'number of n codes: %i' % len(self.n_code_list),
2351 'number of ns codes: %i' % len(self.ns_code_list),
2352 'number of nsl codes: %i' % len(self.nsl_code_list),
2353 'number of nslc codes: %i' % len(self.nslc_code_list)
2354 ]
2355 return '\n'.join(lst)
2357 def summary_stages(self):
2358 data = []
2359 for network, station, channel in self.iter_network_station_channels():
2360 nslc = (network.code, station.code, channel.location_code,
2361 channel.code)
2363 stages = []
2364 in_units = '?'
2365 out_units = '?'
2366 if channel.response:
2367 sens = channel.response.instrument_sensitivity
2368 if sens:
2369 in_units = sanitize_units(sens.input_units.name)
2370 out_units = sanitize_units(sens.output_units.name)
2372 for stage in channel.response.stage_list:
2373 stages.append(stage.summary())
2375 data.append(
2376 (nslc, tts(channel.start_date), tts(channel.end_date),
2377 in_units, out_units, stages))
2379 data.sort()
2381 lst = []
2382 for nslc, stmin, stmax, in_units, out_units, stages in data:
2383 lst.append(' %s: %s - %s, %s -> %s' % (
2384 '.'.join(nslc), stmin, stmax, in_units, out_units))
2385 for stage in stages:
2386 lst.append(' %s' % stage)
2388 return '\n'.join(lst)
2390 def _check_overlaps(self):
2391 by_nslc = {}
2392 for network in self.network_list:
2393 for station in network.station_list:
2394 for channel in station.channel_list:
2395 nslc = (network.code, station.code, channel.location_code,
2396 channel.code)
2397 if nslc not in by_nslc:
2398 by_nslc[nslc] = []
2400 by_nslc[nslc].append(channel)
2402 errors = []
2403 for nslc, channels in by_nslc.items():
2404 errors.extend(check_overlaps('Channel', nslc, channels))
2406 return errors
2408 def check(self):
2409 errors = []
2410 for meth in [self._check_overlaps]:
2411 errors.extend(meth())
2413 if errors:
2414 raise Inconsistencies(
2415 'Inconsistencies found in StationXML:\n '
2416 + '\n '.join(errors))
2419def load_channel_table(stream):
2421 networks = {}
2422 stations = {}
2424 for line in stream:
2425 line = str(line.decode('ascii'))
2426 if line.startswith('#'):
2427 continue
2429 t = line.rstrip().split('|')
2431 if len(t) != 17:
2432 logger.warning('Invalid channel record: %s' % line)
2433 continue
2435 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2436 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2438 try:
2439 scale = float(scale)
2440 except ValueError:
2441 scale = None
2443 try:
2444 scale_freq = float(scale_freq)
2445 except ValueError:
2446 scale_freq = None
2448 try:
2449 depth = float(dep)
2450 except ValueError:
2451 depth = 0.0
2453 try:
2454 azi = float(azi)
2455 dip = float(dip)
2456 except ValueError:
2457 azi = None
2458 dip = None
2460 try:
2461 if net not in networks:
2462 network = Network(code=net)
2463 else:
2464 network = networks[net]
2466 if (net, sta) not in stations:
2467 station = Station(
2468 code=sta, latitude=lat, longitude=lon, elevation=ele)
2470 station.regularize()
2471 else:
2472 station = stations[net, sta]
2474 if scale:
2475 resp = Response(
2476 instrument_sensitivity=Sensitivity(
2477 value=scale,
2478 frequency=scale_freq,
2479 input_units=scale_units))
2480 else:
2481 resp = None
2483 channel = Channel(
2484 code=cha,
2485 location_code=loc.strip(),
2486 latitude=lat,
2487 longitude=lon,
2488 elevation=ele,
2489 depth=depth,
2490 azimuth=azi,
2491 dip=dip,
2492 sensor=Equipment(description=sens),
2493 response=resp,
2494 sample_rate=sample_rate,
2495 start_date=start_date,
2496 end_date=end_date or None)
2498 channel.regularize()
2500 except ValidationError:
2501 raise InvalidRecord(line)
2503 if net not in networks:
2504 networks[net] = network
2506 if (net, sta) not in stations:
2507 stations[net, sta] = station
2508 network.station_list.append(station)
2510 station.channel_list.append(channel)
2512 return FDSNStationXML(
2513 source='created from table input',
2514 created=time.time(),
2515 network_list=sorted(networks.values(), key=lambda x: x.code))
2518def primitive_merge(sxs):
2519 networks = []
2520 for sx in sxs:
2521 networks.extend(sx.network_list)
2523 return FDSNStationXML(
2524 source='merged from different sources',
2525 created=time.time(),
2526 network_list=copy.deepcopy(
2527 sorted(networks, key=lambda x: x.code)))
2530def split_channels(sx):
2531 for nslc in sx.nslc_code_list:
2532 network_list = sx.network_list
2533 network_list_filtered = [
2534 network for network in network_list
2535 if network.code == nslc[0]]
2537 for network in network_list_filtered:
2538 sx.network_list = [network]
2539 station_list = network.station_list
2540 station_list_filtered = [
2541 station for station in station_list
2542 if station.code == nslc[1]]
2544 for station in station_list_filtered:
2545 network.station_list = [station]
2546 channel_list = station.channel_list
2547 station.channel_list = [
2548 channel for channel in channel_list
2549 if (channel.location_code, channel.code) == nslc[2:4]]
2551 yield nslc, copy.deepcopy(sx)
2553 station.channel_list = channel_list
2555 network.station_list = station_list
2557 sx.network_list = network_list
2560if __name__ == '__main__':
2561 from optparse import OptionParser
2563 util.setup_logging('pyrocko.io.stationxml', 'warning')
2565 usage = \
2566 'python -m pyrocko.io.stationxml check|stats|stages ' \
2567 '<filename> [options]'
2569 description = '''Torture StationXML file.'''
2571 parser = OptionParser(
2572 usage=usage,
2573 description=description,
2574 formatter=util.BetterHelpFormatter())
2576 (options, args) = parser.parse_args(sys.argv[1:])
2578 if len(args) != 2:
2579 parser.print_help()
2580 sys.exit(1)
2582 action, path = args
2584 sx = load_xml(filename=path)
2585 if action == 'check':
2586 try:
2587 sx.check()
2588 except Inconsistencies as e:
2589 logger.error(e)
2590 sys.exit(1)
2592 elif action == 'stats':
2593 print(sx.summary())
2595 elif action == 'stages':
2596 print(sx.summary_stages())
2598 else:
2599 parser.print_help()
2600 sys.exit('unknown action: %s' % action)