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}
49unit_to_quantity = {
50 'M/S': 'velocity',
51 'M': 'displacement',
52 'M/S**2': 'acceleration',
53 'V': 'voltage',
54 'COUNTS': 'counts',
55 'COUNT': 'counts',
56 'PA': 'pressure',
57 'RAD/S': 'rotation'}
60def to_quantity(unit, context, delivery):
62 if unit is None:
63 return None
65 name = unit.name.upper()
66 if name in unit_to_quantity:
67 return unit_to_quantity[name]
68 else:
69 delivery.log.append((
70 'warning',
71 'Units not supported by Squirrel framework: %s' % (
72 unit.name.upper() + (
73 ' (%s)' % unit.description
74 if unit.description else '')),
75 context))
77 return 'unsupported_quantity(%s)' % unit
80class StationXMLError(Exception):
81 pass
84class Inconsistencies(StationXMLError):
85 pass
88class NoResponseInformation(StationXMLError):
89 pass
92class MultipleResponseInformation(StationXMLError):
93 pass
96class InconsistentResponseInformation(StationXMLError):
97 pass
100class InconsistentChannelLocations(StationXMLError):
101 pass
104class InvalidRecord(StationXMLError):
105 def __init__(self, line):
106 StationXMLError.__init__(self)
107 self._line = line
109 def __str__(self):
110 return 'Invalid record: "%s"' % self._line
113_exceptions = dict(
114 Inconsistencies=Inconsistencies,
115 NoResponseInformation=NoResponseInformation,
116 MultipleResponseInformation=MultipleResponseInformation,
117 InconsistentResponseInformation=InconsistentResponseInformation,
118 InconsistentChannelLocations=InconsistentChannelLocations,
119 InvalidRecord=InvalidRecord,
120 ValueError=ValueError)
123_logs = dict(
124 info=logger.info,
125 warning=logger.warning,
126 error=logger.error)
129class DeliveryError(StationXMLError):
130 pass
133class Delivery(Object):
135 def __init__(self, payload=None, log=None, errors=None, error=None):
136 if payload is None:
137 payload = []
139 if log is None:
140 log = []
142 if errors is None:
143 errors = []
145 if error is not None:
146 errors.append(error)
148 Object.__init__(self, payload=payload, log=log, errors=errors)
150 payload = List.T(Any.T())
151 log = List.T(Tuple.T(3, String.T()))
152 errors = List.T(Tuple.T(3, String.T()))
154 def extend(self, other):
155 self.payload.extend(other.payload)
156 self.log.extend(other.log)
157 self.errors.extend(other.errors)
159 def extend_without_payload(self, other):
160 self.log.extend(other.log)
161 self.errors.extend(other.errors)
162 return other.payload
164 def emit_log(self):
165 for name, message, context in self.log:
166 message = '%s: %s' % (context, message)
167 _logs[name](message)
169 def expect(self, quiet=False):
170 if not quiet:
171 self.emit_log()
173 if self.errors:
174 name, message, context = self.errors[0]
175 if context:
176 message += ' (%s)' % context
178 if len(self.errors) > 1:
179 message += ' Additional errors pending.'
181 raise _exceptions[name](message)
183 return self.payload
185 def expect_one(self, quiet=False):
186 payload = self.expect(quiet=quiet)
187 if len(payload) != 1:
188 raise DeliveryError(
189 'Expected 1 element but got %i.' % len(payload))
191 return payload[0]
194def wrap(s, width=80, indent=4):
195 words = s.split()
196 lines = []
197 t = []
198 n = 0
199 for w in words:
200 if n + len(w) >= width:
201 lines.append(' '.join(t))
202 n = indent
203 t = [' '*(indent-1)]
205 t.append(w)
206 n += len(w) + 1
208 lines.append(' '.join(t))
209 return '\n'.join(lines)
212def same(x, eps=0.0):
213 if any(type(x[0]) != type(r) for r in x):
214 return False
216 if isinstance(x[0], float):
217 return all(abs(r-x[0]) <= eps for r in x)
218 else:
219 return all(r == x[0] for r in x)
222def same_sample_rate(a, b, eps=1.0e-6):
223 return abs(a - b) < min(a, b)*eps
226def evaluate1(resp, f):
227 return resp.evaluate(num.array([f], dtype=float))[0]
230def check_resp(resp, value, frequency, limit_db, prelude, context):
232 try:
233 value_resp = num.abs(evaluate1(resp, frequency))
234 except response.InvalidResponseError as e:
235 return Delivery(log=[(
236 'warning',
237 'Could not check response: %s' % str(e),
238 context)])
240 if value_resp == 0.0:
241 return Delivery(log=[(
242 'warning',
243 '%s\n'
244 ' computed response is zero' % prelude,
245 context)])
247 diff_db = 20.0 * num.log10(value_resp/value)
249 if num.abs(diff_db) > limit_db:
250 return Delivery(log=[(
251 'warning',
252 '%s\n'
253 ' reported value: %g\n'
254 ' computed value: %g\n'
255 ' at frequency [Hz]: %g\n'
256 ' factor: %g\n'
257 ' difference [dB]: %g\n'
258 ' limit [dB]: %g' % (
259 prelude,
260 value,
261 value_resp,
262 frequency,
263 value_resp/value,
264 diff_db,
265 limit_db),
266 context)])
268 return Delivery()
271def tts(t):
272 if t is None:
273 return '?'
274 else:
275 return util.tts(t, format='%Y-%m-%d')
278this_year = time.gmtime()[0]
281class DummyAwareOptionalTimestamp(Object):
282 '''
283 Optional timestamp with support for some common placeholder values.
285 Some StationXML files contain arbitrary placeholder values for open end
286 intervals, like "2100-01-01". Depending on the time range supported by the
287 system, these dates are translated into ``None`` to prevent crashes with
288 this type.
289 '''
290 dummy_for = (hpfloat, float)
291 dummy_for_description = 'time_float'
293 class __T(TBase):
295 def regularize_extra(self, val):
296 time_float = get_time_float()
298 if isinstance(val, datetime.datetime):
299 tt = val.utctimetuple()
300 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
302 elif isinstance(val, datetime.date):
303 tt = val.timetuple()
304 val = time_float(calendar.timegm(tt))
306 elif isinstance(val, (str, newstr)):
307 val = val.strip()
309 tz_offset = 0
311 m = re_tz.search(val)
312 if m:
313 sh = m.group(2)
314 sm = m.group(4)
315 tz_offset = (int(sh)*3600 if sh else 0) \
316 + (int(sm)*60 if sm else 0)
318 val = re_tz.sub('', val)
320 if len(val) > 10 and val[10] == 'T':
321 val = val.replace('T', ' ', 1)
323 try:
324 val = util.str_to_time(val) - tz_offset
326 except util.TimeStrError:
327 year = int(val[:4])
328 if sys.maxsize > 2**32: # if we're on 64bit
329 if year > this_year + 100:
330 return None # StationXML contained a dummy date
332 if year < 1903: # for macOS, 1900-01-01 dummy dates
333 return None
335 else: # 32bit end of time is in 2038
336 if this_year < 2037 and year > 2037 or year < 1903:
337 return None # StationXML contained a dummy date
339 raise
341 elif isinstance(val, (int, float)):
342 val = time_float(val)
344 else:
345 raise ValidationError(
346 '%s: cannot convert "%s" to type %s' % (
347 self.xname(), val, time_float))
349 return val
351 def to_save(self, val):
352 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
353 .rstrip('0').rstrip('.')
355 def to_save_xml(self, val):
356 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
357 .rstrip('0').rstrip('.') + 'Z'
360class Nominal(StringChoice):
361 choices = [
362 'NOMINAL',
363 'CALCULATED']
366class Email(UnicodePattern):
367 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
370class RestrictedStatus(StringChoice):
371 choices = [
372 'open',
373 'closed',
374 'partial']
377class Type(StringChoice):
378 choices = [
379 'TRIGGERED',
380 'CONTINUOUS',
381 'HEALTH',
382 'GEOPHYSICAL',
383 'WEATHER',
384 'FLAG',
385 'SYNTHESIZED',
386 'INPUT',
387 'EXPERIMENTAL',
388 'MAINTENANCE',
389 'BEAM']
391 class __T(StringChoice.T):
392 def validate_extra(self, val):
393 if val not in self.choices:
394 logger.warning(
395 'channel type: "%s" is not a valid choice out of %s' %
396 (val, repr(self.choices)))
399class PzTransferFunction(StringChoice):
400 choices = [
401 'LAPLACE (RADIANS/SECOND)',
402 'LAPLACE (HERTZ)',
403 'DIGITAL (Z-TRANSFORM)']
406class Symmetry(StringChoice):
407 choices = [
408 'NONE',
409 'EVEN',
410 'ODD']
413class CfTransferFunction(StringChoice):
415 class __T(StringChoice.T):
416 def validate(self, val, regularize=False, depth=-1):
417 if regularize:
418 try:
419 val = str(val)
420 except ValueError:
421 raise ValidationError(
422 '%s: cannot convert to string %s' % (self.xname,
423 repr(val)))
425 val = self.dummy_cls.replacements.get(val, val)
427 self.validate_extra(val)
428 return val
430 choices = [
431 'ANALOG (RADIANS/SECOND)',
432 'ANALOG (HERTZ)',
433 'DIGITAL']
435 replacements = {
436 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
437 'ANALOG (HZ)': 'ANALOG (HERTZ)',
438 }
441class Approximation(StringChoice):
442 choices = [
443 'MACLAURIN']
446class PhoneNumber(StringPattern):
447 pattern = '[0-9]+-[0-9]+'
450class Site(Object):
451 '''
452 Description of a site location using name and optional geopolitical
453 boundaries (country, city, etc.).
454 '''
456 name = Unicode.T(default='', xmltagname='Name')
457 description = Unicode.T(optional=True, xmltagname='Description')
458 town = Unicode.T(optional=True, xmltagname='Town')
459 county = Unicode.T(optional=True, xmltagname='County')
460 region = Unicode.T(optional=True, xmltagname='Region')
461 country = Unicode.T(optional=True, xmltagname='Country')
464class ExternalReference(Object):
465 '''
466 This type contains a URI and description for external data that users may
467 want to reference in StationXML.
468 '''
470 uri = String.T(xmltagname='URI')
471 description = Unicode.T(xmltagname='Description')
474class Units(Object):
475 '''
476 A type to document units. Corresponds to SEED blockette 34.
477 '''
479 def __init__(self, name=None, **kwargs):
480 Object.__init__(self, name=name, **kwargs)
482 name = String.T(xmltagname='Name')
483 description = Unicode.T(optional=True, xmltagname='Description')
486class Counter(Int):
487 pass
490class SampleRateRatio(Object):
491 '''
492 Sample rate expressed as number of samples in a number of seconds.
493 '''
495 number_samples = Int.T(xmltagname='NumberSamples')
496 number_seconds = Int.T(xmltagname='NumberSeconds')
499class Gain(Object):
500 '''
501 Complex type for sensitivity and frequency ranges. This complex type can be
502 used to represent both overall sensitivities and individual stage gains.
503 The FrequencyRangeGroup is an optional construct that defines a pass band
504 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
505 is valid within the number of decibels specified in FrequencyDBVariation.
506 '''
508 def __init__(self, value=None, **kwargs):
509 Object.__init__(self, value=value, **kwargs)
511 value = Float.T(optional=True, xmltagname='Value')
512 frequency = Float.T(optional=True, xmltagname='Frequency')
514 def summary(self):
515 return 'gain(%g @ %g)' % (self.value, self.frequency)
518class NumeratorCoefficient(Object):
519 i = Int.T(optional=True, xmlstyle='attribute')
520 value = Float.T(xmlstyle='content')
523class FloatNoUnit(Object):
524 def __init__(self, value=None, **kwargs):
525 Object.__init__(self, value=value, **kwargs)
527 plus_error = Float.T(optional=True, xmlstyle='attribute')
528 minus_error = Float.T(optional=True, xmlstyle='attribute')
529 value = Float.T(xmlstyle='content')
532class FloatWithUnit(FloatNoUnit):
533 unit = String.T(optional=True, xmlstyle='attribute')
536class Equipment(Object):
537 resource_id = String.T(optional=True, xmlstyle='attribute')
538 type = String.T(optional=True, xmltagname='Type')
539 description = Unicode.T(optional=True, xmltagname='Description')
540 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
541 vendor = Unicode.T(optional=True, xmltagname='Vendor')
542 model = Unicode.T(optional=True, xmltagname='Model')
543 serial_number = String.T(optional=True, xmltagname='SerialNumber')
544 installation_date = DummyAwareOptionalTimestamp.T(
545 optional=True,
546 xmltagname='InstallationDate')
547 removal_date = DummyAwareOptionalTimestamp.T(
548 optional=True,
549 xmltagname='RemovalDate')
550 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
553class PhoneNumber(Object):
554 description = Unicode.T(optional=True, xmlstyle='attribute')
555 country_code = Int.T(optional=True, xmltagname='CountryCode')
556 area_code = Int.T(xmltagname='AreaCode')
557 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
560class BaseFilter(Object):
561 '''
562 The BaseFilter is derived by all filters.
563 '''
565 resource_id = String.T(optional=True, xmlstyle='attribute')
566 name = String.T(optional=True, xmlstyle='attribute')
567 description = Unicode.T(optional=True, xmltagname='Description')
568 input_units = Units.T(optional=True, xmltagname='InputUnits')
569 output_units = Units.T(optional=True, xmltagname='OutputUnits')
572class Sensitivity(Gain):
573 '''
574 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
575 construct that defines a pass band in Hertz (FrequencyStart and
576 FrequencyEnd) in which the SensitivityValue is valid within the number of
577 decibels specified in FrequencyDBVariation.
578 '''
580 input_units = Units.T(optional=True, xmltagname='InputUnits')
581 output_units = Units.T(optional=True, xmltagname='OutputUnits')
582 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
583 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
584 frequency_db_variation = Float.T(optional=True,
585 xmltagname='FrequencyDBVariation')
587 def get_pyrocko_response(self):
588 return Delivery(
589 [response.PoleZeroResponse(constant=self.value)])
592class Coefficient(FloatNoUnit):
593 number = Counter.T(optional=True, xmlstyle='attribute')
596class PoleZero(Object):
597 '''
598 Complex numbers used as poles or zeros in channel response.
599 '''
601 number = Int.T(optional=True, xmlstyle='attribute')
602 real = FloatNoUnit.T(xmltagname='Real')
603 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
605 def value(self):
606 return self.real.value + 1J * self.imaginary.value
609class ClockDrift(FloatWithUnit):
610 unit = String.T(default='SECONDS/SAMPLE', optional=True,
611 xmlstyle='attribute') # fixed
614class Second(FloatWithUnit):
615 '''
616 A time value in seconds.
617 '''
619 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
620 # fixed unit
623class Voltage(FloatWithUnit):
624 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
625 # fixed unit
628class Angle(FloatWithUnit):
629 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
630 # fixed unit
633class Azimuth(FloatWithUnit):
634 '''
635 Instrument azimuth, degrees clockwise from North.
636 '''
638 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
639 # fixed unit
642class Dip(FloatWithUnit):
643 '''
644 Instrument dip in degrees down from horizontal. Together azimuth and dip
645 describe the direction of the sensitive axis of the instrument.
646 '''
648 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
649 # fixed unit
652class Distance(FloatWithUnit):
653 '''
654 Extension of FloatWithUnit for distances, elevations, and depths.
655 '''
657 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
658 # NOT fixed unit!
661class Frequency(FloatWithUnit):
662 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
663 # fixed unit
666class SampleRate(FloatWithUnit):
667 '''
668 Sample rate in samples per second.
669 '''
671 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
672 # fixed unit
675class Person(Object):
676 '''
677 Representation of a person's contact information. A person can belong to
678 multiple agencies and have multiple email addresses and phone numbers.
679 '''
681 name_list = List.T(Unicode.T(xmltagname='Name'))
682 agency_list = List.T(Unicode.T(xmltagname='Agency'))
683 email_list = List.T(Email.T(xmltagname='Email'))
684 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
687class FIR(BaseFilter):
688 '''
689 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
690 also commonly documented using the Coefficients element.
691 '''
693 symmetry = Symmetry.T(xmltagname='Symmetry')
694 numerator_coefficient_list = List.T(
695 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
697 def summary(self):
698 return 'fir(%i%s)' % (
699 self.get_ncoefficiencs(),
700 ',sym' if self.get_effective_symmetry() != 'NONE' else '')
702 def get_effective_coefficients(self):
703 b = num.array(
704 [v.value for v in self.numerator_coefficient_list],
705 dtype=float)
707 if self.symmetry == 'ODD':
708 b = num.concatenate((b, b[-2::-1]))
709 elif self.symmetry == 'EVEN':
710 b = num.concatenate((b, b[::-1]))
712 return b
714 def get_effective_symmetry(self):
715 if self.symmetry == 'NONE':
716 b = self.get_effective_coefficients()
717 if num.all(b - b[::-1] == 0):
718 return ['EVEN', 'ODD'][b.size % 2]
720 return self.symmetry
722 def get_ncoefficiencs(self):
723 nf = len(self.numerator_coefficient_list)
724 if self.symmetry == 'ODD':
725 nc = nf * 2 + 1
726 elif self.symmetry == 'EVEN':
727 nc = nf * 2
728 else:
729 nc = nf
731 return nc
733 def estimate_delay(self, deltat):
734 nc = self.get_ncoefficiencs()
735 if nc > 0:
736 return deltat * (nc - 1) / 2.0
737 else:
738 return 0.0
740 def get_pyrocko_response(
741 self, context, deltat, delay_responses, normalization_frequency):
743 context += self.summary()
745 if not self.numerator_coefficient_list:
746 return Delivery([])
748 b = self.get_effective_coefficients()
750 log = []
751 drop_phase = self.get_effective_symmetry() != 'NONE'
753 if not deltat:
754 log.append((
755 'error',
756 'Digital response requires knowledge about sampling '
757 'interval. Response will be unusable.',
758 context))
760 resp = response.DigitalFilterResponse(
761 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase)
763 if normalization_frequency is not None and deltat is not None:
764 normalization_frequency = 0.0
765 normalization = num.abs(evaluate1(resp, normalization_frequency))
767 if num.abs(normalization - 1.0) > 1e-2:
768 log.append((
769 'warning',
770 'FIR filter coefficients are not normalized. Normalizing '
771 'them. Factor: %g' % normalization, context))
773 resp = response.DigitalFilterResponse(
774 (b/normalization).tolist(), [1.0], deltat,
775 drop_phase=drop_phase)
777 resps = [resp]
779 if not drop_phase:
780 resps.extend(delay_responses)
782 return Delivery(resps, log=log)
785class Coefficients(BaseFilter):
786 '''
787 Response: coefficients for FIR filter. Laplace transforms or IIR filters
788 can be expressed using type as well but the PolesAndZeros should be used
789 instead. Corresponds to SEED blockette 54.
790 '''
792 cf_transfer_function_type = CfTransferFunction.T(
793 xmltagname='CfTransferFunctionType')
794 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
795 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
797 def summary(self):
798 return 'coef_%s(%i,%i%s)' % (
799 'ABC?'[
800 CfTransferFunction.choices.index(
801 self.cf_transfer_function_type)],
802 len(self.numerator_list),
803 len(self.denominator_list),
804 ',sym' if self.is_symmetric_fir else '')
806 def estimate_delay(self, deltat):
807 nc = len(self.numerator_list)
808 if nc > 0:
809 return deltat * (len(self.numerator_list) - 1) / 2.0
810 else:
811 return 0.0
813 def is_symmetric_fir(self):
814 if len(self.denominator_list) != 0:
815 return False
816 b = [v.value for v in self.numerator_list]
817 return b == b[::-1]
819 def get_pyrocko_response(
820 self, context, deltat, delay_responses, normalization_frequency):
822 context += self.summary()
824 factor = 1.0
825 if self.cf_transfer_function_type == 'ANALOG (HERTZ)':
826 factor = 2.0*math.pi
828 if not self.numerator_list and not self.denominator_list:
829 return Delivery(payload=[])
831 b = num.array(
832 [v.value*factor for v in self.numerator_list], dtype=float)
834 a = num.array(
835 [1.0] + [v.value*factor for v in self.denominator_list],
836 dtype=float)
838 log = []
839 resps = []
840 if self.cf_transfer_function_type in [
841 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']:
842 resps.append(response.AnalogFilterResponse(b, a))
844 elif self.cf_transfer_function_type == 'DIGITAL':
845 if not deltat:
846 log.append((
847 'error',
848 'Digital response requires knowledge about sampling '
849 'interval. Response will be unusable.',
850 context))
852 drop_phase = self.is_symmetric_fir()
853 resp = response.DigitalFilterResponse(
854 b, a, deltat or 0.0, drop_phase=drop_phase)
856 if normalization_frequency is not None and deltat is not None:
857 normalization = num.abs(evaluate1(
858 resp, normalization_frequency))
860 if num.abs(normalization - 1.0) > 1e-2:
861 log.append((
862 'warning',
863 'FIR filter coefficients are not normalized. '
864 'Normalizing them. Factor: %g' % normalization,
865 context))
867 resp = response.DigitalFilterResponse(
868 (b/normalization).tolist(), [1.0], deltat,
869 drop_phase=drop_phase)
871 resps.append(resp)
873 if not drop_phase:
874 resps.extend(delay_responses)
876 else:
877 return Delivery(error=(
878 'ValueError',
879 'Unknown transfer function type: %s' % (
880 self.cf_transfer_function_type)))
882 return Delivery(payload=resps, log=log)
885class Latitude(FloatWithUnit):
886 '''
887 Type for latitude coordinate.
888 '''
890 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
891 # fixed unit
892 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
895class Longitude(FloatWithUnit):
896 '''
897 Type for longitude coordinate.
898 '''
900 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
901 # fixed unit
902 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
905class PolesZeros(BaseFilter):
906 '''
907 Response: complex poles and zeros. Corresponds to SEED blockette 53.
908 '''
910 pz_transfer_function_type = PzTransferFunction.T(
911 xmltagname='PzTransferFunctionType')
912 normalization_factor = Float.T(default=1.0,
913 xmltagname='NormalizationFactor')
914 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
915 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
916 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
918 def summary(self):
919 return 'pz_%s(%i,%i)' % (
920 'ABC?'[
921 PzTransferFunction.choices.index(
922 self.pz_transfer_function_type)],
923 len(self.pole_list),
924 len(self.zero_list))
926 def get_pyrocko_response(self, context='', deltat=None):
928 context += self.summary()
930 factor = 1.0
931 cfactor = 1.0
932 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
933 factor = 2. * math.pi
934 cfactor = (2. * math.pi)**(
935 len(self.pole_list) - len(self.zero_list))
937 log = []
938 if self.normalization_factor is None \
939 or self.normalization_factor == 0.0:
941 log.append((
942 'warning',
943 'No pole-zero normalization factor given. '
944 'Assuming a value of 1.0',
945 context))
947 nfactor = 1.0
948 else:
949 nfactor = self.normalization_factor
951 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)'
952 if not is_digital:
953 resp = response.PoleZeroResponse(
954 constant=nfactor*cfactor,
955 zeros=[z.value()*factor for z in self.zero_list],
956 poles=[p.value()*factor for p in self.pole_list])
957 else:
958 if not deltat:
959 log.append((
960 'error',
961 'Digital response requires knowledge about sampling '
962 'interval. Response will be unusable.',
963 context))
965 resp = response.DigitalPoleZeroResponse(
966 constant=nfactor*cfactor,
967 zeros=[z.value()*factor for z in self.zero_list],
968 poles=[p.value()*factor for p in self.pole_list],
969 deltat=deltat or 0.0)
971 if not self.normalization_frequency.value:
972 log.append((
973 'warning',
974 'Cannot check pole-zero normalization factor: '
975 'No normalization frequency given.',
976 context))
978 else:
979 if is_digital and not deltat:
980 log.append((
981 'warning',
982 'Cannot check computed vs reported normalization '
983 'factor without knowing the sampling interval.',
984 context))
985 else:
986 computed_normalization_factor = nfactor / abs(evaluate1(
987 resp, self.normalization_frequency.value))
989 db = 20.0 * num.log10(
990 computed_normalization_factor / nfactor)
992 if abs(db) > 0.17:
993 log.append((
994 'warning',
995 'Computed and reported normalization factors differ '
996 'by %g dB: computed: %g, reported: %g' % (
997 db,
998 computed_normalization_factor,
999 nfactor),
1000 context))
1002 return Delivery([resp], log)
1005class ResponseListElement(Object):
1006 frequency = Frequency.T(xmltagname='Frequency')
1007 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
1008 phase = Angle.T(xmltagname='Phase')
1011class Polynomial(BaseFilter):
1012 '''
1013 Response: expressed as a polynomial (allows non-linear sensors to be
1014 described). Corresponds to SEED blockette 62. Can be used to describe a
1015 stage of acquisition or a complete system.
1016 '''
1018 approximation_type = Approximation.T(default='MACLAURIN',
1019 xmltagname='ApproximationType')
1020 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
1021 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
1022 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
1023 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
1024 maximum_error = Float.T(xmltagname='MaximumError')
1025 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
1027 def summary(self):
1028 return 'poly(%i)' % len(self.coefficient_list)
1031class Decimation(Object):
1032 '''
1033 Corresponds to SEED blockette 57.
1034 '''
1036 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
1037 factor = Int.T(xmltagname='Factor')
1038 offset = Int.T(xmltagname='Offset')
1039 delay = FloatWithUnit.T(xmltagname='Delay')
1040 correction = FloatWithUnit.T(xmltagname='Correction')
1042 def summary(self):
1043 return 'deci(%i, %g -> %g, %g)' % (
1044 self.factor,
1045 self.input_sample_rate.value,
1046 self.input_sample_rate.value / self.factor,
1047 self.delay.value)
1049 def get_pyrocko_response(self):
1050 if self.delay and self.delay.value != 0.0:
1051 return Delivery([response.DelayResponse(delay=-self.delay.value)])
1052 else:
1053 return Delivery([])
1056class Operator(Object):
1057 agency_list = List.T(Unicode.T(xmltagname='Agency'))
1058 contact_list = List.T(Person.T(xmltagname='Contact'))
1059 web_site = String.T(optional=True, xmltagname='WebSite')
1062class Comment(Object):
1063 '''
1064 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
1065 and 59.
1066 '''
1068 id = Counter.T(optional=True, xmlstyle='attribute')
1069 value = Unicode.T(xmltagname='Value')
1070 begin_effective_time = DummyAwareOptionalTimestamp.T(
1071 optional=True,
1072 xmltagname='BeginEffectiveTime')
1073 end_effective_time = DummyAwareOptionalTimestamp.T(
1074 optional=True,
1075 xmltagname='EndEffectiveTime')
1076 author_list = List.T(Person.T(xmltagname='Author'))
1079class ResponseList(BaseFilter):
1080 '''
1081 Response: list of frequency, amplitude and phase values. Corresponds to
1082 SEED blockette 55.
1083 '''
1085 response_list_element_list = List.T(
1086 ResponseListElement.T(xmltagname='ResponseListElement'))
1088 def summary(self):
1089 return 'list(%i)' % len(self.response_list_element_list)
1092class Log(Object):
1093 '''
1094 Container for log entries.
1095 '''
1097 entry_list = List.T(Comment.T(xmltagname='Entry'))
1100class ResponseStage(Object):
1101 '''
1102 This complex type represents channel response and covers SEED blockettes 53
1103 to 56.
1104 '''
1106 number = Counter.T(xmlstyle='attribute')
1107 resource_id = String.T(optional=True, xmlstyle='attribute')
1108 poles_zeros_list = List.T(
1109 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
1110 coefficients_list = List.T(
1111 Coefficients.T(optional=True, xmltagname='Coefficients'))
1112 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
1113 fir = FIR.T(optional=True, xmltagname='FIR')
1114 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
1115 decimation = Decimation.T(optional=True, xmltagname='Decimation')
1116 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
1118 def summary(self):
1119 elements = []
1121 for stuff in [
1122 self.poles_zeros_list, self.coefficients_list,
1123 self.response_list, self.fir, self.polynomial,
1124 self.decimation, self.stage_gain]:
1126 if stuff:
1127 if isinstance(stuff, Object):
1128 elements.append(stuff.summary())
1129 else:
1130 elements.extend(obj.summary() for obj in stuff)
1132 return '%i: %s %s -> %s' % (
1133 self.number,
1134 ', '.join(elements),
1135 self.input_units.name.upper() if self.input_units else '?',
1136 self.output_units.name.upper() if self.output_units else '?')
1138 def get_squirrel_response_stage(self, context):
1139 from pyrocko.squirrel.model import ResponseStage
1141 delivery = Delivery()
1142 delivery_pr = self.get_pyrocko_response(context)
1143 log = delivery_pr.log
1144 delivery_pr.log = []
1145 elements = delivery.extend_without_payload(delivery_pr)
1147 delivery.payload = [ResponseStage(
1148 input_quantity=to_quantity(self.input_units, context, delivery),
1149 output_quantity=to_quantity(self.output_units, context, delivery),
1150 input_sample_rate=self.input_sample_rate,
1151 output_sample_rate=self.output_sample_rate,
1152 elements=elements,
1153 log=log)]
1155 return delivery
1157 def get_pyrocko_response(self, context, gain_only=False):
1159 context = context + ', stage %i' % self.number
1161 responses = []
1162 log = []
1163 if self.stage_gain:
1164 normalization_frequency = self.stage_gain.frequency or 0.0
1165 else:
1166 normalization_frequency = 0.0
1168 if not gain_only:
1169 deltat = None
1170 delay_responses = []
1171 if self.decimation:
1172 rate = self.decimation.input_sample_rate.value
1173 if rate > 0.0:
1174 deltat = 1.0 / rate
1175 delivery = self.decimation.get_pyrocko_response()
1176 if delivery.errors:
1177 return delivery
1179 delay_responses = delivery.payload
1180 log.extend(delivery.log)
1182 for pzs in self.poles_zeros_list:
1183 delivery = pzs.get_pyrocko_response(context, deltat)
1184 if delivery.errors:
1185 return delivery
1187 pz_resps = delivery.payload
1188 log.extend(delivery.log)
1189 responses.extend(pz_resps)
1191 # emulate incorrect? evalresp behaviour
1192 if pzs.normalization_frequency != normalization_frequency \
1193 and normalization_frequency != 0.0:
1195 try:
1196 trial = response.MultiplyResponse(pz_resps)
1197 anorm = num.abs(evaluate1(
1198 trial, pzs.normalization_frequency.value))
1199 asens = num.abs(
1200 evaluate1(trial, normalization_frequency))
1202 factor = anorm/asens
1204 if abs(factor - 1.0) > 0.01:
1205 log.append((
1206 'warning',
1207 'PZ normalization frequency (%g) is different '
1208 'from stage gain frequency (%s) -> Emulating '
1209 'possibly incorrect evalresp behaviour. '
1210 'Correction factor: %g' % (
1211 pzs.normalization_frequency.value,
1212 normalization_frequency,
1213 factor),
1214 context))
1216 responses.append(
1217 response.PoleZeroResponse(constant=factor))
1218 except response.InvalidResponseError as e:
1219 log.append((
1220 'warning',
1221 'Could not check response: %s' % str(e),
1222 context))
1224 if len(self.poles_zeros_list) > 1:
1225 log.append((
1226 'warning',
1227 'Multiple poles and zeros records in single response '
1228 'stage.',
1229 context))
1231 for cfs in self.coefficients_list + (
1232 [self.fir] if self.fir else []):
1234 delivery = cfs.get_pyrocko_response(
1235 context, deltat, delay_responses,
1236 normalization_frequency)
1238 if delivery.errors:
1239 return delivery
1241 responses.extend(delivery.payload)
1242 log.extend(delivery.log)
1244 if len(self.coefficients_list) > 1:
1245 log.append((
1246 'warning',
1247 'Multiple filter coefficients lists in single response '
1248 'stage.',
1249 context))
1251 if self.response_list:
1252 log.append((
1253 'warning',
1254 'Unhandled response element of type: ResponseList',
1255 context))
1257 if self.polynomial:
1258 log.append((
1259 'warning',
1260 'Unhandled response element of type: Polynomial',
1261 context))
1263 if self.stage_gain:
1264 responses.append(
1265 response.PoleZeroResponse(constant=self.stage_gain.value))
1267 return Delivery(responses, log)
1269 @property
1270 def input_units(self):
1271 for e in (self.poles_zeros_list + self.coefficients_list +
1272 [self.response_list, self.fir, self.polynomial]):
1273 if e is not None:
1274 return e.input_units
1276 return None
1278 @property
1279 def output_units(self):
1280 for e in (self.poles_zeros_list + self.coefficients_list +
1281 [self.response_list, self.fir, self.polynomial]):
1282 if e is not None:
1283 return e.output_units
1285 return None
1287 @property
1288 def input_sample_rate(self):
1289 if self.decimation:
1290 return self.decimation.input_sample_rate.value
1292 return None
1294 @property
1295 def output_sample_rate(self):
1296 if self.decimation:
1297 return self.decimation.input_sample_rate.value \
1298 / self.decimation.factor
1300 return None
1303class Response(Object):
1304 resource_id = String.T(optional=True, xmlstyle='attribute')
1305 instrument_sensitivity = Sensitivity.T(optional=True,
1306 xmltagname='InstrumentSensitivity')
1307 instrument_polynomial = Polynomial.T(optional=True,
1308 xmltagname='InstrumentPolynomial')
1309 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
1311 def check_sample_rates(self, channel):
1313 if self.stage_list:
1314 sample_rate = None
1316 for stage in self.stage_list:
1317 if stage.decimation:
1318 input_sample_rate = \
1319 stage.decimation.input_sample_rate.value
1321 if sample_rate is not None and not same_sample_rate(
1322 sample_rate, input_sample_rate):
1324 logger.warning(
1325 'Response stage %i has unexpected input sample '
1326 'rate: %g Hz (expected: %g Hz)' % (
1327 stage.number,
1328 input_sample_rate,
1329 sample_rate))
1331 sample_rate = input_sample_rate / stage.decimation.factor
1333 if sample_rate is not None and channel.sample_rate \
1334 and not same_sample_rate(
1335 sample_rate, channel.sample_rate.value):
1337 logger.warning(
1338 'Channel sample rate (%g Hz) does not match sample rate '
1339 'deducted from response stages information (%g Hz).' % (
1340 channel.sample_rate.value,
1341 sample_rate))
1343 def check_units(self):
1345 if self.instrument_sensitivity \
1346 and self.instrument_sensitivity.input_units:
1348 units = self.instrument_sensitivity.input_units.name.upper()
1350 if self.stage_list:
1351 for stage in self.stage_list:
1352 if units and stage.input_units \
1353 and stage.input_units.name.upper() != units:
1355 logger.warning(
1356 'Input units of stage %i (%s) do not match %s (%s).'
1357 % (
1358 stage.number,
1359 units,
1360 'output units of stage %i'
1361 if stage.number == 0
1362 else 'sensitivity input units',
1363 units))
1365 if stage.output_units:
1366 units = stage.output_units.name.upper()
1367 else:
1368 units = None
1370 sout_units = self.instrument_sensitivity.output_units
1371 if self.instrument_sensitivity and sout_units:
1372 if units is not None and units != sout_units.name.upper():
1373 logger.warning(
1374 'Output units of stage %i (%s) do not match %s (%s).'
1375 % (
1376 stage.number,
1377 units,
1378 'sensitivity output units',
1379 sout_units.name.upper()))
1381 def _sensitivity_checkpoints(self, responses, context):
1382 delivery = Delivery()
1384 if self.instrument_sensitivity:
1385 sval = self.instrument_sensitivity.value
1386 sfreq = self.instrument_sensitivity.frequency
1387 if sval is None:
1388 delivery.log.append((
1389 'warning',
1390 'No sensitivity value given.',
1391 context))
1393 elif sval is None:
1394 delivery.log.append((
1395 'warning',
1396 'Reported sensitivity value is zero.',
1397 context))
1399 elif sfreq is None:
1400 delivery.log.append((
1401 'warning',
1402 'Sensitivity frequency not given.',
1403 context))
1405 else:
1406 trial = response.MultiplyResponse(responses)
1408 delivery.extend(
1409 check_resp(
1410 trial, sval, sfreq, 0.1,
1411 'Instrument sensitivity value inconsistent with '
1412 'sensitivity computed from complete response',
1413 context))
1415 delivery.payload.append(response.FrequencyResponseCheckpoint(
1416 frequency=sfreq, value=sval))
1418 return delivery
1420 def get_squirrel_response(self, context, **kwargs):
1421 from pyrocko.squirrel.model import Response
1423 if self.stage_list:
1424 delivery = Delivery()
1425 for istage, stage in enumerate(self.stage_list):
1426 delivery.extend(stage.get_squirrel_response_stage(context))
1428 checkpoints = []
1429 if not delivery.errors:
1430 all_responses = []
1431 for sq_stage in delivery.payload:
1432 all_responses.extend(sq_stage.elements)
1434 checkpoints.extend(delivery.extend_without_payload(
1435 self._sensitivity_checkpoints(all_responses, context)))
1437 sq_stages = delivery.payload
1438 if sq_stages:
1439 if sq_stages[0].input_quantity is None \
1440 and self.instrument_sensitivity is not None:
1442 sq_stages[0].input_quantity = to_quantity(
1443 self.instrument_sensitivity.input_units,
1444 context, delivery)
1445 sq_stages[-1].output_quantity = to_quantity(
1446 self.instrument_sensitivity.output_units,
1447 context, delivery)
1449 sq_stages = delivery.expect()
1451 return Response(
1452 stages=sq_stages,
1453 log=delivery.log,
1454 checkpoints=checkpoints,
1455 **kwargs)
1457 elif self.instrument_sensitivity:
1458 raise NoResponseInformation(
1459 'Only instrument sensitivity given (won\'t use it). (%s).'
1460 % context)
1461 else:
1462 raise NoResponseInformation(
1463 'Empty instrument response (%s).'
1464 % context)
1466 def get_pyrocko_response(
1467 self, context, fake_input_units=None, stages=(0, 1)):
1469 delivery = Delivery()
1470 if self.stage_list:
1471 for istage, stage in enumerate(self.stage_list):
1472 delivery.extend(stage.get_pyrocko_response(
1473 context, gain_only=not (
1474 stages is None or stages[0] <= istage < stages[1])))
1476 elif self.instrument_sensitivity:
1477 delivery.extend(self.instrument_sensitivity.get_pyrocko_response())
1479 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context)
1480 checkpoints = delivery.extend_without_payload(delivery_cp)
1481 if delivery.errors:
1482 return delivery
1484 if fake_input_units is not None:
1485 if not self.instrument_sensitivity or \
1486 self.instrument_sensitivity.input_units is None:
1488 delivery.errors.append((
1489 'NoResponseInformation',
1490 'No input units given, so cannot convert to requested '
1491 'units: %s' % fake_input_units.upper(),
1492 context))
1494 return delivery
1496 input_units = self.instrument_sensitivity.input_units.name.upper()
1498 conresp = None
1499 try:
1500 conresp = conversion[
1501 fake_input_units.upper(), input_units]
1503 except KeyError:
1504 delivery.errors.append((
1505 'NoResponseInformation',
1506 'Cannot convert between units: %s, %s'
1507 % (fake_input_units, input_units),
1508 context))
1510 if conresp is not None:
1511 delivery.payload.append(conresp)
1512 for checkpoint in checkpoints:
1513 checkpoint.value *= num.abs(evaluate1(
1514 conresp, checkpoint.frequency))
1516 delivery.payload = [
1517 response.MultiplyResponse(
1518 delivery.payload,
1519 checkpoints=checkpoints)]
1521 return delivery
1523 @classmethod
1524 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
1525 normalization_frequency=1.0):
1526 '''
1527 Convert Pyrocko pole-zero response to StationXML response.
1529 :param presponse: Pyrocko pole-zero response
1530 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse`
1531 :param input_unit: Input unit to be reported in the StationXML
1532 response.
1533 :type input_unit: str
1534 :param output_unit: Output unit to be reported in the StationXML
1535 response.
1536 :type output_unit: str
1537 :param normalization_frequency: Frequency where the normalization
1538 factor for the StationXML response should be computed.
1539 :type normalization_frequency: float
1540 '''
1542 norm_factor = 1.0/float(abs(
1543 evaluate1(presponse, normalization_frequency)
1544 / presponse.constant))
1546 pzs = PolesZeros(
1547 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
1548 normalization_factor=norm_factor,
1549 normalization_frequency=Frequency(normalization_frequency),
1550 zero_list=[PoleZero(real=FloatNoUnit(z.real),
1551 imaginary=FloatNoUnit(z.imag))
1552 for z in presponse.zeros],
1553 pole_list=[PoleZero(real=FloatNoUnit(z.real),
1554 imaginary=FloatNoUnit(z.imag))
1555 for z in presponse.poles])
1557 pzs.validate()
1559 stage = ResponseStage(
1560 number=1,
1561 poles_zeros_list=[pzs],
1562 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
1564 resp = Response(
1565 instrument_sensitivity=Sensitivity(
1566 value=stage.stage_gain.value,
1567 frequency=normalization_frequency,
1568 input_units=Units(input_unit),
1569 output_units=Units(output_unit)),
1571 stage_list=[stage])
1573 return resp
1576class BaseNode(Object):
1577 '''
1578 A base node type for derivation from: Network, Station and Channel types.
1579 '''
1581 code = String.T(xmlstyle='attribute')
1582 start_date = DummyAwareOptionalTimestamp.T(optional=True,
1583 xmlstyle='attribute')
1584 end_date = DummyAwareOptionalTimestamp.T(optional=True,
1585 xmlstyle='attribute')
1586 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
1587 alternate_code = String.T(optional=True, xmlstyle='attribute')
1588 historical_code = String.T(optional=True, xmlstyle='attribute')
1589 description = Unicode.T(optional=True, xmltagname='Description')
1590 comment_list = List.T(Comment.T(xmltagname='Comment'))
1592 def spans(self, *args):
1593 if len(args) == 0:
1594 return True
1595 elif len(args) == 1:
1596 return ((self.start_date is None or
1597 self.start_date <= args[0]) and
1598 (self.end_date is None or
1599 args[0] <= self.end_date))
1601 elif len(args) == 2:
1602 return ((self.start_date is None or
1603 args[1] >= self.start_date) and
1604 (self.end_date is None or
1605 self.end_date >= args[0]))
1608def overlaps(a, b):
1609 return (
1610 a.start_date is None or b.end_date is None
1611 or a.start_date < b.end_date
1612 ) and (
1613 b.start_date is None or a.end_date is None
1614 or b.start_date < a.end_date
1615 )
1618class Channel(BaseNode):
1619 '''
1620 Equivalent to SEED blockette 52 and parent element for the related the
1621 response blockettes.
1622 '''
1624 location_code = String.T(xmlstyle='attribute')
1625 external_reference_list = List.T(
1626 ExternalReference.T(xmltagname='ExternalReference'))
1627 latitude = Latitude.T(xmltagname='Latitude')
1628 longitude = Longitude.T(xmltagname='Longitude')
1629 elevation = Distance.T(xmltagname='Elevation')
1630 depth = Distance.T(xmltagname='Depth')
1631 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
1632 dip = Dip.T(optional=True, xmltagname='Dip')
1633 type_list = List.T(Type.T(xmltagname='Type'))
1634 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
1635 sample_rate_ratio = SampleRateRatio.T(optional=True,
1636 xmltagname='SampleRateRatio')
1637 storage_format = String.T(optional=True, xmltagname='StorageFormat')
1638 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
1639 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
1640 sensor = Equipment.T(optional=True, xmltagname='Sensor')
1641 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
1642 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
1643 equipment = Equipment.T(optional=True, xmltagname='Equipment')
1644 response = Response.T(optional=True, xmltagname='Response')
1646 @property
1647 def position_values(self):
1648 lat = self.latitude.value
1649 lon = self.longitude.value
1650 elevation = value_or_none(self.elevation)
1651 depth = value_or_none(self.depth)
1652 return lat, lon, elevation, depth
1655class Station(BaseNode):
1656 '''
1657 This type represents a Station epoch. It is common to only have a single
1658 station epoch with the station's creation and termination dates as the
1659 epoch start and end dates.
1660 '''
1662 latitude = Latitude.T(xmltagname='Latitude')
1663 longitude = Longitude.T(xmltagname='Longitude')
1664 elevation = Distance.T(xmltagname='Elevation')
1665 site = Site.T(optional=True, xmltagname='Site')
1666 vault = Unicode.T(optional=True, xmltagname='Vault')
1667 geology = Unicode.T(optional=True, xmltagname='Geology')
1668 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
1669 operator_list = List.T(Operator.T(xmltagname='Operator'))
1670 creation_date = DummyAwareOptionalTimestamp.T(
1671 optional=True, xmltagname='CreationDate')
1672 termination_date = DummyAwareOptionalTimestamp.T(
1673 optional=True, xmltagname='TerminationDate')
1674 total_number_channels = Counter.T(
1675 optional=True, xmltagname='TotalNumberChannels')
1676 selected_number_channels = Counter.T(
1677 optional=True, xmltagname='SelectedNumberChannels')
1678 external_reference_list = List.T(
1679 ExternalReference.T(xmltagname='ExternalReference'))
1680 channel_list = List.T(Channel.T(xmltagname='Channel'))
1682 @property
1683 def position_values(self):
1684 lat = self.latitude.value
1685 lon = self.longitude.value
1686 elevation = value_or_none(self.elevation)
1687 return lat, lon, elevation
1690class Network(BaseNode):
1691 '''
1692 This type represents the Network layer, all station metadata is contained
1693 within this element. The official name of the network or other descriptive
1694 information can be included in the Description element. The Network can
1695 contain 0 or more Stations.
1696 '''
1698 total_number_stations = Counter.T(optional=True,
1699 xmltagname='TotalNumberStations')
1700 selected_number_stations = Counter.T(optional=True,
1701 xmltagname='SelectedNumberStations')
1702 station_list = List.T(Station.T(xmltagname='Station'))
1704 @property
1705 def station_code_list(self):
1706 return sorted(set(s.code for s in self.station_list))
1708 @property
1709 def sl_code_list(self):
1710 sls = set()
1711 for station in self.station_list:
1712 for channel in station.channel_list:
1713 sls.add((station.code, channel.location_code))
1715 return sorted(sls)
1717 def summary(self, width=80, indent=4):
1718 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1719 lines = ['%s (%i):' % (self.code, len(sls))]
1720 if sls:
1721 ssls = ['.'.join(x for x in c if x) for c in sls]
1722 w = max(len(x) for x in ssls)
1723 n = (width - indent) / (w+1)
1724 while ssls:
1725 lines.append(
1726 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1728 ssls[:n] = []
1730 return '\n'.join(lines)
1733def value_or_none(x):
1734 if x is not None:
1735 return x.value
1736 else:
1737 return None
1740def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1742 pos = lat, lon, elevation, depth = \
1743 channels[0].position_values
1745 if not all(pos == x.position_values for x in channels):
1746 info = '\n'.join(
1747 ' %s: %s' % (x.code, x.position_values) for
1748 x in channels)
1750 mess = 'encountered inconsistencies in channel ' \
1751 'lat/lon/elevation/depth ' \
1752 'for %s.%s.%s: \n%s' % (nsl + (info,))
1754 if inconsistencies == 'raise':
1755 raise InconsistentChannelLocations(mess)
1757 elif inconsistencies == 'warn':
1758 logger.warning(mess)
1759 logger.warning(' -> using mean values')
1761 apos = num.array([x.position_values for x in channels], dtype=float)
1762 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1763 / num.sum(num.isfinite(apos), axis=0)
1765 groups = {}
1766 for channel in channels:
1767 if channel.code not in groups:
1768 groups[channel.code] = []
1770 groups[channel.code].append(channel)
1772 pchannels = []
1773 for code in sorted(groups.keys()):
1774 data = [
1775 (channel.code, value_or_none(channel.azimuth),
1776 value_or_none(channel.dip))
1777 for channel in groups[code]]
1779 azimuth, dip = util.consistency_merge(
1780 data,
1781 message='channel orientation values differ:',
1782 error=inconsistencies)
1784 pchannels.append(
1785 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1787 return pyrocko.model.Station(
1788 *nsl,
1789 lat=mlat,
1790 lon=mlon,
1791 elevation=mele,
1792 depth=mdep,
1793 channels=pchannels)
1796class FDSNStationXML(Object):
1797 '''
1798 Top-level type for Station XML. Required field are Source (network ID of
1799 the institution sending the message) and one or more Network containers or
1800 one or more Station containers.
1801 '''
1803 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1804 source = String.T(xmltagname='Source')
1805 sender = String.T(optional=True, xmltagname='Sender')
1806 module = String.T(optional=True, xmltagname='Module')
1807 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1808 created = Timestamp.T(optional=True, xmltagname='Created')
1809 network_list = List.T(Network.T(xmltagname='Network'))
1811 xmltagname = 'FDSNStationXML'
1812 guessable_xmlns = [guts_xmlns]
1814 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1815 time=None, timespan=None,
1816 inconsistencies='warn'):
1818 assert inconsistencies in ('raise', 'warn')
1820 if nslcs is not None:
1821 nslcs = set(nslcs)
1823 if nsls is not None:
1824 nsls = set(nsls)
1826 tt = ()
1827 if time is not None:
1828 tt = (time,)
1829 elif timespan is not None:
1830 tt = timespan
1832 pstations = []
1833 for network in self.network_list:
1834 if not network.spans(*tt):
1835 continue
1837 for station in network.station_list:
1838 if not station.spans(*tt):
1839 continue
1841 if station.channel_list:
1842 loc_to_channels = {}
1843 for channel in station.channel_list:
1844 if not channel.spans(*tt):
1845 continue
1847 loc = channel.location_code.strip()
1848 if loc not in loc_to_channels:
1849 loc_to_channels[loc] = []
1851 loc_to_channels[loc].append(channel)
1853 for loc in sorted(loc_to_channels.keys()):
1854 channels = loc_to_channels[loc]
1855 if nslcs is not None:
1856 channels = [channel for channel in channels
1857 if (network.code, station.code, loc,
1858 channel.code) in nslcs]
1860 if not channels:
1861 continue
1863 nsl = network.code, station.code, loc
1864 if nsls is not None and nsl not in nsls:
1865 continue
1867 pstations.append(
1868 pyrocko_station_from_channels(
1869 nsl,
1870 channels,
1871 inconsistencies=inconsistencies))
1872 else:
1873 pstations.append(pyrocko.model.Station(
1874 network.code, station.code, '*',
1875 lat=station.latitude.value,
1876 lon=station.longitude.value,
1877 elevation=value_or_none(station.elevation),
1878 name=station.description or ''))
1880 return pstations
1882 @classmethod
1883 def from_pyrocko_stations(
1884 cls, pyrocko_stations, add_flat_responses_from=None):
1886 '''
1887 Generate :py:class:`FDSNStationXML` from list of
1888 :py:class;`pyrocko.model.Station` instances.
1890 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1891 instances.
1892 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1893 '''
1894 from collections import defaultdict
1895 network_dict = defaultdict(list)
1897 if add_flat_responses_from:
1898 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1899 extra = dict(
1900 response=Response(
1901 instrument_sensitivity=Sensitivity(
1902 value=1.0,
1903 frequency=1.0,
1904 input_units=Units(name=add_flat_responses_from))))
1905 else:
1906 extra = {}
1908 have_offsets = set()
1909 for s in pyrocko_stations:
1911 if s.north_shift != 0.0 or s.east_shift != 0.0:
1912 have_offsets.add(s.nsl())
1914 network, station, location = s.nsl()
1915 channel_list = []
1916 for c in s.channels:
1917 channel_list.append(
1918 Channel(
1919 location_code=location,
1920 code=c.name,
1921 latitude=Latitude(value=s.effective_lat),
1922 longitude=Longitude(value=s.effective_lon),
1923 elevation=Distance(value=s.elevation),
1924 depth=Distance(value=s.depth),
1925 azimuth=Azimuth(value=c.azimuth),
1926 dip=Dip(value=c.dip),
1927 **extra
1928 )
1929 )
1931 network_dict[network].append(
1932 Station(
1933 code=station,
1934 latitude=Latitude(value=s.effective_lat),
1935 longitude=Longitude(value=s.effective_lon),
1936 elevation=Distance(value=s.elevation),
1937 channel_list=channel_list)
1938 )
1940 if have_offsets:
1941 logger.warning(
1942 'StationXML does not support Cartesian offsets in '
1943 'coordinates. Storing effective lat/lon for stations: %s' %
1944 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1946 timestamp = util.to_time_float(time.time())
1947 network_list = []
1948 for k, station_list in network_dict.items():
1950 network_list.append(
1951 Network(
1952 code=k, station_list=station_list,
1953 total_number_stations=len(station_list)))
1955 sxml = FDSNStationXML(
1956 source='from pyrocko stations list', created=timestamp,
1957 network_list=network_list)
1959 sxml.validate()
1960 return sxml
1962 def iter_network_stations(
1963 self, net=None, sta=None, time=None, timespan=None):
1965 tt = ()
1966 if time is not None:
1967 tt = (time,)
1968 elif timespan is not None:
1969 tt = timespan
1971 for network in self.network_list:
1972 if not network.spans(*tt) or (
1973 net is not None and network.code != net):
1974 continue
1976 for station in network.station_list:
1977 if not station.spans(*tt) or (
1978 sta is not None and station.code != sta):
1979 continue
1981 yield (network, station)
1983 def iter_network_station_channels(
1984 self, net=None, sta=None, loc=None, cha=None,
1985 time=None, timespan=None):
1987 if loc is not None:
1988 loc = loc.strip()
1990 tt = ()
1991 if time is not None:
1992 tt = (time,)
1993 elif timespan is not None:
1994 tt = timespan
1996 for network in self.network_list:
1997 if not network.spans(*tt) or (
1998 net is not None and network.code != net):
1999 continue
2001 for station in network.station_list:
2002 if not station.spans(*tt) or (
2003 sta is not None and station.code != sta):
2004 continue
2006 if station.channel_list:
2007 for channel in station.channel_list:
2008 if (not channel.spans(*tt) or
2009 (cha is not None and channel.code != cha) or
2010 (loc is not None and
2011 channel.location_code.strip() != loc)):
2012 continue
2014 yield (network, station, channel)
2016 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
2017 time=None, timespan=None):
2019 groups = {}
2020 for network, station, channel in self.iter_network_station_channels(
2021 net, sta, loc, cha, time=time, timespan=timespan):
2023 net = network.code
2024 sta = station.code
2025 cha = channel.code
2026 loc = channel.location_code.strip()
2027 if len(cha) == 3:
2028 bic = cha[:2] # band and intrument code according to SEED
2029 elif len(cha) == 1:
2030 bic = ''
2031 else:
2032 bic = cha
2034 if channel.response and \
2035 channel.response.instrument_sensitivity and \
2036 channel.response.instrument_sensitivity.input_units:
2038 unit = channel.response.instrument_sensitivity\
2039 .input_units.name.upper()
2040 else:
2041 unit = None
2043 bic = (bic, unit)
2045 k = net, sta, loc
2046 if k not in groups:
2047 groups[k] = {}
2049 if bic not in groups[k]:
2050 groups[k][bic] = []
2052 groups[k][bic].append(channel)
2054 for nsl, bic_to_channels in groups.items():
2055 bad_bics = []
2056 for bic, channels in bic_to_channels.items():
2057 sample_rates = []
2058 for channel in channels:
2059 sample_rates.append(channel.sample_rate.value)
2061 if not same(sample_rates):
2062 scs = ','.join(channel.code for channel in channels)
2063 srs = ', '.join('%e' % x for x in sample_rates)
2064 err = 'ignoring channels with inconsistent sampling ' + \
2065 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
2067 logger.warning(err)
2068 bad_bics.append(bic)
2070 for bic in bad_bics:
2071 del bic_to_channels[bic]
2073 return groups
2075 def choose_channels(
2076 self,
2077 target_sample_rate=None,
2078 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
2079 priority_units=['M/S', 'M/S**2'],
2080 priority_instrument_code=['H', 'L'],
2081 time=None,
2082 timespan=None):
2084 nslcs = {}
2085 for nsl, bic_to_channels in self.get_channel_groups(
2086 time=time, timespan=timespan).items():
2088 useful_bics = []
2089 for bic, channels in bic_to_channels.items():
2090 rate = channels[0].sample_rate.value
2092 if target_sample_rate is not None and \
2093 rate < target_sample_rate*0.99999:
2094 continue
2096 if len(bic[0]) == 2:
2097 if bic[0][0] not in priority_band_code:
2098 continue
2100 if bic[0][1] not in priority_instrument_code:
2101 continue
2103 unit = bic[1]
2105 prio_unit = len(priority_units)
2106 try:
2107 prio_unit = priority_units.index(unit)
2108 except ValueError:
2109 pass
2111 prio_inst = len(priority_instrument_code)
2112 prio_band = len(priority_band_code)
2113 if len(channels[0].code) == 3:
2114 try:
2115 prio_inst = priority_instrument_code.index(
2116 channels[0].code[1])
2117 except ValueError:
2118 pass
2120 try:
2121 prio_band = priority_band_code.index(
2122 channels[0].code[0])
2123 except ValueError:
2124 pass
2126 if target_sample_rate is None:
2127 rate = -rate
2129 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
2130 prio_inst, bic))
2132 useful_bics.sort()
2134 for _, _, rate, _, _, bic in useful_bics:
2135 channels = sorted(
2136 bic_to_channels[bic],
2137 key=lambda channel: channel.code)
2139 if channels:
2140 for channel in channels:
2141 nslcs[nsl + (channel.code,)] = channel
2143 break
2145 return nslcs
2147 def get_pyrocko_response(
2148 self, nslc,
2149 time=None, timespan=None, fake_input_units=None, stages=(0, 1)):
2151 net, sta, loc, cha = nslc
2152 resps = []
2153 for _, _, channel in self.iter_network_station_channels(
2154 net, sta, loc, cha, time=time, timespan=timespan):
2155 resp = channel.response
2156 if resp:
2157 resp.check_sample_rates(channel)
2158 resp.check_units()
2159 resps.append(resp.get_pyrocko_response(
2160 '.'.join(nslc),
2161 fake_input_units=fake_input_units,
2162 stages=stages).expect_one())
2164 if not resps:
2165 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
2166 elif len(resps) > 1:
2167 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
2169 return resps[0]
2171 @property
2172 def n_code_list(self):
2173 return sorted(set(x.code for x in self.network_list))
2175 @property
2176 def ns_code_list(self):
2177 nss = set()
2178 for network in self.network_list:
2179 for station in network.station_list:
2180 nss.add((network.code, station.code))
2182 return sorted(nss)
2184 @property
2185 def nsl_code_list(self):
2186 nsls = set()
2187 for network in self.network_list:
2188 for station in network.station_list:
2189 for channel in station.channel_list:
2190 nsls.add(
2191 (network.code, station.code, channel.location_code))
2193 return sorted(nsls)
2195 @property
2196 def nslc_code_list(self):
2197 nslcs = set()
2198 for network in self.network_list:
2199 for station in network.station_list:
2200 for channel in station.channel_list:
2201 nslcs.add(
2202 (network.code, station.code, channel.location_code,
2203 channel.code))
2205 return sorted(nslcs)
2207 def summary(self):
2208 lst = [
2209 'number of n codes: %i' % len(self.n_code_list),
2210 'number of ns codes: %i' % len(self.ns_code_list),
2211 'number of nsl codes: %i' % len(self.nsl_code_list),
2212 'number of nslc codes: %i' % len(self.nslc_code_list)
2213 ]
2214 return '\n'.join(lst)
2216 def summary_stages(self):
2217 data = []
2218 for network, station, channel in self.iter_network_station_channels():
2219 nslc = (network.code, station.code, channel.location_code,
2220 channel.code)
2222 stages = []
2223 in_units = '?'
2224 out_units = '?'
2225 if channel.response:
2226 sens = channel.response.instrument_sensitivity
2227 if sens:
2228 in_units = sens.input_units.name.upper()
2229 out_units = sens.output_units.name.upper()
2231 for stage in channel.response.stage_list:
2232 stages.append(stage.summary())
2234 data.append(
2235 (nslc, tts(channel.start_date), tts(channel.end_date),
2236 in_units, out_units, stages))
2238 data.sort()
2240 lst = []
2241 for nslc, stmin, stmax, in_units, out_units, stages in data:
2242 lst.append(' %s: %s - %s, %s -> %s' % (
2243 '.'.join(nslc), stmin, stmax, in_units, out_units))
2244 for stage in stages:
2245 lst.append(' %s' % stage)
2247 return '\n'.join(lst)
2249 def _check_overlaps(self):
2250 by_nslc = {}
2251 for network in self.network_list:
2252 for station in network.station_list:
2253 for channel in station.channel_list:
2254 nslc = (network.code, station.code, channel.location_code,
2255 channel.code)
2256 if nslc not in by_nslc:
2257 by_nslc[nslc] = []
2259 by_nslc[nslc].append(channel)
2261 errors = []
2262 for nslc, channels in by_nslc.items():
2263 for ia, a in enumerate(channels):
2264 for b in channels[ia+1:]:
2265 if overlaps(a, b):
2266 errors.append(
2267 'Channel epochs overlap for %s:\n'
2268 ' %s - %s\n %s - %s' % (
2269 '.'.join(nslc),
2270 tts(a.start_date), tts(a.end_date),
2271 tts(b.start_date), tts(b.end_date)))
2272 return errors
2274 def check(self):
2275 errors = []
2276 for meth in [self._check_overlaps]:
2277 errors.extend(meth())
2279 if errors:
2280 raise Inconsistencies(
2281 'Inconsistencies found in StationXML:\n '
2282 + '\n '.join(errors))
2285def load_channel_table(stream):
2287 networks = {}
2288 stations = {}
2290 for line in stream:
2291 line = str(line.decode('ascii'))
2292 if line.startswith('#'):
2293 continue
2295 t = line.rstrip().split('|')
2297 if len(t) != 17:
2298 logger.warning('Invalid channel record: %s' % line)
2299 continue
2301 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
2302 scale_freq, scale_units, sample_rate, start_date, end_date) = t
2304 try:
2305 scale = float(scale)
2306 except ValueError:
2307 scale = None
2309 try:
2310 scale_freq = float(scale_freq)
2311 except ValueError:
2312 scale_freq = None
2314 try:
2315 depth = float(dep)
2316 except ValueError:
2317 depth = 0.0
2319 try:
2320 azi = float(azi)
2321 dip = float(dip)
2322 except ValueError:
2323 azi = None
2324 dip = None
2326 try:
2327 if net not in networks:
2328 network = Network(code=net)
2329 else:
2330 network = networks[net]
2332 if (net, sta) not in stations:
2333 station = Station(
2334 code=sta, latitude=lat, longitude=lon, elevation=ele)
2336 station.regularize()
2337 else:
2338 station = stations[net, sta]
2340 if scale:
2341 resp = Response(
2342 instrument_sensitivity=Sensitivity(
2343 value=scale,
2344 frequency=scale_freq,
2345 input_units=scale_units))
2346 else:
2347 resp = None
2349 channel = Channel(
2350 code=cha,
2351 location_code=loc.strip(),
2352 latitude=lat,
2353 longitude=lon,
2354 elevation=ele,
2355 depth=depth,
2356 azimuth=azi,
2357 dip=dip,
2358 sensor=Equipment(description=sens),
2359 response=resp,
2360 sample_rate=sample_rate,
2361 start_date=start_date,
2362 end_date=end_date or None)
2364 channel.regularize()
2366 except ValidationError:
2367 raise InvalidRecord(line)
2369 if net not in networks:
2370 networks[net] = network
2372 if (net, sta) not in stations:
2373 stations[net, sta] = station
2374 network.station_list.append(station)
2376 station.channel_list.append(channel)
2378 return FDSNStationXML(
2379 source='created from table input',
2380 created=time.time(),
2381 network_list=sorted(networks.values(), key=lambda x: x.code))
2384def primitive_merge(sxs):
2385 networks = []
2386 for sx in sxs:
2387 networks.extend(sx.network_list)
2389 return FDSNStationXML(
2390 source='merged from different sources',
2391 created=time.time(),
2392 network_list=copy.deepcopy(
2393 sorted(networks, key=lambda x: x.code)))
2396def split_channels(sx):
2397 for nslc in sx.nslc_code_list:
2398 network_list = sx.network_list
2399 network_list_filtered = [
2400 network for network in network_list
2401 if network.code == nslc[0]]
2403 for network in network_list_filtered:
2404 sx.network_list = [network]
2405 station_list = network.station_list
2406 station_list_filtered = [
2407 station for station in station_list
2408 if station.code == nslc[1]]
2410 for station in station_list_filtered:
2411 network.station_list = [station]
2412 channel_list = station.channel_list
2413 station.channel_list = [
2414 channel for channel in channel_list
2415 if (channel.location_code, channel.code) == nslc[2:4]]
2417 yield nslc, copy.deepcopy(sx)
2419 station.channel_list = channel_list
2421 network.station_list = station_list
2423 sx.network_list = network_list
2426if __name__ == '__main__':
2427 from optparse import OptionParser
2429 util.setup_logging('pyrocko.io.stationxml', 'warning')
2431 usage = \
2432 'python -m pyrocko.io.stationxml check|stats|stages ' \
2433 '<filename> [options]'
2435 description = '''Torture StationXML file.'''
2437 parser = OptionParser(
2438 usage=usage,
2439 description=description,
2440 formatter=util.BetterHelpFormatter())
2442 (options, args) = parser.parse_args(sys.argv[1:])
2444 if len(args) != 2:
2445 parser.print_help()
2446 sys.exit(1)
2448 action, path = args
2450 sx = load_xml(filename=path)
2451 if action == 'check':
2452 try:
2453 sx.check()
2454 except Inconsistencies as e:
2455 logger.error(e)
2456 sys.exit(1)
2458 elif action == 'stats':
2459 print(sx.summary())
2461 elif action == 'stages':
2462 print(sx.summary_stages())
2464 else:
2465 parser.print_help()
2466 sys.exit('unknown action: %s' % action)