1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7QuakeML 1.2 input, output, and data model. 

8 

9This modules provides support to read and write `QuakeML version 1.2 

10<https://quake.ethz.ch/quakeml>`_. It defines a hierarchy of Python objects, 

11closely following the QuakeML data model. 

12 

13QuakeML is a flexible, extensible and modular XML representation of 

14seismological data which is intended to cover a broad range of fields of 

15application in modern seismology. It covers basic seismic event description, 

16including moment tensors. 

17 

18For convenience and ease of use, this documentation contains excerpts from the 

19`QuakeML Manual 

20<https://quake.ethz.ch/quakeml/docs/REC?action=AttachFile&do=get&target=QuakeML-BED-20130214b.pdf>`_. 

21However, this information is not intended to be complete. Please refer to the 

22QuakeML Manual for details. 

23''' 

24 

25import logging 

26from pyrocko.guts import StringPattern, StringChoice, String, Float, Int,\ 

27 Timestamp, Object, List, Union, Bool, Unicode 

28from pyrocko.model import event 

29from pyrocko.gui import marker 

30from pyrocko import moment_tensor 

31import numpy as num 

32 

33logger = logging.getLogger('pyrocko.io.quakeml') 

34 

35guts_prefix = 'quakeml' 

36guts_xmlns = 'http://quakeml.org/xmlns/bed/1.2' 

37polarity_choices = {'positive': 1, 'negative': -1, 'undecidable': None} 

38 

39 

40class QuakeMLError(Exception): 

41 pass 

42 

43 

44class NoPreferredOriginSet(QuakeMLError): 

45 pass 

46 

47 

48def one_element_or_none(li): 

49 if len(li) == 1: 

50 return li[0] 

51 elif len(li) == 0: 

52 return None 

53 else: 

54 logger.warning('More than one element in list: {}'.format(li)) 

55 return None 

56 

57 

58class ResourceIdentifier(StringPattern): 

59 ''' 

60 Identifies resource origin. 

61 

62 They consist of an authority identifier, a unique resource identifier, and 

63 an optional local identifier. The URI schema name smi stands for 

64 seismological meta-information, thus indicating a connection to a set of 

65 metadata associated with the resource. 

66 ''' 

67 pattern = "^(smi|quakeml):[\\w\\d][\\w\\d\\-\\.\\*\\(\\)_~']{2,}/[\\w" +\ 

68 "\\d\\-\\.\\*\\(\\)_~'][\\w\\d\\-\\.\\*\\(\\)\\+\\?_~'=,;#/&]*$" 

69 

70 

71class WhitespaceOrEmptyStringType(StringPattern): 

72 pattern = '^\\s*$' 

73 

74 

75class OriginUncertaintyDescription(StringChoice): 

76 ''' 

77 Preferred uncertainty description. 

78 ''' 

79 choices = [ 

80 'horizontal uncertainty', 

81 'uncertainty ellipse', 

82 'confidence ellipsoid'] 

83 

84 

85class AmplitudeCategory(StringChoice): 

86 ''' 

87 Description of the way the waveform trace is evaluated to get an amplitude 

88 value. 

89 

90 This can be just reading a single value for a given point in time (point), 

91 taking a mean value over a time interval (mean), integrating the trace 

92 over a time interval (integral), specifying just a time interval 

93 (duration), or evaluating a period (period). 

94 ''' 

95 choices = ['point', 'mean', 'duration', 'period', 'integral', 'other'] 

96 

97 

98class OriginDepthType(StringChoice): 

99 ''' 

100 Type of depth determination. 

101 ''' 

102 choices = [ 

103 'from location', 

104 'from moment tensor inversion', 

105 'from modeling of broad-band P waveforms', 

106 'constrained by depth phases', 

107 'constrained by direct phases', 

108 'constrained by depth and direct phases', 

109 'operator assigned', 

110 'other'] 

111 

112 

113class OriginType(StringChoice): 

114 ''' 

115 Describes the origin type. 

116 ''' 

117 choices = [ 

118 'hypocenter', 

119 'centroid', 

120 'amplitude', 

121 'macroseismic', 

122 'rupture start', 

123 'rupture end'] 

124 

125 

126class MTInversionType(StringChoice): 

127 ''' 

128 Type of moment tensor inversion. 

129 

130 Users should avoid to give contradictory information in 

131 :py:class:`MTInversionType` and :py:gattr:`MomentTensor.method_id`. 

132 ''' 

133 choices = ['general', 'zero trace', 'double couple'] 

134 

135 

136class EvaluationMode(StringChoice): 

137 ''' 

138 Mode of an evaluation. 

139 

140 Used in :py:class:`Pick`, :py:class:`Amplitude`, :py:class:`Magnitude`, 

141 :py:class:`Origin`, :py:class:`FocalMechanism`. 

142 ''' 

143 choices = ['manual', 'automatic'] 

144 

145 

146class EvaluationStatus(StringChoice): 

147 ''' 

148 Status of an evaluation. 

149 

150 Used in :py:class:`Pick`, :py:class:`Amplitude`, :py:class:`Magnitude`, 

151 :py:class:`Origin`, :py:class:`FocalMechanism`. 

152 ''' 

153 choices = ['preliminary', 'confirmed', 'reviewed', 'final', 'rejected'] 

154 

155 

156class PickOnset(StringChoice): 

157 ''' 

