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

1177 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-23 12:34 +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 

18from collections import defaultdict 

19 

20import numpy as num 

21 

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

23 Unicode, Int, Float, List, Object, Timestamp, 

24 ValidationError, TBase, re_tz, Any, Tuple) 

25from pyrocko.guts import load_xml # noqa 

26from pyrocko.util import hpfloat, time_to_str, get_time_float 

27 

28import pyrocko.model 

29from pyrocko import util, response 

30 

31guts_prefix = 'sx' 

32 

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

34 

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

36 

37conversion = { 

38 ('M', 'M'): None, 

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

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

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

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

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

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

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

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

47 

48 

49unit_to_quantity = { 

50 'M/S': 'velocity', 

51 'M': 'displacement', 

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

53 'V': 'voltage', 

54 'COUNTS': 'counts', 

55 'COUNT': 'counts', 

56 'PA': 'pressure', 

57 'RAD': 'rotation-displacement', 

58 'R': 'rotation-displacement', 

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

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

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

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

63 

64 

65def to_quantity(unit, context, delivery): 

66 

67 if unit is None: 

68 return None 

69 

70 name = unit.name.upper() 

71 if name in unit_to_quantity: 

72 return unit_to_quantity[name] 

73 else: 

74 delivery.log.append(( 

75 'warning', 

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

77 unit.name.upper() + ( 

78 ' (%s)' % unit.description 

79 if unit.description else '')), 

80 context)) 

81 

82 return 'unsupported_quantity(%s)' % unit 

83 

84 

85class StationXMLError(Exception): 

86 pass 

87 

88 

89class Inconsistencies(StationXMLError): 

90 pass 

91 

92 

93class NoResponseInformation(StationXMLError): 

94 pass 

95 

96 

97class MultipleResponseInformation(StationXMLError): 

98 pass 

99 

100 

101class InconsistentResponseInformation(StationXMLError): 

102 pass 

103 

104 

105class InconsistentChannelLocations(StationXMLError): 

106 pass 

107 

108 

109class InvalidRecord(StationXMLError): 

110 def __init__(self, line): 

111 StationXMLError.__init__(self) 

112 self._line = line 

113 

114 def __str__(self): 

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

116 

117 

118_exceptions = dict( 

119 Inconsistencies=Inconsistencies, 

120 NoResponseInformation=NoResponseInformation, 

121 MultipleResponseInformation=MultipleResponseInformation, 

122 InconsistentResponseInformation=InconsistentResponseInformation, 

123 InconsistentChannelLocations=InconsistentChannelLocations, 

124 InvalidRecord=InvalidRecord, 

125 ValueError=ValueError) 

126 

127 

128_logs = dict( 

129 info=logger.info, 

130 warning=logger.warning, 

131 error=logger.error) 

132 

133 

134class DeliveryError(StationXMLError): 

135 pass 

136 

137 

138class Delivery(Object): 

139 

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

141 if payload is None: 

142 payload = [] 

143 

144 if log is None: 

145 log = [] 

146 

147 if errors is None: 

148 errors = [] 

149 

150 if error is not None: 

151 errors.append(error) 

152 

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

154 

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

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

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

158 

159 def extend(self, other): 

160 self.payload.extend(other.payload) 

161 self.log.extend(other.log) 

162 self.errors.extend(other.errors) 

163 

164 def extend_without_payload(self, other): 

165 self.log.extend(other.log) 

166 self.errors.extend(other.errors) 

167 return other.payload 

168 

169 def emit_log(self): 

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

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

172 _logs[name](message) 

173 

174 def expect(self, quiet=False): 

175 if not quiet: 

176 self.emit_log() 

177 

178 if self.errors: 

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

180 if context: 

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

182 

183 if len(self.errors) > 1: 

184 message += ' Additional errors pending.' 

185 

186 raise _exceptions[name](message) 

187 

188 return self.payload 

189 

190 def expect_one(self, quiet=False): 

191 payload = self.expect(quiet=quiet) 

192 if len(payload) != 1: 

193 raise DeliveryError( 

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

195 

196 return payload[0] 

197 

198 

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

200 words = s.split() 

201 lines = [] 

202 t = [] 

203 n = 0 

204 for w in words: 

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

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

207 n = indent 

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

209 

210 t.append(w) 

211 n += len(w) + 1 

212 

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

214 return '\n'.join(lines) 

215 

216 

217def same(x, eps=0.0): 

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

219 return False 

220 

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

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

223 else: 

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

225 

226 

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

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

229 

230 

231def evaluate1(resp, f): 

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

233 

234 

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

236 

237 try: 

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

239 except response.InvalidResponseError as e: 

240 return Delivery(log=[( 

241 'warning', 

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

243 context)]) 

244 

245 if value_resp == 0.0: 

246 return Delivery(log=[( 

247 'warning', 

248 '%s\n' 

249 ' computed response is zero' % prelude, 

250 context)]) 

251 

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

253 

254 if num.abs(diff_db) > limit_db: 

255 return Delivery(log=[( 

256 'warning', 

257 '%s\n' 

258 ' reported value: %g\n' 

259 ' computed value: %g\n' 

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

261 ' factor: %g\n' 

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

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

264 prelude, 

265 value, 

266 value_resp, 

267 frequency, 

268 value_resp/value, 

269 diff_db, 

270 limit_db), 

271 context)]) 

272 

273 return Delivery() 

274 

275 

276def tts(t): 

277 if t is None: 

278 return '<none>' 

279 else: 

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

281 

282 

283def le_open_left(a, b): 

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

285 

286 

287def le_open_right(a, b): 

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

289 

290 

291def eq_open(a, b): 

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

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

294 

295 

296def contains(a, b): 

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

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

299 

300 

301def find_containing(candidates, node): 

302 for candidate in candidates: 

303 if contains(candidate, node): 

304 return candidate 

305 

306 return None 

307 

308 

309this_year = time.gmtime()[0] 

310 

311 

312class DummyAwareOptionalTimestamp(Object): 

313 ''' 

314 Optional timestamp with support for some common placeholder values. 

315 

316 Some StationXML files contain arbitrary placeholder values for open end 

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

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

319 this type. 

320 ''' 

321 dummy_for = (hpfloat, float) 

322 dummy_for_description = 'pyrocko.util.get_time_float' 

323 

324 class __T(TBase): 

325 

326 def regularize_extra(self, val): 

327 time_float = get_time_float() 

328 

329 if isinstance(val, datetime.datetime): 

330 tt = val.utctimetuple() 

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

332 

333 elif isinstance(val, datetime.date): 

334 tt = val.timetuple() 

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

336 

337 elif isinstance(val, str): 

338 val = val.strip() 

339 

340 tz_offset = 0 

341 

342 m = re_tz.search(val) 

343 if m: 

344 sh = m.group(2) 

345 sm = m.group(4) 

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

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

348 

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

350 

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

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

353 

354 try: 

355 val = util.str_to_time(val) - tz_offset 

356 

357 except util.TimeStrError: 

358 year = int(val[:4]) 

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

360 if year > this_year + 100: 

361 return None # StationXML contained a dummy date 

362 

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

364 return None 

365 

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

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

368 return None # StationXML contained a dummy date 

369 

370 raise 

371 

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

373 val = time_float(val) 

374 

375 else: 

376 raise ValidationError( 

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

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

379 

380 return val 

381 

382 def to_save(self, val): 

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

384 .rstrip('0').rstrip('.') 

385 

386 def to_save_xml(self, val): 

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

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

389 

390 

391class Nominal(StringChoice): 

392 choices = [ 

393 'NOMINAL', 

394 'CALCULATED'] 

395 

396 

397class Email(UnicodePattern): 

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

399 

400 

401class RestrictedStatus(StringChoice): 

402 choices = [ 

403 'open', 

404 'closed', 

405 'partial'] 

406 

407 

408class Type(StringChoice): 

409 choices = [ 

410 'TRIGGERED', 

411 'CONTINUOUS', 

412 'HEALTH', 

413 'GEOPHYSICAL', 

414 'WEATHER', 

415 'FLAG', 

416 'SYNTHESIZED', 

417 'INPUT', 

418 'EXPERIMENTAL', 

419 'MAINTENANCE', 

420 'BEAM'] 

421 

422 class __T(StringChoice.T): 

423 def validate_extra(self, val): 

424 if val not in self.choices: 

425 logger.warning( 

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

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

428 

429 

430class PzTransferFunction(StringChoice): 

431 choices = [ 

432 'LAPLACE (RADIANS/SECOND)', 

433 'LAPLACE (HERTZ)', 

434 'DIGITAL (Z-TRANSFORM)'] 

435 

436 

437class Symmetry(StringChoice): 

438 choices = [ 

439 'NONE', 

440 'EVEN', 

441 'ODD'] 

442 

443 

444class CfTransferFunction(StringChoice): 

445 

446 class __T(StringChoice.T): 

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

448 if regularize: 

449 try: 

450 val = str(val) 

451 except ValueError: 

452 raise ValidationError( 

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

454 repr(val))) 

455 

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

457 

458 self.validate_extra(val) 

459 return val 

460 

461 choices = [ 

462 'ANALOG (RADIANS/SECOND)', 

463 'ANALOG (HERTZ)', 

464 'DIGITAL'] 

465 

466 replacements = { 

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

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

469 } 

470 

471 

472class Approximation(StringChoice): 

473 choices = [ 

474 'MACLAURIN'] 

475 

476 

477class PhoneNumber(StringPattern): 

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

479 

480 

481class Site(Object): 

482 ''' 

483 Description of a site location using name and optional geopolitical 

484 boundaries (country, city, etc.). 

485 ''' 

486 

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

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

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

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

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

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

493 

494 

495class ExternalReference(Object): 

496 ''' 

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

498 want to reference in StationXML. 

499 ''' 

500 

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

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

503 

504 

505class Units(Object): 

506 ''' 

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

508 ''' 

509 

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

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

512 

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

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

515 

516 

517class Counter(Int): 

518 pass 

519 

520 

521class SampleRateRatio(Object): 

522 ''' 

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

524 ''' 

525 

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

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

528 

529 

530class Gain(Object): 

531 ''' 

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

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

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

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

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

537 ''' 

538 

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

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

541 

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

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

544 

545 def summary(self): 

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

547 

548 

549class NumeratorCoefficient(Object): 

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

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

552 

553 

554class FloatNoUnit(Object): 

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

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

557 

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

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

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

561 

562 

563class FloatWithUnit(FloatNoUnit): 

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

565 

566 

567class Equipment(Object): 

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

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

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

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

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

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

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

575 installation_date = DummyAwareOptionalTimestamp.T( 

576 optional=True, 

577 xmltagname='InstallationDate') 

578 removal_date = DummyAwareOptionalTimestamp.T( 

579 optional=True, 

580 xmltagname='RemovalDate') 

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

582 

583 

584class PhoneNumber(Object): 

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

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

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

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

589 

590 

591class BaseFilter(Object): 

592 ''' 

593 The BaseFilter is derived by all filters. 

594 ''' 

595 

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

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

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

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

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

601 

602 

603class Sensitivity(Gain): 

604 ''' 

605 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional 

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

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

608 decibels specified in FrequencyDBVariation. 

609 ''' 

610 

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

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

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

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

615 frequency_db_variation = Float.T(optional=True, 

616 xmltagname='FrequencyDBVariation') 

617 

618 def get_pyrocko_response(self): 

619 return Delivery( 

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

621 

622 

623class Coefficient(FloatNoUnit): 

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

625 

626 

627class PoleZero(Object): 

628 ''' 

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

630 ''' 

631 

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

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

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

635 

636 def value(self): 

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

638 

639 

640class ClockDrift(FloatWithUnit): 

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

642 xmlstyle='attribute') # fixed 

643 

644 

645class Second(FloatWithUnit): 

646 ''' 

647 A time value in seconds. 

648 ''' 

649 

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

651 # fixed unit 

652 

653 

654class Voltage(FloatWithUnit): 

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

656 # fixed unit 

657 

658 

659class Angle(FloatWithUnit): 

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

661 # fixed unit 

662 

663 

664class Azimuth(FloatWithUnit): 

665 ''' 

666 Instrument azimuth, degrees clockwise from North. 

667 ''' 

668 

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

670 # fixed unit 

671 

672 

673class Dip(FloatWithUnit): 

674 ''' 

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

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

677 ''' 

678 

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

680 # fixed unit 

681 

682 

683class Distance(FloatWithUnit): 

684 ''' 

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

686 ''' 

687 

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

689 # NOT fixed unit! 

690 

691 

692class Frequency(FloatWithUnit): 

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

694 # fixed unit 

695 

696 

697class SampleRate(FloatWithUnit): 

698 ''' 

699 Sample rate in samples per second. 

700 ''' 

701 

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

703 # fixed unit 

704 

705 

706class Person(Object): 

707 ''' 

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

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

710 ''' 

711 

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

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

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

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

716 

717 

718class FIR(BaseFilter): 

719 ''' 

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

721 also commonly documented using the Coefficients element. 

722 ''' 

723 

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

725 numerator_coefficient_list = List.T( 

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

727 

728 def summary(self): 

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

730 self.get_ncoefficiencs(), 

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

732 

733 def get_effective_coefficients(self): 

734 b = num.array( 

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

736 dtype=float) 

737 

738 if self.symmetry == 'ODD': 

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

740 elif self.symmetry == 'EVEN': 

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

742 

743 return b 

744 

745 def get_effective_symmetry(self): 

746 if self.symmetry == 'NONE': 

747 b = self.get_effective_coefficients() 

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

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

750 

751 return self.symmetry 

752 

753 def get_ncoefficiencs(self): 

754 nf = len(self.numerator_coefficient_list) 

755 if self.symmetry == 'ODD': 

756 nc = nf * 2 + 1 

757 elif self.symmetry == 'EVEN': 

758 nc = nf * 2 

759 else: 

760 nc = nf 

761 

762 return nc 

763 

764 def estimate_delay(self, deltat): 

765 nc = self.get_ncoefficiencs() 

766 if nc > 0: 

767 return deltat * (nc - 1) / 2.0 

768 else: 

769 return 0.0 

770 

771 def get_pyrocko_response( 

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

773 

774 context += self.summary() 

775 

776 if not self.numerator_coefficient_list: 

777 return Delivery([]) 

778 

779 b = self.get_effective_coefficients() 

780 

781 log = [] 

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

783 

784 if not deltat: 

785 log.append(( 

786 'error', 

787 'Digital response requires knowledge about sampling ' 

788 'interval. Response will be unusable.', 

789 context)) 

790 

791 resp = response.DigitalFilterResponse( 

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

793 

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

795 normalization_frequency = 0.0 

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

797 

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

799 log.append(( 

800 'warning', 

801 'FIR filter coefficients are not normalized. Normalizing ' 

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

803 

804 resp = response.DigitalFilterResponse( 

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

806 drop_phase=drop_phase) 

807 

808 resps = [resp] 

809 

810 if not drop_phase: 

811 resps.extend(delay_responses) 

812 

813 return Delivery(resps, log=log) 

814 

815 

816class Coefficients(BaseFilter): 

817 ''' 

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

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

820 instead. Corresponds to SEED blockette 54. 

821 ''' 

822 

823 cf_transfer_function_type = CfTransferFunction.T( 

824 xmltagname='CfTransferFunctionType') 

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

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

827 

828 def summary(self): 

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

830 'ABC?'[ 

831 CfTransferFunction.choices.index( 

832 self.cf_transfer_function_type)], 

833 len(self.numerator_list), 

834 len(self.denominator_list), 

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

836 

837 def estimate_delay(self, deltat): 

838 nc = len(self.numerator_list) 

839 if nc > 0: 

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

841 else: 

842 return 0.0 

843 

844 def is_symmetric_fir(self): 

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

846 return False 

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

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

849 

850 def get_pyrocko_response( 

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

852 

853 context += self.summary() 

854 

855 factor = 1.0 

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

857 factor = 2.0*math.pi 

858 

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

860 return Delivery(payload=[]) 

861 

862 b = num.array( 

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

864 

865 a = num.array( 

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

867 dtype=float) 

868 

869 log = [] 

870 resps = [] 

871 if self.cf_transfer_function_type in [ 

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

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

874 

875 elif self.cf_transfer_function_type == 'DIGITAL': 

876 if not deltat: 

877 log.append(( 

878 'error', 

879 'Digital response requires knowledge about sampling ' 

880 'interval. Response will be unusable.', 

881 context)) 

882 

883 drop_phase = self.is_symmetric_fir() 

884 resp = response.DigitalFilterResponse( 

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

886 

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

888 normalization = num.abs(evaluate1( 

889 resp, normalization_frequency)) 

890 

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

892 log.append(( 

893 'warning', 

894 'FIR filter coefficients are not normalized. ' 

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

896 context)) 

897 

898 resp = response.DigitalFilterResponse( 

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

900 drop_phase=drop_phase) 

901 

902 resps.append(resp) 

903 

904 if not drop_phase: 

905 resps.extend(delay_responses) 

906 

907 else: 

908 return Delivery(error=( 

909 'ValueError', 

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

911 self.cf_transfer_function_type))) 

912 

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

914 

915 

916class Latitude(FloatWithUnit): 

917 ''' 

918 Type for latitude coordinate. 

919 ''' 

920 

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

922 # fixed unit 

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

924 

925 

926class Longitude(FloatWithUnit): 

927 ''' 

928 Type for longitude coordinate. 

929 ''' 

930 

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

932 # fixed unit 

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

934 

935 

936class PolesZeros(BaseFilter): 

937 ''' 

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

939 ''' 

940 

941 pz_transfer_function_type = PzTransferFunction.T( 

942 xmltagname='PzTransferFunctionType') 

943 normalization_factor = Float.T( 

944 default=1.0, 

945 xmltagname='NormalizationFactor') 

946 normalization_frequency = Frequency.T( 

947 optional=True, # but required by standard 

948 xmltagname='NormalizationFrequency') 

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

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

951 

952 def summary(self): 

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

954 'ABC?'[ 

955 PzTransferFunction.choices.index( 

956 self.pz_transfer_function_type)], 

957 len(self.pole_list), 

958 len(self.zero_list)) 

959 

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

961 

962 context += self.summary() 

963 

964 factor = 1.0 

965 cfactor = 1.0 

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

967 factor = 2. * math.pi 

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

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

970 

971 log = [] 

972 if self.normalization_factor is None \ 

973 or self.normalization_factor == 0.0: 

974 

975 log.append(( 

976 'warning', 

977 'No pole-zero normalization factor given. ' 

978 'Assuming a value of 1.0', 

979 context)) 

980 

981 nfactor = 1.0 

982 else: 

983 nfactor = self.normalization_factor 

984 

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

986 if not is_digital: 

987 resp = response.PoleZeroResponse( 

988 constant=nfactor*cfactor, 

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

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

991 else: 

992 if not deltat: 

993 log.append(( 

994 'error', 

995 'Digital response requires knowledge about sampling ' 

996 'interval. Response will be unusable.', 

997 context)) 

998 

999 resp = response.DigitalPoleZeroResponse( 

1000 constant=nfactor*cfactor, 

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

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

1003 deltat=deltat or 0.0) 

1004 

1005 if not self.normalization_frequency: 

1006 log.append(( 

1007 'warning', 

1008 'Cannot check pole-zero normalization factor: ' 

1009 'No normalization frequency given.', 

1010 context)) 

1011 

1012 else: 

1013 if is_digital and not deltat: 

1014 log.append(( 

1015 'warning', 

1016 'Cannot check computed vs reported normalization ' 

1017 'factor without knowing the sampling interval.', 

1018 context)) 

1019 else: 

1020 computed_normalization_factor = nfactor / abs(evaluate1( 

1021 resp, self.normalization_frequency.value)) 

1022 

1023 db = 20.0 * num.log10( 

1024 computed_normalization_factor / nfactor) 

1025 

1026 if abs(db) > 0.17: 

1027 log.append(( 

1028 'warning', 

1029 'Computed and reported normalization factors differ ' 

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

1031 db, 

1032 computed_normalization_factor, 

1033 nfactor), 

1034 context)) 

1035 

1036 return Delivery([resp], log) 

1037 

1038 

1039class ResponseListElement(Object): 

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

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

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

1043 

1044 

1045class Polynomial(BaseFilter): 

1046 ''' 

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

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

1049 stage of acquisition or a complete system. 

1050 ''' 

1051 

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

1053 xmltagname='ApproximationType') 

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

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

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

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

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

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

1060 

1061 def summary(self): 

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

1063 

1064 

1065class Decimation(Object): 

1066 ''' 

1067 Corresponds to SEED blockette 57. 

1068 ''' 

1069 

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

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

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

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

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

1075 

1076 def summary(self): 

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

1078 self.factor, 

1079 self.input_sample_rate.value, 

1080 self.input_sample_rate.value / self.factor, 

1081 self.delay.value) 

1082 

1083 def get_pyrocko_response(self): 

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

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

1086 else: 

1087 return Delivery([]) 

1088 

1089 

1090class Operator(Object): 

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

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

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

1094 

1095 

1096class Comment(Object): 

1097 ''' 

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

1099 and 59. 

1100 ''' 

1101 

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

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

1104 begin_effective_time = DummyAwareOptionalTimestamp.T( 

1105 optional=True, 

1106 xmltagname='BeginEffectiveTime') 

1107 end_effective_time = DummyAwareOptionalTimestamp.T( 

1108 optional=True, 

1109 xmltagname='EndEffectiveTime') 

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

1111 

1112 

1113class ResponseList(BaseFilter): 

1114 ''' 

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

1116 SEED blockette 55. 

1117 ''' 

1118 

1119 response_list_element_list = List.T( 

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

1121 

1122 def summary(self): 

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

1124 

1125 

1126class Log(Object): 

1127 ''' 

1128 Container for log entries. 

1129 ''' 

1130 

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

1132 

1133 

1134class ResponseStage(Object): 

1135 ''' 

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

1137 to 56. 

1138 ''' 

1139 

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

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

1142 poles_zeros_list = List.T( 

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

1144 coefficients_list = List.T( 

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

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

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

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

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

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

1151 

1152 def summary(self): 

1153 elements = [] 

1154 

1155 for stuff in [ 

1156 self.poles_zeros_list, self.coefficients_list, 

1157 self.response_list, self.fir, self.polynomial, 

1158 self.decimation, self.stage_gain]: 

1159 

1160 if stuff: 

1161 if isinstance(stuff, Object): 

1162 elements.append(stuff.summary()) 

1163 else: 

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

1165 

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

1167 self.number, 

1168 ', '.join(elements), 

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

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

1171 

1172 def get_squirrel_response_stage(self, context): 

1173 from pyrocko.squirrel.model import ResponseStage 

1174 

1175 delivery = Delivery() 

1176 delivery_pr = self.get_pyrocko_response(context) 

1177 log = delivery_pr.log 

1178 delivery_pr.log = [] 

1179 elements = delivery.extend_without_payload(delivery_pr) 

1180 

1181 delivery.payload = [ResponseStage( 

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

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

1184 input_sample_rate=self.input_sample_rate, 

1185 output_sample_rate=self.output_sample_rate, 

1186 elements=elements, 

1187 log=log)] 

1188 

1189 return delivery 

1190 

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

1192 

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

1194 

1195 responses = [] 

1196 log = [] 

1197 if self.stage_gain: 

1198 normalization_frequency = self.stage_gain.frequency or 0.0 

1199 else: 

1200 normalization_frequency = 0.0 

1201 

1202 if not gain_only: 

1203 deltat = None 

1204 delay_responses = [] 

1205 if self.decimation: 

1206 rate = self.decimation.input_sample_rate.value 

1207 if rate > 0.0: 

1208 deltat = 1.0 / rate 

1209 delivery = self.decimation.get_pyrocko_response() 

1210 if delivery.errors: 

1211 return delivery 

1212 

1213 delay_responses = delivery.payload 

1214 log.extend(delivery.log) 

1215 

1216 for pzs in self.poles_zeros_list: 

1217 delivery = pzs.get_pyrocko_response(context, deltat) 

1218 if delivery.errors: 

1219 return delivery 

1220 

1221 pz_resps = delivery.payload 

1222 log.extend(delivery.log) 

1223 responses.extend(pz_resps) 

1224 

1225 # emulate incorrect? evalresp behaviour 

1226 if pzs.normalization_frequency is not None and \ 

1227 pzs.normalization_frequency.value \ 

1228 != normalization_frequency \ 

1229 and normalization_frequency != 0.0: 

1230 

1231 try: 

1232 trial = response.MultiplyResponse(pz_resps) 

1233 anorm = num.abs(evaluate1( 

1234 trial, pzs.normalization_frequency.value)) 

1235 asens = num.abs( 

1236 evaluate1(trial, normalization_frequency)) 

1237 

1238 factor = anorm/asens 

1239 

1240 if abs(factor - 1.0) > 0.01: 

1241 log.append(( 

1242 'warning', 

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

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

1245 'possibly incorrect evalresp behaviour. ' 

1246 'Correction factor: %g' % ( 

1247 pzs.normalization_frequency.value, 

1248 normalization_frequency, 

1249 factor), 

1250 context)) 

1251 

1252 responses.append( 

1253 response.PoleZeroResponse(constant=factor)) 

1254 except response.InvalidResponseError as e: 

1255 log.append(( 

1256 'warning', 

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

1258 context)) 

1259 

1260 if len(self.poles_zeros_list) > 1: 

1261 log.append(( 

1262 'warning', 

1263 'Multiple poles and zeros records in single response ' 

1264 'stage.', 

1265 context)) 

1266 

1267 for cfs in self.coefficients_list + ( 

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

1269 

1270 delivery = cfs.get_pyrocko_response( 

1271 context, deltat, delay_responses, 

1272 normalization_frequency) 

1273 

1274 if delivery.errors: 

1275 return delivery 

1276 

1277 responses.extend(delivery.payload) 

1278 log.extend(delivery.log) 

1279 

1280 if len(self.coefficients_list) > 1: 

1281 log.append(( 

1282 'warning', 

1283 'Multiple filter coefficients lists in single response ' 

1284 'stage.', 

1285 context)) 

1286 

1287 if self.response_list: 

1288 log.append(( 

1289 'warning', 

1290 'Unhandled response element of type: ResponseList', 

1291 context)) 

1292 

1293 if self.polynomial: 

1294 log.append(( 

1295 'warning', 

1296 'Unhandled response element of type: Polynomial', 

1297 context)) 

1298 

1299 if self.stage_gain: 

1300 responses.append( 

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

1302 

1303 return Delivery(responses, log) 

1304 

1305 @property 

1306 def input_units(self): 

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

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

1309 if e is not None: 

1310 return e.input_units 

1311 

1312 return None 

1313 

1314 @property 

1315 def output_units(self): 

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

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

1318 if e is not None: 

1319 return e.output_units 

1320 

1321 return None 

1322 

1323 @property 

1324 def input_sample_rate(self): 

1325 if self.decimation: 

1326 return self.decimation.input_sample_rate.value 

1327 

1328 return None 

1329 

1330 @property 

1331 def output_sample_rate(self): 

1332 if self.decimation: 

1333 return self.decimation.input_sample_rate.value \ 

1334 / self.decimation.factor 

1335 

1336 return None 

1337 

1338 

1339class Response(Object): 

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

1341 instrument_sensitivity = Sensitivity.T(optional=True, 

1342 xmltagname='InstrumentSensitivity') 

1343 instrument_polynomial = Polynomial.T(optional=True, 

1344 xmltagname='InstrumentPolynomial') 

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

1346 

1347 @property 

1348 def output_sample_rate(self): 

1349 if self.stage_list: 

1350 return self.stage_list[-1].output_sample_rate 

1351 else: 

1352 return None 

1353 

1354 def check_sample_rates(self, channel): 

1355 

1356 if self.stage_list: 

1357 sample_rate = None 

1358 

1359 for stage in self.stage_list: 

1360 if stage.decimation: 

1361 input_sample_rate = \ 

1362 stage.decimation.input_sample_rate.value 

1363 

1364 if sample_rate is not None and not same_sample_rate( 

1365 sample_rate, input_sample_rate): 

1366 

1367 logger.warning( 

1368 'Response stage %i has unexpected input sample ' 

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

1370 stage.number, 

1371 input_sample_rate, 

1372 sample_rate)) 

1373 

1374 sample_rate = input_sample_rate / stage.decimation.factor 

1375 

1376 if sample_rate is not None and channel.sample_rate \ 

1377 and not same_sample_rate( 

1378 sample_rate, channel.sample_rate.value): 

1379 

1380 logger.warning( 

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

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

1383 channel.sample_rate.value, 

1384 sample_rate)) 

1385 

1386 def check_units(self): 

1387 

1388 if self.instrument_sensitivity \ 

1389 and self.instrument_sensitivity.input_units: 

1390 

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

1392 

1393 if self.stage_list: 

1394 for stage in self.stage_list: 

1395 if units and stage.input_units \ 

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

1397 

1398 logger.warning( 

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

1400 % ( 

1401 stage.number, 

1402 units, 

1403 'output units of stage %i' 

1404 if stage.number == 0 

1405 else 'sensitivity input units', 

1406 units)) 

1407 

1408 if stage.output_units: 

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

1410 else: 

1411 units = None 

1412 

1413 sout_units = self.instrument_sensitivity.output_units 

1414 if self.instrument_sensitivity and sout_units: 

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

1416 logger.warning( 

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

1418 % ( 

1419 stage.number, 

1420 units, 

1421 'sensitivity output units', 

1422 sout_units.name.upper())) 

1423 

1424 def _sensitivity_checkpoints(self, responses, context): 

1425 delivery = Delivery() 

1426 

1427 if self.instrument_sensitivity: 

1428 sval = self.instrument_sensitivity.value 

1429 sfreq = self.instrument_sensitivity.frequency 

1430 if sval is None: 

1431 delivery.log.append(( 

1432 'warning', 

1433 'No sensitivity value given.', 

1434 context)) 

1435 

1436 elif sval is None: 

1437 delivery.log.append(( 

1438 'warning', 

1439 'Reported sensitivity value is zero.', 

1440 context)) 

1441 

1442 elif sfreq is None: 

1443 delivery.log.append(( 

1444 'warning', 

1445 'Sensitivity frequency not given.', 

1446 context)) 

1447 

1448 else: 

1449 trial = response.MultiplyResponse(responses) 

1450 

1451 delivery.extend( 

1452 check_resp( 

1453 trial, sval, sfreq, 0.1, 

1454 'Instrument sensitivity value inconsistent with ' 

1455 'sensitivity computed from complete response.', 

1456 context)) 

1457 

1458 delivery.payload.append(response.FrequencyResponseCheckpoint( 

1459 frequency=sfreq, value=sval)) 

1460 

1461 return delivery 

1462 

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

1464 from pyrocko.squirrel.model import Response 

1465 

1466 if self.stage_list: 

1467 delivery = Delivery() 

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

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

1470 

1471 checkpoints = [] 

1472 if not delivery.errors: 

1473 all_responses = [] 

1474 for sq_stage in delivery.payload: 

1475 all_responses.extend(sq_stage.elements) 

1476 

1477 checkpoints.extend(delivery.extend_without_payload( 

1478 self._sensitivity_checkpoints(all_responses, context))) 

1479 

1480 sq_stages = delivery.payload 

1481 if sq_stages: 

1482 if sq_stages[0].input_quantity is None \ 

1483 and self.instrument_sensitivity is not None: 

1484 

1485 sq_stages[0].input_quantity = to_quantity( 

1486 self.instrument_sensitivity.input_units, 

1487 context, delivery) 

1488 sq_stages[-1].output_quantity = to_quantity( 

1489 self.instrument_sensitivity.output_units, 

1490 context, delivery) 

1491 

1492 sq_stages = delivery.expect() 

1493 

1494 return Response( 

1495 stages=sq_stages, 

1496 log=delivery.log, 

1497 checkpoints=checkpoints, 

1498 **kwargs) 

1499 

1500 elif self.instrument_sensitivity: 

1501 raise NoResponseInformation( 

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

1503 % context) 

1504 else: 

1505 raise NoResponseInformation( 

1506 'Empty instrument response (%s).' 

1507 % context) 

1508 

1509 def get_pyrocko_response( 

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

1511 

1512 delivery = Delivery() 

1513 if self.stage_list: 

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

1515 delivery.extend(stage.get_pyrocko_response( 

1516 context, gain_only=not ( 

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

1518 

1519 elif self.instrument_sensitivity: 

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

1521 

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

1523 checkpoints = delivery.extend_without_payload(delivery_cp) 

1524 if delivery.errors: 

1525 return delivery 

1526 

1527 if fake_input_units is not None: 

1528 if not self.instrument_sensitivity or \ 

1529 self.instrument_sensitivity.input_units is None: 

1530 

1531 delivery.errors.append(( 

1532 'NoResponseInformation', 

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

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

1535 context)) 

1536 

1537 return delivery 

1538 

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

1540 

1541 conresp = None 

1542 try: 

1543 conresp = conversion[ 

1544 fake_input_units.upper(), input_units] 

1545 

1546 except KeyError: 

1547 delivery.errors.append(( 

1548 'NoResponseInformation', 

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

1550 % (fake_input_units, input_units), 

1551 context)) 

1552 

1553 if conresp is not None: 

1554 delivery.payload.append(conresp) 

1555 for checkpoint in checkpoints: 

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

1557 conresp, checkpoint.frequency)) 

1558 

1559 delivery.payload = [ 

1560 response.MultiplyResponse( 

1561 delivery.payload, 

1562 checkpoints=checkpoints)] 

1563 

1564 return delivery 

1565 

1566 @classmethod 

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

1568 normalization_frequency=1.0): 

1569 ''' 

1570 Convert Pyrocko pole-zero response to StationXML response. 

1571 

1572 :param presponse: Pyrocko pole-zero response 

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

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

1575 response. 

1576 :type input_unit: str 

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

1578 response. 

1579 :type output_unit: str 

1580 :param normalization_frequency: Frequency where the normalization 

1581 factor for the StationXML response should be computed. 

1582 :type normalization_frequency: float 

1583 ''' 

1584 

1585 norm_factor = 1.0/float(abs( 

1586 evaluate1(presponse, normalization_frequency) 

1587 / presponse.constant)) 

1588 

1589 pzs = PolesZeros( 

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

1591 normalization_factor=norm_factor, 

1592 normalization_frequency=Frequency(normalization_frequency), 

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

1594 imaginary=FloatNoUnit(z.imag)) 

1595 for z in presponse.zeros], 

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

1597 imaginary=FloatNoUnit(z.imag)) 

1598 for z in presponse.poles]) 

1599 

1600 pzs.validate() 

1601 

1602 stage = ResponseStage( 

1603 number=1, 

1604 poles_zeros_list=[pzs], 

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

1606 

1607 resp = Response( 

1608 instrument_sensitivity=Sensitivity( 

1609 value=stage.stage_gain.value, 

1610 frequency=normalization_frequency, 

1611 input_units=Units(input_unit), 

1612 output_units=Units(output_unit)), 

1613 

1614 stage_list=[stage]) 

1615 

1616 return resp 

1617 

1618 

1619class BaseNode(Object): 

1620 ''' 

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

1622 ''' 

1623 

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

1625 start_date = DummyAwareOptionalTimestamp.T(optional=True, 

1626 xmlstyle='attribute') 

1627 end_date = DummyAwareOptionalTimestamp.T(optional=True, 

1628 xmlstyle='attribute') 

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

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

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

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

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

1634 

1635 def spans(self, *args): 

1636 if len(args) == 0: 

1637 return True 

1638 elif len(args) == 1: 

1639 return ((self.start_date is None or 

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

1641 (self.end_date is None or 

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

1643 

1644 elif len(args) == 2: 

1645 return ((self.start_date is None or 

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

1647 (self.end_date is None or 

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

1649 

1650 

1651def overlaps(a, b): 

1652 return ( 

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

1654 or a.start_date < b.end_date 

1655 ) and ( 

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

1657 or b.start_date < a.end_date 

1658 ) 

1659 

1660 

1661def check_overlaps(node_type_name, codes, nodes): 

1662 errors = [] 

1663 for ia, a in enumerate(nodes): 

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

1665 if overlaps(a, b): 

1666 errors.append( 

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

1668 ' %s - %s\n %s - %s' % ( 

1669 node_type_name, 

1670 '.'.join(codes), 

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

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

1673 

1674 return errors 

1675 

1676 

1677class Channel(BaseNode): 

1678 ''' 

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

1680 response blockettes. 

1681 ''' 

1682 

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

1684 external_reference_list = List.T( 

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

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

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

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

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

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

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

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

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

1694 sample_rate_ratio = SampleRateRatio.T(optional=True, 

1695 xmltagname='SampleRateRatio') 

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

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

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

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

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

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

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

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

1704 

1705 @property 

1706 def position_values(self): 

1707 lat = self.latitude.value 

1708 lon = self.longitude.value 

1709 elevation = value_or_none(self.elevation) 

1710 depth = value_or_none(self.depth) 

1711 return lat, lon, elevation, depth 

1712 

1713 

1714class Station(BaseNode): 

1715 ''' 

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

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

1718 epoch start and end dates. 

1719 ''' 

1720 

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

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

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

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

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

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

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

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

1729 creation_date = DummyAwareOptionalTimestamp.T( 

1730 optional=True, xmltagname='CreationDate') 

1731 termination_date = DummyAwareOptionalTimestamp.T( 

1732 optional=True, xmltagname='TerminationDate') 

1733 total_number_channels = Counter.T( 

1734 optional=True, xmltagname='TotalNumberChannels') 

1735 selected_number_channels = Counter.T( 

1736 optional=True, xmltagname='SelectedNumberChannels') 

1737 external_reference_list = List.T( 

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

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

1740 

1741 @property 

1742 def position_values(self): 

1743 lat = self.latitude.value 

1744 lon = self.longitude.value 

1745 elevation = value_or_none(self.elevation) 

1746 return lat, lon, elevation 

1747 

1748 

1749class Network(BaseNode): 

1750 ''' 

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

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

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

1754 contain 0 or more Stations. 

1755 ''' 

1756 

1757 total_number_stations = Counter.T(optional=True, 

1758 xmltagname='TotalNumberStations') 

1759 selected_number_stations = Counter.T(optional=True, 

1760 xmltagname='SelectedNumberStations') 

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

1762 

1763 @property 

1764 def station_code_list(self): 

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

1766 

1767 @property 

1768 def sl_code_list(self): 

1769 sls = set() 

1770 for station in self.station_list: 

1771 for channel in station.channel_list: 

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

1773 

1774 return sorted(sls) 

1775 

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

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

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

1779 if sls: 

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

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

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

1783 while ssls: 

1784 lines.append( 

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

1786 

1787 ssls[:n] = [] 

1788 

1789 return '\n'.join(lines) 

1790 

1791 

1792def value_or_none(x): 

1793 if x is not None: 

1794 return x.value 

1795 else: 

1796 return None 

1797 

1798 

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

1800 

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

1802 channels[0].position_values 

1803 

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

1805 info = '\n'.join( 

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

1807 x in channels) 

1808 

1809 mess = 'encountered inconsistencies in channel ' \ 

1810 'lat/lon/elevation/depth ' \ 

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

1812 

1813 if inconsistencies == 'raise': 

1814 raise InconsistentChannelLocations(mess) 

1815 

1816 elif inconsistencies == 'warn': 

1817 logger.warning(mess) 

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

1819 

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

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

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

1823 

1824 groups = {} 

1825 for channel in channels: 

1826 if channel.code not in groups: 

1827 groups[channel.code] = [] 

1828 

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

1830 

1831 pchannels = [] 

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

1833 data = [ 

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

1835 value_or_none(channel.dip)) 

1836 for channel in groups[code]] 

1837 

1838 azimuth, dip = util.consistency_merge( 

1839 data, 

1840 message='channel orientation values differ:', 

1841 error=inconsistencies) 

1842 

1843 pchannels.append( 

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

1845 

1846 return pyrocko.model.Station( 

1847 *nsl, 

1848 lat=mlat, 

1849 lon=mlon, 

1850 elevation=mele, 

1851 depth=mdep, 

1852 channels=pchannels) 

1853 

1854 

1855class FDSNStationXML(Object): 

1856 ''' 

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

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

1859 one or more Station containers. 

1860 ''' 

1861 

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

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

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

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

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

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

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

1869 

1870 xmltagname = 'FDSNStationXML' 

1871 guessable_xmlns = [guts_xmlns] 

1872 

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

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

1875 if self.created is None: 

1876 self.created = time.time() 

1877 

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

1879 time=None, timespan=None, 

1880 inconsistencies='warn'): 

1881 

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

1883 

1884 if nslcs is not None: 

1885 nslcs = set(nslcs) 

1886 

1887 if nsls is not None: 

1888 nsls = set(nsls) 

1889 

1890 tt = () 

1891 if time is not None: 

1892 tt = (time,) 

1893 elif timespan is not None: 

1894 tt = timespan 

1895 

1896 ns_have = set() 

1897 pstations = [] 

1898 sensor_to_channels = defaultdict(list) 

1899 for network in self.network_list: 

1900 if not network.spans(*tt): 

1901 continue 

1902 

1903 for station in network.station_list: 

1904 if not station.spans(*tt): 

1905 continue 

1906 

1907 if station.channel_list: 

1908 loc_to_channels = {} 

1909 for channel in station.channel_list: 

1910 if not channel.spans(*tt): 

1911 continue 

1912 

1913 loc = channel.location_code.strip() 

1914 if loc not in loc_to_channels: 

1915 loc_to_channels[loc] = [] 

1916 

1917 loc_to_channels[loc].append(channel) 

1918 

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

1920 channels = loc_to_channels[loc] 

1921 if nslcs is not None: 

1922 channels = [channel for channel in channels 

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

1924 channel.code) in nslcs] 

1925 

1926 if not channels: 

1927 continue 

1928 

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

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

1931 continue 

1932 

1933 for channel in channels: 

1934 k = nsl, channel.code[:-1], \ 

1935 channel.start_date, channel.end_date 

1936 sensor_to_channels[k].append(channel) 

1937 

1938 else: 

1939 ns = (network.code, station.code) 

1940 if ns in ns_have: 

1941 message = 'Duplicate station ' \ 

1942 '(multiple epochs match): %s.%s ' % ns 

1943 

1944 if inconsistencies == 'raise': 

1945 raise Inconsistencies(message) 

1946 else: 

1947 logger.warning(message) 

1948 

1949 ns_have.add(ns) 

1950 

1951 pstations.append(pyrocko.model.Station( 

1952 network.code, station.code, '*', 

1953 lat=station.latitude.value, 

1954 lon=station.longitude.value, 

1955 elevation=value_or_none(station.elevation), 

1956 name=station.description or '')) 

1957 

1958 sensor_have = set() 

1959 for (nsl, bi, _, _), channels in sensor_to_channels.items(): 

1960 if (nsl, bi) in sensor_have: 

1961 message = 'Duplicate station ' \ 

1962 '(multiple epochs match): %s.%s.%s' % nsl 

1963 

1964 if inconsistencies == 'raise': 

1965 raise Inconsistencies(message) 

1966 else: 

1967 logger.warning(message) 

1968 

1969 sensor_have.add((nsl, bi)) 

1970 

1971 pstations.append( 

1972 pyrocko_station_from_channels( 

1973 nsl, 

1974 channels, 

1975 inconsistencies=inconsistencies)) 

1976 

1977 return pstations 

1978 

1979 @classmethod 

1980 def from_pyrocko_stations( 

1981 cls, pyrocko_stations, add_flat_responses_from=None): 

1982 

1983 ''' 

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

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

1986 

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

1988 instances. 

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

1990 ''' 

1991 network_dict = defaultdict(list) 

1992 

1993 if add_flat_responses_from: 

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

1995 extra = dict( 

1996 response=Response( 

1997 instrument_sensitivity=Sensitivity( 

1998 value=1.0, 

1999 frequency=1.0, 

2000 input_units=Units(name=add_flat_responses_from)))) 

2001 else: 

2002 extra = {} 

2003 

2004 have_offsets = set() 

2005 for s in pyrocko_stations: 

2006 

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

2008 have_offsets.add(s.nsl()) 

2009 

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

2011 channel_list = [] 

2012 for c in s.channels: 

2013 channel_list.append( 

2014 Channel( 

2015 location_code=location, 

2016 code=c.name, 

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

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

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

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

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

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

2023 **extra 

2024 ) 

2025 ) 

2026 

2027 network_dict[network].append( 

2028 Station( 

2029 code=station, 

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

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

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

2033 channel_list=channel_list) 

2034 ) 

2035 

2036 if have_offsets: 

2037 logger.warning( 

2038 'StationXML does not support Cartesian offsets in ' 

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

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

2041 

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

2043 network_list = [] 

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

2045 

2046 network_list.append( 

2047 Network( 

2048 code=k, station_list=station_list, 

2049 total_number_stations=len(station_list))) 

2050 

2051 sxml = FDSNStationXML( 

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

2053 network_list=network_list) 

2054 

2055 sxml.validate() 

2056 return sxml 

2057 

2058 def iter_network_stations( 

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

2060 

2061 tt = () 

2062 if time is not None: 

2063 tt = (time,) 

2064 elif timespan is not None: 

2065 tt = timespan 

2066 

2067 for network in self.network_list: 

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

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

2070 continue 

2071 

2072 for station in network.station_list: 

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

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

2075 continue 

2076 

2077 yield (network, station) 

2078 

2079 def iter_network_station_channels( 

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

2081 time=None, timespan=None): 

2082 

2083 if loc is not None: 

2084 loc = loc.strip() 

2085 

2086 tt = () 

2087 if time is not None: 

2088 tt = (time,) 

2089 elif timespan is not None: 

2090 tt = timespan 

2091 

2092 for network in self.network_list: 

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

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

2095 continue 

2096 

2097 for station in network.station_list: 

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

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

2100 continue 

2101 

2102 if station.channel_list: 

2103 for channel in station.channel_list: 

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

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

2106 (loc is not None and 

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

2108 continue 

2109 

2110 yield (network, station, channel) 

2111 

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

2113 time=None, timespan=None): 

2114 

2115 groups = {} 

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

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

2118 

2119 net = network.code 

2120 sta = station.code 

2121 cha = channel.code 

2122 loc = channel.location_code.strip() 

2123 if len(cha) == 3: 

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

2125 elif len(cha) == 1: 

2126 bic = '' 

2127 else: 

2128 bic = cha 

2129 

2130 if channel.response and \ 

2131 channel.response.instrument_sensitivity and \ 

2132 channel.response.instrument_sensitivity.input_units: 

2133 

2134 unit = channel.response.instrument_sensitivity\ 

2135 .input_units.name.upper() 

2136 else: 

2137 unit = None 

2138 

2139 bic = (bic, unit) 

2140 

2141 k = net, sta, loc 

2142 if k not in groups: 

2143 groups[k] = {} 

2144 

2145 if bic not in groups[k]: 

2146 groups[k][bic] = [] 

2147 

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

2149 

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

2151 bad_bics = [] 

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

2153 sample_rates = [] 

2154 for channel in channels: 

2155 sample_rates.append(channel.sample_rate.value) 

2156 

2157 if not same(sample_rates): 

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

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

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

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

2162 

2163 logger.warning(err) 

2164 bad_bics.append(bic) 

2165 

2166 for bic in bad_bics: 

2167 del bic_to_channels[bic] 

2168 

2169 return groups 

2170 

2171 def choose_channels( 

2172 self, 

2173 target_sample_rate=None, 

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

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

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

2177 time=None, 

2178 timespan=None): 

2179 

2180 nslcs = {} 

2181 for nsl, bic_to_channels in self.get_channel_groups( 

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

2183 

2184 useful_bics = [] 

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

2186 rate = channels[0].sample_rate.value 

2187 

2188 if target_sample_rate is not None and \ 

2189 rate < target_sample_rate*0.99999: 

2190 continue 

2191 

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

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

2194 continue 

2195 

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

2197 continue 

2198 

2199 unit = bic[1] 

2200 

2201 prio_unit = len(priority_units) 

2202 try: 

2203 prio_unit = priority_units.index(unit) 

2204 except ValueError: 

2205 pass 

2206 

2207 prio_inst = len(priority_instrument_code) 

2208 prio_band = len(priority_band_code) 

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

2210 try: 

2211 prio_inst = priority_instrument_code.index( 

2212 channels[0].code[1]) 

2213 except ValueError: 

2214 pass 

2215 

2216 try: 

2217 prio_band = priority_band_code.index( 

2218 channels[0].code[0]) 

2219 except ValueError: 

2220 pass 

2221 

2222 if target_sample_rate is None: 

2223 rate = -rate 

2224 

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

2226 prio_inst, bic)) 

2227 

2228 useful_bics.sort() 

2229 

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

2231 channels = sorted( 

2232 bic_to_channels[bic], 

2233 key=lambda channel: channel.code) 

2234 

2235 if channels: 

2236 for channel in channels: 

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

2238 

2239 break 

2240 

2241 return nslcs 

2242 

2243 def get_pyrocko_response( 

2244 self, nslc, 

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

2246 

2247 net, sta, loc, cha = nslc 

2248 resps = [] 

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

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

2251 resp = channel.response 

2252 if resp: 

2253 resp.check_sample_rates(channel) 

2254 resp.check_units() 

2255 resps.append(resp.get_pyrocko_response( 

2256 '.'.join(nslc), 

2257 fake_input_units=fake_input_units, 

2258 stages=stages).expect_one()) 

2259 

2260 if not resps: 

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

2262 elif len(resps) > 1: 

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

2264 

2265 return resps[0] 

2266 

2267 @property 

2268 def n_code_list(self): 

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

2270 

2271 @property 

2272 def ns_code_list(self): 

2273 nss = set() 

2274 for network in self.network_list: 

2275 for station in network.station_list: 

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

2277 

2278 return sorted(nss) 

2279 

2280 @property 

2281 def nsl_code_list(self): 

2282 nsls = set() 

2283 for network in self.network_list: 

2284 for station in network.station_list: 

2285 for channel in station.channel_list: 

2286 nsls.add( 

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

2288 

2289 return sorted(nsls) 

2290 

2291 @property 

2292 def nslc_code_list(self): 

2293 nslcs = set() 

2294 for network in self.network_list: 

2295 for station in network.station_list: 

2296 for channel in station.channel_list: 

2297 nslcs.add( 

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

2299 channel.code)) 

2300 

2301 return sorted(nslcs) 

2302 

2303 def summary(self): 

2304 lst = [ 

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

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

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

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

2309 ] 

2310 return '\n'.join(lst) 

2311 

2312 def summary_stages(self): 

2313 data = [] 

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

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

2316 channel.code) 

2317 

2318 stages = [] 

2319 in_units = '?' 

2320 out_units = '?' 

2321 if channel.response: 

2322 sens = channel.response.instrument_sensitivity 

2323 if sens: 

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

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

2326 

2327 for stage in channel.response.stage_list: 

2328 stages.append(stage.summary()) 

2329 

2330 data.append( 

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

2332 in_units, out_units, stages)) 

2333 

2334 data.sort() 

2335 

2336 lst = [] 

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

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

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

2340 for stage in stages: 

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

2342 

2343 return '\n'.join(lst) 

2344 

2345 def _check_overlaps(self): 

2346 by_nslc = {} 

2347 for network in self.network_list: 

2348 for station in network.station_list: 

2349 for channel in station.channel_list: 

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

2351 channel.code) 

2352 if nslc not in by_nslc: 

2353 by_nslc[nslc] = [] 

2354 

2355 by_nslc[nslc].append(channel) 

2356 

2357 errors = [] 

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

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

2360 

2361 return errors 

2362 

2363 def check(self): 

2364 errors = [] 

2365 for meth in [self._check_overlaps]: 

2366 errors.extend(meth()) 

2367 

2368 if errors: 

2369 raise Inconsistencies( 

2370 'Inconsistencies found in StationXML:\n ' 

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

2372 

2373 

2374def load_channel_table(stream): 

2375 

2376 networks = {} 

2377 stations = {} 

2378 

2379 for line in stream: 

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

2381 if line.startswith('#'): 

2382 continue 

2383 

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

2385 

2386 if len(t) != 17: 

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

2388 continue 

2389 

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

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

2392 

2393 try: 

2394 scale = float(scale) 

2395 except ValueError: 

2396 scale = None 

2397 

2398 try: 

2399 scale_freq = float(scale_freq) 

2400 except ValueError: 

2401 scale_freq = None 

2402 

2403 try: 

2404 depth = float(dep) 

2405 except ValueError: 

2406 depth = 0.0 

2407 

2408 try: 

2409 azi = float(azi) 

2410 dip = float(dip) 

2411 except ValueError: 

2412 azi = None 

2413 dip = None 

2414 

2415 try: 

2416 if net not in networks: 

2417 network = Network(code=net) 

2418 else: 

2419 network = networks[net] 

2420 

2421 if (net, sta) not in stations: 

2422 station = Station( 

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

2424 

2425 station.regularize() 

2426 else: 

2427 station = stations[net, sta] 

2428 

2429 if scale: 

2430 resp = Response( 

2431 instrument_sensitivity=Sensitivity( 

2432 value=scale, 

2433 frequency=scale_freq, 

2434 input_units=scale_units)) 

2435 else: 

2436 resp = None 

2437 

2438 channel = Channel( 

2439 code=cha, 

2440 location_code=loc.strip(), 

2441 latitude=lat, 

2442 longitude=lon, 

2443 elevation=ele, 

2444 depth=depth, 

2445 azimuth=azi, 

2446 dip=dip, 

2447 sensor=Equipment(description=sens), 

2448 response=resp, 

2449 sample_rate=sample_rate, 

2450 start_date=start_date, 

2451 end_date=end_date or None) 

2452 

2453 channel.regularize() 

2454 

2455 except ValidationError: 

2456 raise InvalidRecord(line) 

2457 

2458 if net not in networks: 

2459 networks[net] = network 

2460 

2461 if (net, sta) not in stations: 

2462 stations[net, sta] = station 

2463 network.station_list.append(station) 

2464 

2465 station.channel_list.append(channel) 

2466 

2467 return FDSNStationXML( 

2468 source='created from table input', 

2469 created=time.time(), 

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

2471 

2472 

2473def primitive_merge(sxs): 

2474 networks = [] 

2475 for sx in sxs: 

2476 networks.extend(sx.network_list) 

2477 

2478 return FDSNStationXML( 

2479 source='merged from different sources', 

2480 created=time.time(), 

2481 network_list=copy.deepcopy( 

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

2483 

2484 

2485def split_channels(sx): 

2486 for nslc in sx.nslc_code_list: 

2487 network_list = sx.network_list 

2488 network_list_filtered = [ 

2489 network for network in network_list 

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

2491 

2492 for network in network_list_filtered: 

2493 sx.network_list = [network] 

2494 station_list = network.station_list 

2495 station_list_filtered = [ 

2496 station for station in station_list 

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

2498 

2499 for station in station_list_filtered: 

2500 network.station_list = [station] 

2501 channel_list = station.channel_list 

2502 station.channel_list = [ 

2503 channel for channel in channel_list 

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

2505 

2506 yield nslc, copy.deepcopy(sx) 

2507 

2508 station.channel_list = channel_list 

2509 

2510 network.station_list = station_list 

2511 

2512 sx.network_list = network_list 

2513 

2514 

2515if __name__ == '__main__': 

2516 from optparse import OptionParser 

2517 

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

2519 

2520 usage = \ 

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

2522 '<filename> [options]' 

2523 

2524 description = '''Torture StationXML file.''' 

2525 

2526 parser = OptionParser( 

2527 usage=usage, 

2528 description=description, 

2529 formatter=util.BetterHelpFormatter()) 

2530 

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

2532 

2533 if len(args) != 2: 

2534 parser.print_help() 

2535 sys.exit(1) 

2536 

2537 action, path = args 

2538 

2539 sx = load_xml(filename=path) 

2540 if action == 'check': 

2541 try: 

2542 sx.check() 

2543 except Inconsistencies as e: 

2544 logger.error(e) 

2545 sys.exit(1) 

2546 

2547 elif action == 'stats': 

2548 print(sx.summary()) 

2549 

2550 elif action == 'stages': 

2551 print(sx.summary_stages()) 

2552 

2553 else: 

2554 parser.print_help() 

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