1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import sys
8import time
9import logging
10import datetime
11import calendar
12import math
13import copy
15import numpy as num
17from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String,
18 Unicode, Int, Float, List, Object, Timestamp,
19 ValidationError, TBase, re_tz, Any, Tuple)
20from pyrocko.guts import load_xml # noqa
21from pyrocko.util import hpfloat, time_to_str, get_time_float
23import pyrocko.model
24from pyrocko import util, response
26try:
27 newstr = unicode
28except NameError:
29 newstr = str
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}
49class StationXMLError(Exception):
50 pass
53class Inconsistencies(StationXMLError):
54 pass
57class NoResponseInformation(StationXMLError):
58 pass
61class MultipleResponseInformation(StationXMLError):
62 pass
65class InconsistentResponseInformation(StationXMLError):
66 pass
69class InconsistentChannelLocations(StationXMLError):
70 pass
73class InvalidRecord(StationXMLError):
74 def __init__(self, line):
75 StationXMLError.__init__(self)
76 self._line = line
78 def __str__(self):
79 return 'Invalid record: "%s"' % self._line
82_exceptions = dict(
83 Inconsistencies=Inconsistencies,
84 NoResponseInformation=NoResponseInformation,
85 MultipleResponseInformation=MultipleResponseInformation,
86 InconsistentResponseInformation=InconsistentResponseInformation,
87 InconsistentChannelLocations=InconsistentChannelLocations,
88 InvalidRecord=InvalidRecord,
89 ValueError=ValueError)
92_logs = dict(
93 info=logger.info,
94 warning=logger.warning,
95 error=logger.error)
98class DeliveryError(StationXMLError):
99 pass
102class Delivery(Object):
104 def __init__(self, payload=None, log=None, errors=None, error=None):
105 if payload is None:
106 payload = []
108 if log is None:
109 log = []
111 if errors is None:
112 errors = []
114 if error is not None:
115 errors.append(error)
117 Object.__init__(self, payload=payload, log=log, errors=errors)
119 payload = List.T(Any.T())
120 log = List.T(Tuple.T(3, String.T()))
121 errors = List.T(Tuple.T(3, String.T()))
123 def extend(self, other):
124 self.payload.extend(other.payload)
125 self.log.extend(other.log)
126 self.errors.extend(other.errors)
128 def extend_without_payload(self, other):
129 self.log.extend(other.log)
130 self.errors.extend(other.errors)
131 return other.payload
133 def emit_log(self):
134 for name, message, context in self.log:
135 message = '%s: %s' % (context, message)
136 _logs[name](message)
138 def expect(self, quiet=False):
139 if not quiet:
140 self.emit_log()
142 if self.errors:
143 name, message, context = self.errors[0]
144 if context:
145 message += ' (%s)' % context
147 if len(self.errors) > 1:
148 message += ' Additional errors pending.'
150 raise _exceptions[name](message)
152 return self.payload
154 def expect_one(self, quiet=False):
155 payload = self.expect(quiet=quiet)
156 if len(payload) != 1:
157 raise DeliveryError(
158 'Expected 1 element but got %i.' % len(payload))
160 return payload[0]
163def wrap(s, width=80, indent=4):
164 words = s.split()
165 lines = []
166 t = []
167 n = 0
168 for w in words:
169 if n + len(w) >= width:
170 lines.append(' '.join(t))
171 n = indent
172 t = [' '*(indent-1)]
174 t.append(w)
175 n += len(w) + 1
177 lines.append(' '.join(t))
178 return '\n'.join(lines)
181def same(x, eps=0.0):
182 if any(type(x[0]) != type(r) for r in x):
183 return False
185 if isinstance(x[0], float):
186 return all(abs(r-x[0]) <= eps for r in x)
187 else:
188 return all(r == x[0] for r in x)
191def same_sample_rate(a, b, eps=1.0e-6):
192 return abs(a - b) < min(a, b)*eps
195def evaluate1(resp, f):
196 return resp.evaluate(num.array([f], dtype=float))[0]
199def check_resp(resp, value, frequency, limit_db, prelude, context):
201 try:
202 value_resp = num.abs(evaluate1(resp, frequency))
203 except response.InvalidResponseError as e:
204 return Delivery(log=[(
205 'warning',
206 'Could not check response: %s' % str(e),
207 context)])
209 if value_resp == 0.0:
210 return Delivery(log=[(
211 'warning',
212 '%s\n'
213 ' computed response is zero' % prelude,
214 context)])
216 diff_db = 20.0 * num.log10(value_resp/value)
218 if num.abs(diff_db) > limit_db:
219 return Delivery(log=[(
220 'warning',
221 '%s\n'
222 ' reported value: %g\n'
223 ' computed value: %g\n'
224 ' at frequency [Hz]: %g\n'
225 ' factor: %g\n'
226 ' difference [dB]: %g\n'
227 ' limit [dB]: %g' % (
228 prelude,
229 value,
230 value_resp,
231 frequency,
232 value_resp/value,
233 diff_db,
234 limit_db),
235 context)])
237 return Delivery()
240def tts(t):
241 if t is None:
242 return '?'
243 else:
244 return util.tts(t, format='%Y-%m-%d')
247this_year = time.gmtime()[0]
250class DummyAwareOptionalTimestamp(Object):
251 '''
252 Optional timestamp with support for some common placeholder values.
254 Some StationXML files contain arbitrary placeholder values for open end
255 intervals, like "2100-01-01". Depending on the time range supported by the
256 system, these dates are translated into ``None`` to prevent crashes with
257 this type.
258 '''
259 dummy_for = (hpfloat, float)
260 dummy_for_description = 'time_float'
262 class __T(TBase):
264 def regularize_extra(self, val):
265 time_float = get_time_float()
267 if isinstance(val, datetime.datetime):
268 tt = val.utctimetuple()
269 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
271 elif isinstance(val, datetime.date):
272 tt = val.timetuple()
273 val = time_float(calendar.timegm(tt))
275 elif isinstance(val, (str, newstr)):
276 val = val.strip()
278 tz_offset = 0
280 m = re_tz.search(val)
281 if m:
282 sh = m.group(2)
283 sm = m.group(4)
284 tz_offset = (int(sh)*3600 if sh else 0) \
285 + (int(sm)*60 if sm else 0)
287 val = re_tz.sub('', val)
289 if len(val) > 10 and val[10] == 'T':
290 val = val.replace('T', ' ', 1)
292 try:
293 val = util.str_to_time(val) - tz_offset
295 except util.TimeStrError:
296 year = int(val[:4])
297 if sys.maxsize > 2**32: # if we're on 64bit
298 if year > this_year + 100:
299 return None # StationXML contained a dummy date
301 if year < 1903: # for macOS, 1900-01-01 dummy dates
302 return None
304 else: # 32bit end of time is in 2038
305 if this_year < 2037 and year > 2037 or year < 1903:
306 return None # StationXML contained a dummy date
308 raise
310 elif isinstance(val, (int, float)):
311 val = time_float(val)
313 else:
314 raise ValidationError(
315 '%s: cannot convert "%s" to type %s' % (
316 self.xname(), val, time_float))
318 return val
320 def to_save(self, val):
321 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
322 .rstrip('0').rstrip('.')
324 def to_save_xml(self, val):
325 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
326 .rstrip('0').rstrip('.') + 'Z'
329class Nominal(StringChoice):
330 choices = [
331 'NOMINAL',
332 'CALCULATED']
335class Email(UnicodePattern):
336 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
339class RestrictedStatus(StringChoice):
340 choices = [
341 'open',
342 'closed',
343 'partial']
346class Type(StringChoice):
347 choices = [
348 'TRIGGERED',
349 'CONTINUOUS',
350 'HEALTH',
351 'GEOPHYSICAL',
352 'WEATHER',
353 'FLAG',
354 'SYNTHESIZED',
355 'INPUT',
356 'EXPERIMENTAL',
357 'MAINTENANCE',
358 'BEAM']
360 class __T(StringChoice.T):
361 def validate_extra(self, val):
362 if val not in self.choices:
363 logger.warning(
364 'channel type: "%s" is not a valid choice out of %s' %
365 (val, repr(self.choices)))
368class PzTransferFunction(StringChoice):
369 choices = [
370 'LAPLACE (RADIANS/SECOND)',
371 'LAPLACE (HERTZ)',
372 'DIGITAL (Z-TRANSFORM)']
375class Symmetry(StringChoice):
376 choices = [
377 'NONE',
378 'EVEN',
379 'ODD']
382class CfTransferFunction(StringChoice):
384 class __T(StringChoice.T):
385 def validate(self, val, regularize=False, depth=-1):
386 if regularize:
387 try:
388 val = str(val)
389 except ValueError:
390 raise ValidationError(
391 '%s: cannot convert to string %s' % (self.xname,
392 repr(val)))
394 val = self.dummy_cls.replacements.get(val, val)
396 self.validate_extra(val)
397 return val
399 choices = [
400 'ANALOG (RADIANS/SECOND)',
401 'ANALOG (HERTZ)',
402 'DIGITAL']
404 replacements = {
405 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
406 'ANALOG (HZ)': 'ANALOG (HERTZ)',
407 }
410class Approximation(StringChoice):
411 choices = [
412 'MACLAURIN']
415class PhoneNumber(StringPattern):
416 pattern = '[0-9]+-[0-9]+'
419class Site(Object):
420 '''
421 Description of a site location using name and optional geopolitical
422 boundaries (country, city, etc.).
423 '''
425 name = Unicode.T(default='', xmltagname='Name')
426 description = Unicode.T(optional=True, xmltagname='Description')
427 town = Unicode.T(optional=True, xmltagname='Town')
428 county = Unicode.T(optional=True, xmltagname='County')
429 region = Unicode.T(optional=True, xmltagname='Region')
430 country = Unicode.T(optional=True, xmltagname='Country')
433class ExternalReference(Object):
434 '''
435 This type contains a URI and description for external data that users may
436 want to reference in StationXML.
437 '''
439 uri = String.T(xmltagname='URI')
440 description = Unicode.T(xmltagname='Description')
443class Units(Object):
444 '''
445 A type to document units. Corresponds to SEED blockette 34.
446 '''
448 def __init__(self, name=None, **kwargs):
449 Object.__init__(self, name=name, **kwargs)
451 name = String.T(xmltagname='Name')
452 description = Unicode.T(optional=True, xmltagname='Description')
455class Counter(Int):
456 pass
459class SampleRateRatio(Object):
460 '''
461 Sample rate expressed as number of samples in a number of seconds.
462 '''
464 number_samples = Int.T(xmltagname='NumberSamples')
465 number_seconds = Int.T(xmltagname='NumberSeconds')
468class Gain(Object):
469 '''
470 Complex type for sensitivity and frequency ranges. This complex type can be
471 used to represent both overall sensitivities and individual stage gains.
472 The FrequencyRangeGroup is an optional construct that defines a pass band
473 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
474 is valid within the number of decibels specified in FrequencyDBVariation.
475 '''
477 def __init__(self, value=None, **kwargs):
478 Object.__init__(self, value=value, **kwargs)
480 value = Float.T(optional=True, xmltagname='Value')
481 frequency = Float.T(optional=True, xmltagname='Frequency')
483 def summary(self):
484 return 'gain(%g @ %g)' % (self.value, self.frequency)
487class NumeratorCoefficient(Object):
488 i = Int.T(optional=True, xmlstyle='attribute')
489 value = Float.T(xmlstyle='content')
492class FloatNoUnit(Object):
493 def __init__(self, value=None, **kwargs):
494 Object.__init__(self, value=value, **kwargs)
496 plus_error = Float.T(optional=True, xmlstyle='attribute')
497 minus_error = Float.T(optional=True, xmlstyle='attribute')
498 value = Float.T(xmlstyle='content')
501class FloatWithUnit(FloatNoUnit):
502 unit = String.T(optional=True, xmlstyle='attribute')
505class Equipment(Object):
506 resource_id = String.T(optional=True, xmlstyle='attribute')
507 type = String.T(optional=True, xmltagname='Type')
508 description = Unicode.T(optional=True, xmltagname='Description')
509 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
510 vendor = Unicode.T(optional=True, xmltagname='Vendor')
511 model = Unicode.T(optional=True, xmltagname='Model')
512 serial_number = String.T(optional=True, xmltagname='SerialNumber')
513 installation_date = DummyAwareOptionalTimestamp.T(
514 optional=True,
515 xmltagname='InstallationDate')
516 removal_date = DummyAwareOptionalTimestamp.T(
517 optional=True,
518 xmltagname='RemovalDate')
519 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
522class PhoneNumber(Object):
523 description = Unicode.T(optional=True, xmlstyle='attribute')
524 country_code = Int.T(optional=True, xmltagname='CountryCode')
525 area_code = Int.T(xmltagname='AreaCode')
526 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
529class BaseFilter(Object):
530 '''
531 The BaseFilter is derived by all filters.
532 '''
534 resource_id = String.T(optional=True, xmlstyle='attribute')
535 name = String.T(optional=True, xmlstyle='attribute')
536 description = Unicode.T(optional=True, xmltagname='Description')
537 input_units = Units.T(optional=True, xmltagname='InputUnits')
538 output_units = Units.T(optional=True, xmltagname='OutputUnits')
541class Sensitivity(Gain):
542 '''
543 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
544 construct that defines a pass band in Hertz (FrequencyStart and
545 FrequencyEnd) in which the SensitivityValue is valid within the number of
546 decibels specified in FrequencyDBVariation.
547 '''
549 input_units = Units.T(optional=True, xmltagname='InputUnits')
550 output_units = Units.T(optional=True, xmltagname='OutputUnits')
551 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
552 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
553 frequency_db_variation = Float.T(optional=True,
554 xmltagname='FrequencyDBVariation')
556 def get_pyrocko_response(self):
557 return Delivery(
558 [response.PoleZeroResponse(constant=self.value)])
561class Coefficient(FloatNoUnit):
562 number = Counter.T(optional=True, xmlstyle='attribute')
565class PoleZero(Object):
566 '''
567 Complex numbers used as poles or zeros in channel response.
568 '''
570 number = Int.T(optional=True, xmlstyle='attribute')
571 real = FloatNoUnit.T(xmltagname='Real')
572 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
574 def value(self):
575 return self.real.value + 1J * self.imaginary.value
578class ClockDrift(FloatWithUnit):
579 unit = String.T(default='SECONDS/SAMPLE', optional=True,
580 xmlstyle='attribute') # fixed
583class Second(FloatWithUnit):
584 '''
585 A time value in seconds.
586 '''
588 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
589 # fixed unit
592class Voltage(FloatWithUnit):
593 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
594 # fixed unit
597class Angle(FloatWithUnit):
598 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
599 # fixed unit
602class Azimuth(FloatWithUnit):
603 '''
604 Instrument azimuth, degrees clockwise from North.
605 '''
607 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
608 # fixed unit
611class Dip(FloatWithUnit):
612 '''
613 Instrument dip in degrees down from horizontal. Together azimuth and dip
614 describe the direction of the sensitive axis of the instrument.
615 '''
617 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
618 # fixed unit
621class Distance(FloatWithUnit):
622 '''
623 Extension of FloatWithUnit for distances, elevations, and depths.
624 '''
626 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
627 # NOT fixed unit!
630class Frequency(FloatWithUnit):
631 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
632 # fixed unit
635class SampleRate(FloatWithUnit):
636 '''
637 Sample rate in samples per second.
638 '''
640 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
641 # fixed unit
644class Person(Object):
645 '''
646 Representation of a person's contact information. A person can belong to
647 multiple agencies and have multiple email addresses and phone numbers.
648 '''
650 name_list = List.T(Unicode.T(xmltagname='Name'))
651 agency_list = List.T(Unicode.T(xmltagname='Agency'))
652 email_list = List.T(Email.T(xmltagname='Email'))
653 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
656class FIR(BaseFilter):
657 '''
658 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
659 also commonly documented using the Coefficients element.
660 '''
662 symmetry = Symmetry.T(xmltagname='Symmetry')
663 numerator_coefficient_list = List.T(
664 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
666 def summary(self):
667 return 'fir(%i%s)' % (
668 self.get_ncoefficiencs(),
669 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
671 def get_effective_coefficients(self):
672 b = num.array(
673 [v.value for v in self.numerator_coefficient_list],
674 dtype=float)
676 if self.symmetry == 'ODD':
677 b = num.concatenate((b, b[-2::-1]))
678 elif self.symmetry == 'EVEN':
679 b = num.concatenate((b, b[::-1]))
681 return b
683 def get_effective_symmetry(self):
684 if self.symmetry == 'NONE':
685 b = self.get_effective_coefficients()
686 if num.all(b - b[::-1] == 0):
687 return ['EVEN', 'ODD'][b.size % 2]
689 return self.symmetry
691 def get_ncoefficiencs(self):
692 nf = len(self.numerator_coefficient_list)
693 if self.symmetry == 'ODD':
694 nc = nf * 2 + 1
695 elif self.symmetry == 'EVEN':
696 nc = nf * 2
697 else:
698 nc = nf
700 return nc
702 def estimate_delay(self, deltat):
703 nc = self.get_ncoefficiencs()
704 if nc > 0:
705 return deltat * (nc - 1) / 2.0
706 else:
707 return 0.0
709 def get_pyrocko_response(
710 self, context, deltat, delay_responses, normalization_frequency):
712 context += self.summary()
714 if not self.numerator_coefficient_list:
715 return Delivery([])
717 b = self.get_effective_coefficients()
719 log = []
720 drop_phase = self.get_effective_symmetry() != 'NONE'
722 if not deltat:
723 log.append((
724 'error',
725 'Digital response requires knowledge about sampling '
726 'interval. Response will be unusable.',
727 context))
729 resp = response.DigitalFilterResponse(
730 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
732 if normalization_frequency is not None and deltat is not None:
733 normalization_frequency = 0.0
734 normalization = num.abs(evaluate1(resp, normalization_frequency))
736 if num.abs(normalization - 1.0) > 1e-2:
737 log.append((
738 'warning',
739 'FIR filter coefficients are not normalized. Normalizing '
740 'them. Factor: %g' % normalization, context))
742 resp = response.DigitalFilterResponse(
743 (b/normalization).tolist(), [1.0], deltat,
744 drop_phase=drop_phase)
746 resps = [resp]
748 if not drop_phase:
749 resps.extend(delay_responses)
751 return Delivery(resps, log=log)
754class Coefficients(BaseFilter):
755 '''
756 Response: coefficients for FIR filter. Laplace transforms or IIR filters
757 can be expressed using type as well but the PolesAndZeros should be used
758 instead. Corresponds to SEED blockette 54.
759 '''
761 cf_transfer_function_type = CfTransferFunction.T(
762 xmltagname='CfTransferFunctionType')
763 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
764 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
766 def summary(self):
767 return 'coef_%s(%i,%i%s)' % (
768 'ABC?'[
769 CfTransferFunction.choices.index(
770 self.cf_transfer_function_type)],
771 len(self.numerator_list),
772 len(self.denominator_list),
773 ',sym' if self.is_symmetric_fir else '')
775 def estimate_delay(self, deltat):
776 nc = len(self.numerator_list)
777 if nc > 0:
778 return deltat * (len(self.numerator_list) - 1) / 2.0
779 else:
780 return 0.0
782 def is_symmetric_fir(self):
783 if len(self.denominator_list) != 0:
784 return False
785 b = [v.value for v in self.numerator_list]
786 return b == b[::-1]
788 def get_pyrocko_response(
789 self, context, deltat, delay_responses, normalization_frequency):
791 context += self.summary()
793 factor = 1.0
794 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
795 factor = 2.0*math.pi
797 if not self.numerator_list and not self.denominator_list:
798 return Delivery(payload=[])
800 b = num.array(
801 [v.value*factor for v in self.numerator_list], dtype=float)
803 a = num.array(
804 [1.0] + [v.value*factor for v in self.denominator_list],
805 dtype=float)
807 log = []
808 resps = []
809 if self.cf_transfer_function_type in [
810 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
811 resps.append(response.AnalogFilterResponse(b, a))
813 elif self.cf_transfer_function_type == 'DIGITAL':
814 if not deltat:
815 log.append((
816 'error',
817 'Digital response requires knowledge about sampling '
818 'interval. Response will be unusable.',
819 context))
821 drop_phase = self.is_symmetric_fir()
822 resp = response.DigitalFilterResponse(
823 b, a, deltat or 0.0, drop_phase=drop_phase)
825 if normalization_frequency is not None and deltat is not None:
826 normalization = num.abs(evaluate1(
827 resp, normalization_frequency))
829 if num.abs(normalization - 1.0) > 1e-2:
830 log.append((
831 'warning',
832 'FIR filter coefficients are not normalized. '
833 'Normalizing them. Factor: %g' % normalization,
834 context))
836 resp = response.DigitalFilterResponse(
837 (b/normalization).tolist(), [1.0], deltat,
838 drop_phase=drop_phase)
840 resps.append(resp)
842 if not drop_phase:
843 resps.extend(delay_responses)
845 else:
846 return Delivery(error=(
847 'ValueError',
848 'Unknown transfer function type: %s' % (
849 self.cf_transfer_function_type)))
851 return Delivery(payload=resps, log=log)
854class Latitude(FloatWithUnit):
855 '''
856 Type for latitude coordinate.
857 '''
859 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
860 # fixed unit
861 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
864class Longitude(FloatWithUnit):
865 '''
866 Type for longitude coordinate.
867 '''
869 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
870 # fixed unit
871 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
874class PolesZeros(BaseFilter):
875 '''
876 Response: complex poles and zeros. Corresponds to SEED blockette 53.
877 '''
879 pz_transfer_function_type = PzTransferFunction.T(
880 xmltagname='PzTransferFunctionType')
881 normalization_factor = Float.T(default=1.0,
882 xmltagname='NormalizationFactor')
883 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
884 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
885 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
887 def summary(self):
888 return 'pz_%s(%i,%i)' % (
889 'ABC?'[
890 PzTransferFunction.choices.index(
891 self.pz_transfer_function_type)],
892 len(self.pole_list),
893 len(self.zero_list))
895 def get_pyrocko_response(self, context='', deltat=None):
897 context += self.summary()
899 factor = 1.0
900 cfactor = 1.0
901 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
902 factor = 2. * math.pi
903 cfactor = (2. * math.pi)**(
904 len(self.pole_list) - len(self.zero_list))
906 log = []
907 if self.normalization_factor is None \
908 or self.normalization_factor == 0.0:
910 log.append((
911 'warning',
912 'No pole-zero normalization factor given. '
913 'Assuming a value of 1.0',
914 context))
916 nfactor = 1.0
917 else:
918 nfactor = self.normalization_factor
920 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
921 if not is_digital:
922 resp = response.PoleZeroResponse(
923 constant=nfactor*cfactor,
924 zeros=[z.value()*factor for z in self.zero_list],
925 poles=[p.value()*factor for p in self.pole_list])
926 else:
927 if not deltat:
928 log.append((
929 'error',
930 'Digital response requires knowledge about sampling '
931 'interval. Response will be unusable.',
932 context))
934 resp = response.DigitalPoleZeroResponse(
935 constant=nfactor*cfactor,
936 zeros=[z.value()*factor for z in self.zero_list],
937 poles=[p.value()*factor for p in self.pole_list],
938 deltat=deltat or 0.0)
940 if not self.normalization_frequency.value:
941 log.append((
942 'warning',
943 'Cannot check pole-zero normalization factor: '
944 'No normalization frequency given.',
945 context))
947 else:
948 if is_digital and not deltat:
949 log.append((
950 'warning',
951 'Cannot check computed vs reported normalization '
952 'factor without knowing the sampling interval.',
953 context))
954 else:
955 computed_normalization_factor = nfactor / abs(evaluate1(
956 resp, self.normalization_frequency.value))
958 db = 20.0 * num.log10(
959 computed_normalization_factor / nfactor)
961 if abs(db) > 0.17:
962 log.append((
963 'warning',
964 'Computed and reported normalization factors differ '
965 'by %g dB: computed: %g, reported: %g' % (
966 db,
967 computed_normalization_factor,
968 nfactor),
969 context))
971 return Delivery([resp], log)
974class ResponseListElement(Object):
975 frequency = Frequency.T(xmltagname='Frequency')
976 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
977 phase = Angle.T(xmltagname='Phase')
980class Polynomial(BaseFilter):
981 '''
982 Response: expressed as a polynomial (allows non-linear sensors to be
983 described). Corresponds to SEED blockette 62. Can be used to describe a
984 stage of acquisition or a complete system.
985 '''
987 approximation_type = Approximation.T(default='MACLAURIN',
988 xmltagname='ApproximationType')
989 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
990 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
991 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
992 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
993 maximum_error = Float.T(xmltagname='MaximumError')
994 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
996 def summary(self):
997 return 'poly(%i)' % len(self.coefficient_list)
1000class Decimation(Object):
1001 '''
1002 Corresponds to SEED blockette 57.
1003 '''
1005 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1006 factor = Int.T(xmltagname='Factor')
1007 offset = Int.T(xmltagname='Offset')
1008 delay = FloatWithUnit.T(xmltagname='Delay')
1009 correction = FloatWithUnit.T(xmltagname='Correction')
1011 def summary(self):
1012 return 'deci(%i, %g -> %g, %g)' % (
1013 self.factor,
1014 self.input_sample_rate.value,
1015 self.input_sample_rate.value / self.factor,
1016 self.delay.value)
1018 def get_pyrocko_response(self):
1019 if self.delay and self.delay.value != 0.0:
1020 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1021 else:
1022 return Delivery([])
1025class Operator(Object):
1026 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1027 contact_list = List.T(Person.T(xmltagname='Contact'))
1028 web_site = String.T(optional=True, xmltagname='WebSite')
1031class Comment(Object):
1032 '''
1033 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1034 and 59.
1035 '''
1037 id = Counter.T(optional=True, xmlstyle='attribute')
1038 value = Unicode.T(xmltagname='Value')
1039 begin_effective_time = DummyAwareOptionalTimestamp.T(
1040 optional=True,
1041 xmltagname='BeginEffectiveTime')
1042 end_effective_time = DummyAwareOptionalTimestamp.T(
1043 optional=True,
1044 xmltagname='EndEffectiveTime')
1045 author_list = List.T(Person.T(xmltagname='Author'))
1048class ResponseList(BaseFilter):
1049 '''
1050 Response: list of frequency, amplitude and phase values. Corresponds to
1051 SEED blockette 55.
1052 '''
1054 response_list_element_list = List.T(
1055 ResponseListElement.T(xmltagname='ResponseListElement'))
1057 def summary(self):
1058 return 'list(%i)' % len(self.response_list_element_list)
1061class Log(Object):
1062 '''
1063 Container for log entries.
1064 '''
1066 entry_list = List.T(Comment.T(xmltagname='Entry'))
1069class ResponseStage(Object):
1070 '''
1071 This complex type represents channel response and covers SEED blockettes 53
1072 to 56.
1073 '''
1075 number = Counter.T(xmlstyle='attribute')
1076 resource_id = String.T(optional=True, xmlstyle='attribute')
1077 poles_zeros_list = List.T(
1078 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1079 coefficients_list = List.T(
1080 Coefficients.T(optional=True, xmltagname='Coefficients'))
1081 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1082 fir = FIR.T(optional=True, xmltagname='FIR')
1083 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1084 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1085 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1087 def summary(self):
1088 elements = []
1090 for stuff in [
1091 self.poles_zeros_list, self.coefficients_list,
1092 self.response_list, self.fir, self.polynomial,
1093 self.decimation, self.stage_gain]:
1095 if stuff:
1096 if isinstance(stuff, Object):
1097 elements.append(stuff.summary())
1098 else:
1099 elements.extend(obj.summary() for obj in stuff)
1101 return '%i: %s %s -> %s' % (
1102 self.number,
1103 ', '.join(elements),
1104 self.input_units.name.upper() if self.input_units else '?',
1105 self.output_units.name.upper() if self.output_units else '?')
1107 def get_squirrel_response_stage(self, context):
1108 from pyrocko.squirrel.model import ResponseStage
1110 delivery = Delivery()
1111 delivery_pr = self.get_pyrocko_response(context)
1112 log = delivery_pr.log
1113 delivery_pr.log = []
1114 elements = delivery.extend_without_payload(delivery_pr)
1116 def to_quantity(unit):
1117 trans = {
1118 'M/S': 'velocity',
1119 'M': 'displacement',
1120 'M/S**2': 'acceleration',
1121 'V': 'voltage',
1122 'COUNTS': 'counts',
1123 'COUNT': 'counts',
1124 'PA': 'pressure',
1125 'RAD/S': 'rotation'}
1127 if unit is None:
1128 return None
1130 name = unit.name.upper()
1131 if name in trans:
1132 return trans[name]
1133 else:
1134 delivery.log.append((
1135 'warning',
1136 'Units not supported by Squirrel framework: %s' % (
1137 unit.name.upper() + (
1138 ' (%s)' % unit.description
1139 if unit.description else '')),
1140 context))
1142 return 'unsupported_quantity(%s)' % unit
1144 delivery.payload = [ResponseStage(
1145 input_quantity=to_quantity(self.input_units),
1146 output_quantity=to_quantity(self.output_units),
1147 input_sample_rate=self.input_sample_rate,
1148 output_sample_rate=self.output_sample_rate,
1149 elements=elements,
1150 log=log)]
1152 return delivery
1154 def get_pyrocko_response(self, context, gain_only=False):
1156 context = context + ', stage %i' % self.number
1158 responses = []
1159 log = []
1160 if self.stage_gain:
1161 normalization_frequency = self.stage_gain.frequency or 0.0
1162 else:
1163 normalization_frequency = 0.0
1165 if not gain_only:
1166 deltat = None
1167 delay_responses = []
1168 if self.decimation:
1169 rate = self.decimation.input_sample_rate.value
1170 if rate > 0.0:
1171 deltat = 1.0 / rate
1172 delivery = self.decimation.get_pyrocko_response()
1173 if delivery.errors:
1174 return delivery
1176 delay_responses = delivery.payload
1177 log.extend(delivery.log)
1179 for pzs in self.poles_zeros_list:
1180 delivery = pzs.get_pyrocko_response(context, deltat)
1181 if delivery.errors:
1182 return delivery
1184 pz_resps = delivery.payload
1185 log.extend(delivery.log)
1186 responses.extend(pz_resps)
1188 # emulate incorrect? evalresp behaviour
1189 if pzs.normalization_frequency != normalization_frequency \
1190 and normalization_frequency != 0.0:
1192 try:
1193 trial = response.MultiplyResponse(pz_resps)
1194 anorm = num.abs(evaluate1(
1195 trial, pzs.normalization_frequency.value))
1196 asens = num.abs(
1197 evaluate1(trial, normalization_frequency))
1199 factor = anorm/asens
1201 if abs(factor - 1.0) > 0.01:
1202 log.append((
1203 'warning',
1204 'PZ normalization frequency (%g) is different '
1205 'from stage gain frequency (%s) -> Emulating '
1206 'possibly incorrect evalresp behaviour. '
1207 'Correction factor: %g' % (
1208 pzs.normalization_frequency.value,
1209 normalization_frequency,
1210 factor),
1211 context))
1213 responses.append(
1214 response.PoleZeroResponse(constant=factor))
1215 except response.InvalidResponseError as e:
1216 log.append((
1217 'warning',
1218 'Could not check response: %s' % str(e),
1219 context))
1221 if len(self.poles_zeros_list) > 1:
1222 log.append((
1223 'warning',
1224 'Multiple poles and zeros records in single response '
1225 'stage.',
1226 context))
1228 for cfs in self.coefficients_list + (
1229 [self.fir] if self.fir else []):
1231 delivery = cfs.get_pyrocko_response(
1232 context, deltat, delay_responses,
1233 normalization_frequency)
1235 if delivery.errors:
1236 return delivery
1238 responses.extend(delivery.payload)
1239 log.extend(delivery.log)
1241 if len(self.coefficients_list) > 1:
1242 log.append((
1243 'warning',
1244 'Multiple filter coefficients lists in single response '
1245 'stage.',
1246 context))
1248 if self.response_list:
1249 log.append((
1250 'warning',
1251 'Unhandled response element of type: ResponseList',
1252 context))
1254 if self.polynomial:
1255 log.append((
1256 'warning',
1257 'Unhandled response element of type: Polynomial',
1258 context))
1260 if self.stage_gain:
1261 responses.append(
1262 response.PoleZeroResponse(constant=self.stage_gain.value))
1264 return Delivery(responses, log)
1266 @property
1267 def input_units(self):
1268 for e in (self.poles_zeros_list + self.coefficients_list +
1269 [self.response_list, self.fir, self.polynomial]):
1270 if e is not None:
1271 return e.input_units
1273 return None
1275 @property
1276 def output_units(self):
1277 for e in (self.poles_zeros_list + self.coefficients_list +
1278 [self.response_list, self.fir, self.polynomial]):
1279 if e is not None:
1280 return e.output_units
1282 return None
1284 @property
1285 def input_sample_rate(self):
1286 if self.decimation:
1287 return self.decimation.input_sample_rate.value
1289 return None
1291 @property
1292 def output_sample_rate(self):
1293 if self.decimation:
1294 return self.decimation.input_sample_rate.value \
1295 / self.decimation.factor
1297 return None
1300class Response(Object):
1301 resource_id = String.T(optional=True, xmlstyle='attribute')
1302 instrument_sensitivity = Sensitivity.T(optional=True,
1303 xmltagname='InstrumentSensitivity')
1304 instrument_polynomial = Polynomial.T(optional=True,
1305 xmltagname='InstrumentPolynomial')
1306 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1308 def check_sample_rates(self, channel):
1310 if self.stage_list:
1311 sample_rate = None
1313 for stage in self.stage_list:
1314 if stage.decimation:
1315 input_sample_rate = \
1316 stage.decimation.input_sample_rate.value
1318 if sample_rate is not None and not same_sample_rate(
1319 sample_rate, input_sample_rate):
1321 logger.warning(
1322 'Response stage %i has unexpected input sample '
1323 'rate: %g Hz (expected: %g Hz)' % (
1324 stage.number,
1325 input_sample_rate,
1326 sample_rate))
1328 sample_rate = input_sample_rate / stage.decimation.factor
1330 if sample_rate is not None and channel.sample_rate \
1331 and not same_sample_rate(
1332 sample_rate, channel.sample_rate.value):
1334 logger.warning(
1335 'Channel sample rate (%g Hz) does not match sample rate '
1336 'deducted from response stages information (%g Hz).' % (
1337 channel.sample_rate.value,
1338 sample_rate))
1340 def check_units(self):
1342 if self.instrument_sensitivity \
1343 and self.instrument_sensitivity.input_units:
1345 units = self.instrument_sensitivity.input_units.name.upper()
1347 if self.stage_list:
1348 for stage in self.stage_list:
1349 if units and stage.input_units \
1350 and stage.input_units.name.upper() != units:
1352 logger.warning(
1353 'Input units of stage %i (%s) do not match %s (%s).'
1354 % (
1355 stage.number,
1356 units,
1357 'output units of stage %i'
1358 if stage.number == 0
1359 else 'sensitivity input units',
1360 units))
1362 if stage.output_units:
1363 units = stage.output_units.name.upper()
1364 else:
1365 units = None
1367 sout_units = self.instrument_sensitivity.output_units
1368 if self.instrument_sensitivity and sout_units:
1369 if units is not None and units != sout_units.name.upper():
1370 logger.warning(
1371 'Output units of stage %i (%s) do not match %s (%s).'
1372 % (
1373 stage.number,
1374 units,
1375 'sensitivity output units',
1376 sout_units.name.upper()))
1378 def _sensitivity_checkpoints(self, responses, context):
1379 delivery = Delivery()
1381 if self.instrument_sensitivity:
1382 sval = self.instrument_sensitivity.value
1383 sfreq = self.instrument_sensitivity.frequency
1384 if sval is None:
1385 delivery.log.append((
1386 'warning',
1387 'No sensitivity value given.',
1388 context))
1390 elif sval is None:
1391 delivery.log.append((
1392 'warning',
1393 'Reported sensitivity value is zero.',
1394 context))
1396 elif sfreq is None:
1397 delivery.log.append((
1398 'warning',
1399 'Sensitivity frequency not given.',
1400 context))
1402 else:
1403 trial = response.MultiplyResponse(responses)
1405 delivery.extend(
1406 check_resp(
1407 trial, sval, sfreq, 0.1,
1408 'Instrument sensitivity value inconsistent with '
1409 'sensitivity computed from complete response',
1410 context))
1412 delivery.payload.append(response.FrequencyResponseCheckpoint(
1413 frequency=sfreq, value=sval))
1415 return delivery
1417 def get_squirrel_response(self, context, **kwargs):
1418 from pyrocko.squirrel.model import Response
1420 if self.stage_list:
1421 delivery = Delivery()
1422 for istage, stage in enumerate(self.stage_list):
1423 delivery.extend(stage.get_squirrel_response_stage(context))
1425 checkpoints = []
1426 if not delivery.errors:
1427 all_responses = []
1428 for sq_stage in delivery.payload:
1429 all_responses.extend(sq_stage.elements)
1431 checkpoints.extend(delivery.extend_without_payload(
1432 self._sensitivity_checkpoints(all_responses, context)))
1434 sq_stages = delivery.expect()
1436 return Response(
1437 stages=sq_stages,
1438 log=delivery.log,
1439 checkpoints=checkpoints,
1440 **kwargs)
1442 elif self.instrument_sensitivity:
1443 raise NoResponseInformation(
1444 'Only instrument sensitivity given (won\'t use it). (%s).'
1445 % context)
1446 else:
1447 raise NoResponseInformation(
1448 'Empty instrument response (%s).'
1449 % context)
1451 def get_pyrocko_response(
1452 self, context, fake_input_units=None, stages=(0, 1)):
1454 delivery = Delivery()
1455 if self.stage_list:
1456 for istage, stage in enumerate(self.stage_list):
1457 delivery.extend(stage.get_pyrocko_response(
1458 context, gain_only=not (
1459 stages is None or stages[0] <= istage < stages[1])))
1461 elif self.instrument_sensitivity:
1462 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1464 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1465 checkpoints = delivery.extend_without_payload(delivery_cp)
1466 if delivery.errors:
1467 return delivery
1469 if fake_input_units is not None:
1470 if not self.instrument_sensitivity or \
1471 self.instrument_sensitivity.input_units is None:
1473 delivery.errors.append((
1474 'NoResponseInformation',
1475 'No input units given, so cannot convert to requested '
1476 'units: %s' % fake_input_units.upper(),
1477 context))
1479 return delivery
1481 input_units = self.instrument_sensitivity.input_units.name.upper()
1483 conresp = None
1484 try:
1485 conresp = conversion[
1486 fake_input_units.upper(), input_units]
1488 except KeyError:
1489 delivery.errors.append((
1490 'NoResponseInformation',
1491 'Cannot convert between units: %s, %s'
1492 % (fake_input_units, input_units),
1493 context))
1495 if conresp is not None:
1496 delivery.payload.append(conresp)
1497 for checkpoint in checkpoints:
1498 checkpoint.value *= num.abs(evaluate1(
1499 conresp, checkpoint.frequency))
1501 delivery.payload = [
1502 response.MultiplyResponse(
1503 delivery.payload,
1504 checkpoints=checkpoints)]
1506 return delivery
1508 @classmethod
1509 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1510 normalization_frequency=1.0):
1511 '''
1512 Convert Pyrocko pole-zero response to StationXML response.
1514 :param presponse: Pyrocko pole-zero response
1515 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1516 :param input_unit: Input unit to be reported in the StationXML
1517 response.
1518 :type input_unit: str
1519 :param output_unit: Output unit to be reported in the StationXML
1520 response.
1521 :type output_unit: str
1522 :param normalization_frequency: Frequency where the normalization
1523 factor for the StationXML response should be computed.
1524 :type normalization_frequency: float
1525 '''
1527 norm_factor = 1.0/float(abs(
1528 evaluate1(presponse, normalization_frequency)
1529 / presponse.constant))
1531 pzs = PolesZeros(
1532 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1533 normalization_factor=norm_factor,
1534 normalization_frequency=Frequency(normalization_frequency),
1535 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1536 imaginary=FloatNoUnit(z.imag))
1537 for z in presponse.zeros],
1538 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1539 imaginary=FloatNoUnit(z.imag))
1540 for z in presponse.poles])
1542 pzs.validate()
1544 stage = ResponseStage(
1545 number=1,
1546 poles_zeros_list=[pzs],
1547 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1549 resp = Response(
1550 instrument_sensitivity=Sensitivity(
1551 value=stage.stage_gain.value,
1552 frequency=normalization_frequency,
1553 input_units=Units(input_unit),
1554 output_units=Units(output_unit)),
1556 stage_list=[stage])
1558 return resp
1561class BaseNode(Object):
1562 '''
1563 A base node type for derivation from: Network, Station and Channel types.
1564 '''
1566 code = String.T(xmlstyle='attribute')
1567 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1568 xmlstyle='attribute')
1569 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1570 xmlstyle='attribute')
1571 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1572 alternate_code = String.T(optional=True, xmlstyle='attribute')
1573 historical_code = String.T(optional=True, xmlstyle='attribute')
1574 description = Unicode.T(optional=True, xmltagname='Description')
1575 comment_list = List.T(Comment.T(xmltagname='Comment'))
1577 def spans(self, *args):
1578 if len(args) == 0:
1579 return True
1580 elif len(args) == 1:
1581 return ((self.start_date is None or
1582 self.start_date <= args[0]) and
1583 (self.end_date is None or
1584 args[0] <= self.end_date))
1586 elif len(args) == 2:
1587 return ((self.start_date is None or
1588 args[1] >= self.start_date) and
1589 (self.end_date is None or
1590 self.end_date >= args[0]))
1593def overlaps(a, b):
1594 return (
1595 a.start_date is None or b.end_date is None
1596 or a.start_date < b.end_date
1597 ) and (
1598 b.start_date is None or a.end_date is None
1599 or b.start_date < a.end_date
1600 )
1603class Channel(BaseNode):
1604 '''
1605 Equivalent to SEED blockette 52 and parent element for the related the
1606 response blockettes.
1607 '''
1609 location_code = String.T(xmlstyle='attribute')
1610 external_reference_list = List.T(
1611 ExternalReference.T(xmltagname='ExternalReference'))
1612 latitude = Latitude.T(xmltagname='Latitude')
1613 longitude = Longitude.T(xmltagname='Longitude')
1614 elevation = Distance.T(xmltagname='Elevation')
1615 depth = Distance.T(xmltagname='Depth')
1616 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1617 dip = Dip.T(optional=True, xmltagname='Dip')
1618 type_list = List.T(Type.T(xmltagname='Type'))
1619 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1620 sample_rate_ratio = SampleRateRatio.T(optional=True,
1621 xmltagname='SampleRateRatio')
1622 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1623 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1624 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1625 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1626 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1627 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1628 equipment = Equipment.T(optional=True, xmltagname='Equipment')
1629 response = Response.T(optional=True, xmltagname='Response')
1631 @property
1632 def position_values(self):
1633 lat = self.latitude.value
1634 lon = self.longitude.value
1635 elevation = value_or_none(self.elevation)
1636 depth = value_or_none(self.depth)
1637 return lat, lon, elevation, depth
1640class Station(BaseNode):
1641 '''
1642 This type represents a Station epoch. It is common to only have a single
1643 station epoch with the station's creation and termination dates as the
1644 epoch start and end dates.
1645 '''
1647 latitude = Latitude.T(xmltagname='Latitude')
1648 longitude = Longitude.T(xmltagname='Longitude')
1649 elevation = Distance.T(xmltagname='Elevation')
1650 site = Site.T(optional=True, xmltagname='Site')
1651 vault = Unicode.T(optional=True, xmltagname='Vault')
1652 geology = Unicode.T(optional=True, xmltagname='Geology')
1653 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1654 operator_list = List.T(Operator.T(xmltagname='Operator'))
1655 creation_date = DummyAwareOptionalTimestamp.T(
1656 optional=True, xmltagname='CreationDate')
1657 termination_date = DummyAwareOptionalTimestamp.T(
1658 optional=True, xmltagname='TerminationDate')
1659 total_number_channels = Counter.T(
1660 optional=True, xmltagname='TotalNumberChannels')
1661 selected_number_channels = Counter.T(
1662 optional=True, xmltagname='SelectedNumberChannels')
1663 external_reference_list = List.T(
1664 ExternalReference.T(xmltagname='ExternalReference'))
1665 channel_list = List.T(Channel.T(xmltagname='Channel'))
1667 @property
1668 def position_values(self):
1669 lat = self.latitude.value
1670 lon = self.longitude.value
1671 elevation = value_or_none(self.elevation)
1672 return lat, lon, elevation
1675class Network(BaseNode):
1676 '''
1677 This type represents the Network layer, all station metadata is contained
1678 within this element. The official name of the network or other descriptive
1679 information can be included in the Description element. The Network can
1680 contain 0 or more Stations.
1681 '''
1683 total_number_stations = Counter.T(optional=True,
1684 xmltagname='TotalNumberStations')
1685 selected_number_stations = Counter.T(optional=True,
1686 xmltagname='SelectedNumberStations')
1687 station_list = List.T(Station.T(xmltagname='Station'))
1689 @property
1690 def station_code_list(self):
1691 return sorted(set(s.code for s in self.station_list))
1693 @property
1694 def sl_code_list(self):
1695 sls = set()
1696 for station in self.station_list:
1697 for channel in station.channel_list:
1698 sls.add((station.code, channel.location_code))
1700 return sorted(sls)
1702 def summary(self, width=80, indent=4):
1703 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1704 lines = ['%s (%i):' % (self.code, len(sls))]
1705 if sls:
1706 ssls = ['.'.join(x for x in c if x) for c in sls]
1707 w = max(len(x) for x in ssls)
1708 n = (width - indent) / (w+1)
1709 while ssls:
1710 lines.append(
1711 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1713 ssls[:n] = []
1715 return '\n'.join(lines)
1718def value_or_none(x):
1719 if x is not None:
1720 return x.value
1721 else:
1722 return None
1725def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1727 pos = lat, lon, elevation, depth = \
1728 channels[0].position_values
1730 if not all(pos == x.position_values for x in channels):
1731 info = '\n'.join(
1732 ' %s: %s' % (x.code, x.position_values) for
1733 x in channels)
1735 mess = 'encountered inconsistencies in channel ' \
1736 'lat/lon/elevation/depth ' \
1737 'for %s.%s.%s: \n%s' % (nsl + (info,))
1739 if inconsistencies == 'raise':
1740 raise InconsistentChannelLocations(mess)
1742 elif inconsistencies == 'warn':
1743 logger.warning(mess)
1744 logger.warning(' -> using mean values')
1746 apos = num.array([x.position_values for x in channels], dtype=float)
1747 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1748 / num.sum(num.isfinite(apos), axis=0)
1750 groups = {}
1751 for channel in channels:
1752 if channel.code not in groups:
1753 groups[channel.code] = []
1755 groups[channel.code].append(channel)
1757 pchannels = []
1758 for code in sorted(groups.keys()):
1759 data = [
1760 (channel.code, value_or_none(channel.azimuth),
1761 value_or_none(channel.dip))
1762 for channel in groups[code]]
1764 azimuth, dip = util.consistency_merge(
1765 data,
1766 message='channel orientation values differ:',
1767 error=inconsistencies)
1769 pchannels.append(
1770 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1772 return pyrocko.model.Station(
1773 *nsl,
1774 lat=mlat,
1775 lon=mlon,
1776 elevation=mele,
1777 depth=mdep,
1778 channels=pchannels)
1781class FDSNStationXML(Object):
1782 '''
1783 Top-level type for Station XML. Required field are Source (network ID of
1784 the institution sending the message) and one or more Network containers or
1785 one or more Station containers.
1786 '''
1788 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1789 source = String.T(xmltagname='Source')
1790 sender = String.T(optional=True, xmltagname='Sender')
1791 module = String.T(optional=True, xmltagname='Module')
1792 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1793 created = Timestamp.T(optional=True, xmltagname='Created')
1794 network_list = List.T(Network.T(xmltagname='Network'))
1796 xmltagname = 'FDSNStationXML'
1797 guessable_xmlns = [guts_xmlns]
1799 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1800 time=None, timespan=None,
1801 inconsistencies='warn'):
1803 assert inconsistencies in ('raise', 'warn')
1805 if nslcs is not None:
1806 nslcs = set(nslcs)
1808 if nsls is not None:
1809 nsls = set(nsls)
1811 tt = ()
1812 if time is not None:
1813 tt = (time,)
1814 elif timespan is not None:
1815 tt = timespan
1817 pstations = []
1818 for network in self.network_list:
1819 if not network.spans(*tt):
1820 continue
1822 for station in network.station_list:
1823 if not station.spans(*tt):
1824 continue
1826 if station.channel_list:
1827 loc_to_channels = {}
1828 for channel in station.channel_list:
1829 if not channel.spans(*tt):
1830 continue
1832 loc = channel.location_code.strip()
1833 if loc not in loc_to_channels:
1834 loc_to_channels[loc] = []
1836 loc_to_channels[loc].append(channel)
1838 for loc in sorted(loc_to_channels.keys()):
1839 channels = loc_to_channels[loc]
1840 if nslcs is not None:
1841 channels = [channel for channel in channels
1842 if (network.code, station.code, loc,
1843 channel.code) in nslcs]
1845 if not channels:
1846 continue
1848 nsl = network.code, station.code, loc
1849 if nsls is not None and nsl not in nsls:
1850 continue
1852 pstations.append(
1853 pyrocko_station_from_channels(
1854 nsl,
1855 channels,
1856 inconsistencies=inconsistencies))
1857 else:
1858 pstations.append(pyrocko.model.Station(
1859 network.code, station.code, '*',
1860 lat=station.latitude.value,
1861 lon=station.longitude.value,
1862 elevation=value_or_none(station.elevation),
1863 name=station.description or ''))
1865 return pstations
1867 @classmethod
1868 def from_pyrocko_stations(
1869 cls, pyrocko_stations, add_flat_responses_from=None):
1871 '''
1872 Generate :py:class:`FDSNStationXML` from list of
1873 :py:class;`pyrocko.model.Station` instances.
1875 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1876 instances.
1877 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1878 '''
1879 from collections import defaultdict
1880 network_dict = defaultdict(list)
1882 if add_flat_responses_from:
1883 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1884 extra = dict(
1885 response=Response(
1886 instrument_sensitivity=Sensitivity(
1887 value=1.0,
1888 frequency=1.0,
1889 input_units=Units(name=add_flat_responses_from))))
1890 else:
1891 extra = {}
1893 have_offsets = set()
1894 for s in pyrocko_stations:
1896 if s.north_shift != 0.0 or s.east_shift != 0.0:
1897 have_offsets.add(s.nsl())
1899 network, station, location = s.nsl()
1900 channel_list = []
1901 for c in s.channels:
1902 channel_list.append(
1903 Channel(
1904 location_code=location,
1905 code=c.name,
1906 latitude=Latitude(value=s.effective_lat),
1907 longitude=Longitude(value=s.effective_lon),
1908 elevation=Distance(value=s.elevation),
1909 depth=Distance(value=s.depth),
1910 azimuth=Azimuth(value=c.azimuth),
1911 dip=Dip(value=c.dip),
1912 **extra
1913 )
1914 )
1916 network_dict[network].append(
1917 Station(
1918 code=station,
1919 latitude=Latitude(value=s.effective_lat),
1920 longitude=Longitude(value=s.effective_lon),
1921 elevation=Distance(value=s.elevation),
1922 channel_list=channel_list)
1923 )
1925 if have_offsets:
1926 logger.warning(
1927 'StationXML does not support Cartesian offsets in '
1928 'coordinates. Storing effective lat/lon for stations: %s' %
1929 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1931 timestamp = util.to_time_float(time.time())
1932 network_list = []
1933 for k, station_list in network_dict.items():
1935 network_list.append(
1936 Network(
1937 code=k, station_list=station_list,
1938 total_number_stations=len(station_list)))
1940 sxml = FDSNStationXML(
1941 source='from pyrocko stations list', created=timestamp,
1942 network_list=network_list)
1944 sxml.validate()
1945 return sxml
1947 def iter_network_stations(
1948 self, net=None, sta=None, time=None, timespan=None):
1950 tt = ()
1951 if time is not None:
1952 tt = (time,)
1953 elif timespan is not None:
1954 tt = timespan
1956 for network in self.network_list:
1957 if not network.spans(*tt) or (
1958 net is not None and network.code != net):
1959 continue
1961 for station in network.station_list:
1962 if not station.spans(*tt) or (
1963 sta is not None and station.code != sta):
1964 continue
1966 yield (network, station)
1968 def iter_network_station_channels(
1969 self, net=None, sta=None, loc=None, cha=None,
1970 time=None, timespan=None):
1972 if loc is not None:
1973 loc = loc.strip()
1975 tt = ()
1976 if time is not None:
1977 tt = (time,)
1978 elif timespan is not None:
1979 tt = timespan
1981 for network in self.network_list:
1982 if not network.spans(*tt) or (
1983 net is not None and network.code != net):
1984 continue
1986 for station in network.station_list:
1987 if not station.spans(*tt) or (
1988 sta is not None and station.code != sta):
1989 continue
1991 if station.channel_list:
1992 for channel in station.channel_list:
1993 if (not channel.spans(*tt) or
1994 (cha is not None and channel.code != cha) or
1995 (loc is not None and
1996 channel.location_code.strip() != loc)):
1997 continue
1999 yield (network, station, channel)
2001 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2002 time=None, timespan=None):
2004 groups = {}
2005 for network, station, channel in self.iter_network_station_channels(
2006 net, sta, loc, cha, time=time, timespan=timespan):
2008 net = network.code
2009 sta = station.code
2010 cha = channel.code
2011 loc = channel.location_code.strip()
2012 if len(cha) == 3:
2013 bic = cha[:2] # band and intrument code according to SEED
2014 elif len(cha) == 1:
2015 bic = ''
2016 else:
2017 bic = cha
2019 if channel.response and \
2020 channel.response.instrument_sensitivity and \
2021 channel.response.instrument_sensitivity.input_units:
2023 unit = channel.response.instrument_sensitivity\
2024 .input_units.name.upper()
2025 else:
2026 unit = None
2028 bic = (bic, unit)
2030 k = net, sta, loc
2031 if k not in groups:
2032 groups[k] = {}
2034 if bic not in groups[k]:
2035 groups[k][bic] = []
2037 groups[k][bic].append(channel)
2039 for nsl, bic_to_channels in groups.items():
2040 bad_bics = []
2041 for bic, channels in bic_to_channels.items():
2042 sample_rates = []
2043 for channel in channels:
2044 sample_rates.append(channel.sample_rate.value)
2046 if not same(sample_rates):
2047 scs = ','.join(channel.code for channel in channels)
2048 srs = ', '.join('%e' % x for x in sample_rates)
2049 err = 'ignoring channels with inconsistent sampling ' + \
2050 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2052 logger.warning(err)
2053 bad_bics.append(bic)
2055 for bic in bad_bics:
2056 del bic_to_channels[bic]
2058 return groups
2060 def choose_channels(
2061 self,
2062 target_sample_rate=None,
2063 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2064 priority_units=['M/S', 'M/S**2'],
2065 priority_instrument_code=['H', 'L'],
2066 time=None,
2067 timespan=None):
2069 nslcs = {}
2070 for nsl, bic_to_channels in self.get_channel_groups(
2071 time=time, timespan=timespan).items():
2073 useful_bics = []
2074 for bic, channels in bic_to_channels.items():
2075 rate = channels[0].sample_rate.value
2077 if target_sample_rate is not None and \
2078 rate < target_sample_rate*0.99999:
2079 continue
2081 if len(bic[0]) == 2:
2082 if bic[0][0] not in priority_band_code:
2083 continue
2085 if bic[0][1] not in priority_instrument_code:
2086 continue
2088 unit = bic[1]
2090 prio_unit = len(priority_units)
2091 try:
2092 prio_unit = priority_units.index(unit)
2093 except ValueError:
2094 pass
2096 prio_inst = len(priority_instrument_code)
2097 prio_band = len(priority_band_code)
2098 if len(channels[0].code) == 3:
2099 try:
2100 prio_inst = priority_instrument_code.index(
2101 channels[0].code[1])
2102 except ValueError:
2103 pass
2105 try:
2106 prio_band = priority_band_code.index(
2107 channels[0].code[0])
2108 except ValueError:
2109 pass
2111 if target_sample_rate is None:
2112 rate = -rate
2114 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2115 prio_inst, bic))
2117 useful_bics.sort()
2119 for _, _, rate, _, _, bic in useful_bics:
2120 channels = sorted(
2121 bic_to_channels[bic],
2122 key=lambda channel: channel.code)
2124 if channels:
2125 for channel in channels:
2126 nslcs[nsl + (channel.code,)] = channel
2128 break
2130 return nslcs
2132 def get_pyrocko_response(
2133 self, nslc,
2134 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2136 net, sta, loc, cha = nslc
2137 resps = []
2138 for _, _, channel in self.iter_network_station_channels(
2139 net, sta, loc, cha, time=time, timespan=timespan):
2140 resp = channel.response
2141 if resp:
2142 resp.check_sample_rates(channel)
2143 resp.check_units()
2144 resps.append(resp.get_pyrocko_response(
2145 '.'.join(nslc),
2146 fake_input_units=fake_input_units,
2147 stages=stages).expect_one())
2149 if not resps:
2150 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2151 elif len(resps) > 1:
2152 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2154 return resps[0]
2156 @property
2157 def n_code_list(self):
2158 return sorted(set(x.code for x in self.network_list))
2160 @property
2161 def ns_code_list(self):
2162 nss = set()
2163 for network in self.network_list:
2164 for station in network.station_list:
2165 nss.add((network.code, station.code))
2167 return sorted(nss)
2169 @property
2170 def nsl_code_list(self):
2171 nsls = set()
2172 for network in self.network_list:
2173 for station in network.station_list:
2174 for channel in station.channel_list:
2175 nsls.add(
2176 (network.code, station.code, channel.location_code))
2178 return sorted(nsls)
2180 @property
2181 def nslc_code_list(self):
2182 nslcs = set()
2183 for network in self.network_list:
2184 for station in network.station_list:
2185 for channel in station.channel_list:
2186 nslcs.add(
2187 (network.code, station.code, channel.location_code,
2188 channel.code))
2190 return sorted(nslcs)
2192 def summary(self):
2193 lst = [
2194 'number of n codes: %i' % len(self.n_code_list),
2195 'number of ns codes: %i' % len(self.ns_code_list),
2196 'number of nsl codes: %i' % len(self.nsl_code_list),
2197 'number of nslc codes: %i' % len(self.nslc_code_list)
2198 ]
2199 return '\n'.join(lst)
2201 def summary_stages(self):
2202 data = []
2203 for network, station, channel in self.iter_network_station_channels():
2204 nslc = (network.code, station.code, channel.location_code,
2205 channel.code)
2207 stages = []
2208 in_units = '?'
2209 out_units = '?'
2210 if channel.response:
2211 sens = channel.response.instrument_sensitivity
2212 if sens:
2213 in_units = sens.input_units.name.upper()
2214 out_units = sens.output_units.name.upper()
2216 for stage in channel.response.stage_list:
2217 stages.append(stage.summary())
2219 data.append(
2220 (nslc, tts(channel.start_date), tts(channel.end_date),
2221 in_units, out_units, stages))
2223 data.sort()
2225 lst = []
2226 for nslc, stmin, stmax, in_units, out_units, stages in data:
2227 lst.append(' %s: %s - %s, %s -> %s' % (
2228 '.'.join(nslc), stmin, stmax, in_units, out_units))
2229 for stage in stages:
2230 lst.append(' %s' % stage)
2232 return '\n'.join(lst)
2234 def _check_overlaps(self):
2235 by_nslc = {}
2236 for network in self.network_list:
2237 for station in network.station_list:
2238 for channel in station.channel_list:
2239 nslc = (network.code, station.code, channel.location_code,
2240 channel.code)
2241 if nslc not in by_nslc:
2242 by_nslc[nslc] = []
2244 by_nslc[nslc].append(channel)
2246 errors = []
2247 for nslc, channels in by_nslc.items():
2248 for ia, a in enumerate(channels):
2249 for b in channels[ia+1:]:
2250 if overlaps(a, b):
2251 errors.append(
2252 'Channel epochs overlap for %s:\n'
2253 ' %s - %s\n %s - %s' % (
2254 '.'.join(nslc),
2255 tts(a.start_date), tts(a.end_date),
2256 tts(b.start_date), tts(b.end_date)))
2257 return errors
2259 def check(self):
2260 errors = []
2261 for meth in [self._check_overlaps]:
2262 errors.extend(meth())
2264 if errors:
2265 raise Inconsistencies(
2266 'Inconsistencies found in StationXML:\n '
2267 + '\n '.join(errors))
2270def load_channel_table(stream):
2272 networks = {}
2273 stations = {}
2275 for line in stream:
2276 line = str(line.decode('ascii'))
2277 if line.startswith('#'):
2278 continue
2280 t = line.rstrip().split('|')
2282 if len(t) != 17:
2283 logger.warning('Invalid channel record: %s' % line)
2284 continue
2286 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2287 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2289 try:
2290 scale = float(scale)
2291 except ValueError:
2292 scale = None
2294 try:
2295 scale_freq = float(scale_freq)
2296 except ValueError:
2297 scale_freq = None
2299 try:
2300 depth = float(dep)
2301 except ValueError:
2302 depth = 0.0
2304 try:
2305 azi = float(azi)
2306 dip = float(dip)
2307 except ValueError:
2308 azi = None
2309 dip = None
2311 try:
2312 if net not in networks:
2313 network = Network(code=net)
2314 else:
2315 network = networks[net]
2317 if (net, sta) not in stations:
2318 station = Station(
2319 code=sta, latitude=lat, longitude=lon, elevation=ele)
2321 station.regularize()
2322 else:
2323 station = stations[net, sta]
2325 if scale:
2326 resp = Response(
2327 instrument_sensitivity=Sensitivity(
2328 value=scale,
2329 frequency=scale_freq,
2330 input_units=scale_units))
2331 else:
2332 resp = None
2334 channel = Channel(
2335 code=cha,
2336 location_code=loc.strip(),
2337 latitude=lat,
2338 longitude=lon,
2339 elevation=ele,
2340 depth=depth,
2341 azimuth=azi,
2342 dip=dip,
2343 sensor=Equipment(description=sens),
2344 response=resp,
2345 sample_rate=sample_rate,
2346 start_date=start_date,
2347 end_date=end_date or None)
2349 channel.regularize()
2351 except ValidationError:
2352 raise InvalidRecord(line)
2354 if net not in networks:
2355 networks[net] = network
2357 if (net, sta) not in stations:
2358 stations[net, sta] = station
2359 network.station_list.append(station)
2361 station.channel_list.append(channel)
2363 return FDSNStationXML(
2364 source='created from table input',
2365 created=time.time(),
2366 network_list=sorted(networks.values(), key=lambda x: x.code))
2369def primitive_merge(sxs):
2370 networks = []
2371 for sx in sxs:
2372 networks.extend(sx.network_list)
2374 return FDSNStationXML(
2375 source='merged from different sources',
2376 created=time.time(),
2377 network_list=copy.deepcopy(
2378 sorted(networks, key=lambda x: x.code)))
2381def split_channels(sx):
2382 for nslc in sx.nslc_code_list:
2383 network_list = sx.network_list
2384 network_list_filtered = [
2385 network for network in network_list
2386 if network.code == nslc[0]]
2388 for network in network_list_filtered:
2389 sx.network_list = [network]
2390 station_list = network.station_list
2391 station_list_filtered = [
2392 station for station in station_list
2393 if station.code == nslc[1]]
2395 for station in station_list_filtered:
2396 network.station_list = [station]
2397 channel_list = station.channel_list
2398 station.channel_list = [
2399 channel for channel in channel_list
2400 if (channel.location_code, channel.code) == nslc[2:4]]
2402 yield nslc, copy.deepcopy(sx)
2404 station.channel_list = channel_list
2406 network.station_list = station_list
2408 sx.network_list = network_list
2411if __name__ == '__main__':
2412 from optparse import OptionParser
2414 util.setup_logging('pyrocko.io.stationxml', 'warning')
2416 usage = \
2417 'python -m pyrocko.io.stationxml check|stats|stages ' \
2418 '<filename> [options]'
2420 description = '''Torture StationXML file.'''
2422 parser = OptionParser(
2423 usage=usage,
2424 description=description,
2425 formatter=util.BetterHelpFormatter())
2427 (options, args) = parser.parse_args(sys.argv[1:])
2429 if len(args) != 2:
2430 parser.print_help()
2431 sys.exit(1)
2433 action, path = args
2435 sx = load_xml(filename=path)
2436 if action == 'check':
2437 try:
2438 sx.check()
2439 except Inconsistencies as e:
2440 logger.error(e)
2441 sys.exit(1)
2443 elif action == 'stats':
2444 print(sx.summary())
2446 elif action == 'stages':
2447 print(sx.summary_stages())
2449 else:
2450 parser.print_help()
2451 sys.exit('unknown action: %s' % action)