158 Flag that roughly categorizes the sharpness of the onset. 

159 ''' 

160 choices = ['emergent', 'impulsive', 'questionable'] 

161 

162 

163class EventType(StringChoice): 

164 ''' 

165 Describes the type of an event. 

166 ''' 

167 choices = [ 

168 'not existing', 

169 'not reported', 

170 'earthquake', 

171 'anthropogenic event', 

172 'collapse', 

173 'cavity collapse', 

174 'mine collapse', 

175 'building collapse', 

176 'explosion', 

177 'accidental explosion', 

178 'chemical explosion', 

179 'controlled explosion', 

180 'experimental explosion', 

181 'industrial explosion', 

182 'mining explosion', 

183 'quarry blast', 

184 'road cut', 

185 'blasting levee', 

186 'nuclear explosion', 

187 'induced or triggered event', 

188 'rock burst', 

189 'reservoir loading', 

190 'fluid injection', 

191 'fluid extraction', 

192 'crash', 

193 'plane crash', 

194 'train crash', 

195 'boat crash', 

196 'other event', 

197 'atmospheric event', 

198 'sonic boom', 

199 'sonic blast', 

200 'acoustic noise', 

201 'thunder', 

202 'avalanche', 

203 'snow avalanche', 

204 'debris avalanche', 

205 'hydroacoustic event', 

206 'ice quake', 

207 'slide', 

208 'landslide', 

209 'rockslide', 

210 'meteorite', 

211 'volcanic eruption', 

212 'duplicate earthquake', 

213 'rockburst'] 

214 

215 

216class DataUsedWaveType(StringChoice): 

217 ''' 

218 Type of waveform data. 

219 ''' 

220 choices = [ 

221 'P waves', 

222 'body waves', 

223 'surface waves', 

224 'mantle waves', 

225 'combined', 

226 'unknown'] 

227 

228 

229class AmplitudeUnit(StringChoice): 

230 ''' 

231 Provides the most likely measurement units. 

232 

233 The measurement units for physical quantity are described in the 

234 :py:gattr:`Amplitude.generic_amplitude` attribute. Possible values are 

235 specified as combination of SI base units. 

236 ''' 

237 choices = ['m', 's', 'm/s', 'm/(s*s)', 'm*s', 'dimensionless', 'other'] 

238 

239 

240class EventDescriptionType(StringChoice): 

241 ''' 

242 Category of earthquake description. 

243 ''' 

244 choices = [ 

245 'felt report', 

246 'Flinn-Engdahl region', 

247 'local time', 

248 'tectonic summary', 

249 'nearest cities', 

250 'earthquake name', 

251 'region name'] 

252 

253 

254class MomentTensorCategory(StringChoice): 

255 ''' 

256 Category of moment tensor inversion. 

257 ''' 

258 choices = ['teleseismic', 'regional'] 

259 

260 

261class EventTypeCertainty(StringChoice): 

262 ''' 

263 Denotes how certain the information on event type is. 

264 ''' 

265 choices = ['known', 'suspected'] 

266 

267 

268class SourceTimeFunctionType(StringChoice): 

269 ''' 

270 Type of source time function. 

271 ''' 

272 choices = ['box car', 'triangle', 'trapezoid', 'unknown'] 

273 

274 

275class PickPolarity(StringChoice): 

276 choices = list(polarity_choices.keys()) 

277 

278 

279class AgencyID(String): 

280 pass 

281 

282 

283class Author(Unicode): 

284 pass 

285 

286 

287class Version(String): 

288 pass 

289 

290 

291class Phase(Object): 

292 value = String.T(xmlstyle='content') 

293 

294 

295class GroundTruthLevel(String): 

296 pass 

297 

298 

299class AnonymousNetworkCode(String): 

300 pass 

301 

302 

303class AnonymousStationCode(String): 

304 pass 

305 

306 

307class AnonymousChannelCode(String): 

308 pass 

309 

310 

311class AnonymousLocationCode(String): 

312 pass 

313 

314 

315class Type(String): 

316 pass 

317 

318 

319class MagnitudeHint(String): 

320 pass 

321 

322 

323class Region(Unicode): 

324 pass 

325 

326 

327class RealQuantity(Object): 

328 value = Float.T() 

329 uncertainty = Float.T(optional=True) 

330 lower_uncertainty = Float.T(optional=True) 

331 upper_uncertainty = Float.T(optional=True) 

332 confidence_level = Float.T(optional=True) 

333 

334 

335class IntegerQuantity(Object): 

336 value = Int.T() 

337 uncertainty = Int.T(optional=True) 

338 lower_uncertainty = Int.T(optional=True) 

339 upper_uncertainty = Int.T(optional=True) 

340 confidence_level = Float.T(optional=True) 

341 

342 

343class ConfidenceEllipsoid(Object): 

344 ''' 

345 Representation of the location uncertainty as a confidence ellipsoid with 

346 arbitrary orientation in space. 

347 ''' 

348 semi_major_axis_length = Float.T() 

349 semi_minor_axis_length = Float.T() 

350 semi_intermediate_axis_length = Float.T() 

351 major_axis_plunge = Float.T() 

352 major_axis_azimuth = Float.T() 

353 major_axis_rotation = Float.T() 

354 

355 

356class TimeQuantity(Object): 

357 ''' 

358 Representation of a point in time. 

359 

360 It's given in ISO 8601 format, with optional symmetric or asymmetric 

