1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
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)
20from pyrocko.guts import load_xml # noqa
21from pyrocko.util import hpfloat, time_to_str, get_time_float
23import pyrocko.model
24from pyrocko import trace, util
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'): trace.IntegrationResponse(1),
40 ('M/S**2', 'M'): trace.IntegrationResponse(2),
41 ('M', 'M/S'): trace.DifferentiationResponse(1),
42 ('M/S', 'M/S'): None,
43 ('M/S**2', 'M/S'): trace.IntegrationResponse(1),
44 ('M', 'M/S**2'): trace.DifferentiationResponse(2),
45 ('M/S', 'M/S**2'): trace.DifferentiationResponse(1),
46 ('M/S**2', 'M/S**2'): None}
49class NoResponseInformation(Exception):
50 pass
53class MultipleResponseInformation(Exception):
54 pass
57def wrap(s, width=80, indent=4):
58 words = s.split()
59 lines = []
60 t = []
61 n = 0
62 for w in words:
63 if n + len(w) >= width:
64 lines.append(' '.join(t))
65 n = indent
66 t = [' '*(indent-1)]
68 t.append(w)
69 n += len(w) + 1
71 lines.append(' '.join(t))
72 return '\n'.join(lines)
75def same(x, eps=0.0):
76 if any(type(x[0]) != type(r) for r in x):
77 return False
79 if isinstance(x[0], float):
80 return all(abs(r-x[0]) <= eps for r in x)
81 else:
82 return all(r == x[0] for r in x)
85class InconsistentResponseInformation(Exception):
86 pass
89def check_resp(resp, value, frequency, limit_db, prelude=''):
90 if not value or not frequency:
91 logger.warn('Cannot validate frequency response')
92 return
94 value_resp = num.abs(
95 resp.evaluate(num.array([frequency], dtype=float)))[0]
97 if value_resp == 0.0:
98 raise InconsistentResponseInformation(
99 '%s\n'
100 ' computed response is zero' % prelude)
102 diff_db = 20.0 * num.log10(value_resp/value)
104 if num.abs(diff_db) > limit_db:
105 raise InconsistentResponseInformation(
106 '%s\n'
107 ' reported value: %g\n'
108 ' computed value: %g\n'
109 ' at frequency [Hz]: %g\n'
110 ' difference [dB]: %g\n'
111 ' limit [dB]: %g' % (
112 prelude,
113 value,
114 value_resp,
115 frequency,
116 diff_db,
117 limit_db))
120this_year = time.gmtime()[0]
123class DummyAwareOptionalTimestamp(Object):
124 dummy_for = (hpfloat, float)
125 dummy_for_description = 'time_float'
127 class __T(TBase):
129 def regularize_extra(self, val):
130 time_float = get_time_float()
132 if isinstance(val, datetime.datetime):
133 tt = val.utctimetuple()
134 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6
136 elif isinstance(val, datetime.date):
137 tt = val.timetuple()
138 val = time_float(calendar.timegm(tt))
140 elif isinstance(val, (str, newstr)):
141 val = val.strip()
143 tz_offset = 0
145 m = re_tz.search(val)
146 if m:
147 sh = m.group(2)
148 sm = m.group(4)
149 tz_offset = (int(sh)*3600 if sh else 0) \
150 + (int(sm)*60 if sm else 0)
152 val = re_tz.sub('', val)
154 if len(val) > 10 and val[10] == 'T':
155 val = val.replace('T', ' ', 1)
157 try:
158 val = util.str_to_time(val) - tz_offset
160 except util.TimeStrError:
161 year = int(val[:4])
162 if sys.maxsize > 2**32: # if we're on 64bit
163 if year > this_year + 100:
164 return None # StationXML contained a dummy date
166 if year < 1903: # for macOS, 1900-01-01 dummy dates
167 return None
169 else: # 32bit end of time is in 2038
170 if this_year < 2037 and year > 2037 or year < 1903:
171 return None # StationXML contained a dummy date
173 raise
175 elif isinstance(val, (int, float)):
176 val = time_float(val)
178 else:
179 raise ValidationError(
180 '%s: cannot convert "%s" to type %s' % (
181 self.xname(), val, time_float))
183 return val
185 def to_save(self, val):
186 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\
187 .rstrip('0').rstrip('.')
189 def to_save_xml(self, val):
190 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\
191 .rstrip('0').rstrip('.') + 'Z'
194class Nominal(StringChoice):
195 choices = [
196 'NOMINAL',
197 'CALCULATED']
200class Email(UnicodePattern):
201 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+'
204class RestrictedStatus(StringChoice):
205 choices = [
206 'open',
207 'closed',
208 'partial']
211class Type(StringChoice):
212 choices = [
213 'TRIGGERED',
214 'CONTINUOUS',
215 'HEALTH',
216 'GEOPHYSICAL',
217 'WEATHER',
218 'FLAG',
219 'SYNTHESIZED',
220 'INPUT',
221 'EXPERIMENTAL',
222 'MAINTENANCE',
223 'BEAM']
225 class __T(StringChoice.T):
226 def validate_extra(self, val):
227 if val not in self.choices:
228 logger.warn(
229 'channel type: "%s" is not a valid choice out of %s' %
230 (val, repr(self.choices)))
233class PzTransferFunction(StringChoice):
234 choices = [
235 'LAPLACE (RADIANS/SECOND)',
236 'LAPLACE (HERTZ)',
237 'DIGITAL (Z-TRANSFORM)']
240class Symmetry(StringChoice):
241 choices = [
242 'NONE',
243 'EVEN',
244 'ODD']
247class CfTransferFunction(StringChoice):
249 class __T(StringChoice.T):
250 def validate(self, val, regularize=False, depth=-1):
251 if regularize:
252 try:
253 val = str(val)
254 except ValueError:
255 raise ValidationError(
256 '%s: cannot convert to string %s' % (self.xname,
257 repr(val)))
259 val = self.dummy_cls.replacements.get(val, val)
261 self.validate_extra(val)
262 return val
264 choices = [
265 'ANALOG (RADIANS/SECOND)',
266 'ANALOG (HERTZ)',
267 'DIGITAL']
269 replacements = {
270 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)',
271 'ANALOG (HZ)': 'ANALOG (HERTZ)',
272 }
275class Approximation(StringChoice):
276 choices = [
277 'MACLAURIN']
280class PhoneNumber(StringPattern):
281 pattern = '[0-9]+-[0-9]+'
284class Site(Object):
285 '''
286 Description of a site location using name and optional geopolitical
287 boundaries (country, city, etc.).
288 '''
290 name = Unicode.T(xmltagname='Name')
291 description = Unicode.T(optional=True, xmltagname='Description')
292 town = Unicode.T(optional=True, xmltagname='Town')
293 county = Unicode.T(optional=True, xmltagname='County')
294 region = Unicode.T(optional=True, xmltagname='Region')
295 country = Unicode.T(optional=True, xmltagname='Country')
298class ExternalReference(Object):
299 '''
300 This type contains a URI and description for external data that users may
301 want to reference in StationXML.
302 '''
304 uri = String.T(xmltagname='URI')
305 description = Unicode.T(xmltagname='Description')
308class Units(Object):
309 '''
310 A type to document units. Corresponds to SEED blockette 34.
311 '''
313 def __init__(self, name=None, **kwargs):
314 Object.__init__(self, name=name, **kwargs)
316 name = String.T(xmltagname='Name')
317 description = Unicode.T(optional=True, xmltagname='Description')
320class Counter(Int):
321 pass
324class SampleRateRatio(Object):
325 '''
326 Sample rate expressed as number of samples in a number of seconds.
327 '''
329 number_samples = Int.T(xmltagname='NumberSamples')
330 number_seconds = Int.T(xmltagname='NumberSeconds')
333class Gain(Object):
334 '''
335 Complex type for sensitivity and frequency ranges. This complex type can be
336 used to represent both overall sensitivities and individual stage gains.
337 The FrequencyRangeGroup is an optional construct that defines a pass band
338 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue
339 is valid within the number of decibels specified in FrequencyDBVariation.
340 '''
342 def __init__(self, value=None, **kwargs):
343 Object.__init__(self, value=value, **kwargs)
345 value = Float.T(optional=True, xmltagname='Value')
346 frequency = Float.T(optional=True, xmltagname='Frequency')
349class NumeratorCoefficient(Object):
350 i = Int.T(optional=True, xmlstyle='attribute')
351 value = Float.T(xmlstyle='content')
354class FloatNoUnit(Object):
355 def __init__(self, value=None, **kwargs):
356 Object.__init__(self, value=value, **kwargs)
358 plus_error = Float.T(optional=True, xmlstyle='attribute')
359 minus_error = Float.T(optional=True, xmlstyle='attribute')
360 value = Float.T(xmlstyle='content')
363class FloatWithUnit(FloatNoUnit):
364 unit = String.T(optional=True, xmlstyle='attribute')
367class Equipment(Object):
368 resource_id = String.T(optional=True, xmlstyle='attribute')
369 type = String.T(optional=True, xmltagname='Type')
370 description = Unicode.T(optional=True, xmltagname='Description')
371 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer')
372 vendor = Unicode.T(optional=True, xmltagname='Vendor')
373 model = Unicode.T(optional=True, xmltagname='Model')
374 serial_number = String.T(optional=True, xmltagname='SerialNumber')
375 installation_date = DummyAwareOptionalTimestamp.T(
376 optional=True,
377 xmltagname='InstallationDate')
378 removal_date = DummyAwareOptionalTimestamp.T(
379 optional=True,
380 xmltagname='RemovalDate')
381 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate'))
384class PhoneNumber(Object):
385 description = Unicode.T(optional=True, xmlstyle='attribute')
386 country_code = Int.T(optional=True, xmltagname='CountryCode')
387 area_code = Int.T(xmltagname='AreaCode')
388 phone_number = PhoneNumber.T(xmltagname='PhoneNumber')
391class BaseFilter(Object):
392 '''
393 The BaseFilter is derived by all filters.
394 '''
396 resource_id = String.T(optional=True, xmlstyle='attribute')
397 name = String.T(optional=True, xmlstyle='attribute')
398 description = Unicode.T(optional=True, xmltagname='Description')
399 input_units = Units.T(optional=True, xmltagname='InputUnits')
400 output_units = Units.T(optional=True, xmltagname='OutputUnits')
403class Sensitivity(Gain):
404 '''
405 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional
406 construct that defines a pass band in Hertz (FrequencyStart and
407 FrequencyEnd) in which the SensitivityValue is valid within the number of
408 decibels specified in FrequencyDBVariation.
409 '''
411 input_units = Units.T(optional=True, xmltagname='InputUnits')
412 output_units = Units.T(optional=True, xmltagname='OutputUnits')
413 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart')
414 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd')
415 frequency_db_variation = Float.T(optional=True,
416 xmltagname='FrequencyDBVariation')
419class Coefficient(FloatNoUnit):
420 number = Counter.T(optional=True, xmlstyle='attribute')
423class PoleZero(Object):
424 '''
425 Complex numbers used as poles or zeros in channel response.
426 '''
428 number = Int.T(optional=True, xmlstyle='attribute')
429 real = FloatNoUnit.T(xmltagname='Real')
430 imaginary = FloatNoUnit.T(xmltagname='Imaginary')
432 def value(self):
433 return self.real.value + 1J * self.imaginary.value
436class ClockDrift(FloatWithUnit):
437 unit = String.T(default='SECONDS/SAMPLE', optional=True,
438 xmlstyle='attribute') # fixed
441class Second(FloatWithUnit):
442 '''
443 A time value in seconds.
444 '''
446 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute')
447 # fixed unit
450class Voltage(FloatWithUnit):
451 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute')
452 # fixed unit
455class Angle(FloatWithUnit):
456 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
457 # fixed unit
460class Azimuth(FloatWithUnit):
461 '''
462 Instrument azimuth, degrees clockwise from North.
463 '''
465 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
466 # fixed unit
469class Dip(FloatWithUnit):
470 '''
471 Instrument dip in degrees down from horizontal. Together azimuth and dip
472 describe the direction of the sensitive axis of the instrument.
473 '''
475 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
476 # fixed unit
479class Distance(FloatWithUnit):
480 '''
481 Extension of FloatWithUnit for distances, elevations, and depths.
482 '''
484 unit = String.T(default='METERS', optional=True, xmlstyle='attribute')
485 # NOT fixed unit!
488class Frequency(FloatWithUnit):
489 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute')
490 # fixed unit
493class SampleRate(FloatWithUnit):
494 '''
495 Sample rate in samples per second.
496 '''
498 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute')
499 # fixed unit
502class Person(Object):
503 '''
504 Representation of a person's contact information. A person can belong to
505 multiple agencies and have multiple email addresses and phone numbers.
506 '''
508 name_list = List.T(Unicode.T(xmltagname='Name'))
509 agency_list = List.T(Unicode.T(xmltagname='Agency'))
510 email_list = List.T(Email.T(xmltagname='Email'))
511 phone_list = List.T(PhoneNumber.T(xmltagname='Phone'))
514class FIR(BaseFilter):
515 '''
516 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are
517 also commonly documented using the Coefficients element.
518 '''
520 symmetry = Symmetry.T(xmltagname='Symmetry')
521 numerator_coefficient_list = List.T(
522 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient'))
525class Coefficients(BaseFilter):
526 '''
527 Response: coefficients for FIR filter. Laplace transforms or IIR filters
528 can be expressed using type as well but the PolesAndZeros should be used
529 instead. Corresponds to SEED blockette 54.
530 '''
532 cf_transfer_function_type = CfTransferFunction.T(
533 xmltagname='CfTransferFunctionType')
534 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator'))
535 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator'))
538class Latitude(FloatWithUnit):
539 '''
540 Type for latitude coordinate.
541 '''
543 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
544 # fixed unit
545 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
548class Longitude(FloatWithUnit):
549 '''
550 Type for longitude coordinate.
551 '''
553 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute')
554 # fixed unit
555 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute')
558class PolesZeros(BaseFilter):
559 '''
560 Response: complex poles and zeros. Corresponds to SEED blockette 53.
561 '''
563 pz_transfer_function_type = PzTransferFunction.T(
564 xmltagname='PzTransferFunctionType')
565 normalization_factor = Float.T(default=1.0,
566 xmltagname='NormalizationFactor')
567 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency')
568 zero_list = List.T(PoleZero.T(xmltagname='Zero'))
569 pole_list = List.T(PoleZero.T(xmltagname='Pole'))
571 def get_pyrocko_response(self, nslc):
572 if self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)':
573 logger.warn(
574 'unhandled pole-zero response of type "DIGITAL (Z-TRANSFORM)" '
575 '(%s)' % '.'.join(nslc))
577 return []
579 if self.pz_transfer_function_type not in (
580 'LAPLACE (RADIANS/SECOND)',
581 'LAPLACE (HERTZ)'):
583 raise NoResponseInformation(
584 'cannot convert PoleZero response of type %s (%s)' %
585 (self.pz_transfer_function_type, '.'.join(nslc)))
587 factor = 1.0
588 cfactor = 1.0
589 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
590 factor = 2. * math.pi
591 cfactor = (2. * math.pi)**(
592 len(self.pole_list) - len(self.zero_list))
594 if self.normalization_factor is None \
595 or self.normalization_factor == 0.0:
597 logger.warn(
598 'no pole-zero normalization factor given. Assuming a value of '
599 '1.0 (%s)' % '.'.join(nslc))
601 nfactor = 1.0
602 else:
603 nfactor = self.normalization_factor
605 resp = trace.PoleZeroResponse(
606 constant=nfactor*cfactor,
607 zeros=[z.value()*factor for z in self.zero_list],
608 poles=[p.value()*factor for p in self.pole_list])
610 if not self.normalization_frequency.value:
611 logger.warn(
612 'cannot check pole-zero normalization factor (%s)'
613 % '.'.join(nslc))
615 else:
616 computed_normalization_factor = nfactor / abs(
617 resp.evaluate(
618 num.array([self.normalization_frequency.value]))[0])
620 db = 20.0 * num.log10(
621 computed_normalization_factor / nfactor)
623 if abs(db) > 0.17:
624 logger.warn(
625 'computed and reported normalization factors differ by '
626 '%g dB: computed: %g, reported: %g (%s)' % (
627 db,
628 computed_normalization_factor,
629 nfactor,
630 '.'.join(nslc)))
632 return [resp]
635class ResponseListElement(Object):
636 frequency = Frequency.T(xmltagname='Frequency')
637 amplitude = FloatWithUnit.T(xmltagname='Amplitude')
638 phase = Angle.T(xmltagname='Phase')
641class Polynomial(BaseFilter):
642 '''
643 Response: expressed as a polynomial (allows non-linear sensors to be
644 described). Corresponds to SEED blockette 62. Can be used to describe a
645 stage of acquisition or a complete system.
646 '''
648 approximation_type = Approximation.T(default='MACLAURIN',
649 xmltagname='ApproximationType')
650 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound')
651 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound')
652 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound')
653 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound')
654 maximum_error = Float.T(xmltagname='MaximumError')
655 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient'))
658class Decimation(Object):
659 '''
660 Corresponds to SEED blockette 57.
661 '''
663 input_sample_rate = Frequency.T(xmltagname='InputSampleRate')
664 factor = Int.T(xmltagname='Factor')
665 offset = Int.T(xmltagname='Offset')
666 delay = FloatWithUnit.T(xmltagname='Delay')
667 correction = FloatWithUnit.T(xmltagname='Correction')
670class Operator(Object):
671 agency_list = List.T(Unicode.T(xmltagname='Agency'))
672 contact_list = List.T(Person.T(xmltagname='Contact'))
673 web_site = String.T(optional=True, xmltagname='WebSite')
676class Comment(Object):
677 '''
678 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51
679 and 59.
680 '''
682 id = Counter.T(optional=True, xmlstyle='attribute')
683 value = Unicode.T(xmltagname='Value')
684 begin_effective_time = DummyAwareOptionalTimestamp.T(
685 optional=True,
686 xmltagname='BeginEffectiveTime')
687 end_effective_time = DummyAwareOptionalTimestamp.T(
688 optional=True,
689 xmltagname='EndEffectiveTime')
690 author_list = List.T(Person.T(xmltagname='Author'))
693class ResponseList(BaseFilter):
694 '''
695 Response: list of frequency, amplitude and phase values. Corresponds to
696 SEED blockette 55.
697 '''
699 response_list_element_list = List.T(
700 ResponseListElement.T(xmltagname='ResponseListElement'))
703class Log(Object):
704 '''
705 Container for log entries.
706 '''
708 entry_list = List.T(Comment.T(xmltagname='Entry'))
711class ResponseStage(Object):
712 '''
713 This complex type represents channel response and covers SEED blockettes 53
714 to 56.
715 '''
717 number = Counter.T(xmlstyle='attribute')
718 resource_id = String.T(optional=True, xmlstyle='attribute')
719 poles_zeros_list = List.T(
720 PolesZeros.T(optional=True, xmltagname='PolesZeros'))
721 coefficients_list = List.T(
722 Coefficients.T(optional=True, xmltagname='Coefficients'))
723 response_list = ResponseList.T(optional=True, xmltagname='ResponseList')
724 fir = FIR.T(optional=True, xmltagname='FIR')
725 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial')
726 decimation = Decimation.T(optional=True, xmltagname='Decimation')
727 stage_gain = Gain.T(optional=True, xmltagname='StageGain')
729 def get_pyrocko_response(self, nslc):
730 responses = []
731 for pzs in self.poles_zeros_list:
732 responses.extend(pzs.get_pyrocko_response(nslc))
734 if len(self.poles_zeros_list) > 1:
735 logger.warn(
736 'multiple poles and zeros records in single response stage '
737 '(%s.%s.%s.%s)' % nslc)
739 if (self.coefficients_list or
740 self.response_list or
741 self.fir or
742 self.polynomial):
744 logger.debug('unhandled response at stage %i' % self.number)
746 if self.stage_gain:
747 responses.append(
748 trace.PoleZeroResponse(constant=self.stage_gain.value))
750 return responses
752 @property
753 def input_units(self):
754 for e in (self.poles_zeros_list + self.coefficients_list +
755 [self.response_list, self.fir, self.polynomial]):
756 if e is not None:
757 return e.input_units
759 @property
760 def output_units(self):
761 for e in (self.poles_zeros_list + self.coefficients_list +
762 [self.response_list, self.fir, self.polynomial]):
763 if e is not None:
764 return e.output_units
767class Response(Object):
768 resource_id = String.T(optional=True, xmlstyle='attribute')
769 instrument_sensitivity = Sensitivity.T(optional=True,
770 xmltagname='InstrumentSensitivity')
771 instrument_polynomial = Polynomial.T(optional=True,
772 xmltagname='InstrumentPolynomial')
773 stage_list = List.T(ResponseStage.T(xmltagname='Stage'))
775 def get_pyrocko_response(self, nslc, fake_input_units=None):
776 responses = []
777 for stage in self.stage_list:
778 responses.extend(stage.get_pyrocko_response(nslc))
780 if not self.stage_list and self.instrument_sensitivity:
781 responses.append(
782 trace.PoleZeroResponse(
783 constant=self.instrument_sensitivity.value))
785 if self.instrument_sensitivity:
786 trial = trace.MultiplyResponse(responses)
787 sval = self.instrument_sensitivity.value
788 sfreq = self.instrument_sensitivity.frequency
789 try:
790 check_resp(
791 trial, sval, sfreq, 0.1,
792 prelude='Instrument sensitivity value inconsistent with '
793 'sensitivity computed from complete response\n'
794 ' channel: %s' % '.'.join(nslc))
795 except InconsistentResponseInformation as e:
796 logger.warn(str(e))
798 if fake_input_units is not None:
799 if not self.instrument_sensitivity or \
800 self.instrument_sensitivity.input_units is None:
802 raise NoResponseInformation('no input units given')
804 input_units = self.instrument_sensitivity.input_units.name
806 try:
807 conresp = conversion[
808 fake_input_units.upper(), input_units.upper()]
810 except KeyError:
811 raise NoResponseInformation(
812 'cannot convert between units: %s, %s'
813 % (fake_input_units, input_units))
815 if conresp is not None:
816 responses.append(conresp)
818 return trace.MultiplyResponse(responses)
820 @classmethod
821 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit,
822 normalization_frequency=1.0):
823 '''
824 Convert Pyrocko pole-zero response to StationXML response.
826 :param presponse: Pyrocko pole-zero response
827 :type presponse: :py:class:`~pyrocko.trace.PoleZeroResponse`
828 :param input_unit: Input unit to be reported in the StationXML
829 response.
830 :type input_unit: str
831 :param output_unit: Output unit to be reported in the StationXML
832 response.
833 :type output_unit: str
834 :param normalization_frequency: Frequency where the normalization
835 factor for the StationXML response should be computed.
836 :type normalization_frequency: float
837 '''
839 norm_factor = 1.0/float(abs(
840 presponse.evaluate(num.array([normalization_frequency]))[0]
841 / presponse.constant))
843 pzs = PolesZeros(
844 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)',
845 normalization_factor=norm_factor,
846 normalization_frequency=Frequency(normalization_frequency),
847 zero_list=[PoleZero(real=FloatNoUnit(z.real),
848 imaginary=FloatNoUnit(z.imag))
849 for z in presponse.zeros],
850 pole_list=[PoleZero(real=FloatNoUnit(z.real),
851 imaginary=FloatNoUnit(z.imag))
852 for z in presponse.poles])
854 pzs.validate()
856 stage = ResponseStage(
857 number=1,
858 poles_zeros_list=[pzs],
859 stage_gain=Gain(float(abs(presponse.constant))/norm_factor))
861 resp = Response(
862 instrument_sensitivity=Sensitivity(
863 value=stage.stage_gain.value,
864 input_units=Units(input_unit),
865 output_units=Units(output_unit)),
867 stage_list=[stage])
869 return resp
872class BaseNode(Object):
873 '''
874 A base node type for derivation from: Network, Station and Channel types.
875 '''
877 code = String.T(xmlstyle='attribute')
878 start_date = DummyAwareOptionalTimestamp.T(optional=True,
879 xmlstyle='attribute')
880 end_date = DummyAwareOptionalTimestamp.T(optional=True,
881 xmlstyle='attribute')
882 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute')
883 alternate_code = String.T(optional=True, xmlstyle='attribute')
884 historical_code = String.T(optional=True, xmlstyle='attribute')
885 description = Unicode.T(optional=True, xmltagname='Description')
886 comment_list = List.T(Comment.T(xmltagname='Comment'))
888 def spans(self, *args):
889 if len(args) == 0:
890 return True
891 elif len(args) == 1:
892 return ((self.start_date is None or
893 self.start_date <= args[0]) and
894 (self.end_date is None or
895 args[0] <= self.end_date))
897 elif len(args) == 2:
898 return ((self.start_date is None or
899 args[1] >= self.start_date) and
900 (self.end_date is None or
901 self.end_date >= args[0]))
904class Channel(BaseNode):
905 '''
906 Equivalent to SEED blockette 52 and parent element for the related the
907 response blockettes.
908 '''
910 location_code = String.T(xmlstyle='attribute')
911 external_reference_list = List.T(
912 ExternalReference.T(xmltagname='ExternalReference'))
913 latitude = Latitude.T(xmltagname='Latitude')
914 longitude = Longitude.T(xmltagname='Longitude')
915 elevation = Distance.T(xmltagname='Elevation')
916 depth = Distance.T(xmltagname='Depth')
917 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth')
918 dip = Dip.T(optional=True, xmltagname='Dip')
919 type_list = List.T(Type.T(xmltagname='Type'))
920 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate')
921 sample_rate_ratio = SampleRateRatio.T(optional=True,
922 xmltagname='SampleRateRatio')
923 storage_format = String.T(optional=True, xmltagname='StorageFormat')
924 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift')
925 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits')
926 sensor = Equipment.T(optional=True, xmltagname='Sensor')
927 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier')
928 data_logger = Equipment.T(optional=True, xmltagname='DataLogger')
929 equipment = Equipment.T(optional=True, xmltagname='Equipment')
930 response = Response.T(optional=True, xmltagname='Response')
932 @property
933 def position_values(self):
934 lat = self.latitude.value
935 lon = self.longitude.value
936 elevation = value_or_none(self.elevation)
937 depth = value_or_none(self.depth)
938 return lat, lon, elevation, depth
941class Station(BaseNode):
942 '''
943 This type represents a Station epoch. It is common to only have a single
944 station epoch with the station's creation and termination dates as the
945 epoch start and end dates.
946 '''
948 latitude = Latitude.T(xmltagname='Latitude')
949 longitude = Longitude.T(xmltagname='Longitude')
950 elevation = Distance.T(xmltagname='Elevation')
951 site = Site.T(optional=True, xmltagname='Site')
952 vault = Unicode.T(optional=True, xmltagname='Vault')
953 geology = Unicode.T(optional=True, xmltagname='Geology')
954 equipment_list = List.T(Equipment.T(xmltagname='Equipment'))
955 operator_list = List.T(Operator.T(xmltagname='Operator'))
956 creation_date = DummyAwareOptionalTimestamp.T(
957 optional=True, xmltagname='CreationDate')
958 termination_date = DummyAwareOptionalTimestamp.T(
959 optional=True, xmltagname='TerminationDate')
960 total_number_channels = Counter.T(
961 optional=True, xmltagname='TotalNumberChannels')
962 selected_number_channels = Counter.T(
963 optional=True, xmltagname='SelectedNumberChannels')
964 external_reference_list = List.T(
965 ExternalReference.T(xmltagname='ExternalReference'))
966 channel_list = List.T(Channel.T(xmltagname='Channel'))
968 @property
969 def position_values(self):
970 lat = self.latitude.value
971 lon = self.longitude.value
972 elevation = value_or_none(self.elevation)
973 return lat, lon, elevation
976class Network(BaseNode):
977 '''
978 This type represents the Network layer, all station metadata is contained
979 within this element. The official name of the network or other descriptive
980 information can be included in the Description element. The Network can
981 contain 0 or more Stations.
982 '''
984 total_number_stations = Counter.T(optional=True,
985 xmltagname='TotalNumberStations')
986 selected_number_stations = Counter.T(optional=True,
987 xmltagname='SelectedNumberStations')
988 station_list = List.T(Station.T(xmltagname='Station'))
990 @property
991 def station_code_list(self):
992 return sorted(set(s.code for s in self.station_list))
994 @property
995 def sl_code_list(self):
996 sls = set()
997 for station in self.station_list:
998 for channel in station.channel_list:
999 sls.add((station.code, channel.location_code))
1001 return sorted(sls)
1003 def summary(self, width=80, indent=4):
1004 sls = self.sl_code_list or [(x,) for x in self.station_code_list]
1005 lines = ['%s (%i):' % (self.code, len(sls))]
1006 if sls:
1007 ssls = ['.'.join(x for x in c if x) for c in sls]
1008 w = max(len(x) for x in ssls)
1009 n = (width - indent) / (w+1)
1010 while ssls:
1011 lines.append(
1012 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n]))
1014 ssls[:n] = []
1016 return '\n'.join(lines)
1019def value_or_none(x):
1020 if x is not None:
1021 return x.value
1022 else:
1023 return None
1026def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'):
1028 pos = lat, lon, elevation, depth = \
1029 channels[0].position_values
1031 if not all(pos == x.position_values for x in channels):
1032 info = '\n'.join(
1033 ' %s: %s' % (x.code, x.position_values) for
1034 x in channels)
1036 mess = 'encountered inconsistencies in channel ' \
1037 'lat/lon/elevation/depth ' \
1038 'for %s.%s.%s: \n%s' % (nsl + (info,))
1040 if inconsistencies == 'raise':
1041 raise InconsistentChannelLocations(mess)
1043 elif inconsistencies == 'warn':
1044 logger.warn(mess)
1045 logger.warn(' -> using mean values')
1047 apos = num.array([x.position_values for x in channels], dtype=float)
1048 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \
1049 / num.sum(num.isfinite(apos), axis=0)
1051 groups = {}
1052 for channel in channels:
1053 if channel.code not in groups:
1054 groups[channel.code] = []
1056 groups[channel.code].append(channel)
1058 pchannels = []
1059 for code in sorted(groups.keys()):
1060 data = [
1061 (channel.code, value_or_none(channel.azimuth),
1062 value_or_none(channel.dip))
1063 for channel in groups[code]]
1065 azimuth, dip = util.consistency_merge(
1066 data,
1067 message='channel orientation values differ:',
1068 error=inconsistencies)
1070 pchannels.append(
1071 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip))
1073 return pyrocko.model.Station(
1074 *nsl,
1075 lat=mlat,
1076 lon=mlon,
1077 elevation=mele,
1078 depth=mdep,
1079 channels=pchannels)
1082class FDSNStationXML(Object):
1083 '''
1084 Top-level type for Station XML. Required field are Source (network ID of
1085 the institution sending the message) and one or more Network containers or
1086 one or more Station containers.
1087 '''
1089 schema_version = Float.T(default=1.0, xmlstyle='attribute')
1090 source = String.T(xmltagname='Source')
1091 sender = String.T(optional=True, xmltagname='Sender')
1092 module = String.T(optional=True, xmltagname='Module')
1093 module_uri = String.T(optional=True, xmltagname='ModuleURI')
1094 created = Timestamp.T(optional=True, xmltagname='Created')
1095 network_list = List.T(Network.T(xmltagname='Network'))
1097 xmltagname = 'FDSNStationXML'
1098 guessable_xmlns = [guts_xmlns]
1100 def get_pyrocko_stations(self, nslcs=None, nsls=None,
1101 time=None, timespan=None,
1102 inconsistencies='warn'):
1104 assert inconsistencies in ('raise', 'warn')
1106 if nslcs is not None:
1107 nslcs = set(nslcs)
1109 if nsls is not None:
1110 nsls = set(nsls)
1112 tt = ()
1113 if time is not None:
1114 tt = (time,)
1115 elif timespan is not None:
1116 tt = timespan
1118 pstations = []
1119 for network in self.network_list:
1120 if not network.spans(*tt):
1121 continue
1123 for station in network.station_list:
1124 if not station.spans(*tt):
1125 continue
1127 if station.channel_list:
1128 loc_to_channels = {}
1129 for channel in station.channel_list:
1130 if not channel.spans(*tt):
1131 continue
1133 loc = channel.location_code.strip()
1134 if loc not in loc_to_channels:
1135 loc_to_channels[loc] = []
1137 loc_to_channels[loc].append(channel)
1139 for loc in sorted(loc_to_channels.keys()):
1140 channels = loc_to_channels[loc]
1141 if nslcs is not None:
1142 channels = [channel for channel in channels
1143 if (network.code, station.code, loc,
1144 channel.code) in nslcs]
1146 if not channels:
1147 continue
1149 nsl = network.code, station.code, loc
1150 if nsls is not None and nsl not in nsls:
1151 continue
1153 pstations.append(
1154 pyrocko_station_from_channels(
1155 nsl,
1156 channels,
1157 inconsistencies=inconsistencies))
1158 else:
1159 pstations.append(pyrocko.model.Station(
1160 network.code, station.code, '*',
1161 lat=station.latitude.value,
1162 lon=station.longitude.value,
1163 elevation=value_or_none(station.elevation),
1164 name=station.description or ''))
1166 return pstations
1168 @classmethod
1169 def from_pyrocko_stations(
1170 cls, pyrocko_stations, add_flat_responses_from=None):
1172 '''
1173 Generate :py:class:`FDSNStationXML` from list of
1174 :py:class;`pyrocko.model.Station` instances.
1176 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station`
1177 instances.
1178 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2'
1179 '''
1180 from collections import defaultdict
1181 network_dict = defaultdict(list)
1183 if add_flat_responses_from:
1184 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2')
1185 extra = dict(
1186 response=Response(
1187 instrument_sensitivity=Sensitivity(
1188 value=1.0,
1189 frequency=1.0,
1190 input_units=Units(name=add_flat_responses_from))))
1191 else:
1192 extra = {}
1194 have_offsets = set()
1195 for s in pyrocko_stations:
1197 if s.north_shift != 0.0 or s.east_shift != 0.0:
1198 have_offsets.add(s.nsl())
1200 network, station, location = s.nsl()
1201 channel_list = []
1202 for c in s.channels:
1203 channel_list.append(
1204 Channel(
1205 location_code=location,
1206 code=c.name,
1207 latitude=Latitude(value=s.effective_lat),
1208 longitude=Longitude(value=s.effective_lon),
1209 elevation=Distance(value=s.elevation),
1210 depth=Distance(value=s.depth),
1211 azimuth=Azimuth(value=c.azimuth),
1212 dip=Dip(value=c.dip),
1213 **extra
1214 )
1215 )
1217 network_dict[network].append(
1218 Station(
1219 code=station,
1220 latitude=Latitude(value=s.effective_lat),
1221 longitude=Longitude(value=s.effective_lon),
1222 elevation=Distance(value=s.elevation),
1223 channel_list=channel_list)
1224 )
1226 if have_offsets:
1227 logger.warn(
1228 'StationXML does not support Cartesian offsets in '
1229 'coordinates. Storing effective lat/lon for stations: %s' %
1230 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets)))
1232 timestamp = util.to_time_float(time.time())
1233 network_list = []
1234 for k, station_list in network_dict.items():
1236 network_list.append(
1237 Network(
1238 code=k, station_list=station_list,
1239 total_number_stations=len(station_list)))
1241 sxml = FDSNStationXML(
1242 source='from pyrocko stations list', created=timestamp,
1243 network_list=network_list)
1245 sxml.validate()
1246 return sxml
1248 def iter_network_stations(
1249 self, net=None, sta=None, time=None, timespan=None):
1251 tt = ()
1252 if time is not None:
1253 tt = (time,)
1254 elif timespan is not None:
1255 tt = timespan
1257 for network in self.network_list:
1258 if not network.spans(*tt) or (
1259 net is not None and network.code != net):
1260 continue
1262 for station in network.station_list:
1263 if not station.spans(*tt) or (
1264 sta is not None and station.code != sta):
1265 continue
1267 yield (network, station)
1269 def iter_network_station_channels(
1270 self, net=None, sta=None, loc=None, cha=None,
1271 time=None, timespan=None):
1273 if loc is not None:
1274 loc = loc.strip()
1276 tt = ()
1277 if time is not None:
1278 tt = (time,)
1279 elif timespan is not None:
1280 tt = timespan
1282 for network in self.network_list:
1283 if not network.spans(*tt) or (
1284 net is not None and network.code != net):
1285 continue
1287 for station in network.station_list:
1288 if not station.spans(*tt) or (
1289 sta is not None and station.code != sta):
1290 continue
1292 if station.channel_list:
1293 for channel in station.channel_list:
1294 if (not channel.spans(*tt) or
1295 (cha is not None and channel.code != cha) or
1296 (loc is not None and
1297 channel.location_code.strip() != loc)):
1298 continue
1300 yield (network, station, channel)
1302 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None,
1303 time=None, timespan=None):
1305 groups = {}
1306 for network, station, channel in self.iter_network_station_channels(
1307 net, sta, loc, cha, time=time, timespan=timespan):
1309 net = network.code
1310 sta = station.code
1311 cha = channel.code
1312 loc = channel.location_code.strip()
1313 if len(cha) == 3:
1314 bic = cha[:2] # band and intrument code according to SEED
1315 elif len(cha) == 1:
1316 bic = ''
1317 else:
1318 bic = cha
1320 if channel.response and \
1321 channel.response.instrument_sensitivity and \
1322 channel.response.instrument_sensitivity.input_units:
1324 unit = channel.response.instrument_sensitivity.input_units.name
1325 else:
1326 unit = None
1328 bic = (bic, unit)
1330 k = net, sta, loc
1331 if k not in groups:
1332 groups[k] = {}
1334 if bic not in groups[k]:
1335 groups[k][bic] = []
1337 groups[k][bic].append(channel)
1339 for nsl, bic_to_channels in groups.items():
1340 bad_bics = []
1341 for bic, channels in bic_to_channels.items():
1342 sample_rates = []
1343 for channel in channels:
1344 sample_rates.append(channel.sample_rate.value)
1346 if not same(sample_rates):
1347 scs = ','.join(channel.code for channel in channels)
1348 srs = ', '.join('%e' % x for x in sample_rates)
1349 err = 'ignoring channels with inconsistent sampling ' + \
1350 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs))
1352 logger.warn(err)
1353 bad_bics.append(bic)
1355 for bic in bad_bics:
1356 del bic_to_channels[bic]
1358 return groups
1360 def choose_channels(
1361 self,
1362 target_sample_rate=None,
1363 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'],
1364 priority_units=['M/S', 'M/S**2'],
1365 priority_instrument_code=['H', 'L'],
1366 time=None,
1367 timespan=None):
1369 nslcs = {}
1370 for nsl, bic_to_channels in self.get_channel_groups(
1371 time=time, timespan=timespan).items():
1373 useful_bics = []
1374 for bic, channels in bic_to_channels.items():
1375 rate = channels[0].sample_rate.value
1377 if target_sample_rate is not None and \
1378 rate < target_sample_rate*0.99999:
1379 continue
1381 if len(bic[0]) == 2:
1382 if bic[0][0] not in priority_band_code:
1383 continue
1385 if bic[0][1] not in priority_instrument_code:
1386 continue
1388 unit = bic[1]
1390 prio_unit = len(priority_units)
1391 try:
1392 prio_unit = priority_units.index(unit)
1393 except ValueError:
1394 pass
1396 prio_inst = len(priority_instrument_code)
1397 prio_band = len(priority_band_code)
1398 if len(channels[0].code) == 3:
1399 try:
1400 prio_inst = priority_instrument_code.index(
1401 channels[0].code[1])
1402 except ValueError:
1403 pass
1405 try:
1406 prio_band = priority_band_code.index(
1407 channels[0].code[0])
1408 except ValueError:
1409 pass
1411 if target_sample_rate is None:
1412 rate = -rate
1414 useful_bics.append((-len(channels), prio_band, rate, prio_unit,
1415 prio_inst, bic))
1417 useful_bics.sort()
1419 for _, _, rate, _, _, bic in useful_bics:
1420 channels = sorted(
1421 bic_to_channels[bic],
1422 key=lambda channel: channel.code)
1424 if channels:
1425 for channel in channels:
1426 nslcs[nsl + (channel.code,)] = channel
1428 break
1430 return nslcs
1432 def get_pyrocko_response(
1433 self, nslc, time=None, timespan=None, fake_input_units=None):
1435 net, sta, loc, cha = nslc
1436 resps = []
1437 for _, _, channel in self.iter_network_station_channels(
1438 net, sta, loc, cha, time=time, timespan=timespan):
1439 resp = channel.response
1440 if resp:
1441 resps.append(resp.get_pyrocko_response(
1442 nslc, fake_input_units=fake_input_units))
1444 if not resps:
1445 raise NoResponseInformation('%s.%s.%s.%s' % nslc)
1446 elif len(resps) > 1:
1447 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc)
1449 return resps[0]
1451 @property
1452 def n_code_list(self):
1453 return sorted(set(x.code for x in self.network_list))
1455 @property
1456 def ns_code_list(self):
1457 nss = set()
1458 for network in self.network_list:
1459 for station in network.station_list:
1460 nss.add((network.code, station.code))
1462 return sorted(nss)
1464 @property
1465 def nsl_code_list(self):
1466 nsls = set()
1467 for network in self.network_list:
1468 for station in network.station_list:
1469 for channel in station.channel_list:
1470 nsls.add(
1471 (network.code, station.code, channel.location_code))
1473 return sorted(nsls)
1475 @property
1476 def nslc_code_list(self):
1477 nslcs = set()
1478 for network in self.network_list:
1479 for station in network.station_list:
1480 for channel in station.channel_list:
1481 nslcs.add(
1482 (network.code, station.code, channel.location_code,
1483 channel.code))
1485 return sorted(nslcs)
1487 def summary(self):
1488 lst = [
1489 'number of n codes: %i' % len(self.n_code_list),
1490 'number of ns codes: %i' % len(self.ns_code_list),
1491 'number of nsl codes: %i' % len(self.nsl_code_list),
1492 'number of nslc codes: %i' % len(self.nslc_code_list)
1493 ]
1495 return '\n'.join(lst)
1498class InconsistentChannelLocations(Exception):
1499 pass
1502class InvalidRecord(Exception):
1503 def __init__(self, line):
1504 Exception.__init__(self)
1505 self._line = line
1507 def __str__(self):
1508 return 'Invalid record: "%s"' % self._line
1511def load_channel_table(stream):
1513 networks = {}
1514 stations = {}
1516 for line in stream:
1517 line = str(line.decode('ascii'))
1518 if line.startswith('#'):
1519 continue
1521 t = line.rstrip().split('|')
1523 if len(t) != 17:
1524 logger.warn('Invalid channel record: %s' % line)
1525 continue
1527 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale,
1528 scale_freq, scale_units, sample_rate, start_date, end_date) = t
1530 try:
1531 scale = float(scale)
1532 except ValueError:
1533 scale = None
1535 try:
1536 scale_freq = float(scale_freq)
1537 except ValueError:
1538 scale_freq = None
1540 try:
1541 depth = float(dep)
1542 except ValueError:
1543 depth = 0.0
1545 try:
1546 azi = float(azi)
1547 dip = float(dip)
1548 except ValueError:
1549 azi = None
1550 dip = None
1552 try:
1553 if net not in networks:
1554 network = Network(code=net)
1555 else:
1556 network = networks[net]
1558 if (net, sta) not in stations:
1559 station = Station(
1560 code=sta, latitude=lat, longitude=lon, elevation=ele)
1562 station.regularize()
1563 else:
1564 station = stations[net, sta]
1566 if scale:
1567 resp = Response(
1568 instrument_sensitivity=Sensitivity(
1569 value=scale,
1570 frequency=scale_freq,
1571 input_units=scale_units))
1572 else:
1573 resp = None
1575 channel = Channel(
1576 code=cha,
1577 location_code=loc.strip(),
1578 latitude=lat,
1579 longitude=lon,
1580 elevation=ele,
1581 depth=depth,
1582 azimuth=azi,
1583 dip=dip,
1584 sensor=Equipment(description=sens),
1585 response=resp,
1586 sample_rate=sample_rate,
1587 start_date=start_date,
1588 end_date=end_date or None)
1590 channel.regularize()
1592 except ValidationError:
1593 raise InvalidRecord(line)
1595 if net not in networks:
1596 networks[net] = network
1598 if (net, sta) not in stations:
1599 stations[net, sta] = station
1600 network.station_list.append(station)
1602 station.channel_list.append(channel)
1604 return FDSNStationXML(
1605 source='created from table input',
1606 created=time.time(),
1607 network_list=sorted(networks.values(), key=lambda x: x.code))
1610def primitive_merge(sxs):
1611 networks = []
1612 for sx in sxs:
1613 networks.extend(sx.network_list)
1615 return FDSNStationXML(
1616 source='merged from different sources',
1617 created=time.time(),
1618 network_list=copy.deepcopy(
1619 sorted(networks, key=lambda x: x.code)))