Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/stationxml.py: 70%

1161 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 06:59 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7`FDSN StationXML <https://www.fdsn.org/xml/station/>`_ input, output and data 

8model. 

9''' 

10 

11import sys 

12import time 

13import logging 

14import datetime 

15import calendar 

16import math 

17import copy 

18 

19import numpy as num 

20 

21from pyrocko.guts import (StringChoice, StringPattern, UnicodePattern, String, 

22 Unicode, Int, Float, List, Object, Timestamp, 

23 ValidationError, TBase, re_tz, Any, Tuple) 

24from pyrocko.guts import load_xml # noqa 

25from pyrocko.util import hpfloat, time_to_str, get_time_float 

26 

27import pyrocko.model 

28from pyrocko import util, response 

29 

30guts_prefix = 'sx' 

31 

32guts_xmlns = 'http://www.fdsn.org/xml/station/1' 

33 

34logger = logging.getLogger('pyrocko.io.stationxml') 

35 

36conversion = { 

37 ('M', 'M'): None, 

38 ('M/S', 'M'): response.IntegrationResponse(1), 

39 ('M/S**2', 'M'): response.IntegrationResponse(2), 

40 ('M', 'M/S'): response.DifferentiationResponse(1), 

41 ('M/S', 'M/S'): None, 

42 ('M/S**2', 'M/S'): response.IntegrationResponse(1), 

43 ('M', 'M/S**2'): response.DifferentiationResponse(2), 

44 ('M/S', 'M/S**2'): response.DifferentiationResponse(1), 

45 ('M/S**2', 'M/S**2'): None} 

46 

47 

48unit_to_quantity = { 

49 'M/S': 'velocity', 

50 'M': 'displacement', 

51 'M/S**2': 'acceleration', 

52 'V': 'voltage', 

53 'COUNTS': 'counts', 

54 'COUNT': 'counts', 

55 'PA': 'pressure', 

56 'RAD': 'rotation-displacement', 

57 'R': 'rotation-displacement', 

58 'RAD/S': 'rotation-velocity', 

59 'R/S': 'rotation-velocity', 

60 'RAD/S**2': 'rotation-acceleration', 

61 'R/S**2': 'rotation-acceleration'} 

62 

63 

64def to_quantity(unit, context, delivery): 

65 

66 if unit is None: 

67 return None 

68 

69 name = unit.name.upper() 

70 if name in unit_to_quantity: 

71 return unit_to_quantity[name] 

72 else: 

73 delivery.log.append(( 

74 'warning', 

75 'Units not supported by Squirrel framework: %s' % ( 

76 unit.name.upper() + ( 

77 ' (%s)' % unit.description 

78 if unit.description else '')), 

79 context)) 

80 

81 return 'unsupported_quantity(%s)' % unit 

82 

83 

84class StationXMLError(Exception): 

85 pass 

86 

87 

88class Inconsistencies(StationXMLError): 

89 pass 

90 

91 

92class NoResponseInformation(StationXMLError): 

93 pass 

94 

95 

96class MultipleResponseInformation(StationXMLError): 

97 pass 

98 

99 

100class InconsistentResponseInformation(StationXMLError): 

101 pass 

102 

103 

104class InconsistentChannelLocations(StationXMLError): 

105 pass 

106 

107 

108class InvalidRecord(StationXMLError): 

109 def __init__(self, line): 

110 StationXMLError.__init__(self) 

111 self._line = line 

112 

113 def __str__(self): 

114 return 'Invalid record: "%s"' % self._line 

115 

116 

117_exceptions = dict( 

118 Inconsistencies=Inconsistencies, 

119 NoResponseInformation=NoResponseInformation, 

120 MultipleResponseInformation=MultipleResponseInformation, 

121 InconsistentResponseInformation=InconsistentResponseInformation, 

122 InconsistentChannelLocations=InconsistentChannelLocations, 

123 InvalidRecord=InvalidRecord, 

124 ValueError=ValueError) 

125 

126 

127_logs = dict( 

128 info=logger.info, 

129 warning=logger.warning, 

130 error=logger.error) 

131 

132 

133class DeliveryError(StationXMLError): 

134 pass 

135 

136 

137class Delivery(Object): 

138 

139 def __init__(self, payload=None, log=None, errors=None, error=None): 

140 if payload is None: 

141 payload = [] 

142 

143 if log is None: 

144 log = [] 

145 

146 if errors is None: 

147 errors = [] 

148 

149 if error is not None: 

150 errors.append(error) 

151 

152 Object.__init__(self, payload=payload, log=log, errors=errors) 

153 

154 payload = List.T(Any.T()) 

155 log = List.T(Tuple.T(3, String.T())) 

156 errors = List.T(Tuple.T(3, String.T())) 

157 

158 def extend(self, other): 

159 self.payload.extend(other.payload) 

160 self.log.extend(other.log) 

161 self.errors.extend(other.errors) 

162 

163 def extend_without_payload(self, other): 

164 self.log.extend(other.log) 

165 self.errors.extend(other.errors) 

166 return other.payload 

167 

168 def emit_log(self): 

169 for name, message, context in self.log: 

170 message = '%s: %s' % (context, message) 

171 _logs[name](message) 

172 

173 def expect(self, quiet=False): 

174 if not quiet: 

175 self.emit_log() 

176 

177 if self.errors: 

178 name, message, context = self.errors[0] 

179 if context: 

180 message += ' (%s)' % context 

181 

182 if len(self.errors) > 1: 

183 message += ' Additional errors pending.' 

184 

185 raise _exceptions[name](message) 

186 

187 return self.payload 

188 

189 def expect_one(self, quiet=False): 

190 payload = self.expect(quiet=quiet) 

191 if len(payload) != 1: 

192 raise DeliveryError( 

193 'Expected 1 element but got %i.' % len(payload)) 

194 

195 return payload[0] 

196 

197 

198def wrap(s, width=80, indent=4): 

199 words = s.split() 

200 lines = [] 

201 t = [] 

202 n = 0 

203 for w in words: 

204 if n + len(w) >= width: 

205 lines.append(' '.join(t)) 

206 n = indent 

207 t = [' '*(indent-1)] 

208 

209 t.append(w) 

210 n += len(w) + 1 

211 

212 lines.append(' '.join(t)) 

213 return '\n'.join(lines) 

214 

215 

216def same(x, eps=0.0): 

217 if any(type(x[0]) != type(r) for r in x): 

218 return False 

219 

220 if isinstance(x[0], float): 

221 return all(abs(r-x[0]) <= eps for r in x) 

222 else: 

223 return all(r == x[0] for r in x) 

224 

225 

226def same_sample_rate(a, b, eps=1.0e-6): 

227 return abs(a - b) < min(a, b)*eps 

228 

229 

230def evaluate1(resp, f): 

231 return resp.evaluate(num.array([f], dtype=float))[0] 

232 

233 

234def check_resp(resp, value, frequency, limit_db, prelude, context): 

235 

236 try: 

237 value_resp = num.abs(evaluate1(resp, frequency)) 

238 except response.InvalidResponseError as e: 

239 return Delivery(log=[( 

240 'warning', 

241 'Could not check response: %s' % str(e), 

242 context)]) 

243 

244 if value_resp == 0.0: 

245 return Delivery(log=[( 

246 'warning', 

247 '%s\n' 

248 ' computed response is zero' % prelude, 

249 context)]) 

250 

251 diff_db = 20.0 * num.log10(value_resp/value) 

252 

253 if num.abs(diff_db) > limit_db: 

254 return Delivery(log=[( 

255 'warning', 

256 '%s\n' 

257 ' reported value: %g\n' 

258 ' computed value: %g\n' 

259 ' at frequency [Hz]: %g\n' 

260 ' factor: %g\n' 

261 ' difference [dB]: %g\n' 

262 ' limit [dB]: %g' % ( 

263 prelude, 

264 value, 

265 value_resp, 

266 frequency, 

267 value_resp/value, 

268 diff_db, 

269 limit_db), 

270 context)]) 

271 

272 return Delivery() 

273 

274 

275def tts(t): 

276 if t is None: 

277 return '<none>' 

278 else: 

279 return util.tts(t, format='%Y-%m-%d %H:%M:%S') 

280 

281 

282def le_open_left(a, b): 

283 return a is None or (b is not None and a <= b) 

284 

285 

286def le_open_right(a, b): 

287 return b is None or (a is not None and a <= b) 

288 

289 

290def eq_open(a, b): 

291 return (a is None and b is None) \ 

292 or (a is not None and b is not None and a == b) 

293 

294 

295def contains(a, b): 

296 return le_open_left(a.start_date, b.start_date) \ 

297 and le_open_right(b.end_date, a.end_date) 

298 

299 

300def find_containing(candidates, node): 

301 for candidate in candidates: 

302 if contains(candidate, node): 

303 return candidate 

304 

305 return None 

306 

307 

308this_year = time.gmtime()[0] 

309 

310 

311class DummyAwareOptionalTimestamp(Object): 

312 ''' 

313 Optional timestamp with support for some common placeholder values. 

314 

315 Some StationXML files contain arbitrary placeholder values for open end 

316 intervals, like "2100-01-01". Depending on the time range supported by the 

317 system, these dates are translated into ``None`` to prevent crashes with 

318 this type. 