361 uncertainties given in seconds. The time has to be specified in UTC. 

362 ''' 

363 value = Timestamp.T() 

364 uncertainty = Float.T(optional=True) 

365 lower_uncertainty = Float.T(optional=True) 

366 upper_uncertainty = Float.T(optional=True) 

367 confidence_level = Float.T(optional=True) 

368 

369 

370class TimeWindow(Object): 

371 ''' 

372 Representation of a time window for amplitude measurements. 

373 

374 Which is given by a central point in time, and points in time before and 

375 after this central point. Both points before and after may coincide with 

376 the central point. 

377 ''' 

378 begin = Float.T() 

379 end = Float.T() 

380 reference = Timestamp.T() 

381 

382 

383class ResourceReference(ResourceIdentifier): 

384 ''' 

385 This type is used to refer to QuakeML resources as described in Sect. 3.1 

386 in the `QuakeML manual <https://quake.ethz.ch/quakeml/docs/REC?action=Att\ 

387 achFile&do=get&target=QuakeML-BED-20130214b.pdf>`_. 

388 ''' 

389 pass 

390 

391 

392class DataUsed(Object): 

393 ''' 

394 Description of the type of data used in a moment-tensor inversion. 

395 ''' 

396 wave_type = DataUsedWaveType.T() 

397 station_count = Int.T(optional=True) 

398 component_count = Int.T(optional=True) 

399 shortest_period = Float.T(optional=True) 

400 longest_period = Float.T(optional=True) 

401 

402 

403class EventDescription(Object): 

404 ''' 

405 Free-form string with additional event description. 

406 

407 This can be a well-known name, like 1906 San Francisco Earthquake. 

408 A number of categories can be given in :py:gattr:`type`. 

409 ''' 

410 text = Unicode.T() 

411 type = EventDescriptionType.T(optional=True) 

412 

413 

414class SourceTimeFunction(Object): 

415 ''' 

416 Source time function used in moment-tensor inversion. 

417 ''' 

418 type = SourceTimeFunctionType.T() 

419 duration = Float.T() 

420 rise_time = Float.T(optional=True) 

421 decay_time = Float.T(optional=True) 

422 

423 

424class OriginQuality(Object): 

425 ''' 

426 Description of an origin's quality. 

427 

428 It contains various attributes commonly used to describe the 

429 quality of an origin, e. g., errors, azimuthal coverage, etc. 

430 :py:class:`Origin` objects have an optional attribute of the type 

431 :py:gattr:`OriginQuality`. 

432 ''' 

433 associated_phase_count = Int.T(optional=True) 

434 used_phase_count = Int.T(optional=True) 

435 associated_station_count = Int.T(optional=True) 

436 used_station_count = Int.T(optional=True) 

437 depth_phase_count = Int.T(optional=True) 

438 standard_error = Float.T(optional=True) 

439 azimuthal_gap = Float.T(optional=True) 

440 secondary_azimuthal_gap = Float.T(optional=True) 

441 ground_truth_level = GroundTruthLevel.T(optional=True) 

442 maximum_distance = Float.T(optional=True) 

443 minimum_distance = Float.T(optional=True) 

444 median_distance = Float.T(optional=True) 

445 

446 

447class Axis(Object): 

448 ''' 

449 Representation of an eigenvector of a moment tensor. 

450 

451 Which is expressed in its principal-axes system and uses the angles 

452 :py:gattr:`azimuth`, :py:gattr:`plunge`, and the eigenvalue 

453 :py:gattr:`length`. 

454 ''' 

455 azimuth = RealQuantity.T() 

456 plunge = RealQuantity.T() 

457 length = RealQuantity.T() 

458 

459 

460class Tensor(Object): 

461 ''' 

462 Representation of the six independent moment-tensor elements in spherical 

463 coordinates. 

464 

465 These are the moment-tensor elements Mrr, Mtt, Mpp, Mrt, Mrp, Mtp in the 

466 spherical coordinate system defined by local upward vertical (r), 

467 North-South (t), and West-East (p) directions. 

468 ''' 

469 mrr = RealQuantity.T(xmltagname='Mrr') 

470 mtt = RealQuantity.T(xmltagname='Mtt') 

471 mpp = RealQuantity.T(xmltagname='Mpp') 

472 mrt = RealQuantity.T(xmltagname='Mrt') 

473 mrp = RealQuantity.T(xmltagname='Mrp') 

474 mtp = RealQuantity.T(xmltagname='Mtp') 

475 

476 

477class NodalPlane(Object): 

478 ''' 

479 Description of a nodal plane of a focal mechanism. 

480 ''' 

481 strike = RealQuantity.T() 

482 dip = RealQuantity.T() 

483 rake = RealQuantity.T() 

484 

485 

486class CompositeTime(Object): 

487 ''' 

488 Representation of a time instant. 

489 

490 If the specification is given with no greater accuracy than days (i.e., no 

491 time components are given), the date refers to local time. However, if 

492 time components are given, they have to refer to UTC. 

493 ''' 

494 year = IntegerQuantity.T(optional=True) 

495 month = IntegerQuantity.T(optional=True) 

496 day = IntegerQuantity.T(optional=True) 

497 hour = IntegerQuantity.T(optional=True) 

498 minute = IntegerQuantity.T(optional=True) 

499 second = RealQuantity.T(optional=True) 

500 

501 

502class OriginUncertainty(Object): 

503 ''' 

