1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import sys 

7import time 

8import logging 

9import datetime 

10import calendar 

11import math 

12import copy 

13 

14import numpy as num 

15 

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

17 Unicode, Int, Float, List, Object, Timestamp, 

18 ValidationError, TBase, re_tz, Any, Tuple) 

19from pyrocko.guts import load_xml # noqa 

20from pyrocko.util import hpfloat, time_to_str, get_time_float 

21 

22import pyrocko.model 

23from pyrocko import util, response 

24 

25guts_prefix = 'sx' 

26 

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

28 

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

30 

31conversion = { 

32 ('M', 'M'): None, 

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

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

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

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

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

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

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

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

41 

42 

43unit_to_quantity = { 

44 'M/S': 'velocity', 

45 'M': 'displacement', 

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

47 'V': 'voltage', 

48 'COUNTS': 'counts', 

49 'COUNT': 'counts', 

50 'PA': 'pressure', 

51 'RAD': 'rotation-displacement', 

52 'R': 'rotation-displacement', 

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

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

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

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

57 

58 

59def to_quantity(unit, context, delivery): 

60 

61 if unit is None: 

62 return None 

63 

64 name = unit.name.upper() 

65 if name in unit_to_quantity: 

66 return unit_to_quantity[name] 

67 else: 

68 delivery.log.append(( 

69 'warning', 

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

71 unit.name.upper() + ( 

72 ' (%s)' % unit.description 

73 if unit.description else '')), 

74 context)) 

75 

76 return 'unsupported_quantity(%s)' % unit 

77 

78 

79class StationXMLError(Exception): 

80 pass 

81 

82 

83class Inconsistencies(StationXMLError): 

84 pass 

85 

86 

87class NoResponseInformation(StationXMLError): 

88 pass 

89 

90 

91class MultipleResponseInformation(StationXMLError): 

92 pass 

93 

94 

95class InconsistentResponseInformation(StationXMLError): 

96 pass 

97 

98 

99class InconsistentChannelLocations(StationXMLError): 

100 pass 

101 

102 

103class InvalidRecord(StationXMLError): 

104 def __init__(self, line): 

105 StationXMLError.__init__(self) 

106 self._line = line 

107 

108 def __str__(self): 

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

110 

111 

112_exceptions = dict( 

113 Inconsistencies=Inconsistencies, 

114 NoResponseInformation=NoResponseInformation, 

115 MultipleResponseInformation=MultipleResponseInformation, 

116 InconsistentResponseInformation=InconsistentResponseInformation, 

117 InconsistentChannelLocations=InconsistentChannelLocations, 

118 InvalidRecord=InvalidRecord, 

119 ValueError=ValueError) 

120 

121 

122_logs = dict( 

123 info=logger.info, 

124 warning=logger.warning, 

125 error=logger.error) 

126 

127 

128class DeliveryError(StationXMLError): 

129 pass 

130 

131 

132class Delivery(Object): 

133 

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

135 if payload is None: 

136 payload = [] 

137 

138 if log is None: 

139 log = [] 

140 

141 if errors is None: 

142 errors = [] 

143 

144 if error is not None: 

145 errors.append(error) 

146 

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

148 

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

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

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

152 

153 def extend(self, other): 

154 self.payload.extend(other.payload) 

155 self.log.extend(other.log) 

156 self.errors.extend(other.errors) 

157 

158 def extend_without_payload(self, other): 

159 self.log.extend(other.log) 

160 self.errors.extend(other.errors) 

161 return other.payload 

162 

163 def emit_log(self): 

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

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

166 _logs[name](message) 

167 

168 def expect(self, quiet=False): 

169 if not quiet: 

170 self.emit_log() 

171 

172 if self.errors: 

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

174 if context: 

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

176 

177 if len(self.errors) > 1: 

178 message += ' Additional errors pending.' 

179 

180 raise _exceptions[name](message) 

181 

182 return self.payload 

183 

184 def expect_one(self, quiet=False): 

185 payload = self.expect(quiet=quiet) 

186 if len(payload) != 1: 

187 raise DeliveryError( 

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

189 

190 return payload[0] 

191 

192 

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

194 words = s.split() 

195 lines = [] 

196 t = [] 

197 n = 0 

198 for w in words: 

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

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

201 n = indent 

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

203 

204 t.append(w) 

205 n += len(w) + 1 

206 

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

208 return '\n'.join(lines) 

209 

210 

211def same(x, eps=0.0): 

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

213 return False 

214 

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

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

217 else: 

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

219 

220 

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

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

223 

224 

225def evaluate1(resp, f): 

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

227 

228 

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

230 

231 try: 

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

233 except response.InvalidResponseError as e: 

234 return Delivery(log=[( 

235 'warning', 

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

237 context)]) 

238 

239 if value_resp == 0.0: 

240 return Delivery(log=[( 

241 'warning', 

242 '%s\n' 

243 ' computed response is zero' % prelude, 

244 context)]) 

245 

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

247 

248 if num.abs(diff_db) > limit_db: 

249 return Delivery(log=[( 

250 'warning', 

251 '%s\n' 

252 ' reported value: %g\n' 

253 ' computed value: %g\n' 

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

255 ' factor: %g\n' 

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

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

258 prelude, 

259 value, 

260 value_resp, 

261 frequency, 

262 value_resp/value, 

263 diff_db, 

264 limit_db), 

265 context)]) 

266 

267 return Delivery() 

268 

269 

270def tts(t): 

271 if t is None: 

272 return '?' 

273 else: 

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

275 

276 

277def le_open_left(a, b): 

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

279 

280 

281def le_open_right(a, b): 

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

283 

284 

285def eq_open(a, b): 

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

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

288 

289 

290def contains(a, b): 

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

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

293 

294 

295def find_containing(candidates, node): 

296 for candidate in candidates: 

297 if contains(candidate, node): 

298 return candidate 

299 

300 return None 

301 

302 

303this_year = time.gmtime()[0] 

304 

305 

306class DummyAwareOptionalTimestamp(Object): 

307 ''' 

308 Optional timestamp with support for some common placeholder values. 

309 

310 Some StationXML files contain arbitrary placeholder values for open end 

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

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

313 this type. 

314 ''' 

315 dummy_for = (hpfloat, float) 

316 dummy_for_description = 'time_float' 

317 

318 class __T(TBase): 

319 

320 def regularize_extra(self, val): 

321 time_float = get_time_float() 

322 

323 if isinstance(val, datetime.datetime): 

324 tt = val.utctimetuple() 

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

326 

327 elif isinstance(val, datetime.date): 

328 tt = val.timetuple() 

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

330 

331 elif isinstance(val, str): 

332 val = val.strip() 

333 

334 tz_offset = 0 

335 

336 m = re_tz.search(val) 

337 if m: 

338 sh = m.group(2) 

339 sm = m.group(4) 

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

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

342 

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

344 

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

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

347 

348 try: 

349 val = util.str_to_time(val) - tz_offset 

350 

351 except util.TimeStrError: 

352 year = int(val[:4]) 

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

354 if year > this_year + 100: 

355 return None # StationXML contained a dummy date 

356 

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

358 return None 

359 

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

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

362 return None # StationXML contained a dummy date 

363 

364 raise 

365 

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

367 val = time_float(val) 

368 

369 else: 

370 raise ValidationError( 

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

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

373 

374 return val 

375 

376 def to_save(self, val): 

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

378 .rstrip('0').rstrip('.') 

379 

380 def to_save_xml(self, val): 

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

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

383 

384 

385class Nominal(StringChoice): 

386 choices = [ 

387 'NOMINAL', 

388 'CALCULATED'] 

389 

390 

391class Email(UnicodePattern): 

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

393 

394 

395class RestrictedStatus(StringChoice): 

396 choices = [ 

397 'open', 

398 'closed', 

399 'partial'] 

400 

401 

402class Type(StringChoice): 

403 choices = [ 

404 'TRIGGERED', 

405 'CONTINUOUS', 

406 'HEALTH', 

407 'GEOPHYSICAL', 

408 'WEATHER', 

409 'FLAG', 

410 'SYNTHESIZED', 

411 'INPUT', 

412 'EXPERIMENTAL', 

413 'MAINTENANCE', 

414 'BEAM'] 

415 

416 class __T(StringChoice.T): 

417 def validate_extra(self, val): 

418 if val not in self.choices: 

419 logger.warning( 

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

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

422 

423 

424class PzTransferFunction(StringChoice): 

425 choices = [ 

426 'LAPLACE (RADIANS/SECOND)', 

427 'LAPLACE (HERTZ)', 

428 'DIGITAL (Z-TRANSFORM)'] 

429 

430 

431class Symmetry(StringChoice): 

432 choices = [ 

433 'NONE', 

434 'EVEN', 

435 'ODD'] 

436 

437 

438class CfTransferFunction(StringChoice): 

439 

440 class __T(StringChoice.T): 

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

442 if regularize: 

443 try: 

444 val = str(val) 

445 except ValueError: 

446 raise ValidationError( 

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

448 repr(val))) 

449 

450 val = self.dummy_cls.replacements.get(val, val) 

451 

452 self.validate_extra(val) 

453 return val 

454 

455 choices = [ 

456 'ANALOG (RADIANS/SECOND)', 

457 'ANALOG (HERTZ)', 

458 'DIGITAL'] 

459 

460 replacements = { 

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

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

463 } 

464 

465 

466class Approximation(StringChoice): 

467 choices = [ 

468 'MACLAURIN'] 

469 

470 

471class PhoneNumber(StringPattern): 

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

473 

474 

475class Site(Object): 

476 ''' 

477 Description of a site location using name and optional geopolitical 

478 boundaries (country, city, etc.). 

479 ''' 

480 

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

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

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

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

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

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

487 

488 

489class ExternalReference(Object): 

490 ''' 

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

492 want to reference in StationXML. 

493 ''' 

494 

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

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

497 

498 

499class Units(Object): 

500 ''' 

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

502 ''' 

503 

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

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

506 

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

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

509 

510 

511class Counter(Int): 

512 pass 

513 

514 

515class SampleRateRatio(Object): 

516 ''' 

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

518 ''' 

519 

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

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

522 

523 

524class Gain(Object): 

525 ''' 

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

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

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

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

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

531 ''' 

532 

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

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

535 

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

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

538 

539 def summary(self): 

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

541 

542 

543class NumeratorCoefficient(Object): 

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

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

546 

547 

548class FloatNoUnit(Object): 

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

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

551 

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

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

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

555 

556 

557class FloatWithUnit(FloatNoUnit): 

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

559 

560 

561class Equipment(Object): 

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

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

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

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

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

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

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

569 installation_date = DummyAwareOptionalTimestamp.T( 

570 optional=True, 

571 xmltagname='InstallationDate') 

572 removal_date = DummyAwareOptionalTimestamp.T( 

573 optional=True, 

574 xmltagname='RemovalDate') 

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

576 

577 

578class PhoneNumber(Object): 

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

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

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

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

583 

584 

585class BaseFilter(Object): 

586 ''' 

587 The BaseFilter is derived by all filters. 

588 ''' 

589 

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

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

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

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

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

595 

596 

597class Sensitivity(Gain): 

598 ''' 

599 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional 

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

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

602 decibels specified in FrequencyDBVariation. 

603 ''' 

604 

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

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

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

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

609 frequency_db_variation = Float.T(optional=True, 

610 xmltagname='FrequencyDBVariation') 

611 

612 def get_pyrocko_response(self): 

613 return Delivery( 

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

615 

616 

617class Coefficient(FloatNoUnit): 

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

619 

620 

621class PoleZero(Object): 

622 ''' 

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

624 ''' 

625 

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

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

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

629 

630 def value(self): 

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

632 

633 

634class ClockDrift(FloatWithUnit): 

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

636 xmlstyle='attribute') # fixed 

637 

638 

639class Second(FloatWithUnit): 

640 ''' 

641 A time value in seconds. 

642 ''' 

643 

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

645 # fixed unit 

646 

647 

648class Voltage(FloatWithUnit): 

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

650 # fixed unit 

651 

652 

653class Angle(FloatWithUnit): 

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

655 # fixed unit 

656 

657 

658class Azimuth(FloatWithUnit): 

659 ''' 

660 Instrument azimuth, degrees clockwise from North. 

661 ''' 

662 

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

664 # fixed unit 

665 

666 

667class Dip(FloatWithUnit): 

668 ''' 

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

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

671 ''' 

672 

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

674 # fixed unit 

675 

676 

677class Distance(FloatWithUnit): 

678 ''' 

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

680 ''' 

681 

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

683 # NOT fixed unit! 

684 

685 

686class Frequency(FloatWithUnit): 

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

688 # fixed unit 

689 

690 

691class SampleRate(FloatWithUnit): 

692 ''' 

693 Sample rate in samples per second. 

694 ''' 

695 

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

697 # fixed unit 

698 

699 

700class Person(Object): 

701 ''' 

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

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

704 ''' 

705 

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

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

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

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

710 

711 

712class FIR(BaseFilter): 

713 ''' 

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

715 also commonly documented using the Coefficients element. 

716 ''' 

717 

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

719 numerator_coefficient_list = List.T( 

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

721 

722 def summary(self): 

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

724 self.get_ncoefficiencs(), 

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

726 

727 def get_effective_coefficients(self): 

728 b = num.array( 

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

730 dtype=float) 

731 

732 if self.symmetry == 'ODD': 

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

734 elif self.symmetry == 'EVEN': 

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

736 

737 return b 

738 

739 def get_effective_symmetry(self): 

740 if self.symmetry == 'NONE': 

741 b = self.get_effective_coefficients() 

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

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

744 

745 return self.symmetry 

746 

747 def get_ncoefficiencs(self): 

748 nf = len(self.numerator_coefficient_list) 

749 if self.symmetry == 'ODD': 

750 nc = nf * 2 + 1 

751 elif self.symmetry == 'EVEN': 

752 nc = nf * 2 

753 else: 

754 nc = nf 

755 

756 return nc 

757 

758 def estimate_delay(self, deltat): 

759 nc = self.get_ncoefficiencs() 

760 if nc > 0: 

761 return deltat * (nc - 1) / 2.0 

762 else: 

763 return 0.0 

764 

765 def get_pyrocko_response( 

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

767 

768 context += self.summary() 

769 

770 if not self.numerator_coefficient_list: 

771 return Delivery([]) 

772 

773 b = self.get_effective_coefficients() 

774 

775 log = [] 

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

777 

778 if not deltat: 

779 log.append(( 

780 'error', 

781 'Digital response requires knowledge about sampling ' 

782 'interval. Response will be unusable.', 

783 context)) 

784 

785 resp = response.DigitalFilterResponse( 

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

787 

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

789 normalization_frequency = 0.0 

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

791 

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

793 log.append(( 

794 'warning', 

795 'FIR filter coefficients are not normalized. Normalizing ' 

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

797 

798 resp = response.DigitalFilterResponse( 

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

800 drop_phase=drop_phase) 

801 

802 resps = [resp] 

803 

804 if not drop_phase: 

805 resps.extend(delay_responses) 

806 

807 return Delivery(resps, log=log) 

808 

809 

810class Coefficients(BaseFilter): 

811 ''' 

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

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

814 instead. Corresponds to SEED blockette 54. 

815 ''' 

816 

817 cf_transfer_function_type = CfTransferFunction.T( 

818 xmltagname='CfTransferFunctionType') 

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

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

821 

822 def summary(self): 

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

824 'ABC?'[ 

825 CfTransferFunction.choices.index( 

826 self.cf_transfer_function_type)], 

827 len(self.numerator_list), 

828 len(self.denominator_list), 

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

830 

831 def estimate_delay(self, deltat): 

832 nc = len(self.numerator_list) 

833 if nc > 0: 

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

835 else: 

836 return 0.0 

837 

838 def is_symmetric_fir(self): 

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

840 return False 

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

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

843 

844 def get_pyrocko_response( 

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

846 

847 context += self.summary() 

848 

849 factor = 1.0 

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

851 factor = 2.0*math.pi 

852 

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

854 return Delivery(payload=[]) 

855 

856 b = num.array( 

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

858 

859 a = num.array( 

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

861 dtype=float) 

862 

863 log = [] 

864 resps = [] 

865 if self.cf_transfer_function_type in [ 

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

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

868 

869 elif self.cf_transfer_function_type == 'DIGITAL': 

870 if not deltat: 

871 log.append(( 

872 'error', 

873 'Digital response requires knowledge about sampling ' 

874 'interval. Response will be unusable.', 

875 context)) 

876 

877 drop_phase = self.is_symmetric_fir() 

878 resp = response.DigitalFilterResponse( 

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

880 

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

882 normalization = num.abs(evaluate1( 

883 resp, normalization_frequency)) 

884 

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

886 log.append(( 

887 'warning', 

888 'FIR filter coefficients are not normalized. ' 

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

890 context)) 

891 

892 resp = response.DigitalFilterResponse( 

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

894 drop_phase=drop_phase) 

895 

896 resps.append(resp) 

897 

898 if not drop_phase: 

899 resps.extend(delay_responses) 

900 

901 else: 

902 return Delivery(error=( 

903 'ValueError', 

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

905 self.cf_transfer_function_type))) 

906 

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

908 

909 

910class Latitude(FloatWithUnit): 

911 ''' 

912 Type for latitude coordinate. 

913 ''' 

914 

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

916 # fixed unit 

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

918 

919 

920class Longitude(FloatWithUnit): 

921 ''' 

922 Type for longitude coordinate. 

923 ''' 

924 

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

926 # fixed unit 

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

928 

929 

930class PolesZeros(BaseFilter): 

931 ''' 

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

933 ''' 

934 

935 pz_transfer_function_type = PzTransferFunction.T( 

936 xmltagname='PzTransferFunctionType') 

937 normalization_factor = Float.T(default=1.0, 

938 xmltagname='NormalizationFactor') 

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

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

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

942 

943 def summary(self): 

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

945 'ABC?'[ 

946 PzTransferFunction.choices.index( 

947 self.pz_transfer_function_type)], 

948 len(self.pole_list), 

949 len(self.zero_list)) 

950 

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

952 

953 context += self.summary() 

954 

955 factor = 1.0 

956 cfactor = 1.0 

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

958 factor = 2. * math.pi 

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

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

961 

962 log = [] 

963 if self.normalization_factor is None \ 

964 or self.normalization_factor == 0.0: 

965 

966 log.append(( 

967 'warning', 

968 'No pole-zero normalization factor given. ' 

969 'Assuming a value of 1.0', 

970 context)) 

971 

972 nfactor = 1.0 

973 else: 

974 nfactor = self.normalization_factor 

975 

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

977 if not is_digital: 

978 resp = response.PoleZeroResponse( 

979 constant=nfactor*cfactor, 

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

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

982 else: 

983 if not deltat: 

984 log.append(( 

985 'error', 

986 'Digital response requires knowledge about sampling ' 

987 'interval. Response will be unusable.', 

988 context)) 

989 

990 resp = response.DigitalPoleZeroResponse( 

991 constant=nfactor*cfactor, 

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

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

994 deltat=deltat or 0.0) 

995 

996 if not self.normalization_frequency.value: 

997 log.append(( 

998 'warning', 

999 'Cannot check pole-zero normalization factor: ' 

1000 'No normalization frequency given.', 

1001 context)) 

1002 

1003 else: 

1004 if is_digital and not deltat: 

1005 log.append(( 

1006 'warning', 

1007 'Cannot check computed vs reported normalization ' 

1008 'factor without knowing the sampling interval.', 

1009 context)) 

1010 else: 

1011 computed_normalization_factor = nfactor / abs(evaluate1( 

1012 resp, self.normalization_frequency.value)) 

1013 

1014 db = 20.0 * num.log10( 

1015 computed_normalization_factor / nfactor) 

1016 

1017 if abs(db) > 0.17: 

1018 log.append(( 

1019 'warning', 

1020 'Computed and reported normalization factors differ ' 

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

1022 db, 

1023 computed_normalization_factor, 

1024 nfactor), 

1025 context)) 

1026 

1027 return Delivery([resp], log) 

1028 

1029 

1030class ResponseListElement(Object): 

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

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

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

1034 

1035 

1036class Polynomial(BaseFilter): 

1037 ''' 

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

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

1040 stage of acquisition or a complete system. 

1041 ''' 

1042 

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

1044 xmltagname='ApproximationType') 

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

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

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

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

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

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

1051 

1052 def summary(self): 

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

1054 

1055 

1056class Decimation(Object): 

1057 ''' 

1058 Corresponds to SEED blockette 57. 

1059 ''' 

1060 

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

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

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

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

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

1066 

1067 def summary(self): 

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

1069 self.factor, 

1070 self.input_sample_rate.value, 

1071 self.input_sample_rate.value / self.factor, 

1072 self.delay.value) 

1073 

1074 def get_pyrocko_response(self): 

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

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

1077 else: 

1078 return Delivery([]) 

1079 

1080 

1081class Operator(Object): 

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

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

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

1085 

1086 

1087class Comment(Object): 

1088 ''' 

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

1090 and 59. 

1091 ''' 

1092 

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

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

1095 begin_effective_time = DummyAwareOptionalTimestamp.T( 

1096 optional=True, 

1097 xmltagname='BeginEffectiveTime') 

1098 end_effective_time = DummyAwareOptionalTimestamp.T( 

1099 optional=True, 

1100 xmltagname='EndEffectiveTime') 

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

1102 

1103 

1104class ResponseList(BaseFilter): 

1105 ''' 

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

1107 SEED blockette 55. 

1108 ''' 

1109 

1110 response_list_element_list = List.T( 

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

1112 

1113 def summary(self): 

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

1115 

1116 

1117class Log(Object): 

1118 ''' 

1119 Container for log entries. 

1120 ''' 

1121 

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

1123 

1124 

1125class ResponseStage(Object): 

1126 ''' 

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

1128 to 56. 

1129 ''' 

1130 

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

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

1133 poles_zeros_list = List.T( 

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

1135 coefficients_list = List.T( 

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

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

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

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

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

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

1142 

1143 def summary(self): 

1144 elements = [] 

1145 

1146 for stuff in [ 

1147 self.poles_zeros_list, self.coefficients_list, 

1148 self.response_list, self.fir, self.polynomial, 

1149 self.decimation, self.stage_gain]: 

1150 

1151 if stuff: 

1152 if isinstance(stuff, Object): 

1153 elements.append(stuff.summary()) 

1154 else: 

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

1156 

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

1158 self.number, 

1159 ', '.join(elements), 

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

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

1162 

1163 def get_squirrel_response_stage(self, context): 

1164 from pyrocko.squirrel.model import ResponseStage 

1165 

1166 delivery = Delivery() 

1167 delivery_pr = self.get_pyrocko_response(context) 

1168 log = delivery_pr.log 

1169 delivery_pr.log = [] 

1170 elements = delivery.extend_without_payload(delivery_pr) 

1171 

1172 delivery.payload = [ResponseStage( 

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

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

1175 input_sample_rate=self.input_sample_rate, 

1176 output_sample_rate=self.output_sample_rate, 

1177 elements=elements, 

1178 log=log)] 

1179 

1180 return delivery 

1181 

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

1183 

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

1185 

1186 responses = [] 

1187 log = [] 

1188 if self.stage_gain: 

1189 normalization_frequency = self.stage_gain.frequency or 0.0 

1190 else: 

1191 normalization_frequency = 0.0 

1192 

1193 if not gain_only: 

1194 deltat = None 

1195 delay_responses = [] 

1196 if self.decimation: 

1197 rate = self.decimation.input_sample_rate.value 

1198 if rate > 0.0: 

1199 deltat = 1.0 / rate 

1200 delivery = self.decimation.get_pyrocko_response() 

1201 if delivery.errors: 

1202 return delivery 

1203 

1204 delay_responses = delivery.payload 

1205 log.extend(delivery.log) 

1206 

1207 for pzs in self.poles_zeros_list: 

1208 delivery = pzs.get_pyrocko_response(context, deltat) 

1209 if delivery.errors: 

1210 return delivery 

1211 

1212 pz_resps = delivery.payload 

1213 log.extend(delivery.log) 

1214 responses.extend(pz_resps) 

1215 

1216 # emulate incorrect? evalresp behaviour 

1217 if pzs.normalization_frequency != normalization_frequency \ 

1218 and normalization_frequency != 0.0: 

1219 

1220 try: 

1221 trial = response.MultiplyResponse(pz_resps) 

1222 anorm = num.abs(evaluate1( 

1223 trial, pzs.normalization_frequency.value)) 

1224 asens = num.abs( 

1225 evaluate1(trial, normalization_frequency)) 

1226 

1227 factor = anorm/asens 

1228 

1229 if abs(factor - 1.0) > 0.01: 

1230 log.append(( 

1231 'warning', 

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

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

1234 'possibly incorrect evalresp behaviour. ' 

1235 'Correction factor: %g' % ( 

1236 pzs.normalization_frequency.value, 

1237 normalization_frequency, 

1238 factor), 

1239 context)) 

1240 

1241 responses.append( 

1242 response.PoleZeroResponse(constant=factor)) 

1243 except response.InvalidResponseError as e: 

1244 log.append(( 

1245 'warning', 

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

1247 context)) 

1248 

1249 if len(self.poles_zeros_list) > 1: 

1250 log.append(( 

1251 'warning', 

1252 'Multiple poles and zeros records in single response ' 

1253 'stage.', 

1254 context)) 

1255 

1256 for cfs in self.coefficients_list + ( 

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

1258 

1259 delivery = cfs.get_pyrocko_response( 

1260 context, deltat, delay_responses, 

1261 normalization_frequency) 

1262 

1263 if delivery.errors: 

1264 return delivery 

1265 

1266 responses.extend(delivery.payload) 

1267 log.extend(delivery.log) 

1268 

1269 if len(self.coefficients_list) > 1: 

1270 log.append(( 

1271 'warning', 

1272 'Multiple filter coefficients lists in single response ' 

1273 'stage.', 

1274 context)) 

1275 

1276 if self.response_list: 

1277 log.append(( 

1278 'warning', 

1279 'Unhandled response element of type: ResponseList', 

1280 context)) 

1281 

1282 if self.polynomial: 

1283 log.append(( 

1284 'warning', 

1285 'Unhandled response element of type: Polynomial', 

1286 context)) 

1287 

1288 if self.stage_gain: 

1289 responses.append( 

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

1291 

1292 return Delivery(responses, log) 

1293 

1294 @property 

1295 def input_units(self): 

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

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

1298 if e is not None: 

1299 return e.input_units 

1300 

1301 return None 

1302 

1303 @property 

1304 def output_units(self): 

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

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

1307 if e is not None: 

1308 return e.output_units 

1309 

1310 return None 

1311 

1312 @property 

1313 def input_sample_rate(self): 

1314 if self.decimation: 

1315 return self.decimation.input_sample_rate.value 

1316 

1317 return None 

1318 

1319 @property 

1320 def output_sample_rate(self): 

1321 if self.decimation: 

1322 return self.decimation.input_sample_rate.value \ 

1323 / self.decimation.factor 

1324 

1325 return None 

1326 

1327 

1328class Response(Object): 

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

1330 instrument_sensitivity = Sensitivity.T(optional=True, 

1331 xmltagname='InstrumentSensitivity') 

1332 instrument_polynomial = Polynomial.T(optional=True, 

1333 xmltagname='InstrumentPolynomial') 

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

1335 

1336 @property 

1337 def output_sample_rate(self): 

1338 if self.stage_list: 

1339 return self.stage_list[-1].output_sample_rate 

1340 else: 

1341 return None 

1342 

1343 def check_sample_rates(self, channel): 

1344 

1345 if self.stage_list: 

1346 sample_rate = None 

1347 

1348 for stage in self.stage_list: 

1349 if stage.decimation: 

1350 input_sample_rate = \ 

1351 stage.decimation.input_sample_rate.value 

1352 

1353 if sample_rate is not None and not same_sample_rate( 

1354 sample_rate, input_sample_rate): 

1355 

1356 logger.warning( 

1357 'Response stage %i has unexpected input sample ' 

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

1359 stage.number, 

1360 input_sample_rate, 

1361 sample_rate)) 

1362 

1363 sample_rate = input_sample_rate / stage.decimation.factor 

1364 

1365 if sample_rate is not None and channel.sample_rate \ 

1366 and not same_sample_rate( 

1367 sample_rate, channel.sample_rate.value): 

1368 

1369 logger.warning( 

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

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

1372 channel.sample_rate.value, 

1373 sample_rate)) 

1374 

1375 def check_units(self): 

1376 

1377 if self.instrument_sensitivity \ 

1378 and self.instrument_sensitivity.input_units: 

1379 

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

1381 

1382 if self.stage_list: 

1383 for stage in self.stage_list: 

1384 if units and stage.input_units \ 

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

1386 

1387 logger.warning( 

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

1389 % ( 

1390 stage.number, 

1391 units, 

1392 'output units of stage %i' 

1393 if stage.number == 0 

1394 else 'sensitivity input units', 

1395 units)) 

1396 

1397 if stage.output_units: 

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

1399 else: 

1400 units = None 

1401 

1402 sout_units = self.instrument_sensitivity.output_units 

1403 if self.instrument_sensitivity and sout_units: 

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

1405 logger.warning( 

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

1407 % ( 

1408 stage.number, 

1409 units, 

1410 'sensitivity output units', 

1411 sout_units.name.upper())) 

1412 

1413 def _sensitivity_checkpoints(self, responses, context): 

1414 delivery = Delivery() 

1415 

1416 if self.instrument_sensitivity: 

1417 sval = self.instrument_sensitivity.value 

1418 sfreq = self.instrument_sensitivity.frequency 

1419 if sval is None: 

1420 delivery.log.append(( 

1421 'warning', 

1422 'No sensitivity value given.', 

1423 context)) 

1424 

1425 elif sval is None: 

1426 delivery.log.append(( 

1427 'warning', 

1428 'Reported sensitivity value is zero.', 

1429 context)) 

1430 

1431 elif sfreq is None: 

1432 delivery.log.append(( 

1433 'warning', 

1434 'Sensitivity frequency not given.', 

1435 context)) 

1436 

1437 else: 

1438 trial = response.MultiplyResponse(responses) 

1439 

1440 delivery.extend( 

1441 check_resp( 

1442 trial, sval, sfreq, 0.1, 

1443 'Instrument sensitivity value inconsistent with ' 

1444 'sensitivity computed from complete response', 

1445 context)) 

1446 

1447 delivery.payload.append(response.FrequencyResponseCheckpoint( 

1448 frequency=sfreq, value=sval)) 

1449 

1450 return delivery 

1451 

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

1453 from pyrocko.squirrel.model import Response 

1454 

1455 if self.stage_list: 

1456 delivery = Delivery() 

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

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

1459 

1460 checkpoints = [] 

1461 if not delivery.errors: 

1462 all_responses = [] 

1463 for sq_stage in delivery.payload: 

1464 all_responses.extend(sq_stage.elements) 

1465 

1466 checkpoints.extend(delivery.extend_without_payload( 

1467 self._sensitivity_checkpoints(all_responses, context))) 

1468 

1469 sq_stages = delivery.payload 

1470 if sq_stages: 

1471 if sq_stages[0].input_quantity is None \ 

1472 and self.instrument_sensitivity is not None: 

1473 

1474 sq_stages[0].input_quantity = to_quantity( 

1475 self.instrument_sensitivity.input_units, 

1476 context, delivery) 

1477 sq_stages[-1].output_quantity = to_quantity( 

1478 self.instrument_sensitivity.output_units, 

1479 context, delivery) 

1480 

1481 sq_stages = delivery.expect() 

1482 

1483 return Response( 

1484 stages=sq_stages, 

1485 log=delivery.log, 

1486 checkpoints=checkpoints, 

1487 **kwargs) 

1488 

1489 elif self.instrument_sensitivity: 

1490 raise NoResponseInformation( 

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

1492 % context) 

1493 else: 

1494 raise NoResponseInformation( 

1495 'Empty instrument response (%s).' 

1496 % context) 

1497 

1498 def get_pyrocko_response( 

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

1500 

1501 delivery = Delivery() 

1502 if self.stage_list: 

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

1504 delivery.extend(stage.get_pyrocko_response( 

1505 context, gain_only=not ( 

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

1507 

1508 elif self.instrument_sensitivity: 

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

1510 

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

1512 checkpoints = delivery.extend_without_payload(delivery_cp) 

1513 if delivery.errors: 

1514 return delivery 

1515 

1516 if fake_input_units is not None: 

1517 if not self.instrument_sensitivity or \ 

1518 self.instrument_sensitivity.input_units is None: 

1519 

1520 delivery.errors.append(( 

1521 'NoResponseInformation', 

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

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

1524 context)) 

1525 

1526 return delivery 

1527 

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

1529 

1530 conresp = None 

1531 try: 

1532 conresp = conversion[ 

1533 fake_input_units.upper(), input_units] 

1534 

1535 except KeyError: 

1536 delivery.errors.append(( 

1537 'NoResponseInformation', 

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

1539 % (fake_input_units, input_units), 

1540 context)) 

1541 

1542 if conresp is not None: 

1543 delivery.payload.append(conresp) 

1544 for checkpoint in checkpoints: 

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

1546 conresp, checkpoint.frequency)) 

1547 

1548 delivery.payload = [ 

1549 response.MultiplyResponse( 

1550 delivery.payload, 

1551 checkpoints=checkpoints)] 

1552 

1553 return delivery 

1554 

1555 @classmethod 

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

1557 normalization_frequency=1.0): 

1558 ''' 

1559 Convert Pyrocko pole-zero response to StationXML response. 

1560 

1561 :param presponse: Pyrocko pole-zero response 

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

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

1564 response. 

1565 :type input_unit: str 

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

1567 response. 

1568 :type output_unit: str 

1569 :param normalization_frequency: Frequency where the normalization 

1570 factor for the StationXML response should be computed. 

1571 :type normalization_frequency: float 

1572 ''' 

1573 

1574 norm_factor = 1.0/float(abs( 

1575 evaluate1(presponse, normalization_frequency) 

1576 / presponse.constant)) 

1577 

1578 pzs = PolesZeros( 

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

1580 normalization_factor=norm_factor, 

1581 normalization_frequency=Frequency(normalization_frequency), 

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

1583 imaginary=FloatNoUnit(z.imag)) 

1584 for z in presponse.zeros], 

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

1586 imaginary=FloatNoUnit(z.imag)) 

1587 for z in presponse.poles]) 

1588 

1589 pzs.validate() 

1590 

1591 stage = ResponseStage( 

1592 number=1, 

1593 poles_zeros_list=[pzs], 

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

1595 

1596 resp = Response( 

1597 instrument_sensitivity=Sensitivity( 

1598 value=stage.stage_gain.value, 

1599 frequency=normalization_frequency, 

1600 input_units=Units(input_unit), 

1601 output_units=Units(output_unit)), 

1602 

1603 stage_list=[stage]) 

1604 

1605 return resp 

1606 

1607 

1608class BaseNode(Object): 

1609 ''' 

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

1611 ''' 

1612 

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

1614 start_date = DummyAwareOptionalTimestamp.T(optional=True, 

1615 xmlstyle='attribute') 

1616 end_date = DummyAwareOptionalTimestamp.T(optional=True, 

1617 xmlstyle='attribute') 

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

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

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

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

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

1623 

1624 def spans(self, *args): 

1625 if len(args) == 0: 

1626 return True 

1627 elif len(args) == 1: 

1628 return ((self.start_date is None or 

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

1630 (self.end_date is None or 

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

1632 

1633 elif len(args) == 2: 

1634 return ((self.start_date is None or 

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

1636 (self.end_date is None or 

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

1638 

1639 

1640def overlaps(a, b): 

1641 return ( 

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

1643 or a.start_date < b.end_date 

1644 ) and ( 

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

1646 or b.start_date < a.end_date 

1647 ) 

1648 

1649 

1650def check_overlaps(node_type_name, codes, nodes): 

1651 errors = [] 

1652 for ia, a in enumerate(nodes): 

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

1654 if overlaps(a, b): 

1655 errors.append( 

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

1657 ' %s - %s\n %s - %s' % ( 

1658 node_type_name, 

1659 '.'.join(codes), 

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

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

1662 

1663 return errors 

1664 

1665 

1666class Channel(BaseNode): 

1667 ''' 

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

1669 response blockettes. 

1670 ''' 

1671 

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

1673 external_reference_list = List.T( 

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

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

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

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

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

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

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

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

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

1683 sample_rate_ratio = SampleRateRatio.T(optional=True, 

1684 xmltagname='SampleRateRatio') 

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

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

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

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

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

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

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

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

1693 

1694 @property 

1695 def position_values(self): 

1696 lat = self.latitude.value 

1697 lon = self.longitude.value 

1698 elevation = value_or_none(self.elevation) 

1699 depth = value_or_none(self.depth) 

1700 return lat, lon, elevation, depth 

1701 

1702 

1703class Station(BaseNode): 

1704 ''' 

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

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

1707 epoch start and end dates. 

1708 ''' 

1709 

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

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

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

1713 site = Site.T(optional=True, xmltagname='Site') 

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

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

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

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

1718 creation_date = DummyAwareOptionalTimestamp.T( 

1719 optional=True, xmltagname='CreationDate') 

1720 termination_date = DummyAwareOptionalTimestamp.T( 

1721 optional=True, xmltagname='TerminationDate') 

1722 total_number_channels = Counter.T( 

1723 optional=True, xmltagname='TotalNumberChannels') 

1724 selected_number_channels = Counter.T( 

1725 optional=True, xmltagname='SelectedNumberChannels') 

1726 external_reference_list = List.T( 

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

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

1729 

1730 @property 

1731 def position_values(self): 

1732 lat = self.latitude.value 

1733 lon = self.longitude.value 

1734 elevation = value_or_none(self.elevation) 

1735 return lat, lon, elevation 

1736 

1737 

1738class Network(BaseNode): 

1739 ''' 

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

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

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

1743 contain 0 or more Stations. 

1744 ''' 

1745 

1746 total_number_stations = Counter.T(optional=True, 

1747 xmltagname='TotalNumberStations') 

1748 selected_number_stations = Counter.T(optional=True, 

1749 xmltagname='SelectedNumberStations') 

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

1751 

1752 @property 

1753 def station_code_list(self): 

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

1755 

1756 @property 

1757 def sl_code_list(self): 

1758 sls = set() 

1759 for station in self.station_list: 

1760 for channel in station.channel_list: 

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

1762 

1763 return sorted(sls) 

1764 

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

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

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

1768 if sls: 

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

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

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

1772 while ssls: 

1773 lines.append( 

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

1775 

1776 ssls[:n] = [] 

1777 

1778 return '\n'.join(lines) 

1779 

1780 

1781def value_or_none(x): 

1782 if x is not None: 

1783 return x.value 

1784 else: 

1785 return None 

1786 

1787 

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

1789 

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

1791 channels[0].position_values 

1792 

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

1794 info = '\n'.join( 

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

1796 x in channels) 

1797 

1798 mess = 'encountered inconsistencies in channel ' \ 

1799 'lat/lon/elevation/depth ' \ 

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

1801 

1802 if inconsistencies == 'raise': 

1803 raise InconsistentChannelLocations(mess) 

1804 

1805 elif inconsistencies == 'warn': 

1806 logger.warning(mess) 

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

1808 

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

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

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

1812 

1813 groups = {} 

1814 for channel in channels: 

1815 if channel.code not in groups: 

1816 groups[channel.code] = [] 

1817 

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

1819 

1820 pchannels = [] 

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

1822 data = [ 

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

1824 value_or_none(channel.dip)) 

1825 for channel in groups[code]] 

1826 

1827 azimuth, dip = util.consistency_merge( 

1828 data, 

1829 message='channel orientation values differ:', 

1830 error=inconsistencies) 

1831 

1832 pchannels.append( 

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

1834 

1835 return pyrocko.model.Station( 

1836 *nsl, 

1837 lat=mlat, 

1838 lon=mlon, 

1839 elevation=mele, 

1840 depth=mdep, 

1841 channels=pchannels) 

1842 

1843 

1844class FDSNStationXML(Object): 

1845 ''' 

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

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

1848 one or more Station containers. 

1849 ''' 

1850 

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

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

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

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

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

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

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

1858 

1859 xmltagname = 'FDSNStationXML' 

1860 guessable_xmlns = [guts_xmlns] 

1861 

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

1863 time=None, timespan=None, 

1864 inconsistencies='warn'): 

1865 

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

1867 

1868 if nslcs is not None: 

1869 nslcs = set(nslcs) 

1870 

1871 if nsls is not None: 

1872 nsls = set(nsls) 

1873 

1874 tt = () 

1875 if time is not None: 

1876 tt = (time,) 

1877 elif timespan is not None: 

1878 tt = timespan 

1879 

1880 pstations = [] 

1881 for network in self.network_list: 

1882 if not network.spans(*tt): 

1883 continue 

1884 

1885 for station in network.station_list: 

1886 if not station.spans(*tt): 

1887 continue 

1888 

1889 if station.channel_list: 

1890 loc_to_channels = {} 

1891 for channel in station.channel_list: 

1892 if not channel.spans(*tt): 

1893 continue 

1894 

1895 loc = channel.location_code.strip() 

1896 if loc not in loc_to_channels: 

1897 loc_to_channels[loc] = [] 

1898 

1899 loc_to_channels[loc].append(channel) 

1900 

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

1902 channels = loc_to_channels[loc] 

1903 if nslcs is not None: 

1904 channels = [channel for channel in channels 

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

1906 channel.code) in nslcs] 

1907 

1908 if not channels: 

1909 continue 

1910 

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

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

1913 continue 

1914 

1915 pstations.append( 

1916 pyrocko_station_from_channels( 

1917 nsl, 

1918 channels, 

1919 inconsistencies=inconsistencies)) 

1920 else: 

1921 pstations.append(pyrocko.model.Station( 

1922 network.code, station.code, '*', 

1923 lat=station.latitude.value, 

1924 lon=station.longitude.value, 

1925 elevation=value_or_none(station.elevation), 

1926 name=station.description or '')) 

1927 

1928 return pstations 

1929 

1930 @classmethod 

1931 def from_pyrocko_stations( 

1932 cls, pyrocko_stations, add_flat_responses_from=None): 

1933 

1934 ''' 

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

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

1937 

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

1939 instances. 

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

1941 ''' 

1942 from collections import defaultdict 

1943 network_dict = defaultdict(list) 

1944 

1945 if add_flat_responses_from: 

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

1947 extra = dict( 

1948 response=Response( 

1949 instrument_sensitivity=Sensitivity( 

1950 value=1.0, 

1951 frequency=1.0, 

1952 input_units=Units(name=add_flat_responses_from)))) 

1953 else: 

1954 extra = {} 

1955 

1956 have_offsets = set() 

1957 for s in pyrocko_stations: 

1958 

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

1960 have_offsets.add(s.nsl()) 

1961 

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

1963 channel_list = [] 

1964 for c in s.channels: 

1965 channel_list.append( 

1966 Channel( 

1967 location_code=location, 

1968 code=c.name, 

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

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

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

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

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

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

1975 **extra 

1976 ) 

1977 ) 

1978 

1979 network_dict[network].append( 

1980 Station( 

1981 code=station, 

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

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

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

1985 channel_list=channel_list) 

1986 ) 

1987 

1988 if have_offsets: 

1989 logger.warning( 

1990 'StationXML does not support Cartesian offsets in ' 

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

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

1993 

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

1995 network_list = [] 

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

1997 

1998 network_list.append( 

1999 Network( 

2000 code=k, station_list=station_list, 

2001 total_number_stations=len(station_list))) 

2002 

2003 sxml = FDSNStationXML( 

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

2005 network_list=network_list) 

2006 

2007 sxml.validate() 

2008 return sxml 

2009 

2010 def iter_network_stations( 

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

2012 

2013 tt = () 

2014 if time is not None: 

2015 tt = (time,) 

2016 elif timespan is not None: 

2017 tt = timespan 

2018 

2019 for network in self.network_list: 

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

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

2022 continue 

2023 

2024 for station in network.station_list: 

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

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

2027 continue 

2028 

2029 yield (network, station) 

2030 

2031 def iter_network_station_channels( 

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

2033 time=None, timespan=None): 

2034 

2035 if loc is not None: 

2036 loc = loc.strip() 

2037 

2038 tt = () 

2039 if time is not None: 

2040 tt = (time,) 

2041 elif timespan is not None: 

2042 tt = timespan 

2043 

2044 for network in self.network_list: 

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

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

2047 continue 

2048 

2049 for station in network.station_list: 

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

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

2052 continue 

2053 

2054 if station.channel_list: 

2055 for channel in station.channel_list: 

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

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

2058 (loc is not None and 

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

2060 continue 

2061 

2062 yield (network, station, channel) 

2063 

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

2065 time=None, timespan=None): 

2066 

2067 groups = {} 

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

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

2070 

2071 net = network.code 

2072 sta = station.code 

2073 cha = channel.code 

2074 loc = channel.location_code.strip() 

2075 if len(cha) == 3: 

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

2077 elif len(cha) == 1: 

2078 bic = '' 

2079 else: 

2080 bic = cha 

2081 

2082 if channel.response and \ 

2083 channel.response.instrument_sensitivity and \ 

2084 channel.response.instrument_sensitivity.input_units: 

2085 

2086 unit = channel.response.instrument_sensitivity\ 

2087 .input_units.name.upper() 

2088 else: 

2089 unit = None 

2090 

2091 bic = (bic, unit) 

2092 

2093 k = net, sta, loc 

2094 if k not in groups: 

2095 groups[k] = {} 

2096 

2097 if bic not in groups[k]: 

2098 groups[k][bic] = [] 

2099 

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

2101 

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

2103 bad_bics = [] 

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

2105 sample_rates = [] 

2106 for channel in channels: 

2107 sample_rates.append(channel.sample_rate.value) 

2108 

2109 if not same(sample_rates): 

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

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

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

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

2114 

2115 logger.warning(err) 

2116 bad_bics.append(bic) 

2117 

2118 for bic in bad_bics: 

2119 del bic_to_channels[bic] 

2120 

2121 return groups 

2122 

2123 def choose_channels( 

2124 self, 

2125 target_sample_rate=None, 

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

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

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

2129 time=None, 

2130 timespan=None): 

2131 

2132 nslcs = {} 

2133 for nsl, bic_to_channels in self.get_channel_groups( 

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

2135 

2136 useful_bics = [] 

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

2138 rate = channels[0].sample_rate.value 

2139 

2140 if target_sample_rate is not None and \ 

2141 rate < target_sample_rate*0.99999: 

2142 continue 

2143 

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

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

2146 continue 

2147 

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

2149 continue 

2150 

2151 unit = bic[1] 

2152 

2153 prio_unit = len(priority_units) 

2154 try: 

2155 prio_unit = priority_units.index(unit) 

2156 except ValueError: 

2157 pass 

2158 

2159 prio_inst = len(priority_instrument_code) 

2160 prio_band = len(priority_band_code) 

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

2162 try: 

2163 prio_inst = priority_instrument_code.index( 

2164 channels[0].code[1]) 

2165 except ValueError: 

2166 pass 

2167 

2168 try: 

2169 prio_band = priority_band_code.index( 

2170 channels[0].code[0]) 

2171 except ValueError: 

2172 pass 

2173 

2174 if target_sample_rate is None: 

2175 rate = -rate 

2176 

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

2178 prio_inst, bic)) 

2179 

2180 useful_bics.sort() 

2181 

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

2183 channels = sorted( 

2184 bic_to_channels[bic], 

2185 key=lambda channel: channel.code) 

2186 

2187 if channels: 

2188 for channel in channels: 

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

2190 

2191 break 

2192 

2193 return nslcs 

2194 

2195 def get_pyrocko_response( 

2196 self, nslc, 

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

2198 

2199 net, sta, loc, cha = nslc 

2200 resps = [] 

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

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

2203 resp = channel.response 

2204 if resp: 

2205 resp.check_sample_rates(channel) 

2206 resp.check_units() 

2207 resps.append(resp.get_pyrocko_response( 

2208 '.'.join(nslc), 

2209 fake_input_units=fake_input_units, 

2210 stages=stages).expect_one()) 

2211 

2212 if not resps: 

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

2214 elif len(resps) > 1: 

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

2216 

2217 return resps[0] 

2218 

2219 @property 

2220 def n_code_list(self): 

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

2222 

2223 @property 

2224 def ns_code_list(self): 

2225 nss = set() 

2226 for network in self.network_list: 

2227 for station in network.station_list: 

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

2229 

2230 return sorted(nss) 

2231 

2232 @property 

2233 def nsl_code_list(self): 

2234 nsls = set() 

2235 for network in self.network_list: 

2236 for station in network.station_list: 

2237 for channel in station.channel_list: 

2238 nsls.add( 

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

2240 

2241 return sorted(nsls) 

2242 

2243 @property 

2244 def nslc_code_list(self): 

2245 nslcs = set() 

2246 for network in self.network_list: 

2247 for station in network.station_list: 

2248 for channel in station.channel_list: 

2249 nslcs.add( 

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

2251 channel.code)) 

2252 

2253 return sorted(nslcs) 

2254 

2255 def summary(self): 

2256 lst = [ 

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

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

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

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

2261 ] 

2262 return '\n'.join(lst) 

2263 

2264 def summary_stages(self): 

2265 data = [] 

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

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

2268 channel.code) 

2269 

2270 stages = [] 

2271 in_units = '?' 

2272 out_units = '?' 

2273 if channel.response: 

2274 sens = channel.response.instrument_sensitivity 

2275 if sens: 

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

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

2278 

2279 for stage in channel.response.stage_list: 

2280 stages.append(stage.summary()) 

2281 

2282 data.append( 

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

2284 in_units, out_units, stages)) 

2285 

2286 data.sort() 

2287 

2288 lst = [] 

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

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

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

2292 for stage in stages: 

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

2294 

2295 return '\n'.join(lst) 

2296 

2297 def _check_overlaps(self): 

2298 by_nslc = {} 

2299 for network in self.network_list: 

2300 for station in network.station_list: 

2301 for channel in station.channel_list: 

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

2303 channel.code) 

2304 if nslc not in by_nslc: 

2305 by_nslc[nslc] = [] 

2306 

2307 by_nslc[nslc].append(channel) 

2308 

2309 errors = [] 

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

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

2312 

2313 return errors 

2314 

2315 def check(self): 

2316 errors = [] 

2317 for meth in [self._check_overlaps]: 

2318 errors.extend(meth()) 

2319 

2320 if errors: 

2321 raise Inconsistencies( 

2322 'Inconsistencies found in StationXML:\n ' 

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

2324 

2325 

2326def load_channel_table(stream): 

2327 

2328 networks = {} 

2329 stations = {} 

2330 

2331 for line in stream: 

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

2333 if line.startswith('#'): 

2334 continue 

2335 

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

2337 

2338 if len(t) != 17: 

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

2340 continue 

2341 

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

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

2344 

2345 try: 

2346 scale = float(scale) 

2347 except ValueError: 

2348 scale = None 

2349 

2350 try: 

2351 scale_freq = float(scale_freq) 

2352 except ValueError: 

2353 scale_freq = None 

2354 

2355 try: 

2356 depth = float(dep) 

2357 except ValueError: 

2358 depth = 0.0 

2359 

2360 try: 

2361 azi = float(azi) 

2362 dip = float(dip) 

2363 except ValueError: 

2364 azi = None 

2365 dip = None 

2366 

2367 try: 

2368 if net not in networks: 

2369 network = Network(code=net) 

2370 else: 

2371 network = networks[net] 

2372 

2373 if (net, sta) not in stations: 

2374 station = Station( 

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

2376 

2377 station.regularize() 

2378 else: 

2379 station = stations[net, sta] 

2380 

2381 if scale: 

2382 resp = Response( 

2383 instrument_sensitivity=Sensitivity( 

2384 value=scale, 

2385 frequency=scale_freq, 

2386 input_units=scale_units)) 

2387 else: 

2388 resp = None 

2389 

2390 channel = Channel( 

2391 code=cha, 

2392 location_code=loc.strip(), 

2393 latitude=lat, 

2394 longitude=lon, 

2395 elevation=ele, 

2396 depth=depth, 

2397 azimuth=azi, 

2398 dip=dip, 

2399 sensor=Equipment(description=sens), 

2400 response=resp, 

2401 sample_rate=sample_rate, 

2402 start_date=start_date, 

2403 end_date=end_date or None) 

2404 

2405 channel.regularize() 

2406 

2407 except ValidationError: 

2408 raise InvalidRecord(line) 

2409 

2410 if net not in networks: 

2411 networks[net] = network 

2412 

2413 if (net, sta) not in stations: 

2414 stations[net, sta] = station 

2415 network.station_list.append(station) 

2416 

2417 station.channel_list.append(channel) 

2418 

2419 return FDSNStationXML( 

2420 source='created from table input', 

2421 created=time.time(), 

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

2423 

2424 

2425def primitive_merge(sxs): 

2426 networks = [] 

2427 for sx in sxs: 

2428 networks.extend(sx.network_list) 

2429 

2430 return FDSNStationXML( 

2431 source='merged from different sources', 

2432 created=time.time(), 

2433 network_list=copy.deepcopy( 

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

2435 

2436 

2437def split_channels(sx): 

2438 for nslc in sx.nslc_code_list: 

2439 network_list = sx.network_list 

2440 network_list_filtered = [ 

2441 network for network in network_list 

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

2443 

2444 for network in network_list_filtered: 

2445 sx.network_list = [network] 

2446 station_list = network.station_list 

2447 station_list_filtered = [ 

2448 station for station in station_list 

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

2450 

2451 for station in station_list_filtered: 

2452 network.station_list = [station] 

2453 channel_list = station.channel_list 

2454 station.channel_list = [ 

2455 channel for channel in channel_list 

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

2457 

2458 yield nslc, copy.deepcopy(sx) 

2459 

2460 station.channel_list = channel_list 

2461 

2462 network.station_list = station_list 

2463 

2464 sx.network_list = network_list 

2465 

2466 

2467if __name__ == '__main__': 

2468 from optparse import OptionParser 

2469 

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

2471 

2472 usage = \ 

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

2474 '<filename> [options]' 

2475 

2476 description = '''Torture StationXML file.''' 

2477 

2478 parser = OptionParser( 

2479 usage=usage, 

2480 description=description, 

2481 formatter=util.BetterHelpFormatter()) 

2482 

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

2484 

2485 if len(args) != 2: 

2486 parser.print_help() 

2487 sys.exit(1) 

2488 

2489 action, path = args 

2490 

2491 sx = load_xml(filename=path) 

2492 if action == 'check': 

2493 try: 

2494 sx.check() 

2495 except Inconsistencies as e: 

2496 logger.error(e) 

2497 sys.exit(1) 

2498 

2499 elif action == 'stats': 

2500 print(sx.summary()) 

2501 

2502 elif action == 'stages': 

2503 print(sx.summary_stages()) 

2504 

2505 else: 

2506 parser.print_help() 

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