319 ''' 

320 dummy_for = (hpfloat, float) 

321 dummy_for_description = 'pyrocko.util.get_time_float' 

322 

323 class __T(TBase): 

324 

325 def regularize_extra(self, val): 

326 time_float = get_time_float() 

327 

328 if isinstance(val, datetime.datetime): 

329 tt = val.utctimetuple() 

330 val = time_float(calendar.timegm(tt)) + val.microsecond * 1e-6 

331 

332 elif isinstance(val, datetime.date): 

333 tt = val.timetuple() 

334 val = time_float(calendar.timegm(tt)) 

335 

336 elif isinstance(val, str): 

337 val = val.strip() 

338 

339 tz_offset = 0 

340 

341 m = re_tz.search(val) 

342 if m: 

343 sh = m.group(2) 

344 sm = m.group(4) 

345 tz_offset = (int(sh)*3600 if sh else 0) \ 

346 + (int(sm)*60 if sm else 0) 

347 

348 val = re_tz.sub('', val) 

349 

350 if len(val) > 10 and val[10] == 'T': 

351 val = val.replace('T', ' ', 1) 

352 

353 try: 

354 val = util.str_to_time(val) - tz_offset 

355 

356 except util.TimeStrError: 

357 year = int(val[:4]) 

358 if sys.maxsize > 2**32: # if we're on 64bit 

359 if year > this_year + 100: 

360 return None # StationXML contained a dummy date 

361 

362 if year < 1903: # for macOS, 1900-01-01 dummy dates 

363 return None 

364 

365 else: # 32bit end of time is in 2038 

366 if this_year < 2037 and year > 2037 or year < 1903: 

367 return None # StationXML contained a dummy date 

368 

369 raise 

370 

371 elif isinstance(val, (int, float)): 

372 val = time_float(val) 

373 

374 else: 

375 raise ValidationError( 

376 '%s: cannot convert "%s" to type %s' % ( 

377 self.xname(), val, time_float)) 

378 

379 return val 

380 

381 def to_save(self, val): 

382 return time_to_str(val, format='%Y-%m-%d %H:%M:%S.9FRAC')\ 

383 .rstrip('0').rstrip('.') 

384 

385 def to_save_xml(self, val): 

386 return time_to_str(val, format='%Y-%m-%dT%H:%M:%S.9FRAC')\ 

387 .rstrip('0').rstrip('.') + 'Z' 

388 

389 

390class Nominal(StringChoice): 

391 choices = [ 

392 'NOMINAL', 

393 'CALCULATED'] 

394 

395 

396class Email(UnicodePattern): 

397 pattern = u'[\\w\\.\\-_]+@[\\w\\.\\-_]+' 

398 

399 

400class RestrictedStatus(StringChoice): 

401 choices = [ 

402 'open', 

403 'closed', 

404 'partial'] 

405 

406 

407class Type(StringChoice): 

408 choices = [ 

409 'TRIGGERED', 

410 'CONTINUOUS', 

411 'HEALTH', 

412 'GEOPHYSICAL', 

413 'WEATHER', 

414 'FLAG', 

415 'SYNTHESIZED', 

416 'INPUT', 

417 'EXPERIMENTAL', 

418 'MAINTENANCE', 

419 'BEAM'] 

420 

421 class __T(StringChoice.T): 

422 def validate_extra(self, val): 

423 if val not in self.choices: 

424 logger.warning( 

425 'channel type: "%s" is not a valid choice out of %s' % 

426 (val, repr(self.choices))) 

427 

428 

429class PzTransferFunction(StringChoice): 

430 choices = [ 

431 'LAPLACE (RADIANS/SECOND)', 

432 'LAPLACE (HERTZ)', 

433 'DIGITAL (Z-TRANSFORM)'] 

434 

435 

436class Symmetry(StringChoice): 

437 choices = [ 

438 'NONE', 

439 'EVEN', 

440 'ODD'] 

441 

442 

443class CfTransferFunction(StringChoice): 

444 

445 class __T(StringChoice.T): 

446 def validate(self, val, regularize=False, depth=-1): 

447 if regularize: 

448 try: 

449 val = str(val) 

450 except ValueError: 

451 raise ValidationError( 

452 '%s: cannot convert to string %s' % (self.xname, 

453 repr(val))) 

454 

455 val = self._dummy_cls.replacements.get(val, val) 

456 

457 self.validate_extra(val) 

458 return val 

459 

460 choices = [ 

461 'ANALOG (RADIANS/SECOND)', 

462 'ANALOG (HERTZ)', 

463 'DIGITAL'] 

464 

465 replacements = { 

466 'ANALOG (RAD/SEC)': 'ANALOG (RADIANS/SECOND)', 

467 'ANALOG (HZ)': 'ANALOG (HERTZ)', 

468 } 

469 

470 

471class Approximation(StringChoice): 

472 choices = [ 

473 'MACLAURIN'] 

474 

475 

476class PhoneNumber(StringPattern): 

477 pattern = '[0-9]+-[0-9]+' 

478 

479 

480class Site(Object): 

481 ''' 

482 Description of a site location using name and optional geopolitical 

483 boundaries (country, city, etc.). 

484 ''' 

485 

486 name = Unicode.T(default='', xmltagname='Name') 

487 description = Unicode.T(optional=True, xmltagname='Description') 

488 town = Unicode.T(optional=True, xmltagname='Town') 

489 county = Unicode.T(optional=True, xmltagname='County') 

490 region = Unicode.T(optional=True, xmltagname='Region') 

491 country = Unicode.T(optional=True, xmltagname='Country') 

492 

493 

494class ExternalReference(Object): 

495 ''' 

496 This type contains a URI and description for external data that users may 

497 want to reference in StationXML. 

498 ''' 

499 

500 uri = String.T(xmltagname='URI') 

501 description = Unicode.T(xmltagname='Description') 

502 

503 

504class Units(Object): 

505 ''' 

506 A type to document units. Corresponds to SEED blockette 34. 

507 ''' 

508 

509 def __init__(self, name=None, **kwargs): 

510 Object.__init__(self, name=name, **kwargs) 

511 

512 name = String.T(xmltagname='Name') 

513 description = Unicode.T(optional=True, xmltagname='Description') 

514 

515 

516class Counter(Int): 

517 pass 

518 

519 

520class SampleRateRatio(Object): 

521 ''' 

522 Sample rate expressed as number of samples in a number of seconds. 

523 ''' 

524 

525 number_samples = Int.T(xmltagname='NumberSamples') 

526 number_seconds = Int.T(xmltagname='NumberSeconds') 

527 

528 

529class Gain(Object): 

530 ''' 

531 Complex type for sensitivity and frequency ranges. This complex type can be 

532 used to represent both overall sensitivities and individual stage gains. 

533 The FrequencyRangeGroup is an optional construct that defines a pass band 

534 in Hertz ( FrequencyStart and FrequencyEnd) in which the SensitivityValue 

535 is valid within the number of decibels specified in FrequencyDBVariation. 

536 ''' 

537 

538 def __init__(self, value=None, **kwargs): 

539 Object.__init__(self, value=value, **kwargs) 

540 

541 value = Float.T(optional=True, xmltagname='Value') 

542 frequency = Float.T(optional=True, xmltagname='Frequency') 

543 

544 def summary(self): 

545 return 'gain(%g @ %g)' % (self.value, self.frequency) 

546 

547 

548class NumeratorCoefficient(Object): 

549 i = Int.T(optional=True, xmlstyle='attribute') 

550 value = Float.T(xmlstyle='content') 

551 

552 

553class FloatNoUnit(Object): 

554 def __init__(self, value=None, **kwargs): 

555 Object.__init__(self, value=value, **kwargs) 

556 

557 plus_error = Float.T(optional=True, xmlstyle='attribute') 

558 minus_error = Float.T(optional=True, xmlstyle='attribute') 

559 value = Float.T(xmlstyle='content') 

560 

561 

562class FloatWithUnit(FloatNoUnit): 

563 unit = String.T(optional=True, xmlstyle='attribute') 

564 

565 

566class Equipment(Object): 

567 resource_id = String.T(optional=True, xmlstyle='attribute') 

568 type = String.T(optional=True, xmltagname='Type') 

569 description = Unicode.T(optional=True, xmltagname='Description') 

570 manufacturer = Unicode.T(optional=True, xmltagname='Manufacturer') 

571 vendor = Unicode.T(optional=True, xmltagname='Vendor') 

572 model = Unicode.T(optional=True, xmltagname='Model') 

573 serial_number = String.T(optional=True, xmltagname='SerialNumber') 

574 installation_date = DummyAwareOptionalTimestamp.T( 

575 optional=True, 

576 xmltagname='InstallationDate') 

577 removal_date = DummyAwareOptionalTimestamp.T( 

578 optional=True, 

579 xmltagname='RemovalDate') 

580 calibration_date_list = List.T(Timestamp.T(xmltagname='CalibrationDate')) 

581 

582 

583class PhoneNumber(Object): 

584 description = Unicode.T(optional=True, xmlstyle='attribute') 

585 country_code = Int.T(optional=True, xmltagname='CountryCode') 

586 area_code = Int.T(xmltagname='AreaCode') 

587 phone_number = PhoneNumber.T(xmltagname='PhoneNumber') 

588 

589 

590class BaseFilter(Object): 

591 ''' 

592 The BaseFilter is derived by all filters. 

593 ''' 

594 

595 resource_id = String.T(optional=True, xmlstyle='attribute') 

596 name = String.T(optional=True, xmlstyle='attribute') 

597 description = Unicode.T(optional=True, xmltagname='Description') 

598 input_units = Units.T(optional=True, xmltagname='InputUnits') 

599 output_units = Units.T(optional=True, xmltagname='OutputUnits') 

600 

601 

602class Sensitivity(Gain): 

603 ''' 

604 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional 

605 construct that defines a pass band in Hertz (FrequencyStart and 

606 FrequencyEnd) in which the SensitivityValue is valid within the number of 

607 decibels specified in FrequencyDBVariation. 

608 ''' 

609 

610 input_units = Units.T(optional=True, xmltagname='InputUnits') 

611 output_units = Units.T(optional=True, xmltagname='OutputUnits') 

612 frequency_start = Float.T(optional=True, xmltagname='FrequencyStart') 

613 frequency_end = Float.T(optional=True, xmltagname='FrequencyEnd') 

614 frequency_db_variation = Float.T(optional=True, 

615 xmltagname='FrequencyDBVariation') 

616 

617 def get_pyrocko_response(self): 

618 return Delivery( 

619 [response.PoleZeroResponse(constant=self.value)]) 

620 

621 

622class Coefficient(FloatNoUnit): 

623 number = Counter.T(optional=True, xmlstyle='attribute') 

624 

625 

626class PoleZero(Object): 

627 ''' 

628 Complex numbers used as poles or zeros in channel response. 

629 ''' 

630 

631 number = Int.T(optional=True, xmlstyle='attribute') 

632 real = FloatNoUnit.T(xmltagname='Real') 

633 imaginary = FloatNoUnit.T(xmltagname='Imaginary') 

634 

635 def value(self): 

636 return self.real.value + 1J * self.imaginary.value 

637 

638 

639class ClockDrift(FloatWithUnit): 

640 unit = String.T(default='SECONDS/SAMPLE', optional=True, 

641 xmlstyle='attribute') # fixed 

642 

643 

644class Second(FloatWithUnit): 

645 ''' 

646 A time value in seconds. 

647 ''' 

648 

649 unit = String.T(default='SECONDS', optional=True, xmlstyle='attribute') 

650 # fixed unit 

651 

652 

653class Voltage(FloatWithUnit): 

654 unit = String.T(default='VOLTS', optional=True, xmlstyle='attribute') 

655 # fixed unit 

656 

657 

658class Angle(FloatWithUnit): 

659 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute') 

660 # fixed unit 

661 

662 

663class Azimuth(FloatWithUnit): 

664 ''' 

665 Instrument azimuth, degrees clockwise from North. 

666 ''' 

667 

668 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute') 

669 # fixed unit 

670 

671 

672class Dip(FloatWithUnit): 

673 ''' 

674 Instrument dip in degrees down from horizontal. Together azimuth and dip 

675 describe the direction of the sensitive axis of the instrument. 

676 ''' 

677 

678 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute') 

679 # fixed unit 

680 

681 

682class Distance(FloatWithUnit): 

683 ''' 