504 Description of the location uncertainties of an origin. 

505 

506 The uncertainty can be described either as a simple circular horizontal 

507 uncertainty, an uncertainty ellipse according to IMS1.0, or a confidence 

508 ellipsoid. If multiple uncertainty models are given, the preferred variant 

509 can be specified in the attribute :py:gattr:`preferred_description`. 

510 ''' 

511 horizontal_uncertainty = Float.T(optional=True) 

512 min_horizontal_uncertainty = Float.T(optional=True) 

513 max_horizontal_uncertainty = Float.T(optional=True) 

514 azimuth_max_horizontal_uncertainty = Float.T(optional=True) 

515 confidence_ellipsoid = ConfidenceEllipsoid.T(optional=True) 

516 preferred_description = OriginUncertaintyDescription.T(optional=True) 

517 confidence_level = Float.T(optional=True) 

518 

519 

520class ResourceReferenceOptional(Union): 

521 members = [ResourceReference.T(), WhitespaceOrEmptyStringType.T()] 

522 

523 

524class CreationInfo(Object): 

525 ''' 

526 Description of creation metadata (author, version, and creation time) of a 

527 resource. 

528 ''' 

529 agency_id = AgencyID.T(optional=True, xmltagname='agencyID') 

530 agency_uri = ResourceReference.T(optional=True, xmltagname='agencyURI') 

531 author = Author.T(optional=True) 

532 author_uri = ResourceReference.T(optional=True, xmltagname='authorURI') 

533 creation_time = Timestamp.T(optional=True) 

534 version = Version.T(optional=True) 

535 

536 

537class StationMagnitudeContribution(Object): 

538 ''' 

539 Description of the weighting of magnitude values from several 

540 :py:class:`StationMagnitude` objects for computing a network magnitude 

541 estimation. 

542 ''' 

543 station_magnitude_id = ResourceReference.T(xmltagname='stationMagnitudeID') 

544 residual = Float.T(optional=True) 

545 weight = Float.T(optional=True) 

546 

547 

548class PrincipalAxes(Object): 

549 ''' 

550 Representation of the principal axes of a moment tensor. 

551 

552 :py:gattr:`t_axis` and :py:gattr:`p_axis` are required, while 

553 :py:gattr:`n_axis` is optional. 

554 ''' 

555 t_axis = Axis.T() 

556 p_axis = Axis.T() 

557 n_axis = Axis.T(optional=True) 

558 

559 

560class NodalPlanes(Object): 

561 ''' 

562 Representation of the nodal planes of a moment tensor. 

563 

564 The attribute :py:gattr:`preferred_plane` can be used to define which plane 

565 is the preferred one. 

566 ''' 

567 preferred_plane = Int.T(optional=True, xmlstyle='attribute') 

568 nodal_plane1 = NodalPlane.T(optional=True) 

569 nodal_plane2 = NodalPlane.T(optional=True) 

570 

571 

572class WaveformStreamID(Object): 

573 ''' 

574 Reference to a stream description in an inventory. 

575 

576 This is mostly equivalent to the combination of networkCode, stationCode, 

577 locationCode, and channelCode. However, additional information, e.g., 

578 sampling rate, can be referenced by the resourceURI. It is recommended to 

579 use resourceURI as a flexible, abstract, and unique stream ID that allows 

580 to describe different processing levels, or resampled/filtered products of 

581 the same initial stream, without violating the intrinsic meaning of the 

582 legacy identifiers (:py:gattr:`network_code`, :py:gattr:`station_code`, 

583 :py:gattr:`channel_code`, and :py:gattr:`location_code`). However, for 

584 operation in the context of legacy systems, the classical identifier 

585 components are supported. 

586 ''' 

587 value = ResourceReferenceOptional.T(xmlstyle='content') 

588 network_code = AnonymousNetworkCode.T(xmlstyle='attribute') 

589 station_code = AnonymousStationCode.T(xmlstyle='attribute') 

590 channel_code = AnonymousChannelCode.T(optional=True, xmlstyle='attribute') 

591 location_code = AnonymousLocationCode.T( 

592 optional=True, xmlstyle='attribute') 

593 

594 @property 

595 def nslc_id(self): 

596 return (self.network_code, self.station_code, self.location_code, 

597 self.channel_code) 

598 

599 

600class Comment(Object): 

601 ''' 

602 Comment to a resource together with author and creation time information. 

603 ''' 

604 id = ResourceReference.T(optional=True, xmlstyle='attribute') 

605 text = Unicode.T() 

606 creation_info = CreationInfo.T(optional=True) 

607 

608 

609class MomentTensor(Object): 

610 ''' 

611 Representation of a moment tensor solution for an event. 

612 

613 It is an optional part of a :py:class:`FocalMechanism` description. 

