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 dummy_for = (hpfloat, float)
252 dummy_for_description = 'time_float'
254 class __T(TBase):
256 def regularize_extra(self, val):
257 time_float = get_time_float()
259 if isinstance(val, datetime.datetime):
260 tt = val.utctimetuple()
261 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
263 elif isinstance(val, datetime.date):
264 tt = val.timetuple()
265 val = time_float(calendar.timegm(tt))
267 elif isinstance(val, (str, newstr)):
268 val = val.strip()
270 tz_offset = 0
272 m = re_tz.search(val)
273 if m:
274 sh = m.group(2)
275 sm = m.group(4)
276 tz_offset = (int(sh)*3600 if sh else 0) \
277 + (int(sm)*60 if sm else 0)
279 val = re_tz.sub('', val)
281 if len(val) > 10 and val[10] == 'T':
282 val = val.replace('T', ' ', 1)
284 try:
285 val = util.str_to_time(val) - tz_offset
287 except util.TimeStrError:
288 year = int(val[:4])
289 if sys.maxsize > 2**32: # if we're on 64bit
290 if year > this_year + 100:
291 return None # StationXML contained a dummy date
293 if year < 1903: # for macOS, 1900-01-01 dummy dates
294 return None
296 else: # 32bit end of time is in 2038
297 if this_year < 2037 and year > 2037 or year < 1903:
298 return None # StationXML contained a dummy date
300 raise
302 elif isinstance(val, (int, float)):
303 val = time_float(val)
305 else:
306 raise ValidationError(
307 '%s: cannot convert "%s" to type %s' % (
308 self.xname(), val, time_float))
310 return val
312 def to_save(self, val):
313 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
314 .rstrip('0').rstrip('.')
316 def to_save_xml(self, val):
317 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
318 .rstrip('0').rstrip('.') + 'Z'
321class Nominal(StringChoice):
322 choices = [
323 'NOMINAL',
324 'CALCULATED']
327class Email(UnicodePattern):
328 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
331class RestrictedStatus(StringChoice):
332 choices = [
333 'open',
334 'closed',
335 'partial']
338class Type(StringChoice):
339 choices = [
340 'TRIGGERED',
341 'CONTINUOUS',
342 'HEALTH',
343 'GEOPHYSICAL',
344 'WEATHER',
345 'FLAG',
346 'SYNTHESIZED',
347 'INPUT',
348 'EXPERIMENTAL',
349 'MAINTENANCE',
350 'BEAM']
352 class __T(StringChoice.T):
353 def validate_extra(self, val):
354 if val not in self.choices:
355 logger.warning(
356 'channel type: "%s" is not a valid choice out of %s' %
357 (val, repr(self.choices)))
360class PzTransferFunction(StringChoice):
361 choices = [
362 'LAPLACE (RADIANS/SECOND)',
363 'LAPLACE (HERTZ)',
364 'DIGITAL (Z-TRANSFORM)']
367class Symmetry(StringChoice):
368 choices = [
369 'NONE',
370 'EVEN',
371 'ODD']
374class CfTransferFunction(StringChoice):
376 class __T(StringChoice.T):
377 def validate(self, val, regularize=False, depth=-1):
378 if regularize:
379 try:
380 val = str(val)
381 except ValueError:
382 raise ValidationError(
383 '%s: cannot convert to string %s' % (self.xname,
384 repr(val)))
386 val = self.dummy_cls.replacements.get(val, val)
388 self.validate_extra(val)
389 return val
391 choices = [
392 'ANALOG (RADIANS/SECOND)',
393 'ANALOG (HERTZ)',
394 'DIGITAL']
396 replacements = {
397 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
398 'ANALOG (HZ)': 'ANALOG (HERTZ)',
399 }
402class Approximation(StringChoice):
403 choices = [
404 'MACLAURIN']
407class PhoneNumber(StringPattern):
408 pattern = '[0-9]+-[0-9]+'
411class Site(Object):
412 '''
413 Description of a site location using name and optional geopolitical
414 boundaries (country, city, etc.).
415 '''
417 name = Unicode.T(default='', xmltagname='Name')
418 description = Unicode.T(optional=True, xmltagname='Description')
419 town = Unicode.T(optional=True, xmltagname='Town')
420 county = Unicode.T(optional=True, xmltagname='County')
421 region = Unicode.T(optional=True, xmltagname='Region')
422 country = Unicode.T(optional=True, xmltagname='Country')
425class ExternalReference(Object):
426 '''
427 This type contains a URI and description for external data that users may
428 want to reference in StationXML.
429 '''
431 uri = String.T(xmltagname='URI')
432 description = Unicode.T(xmltagname='Description')
435class Units(Object):
436 '''
437 A type to document units. Corresponds to SEED blockette 34.
438 '''
440 def __init__(self, name=None, **kwargs):
441 Object.__init__(self, name=name, **kwargs)
443 name = String.T(xmltagname='Name')
444 description = Unicode.T(optional=True, xmltagname='Description')
447class Counter(Int):
448 pass
451class SampleRateRatio(Object):
452 '''
453 Sample rate expressed as number of samples in a number of seconds.
454 '''
456 number_samples = Int.T(xmltagname='NumberSamples')
457 number_seconds = Int.T(xmltagname='NumberSeconds')
460class Gain(Object):
461 '''
462 Complex type for sensitivity and frequency ranges. This complex type can be
463 used to represent both overall sensitivities and individual stage gains.
464 The FrequencyRangeGroup is an optional construct that defines a pass band
465 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
466 is valid within the number of decibels specified in FrequencyDBVariation.
467 '''
469 def __init__(self, value=None, **kwargs):
470 Object.__init__(self, value=value, **kwargs)
472 value = Float.T(optional=True, xmltagname='Value')
473 frequency = Float.T(optional=True, xmltagname='Frequency')
475 def summary(self):
476 return 'gain(%g @ %g)' % (self.value, self.frequency)
479class NumeratorCoefficient(Object):
480 i = Int.T(optional=True, xmlstyle='attribute')
481 value = Float.T(xmlstyle='content')
484class FloatNoUnit(Object):
485 def __init__(self, value=None, **kwargs):
486 Object.__init__(self, value=value, **kwargs)
488 plus_error = Float.T(optional=True, xmlstyle='attribute')
489 minus_error = Float.T(optional=True, xmlstyle='attribute')
490 value = Float.T(xmlstyle='content')
493class FloatWithUnit(FloatNoUnit):
494 unit = String.T(optional=True, xmlstyle='attribute')
497class Equipment(Object):
498 resource_id = String.T(optional=True, xmlstyle='attribute')
499 type = String.T(optional=True, xmltagname='Type')
500 description = Unicode.T(optional=True, xmltagname='Description')
501 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
502 vendor = Unicode.T(optional=True, xmltagname='Vendor')
503 model = Unicode.T(optional=True, xmltagname='Model')
504 serial_number = String.T(optional=True, xmltagname='SerialNumber')
505 installation_date = DummyAwareOptionalTimestamp.T(
506 optional=True,
507 xmltagname='InstallationDate')
508 removal_date = DummyAwareOptionalTimestamp.T(
509 optional=True,
510 xmltagname='RemovalDate')
511 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
514class PhoneNumber(Object):
515 description = Unicode.T(optional=True, xmlstyle='attribute')
516 country_code = Int.T(optional=True, xmltagname='CountryCode')
517 area_code = Int.T(xmltagname='AreaCode')
518 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
521class BaseFilter(Object):
522 '''
523 The BaseFilter is derived by all filters.
524 '''
526 resource_id = String.T(optional=True, xmlstyle='attribute')
527 name = String.T(optional=True, xmlstyle='attribute')
528 description = Unicode.T(optional=True, xmltagname='Description')
529 input_units = Units.T(optional=True, xmltagname='InputUnits')
530 output_units = Units.T(optional=True, xmltagname='OutputUnits')
533class Sensitivity(Gain):
534 '''
535 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
536 construct that defines a pass band in Hertz (FrequencyStart and
537 FrequencyEnd) in which the SensitivityValue is valid within the number of
538 decibels specified in FrequencyDBVariation.
539 '''
541 input_units = Units.T(optional=True, xmltagname='InputUnits')
542 output_units = Units.T(optional=True, xmltagname='OutputUnits')
543 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
544 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
545 frequency_db_variation = Float.T(optional=True,
546 xmltagname='FrequencyDBVariation')
548 def get_pyrocko_response(self):
549 return Delivery(
550 [response.PoleZeroResponse(constant=self.value)])
553class Coefficient(FloatNoUnit):
554 number = Counter.T(optional=True, xmlstyle='attribute')
557class PoleZero(Object):
558 '''
559 Complex numbers used as poles or zeros in channel response.
560 '''
562 number = Int.T(optional=True, xmlstyle='attribute')
563 real = FloatNoUnit.T(xmltagname='Real')
564 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
566 def value(self):
567 return self.real.value + 1J * self.imaginary.value
570class ClockDrift(FloatWithUnit):
571 unit = String.T(default='SECONDS/SAMPLE', optional=True,
572 xmlstyle='attribute') # fixed
575class Second(FloatWithUnit):
576 '''
577 A time value in seconds.
578 '''
580 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
581 # fixed unit
584class Voltage(FloatWithUnit):
585 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
586 # fixed unit
589class Angle(FloatWithUnit):
590 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
591 # fixed unit
594class Azimuth(FloatWithUnit):
595 '''
596 Instrument azimuth, degrees clockwise from North.
597 '''
599 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
600 # fixed unit
603class Dip(FloatWithUnit):
604 '''
605 Instrument dip in degrees down from horizontal. Together azimuth and dip
606 describe the direction of the sensitive axis of the instrument.
607 '''
609 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
610 # fixed unit
613class Distance(FloatWithUnit):
614 '''
615 Extension of FloatWithUnit for distances, elevations, and depths.
616 '''
618 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
619 # NOT fixed unit!
622class Frequency(FloatWithUnit):
623 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
624 # fixed unit
627class SampleRate(FloatWithUnit):
628 '''
629 Sample rate in samples per second.
630 '''
632 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
633 # fixed unit
636class Person(Object):
637 '''
638 Representation of a person's contact information. A person can belong to
639 multiple agencies and have multiple email addresses and phone numbers.
640 '''
642 name_list = List.T(Unicode.T(xmltagname='Name'))
643 agency_list = List.T(Unicode.T(xmltagname='Agency'))
644 email_list = List.T(Email.T(xmltagname='Email'))
645 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
648class FIR(BaseFilter):
649 '''
650 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
651 also commonly documented using the Coefficients element.
652 '''
654 symmetry = Symmetry.T(xmltagname='Symmetry')
655 numerator_coefficient_list = List.T(
656 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
658 def summary(self):
659 return 'fir(%i%s)' % (
660 self.get_ncoefficiencs(),
661 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
663 def get_effective_coefficients(self):
664 b = num.array(
665 [v.value for v in self.numerator_coefficient_list],
666 dtype=float)
668 if self.symmetry == 'ODD':
669 b = num.concatenate((b, b[-2::-1]))
670 elif self.symmetry == 'EVEN':
671 b = num.concatenate((b, b[::-1]))
673 return b
675 def get_effective_symmetry(self):
676 if self.symmetry == 'NONE':
677 b = self.get_effective_coefficients()
678 if num.all(b - b[::-1] == 0):
679 return ['EVEN', 'ODD'][b.size % 2]
681 return self.symmetry
683 def get_ncoefficiencs(self):
684 nf = len(self.numerator_coefficient_list)
685 if self.symmetry == 'ODD':
686 nc = nf * 2 + 1
687 elif self.symmetry == 'EVEN':
688 nc = nf * 2
689 else:
690 nc = nf
692 return nc
694 def estimate_delay(self, deltat):
695 nc = self.get_ncoefficiencs()
696 if nc > 0:
697 return deltat * (nc - 1) / 2.0
698 else:
699 return 0.0
701 def get_pyrocko_response(
702 self, context, deltat, delay_responses, normalization_frequency):
704 context += self.summary()
706 if not self.numerator_coefficient_list:
707 return Delivery([])
709 b = self.get_effective_coefficients()
711 log = []
712 drop_phase = self.get_effective_symmetry() != 'NONE'
714 if not deltat:
715 log.append((
716 'error',
717 'Digital response requires knowledge about sampling '
718 'interval. Response will be unusable.',
719 context))
721 resp = response.DigitalFilterResponse(
722 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
724 if normalization_frequency is not None and deltat is not None:
725 normalization_frequency = 0.0
726 normalization = num.abs(evaluate1(resp, normalization_frequency))
728 if num.abs(normalization - 1.0) > 1e-2:
729 log.append((
730 'warning',
731 'FIR filter coefficients are not normalized. Normalizing '
732 'them. Factor: %g' % normalization, context))
734 resp = response.DigitalFilterResponse(
735 (b/normalization).tolist(), [1.0], deltat,
736 drop_phase=drop_phase)
738 resps = [resp]
740 if not drop_phase:
741 resps.extend(delay_responses)
743 return Delivery(resps, log=log)
746class Coefficients(BaseFilter):
747 '''
748 Response: coefficients for FIR filter. Laplace transforms or IIR filters
749 can be expressed using type as well but the PolesAndZeros should be used
750 instead. Corresponds to SEED blockette 54.
751 '''
753 cf_transfer_function_type = CfTransferFunction.T(
754 xmltagname='CfTransferFunctionType')
755 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
756 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
758 def summary(self):
759 return 'coef_%s(%i,%i%s)' % (
760 'ABC?'[
761 CfTransferFunction.choices.index(
762 self.cf_transfer_function_type)],
763 len(self.numerator_list),
764 len(self.denominator_list),
765 ',sym' if self.is_symmetric_fir else '')
767 def estimate_delay(self, deltat):
768 nc = len(self.numerator_list)
769 if nc > 0:
770 return deltat * (len(self.numerator_list) - 1) / 2.0
771 else:
772 return 0.0
774 def is_symmetric_fir(self):
775 if len(self.denominator_list) != 0:
776 return False
777 b = [v.value for v in self.numerator_list]
778 return b == b[::-1]
780 def get_pyrocko_response(
781 self, context, deltat, delay_responses, normalization_frequency):
783 context += self.summary()
785 factor = 1.0
786 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
787 factor = 2.0*math.pi
789 if not self.numerator_list and not self.denominator_list:
790 return Delivery(payload=[])
792 b = num.array(
793 [v.value*factor for v in self.numerator_list], dtype=float)
795 a = num.array(
796 [1.0] + [v.value*factor for v in self.denominator_list],
797 dtype=float)
799 log = []
800 resps = []
801 if self.cf_transfer_function_type in [
802 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
803 resps.append(response.AnalogFilterResponse(b, a))
805 elif self.cf_transfer_function_type == 'DIGITAL':
806 if not deltat:
807 log.append((
808 'error',
809 'Digital response requires knowledge about sampling '
810 'interval. Response will be unusable.',
811 context))
813 drop_phase = self.is_symmetric_fir()
814 resp = response.DigitalFilterResponse(
815 b, a, deltat or 0.0, drop_phase=drop_phase)
817 if normalization_frequency is not None and deltat is not None:
818 normalization = num.abs(evaluate1(
819 resp, normalization_frequency))
821 if num.abs(normalization - 1.0) > 1e-2:
822 log.append((
823 'warning',
824 'FIR filter coefficients are not normalized. '
825 'Normalizing them. Factor: %g' % normalization,
826 context))
828 resp = response.DigitalFilterResponse(
829 (b/normalization).tolist(), [1.0], deltat,
830 drop_phase=drop_phase)
832 resps.append(resp)
834 if not drop_phase:
835 resps.extend(delay_responses)
837 else:
838 return Delivery(error=(
839 'ValueError',
840 'Unknown transfer function type: %s' % (
841 self.cf_transfer_function_type)))
843 return Delivery(payload=resps, log=log)
846class Latitude(FloatWithUnit):
847 '''
848 Type for latitude coordinate.
849 '''
851 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
852 # fixed unit
853 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
856class Longitude(FloatWithUnit):
857 '''
858 Type for longitude coordinate.
859 '''
861 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
862 # fixed unit
863 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
866class PolesZeros(BaseFilter):
867 '''
868 Response: complex poles and zeros. Corresponds to SEED blockette 53.
869 '''
871 pz_transfer_function_type = PzTransferFunction.T(
872 xmltagname='PzTransferFunctionType')
873 normalization_factor = Float.T(default=1.0,
874 xmltagname='NormalizationFactor')
875 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
876 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
877 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
879 def summary(self):
880 return 'pz_%s(%i,%i)' % (
881 'ABC?'[
882 PzTransferFunction.choices.index(
883 self.pz_transfer_function_type)],
884 len(self.pole_list),
885 len(self.zero_list))
887 def get_pyrocko_response(self, context='', deltat=None):
889 context += self.summary()
891 factor = 1.0
892 cfactor = 1.0
893 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
894 factor = 2. * math.pi
895 cfactor = (2. * math.pi)**(
896 len(self.pole_list) - len(self.zero_list))
898 log = []
899 if self.normalization_factor is None \
900 or self.normalization_factor == 0.0:
902 log.append((
903 'warning',
904 'No pole-zero normalization factor given. '
905 'Assuming a value of 1.0',
906 context))
908 nfactor = 1.0
909 else:
910 nfactor = self.normalization_factor
912 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
913 if not is_digital:
914 resp = response.PoleZeroResponse(
915 constant=nfactor*cfactor,
916 zeros=[z.value()*factor for z in self.zero_list],
917 poles=[p.value()*factor for p in self.pole_list])
918 else:
919 if not deltat:
920 log.append((
921 'error',
922 'Digital response requires knowledge about sampling '
923 'interval. Response will be unusable.',
924 context))
926 resp = response.DigitalPoleZeroResponse(
927 constant=nfactor*cfactor,
928 zeros=[z.value()*factor for z in self.zero_list],
929 poles=[p.value()*factor for p in self.pole_list],
930 deltat=deltat or 0.0)
932 if not self.normalization_frequency.value:
933 log.append((
934 'warning',
935 'Cannot check pole-zero normalization factor: '
936 'No normalization frequency given.',
937 context))
939 else:
940 if is_digital and not deltat:
941 log.append((
942 'warning',
943 'Cannot check computed vs reported normalization '
944 'factor without knowing the sampling interval.',
945 context))
946 else:
947 computed_normalization_factor = nfactor / abs(evaluate1(
948 resp, self.normalization_frequency.value))
950 db = 20.0 * num.log10(
951 computed_normalization_factor / nfactor)
953 if abs(db) > 0.17:
954 log.append((
955 'warning',
956 'Computed and reported normalization factors differ '
957 'by %g dB: computed: %g, reported: %g' % (
958 db,
959 computed_normalization_factor,
960 nfactor),
961 context))
963 return Delivery([resp], log)
966class ResponseListElement(Object):
967 frequency = Frequency.T(xmltagname='Frequency')
968 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
969 phase = Angle.T(xmltagname='Phase')
972class Polynomial(BaseFilter):
973 '''
974 Response: expressed as a polynomial (allows non-linear sensors to be
975 described). Corresponds to SEED blockette 62. Can be used to describe a
976 stage of acquisition or a complete system.
977 '''
979 approximation_type = Approximation.T(default='MACLAURIN',
980 xmltagname='ApproximationType')
981 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
982 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
983 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
984 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
985 maximum_error = Float.T(xmltagname='MaximumError')
986 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
988 def summary(self):
989 return 'poly(%i)' % len(self.coefficient_list)
992class Decimation(Object):
993 '''
994 Corresponds to SEED blockette 57.
995 '''
997 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
998 factor = Int.T(xmltagname='Factor')
999 offset = Int.T(xmltagname='Offset')
1000 delay = FloatWithUnit.T(xmltagname='Delay')
1001 correction = FloatWithUnit.T(xmltagname='Correction')
1003 def summary(self):
1004 return 'deci(%i, %g -> %g, %g)' % (
1005 self.factor,
1006 self.input_sample_rate.value,
1007 self.input_sample_rate.value / self.factor,
1008 self.delay.value)
1010 def get_pyrocko_response(self):
1011 if self.delay and self.delay.value != 0.0:
1012 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1013 else:
1014 return Delivery([])
1017class Operator(Object):
1018 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1019 contact_list = List.T(Person.T(xmltagname='Contact'))
1020 web_site = String.T(optional=True, xmltagname='WebSite')
1023class Comment(Object):
1024 '''
1025 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1026 and 59.
1027 '''
1029 id = Counter.T(optional=True, xmlstyle='attribute')
1030 value = Unicode.T(xmltagname='Value')
1031 begin_effective_time = DummyAwareOptionalTimestamp.T(
1032 optional=True,
1033 xmltagname='BeginEffectiveTime')
1034 end_effective_time = DummyAwareOptionalTimestamp.T(
1035 optional=True,
1036 xmltagname='EndEffectiveTime')
1037 author_list = List.T(Person.T(xmltagname='Author'))
1040class ResponseList(BaseFilter):
1041 '''
1042 Response: list of frequency, amplitude and phase values. Corresponds to
1043 SEED blockette 55.
1044 '''
1046 response_list_element_list = List.T(
1047 ResponseListElement.T(xmltagname='ResponseListElement'))
1049 def summary(self):
1050 return 'list(%i)' % len(self.response_list_element_list)
1053class Log(Object):
1054 '''
1055 Container for log entries.
1056 '''
1058 entry_list = List.T(Comment.T(xmltagname='Entry'))
1061class ResponseStage(Object):
1062 '''
1063 This complex type represents channel response and covers SEED blockettes 53
1064 to 56.
1065 '''
1067 number = Counter.T(xmlstyle='attribute')
1068 resource_id = String.T(optional=True, xmlstyle='attribute')
1069 poles_zeros_list = List.T(
1070 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1071 coefficients_list = List.T(
1072 Coefficients.T(optional=True, xmltagname='Coefficients'))
1073 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1074 fir = FIR.T(optional=True, xmltagname='FIR')
1075 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1076 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1077 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1079 def summary(self):
1080 elements = []
1082 for stuff in [
1083 self.poles_zeros_list, self.coefficients_list,
1084 self.response_list, self.fir, self.polynomial,
1085 self.decimation, self.stage_gain]:
1087 if stuff:
1088 if isinstance(stuff, Object):
1089 elements.append(stuff.summary())
1090 else:
1091 elements.extend(obj.summary() for obj in stuff)
1093 return '%i: %s %s -> %s' % (
1094 self.number,
1095 ', '.join(elements),
1096 self.input_units.name.upper() if self.input_units else '?',
1097 self.output_units.name.upper() if self.output_units else '?')
1099 def get_squirrel_response_stage(self, context):
1100 from pyrocko.squirrel.model import ResponseStage
1102 delivery = Delivery()
1103 delivery_pr = self.get_pyrocko_response(context)
1104 log = delivery_pr.log
1105 delivery_pr.log = []
1106 elements = delivery.extend_without_payload(delivery_pr)
1108 def to_quantity(unit):
1109 trans = {
1110 'M/S': 'velocity',
1111 'M': 'displacement',
1112 'M/S**2': 'acceleration',
1113 'V': 'voltage',
1114 'COUNTS': 'counts',
1115 'COUNT': 'counts',
1116 'PA': 'pressure',
1117 'RAD/S': 'rotation'}
1119 if unit is None:
1120 return None
1122 name = unit.name.upper()
1123 if name in trans:
1124 return trans[name]
1125 else:
1126 delivery.log.append((
1127 'warning',
1128 'Units not supported by Squirrel framework: %s' % (
1129 unit.name.upper() + (
1130 ' (%s)' % unit.description
1131 if unit.description else '')),
1132 context))
1134 return 'unsupported_quantity(%s)' % unit
1136 delivery.payload = [ResponseStage(
1137 input_quantity=to_quantity(self.input_units),
1138 output_quantity=to_quantity(self.output_units),
1139 input_sample_rate=self.input_sample_rate,
1140 output_sample_rate=self.output_sample_rate,
1141 elements=elements,
1142 log=log)]
1144 return delivery
1146 def get_pyrocko_response(self, context, gain_only=False):
1148 context = context + ', stage %i' % self.number
1150 responses = []
1151 log = []
1152 if self.stage_gain:
1153 normalization_frequency = self.stage_gain.frequency or 0.0
1154 else:
1155 normalization_frequency = 0.0
1157 if not gain_only:
1158 deltat = None
1159 delay_responses = []
1160 if self.decimation:
1161 rate = self.decimation.input_sample_rate.value
1162 if rate > 0.0:
1163 deltat = 1.0 / rate
1164 delivery = self.decimation.get_pyrocko_response()
1165 if delivery.errors:
1166 return delivery
1168 delay_responses = delivery.payload
1169 log.extend(delivery.log)
1171 for pzs in self.poles_zeros_list:
1172 delivery = pzs.get_pyrocko_response(context, deltat)
1173 if delivery.errors:
1174 return delivery
1176 pz_resps = delivery.payload
1177 log.extend(delivery.log)
1178 responses.extend(pz_resps)
1180 # emulate incorrect? evalresp behaviour
1181 if pzs.normalization_frequency != normalization_frequency \
1182 and normalization_frequency != 0.0:
1184 try:
1185 trial = response.MultiplyResponse(pz_resps)
1186 anorm = num.abs(evaluate1(
1187 trial, pzs.normalization_frequency.value))
1188 asens = num.abs(
1189 evaluate1(trial, normalization_frequency))
1191 factor = anorm/asens
1193 if abs(factor - 1.0) > 0.01:
1194 log.append((
1195 'warning',
1196 'PZ normalization frequency (%g) is different '
1197 'from stage gain frequency (%s) -> Emulating '
1198 'possibly incorrect evalresp behaviour. '
1199 'Correction factor: %g' % (
1200 pzs.normalization_frequency.value,
1201 normalization_frequency,
1202 factor),
1203 context))
1205 responses.append(
1206 response.PoleZeroResponse(constant=factor))
1207 except response.InvalidResponseError as e:
1208 log.append((
1209 'warning',
1210 'Could not check response: %s' % str(e),
1211 context))
1213 if len(self.poles_zeros_list) > 1:
1214 log.append((
1215 'warning',
1216 'Multiple poles and zeros records in single response '
1217 'stage.',
1218 context))
1220 for cfs in self.coefficients_list + (
1221 [self.fir] if self.fir else []):
1223 delivery = cfs.get_pyrocko_response(
1224 context, deltat, delay_responses,
1225 normalization_frequency)
1227 if delivery.errors:
1228 return delivery
1230 responses.extend(delivery.payload)
1231 log.extend(delivery.log)
1233 if len(self.coefficients_list) > 1:
1234 log.append((
1235 'warning',
1236 'Multiple filter coefficients lists in single response '
1237 'stage.',
1238 context))
1240 if self.response_list:
1241 log.append((
1242 'warning',
1243 'Unhandled response element of type: ResponseList',
1244 context))
1246 if self.polynomial:
1247 log.append((
1248 'warning',
1249 'Unhandled response element of type: Polynomial',
1250 context))
1252 if self.stage_gain:
1253 responses.append(
1254 response.PoleZeroResponse(constant=self.stage_gain.value))
1256 return Delivery(responses, log)
1258 @property
1259 def input_units(self):
1260 for e in (self.poles_zeros_list + self.coefficients_list +
1261 [self.response_list, self.fir, self.polynomial]):
1262 if e is not None:
1263 return e.input_units
1265 return None
1267 @property
1268 def output_units(self):
1269 for e in (self.poles_zeros_list + self.coefficients_list +
1270 [self.response_list, self.fir, self.polynomial]):
1271 if e is not None:
1272 return e.output_units
1274 return None
1276 @property
1277 def input_sample_rate(self):
1278 if self.decimation:
1279 return self.decimation.input_sample_rate.value
1281 return None
1283 @property
1284 def output_sample_rate(self):
1285 if self.decimation:
1286 return self.decimation.input_sample_rate.value \
1287 / self.decimation.factor
1289 return None
1292class Response(Object):
1293 resource_id = String.T(optional=True, xmlstyle='attribute')
1294 instrument_sensitivity = Sensitivity.T(optional=True,
1295 xmltagname='InstrumentSensitivity')
1296 instrument_polynomial = Polynomial.T(optional=True,
1297 xmltagname='InstrumentPolynomial')
1298 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1300 def check_sample_rates(self, channel):
1302 if self.stage_list:
1303 sample_rate = None
1305 for stage in self.stage_list:
1306 if stage.decimation:
1307 input_sample_rate = \
1308 stage.decimation.input_sample_rate.value
1310 if sample_rate is not None and not same_sample_rate(
1311 sample_rate, input_sample_rate):
1313 logger.warning(
1314 'Response stage %i has unexpected input sample '
1315 'rate: %g Hz (expected: %g Hz)' % (
1316 stage.number,
1317 input_sample_rate,
1318 sample_rate))
1320 sample_rate = input_sample_rate / stage.decimation.factor
1322 if sample_rate is not None and channel.sample_rate \
1323 and not same_sample_rate(
1324 sample_rate, channel.sample_rate.value):
1326 logger.warning(
1327 'Channel sample rate (%g Hz) does not match sample rate '
1328 'deducted from response stages information (%g Hz).' % (
1329 channel.sample_rate.value,
1330 sample_rate))
1332 def check_units(self):
1334 if self.instrument_sensitivity \
1335 and self.instrument_sensitivity.input_units:
1337 units = self.instrument_sensitivity.input_units.name.upper()
1339 if self.stage_list:
1340 for stage in self.stage_list:
1341 if units and stage.input_units \
1342 and stage.input_units.name.upper() != units:
1344 logger.warning(
1345 'Input units of stage %i (%s) do not match %s (%s).'
1346 % (
1347 stage.number,
1348 units,
1349 'output units of stage %i'
1350 if stage.number == 0
1351 else 'sensitivity input units',
1352 units))
1354 if stage.output_units:
1355 units = stage.output_units.name.upper()
1356 else:
1357 units = None
1359 sout_units = self.instrument_sensitivity.output_units
1360 if self.instrument_sensitivity and sout_units:
1361 if units is not None and units != sout_units.name.upper():
1362 logger.warning(
1363 'Output units of stage %i (%s) do not match %s (%s).'
1364 % (
1365 stage.number,
1366 units,
1367 'sensitivity output units',
1368 sout_units.name.upper()))
1370 def _sensitivity_checkpoints(self, responses, context):
1371 delivery = Delivery()
1373 if self.instrument_sensitivity:
1374 sval = self.instrument_sensitivity.value
1375 sfreq = self.instrument_sensitivity.frequency
1376 if sval is None:
1377 delivery.log.append((
1378 'warning',
1379 'No sensitivity value given.',
1380 context))
1382 elif sval is None:
1383 delivery.log.append((
1384 'warning',
1385 'Reported sensitivity value is zero.',
1386 context))
1388 elif sfreq is None:
1389 delivery.log.append((
1390 'warning',
1391 'Sensitivity frequency not given.',
1392 context))
1394 else:
1395 trial = response.MultiplyResponse(responses)
1397 delivery.extend(
1398 check_resp(
1399 trial, sval, sfreq, 0.1,
1400 'Instrument sensitivity value inconsistent with '
1401 'sensitivity computed from complete response',
1402 context))
1404 delivery.payload.append(response.FrequencyResponseCheckpoint(
1405 frequency=sfreq, value=sval))
1407 return delivery
1409 def get_squirrel_response(self, context, **kwargs):
1410 from pyrocko.squirrel.model import Response
1412 if self.stage_list:
1413 delivery = Delivery()
1414 for istage, stage in enumerate(self.stage_list):
1415 delivery.extend(stage.get_squirrel_response_stage(context))
1417 checkpoints = []
1418 if not delivery.errors:
1419 all_responses = []
1420 for sq_stage in delivery.payload:
1421 all_responses.extend(sq_stage.elements)
1423 checkpoints.extend(delivery.extend_without_payload(
1424 self._sensitivity_checkpoints(all_responses, context)))
1426 sq_stages = delivery.expect()
1428 return Response(
1429 stages=sq_stages,
1430 log=delivery.log,
1431 checkpoints=checkpoints,
1432 **kwargs)
1434 elif self.instrument_sensitivity:
1435 raise NoResponseInformation(
1436 'Only instrument sensitivity given (won\'t use it). (%s).'
1437 % context)
1438 else:
1439 raise NoResponseInformation(
1440 'Empty instrument response (%s).'
1441 % context)
1443 def get_pyrocko_response(
1444 self, context, fake_input_units=None, stages=(0, 1)):
1446 delivery = Delivery()
1447 if self.stage_list:
1448 for istage, stage in enumerate(self.stage_list):
1449 delivery.extend(stage.get_pyrocko_response(
1450 context, gain_only=not (
1451 stages is None or stages[0] <= istage < stages[1])))
1453 elif self.instrument_sensitivity:
1454 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1456 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1457 checkpoints = delivery.extend_without_payload(delivery_cp)
1458 if delivery.errors:
1459 return delivery
1461 if fake_input_units is not None:
1462 if not self.instrument_sensitivity or \
1463 self.instrument_sensitivity.input_units is None:
1465 delivery.errors.append((
1466 'NoResponseInformation',
1467 'No input units given, so cannot convert to requested '
1468 'units: %s' % fake_input_units.upper(),
1469 context))
1471 return delivery
1473 input_units = self.instrument_sensitivity.input_units.name.upper()
1475 conresp = None
1476 try:
1477 conresp = conversion[
1478 fake_input_units.upper(), input_units]
1480 except KeyError:
1481 delivery.errors.append((
1482 'NoResponseInformation',
1483 'Cannot convert between units: %s, %s'
1484 % (fake_input_units, input_units),
1485 context))
1487 if conresp is not None:
1488 delivery.payload.append(conresp)
1489 for checkpoint in checkpoints:
1490 checkpoint.value *= num.abs(evaluate1(
1491 conresp, checkpoint.frequency))
1493 delivery.payload = [
1494 response.MultiplyResponse(
1495 delivery.payload,
1496 checkpoints=checkpoints)]
1498 return delivery
1500 @classmethod
1501 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1502 normalization_frequency=1.0):
1503 '''
1504 Convert Pyrocko pole-zero response to StationXML response.
1506 :param presponse: Pyrocko pole-zero response
1507 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1508 :param input_unit: Input unit to be reported in the StationXML
1509 response.
1510 :type input_unit: str
1511 :param output_unit: Output unit to be reported in the StationXML
1512 response.
1513 :type output_unit: str
1514 :param normalization_frequency: Frequency where the normalization
1515 factor for the StationXML response should be computed.
1516 :type normalization_frequency: float
1517 '''
1519 norm_factor = 1.0/float(abs(
1520 evaluate1(presponse, normalization_frequency)
1521 / presponse.constant))
1523 pzs = PolesZeros(
1524 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1525 normalization_factor=norm_factor,
1526 normalization_frequency=Frequency(normalization_frequency),
1527 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1528 imaginary=FloatNoUnit(z.imag))
1529 for z in presponse.zeros],
1530 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1531 imaginary=FloatNoUnit(z.imag))
1532 for z in presponse.poles])
1534 pzs.validate()
1536 stage = ResponseStage(
1537 number=1,
1538 poles_zeros_list=[pzs],
1539 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1541 resp = Response(
1542 instrument_sensitivity=Sensitivity(
1543 value=stage.stage_gain.value,
1544 frequency=normalization_frequency,
1545 input_units=Units(input_unit),
1546 output_units=Units(output_unit)),
1548 stage_list=[stage])
1550 return resp
1553class BaseNode(Object):
1554 '''
1555 A base node type for derivation from: Network, Station and Channel types.
1556 '''
1558 code = String.T(xmlstyle='attribute')
1559 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1560 xmlstyle='attribute')
1561 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1562 xmlstyle='attribute')
1563 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1564 alternate_code = String.T(optional=True, xmlstyle='attribute')
1565 historical_code = String.T(optional=True, xmlstyle='attribute')
1566 description = Unicode.T(optional=True, xmltagname='Description')
1567 comment_list = List.T(Comment.T(xmltagname='Comment'))
1569 def spans(self, *args):
1570 if len(args) == 0:
1571 return True
1572 elif len(args) == 1:
1573 return ((self.start_date is None or
1574 self.start_date <= args[0]) and
1575 (self.end_date is None or
1576 args[0] <= self.end_date))
1578 elif len(args) == 2:
1579 return ((self.start_date is None or
1580 args[1] >= self.start_date) and
1581 (self.end_date is None or
1582 self.end_date >= args[0]))
1585def overlaps(a, b):
1586 return (
1587 a.start_date is None or b.end_date is None
1588 or a.start_date < b.end_date
1589 ) and (
1590 b.start_date is None or a.end_date is None
1591 or b.start_date < a.end_date
1592 )
1595class Channel(BaseNode):
1596 '''
1597 Equivalent to SEED blockette 52 and parent element for the related the
1598 response blockettes.
1599 '''
1601 location_code = String.T(xmlstyle='attribute')
1602 external_reference_list = List.T(
1603 ExternalReference.T(xmltagname='ExternalReference'))
1604 latitude = Latitude.T(xmltagname='Latitude')
1605 longitude = Longitude.T(xmltagname='Longitude')
1606 elevation = Distance.T(xmltagname='Elevation')
1607 depth = Distance.T(xmltagname='Depth')
1608 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1609 dip = Dip.T(optional=True, xmltagname='Dip')
1610 type_list = List.T(Type.T(xmltagname='Type'))
1611 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1612 sample_rate_ratio = SampleRateRatio.T(optional=True,
1613 xmltagname='SampleRateRatio')
1614 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1615 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1616 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1617 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1618 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1619 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1620 equipment = Equipment.T(optional=True, xmltagname='Equipment')
1621 response = Response.T(optional=True, xmltagname='Response')
1623 @property
1624 def position_values(self):
1625 lat = self.latitude.value
1626 lon = self.longitude.value
1627 elevation = value_or_none(self.elevation)
1628 depth = value_or_none(self.depth)
1629 return lat, lon, elevation, depth
1632class Station(BaseNode):
1633 '''
1634 This type represents a Station epoch. It is common to only have a single
1635 station epoch with the station's creation and termination dates as the
1636 epoch start and end dates.
1637 '''
1639 latitude = Latitude.T(xmltagname='Latitude')
1640 longitude = Longitude.T(xmltagname='Longitude')
1641 elevation = Distance.T(xmltagname='Elevation')
1642 site = Site.T(optional=True, xmltagname='Site')
1643 vault = Unicode.T(optional=True, xmltagname='Vault')
1644 geology = Unicode.T(optional=True, xmltagname='Geology')
1645 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1646 operator_list = List.T(Operator.T(xmltagname='Operator'))
1647 creation_date = DummyAwareOptionalTimestamp.T(
1648 optional=True, xmltagname='CreationDate')
1649 termination_date = DummyAwareOptionalTimestamp.T(
1650 optional=True, xmltagname='TerminationDate')
1651 total_number_channels = Counter.T(
1652 optional=True, xmltagname='TotalNumberChannels')
1653 selected_number_channels = Counter.T(
1654 optional=True, xmltagname='SelectedNumberChannels')
1655 external_reference_list = List.T(
1656 ExternalReference.T(xmltagname='ExternalReference'))
1657 channel_list = List.T(Channel.T(xmltagname='Channel'))
1659 @property
1660 def position_values(self):
1661 lat = self.latitude.value
1662 lon = self.longitude.value
1663 elevation = value_or_none(self.elevation)
1664 return lat, lon, elevation
1667class Network(BaseNode):
1668 '''
1669 This type represents the Network layer, all station metadata is contained
1670 within this element. The official name of the network or other descriptive
1671 information can be included in the Description element. The Network can
1672 contain 0 or more Stations.
1673 '''
1675 total_number_stations = Counter.T(optional=True,
1676 xmltagname='TotalNumberStations')
1677 selected_number_stations = Counter.T(optional=True,
1678 xmltagname='SelectedNumberStations')
1679 station_list = List.T(Station.T(xmltagname='Station'))
1681 @property
1682 def station_code_list(self):
1683 return sorted(set(s.code for s in self.station_list))
1685 @property
1686 def sl_code_list(self):
1687 sls = set()
1688 for station in self.station_list:
1689 for channel in station.channel_list:
1690 sls.add((station.code, channel.location_code))
1692 return sorted(sls)
1694 def summary(self, width=80, indent=4):
1695 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1696 lines = ['%s (%i):' % (self.code, len(sls))]
1697 if sls:
1698 ssls = ['.'.join(x for x in c if x) for c in sls]
1699 w = max(len(x) for x in ssls)
1700 n = (width - indent) / (w+1)
1701 while ssls:
1702 lines.append(
1703 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1705 ssls[:n] = []
1707 return '\n'.join(lines)
1710def value_or_none(x):
1711 if x is not None:
1712 return x.value
1713 else:
1714 return None
1717def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1719 pos = lat, lon, elevation, depth = \
1720 channels[0].position_values
1722 if not all(pos == x.position_values for x in channels):
1723 info = '\n'.join(
1724 ' %s: %s' % (x.code, x.position_values) for
1725 x in channels)
1727 mess = 'encountered inconsistencies in channel ' \
1728 'lat/lon/elevation/depth ' \
1729 'for %s.%s.%s: \n%s' % (nsl + (info,))
1731 if inconsistencies == 'raise':
1732 raise InconsistentChannelLocations(mess)
1734 elif inconsistencies == 'warn':
1735 logger.warning(mess)
1736 logger.warning(' -> using mean values')
1738 apos = num.array([x.position_values for x in channels], dtype=float)
1739 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1740 / num.sum(num.isfinite(apos), axis=0)
1742 groups = {}
1743 for channel in channels:
1744 if channel.code not in groups:
1745 groups[channel.code] = []
1747 groups[channel.code].append(channel)
1749 pchannels = []
1750 for code in sorted(groups.keys()):
1751 data = [
1752 (channel.code, value_or_none(channel.azimuth),
1753 value_or_none(channel.dip))
1754 for channel in groups[code]]
1756 azimuth, dip = util.consistency_merge(
1757 data,
1758 message='channel orientation values differ:',
1759 error=inconsistencies)
1761 pchannels.append(
1762 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1764 return pyrocko.model.Station(
1765 *nsl,
1766 lat=mlat,
1767 lon=mlon,
1768 elevation=mele,
1769 depth=mdep,
1770 channels=pchannels)
1773class FDSNStationXML(Object):
1774 '''
1775 Top-level type for Station XML. Required field are Source (network ID of
1776 the institution sending the message) and one or more Network containers or
1777 one or more Station containers.
1778 '''
1780 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1781 source = String.T(xmltagname='Source')
1782 sender = String.T(optional=True, xmltagname='Sender')
1783 module = String.T(optional=True, xmltagname='Module')
1784 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1785 created = Timestamp.T(optional=True, xmltagname='Created')
1786 network_list = List.T(Network.T(xmltagname='Network'))
1788 xmltagname = 'FDSNStationXML'
1789 guessable_xmlns = [guts_xmlns]
1791 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1792 time=None, timespan=None,
1793 inconsistencies='warn'):
1795 assert inconsistencies in ('raise', 'warn')
1797 if nslcs is not None:
1798 nslcs = set(nslcs)
1800 if nsls is not None:
1801 nsls = set(nsls)
1803 tt = ()
1804 if time is not None:
1805 tt = (time,)
1806 elif timespan is not None:
1807 tt = timespan
1809 pstations = []
1810 for network in self.network_list:
1811 if not network.spans(*tt):
1812 continue
1814 for station in network.station_list:
1815 if not station.spans(*tt):
1816 continue
1818 if station.channel_list:
1819 loc_to_channels = {}
1820 for channel in station.channel_list:
1821 if not channel.spans(*tt):
1822 continue
1824 loc = channel.location_code.strip()
1825 if loc not in loc_to_channels:
1826 loc_to_channels[loc] = []
1828 loc_to_channels[loc].append(channel)
1830 for loc in sorted(loc_to_channels.keys()):
1831 channels = loc_to_channels[loc]
1832 if nslcs is not None:
1833 channels = [channel for channel in channels
1834 if (network.code, station.code, loc,
1835 channel.code) in nslcs]
1837 if not channels:
1838 continue
1840 nsl = network.code, station.code, loc
1841 if nsls is not None and nsl not in nsls:
1842 continue
1844 pstations.append(
1845 pyrocko_station_from_channels(
1846 nsl,
1847 channels,
1848 inconsistencies=inconsistencies))
1849 else:
1850 pstations.append(pyrocko.model.Station(
1851 network.code, station.code, '*',
1852 lat=station.latitude.value,
1853 lon=station.longitude.value,
1854 elevation=value_or_none(station.elevation),
1855 name=station.description or ''))
1857 return pstations
1859 @classmethod
1860 def from_pyrocko_stations(
1861 cls, pyrocko_stations, add_flat_responses_from=None):
1863 '''
1864 Generate :py:class:`FDSNStationXML` from list of
1865 :py:class;`pyrocko.model.Station` instances.
1867 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1868 instances.
1869 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1870 '''
1871 from collections import defaultdict
1872 network_dict = defaultdict(list)
1874 if add_flat_responses_from:
1875 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1876 extra = dict(
1877 response=Response(
1878 instrument_sensitivity=Sensitivity(
1879 value=1.0,
1880 frequency=1.0,
1881 input_units=Units(name=add_flat_responses_from))))
1882 else:
1883 extra = {}
1885 have_offsets = set()
1886 for s in pyrocko_stations:
1888 if s.north_shift != 0.0 or s.east_shift != 0.0:
1889 have_offsets.add(s.nsl())
1891 network, station, location = s.nsl()
1892 channel_list = []
1893 for c in s.channels:
1894 channel_list.append(
1895 Channel(
1896 location_code=location,
1897 code=c.name,
1898 latitude=Latitude(value=s.effective_lat),
1899 longitude=Longitude(value=s.effective_lon),
1900 elevation=Distance(value=s.elevation),
1901 depth=Distance(value=s.depth),
1902 azimuth=Azimuth(value=c.azimuth),
1903 dip=Dip(value=c.dip),
1904 **extra
1905 )
1906 )
1908 network_dict[network].append(
1909 Station(
1910 code=station,
1911 latitude=Latitude(value=s.effective_lat),
1912 longitude=Longitude(value=s.effective_lon),
1913 elevation=Distance(value=s.elevation),
1914 channel_list=channel_list)
1915 )
1917 if have_offsets:
1918 logger.warning(
1919 'StationXML does not support Cartesian offsets in '
1920 'coordinates. Storing effective lat/lon for stations: %s' %
1921 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1923 timestamp = util.to_time_float(time.time())
1924 network_list = []
1925 for k, station_list in network_dict.items():
1927 network_list.append(
1928 Network(
1929 code=k, station_list=station_list,
1930 total_number_stations=len(station_list)))
1932 sxml = FDSNStationXML(
1933 source='from pyrocko stations list', created=timestamp,
1934 network_list=network_list)
1936 sxml.validate()
1937 return sxml
1939 def iter_network_stations(
1940 self, net=None, sta=None, time=None, timespan=None):
1942 tt = ()
1943 if time is not None:
1944 tt = (time,)
1945 elif timespan is not None:
1946 tt = timespan
1948 for network in self.network_list:
1949 if not network.spans(*tt) or (
1950 net is not None and network.code != net):
1951 continue
1953 for station in network.station_list:
1954 if not station.spans(*tt) or (
1955 sta is not None and station.code != sta):
1956 continue
1958 yield (network, station)
1960 def iter_network_station_channels(
1961 self, net=None, sta=None, loc=None, cha=None,
1962 time=None, timespan=None):
1964 if loc is not None:
1965 loc = loc.strip()
1967 tt = ()
1968 if time is not None:
1969 tt = (time,)
1970 elif timespan is not None:
1971 tt = timespan
1973 for network in self.network_list:
1974 if not network.spans(*tt) or (
1975 net is not None and network.code != net):
1976 continue
1978 for station in network.station_list:
1979 if not station.spans(*tt) or (
1980 sta is not None and station.code != sta):
1981 continue
1983 if station.channel_list:
1984 for channel in station.channel_list:
1985 if (not channel.spans(*tt) or
1986 (cha is not None and channel.code != cha) or
1987 (loc is not None and
1988 channel.location_code.strip() != loc)):
1989 continue
1991 yield (network, station, channel)
1993 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
1994 time=None, timespan=None):
1996 groups = {}
1997 for network, station, channel in self.iter_network_station_channels(
1998 net, sta, loc, cha, time=time, timespan=timespan):
2000 net = network.code
2001 sta = station.code
2002 cha = channel.code
2003 loc = channel.location_code.strip()
2004 if len(cha) == 3:
2005 bic = cha[:2] # band and intrument code according to SEED
2006 elif len(cha) == 1:
2007 bic = ''
2008 else:
2009 bic = cha
2011 if channel.response and \
2012 channel.response.instrument_sensitivity and \
2013 channel.response.instrument_sensitivity.input_units:
2015 unit = channel.response.instrument_sensitivity\
2016 .input_units.name.upper()
2017 else:
2018 unit = None
2020 bic = (bic, unit)
2022 k = net, sta, loc
2023 if k not in groups:
2024 groups[k] = {}
2026 if bic not in groups[k]:
2027 groups[k][bic] = []
2029 groups[k][bic].append(channel)
2031 for nsl, bic_to_channels in groups.items():
2032 bad_bics = []
2033 for bic, channels in bic_to_channels.items():
2034 sample_rates = []
2035 for channel in channels:
2036 sample_rates.append(channel.sample_rate.value)
2038 if not same(sample_rates):
2039 scs = ','.join(channel.code for channel in channels)
2040 srs = ', '.join('%e' % x for x in sample_rates)
2041 err = 'ignoring channels with inconsistent sampling ' + \
2042 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2044 logger.warning(err)
2045 bad_bics.append(bic)
2047 for bic in bad_bics:
2048 del bic_to_channels[bic]
2050 return groups
2052 def choose_channels(
2053 self,
2054 target_sample_rate=None,
2055 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2056 priority_units=['M/S', 'M/S**2'],
2057 priority_instrument_code=['H', 'L'],
2058 time=None,
2059 timespan=None):
2061 nslcs = {}
2062 for nsl, bic_to_channels in self.get_channel_groups(
2063 time=time, timespan=timespan).items():
2065 useful_bics = []
2066 for bic, channels in bic_to_channels.items():
2067 rate = channels[0].sample_rate.value
2069 if target_sample_rate is not None and \
2070 rate < target_sample_rate*0.99999:
2071 continue
2073 if len(bic[0]) == 2:
2074 if bic[0][0] not in priority_band_code:
2075 continue
2077 if bic[0][1] not in priority_instrument_code:
2078 continue
2080 unit = bic[1]
2082 prio_unit = len(priority_units)
2083 try:
2084 prio_unit = priority_units.index(unit)
2085 except ValueError:
2086 pass
2088 prio_inst = len(priority_instrument_code)
2089 prio_band = len(priority_band_code)
2090 if len(channels[0].code) == 3:
2091 try:
2092 prio_inst = priority_instrument_code.index(
2093 channels[0].code[1])
2094 except ValueError:
2095 pass
2097 try:
2098 prio_band = priority_band_code.index(
2099 channels[0].code[0])
2100 except ValueError:
2101 pass
2103 if target_sample_rate is None:
2104 rate = -rate
2106 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2107 prio_inst, bic))
2109 useful_bics.sort()
2111 for _, _, rate, _, _, bic in useful_bics:
2112 channels = sorted(
2113 bic_to_channels[bic],
2114 key=lambda channel: channel.code)
2116 if channels:
2117 for channel in channels:
2118 nslcs[nsl + (channel.code,)] = channel
2120 break
2122 return nslcs
2124 def get_pyrocko_response(
2125 self, nslc,
2126 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2128 net, sta, loc, cha = nslc
2129 resps = []
2130 for _, _, channel in self.iter_network_station_channels(
2131 net, sta, loc, cha, time=time, timespan=timespan):
2132 resp = channel.response
2133 if resp:
2134 resp.check_sample_rates(channel)
2135 resp.check_units()
2136 resps.append(resp.get_pyrocko_response(
2137 '.'.join(nslc),
2138 fake_input_units=fake_input_units,
2139 stages=stages).expect_one())
2141 if not resps:
2142 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2143 elif len(resps) > 1:
2144 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2146 return resps[0]
2148 @property
2149 def n_code_list(self):
2150 return sorted(set(x.code for x in self.network_list))
2152 @property
2153 def ns_code_list(self):
2154 nss = set()
2155 for network in self.network_list:
2156 for station in network.station_list:
2157 nss.add((network.code, station.code))
2159 return sorted(nss)
2161 @property
2162 def nsl_code_list(self):
2163 nsls = set()
2164 for network in self.network_list:
2165 for station in network.station_list:
2166 for channel in station.channel_list:
2167 nsls.add(
2168 (network.code, station.code, channel.location_code))
2170 return sorted(nsls)
2172 @property
2173 def nslc_code_list(self):
2174 nslcs = set()
2175 for network in self.network_list:
2176 for station in network.station_list:
2177 for channel in station.channel_list:
2178 nslcs.add(
2179 (network.code, station.code, channel.location_code,
2180 channel.code))
2182 return sorted(nslcs)
2184 def summary(self):
2185 lst = [
2186 'number of n codes: %i' % len(self.n_code_list),
2187 'number of ns codes: %i' % len(self.ns_code_list),
2188 'number of nsl codes: %i' % len(self.nsl_code_list),
2189 'number of nslc codes: %i' % len(self.nslc_code_list)
2190 ]
2191 return '\n'.join(lst)
2193 def summary_stages(self):
2194 data = []
2195 for network, station, channel in self.iter_network_station_channels():
2196 nslc = (network.code, station.code, channel.location_code,
2197 channel.code)
2199 stages = []
2200 in_units = '?'
2201 out_units = '?'
2202 if channel.response:
2203 sens = channel.response.instrument_sensitivity
2204 if sens:
2205 in_units = sens.input_units.name.upper()
2206 out_units = sens.output_units.name.upper()
2208 for stage in channel.response.stage_list:
2209 stages.append(stage.summary())
2211 data.append(
2212 (nslc, tts(channel.start_date), tts(channel.end_date),
2213 in_units, out_units, stages))
2215 data.sort()
2217 lst = []
2218 for nslc, stmin, stmax, in_units, out_units, stages in data:
2219 lst.append(' %s: %s - %s, %s -> %s' % (
2220 '.'.join(nslc), stmin, stmax, in_units, out_units))
2221 for stage in stages:
2222 lst.append(' %s' % stage)
2224 return '\n'.join(lst)
2226 def _check_overlaps(self):
2227 by_nslc = {}
2228 for network in self.network_list:
2229 for station in network.station_list:
2230 for channel in station.channel_list:
2231 nslc = (network.code, station.code, channel.location_code,
2232 channel.code)
2233 if nslc not in by_nslc:
2234 by_nslc[nslc] = []
2236 by_nslc[nslc].append(channel)
2238 errors = []
2239 for nslc, channels in by_nslc.items():
2240 for ia, a in enumerate(channels):
2241 for b in channels[ia+1:]:
2242 if overlaps(a, b):
2243 errors.append(
2244 'Channel epochs overlap for %s:\n'
2245 ' %s - %s\n %s - %s' % (
2246 '.'.join(nslc),
2247 tts(a.start_date), tts(a.end_date),
2248 tts(b.start_date), tts(b.end_date)))
2249 return errors
2251 def check(self):
2252 errors = []
2253 for meth in [self._check_overlaps]:
2254 errors.extend(meth())
2256 if errors:
2257 raise Inconsistencies(
2258 'Inconsistencies found in StationXML:\n '
2259 + '\n '.join(errors))
2262def load_channel_table(stream):
2264 networks = {}
2265 stations = {}
2267 for line in stream:
2268 line = str(line.decode('ascii'))
2269 if line.startswith('#'):
2270 continue
2272 t = line.rstrip().split('|')
2274 if len(t) != 17:
2275 logger.warning('Invalid channel record: %s' % line)
2276 continue
2278 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2279 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2281 try:
2282 scale = float(scale)
2283 except ValueError:
2284 scale = None
2286 try:
2287 scale_freq = float(scale_freq)
2288 except ValueError:
2289 scale_freq = None
2291 try:
2292 depth = float(dep)
2293 except ValueError:
2294 depth = 0.0
2296 try:
2297 azi = float(azi)
2298 dip = float(dip)
2299 except ValueError:
2300 azi = None
2301 dip = None
2303 try:
2304 if net not in networks:
2305 network = Network(code=net)
2306 else:
2307 network = networks[net]
2309 if (net, sta) not in stations:
2310 station = Station(
2311 code=sta, latitude=lat, longitude=lon, elevation=ele)
2313 station.regularize()
2314 else:
2315 station = stations[net, sta]
2317 if scale:
2318 resp = Response(
2319 instrument_sensitivity=Sensitivity(
2320 value=scale,
2321 frequency=scale_freq,
2322 input_units=scale_units))
2323 else:
2324 resp = None
2326 channel = Channel(
2327 code=cha,
2328 location_code=loc.strip(),
2329 latitude=lat,
2330 longitude=lon,
2331 elevation=ele,
2332 depth=depth,
2333 azimuth=azi,
2334 dip=dip,
2335 sensor=Equipment(description=sens),
2336 response=resp,
2337 sample_rate=sample_rate,
2338 start_date=start_date,
2339 end_date=end_date or None)
2341 channel.regularize()
2343 except ValidationError:
2344 raise InvalidRecord(line)
2346 if net not in networks:
2347 networks[net] = network
2349 if (net, sta) not in stations:
2350 stations[net, sta] = station
2351 network.station_list.append(station)
2353 station.channel_list.append(channel)
2355 return FDSNStationXML(
2356 source='created from table input',
2357 created=time.time(),
2358 network_list=sorted(networks.values(), key=lambda x: x.code))
2361def primitive_merge(sxs):
2362 networks = []
2363 for sx in sxs:
2364 networks.extend(sx.network_list)
2366 return FDSNStationXML(
2367 source='merged from different sources',
2368 created=time.time(),
2369 network_list=copy.deepcopy(
2370 sorted(networks, key=lambda x: x.code)))
2373def split_channels(sx):
2374 for nslc in sx.nslc_code_list:
2375 network_list = sx.network_list
2376 network_list_filtered = [
2377 network for network in network_list
2378 if network.code == nslc[0]]
2380 for network in network_list_filtered:
2381 sx.network_list = [network]
2382 station_list = network.station_list
2383 station_list_filtered = [
2384 station for station in station_list
2385 if station.code == nslc[1]]
2387 for station in station_list_filtered:
2388 network.station_list = [station]
2389 channel_list = station.channel_list
2390 station.channel_list = [
2391 channel for channel in channel_list
2392 if (channel.location_code, channel.code) == nslc[2:4]]
2394 yield nslc, copy.deepcopy(sx)
2396 station.channel_list = channel_list
2398 network.station_list = station_list
2400 sx.network_list = network_list
2403if __name__ == '__main__':
2404 from optparse import OptionParser
2406 util.setup_logging('pyrocko.io.stationxml', 'warning')
2408 usage = \
2409 'python -m pyrocko.io.stationxml check|stats|stages ' \
2410 '<filename> [options]'
2412 description = '''Torture StationXML file.'''
2414 parser = OptionParser(
2415 usage=usage,
2416 description=description,
2417 formatter=util.BetterHelpFormatter())
2419 (options, args) = parser.parse_args(sys.argv[1:])
2421 if len(args) != 2:
2422 parser.print_help()
2423 sys.exit(1)
2425 action, path = args
2427 sx = load_xml(filename=path)
2428 if action == 'check':
2429 try:
2430 sx.check()
2431 except Inconsistencies as e:
2432 logger.error(e)
2433 sys.exit(1)
2435 elif action == 'stats':
2436 print(sx.summary())
2438 elif action == 'stages':
2439 print(sx.summary_stages())
2441 else:
2442 parser.print_help()
2443 sys.exit('unknown action: %s' % action)