684 Extension of FloatWithUnit for distances, elevations, and depths. 

685 ''' 

686 

687 unit = String.T(default='METERS', optional=True, xmlstyle='attribute') 

688 # NOT fixed unit! 

689 

690 

691class Frequency(FloatWithUnit): 

692 unit = String.T(default='HERTZ', optional=True, xmlstyle='attribute') 

693 # fixed unit 

694 

695 

696class SampleRate(FloatWithUnit): 

697 ''' 

698 Sample rate in samples per second. 

699 ''' 

700 

701 unit = String.T(default='SAMPLES/S', optional=True, xmlstyle='attribute') 

702 # fixed unit 

703 

704 

705class Person(Object): 

706 ''' 

707 Representation of a person's contact information. A person can belong to 

708 multiple agencies and have multiple email addresses and phone numbers. 

709 ''' 

710 

711 name_list = List.T(Unicode.T(xmltagname='Name')) 

712 agency_list = List.T(Unicode.T(xmltagname='Agency')) 

713 email_list = List.T(Email.T(xmltagname='Email')) 

714 phone_list = List.T(PhoneNumber.T(xmltagname='Phone')) 

715 

716 

717class FIR(BaseFilter): 

718 ''' 

719 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are 

720 also commonly documented using the Coefficients element. 

721 ''' 

722 

723 symmetry = Symmetry.T(xmltagname='Symmetry') 

724 numerator_coefficient_list = List.T( 

725 NumeratorCoefficient.T(xmltagname='NumeratorCoefficient')) 

726 

727 def summary(self): 

728 return 'fir(%i%s)' % ( 

729 self.get_ncoefficiencs(), 

730 ',sym' if self.get_effective_symmetry() != 'NONE' else '') 

731 

732 def get_effective_coefficients(self): 

733 b = num.array( 

734 [v.value for v in self.numerator_coefficient_list], 

735 dtype=float) 

736 

737 if self.symmetry == 'ODD': 

738 b = num.concatenate((b, b[-2::-1])) 

739 elif self.symmetry == 'EVEN': 

740 b = num.concatenate((b, b[::-1])) 

741 

742 return b 

743 

744 def get_effective_symmetry(self): 

745 if self.symmetry == 'NONE': 

746 b = self.get_effective_coefficients() 

747 if num.all(b - b[::-1] == 0): 

748 return ['EVEN', 'ODD'][b.size % 2] 

749 

750 return self.symmetry 

751 

752 def get_ncoefficiencs(self): 

753 nf = len(self.numerator_coefficient_list) 

754 if self.symmetry == 'ODD': 

755 nc = nf * 2 + 1 

756 elif self.symmetry == 'EVEN': 

757 nc = nf * 2 

758 else: 

759 nc = nf 

760 

761 return nc 

762 

763 def estimate_delay(self, deltat): 

764 nc = self.get_ncoefficiencs() 

765 if nc > 0: 

766 return deltat * (nc - 1) / 2.0 

767 else: 

768 return 0.0 

769 

770 def get_pyrocko_response( 

771 self, context, deltat, delay_responses, normalization_frequency): 

772 

773 context += self.summary() 

774 

775 if not self.numerator_coefficient_list: 

776 return Delivery([]) 

777 

778 b = self.get_effective_coefficients() 

779 

780 log = [] 

781 drop_phase = self.get_effective_symmetry() != 'NONE' 

782 

783 if not deltat: 

784 log.append(( 

785 'error', 

786 'Digital response requires knowledge about sampling ' 

787 'interval. Response will be unusable.', 

788 context)) 

789 

790 resp = response.DigitalFilterResponse( 

791 b.tolist(), [1.0], deltat or 0.0, drop_phase=drop_phase) 

792 

793 if normalization_frequency is not None and deltat is not None: 

794 normalization_frequency = 0.0 

795 normalization = num.abs(evaluate1(resp, normalization_frequency)) 

796 

797 if num.abs(normalization - 1.0) > 1e-2: 

798 log.append(( 

799 'warning', 

800 'FIR filter coefficients are not normalized. Normalizing ' 

801 'them. Factor: %g' % normalization, context)) 

802 

803 resp = response.DigitalFilterResponse( 

804 (b/normalization).tolist(), [1.0], deltat, 

805 drop_phase=drop_phase) 

806 

807 resps = [resp] 

808 

809 if not drop_phase: 

810 resps.extend(delay_responses) 

811 

812 return Delivery(resps, log=log) 

813 

814 

815class Coefficients(BaseFilter): 

816 ''' 

817 Response: coefficients for FIR filter. Laplace transforms or IIR filters 

818 can be expressed using type as well but the PolesAndZeros should be used 

819 instead. Corresponds to SEED blockette 54. 

820 ''' 

821 

822 cf_transfer_function_type = CfTransferFunction.T( 

823 xmltagname='CfTransferFunctionType') 

824 numerator_list = List.T(FloatWithUnit.T(xmltagname='Numerator')) 

825 denominator_list = List.T(FloatWithUnit.T(xmltagname='Denominator')) 

826 

827 def summary(self): 

828 return 'coef_%s(%i,%i%s)' % ( 

829 'ABC?'[ 

830 CfTransferFunction.choices.index( 

831 self.cf_transfer_function_type)], 

832 len(self.numerator_list), 

833 len(self.denominator_list), 

834 ',sym' if self.is_symmetric_fir else '') 

835 

836 def estimate_delay(self, deltat): 

837 nc = len(self.numerator_list) 

838 if nc > 0: 

839 return deltat * (len(self.numerator_list) - 1) / 2.0 

840 else: 

841 return 0.0 

842 

843 def is_symmetric_fir(self): 

844 if len(self.denominator_list) != 0: 

845 return False 

846 b = [v.value for v in self.numerator_list] 

847 return b == b[::-1] 

848 

849 def get_pyrocko_response( 

850 self, context, deltat, delay_responses, normalization_frequency): 

851 

852 context += self.summary() 

853 

854 factor = 1.0 

855 if self.cf_transfer_function_type == 'ANALOG (HERTZ)': 

856 factor = 2.0*math.pi 

857 

858 if not self.numerator_list and not self.denominator_list: 

859 return Delivery(payload=[]) 

860 

861 b = num.array( 

862 [v.value*factor for v in self.numerator_list], dtype=float) 

863 

864 a = num.array( 

865 [1.0] + [v.value*factor for v in self.denominator_list], 

866 dtype=float) 

867 

868 log = [] 

869 resps = [] 

870 if self.cf_transfer_function_type in [ 

871 'ANALOG (RADIANS/SECOND)', 'ANALOG (HERTZ)']: 

872 resps.append(response.AnalogFilterResponse(b, a)) 

873 

874 elif self.cf_transfer_function_type == 'DIGITAL': 

875 if not deltat: 

876 log.append(( 

877 'error', 

878 'Digital response requires knowledge about sampling ' 

879 'interval. Response will be unusable.', 

880 context)) 

881 

882 drop_phase = self.is_symmetric_fir() 

883 resp = response.DigitalFilterResponse( 

884 b, a, deltat or 0.0, drop_phase=drop_phase) 

885 

886 if normalization_frequency is not None and deltat is not None: 

887 normalization = num.abs(evaluate1( 

888 resp, normalization_frequency)) 

889 

890 if num.abs(normalization - 1.0) > 1e-2: 

891 log.append(( 

892 'warning', 

893 'FIR filter coefficients are not normalized. ' 

894 'Normalizing them. Factor: %g' % normalization, 

895 context)) 

896 

897 resp = response.DigitalFilterResponse( 

898 (b/normalization).tolist(), [1.0], deltat, 

899 drop_phase=drop_phase) 

900 

901 resps.append(resp) 

902 

903 if not drop_phase: 

904 resps.extend(delay_responses) 

905 

906 else: 

907 return Delivery(error=( 

908 'ValueError', 

909 'Unknown transfer function type: %s' % ( 

910 self.cf_transfer_function_type))) 

911 

912 return Delivery(payload=resps, log=log) 

913 

914 

915class Latitude(FloatWithUnit): 

916 ''' 

917 Type for latitude coordinate. 

918 ''' 

919 

920 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute') 

921 # fixed unit 

922 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute') 

923 

924 

925class Longitude(FloatWithUnit): 

926 ''' 

927 Type for longitude coordinate. 

928 ''' 

929 

930 unit = String.T(default='DEGREES', optional=True, xmlstyle='attribute') 

931 # fixed unit 

932 datum = String.T(default='WGS84', optional=True, xmlstyle='attribute') 

933 

934 

935class PolesZeros(BaseFilter): 

936 ''' 

937 Response: complex poles and zeros. Corresponds to SEED blockette 53. 

938 ''' 

939 

940 pz_transfer_function_type = PzTransferFunction.T( 

941 xmltagname='PzTransferFunctionType') 

942 normalization_factor = Float.T(default=1.0, 

943 xmltagname='NormalizationFactor') 

944 normalization_frequency = Frequency.T(xmltagname='NormalizationFrequency') 

945 zero_list = List.T(PoleZero.T(xmltagname='Zero')) 

946 pole_list = List.T(PoleZero.T(xmltagname='Pole')) 

947 

948 def summary(self): 

949 return 'pz_%s(%i,%i)' % ( 

950 'ABC?'[ 

951 PzTransferFunction.choices.index( 

952 self.pz_transfer_function_type)], 

953 len(self.pole_list), 

954 len(self.zero_list)) 

955 

956 def get_pyrocko_response(self, context='', deltat=None): 

957 

958 context += self.summary() 

959 

960 factor = 1.0 

961 cfactor = 1.0 

962 if self.pz_transfer_function_type == 'LAPLACE (HERTZ)': 

963 factor = 2. * math.pi 

964 cfactor = (2. * math.pi)**( 

965 len(self.pole_list) - len(self.zero_list)) 

966 

967 log = [] 

968 if self.normalization_factor is None \ 

969 or self.normalization_factor == 0.0: 

970 

971 log.append(( 

972 'warning', 

973 'No pole-zero normalization factor given. ' 

974 'Assuming a value of 1.0', 

975 context)) 

976 

977 nfactor = 1.0 

978 else: 

979 nfactor = self.normalization_factor 

980 

981 is_digital = self.pz_transfer_function_type == 'DIGITAL (Z-TRANSFORM)' 

982 if not is_digital: 

983 resp = response.PoleZeroResponse( 

984 constant=nfactor*cfactor, 

985 zeros=[z.value()*factor for z in self.zero_list], 

986 poles=[p.value()*factor for p in self.pole_list]) 

987 else: 

988 if not deltat: 

989 log.append(( 

990 'error', 

991 'Digital response requires knowledge about sampling ' 

992 'interval. Response will be unusable.', 

993 context)) 

994 

995 resp = response.DigitalPoleZeroResponse( 

996 constant=nfactor*cfactor, 

997 zeros=[z.value()*factor for z in self.zero_list], 

998 poles=[p.value()*factor for p in self.pole_list], 

999 deltat=deltat or 0.0) 

1000 

1001 if not self.normalization_frequency.value: 

1002 log.append(( 

1003 'warning', 

1004 'Cannot check pole-zero normalization factor: ' 

1005 'No normalization frequency given.', 

1006 context)) 

1007 

1008 else: 

1009 if is_digital and not deltat: 

1010 log.append(( 

1011 'warning', 

1012 'Cannot check computed vs reported normalization ' 

1013 'factor without knowing the sampling interval.', 

1014 context)) 

1015 else: 

1016 computed_normalization_factor = nfactor / abs(evaluate1( 

1017 resp, self.normalization_frequency.value)) 

1018 

1019 db = 20.0 * num.log10( 

1020 computed_normalization_factor / nfactor) 

1021 

1022 if abs(db) > 0.17: 

1023 log.append(( 

1024 'warning', 

1025 'Computed and reported normalization factors differ ' 

1026 'by %g dB: computed: %g, reported: %g' % ( 

1027 db, 

1028 computed_normalization_factor, 

1029 nfactor), 

1030 context)) 

1031 

1032 return Delivery([resp], log) 

1033 

1034 

1035class ResponseListElement(Object): 

1036 frequency = Frequency.T(xmltagname='Frequency') 

1037 amplitude = FloatWithUnit.T(xmltagname='Amplitude') 

1038 phase = Angle.T(xmltagname='Phase') 

1039 

1040 

1041class Polynomial(BaseFilter): 

1042 ''' 