614 ''' 

615 public_id = ResourceReference.T( 

616 xmlstyle='attribute', xmltagname='publicID') 

617 data_used_list = List.T(DataUsed.T()) 

618 comment_list = List.T(Comment.T()) 

619 derived_origin_id = ResourceReference.T( 

620 optional=True, xmltagname='derivedOriginID') 

621 moment_magnitude_id = ResourceReference.T( 

622 optional=True, xmltagname='momentMagnitudeID') 

623 scalar_moment = RealQuantity.T(optional=True) 

624 tensor = Tensor.T(optional=True) 

625 variance = Float.T(optional=True) 

626 variance_reduction = Float.T(optional=True) 

627 double_couple = Float.T(optional=True) 

628 clvd = Float.T(optional=True) 

629 iso = Float.T(optional=True) 

630 greens_function_id = ResourceReference.T( 

631 optional=True, xmltagname='greensFunctionID') 

632 filter_id = ResourceReference.T(optional=True, xmltagname='filterID') 

633 source_time_function = SourceTimeFunction.T(optional=True) 

634 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

635 category = MomentTensorCategory.T(optional=True) 

636 inversion_type = MTInversionType.T(optional=True) 

637 creation_info = CreationInfo.T(optional=True) 

638 

639 def pyrocko_moment_tensor(self): 

640 mrr = self.tensor.mrr.value 

641 mtt = self.tensor.mtt.value 

642 mpp = self.tensor.mpp.value 

643 mrt = self.tensor.mrt.value 

644 mrp = self.tensor.mrp.value 

645 mtp = self.tensor.mtp.value 

646 mt = moment_tensor.MomentTensor(m_up_south_east=num.array([ 

647 [mrr, mrt, mrp], [mrt, mtt, mtp], [mrp, mtp, mpp]])) 

648 

649 return mt 

650 

651 

652class Amplitude(Object): 

653 ''' 

654 Quantification of the waveform anomaly. 

655 

656 Usually it consists of a single amplitude measurement or a measurement of 

657 the visible signal duration for duration magnitudes. 

658 ''' 

659 public_id = ResourceReference.T( 

660 xmlstyle='attribute', xmltagname='publicID') 

661 comment_list = List.T(Comment.T()) 

662 generic_amplitude = RealQuantity.T() 

663 type = Type.T(optional=True) 

664 category = AmplitudeCategory.T(optional=True) 

665 unit = AmplitudeUnit.T(optional=True) 

666 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

667 period = RealQuantity.T(optional=True) 

668 snr = Float.T(optional=True) 

669 time_window = TimeWindow.T(optional=True) 

670 pick_id = ResourceReference.T(optional=True, xmltagname='pickID') 

671 waveform_id = WaveformStreamID.T(optional=True, xmltagname='waveformID') 

672 filter_id = ResourceReference.T(optional=True, xmltagname='filterID') 

673 scaling_time = TimeQuantity.T(optional=True) 

674 magnitude_hint = MagnitudeHint.T(optional=True) 

675 evaluation_mode = EvaluationMode.T(optional=True) 

676 evaluation_status = EvaluationStatus.T(optional=True) 

677 creation_info = CreationInfo.T(optional=True) 

678 

679 

680class Magnitude(Object): 

681 ''' 

682 Description of a magnitude value. 

683 

684 It can, but does not need to be associated with an origin. Association 

685 with an origin is expressed with the optional attribute 

686 :py:gattr:`origin_id`. It is either a combination of different magnitude 

687 estimations, or it represents the reported magnitude for the given event. 

688 ''' 

689 public_id = ResourceReference.T( 

690 xmlstyle='attribute', xmltagname='publicID') 

691 comment_list = List.T(Comment.T()) 

692 station_magnitude_contribution_list = List.T( 

693 StationMagnitudeContribution.T()) 

694 mag = RealQuantity.T() 

695 type = Type.T(optional=True) 

696 origin_id = ResourceReference.T(optional=True, xmltagname='originID') 

697 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

698 station_count = Int.T(optional=True) 

699 azimuthal_gap = Float.T(optional=True) 

700 evaluation_mode = EvaluationMode.T(optional=True) 

701 evaluation_status = EvaluationStatus.T(optional=True) 

702 creation_info = CreationInfo.T(optional=True) 

703 

704 

705class StationMagnitude(Object): 

706 ''' 

707 Description of a magnitude derived from a single waveform stream. 

708 ''' 

709 public_id = ResourceReference.T( 

710 xmlstyle='attribute', xmltagname='publicID') 

711 comment_list = List.T(Comment.T()) 

712 origin_id = ResourceReference.T(optional=True, xmltagname='originID') 

713 mag = RealQuantity.T() 

714 type = Type.T(optional=True) 

715 amplitude_id = ResourceReference.T(optional=True, xmltagname='amplitudeID') 

716 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

717 waveform_id = WaveformStreamID.T(optional=True, xmltagname='waveformID') 

718 creation_info = CreationInfo.T(optional=True) 

719 

720 

721class Arrival(Object): 

722 ''' 

723 Successful association of a pick with an origin qualifies this pick as an 

724 arrival. 

725 

726 An arrival thus connects a pick with an origin and provides 

727 additional attributes that describe this relationship. Usually 

728 qualification of a pick as an arrival for a given origin is a hypothesis, 

729 which is based on assumptions about the type of arrival (phase) as well as 

730 observed and (on the basis of an earth model) computed arrival times, or 

731 the residual, respectively. Additional pick attributes like the horizontal 

732 slowness and backazimuth of the observed wave—especially if derived from 

733 array data—may further constrain the nature of the arrival. 

734 ''' 

735 public_id = ResourceReference.T( 

736 xmlstyle='attribute', xmltagname='publicID') 

737 comment_list = List.T(Comment.T()) 

738 pick_id = ResourceReference.T(xmltagname='pickID') 

739 phase = Phase.T() 

740 time_correction = Float.T(optional=True) 

741 azimuth = Float.T(optional=True) 

742 distance = Float.T(optional=True) 

743 takeoff_angle = RealQuantity.T(optional=True) 

744 time_residual = Float.T(optional=True) 

745 horizontal_slowness_residual = Float.T(optional=True) 

746 backazimuth_residual = Float.T(optional=True) 

747 time_weight = Float.T(optional=True) 

748 time_used = Int.T(optional=True) 

749 horizontal_slowness_weight = Float.T(optional=True) 

750 backazimuth_weight = Float.T(optional=True) 

751 earth_model_id = ResourceReference.T( 

752 optional=True, xmltagname='earthModelID') 

753 creation_info = CreationInfo.T(optional=True) 

754 

755 

756class Pick(Object): 

757 ''' 

758 A pick is the observation of an amplitude anomaly in a seismogram at a 

759 specific point in time. 

760 

761 It is not necessarily related to a seismic event. 

762 ''' 

763 public_id = ResourceReference.T( 

764 xmlstyle='attribute', xmltagname='publicID') 

765 comment_list = List.T(Comment.T()) 

766 time = TimeQuantity.T() 

767 waveform_id = WaveformStreamID.T(xmltagname='waveformID') 

768 filter_id = ResourceReference.T(optional=True, xmltagname='filterID') 

769 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

770 horizontal_slowness = RealQuantity.T(optional=True) 

771 backazimuth = RealQuantity.T(optional=True) 

772 slowness_method_id = ResourceReference.T( 

773 optional=True, xmltagname='slownessMethodID') 

774 onset = PickOnset.T(optional=True) 

775 phase_hint = Phase.T(optional=True) 

776 polarity = PickPolarity.T(optional=True) 

777 evaluation_mode = EvaluationMode.T(optional=True) 

778 evaluation_status = EvaluationStatus.T(optional=True) 

779 creation_info = CreationInfo.T(optional=True) 

780 

781 @property 

782 def pyrocko_polarity(self): 

783 return polarity_choices.get(self.polarity, None) 

784 

785 def get_pyrocko_phase_marker(self, event=None): 

786 if not self.phase_hint: 

787 logger.warning('Pick %s: phase_hint undefined' % self.public_id) 

788 phasename = 'undefined' 

789 else: 

790 phasename = self.phase_hint.value 

791 

792 return marker.PhaseMarker( 

793 event=event, nslc_ids=[self.waveform_id.nslc_id], 

794 tmin=self.time.value, tmax=self.time.value, 

795 phasename=phasename, 

796 polarity=self.pyrocko_polarity, 

797 automatic=self.evaluation_mode) 

798 

799 

800class FocalMechanism(Object): 

801 ''' 

802 Description of the focal mechanism of an event. 

803 

804 It includes different descriptions like nodal planes, principal axes, and 

805 a moment tensor. The moment tensor description is provided by objects of 

806 the class :py:class:`MomentTensor` which can be specified as 

807 child elements of :py:class:`FocalMechanism`. 

808 ''' 

809 public_id = ResourceReference.T( 

810 xmlstyle='attribute', xmltagname='publicID') 

811 waveform_id_list = List.T(WaveformStreamID.T(xmltagname='waveformID')) 

812 comment_list = List.T(Comment.T()) 

813 moment_tensor_list = List.T(MomentTensor.T()) 

814 triggering_origin_id = ResourceReference.T( 

815 optional=True, xmltagname='triggeringOriginID') 

816 nodal_planes = NodalPlanes.T(optional=True) 

817 principal_axes = PrincipalAxes.T(optional=True) 

818 azimuthal_gap = Float.T(optional=True) 

819 station_polarity_count = Int.T(optional=True) 

820 misfit = Float.T(optional=True) 

821 station_distribution_ratio = Float.T(optional=True) 

822 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

823 evaluation_mode = EvaluationMode.T(optional=True) 

824 evaluation_status = EvaluationStatus.T(optional=True) 

825 creation_info = CreationInfo.T(optional=True) 

826 

827 

828class Origin(Object): 

829 ''' 

830 Representation of the focal time and geographical location of an earthquake 

831 hypocenter, as well as additional meta-information. 

832 

833 :py:class:`Origin` can have objects of type :py:class:`OriginUncertainty` 

834 and :py:class:`Arrival` as child elements. 