1043 Response: expressed as a polynomial (allows non-linear sensors to be 

1044 described). Corresponds to SEED blockette 62. Can be used to describe a 

1045 stage of acquisition or a complete system. 

1046 ''' 

1047 

1048 approximation_type = Approximation.T(default='MACLAURIN', 

1049 xmltagname='ApproximationType') 

1050 frequency_lower_bound = Frequency.T(xmltagname='FrequencyLowerBound') 

1051 frequency_upper_bound = Frequency.T(xmltagname='FrequencyUpperBound') 

1052 approximation_lower_bound = Float.T(xmltagname='ApproximationLowerBound') 

1053 approximation_upper_bound = Float.T(xmltagname='ApproximationUpperBound') 

1054 maximum_error = Float.T(xmltagname='MaximumError') 

1055 coefficient_list = List.T(Coefficient.T(xmltagname='Coefficient')) 

1056 

1057 def summary(self): 

1058 return 'poly(%i)' % len(self.coefficient_list) 

1059 

1060 

1061class Decimation(Object): 

1062 ''' 

1063 Corresponds to SEED blockette 57. 

1064 ''' 

1065 

1066 input_sample_rate = Frequency.T(xmltagname='InputSampleRate') 

1067 factor = Int.T(xmltagname='Factor') 

1068 offset = Int.T(xmltagname='Offset') 

1069 delay = FloatWithUnit.T(xmltagname='Delay') 

1070 correction = FloatWithUnit.T(xmltagname='Correction') 

1071 

1072 def summary(self): 

1073 return 'deci(%i, %g -> %g, %g)' % ( 

1074 self.factor, 

1075 self.input_sample_rate.value, 

1076 self.input_sample_rate.value / self.factor, 

1077 self.delay.value) 

1078 

1079 def get_pyrocko_response(self): 

1080 if self.delay and self.delay.value != 0.0: 

1081 return Delivery([response.DelayResponse(delay=-self.delay.value)]) 

1082 else: 

1083 return Delivery([]) 

1084 

1085 

1086class Operator(Object): 

1087 agency_list = List.T(Unicode.T(xmltagname='Agency')) 

1088 contact_list = List.T(Person.T(xmltagname='Contact')) 

1089 web_site = String.T(optional=True, xmltagname='WebSite') 

1090 

1091 

1092class Comment(Object): 

1093 ''' 

1094 Container for a comment or log entry. Corresponds to SEED blockettes 31, 51 

1095 and 59. 

1096 ''' 

1097 

1098 id = Counter.T(optional=True, xmlstyle='attribute') 

1099 value = Unicode.T(xmltagname='Value') 

1100 begin_effective_time = DummyAwareOptionalTimestamp.T( 

1101 optional=True, 

1102 xmltagname='BeginEffectiveTime') 

1103 end_effective_time = DummyAwareOptionalTimestamp.T( 

1104 optional=True, 

1105 xmltagname='EndEffectiveTime') 

1106 author_list = List.T(Person.T(xmltagname='Author')) 

1107 

1108 

1109class ResponseList(BaseFilter): 

1110 ''' 

1111 Response: list of frequency, amplitude and phase values. Corresponds to 

1112 SEED blockette 55. 

1113 ''' 

1114 

1115 response_list_element_list = List.T( 

1116 ResponseListElement.T(xmltagname='ResponseListElement')) 

1117 

1118 def summary(self): 

1119 return 'list(%i)' % len(self.response_list_element_list) 

1120 

1121 

1122class Log(Object): 

1123 ''' 

1124 Container for log entries. 

1125 ''' 

1126 

1127 entry_list = List.T(Comment.T(xmltagname='Entry')) 

1128 

1129 

1130class ResponseStage(Object): 

1131 ''' 

1132 This complex type represents channel response and covers SEED blockettes 53 

1133 to 56. 