835 ''' 

836 public_id = ResourceReference.T( 

837 xmlstyle='attribute', xmltagname='publicID') 

838 composite_time_list = List.T(CompositeTime.T()) 

839 comment_list = List.T(Comment.T()) 

840 origin_uncertainty_list = List.T(OriginUncertainty.T()) 

841 arrival_list = List.T(Arrival.T()) 

842 time = TimeQuantity.T() 

843 longitude = RealQuantity.T() 

844 latitude = RealQuantity.T() 

845 depth = RealQuantity.T(optional=True) 

846 depth_type = OriginDepthType.T(optional=True) 

847 time_fixed = Bool.T(optional=True) 

848 epicenter_fixed = Bool.T(optional=True) 

849 reference_system_id = ResourceReference.T( 

850 optional=True, xmltagname='referenceSystemID') 

851 method_id = ResourceReference.T(optional=True, xmltagname='methodID') 

852 earth_model_id = ResourceReference.T( 

853 optional=True, xmltagname='earthModelID') 

854 quality = OriginQuality.T(optional=True) 

855 type = OriginType.T(optional=True) 

856 region = Region.T(optional=True) 

857 evaluation_mode = EvaluationMode.T(optional=True) 

858 evaluation_status = EvaluationStatus.T(optional=True) 

859 creation_info = CreationInfo.T(optional=True) 

860 

861 def position_values(self): 

862 lat = self.latitude.value 

863 lon = self.longitude.value 

864 if not self.depth: 

865 logger.warning( 

866 'Origin %s: Depth is undefined. Set to depth=0.' % 

867 self.public_id) 

868 depth = 0. 

869 else: 

870 depth = self.depth.value 

871 

872 return lat, lon, depth 

873 

874 def get_pyrocko_event(self): 

875 lat, lon, depth = self.position_values() 

876 otime = self.time.value 

877 if self.creation_info: 

878 cat = self.creation_info.agency_id 

879 else: 

880 cat = None 

881 

882 return event.Event( 

883 name=self.public_id, 

884 lat=lat, 

885 lon=lon, 

886 time=otime, 

887 depth=depth, 

888 catalog=cat, 

889 region=self.region) 

890 

891 

892class Event(Object): 

893 ''' 

894 Representation of a seismic event. 

895 

896 The Event does not necessarily need to be a tectonic earthquake. An event 

897 is usually associated with one or more origins, which contain information 

898 about focal time and geographical location of the event. Multiple origins 

899 can cover automatic and manual locations, a set of location from different 

900 agencies, locations generated with different location programs and earth 

901 models, etc. Furthermore, an event is usually associated with one or more 

902 magnitudes, and with one or more focal mechanism determinations. In 

903 standard QuakeML-BED, :py:class:`Origin`, :py:class:`Magnitude`, 

904 :py:class:`StationMagnitude`, and :py:class:`FocalMechanism` are child 

905 elements of Event. In BED-RT all these elements are on the same hierarchy 

906 level as child elements of :py:class:`EventParameters`. The association of 

907 origins, magnitudes, and focal mechanisms to a particular event is 

908 expressed using references inside :py:class:`Event`. 

909 ''' 

910 public_id = ResourceReference.T( 

911 xmlstyle='attribute', xmltagname='publicID') 

912 description_list = List.T(EventDescription.T()) 

913 comment_list = List.T(Comment.T()) 

914 focal_mechanism_list = List.T(FocalMechanism.T()) 

915 amplitude_list = List.T(Amplitude.T()) 

916 magnitude_list = List.T(Magnitude.T()) 

917 station_magnitude_list = List.T(StationMagnitude.T()) 

918 origin_list = List.T(Origin.T()) 

919 pick_list = List.T(Pick.T()) 

920 preferred_origin_id = ResourceReference.T( 

921 optional=True, xmltagname='preferredOriginID') 

922 preferred_magnitude_id = ResourceReference.T( 

923 optional=True, xmltagname='preferredMagnitudeID') 

924 preferred_focal_mechanism_id = ResourceReference.T( 

925 optional=True, xmltagname='preferredFocalMechanismID') 

926 type = EventType.T( 

927 optional=True) 

928 type_certainty = EventTypeCertainty.T( 

929 optional=True) 

930 creation_info = CreationInfo.T( 

931 optional=True) 

932 region = Region.T( 

933 optional=True) 

934 

935 def describe(self): 

936 return '''%s: 

937 origins: %i %s 

938 magnitudes: %i %s 

939 focal_machanisms: %i %s 

940 picks: %i 

941 amplitudes: %i 

942 station_magnitudes: %i''' % ( 

943 self.public_id, 

944 len(self.origin_list), 

945 '@' if self.preferred_origin_id else '-', 

946 len(self.magnitude_list), 

947 '@' if self.preferred_magnitude_id else '-', 

948 len(self.focal_mechanism_list), 

949 '@' if self.preferred_focal_mechanism_id else '-', 

950 len(self.pick_list), 

951 len(self.amplitude_list), 

952 len(self.station_magnitude_list)) 

953 

954 def get_pyrocko_phase_markers(self): 

955 event = self.get_pyrocko_event() 

956 return [ 

957 p.get_pyrocko_phase_marker(event=event) for p in self.pick_list] 

958 

959 def get_pyrocko_event(self): 

960 ''' 

961 Convert into Pyrocko event object. 

962 

963 Uses *preferred* origin, magnitude, and moment tensor. If no preferred 

964 item is specified, it picks the first from the list and emits a 

965 warning. 