1134 ''' 

1135 

1136 number = Counter.T(xmlstyle='attribute') 

1137 resource_id = String.T(optional=True, xmlstyle='attribute') 

1138 poles_zeros_list = List.T( 

1139 PolesZeros.T(optional=True, xmltagname='PolesZeros')) 

1140 coefficients_list = List.T( 

1141 Coefficients.T(optional=True, xmltagname='Coefficients')) 

1142 response_list = ResponseList.T(optional=True, xmltagname='ResponseList') 

1143 fir = FIR.T(optional=True, xmltagname='FIR') 

1144 polynomial = Polynomial.T(optional=True, xmltagname='Polynomial') 

1145 decimation = Decimation.T(optional=True, xmltagname='Decimation') 

1146 stage_gain = Gain.T(optional=True, xmltagname='StageGain') 

1147 

1148 def summary(self): 

1149 elements = [] 

1150 

1151 for stuff in [ 

1152 self.poles_zeros_list, self.coefficients_list, 

1153 self.response_list, self.fir, self.polynomial, 

1154 self.decimation, self.stage_gain]: 

1155 

1156 if stuff: 

1157 if isinstance(stuff, Object): 

1158 elements.append(stuff.summary()) 

1159 else: 

1160 elements.extend(obj.summary() for obj in stuff) 

1161 

1162 return '%i: %s %s -> %s' % ( 

1163 self.number, 

1164 ', '.join(elements), 

1165 self.input_units.name.upper() if self.input_units else '?', 

1166 self.output_units.name.upper() if self.output_units else '?') 

1167 

1168 def get_squirrel_response_stage(self, context): 

1169 from pyrocko.squirrel.model import ResponseStage 

1170 

1171 delivery = Delivery() 

1172 delivery_pr = self.get_pyrocko_response(context) 

1173 log = delivery_pr.log 

1174 delivery_pr.log = [] 

1175 elements = delivery.extend_without_payload(delivery_pr) 

1176 

1177 delivery.payload = [ResponseStage( 

1178 input_quantity=to_quantity(self.input_units, context, delivery), 

1179 output_quantity=to_quantity(self.output_units, context, delivery), 

1180 input_sample_rate=self.input_sample_rate, 

1181 output_sample_rate=self.output_sample_rate, 

1182 elements=elements, 

1183 log=log)] 

1184 

1185 return delivery 

1186 

1187 def get_pyrocko_response(self, context, gain_only=False): 

1188 

1189 context = context + ', stage %i' % self.number 

1190 

1191 responses = [] 

1192 log = [] 

1193 if self.stage_gain: 

1194 normalization_frequency = self.stage_gain.frequency or 0.0 

1195 else: 

1196 normalization_frequency = 0.0 

1197 

1198 if not gain_only: 

1199 deltat = None 

1200 delay_responses = [] 

1201 if self.decimation: 

1202 rate = self.decimation.input_sample_rate.value 

1203 if rate > 0.0: 

1204 deltat = 1.0 / rate 

1205 delivery = self.decimation.get_pyrocko_response() 

1206 if delivery.errors: 

1207 return delivery 

1208 

1209 delay_responses = delivery.payload 

1210 log.extend(delivery.log) 

1211 

1212 for pzs in self.poles_zeros_list: 

1213 delivery = pzs.get_pyrocko_response(context, deltat) 

1214 if delivery.errors: 

1215 return delivery 

1216 

1217 pz_resps = delivery.payload 

1218 log.extend(delivery.log) 

1219 responses.extend(pz_resps) 

1220 

1221 # emulate incorrect? evalresp behaviour 

1222 if pzs.normalization_frequency != normalization_frequency \ 

1223 and normalization_frequency != 0.0: 

1224 

1225 try: 

1226 trial = response.MultiplyResponse(pz_resps) 

1227 anorm = num.abs(evaluate1( 

1228 trial, pzs.normalization_frequency.value)) 

1229 asens = num.abs( 

1230 evaluate1(trial, normalization_frequency)) 

1231 

1232 factor = anorm/asens 

1233 

1234 if abs(factor - 1.0) > 0.01: 

1235 log.append(( 

1236 'warning', 

1237 'PZ normalization frequency (%g) is different ' 

1238 'from stage gain frequency (%s) -> Emulating ' 

1239 'possibly incorrect evalresp behaviour. ' 

1240 'Correction factor: %g' % ( 

1241 pzs.normalization_frequency.value, 

1242 normalization_frequency, 

1243 factor), 

1244 context)) 

1245 

1246 responses.append( 

1247 response.PoleZeroResponse(constant=factor)) 

1248 except response.InvalidResponseError as e: 

1249 log.append(( 

1250 'warning', 

1251 'Could not check response: %s' % str(e), 

1252 context)) 

1253 

1254 if len(self.poles_zeros_list) > 1: 

1255 log.append(( 

1256 'warning', 

1257 'Multiple poles and zeros records in single response ' 

1258 'stage.', 

1259 context)) 

1260 

1261 for cfs in self.coefficients_list + ( 

1262 [self.fir] if self.fir else []): 

1263 

1264 delivery = cfs.get_pyrocko_response( 

1265 context, deltat, delay_responses, 

1266 normalization_frequency) 

1267 

1268 if delivery.errors: 

1269 return delivery 

1270 

1271 responses.extend(delivery.payload) 

1272 log.extend(delivery.log) 

1273 

1274 if len(self.coefficients_list) > 1: 

1275 log.append(( 

1276 'warning', 

1277 'Multiple filter coefficients lists in single response ' 

1278 'stage.', 

1279 context)) 

1280 

1281 if self.response_list: 

1282 log.append(( 

1283 'warning', 

1284 'Unhandled response element of type: ResponseList', 

1285 context)) 

1286 

1287 if self.polynomial: 

1288 log.append(( 

1289 'warning', 

1290 'Unhandled response element of type: Polynomial', 

1291 context)) 

1292 

1293 if self.stage_gain: 

1294 responses.append( 

1295 response.PoleZeroResponse(constant=self.stage_gain.value)) 

1296 

1297 return Delivery(responses, log) 

1298 

1299 @property 

1300 def input_units(self): 

1301 for e in (self.poles_zeros_list + self.coefficients_list + 

1302 [self.response_list, self.fir, self.polynomial]): 

1303 if e is not None: 

1304 return e.input_units 

1305 

1306 return None 

1307 

1308 @property 

1309 def output_units(self): 

1310 for e in (self.poles_zeros_list + self.coefficients_list + 

1311 [self.response_list, self.fir, self.polynomial]): 

1312 if e is not None: 

1313 return e.output_units 

1314 

1315 return None 

1316 

1317 @property 

1318 def input_sample_rate(self): 

1319 if self.decimation: 

1320 return self.decimation.input_sample_rate.value 

1321 

1322 return None 

1323 

1324 @property 

1325 def output_sample_rate(self): 

1326 if self.decimation: 

1327 return self.decimation.input_sample_rate.value \ 

1328 / self.decimation.factor 

1329 

1330 return None 

1331 

1332 

1333class Response(Object): 

1334 resource_id = String.T(optional=True, xmlstyle='attribute') 

1335 instrument_sensitivity = Sensitivity.T(optional=True, 

1336 xmltagname='InstrumentSensitivity') 

1337 instrument_polynomial = Polynomial.T(optional=True, 

1338 xmltagname='InstrumentPolynomial') 

1339 stage_list = List.T(ResponseStage.T(xmltagname='Stage')) 

1340 

1341 @property 

1342 def output_sample_rate(self): 

1343 if self.stage_list: 

1344 return self.stage_list[-1].output_sample_rate 

1345 else: 

1346 return None 

1347 

1348 def check_sample_rates(self, channel): 

1349 

1350 if self.stage_list: 

1351 sample_rate = None 

1352 

1353 for stage in self.stage_list: 

1354 if stage.decimation: 

1355 input_sample_rate = \ 

1356 stage.decimation.input_sample_rate.value 

1357 

1358 if sample_rate is not None and not same_sample_rate( 

1359 sample_rate, input_sample_rate): 

1360 

1361 logger.warning( 

1362 'Response stage %i has unexpected input sample ' 

1363 'rate: %g Hz (expected: %g Hz)' % ( 

1364 stage.number, 

1365 input_sample_rate, 

1366 sample_rate)) 

1367 

1368 sample_rate = input_sample_rate / stage.decimation.factor 

1369 

1370 if sample_rate is not None and channel.sample_rate \ 

1371 and not same_sample_rate( 

1372 sample_rate, channel.sample_rate.value): 

1373 

1374 logger.warning( 

1375 'Channel sample rate (%g Hz) does not match sample rate ' 

1376 'deducted from response stages information (%g Hz).' % ( 

1377 channel.sample_rate.value, 

1378 sample_rate)) 

1379 

1380 def check_units(self): 

1381 

1382 if self.instrument_sensitivity \ 

1383 and self.instrument_sensitivity.input_units: 

1384 

1385 units = self.instrument_sensitivity.input_units.name.upper() 

1386 

1387 if self.stage_list: 

1388 for stage in self.stage_list: 

1389 if units and stage.input_units \ 

1390 and stage.input_units.name.upper() != units: 

1391 

1392 logger.warning( 

1393 'Input units of stage %i (%s) do not match %s (%s).' 

1394 % ( 

1395 stage.number, 

1396 units, 

1397 'output units of stage %i' 

1398 if stage.number == 0 

1399 else 'sensitivity input units', 

1400 units)) 

1401 

1402 if stage.output_units: 

1403 units = stage.output_units.name.upper() 

1404 else: 

1405 units = None 

1406 

1407 sout_units = self.instrument_sensitivity.output_units 

1408 if self.instrument_sensitivity and sout_units: 

1409 if units is not None and units != sout_units.name.upper(): 

1410 logger.warning( 

1411 'Output units of stage %i (%s) do not match %s (%s).' 

1412 % ( 

1413 stage.number, 

1414 units, 

1415 'sensitivity output units', 

1416 sout_units.name.upper())) 

1417 

1418 def _sensitivity_checkpoints(self, responses, context): 

1419 delivery = Delivery() 

1420 

1421 if self.instrument_sensitivity: 

1422 sval = self.instrument_sensitivity.value 

1423 sfreq = self.instrument_sensitivity.frequency 

1424 if sval is None: 

1425 delivery.log.append(( 

1426 'warning', 

1427 'No sensitivity value given.', 

1428 context)) 

1429 

1430 elif sval is None: 

1431 delivery.log.append(( 

1432 'warning', 

1433 'Reported sensitivity value is zero.', 

1434 context)) 

1435 

1436 elif sfreq is None: 

1437 delivery.log.append(( 

1438 'warning', 

1439 'Sensitivity frequency not given.', 

1440 context)) 

1441 

1442 else: 

1443 trial = response.MultiplyResponse(responses) 

1444 

1445 delivery.extend( 

1446 check_resp( 

1447 trial, sval, sfreq, 0.1, 

1448 'Instrument sensitivity value inconsistent with ' 

1449 'sensitivity computed from complete response.', 

1450 context)) 

1451 

1452 delivery.payload.append(response.FrequencyResponseCheckpoint( 

1453 frequency=sfreq, value=sval)) 

1454 

1455 return delivery 

1456 

1457 def get_squirrel_response(self, context, **kwargs): 

1458 from pyrocko.squirrel.model import Response 

1459 

1460 if self.stage_list: 

1461 delivery = Delivery() 

1462 for istage, stage in enumerate(self.stage_list): 

1463 delivery.extend(stage.get_squirrel_response_stage(context)) 

1464 

1465 checkpoints = [] 

1466 if not delivery.errors: 

1467 all_responses = [] 

1468 for sq_stage in delivery.payload: 

1469 all_responses.extend(sq_stage.elements) 

1470 

1471 checkpoints.extend(delivery.extend_without_payload( 

1472 self._sensitivity_checkpoints(all_responses, context))) 

1473 

1474 sq_stages = delivery.payload 

1475 if sq_stages: 

1476 if sq_stages[0].input_quantity is None \ 

1477 and self.instrument_sensitivity is not None: 

1478 

1479 sq_stages[0].input_quantity = to_quantity( 

1480 self.instrument_sensitivity.input_units, 

1481 context, delivery) 

1482 sq_stages[-1].output_quantity = to_quantity( 

1483 self.instrument_sensitivity.output_units, 

1484 context, delivery) 

1485 

1486 sq_stages = delivery.expect() 

1487 

1488 return Response( 

1489 stages=sq_stages, 

1490 log=delivery.log, 

1491 checkpoints=checkpoints, 

1492 **kwargs) 

1493 

1494 elif self.instrument_sensitivity: 

1495 raise NoResponseInformation( 

1496 "Only instrument sensitivity given (won't use it). (%s)." 

1497 % context) 

1498 else: 

1499 raise NoResponseInformation( 

1500 'Empty instrument response (%s).' 

1501 % context) 

1502 

1503 def get_pyrocko_response( 

1504 self, context, fake_input_units=None, stages=(0, 1)): 

1505 

1506 delivery = Delivery() 

1507 if self.stage_list: 

1508 for istage, stage in enumerate(self.stage_list): 

1509 delivery.extend(stage.get_pyrocko_response( 

1510 context, gain_only=not ( 

1511 stages is None or stages[0] <= istage < stages[1]))) 

1512 

1513 elif self.instrument_sensitivity: 

1514 delivery.extend(self.instrument_sensitivity.get_pyrocko_response()) 

1515 

1516 delivery_cp = self._sensitivity_checkpoints(delivery.payload, context) 

1517 checkpoints = delivery.extend_without_payload(delivery_cp) 

1518 if delivery.errors: 

1519 return delivery 

1520 

1521 if fake_input_units is not None: 

1522 if not self.instrument_sensitivity or \ 

1523 self.instrument_sensitivity.input_units is None: 

1524 

1525 delivery.errors.append(( 

1526 'NoResponseInformation', 

1527 'No input units given, so cannot convert to requested ' 

1528 'units: %s' % fake_input_units.upper(), 

1529 context)) 

1530 

1531 return delivery 

1532 

1533 input_units = self.instrument_sensitivity.input_units.name.upper() 

1534 

1535 conresp = None 

1536 try: 

1537 conresp = conversion[ 

1538 fake_input_units.upper(), input_units] 

1539 

1540 except KeyError: 

1541 delivery.errors.append(( 

1542 'NoResponseInformation', 

1543 'Cannot convert between units: %s, %s' 

1544 % (fake_input_units, input_units), 

1545 context)) 

1546 

1547 if conresp is not None: 

1548 delivery.payload.append(conresp) 

1549 for checkpoint in checkpoints: 

1550 checkpoint.value *= num.abs(evaluate1( 

1551 conresp, checkpoint.frequency)) 

1552 

1553 delivery.payload = [ 

1554 response.MultiplyResponse( 

1555 delivery.payload, 

1556 checkpoints=checkpoints)] 

1557 

1558 return delivery 

1559 

1560 @classmethod 

1561 def from_pyrocko_pz_response(cls, presponse, input_unit, output_unit, 

1562 normalization_frequency=1.0): 

1563 ''' 

1564 Convert Pyrocko pole-zero response to StationXML response. 

1565 

1566 :param presponse: Pyrocko pole-zero response 

1567 :type presponse: :py:class:`~pyrocko.response.PoleZeroResponse` 

1568 :param input_unit: Input unit to be reported in the StationXML 

1569 response. 

1570 :type input_unit: str 

1571 :param output_unit: Output unit to be reported in the StationXML 

1572 response. 

1573 :type output_unit: str 

1574 :param normalization_frequency: Frequency where the normalization 

1575 factor for the StationXML response should be computed. 

1576 :type normalization_frequency: float 

1577 ''' 

1578 

1579 norm_factor = 1.0/float(abs( 

1580 evaluate1(presponse, normalization_frequency) 

1581 / presponse.constant)) 

1582 

1583 pzs = PolesZeros( 

1584 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)', 

1585 normalization_factor=norm_factor, 

1586 normalization_frequency=Frequency(normalization_frequency), 

1587 zero_list=[PoleZero(real=FloatNoUnit(z.real), 

1588 imaginary=FloatNoUnit(z.imag)) 

1589 for z in presponse.zeros], 

1590 pole_list=[PoleZero(real=FloatNoUnit(z.real), 

1591 imaginary=FloatNoUnit(z.imag)) 

1592 for z in presponse.poles]) 

1593 

1594 pzs.validate() 

1595 

1596 stage = ResponseStage( 

1597 number=1, 

1598 poles_zeros_list=[pzs], 

1599 stage_gain=Gain(float(abs(presponse.constant))/norm_factor)) 

1600 

1601 resp = Response( 

1602 instrument_sensitivity=Sensitivity( 

1603 value=stage.stage_gain.value, 

1604 frequency=normalization_frequency, 

1605 input_units=Units(input_unit), 

1606 output_units=Units(output_unit)), 

1607 

1608 stage_list=[stage]) 

1609 

1610 return resp 

1611 

1612 

1613class BaseNode(Object): 

1614 ''' 

1615 A base node type for derivation from: Network, Station and Channel types. 

1616 ''' 

1617 

1618 code = String.T(xmlstyle='attribute') 

1619 start_date = DummyAwareOptionalTimestamp.T(optional=True, 

1620 xmlstyle='attribute') 

1621 end_date = DummyAwareOptionalTimestamp.T(optional=True, 

1622 xmlstyle='attribute') 

1623 restricted_status = RestrictedStatus.T(optional=True, xmlstyle='attribute') 

1624 alternate_code = String.T(optional=True, xmlstyle='attribute') 

1625 historical_code = String.T(optional=True, xmlstyle='attribute') 

1626 description = Unicode.T(optional=True, xmltagname='Description') 

1627 comment_list = List.T(Comment.T(xmltagname='Comment')) 

1628 

1629 def spans(self, *args): 

1630 if len(args) == 0: 

1631 return True 

1632 elif len(args) == 1: 

1633 return ((self.start_date is None or 

1634 self.start_date <= args[0]) and 

1635 (self.end_date is None or 

1636 args[0] <= self.end_date)) 

1637 

1638 elif len(args) == 2: 

1639 return ((self.start_date is None or 

1640 args[1] >= self.start_date) and 

1641 (self.end_date is None or 

1642 self.end_date >= args[0])) 

1643 

1644 

1645def overlaps(a, b): 

1646 return ( 

1647 a.start_date is None or b.end_date is None 

1648 or a.start_date < b.end_date 

1649 ) and ( 

1650 b.start_date is None or a.end_date is None 

1651 or b.start_date < a.end_date 

1652 ) 

1653 

1654 

1655def check_overlaps(node_type_name, codes, nodes): 

1656 errors = [] 

1657 for ia, a in enumerate(nodes): 

1658 for b in nodes[ia+1:]: 

1659 if overlaps(a, b): 

1660 errors.append( 

1661 '%s epochs overlap for %s:\n' 

1662 ' %s - %s\n %s - %s' % ( 

1663 node_type_name, 

1664 '.'.join(codes), 

1665 tts(a.start_date), tts(a.end_date), 

1666 tts(b.start_date), tts(b.end_date))) 

1667 

1668 return errors 

1669 

1670 

1671class Channel(BaseNode): 

1672 ''' 

1673 Equivalent to SEED blockette 52 and parent element for the related the 

1674 response blockettes. 

1675 ''' 

1676 

1677 location_code = String.T(xmlstyle='attribute') 

1678 external_reference_list = List.T( 

1679 ExternalReference.T(xmltagname='ExternalReference')) 

1680 latitude = Latitude.T(xmltagname='Latitude') 

1681 longitude = Longitude.T(xmltagname='Longitude') 

1682 elevation = Distance.T(xmltagname='Elevation') 

1683 depth = Distance.T(xmltagname='Depth') 

1684 azimuth = Azimuth.T(optional=True, xmltagname='Azimuth') 

1685 dip = Dip.T(optional=True, xmltagname='Dip') 

1686 type_list = List.T(Type.T(xmltagname='Type')) 

1687 sample_rate = SampleRate.T(optional=True, xmltagname='SampleRate') 

1688 sample_rate_ratio = SampleRateRatio.T(optional=True, 

1689 xmltagname='SampleRateRatio') 

1690 storage_format = String.T(optional=True, xmltagname='StorageFormat') 

1691 clock_drift = ClockDrift.T(optional=True, xmltagname='ClockDrift') 

1692 calibration_units = Units.T(optional=True, xmltagname='CalibrationUnits') 

1693 sensor = Equipment.T(optional=True, xmltagname='Sensor') 

1694 pre_amplifier = Equipment.T(optional=True, xmltagname='PreAmplifier') 

1695 data_logger = Equipment.T(optional=True, xmltagname='DataLogger') 

1696 equipment_list = List.T(Equipment.T(xmltagname='Equipment')) 

1697 response = Response.T(optional=True, xmltagname='Response') 

1698 

1699 @property 

1700 def position_values(self): 

1701 lat = self.latitude.value 

1702 lon = self.longitude.value 

1703 elevation = value_or_none(self.elevation) 

1704 depth = value_or_none(self.depth) 

1705 return lat, lon, elevation, depth 

1706 

1707 

1708class Station(BaseNode): 

1709 ''' 

1710 This type represents a Station epoch. It is common to only have a single 

1711 station epoch with the station's creation and termination dates as the 

1712 epoch start and end dates. 

1713 ''' 

1714 

1715 latitude = Latitude.T(xmltagname='Latitude') 

1716 longitude = Longitude.T(xmltagname='Longitude') 

1717 elevation = Distance.T(xmltagname='Elevation') 

1718 site = Site.T(default=Site.D(name=''), xmltagname='Site') 

1719 vault = Unicode.T(optional=True, xmltagname='Vault') 

1720 geology = Unicode.T(optional=True, xmltagname='Geology') 

1721 equipment_list = List.T(Equipment.T(xmltagname='Equipment')) 

1722 operator_list = List.T(Operator.T(xmltagname='Operator')) 

1723 creation_date = DummyAwareOptionalTimestamp.T( 

1724 optional=True, xmltagname='CreationDate') 

1725 termination_date = DummyAwareOptionalTimestamp.T( 

1726 optional=True, xmltagname='TerminationDate') 

1727 total_number_channels = Counter.T( 

1728 optional=True, xmltagname='TotalNumberChannels') 

1729 selected_number_channels = Counter.T( 

1730 optional=True, xmltagname='SelectedNumberChannels') 

1731 external_reference_list = List.T( 

1732 ExternalReference.T(xmltagname='ExternalReference')) 

1733 channel_list = List.T(Channel.T(xmltagname='Channel')) 

1734 

1735 @property 

1736 def position_values(self): 

1737 lat = self.latitude.value 

1738 lon = self.longitude.value 

1739 elevation = value_or_none(self.elevation) 

1740 return lat, lon, elevation 

1741 

1742 

1743class Network(BaseNode): 

1744 ''' 

1745 This type represents the Network layer, all station metadata is contained 

1746 within this element. The official name of the network or other descriptive 

1747 information can be included in the Description element. The Network can 

1748 contain 0 or more Stations. 

1749 ''' 

1750 

1751 total_number_stations = Counter.T(optional=True, 

1752 xmltagname='TotalNumberStations') 

1753 selected_number_stations = Counter.T(optional=True, 

1754 xmltagname='SelectedNumberStations') 

1755 station_list = List.T(Station.T(xmltagname='Station')) 

1756 

1757 @property 

1758 def station_code_list(self): 

1759 return sorted(set(s.code for s in self.station_list)) 

1760 

1761 @property 

1762 def sl_code_list(self): 

1763 sls = set() 

1764 for station in self.station_list: 

1765 for channel in station.channel_list: 

1766 sls.add((station.code, channel.location_code)) 

1767 

1768 return sorted(sls) 

1769 

1770 def summary(self, width=80, indent=4): 

1771 sls = self.sl_code_list or [(x,) for x in self.station_code_list] 

1772 lines = ['%s (%i):' % (self.code, len(sls))] 

1773 if sls: 

1774 ssls = ['.'.join(x for x in c if x) for c in sls] 

1775 w = max(len(x) for x in ssls) 

1776 n = (width - indent) / (w+1) 

1777 while ssls: 

1778 lines.append( 

1779 ' ' * indent + ' '.join(x.ljust(w) for x in ssls[:n])) 

1780 

1781 ssls[:n] = [] 

1782 

1783 return '\n'.join(lines) 

1784 

1785 

1786def value_or_none(x): 

1787 if x is not None: 

1788 return x.value 

1789 else: 

1790 return None 

1791 

1792 

1793def pyrocko_station_from_channels(nsl, channels, inconsistencies='warn'): 

1794 

1795 pos = lat, lon, elevation, depth = \ 

1796 channels[0].position_values 

1797 

1798 if not all(pos == x.position_values for x in channels): 

1799 info = '\n'.join( 

1800 ' %s: %s' % (x.code, x.position_values) for 

1801 x in channels) 

1802 

1803 mess = 'encountered inconsistencies in channel ' \ 

1804 'lat/lon/elevation/depth ' \ 

1805 'for %s.%s.%s: \n%s' % (nsl + (info,)) 

1806 

1807 if inconsistencies == 'raise': 

1808 raise InconsistentChannelLocations(mess) 

1809 

1810 elif inconsistencies == 'warn': 

1811 logger.warning(mess) 

1812 logger.warning(' -> using mean values') 

1813 

1814 apos = num.array([x.position_values for x in channels], dtype=float) 

1815 mlat, mlon, mele, mdep = num.nansum(apos, axis=0) \ 

1816 / num.sum(num.isfinite(apos), axis=0) 

1817 

1818 groups = {} 

1819 for channel in channels: 

1820 if channel.code not in groups: 

1821 groups[channel.code] = [] 

1822 

1823 groups[channel.code].append(channel) 

1824 

1825 pchannels = [] 

1826 for code in sorted(groups.keys()): 

1827 data = [ 

1828 (channel.code, value_or_none(channel.azimuth), 

1829 value_or_none(channel.dip)) 

1830 for channel in groups[code]] 

1831 

1832 azimuth, dip = util.consistency_merge( 

1833 data, 

1834 message='channel orientation values differ:', 

1835 error=inconsistencies) 

1836 

1837 pchannels.append( 

1838 pyrocko.model.Channel(code, azimuth=azimuth, dip=dip)) 

1839 

1840 return pyrocko.model.Station( 

1841 *nsl, 

1842 lat=mlat, 

1843 lon=mlon, 

1844 elevation=mele, 

1845 depth=mdep, 

1846 channels=pchannels) 

1847 

1848 

1849class FDSNStationXML(Object): 

1850 ''' 

1851 Top-level type for Station XML. Required field are Source (network ID of 

1852 the institution sending the message) and one or more Network containers or 

1853 one or more Station containers. 