966 ''' 

967 

968 origin = self.preferred_origin 

969 if not origin and self.origin_list: 

970 origin = self.origin_list[0] 

971 if len(self.origin_list) > 1: 

972 logger.warning( 

973 'Event %s: No preferred origin set, ' 

974 'more than one available, using first' % self.public_id) 

975 

976 if not origin: 

977 raise QuakeMLError( 

978 'No origin available for event: %s' % self.public_id) 

979 

980 ev = origin.get_pyrocko_event() 

981 

982 foc_mech = self.preferred_focal_mechanism 

983 if not foc_mech and self.focal_mechanism_list: 

984 foc_mech = self.focal_mechanism_list[0] 

985 if len(self.focal_mechanism_list) > 1: 

986 logger.warning( 

987 'Event %s: No preferred focal mechanism set, ' 

988 'more than one available, using first' % ev.name) 

989 

990 if foc_mech and foc_mech.moment_tensor_list: 

991 ev.moment_tensor = \ 

992 foc_mech.moment_tensor_list[0].pyrocko_moment_tensor() 

993 

994 if len(foc_mech.moment_tensor_list) > 1: 

995 logger.warning( 

996 'more than one moment tensor available, using first') 

997 

998 mag = None 

999 pref_mag = self.preferred_magnitude 

1000 if pref_mag: 

1001 mag = pref_mag 

1002 elif self.magnitude_list: 

1003 mag = self.magnitude_list[0] 

1004 if len(self.magnitude_list) > 1: 

1005 logger.warning( 

1006 'Event %s: No preferred magnitude set, ' 

1007 'more than one available, using first' % ev.name) 

1008 

1009 if mag: 

1010 ev.magnitude = mag.mag.value 

1011 ev.magnitude_type = mag.type 

1012 

1013 ev.region = self.get_effective_region() 

1014 

1015 return ev 

1016 

1017 def get_effective_region(self): 

1018 if self.region: 

1019 return self.region 

1020 

1021 for desc in self.description_list: 

1022 if desc.type in ('Flinn-Engdahl region', 'region name'): 

1023 return desc.text 

1024 

1025 return None 

1026 

1027 @property 

1028 def preferred_origin(self): 

1029 return one_element_or_none( 

1030 [x for x in self.origin_list 

1031 if x.public_id == self.preferred_origin_id]) 

1032 

1033 @property 

1034 def preferred_magnitude(self): 

1035 return one_element_or_none( 

1036 [x for x in self.magnitude_list 

1037 if x.public_id == self.preferred_magnitude_id]) 

1038 

1039 @property 

1040 def preferred_focal_mechanism(self): 

1041 return one_element_or_none( 

1042 [x for x in self.focal_mechanism_list 

1043 if x.public_id == self.preferred_focal_mechanism_id]) 

1044 

1045 

1046class EventParameters(Object): 

1047 ''' 

1048 In the bulletin-type (non real-time) model, this class serves as a 

1049 container for Event objects. 

1050 

1051 In the real-time version, it can hold objects of type :py:class:`Event`, 

1052 :py:class:`Origin`, :py:class:`Magnitude`, :py:class:`StationMagnitude`, 

1053 :py:class:`FocalMechanism`, Reading, :py:class:`Amplitude` and 

1054 :py:class:`Pick` (real-time mode is not covered by this module at the 

1055 moment). 

1056 ''' 

1057 public_id = ResourceReference.T( 

1058 xmlstyle='attribute', xmltagname='publicID') 

1059 comment_list = List.T(Comment.T()) 

1060 event_list = List.T(Event.T(xmltagname='event')) 

1061 description = Unicode.T(optional=True) 

1062 creation_info = CreationInfo.T(optional=True) 

1063 

1064 

1065class QuakeML(Object): 

1066 ''' 

1067 QuakeML data container. 

1068 ''' 

1069 xmltagname = 'quakeml' 

1070 xmlns = 'http://quakeml.org/xmlns/quakeml/1.2' 

1071 guessable_xmlns = [xmlns, guts_xmlns] 

1072 

1073 event_parameters = EventParameters.T(optional=True) 

1074 

1075 def get_events(self): 

1076 return self.event_parameters.event_list 

1077 

1078 def get_pyrocko_events(self): 

1079 ''' 

1080 Get event information in Pyrocko's basic event format. 

1081 

1082 :rtype: 

1083 List of :py:class:`pyrocko.model.event.Event` objects. 

1084 ''' 

1085 events = [] 

1086 for e in self.event_parameters.event_list: 

1087 events.append(e.get_pyrocko_event()) 

1088 

1089 return events 

1090 

1091 def get_pyrocko_phase_markers(self): 

1092 ''' 

1093 Get pick information in Pyrocko's basic marker format. 

1094 

1095 :rtype: 

1096 List of :py:class:`pyrocko.gui.marker.PhaseMarker` objects. 

1097 ''' 

1098 markers = [] 

1099 for e in self.event_parameters.event_list: 

1100 markers.extend(e.get_pyrocko_phase_markers()) 

1101 

1102 return markers 

1103 

1104 @classmethod 

1105 def load_xml(cls, stream=None, filename=None, string=None): 

1106 ''' 

1107 Load QuakeML data from stream, file or string. 

1108 

1109 :param stream: 

1110 Stream open for reading in binary mode. 

1111 :type stream: 

1112 file-like object, optional 

1113 

1114 :param filename: 

1115 Path to file to be opened for reading. 

1116 :type filename: 

1117 str, optional 

1118 

1119 :param string: 

1120 String with QuakeML data to be deserialized. 

1121 :type string: 

1122 str, optional 

1123 

1124 The arguments ``stream``, ``filename``, and ``string`` are mutually 

1125 exclusive. 

1126 

1127 :returns: 

1128 Parsed QuakeML data structure. 

1129 :rtype: 

1130 :py:class:`QuakeML` object 

1131 

1132 ''' 

1133 

1134 return super(QuakeML, cls).load_xml( 

1135 stream=stream, 

1136 filename=filename, 

1137 string=string, 

1138 ns_hints=[ 

1139 'http://quakeml.org/xmlns/quakeml/1.2', 

1140 'http://quakeml.org/xmlns/bed/1.2'], 

1141 ns_ignore=True)