1854 ''' 

1855 

1856 schema_version = Float.T(default=1.0, xmlstyle='attribute') 

1857 source = String.T(xmltagname='Source') 

1858 sender = String.T(optional=True, xmltagname='Sender') 

1859 module = String.T(optional=True, xmltagname='Module') 

1860 module_uri = String.T(optional=True, xmltagname='ModuleURI') 

1861 created = Timestamp.T(optional=True, xmltagname='Created') 

1862 network_list = List.T(Network.T(xmltagname='Network')) 

1863 

1864 xmltagname = 'FDSNStationXML' 

1865 guessable_xmlns = [guts_xmlns] 

1866 

1867 def __init__(self, *args, **kwargs): 

1868 Object.__init__(self, *args, **kwargs) 

1869 if self.created is None: 

1870 self.created = time.time() 

1871 

1872 def get_pyrocko_stations(self, nslcs=None, nsls=None, 

1873 time=None, timespan=None, 

1874 inconsistencies='warn'): 

1875 

1876 assert inconsistencies in ('raise', 'warn') 

1877 

1878 if nslcs is not None: 

1879 nslcs = set(nslcs) 

1880 

1881 if nsls is not None: 

1882 nsls = set(nsls) 

1883 

1884 tt = () 

1885 if time is not None: 

1886 tt = (time,) 

1887 elif timespan is not None: 

1888 tt = timespan 

1889 

1890 pstations = [] 

1891 for network in self.network_list: 

1892 if not network.spans(*tt): 

1893 continue 

1894 

1895 for station in network.station_list: 

1896 if not station.spans(*tt): 

1897 continue 

1898 

1899 if station.channel_list: 

1900 loc_to_channels = {} 

1901 for channel in station.channel_list: 

1902 if not channel.spans(*tt): 

1903 continue 

1904 

1905 loc = channel.location_code.strip() 

1906 if loc not in loc_to_channels: 

1907 loc_to_channels[loc] = [] 

1908 

1909 loc_to_channels[loc].append(channel) 

1910 

1911 for loc in sorted(loc_to_channels.keys()): 

1912 channels = loc_to_channels[loc] 

1913 if nslcs is not None: 

1914 channels = [channel for channel in channels 

1915 if (network.code, station.code, loc, 

1916 channel.code) in nslcs] 

1917 

1918 if not channels: 

1919 continue 

1920 

1921 nsl = network.code, station.code, loc 

1922 if nsls is not None and nsl not in nsls: 

1923 continue 

1924 

1925 pstations.append( 

1926 pyrocko_station_from_channels( 

1927 nsl, 

1928 channels, 

1929 inconsistencies=inconsistencies)) 

1930 else: 

1931 pstations.append(pyrocko.model.Station( 

1932 network.code, station.code, '*', 

1933 lat=station.latitude.value, 

1934 lon=station.longitude.value, 

1935 elevation=value_or_none(station.elevation), 

1936 name=station.description or '')) 

1937 

1938 return pstations 

1939 

1940 @classmethod 

1941 def from_pyrocko_stations( 

1942 cls, pyrocko_stations, add_flat_responses_from=None): 

1943 

1944 ''' 

1945 Generate :py:class:`FDSNStationXML` from list of 

1946 :py:class;`pyrocko.model.Station` instances. 

1947 

1948 :param pyrocko_stations: list of :py:class;`pyrocko.model.Station` 

1949 instances. 

1950 :param add_flat_responses_from: unit, 'M', 'M/S' or 'M/S**2' 

1951 ''' 

1952 from collections import defaultdict 

1953 network_dict = defaultdict(list) 

1954 

1955 if add_flat_responses_from: 

1956 assert add_flat_responses_from in ('M', 'M/S', 'M/S**2') 

1957 extra = dict( 

1958 response=Response( 

1959 instrument_sensitivity=Sensitivity( 

1960 value=1.0, 

1961 frequency=1.0, 

1962 input_units=Units(name=add_flat_responses_from)))) 

1963 else: 

1964 extra = {} 

1965 

1966 have_offsets = set() 

1967 for s in pyrocko_stations: 

1968 

1969 if s.north_shift != 0.0 or s.east_shift != 0.0: 

1970 have_offsets.add(s.nsl()) 

1971 

1972 network, station, location = s.nsl() 

1973 channel_list = [] 

1974 for c in s.channels: 

1975 channel_list.append( 

1976 Channel( 

1977 location_code=location, 

1978 code=c.name, 

1979 latitude=Latitude(value=s.effective_lat), 

1980 longitude=Longitude(value=s.effective_lon), 

1981 elevation=Distance(value=s.elevation), 

1982 depth=Distance(value=s.depth), 

1983 azimuth=Azimuth(value=c.azimuth), 

1984 dip=Dip(value=c.dip), 

1985 **extra 

1986 ) 

1987 ) 

1988 

1989 network_dict[network].append( 

1990 Station( 

1991 code=station, 

1992 latitude=Latitude(value=s.effective_lat), 

1993 longitude=Longitude(value=s.effective_lon), 

1994 elevation=Distance(value=s.elevation), 

1995 channel_list=channel_list) 

1996 ) 

1997 

1998 if have_offsets: 

1999 logger.warning( 

2000 'StationXML does not support Cartesian offsets in ' 

2001 'coordinates. Storing effective lat/lon for stations: %s' % 

2002 ', '.join('.'.join(nsl) for nsl in sorted(have_offsets))) 

2003 

2004 timestamp = util.to_time_float(time.time()) 

2005 network_list = [] 

2006 for k, station_list in network_dict.items(): 

2007 

2008 network_list.append( 

2009 Network( 

2010 code=k, station_list=station_list, 

2011 total_number_stations=len(station_list))) 

2012 

2013 sxml = FDSNStationXML( 

2014 source='from pyrocko stations list', created=timestamp, 

2015 network_list=network_list) 

2016 

2017 sxml.validate() 

2018 return sxml 

2019 

2020 def iter_network_stations( 

2021 self, net=None, sta=None, time=None, timespan=None): 

2022 

2023 tt = () 

2024 if time is not None: 

2025 tt = (time,) 

2026 elif timespan is not None: 

2027 tt = timespan 

2028 

2029 for network in self.network_list: 

2030 if not network.spans(*tt) or ( 

2031 net is not None and network.code != net): 

2032 continue 

2033 

2034 for station in network.station_list: 

2035 if not station.spans(*tt) or ( 

2036 sta is not None and station.code != sta): 

2037 continue 

2038 

2039 yield (network, station) 

2040 

2041 def iter_network_station_channels( 

2042 self, net=None, sta=None, loc=None, cha=None, 

2043 time=None, timespan=None): 

2044 

2045 if loc is not None: 

2046 loc = loc.strip() 

2047 

2048 tt = () 

2049 if time is not None: 

2050 tt = (time,) 

2051 elif timespan is not None: 

2052 tt = timespan 

2053 

2054 for network in self.network_list: 

2055 if not network.spans(*tt) or ( 

2056 net is not None and network.code != net): 

2057 continue 

2058 

2059 for station in network.station_list: 

2060 if not station.spans(*tt) or ( 

2061 sta is not None and station.code != sta): 

2062 continue 

2063 

2064 if station.channel_list: 

2065 for channel in station.channel_list: 

2066 if (not channel.spans(*tt) or 

2067 (cha is not None and channel.code != cha) or 

2068 (loc is not None and 

2069 channel.location_code.strip() != loc)): 

2070 continue 

2071 

2072 yield (network, station, channel) 

2073 

2074 def get_channel_groups(self, net=None, sta=None, loc=None, cha=None, 

2075 time=None, timespan=None): 

2076 

2077 groups = {} 

2078 for network, station, channel in self.iter_network_station_channels( 

2079 net, sta, loc, cha, time=time, timespan=timespan): 

2080 

2081 net = network.code 

2082 sta = station.code 

2083 cha = channel.code 

2084 loc = channel.location_code.strip() 

2085 if len(cha) == 3: 

2086 bic = cha[:2] # band and intrument code according to SEED 

2087 elif len(cha) == 1: 

2088 bic = '' 

2089 else: 

2090 bic = cha 

2091 

2092 if channel.response and \ 

2093 channel.response.instrument_sensitivity and \ 

2094 channel.response.instrument_sensitivity.input_units: 

2095 

2096 unit = channel.response.instrument_sensitivity\ 

2097 .input_units.name.upper() 

2098 else: 

2099 unit = None 

2100 

2101 bic = (bic, unit) 

2102 

2103 k = net, sta, loc 

2104 if k not in groups: 

2105 groups[k] = {} 

2106 

2107 if bic not in groups[k]: 

2108 groups[k][bic] = [] 

2109 

2110 groups[k][bic].append(channel) 

2111 

2112 for nsl, bic_to_channels in groups.items(): 

2113 bad_bics = [] 

2114 for bic, channels in bic_to_channels.items(): 

2115 sample_rates = [] 

2116 for channel in channels: 

2117 sample_rates.append(channel.sample_rate.value) 

2118 

2119 if not same(sample_rates): 

2120 scs = ','.join(channel.code for channel in channels) 

2121 srs = ', '.join('%e' % x for x in sample_rates) 

2122 err = 'ignoring channels with inconsistent sampling ' + \ 

2123 'rates (%s.%s.%s.%s: %s)' % (nsl + (scs, srs)) 

2124 

2125 logger.warning(err) 

2126 bad_bics.append(bic) 

2127 

2128 for bic in bad_bics: 

2129 del bic_to_channels[bic] 

2130 

2131 return groups 

2132 

2133 def choose_channels( 

2134 self, 

2135 target_sample_rate=None, 

2136 priority_band_code=['H', 'B', 'M', 'L', 'V', 'E', 'S'], 

2137 priority_units=['M/S', 'M/S**2'], 

2138 priority_instrument_code=['H', 'L'], 

2139 time=None, 

2140 timespan=None): 

2141 

2142 nslcs = {} 

2143 for nsl, bic_to_channels in self.get_channel_groups( 

2144 time=time, timespan=timespan).items(): 

2145 

2146 useful_bics = [] 

2147 for bic, channels in bic_to_channels.items(): 

2148 rate = channels[0].sample_rate.value 

2149 

2150 if target_sample_rate is not None and \ 

2151 rate < target_sample_rate*0.99999: 

2152 continue 

2153 

2154 if len(bic[0]) == 2: 

2155 if bic[0][0] not in priority_band_code: 

2156 continue 

2157 

2158 if bic[0][1] not in priority_instrument_code: 

2159 continue 

2160 

2161 unit = bic[1] 

2162 

2163 prio_unit = len(priority_units) 

2164 try: 

2165 prio_unit = priority_units.index(unit) 

2166 except ValueError: 

2167 pass 

2168 

2169 prio_inst = len(priority_instrument_code) 

2170 prio_band = len(priority_band_code) 

2171 if len(channels[0].code) == 3: 

2172 try: 

2173 prio_inst = priority_instrument_code.index( 

2174 channels[0].code[1]) 

2175 except ValueError: 

2176 pass 

2177 

2178 try: 

2179 prio_band = priority_band_code.index( 

2180 channels[0].code[0]) 

2181 except ValueError: 

2182 pass 

2183 

2184 if target_sample_rate is None: 

2185 rate = -rate 

2186 

2187 useful_bics.append((-len(channels), prio_band, rate, prio_unit, 

2188 prio_inst, bic)) 

2189 

2190 useful_bics.sort() 

2191 

2192 for _, _, rate, _, _, bic in useful_bics: 

2193 channels = sorted( 

2194 bic_to_channels[bic], 

2195 key=lambda channel: channel.code) 

2196 

2197 if channels: 

2198 for channel in channels: 

2199 nslcs[nsl + (channel.code,)] = channel 

2200 

2201 break 

2202 

2203 return nslcs 

2204 

2205 def get_pyrocko_response( 

2206 self, nslc, 

2207 time=None, timespan=None, fake_input_units=None, stages=(0, 1)): 

2208 

2209 net, sta, loc, cha = nslc 

2210 resps = [] 

2211 for _, _, channel in self.iter_network_station_channels( 

2212 net, sta, loc, cha, time=time, timespan=timespan): 

2213 resp = channel.response 

2214 if resp: 

2215 resp.check_sample_rates(channel) 

2216 resp.check_units() 

2217 resps.append(resp.get_pyrocko_response( 

2218 '.'.join(nslc), 

2219 fake_input_units=fake_input_units, 

2220 stages=stages).expect_one()) 

2221 

2222 if not resps: 

2223 raise NoResponseInformation('%s.%s.%s.%s' % nslc) 

2224 elif len(resps) > 1: 

2225 raise MultipleResponseInformation('%s.%s.%s.%s' % nslc) 

2226 

2227 return resps[0] 

2228 

2229 @property 

2230 def n_code_list(self): 

2231 return sorted(set(x.code for x in self.network_list)) 

2232 

2233 @property 

2234 def ns_code_list(self): 

2235 nss = set() 

2236 for network in self.network_list: 

2237 for station in network.station_list: 

2238 nss.add((network.code, station.code)) 

2239 

2240 return sorted(nss) 

2241 

2242 @property 

2243 def nsl_code_list(self): 

2244 nsls = set() 

2245 for network in self.network_list: 

2246 for station in network.station_list: 

2247 for channel in station.channel_list: 

2248 nsls.add( 

2249 (network.code, station.code, channel.location_code)) 

2250 

2251 return sorted(nsls) 

2252 

2253 @property 

2254 def nslc_code_list(self): 

2255 nslcs = set() 

2256 for network in self.network_list: 

2257 for station in network.station_list: 

2258 for channel in station.channel_list: 

2259 nslcs.add( 

2260 (network.code, station.code, channel.location_code, 

2261 channel.code)) 

2262 

2263 return sorted(nslcs) 

2264 

2265 def summary(self): 

2266 lst = [ 

2267 'number of n codes: %i' % len(self.n_code_list), 

2268 'number of ns codes: %i' % len(self.ns_code_list), 

2269 'number of nsl codes: %i' % len(self.nsl_code_list), 

2270 'number of nslc codes: %i' % len(self.nslc_code_list) 

2271 ] 

2272 return '\n'.join(lst) 

2273 

2274 def summary_stages(self): 

2275 data = [] 

2276 for network, station, channel in self.iter_network_station_channels(): 

2277 nslc = (network.code, station.code, channel.location_code, 

2278 channel.code) 

2279 

2280 stages = [] 

2281 in_units = '?' 

2282 out_units = '?' 

2283 if channel.response: 

2284 sens = channel.response.instrument_sensitivity 

2285 if sens: 

2286 in_units = sens.input_units.name.upper() 

2287 out_units = sens.output_units.name.upper() 

2288 

2289 for stage in channel.response.stage_list: 

2290 stages.append(stage.summary()) 

2291 

2292 data.append( 

2293 (nslc, tts(channel.start_date), tts(channel.end_date), 

2294 in_units, out_units, stages)) 

2295 

2296 data.sort() 

2297 

2298 lst = [] 

2299 for nslc, stmin, stmax, in_units, out_units, stages in data: 

2300 lst.append(' %s: %s - %s, %s -> %s' % ( 

2301 '.'.join(nslc), stmin, stmax, in_units, out_units)) 

2302 for stage in stages: 

2303 lst.append(' %s' % stage) 

2304 

2305 return '\n'.join(lst) 

2306 

2307 def _check_overlaps(self): 

2308 by_nslc = {} 

2309 for network in self.network_list: 

2310 for station in network.station_list: 

2311 for channel in station.channel_list: 

2312 nslc = (network.code, station.code, channel.location_code, 

2313 channel.code) 

2314 if nslc not in by_nslc: 

2315 by_nslc[nslc] = [] 

2316 

2317 by_nslc[nslc].append(channel) 

2318 

2319 errors = [] 

2320 for nslc, channels in by_nslc.items(): 

2321 errors.extend(check_overlaps('Channel', nslc, channels)) 

2322 

2323 return errors 

2324 

2325 def check(self): 

2326 errors = [] 

2327 for meth in [self._check_overlaps]: 

2328 errors.extend(meth()) 

2329 

2330 if errors: 

2331 raise Inconsistencies( 

2332 'Inconsistencies found in StationXML:\n ' 

2333 + '\n '.join(errors)) 

2334 

2335 

2336def load_channel_table(stream): 

2337 

2338 networks = {} 

2339 stations = {} 

2340 

2341 for line in stream: 

2342 line = str(line.decode('ascii')) 

2343 if line.startswith('#'): 

2344 continue 

2345 

2346 t = line.rstrip().split('|') 

2347 

2348 if len(t) != 17: 

2349 logger.warning('Invalid channel record: %s' % line) 

2350 continue 

2351 

2352 (net, sta, loc, cha, lat, lon, ele, dep, azi, dip, sens, scale, 

2353 scale_freq, scale_units, sample_rate, start_date, end_date) = t 

2354 

2355 try: 

2356 scale = float(scale) 

2357 except ValueError: 

2358 scale = None 

2359 

2360 try: 

2361 scale_freq = float(scale_freq) 

2362 except ValueError: 

2363 scale_freq = None 

2364 

2365 try: 

2366 depth = float(dep) 

2367 except ValueError: 

2368 depth = 0.0 

2369 

2370 try: 

2371 azi = float(azi) 

2372 dip = float(dip) 

2373 except ValueError: 

2374 azi = None 

2375 dip = None 

2376 

2377 try: 

2378 if net not in networks: 

2379 network = Network(code=net) 

2380 else: 

2381 network = networks[net] 

2382 

2383 if (net, sta) not in stations: 

2384 station = Station( 

2385 code=sta, latitude=lat, longitude=lon, elevation=ele) 

2386 

2387 station.regularize() 

2388 else: 

2389 station = stations[net, sta] 

2390 

2391 if scale: 

2392 resp = Response( 

2393 instrument_sensitivity=Sensitivity( 

2394 value=scale, 

2395 frequency=scale_freq, 

2396 input_units=scale_units)) 

2397 else: 

2398 resp = None 

2399 

2400 channel = Channel( 

2401 code=cha, 

2402 location_code=loc.strip(), 

2403 latitude=lat, 

2404 longitude=lon, 

2405 elevation=ele, 

2406 depth=depth, 

2407 azimuth=azi, 

2408 dip=dip, 

2409 sensor=Equipment(description=sens), 

2410 response=resp, 

2411 sample_rate=sample_rate, 

2412 start_date=start_date, 

2413 end_date=end_date or None) 

2414 

2415 channel.regularize() 

2416 

2417 except ValidationError: 

2418 raise InvalidRecord(line) 

2419 

2420 if net not in networks: 

2421 networks[net] = network 

2422 

2423 if (net, sta) not in stations: 

2424 stations[net, sta] = station 

2425 network.station_list.append(station) 

2426 

2427 station.channel_list.append(channel) 

2428 

2429 return FDSNStationXML( 

2430 source='created from table input', 

2431 created=time.time(), 

2432 network_list=sorted(networks.values(), key=lambda x: x.code)) 

2433 

2434 

2435def primitive_merge(sxs): 

2436 networks = [] 

2437 for sx in sxs: 

2438 networks.extend(sx.network_list) 

2439 

2440 return FDSNStationXML( 

2441 source='merged from different sources', 

2442 created=time.time(), 

2443 network_list=copy.deepcopy( 

2444 sorted(networks, key=lambda x: x.code))) 

2445 

2446 

2447def split_channels(sx): 

2448 for nslc in sx.nslc_code_list: 

2449 network_list = sx.network_list 

2450 network_list_filtered = [ 

2451 network for network in network_list 

2452 if network.code == nslc[0]] 

2453 

2454 for network in network_list_filtered: 

2455 sx.network_list = [network] 

2456 station_list = network.station_list 

2457 station_list_filtered = [ 

2458 station for station in station_list 

2459 if station.code == nslc[1]] 

2460 

2461 for station in station_list_filtered: 

2462 network.station_list = [station] 

2463 channel_list = station.channel_list 

2464 station.channel_list = [ 

2465 channel for channel in channel_list 

2466 if (channel.location_code, channel.code) == nslc[2:4]] 

2467 

2468 yield nslc, copy.deepcopy(sx) 

2469 

2470 station.channel_list = channel_list 

2471 

2472 network.station_list = station_list 

2473 

2474 sx.network_list = network_list 

2475 

2476 

2477if __name__ == '__main__': 

2478 from optparse import OptionParser 

2479 

2480 util.setup_logging('pyrocko.io.stationxml', 'warning') 

2481 

2482 usage = \ 

2483 'python -m pyrocko.io.stationxml check|stats|stages ' \ 

2484 '<filename> [options]' 

2485 

2486 description = '''Torture StationXML file.''' 

2487 

2488 parser = OptionParser( 

2489 usage=usage, 

2490 description=description, 

2491 formatter=util.BetterHelpFormatter()) 

2492 

2493 (options, args) = parser.parse_args(sys.argv[1:]) 

2494 

2495 if len(args) != 2: 

2496 parser.print_help() 

2497 sys.exit(1) 

2498 

2499 action, path = args 

2500 

2501 sx = load_xml(filename=path) 

2502 if action == 'check': 

2503 try: 

2504 sx.check() 

2505 except Inconsistencies as e: 

2506 logger.error(e) 

2507 sys.exit(1) 

2508 

2509 elif action == 'stats': 

2510 print(sx.summary()) 

2511 

2512 elif action == 'stages': 

2513 print(sx.summary_stages()) 

2514 

2515 else: 

2516 parser.print_help() 

2517 sys.exit('unknown action: %s' % action)