Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/meta.py: 72%

1445 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-01-10 08:51 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Data model for Green's function store meta-data. 

8''' 

9 

10import math 

11import re 

12import fnmatch 

13import logging 

14 

15import numpy as num 

16from scipy.interpolate import interp1d 

17 

18from pyrocko.guts import (Object, SObject, String, StringChoice, 

19 StringPattern, Unicode, Float, Bool, Int, TBase, 

20 List, ValidationError, Timestamp, Tuple, Dict) 

21from pyrocko.guts import dump, load # noqa 

22from pyrocko.guts_array import literal, Array 

23from pyrocko.model import Location, gnss 

24 

25from pyrocko import cake, orthodrome, spit, moment_tensor, trace 

26 

27from .error import StoreError 

28 

29try: 

30 new_str = unicode 

31except NameError: 

32 new_str = str 

33 

34guts_prefix = 'pf' 

35 

36d2r = math.pi / 180. 

37r2d = 1.0 / d2r 

38km = 1000. 

39vicinity_eps = 1e-5 

40 

41logger = logging.getLogger('pyrocko.gf.meta') 

42 

43 

44def fmt_choices(cls): 

45 return 'choices: %s' % ', '.join( 

46 "``'%s'``" % choice for choice in cls.choices) 

47 

48 

49class InterpolationMethod(StringChoice): 

50 choices = ['nearest_neighbor', 'multilinear'] 

51 

52 

53class SeismosizerTrace(Object): 

54 

55 codes = Tuple.T( 

56 4, String.T(), 

57 default=('', 'STA', '', 'Z'), 

58 help='network, station, location and channel codes') 

59 

60 data = Array.T( 

61 shape=(None,), 

62 dtype=num.float32, 

63 serialize_as='base64', 

64 serialize_dtype=num.dtype('<f4'), 

65 help='numpy array with data samples') 

66 

67 deltat = Float.T( 

68 default=1.0, 

69 help='sampling interval [s]') 

70 

71 tmin = Timestamp.T( 

72 default=Timestamp.D('1970-01-01 00:00:00'), 

73 help='time of first sample as a system timestamp [s]') 

74 

75 def pyrocko_trace(self): 

76 c = self.codes 

77 return trace.Trace(c[0], c[1], c[2], c[3], 

78 ydata=self.data, 

79 deltat=self.deltat, 

80 tmin=self.tmin) 

81 

82 @classmethod 

83 def from_pyrocko_trace(cls, tr, **kwargs): 

84 d = dict( 

85 codes=tr.nslc_id, 

86 tmin=tr.tmin, 

87 deltat=tr.deltat, 

88 data=num.asarray(tr.get_ydata(), dtype=num.float32)) 

89 d.update(kwargs) 

90 return cls(**d) 

91 

92 

93class SeismosizerResult(Object): 

94 n_records_stacked = Int.T(optional=True, default=1) 

95 t_stack = Float.T(optional=True, default=0.) 

96 

97 

98class Result(SeismosizerResult): 

99 trace = SeismosizerTrace.T(optional=True) 

100 n_shared_stacking = Int.T(optional=True, default=1) 

101 t_optimize = Float.T(optional=True, default=0.) 

102 

103 

104class StaticResult(SeismosizerResult): 

105 result = Dict.T( 

106 String.T(), 

107 Array.T(shape=(None,), dtype=float, serialize_as='base64')) 

108 

109 

110class GNSSCampaignResult(StaticResult): 

111 campaign = gnss.GNSSCampaign.T( 

112 optional=True) 

113 

114 

115class SatelliteResult(StaticResult): 

116 

117 theta = Array.T( 

118 optional=True, 

119 shape=(None,), dtype=float, serialize_as='base64') 

120 

121 phi = Array.T( 

122 optional=True, 

123 shape=(None,), dtype=float, serialize_as='base64') 

124 

125 

126class KiteSceneResult(SatelliteResult): 

127 

128 shape = Tuple.T() 

129 

130 def get_scene(self, component='displacement.los'): 

131 try: 

132 from kite import Scene 

133 except ImportError: 

134 raise ImportError('Kite not installed') 

135 

136 def reshape(arr): 

137 return arr.reshape(self.shape) 

138 

139 displacement = self.result[component] 

140 

141 displacement, theta, phi = map( 

142 reshape, (displacement, self.theta, self.phi)) 

143 

144 sc = Scene( 

145 displacement=displacement, 

146 theta=theta, phi=phi, 

147 config=self.config) 

148 

149 return sc 

150 

151 

152class ComponentSchemeDescription(Object): 

153 name = String.T() 

154 source_terms = List.T(String.T()) 

155 ncomponents = Int.T() 

156 default_stored_quantity = String.T(optional=True) 

157 provided_components = List.T(String.T()) 

158 

159 

160component_scheme_descriptions = [ 

161 ComponentSchemeDescription( 

162 name='elastic2', 

163 source_terms=['m0'], 

164 ncomponents=2, 

165 default_stored_quantity='displacement', 

166 provided_components=[ 

167 'n', 'e', 'd']), 

168 ComponentSchemeDescription( 

169 name='elastic8', 

170 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'], 

171 ncomponents=8, 

172 default_stored_quantity='displacement', 

173 provided_components=[ 

174 'n', 'e', 'd']), 

175 ComponentSchemeDescription( 

176 name='elastic10', 

177 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'], 

178 ncomponents=10, 

179 default_stored_quantity='displacement', 

180 provided_components=[ 

181 'n', 'e', 'd']), 

182 ComponentSchemeDescription( 

183 name='elastic10_fd', 

184 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'], 

185 ncomponents=30, 

186 default_stored_quantity='displacement', 

187 provided_components=[ 

188 'n', 'e', 'd']), 

189 ComponentSchemeDescription( 

190 name='elastic18', 

191 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'], 

192 ncomponents=18, 

193 default_stored_quantity='displacement', 

194 provided_components=[ 

195 'n', 'e', 'd']), 

196 ComponentSchemeDescription( 

197 name='elastic5', 

198 source_terms=['fn', 'fe', 'fd'], 

199 ncomponents=5, 

200 default_stored_quantity='displacement', 

201 provided_components=[ 

202 'n', 'e', 'd']), 

203 ComponentSchemeDescription( 

204 name='rotational8', 

205 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'], 

206 ncomponents=8, 

207 default_stored_quantity='rotation_displacement', 

208 provided_components=[ 

209 'n', 'e', 'd']), 

210 ComponentSchemeDescription( 

211 name='poroelastic10', 

212 source_terms=['pore_pressure'], 

213 ncomponents=10, 

214 default_stored_quantity=None, 

215 provided_components=[ 

216 'displacement.n', 'displacement.e', 'displacement.d', 

217 'vertical_tilt.n', 'vertical_tilt.e', 

218 'pore_pressure', 

219 'darcy_velocity.n', 'darcy_velocity.e', 'darcy_velocity.d'])] 

220 

221 

222# new names? 

223# 'mt_to_displacement_1d' 

224# 'mt_to_displacement_1d_ff_only' 

225# 'mt_to_gravity_1d' 

226# 'mt_to_stress_1d' 

227# 'explosion_to_displacement_1d' 

228# 'sf_to_displacement_1d' 

229# 'mt_to_displacement_3d' 

230# 'mt_to_gravity_3d' 

231# 'mt_to_stress_3d' 

232# 'pore_pressure_to_displacement_1d' 

233# 'pore_pressure_to_vertical_tilt_1d' 

234# 'pore_pressure_to_pore_pressure_1d' 

235# 'pore_pressure_to_darcy_velocity_1d' 

236 

237 

238component_schemes = [c.name for c in component_scheme_descriptions] 

239component_scheme_to_description = dict( 

240 (c.name, c) for c in component_scheme_descriptions) 

241 

242 

243class ComponentScheme(StringChoice): 

244 ''' 

245 Different Green's Function component schemes are available: 

246 

247 ================= ========================================================= 

248 Name Description 

249 ================= ========================================================= 

250 ``elastic10`` Elastodynamic for 

251 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and 

252 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT 

253 sources only 

254 ``elastic8`` Elastodynamic for far-field only 

255 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and 

256 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, 

257 MT sources only 

258 ``elastic2`` Elastodynamic for 

259 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and 

260 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, purely 

261 isotropic sources only 

262 ``elastic5`` Elastodynamic for 

263 :py:class:`~pyrocko.gf.meta.ConfigTypeA` 

264 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, SF 

265 sources only 

266 ``elastic18`` Elastodynamic for 

267 :py:class:`~pyrocko.gf.meta.ConfigTypeC` stores, MT 

268 sources only 

269 ``poroelastic10`` Poroelastic for :py:class:`~pyrocko.gf.meta.ConfigTypeA` 

270 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores 

271 ``rotational8`` Elastodynamic rotational motions for 

272 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and 

273 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT 

274 sources only 

275 ================= ========================================================= 

276 ''' 

277 

278 choices = component_schemes 

279 

280 

281class Earthmodel1D(Object): 

282 dummy_for = cake.LayeredModel 

283 dummy_for_description = 'pyrocko.cake.LayeredModel' 

284 

285 class __T(TBase): 

286 def regularize_extra(self, val): 

287 if isinstance(val, str): 

288 val = cake.LayeredModel.from_scanlines( 

289 cake.read_nd_model_str(val)) 

290 

291 return val 

292 

293 def to_save(self, val): 

294 return literal(cake.write_nd_model_str(val)) 

295 

296 

297class StringID(StringPattern): 

298 pattern = r'^[A-Za-z][A-Za-z0-9._]{0,64}$' 

299 

300 

301class ScopeType(StringChoice): 

302 choices = [ 

303 'global', 

304 'regional', 

305 'local', 

306 ] 

307 

308 

309class WaveType(StringChoice): 

310 choices = [ 

311 'full waveform', 

312 'bodywave', 

313 'P wave', 

314 'S wave', 

315 'surface wave', 

316 ] 

317 

318 

319class NearfieldTermsType(StringChoice): 

320 choices = [ 

321 'complete', 

322 'incomplete', 

323 'missing', 

324 ] 

325 

326 

327class QuantityType(StringChoice): 

328 choices = [ 

329 'displacement', 

330 'velocity', 

331 'acceleration', 

332 'rotation_displacement', 

333 'rotation_velocity', 

334 'rotation_acceleration', 

335 'pressure', 

336 'tilt', 

337 'pore_pressure', 

338 'darcy_velocity', 

339 'vertical_tilt'] 

340 

341 

342class Reference(Object): 

343 id = StringID.T() 

344 type = String.T() 

345 title = Unicode.T() 

346 journal = Unicode.T(optional=True) 

347 volume = Unicode.T(optional=True) 

348 number = Unicode.T(optional=True) 

349 pages = Unicode.T(optional=True) 

350 year = String.T() 

351 note = Unicode.T(optional=True) 

352 issn = String.T(optional=True) 

353 doi = String.T(optional=True) 

354 url = String.T(optional=True) 

355 eprint = String.T(optional=True) 

356 authors = List.T(Unicode.T()) 

357 publisher = Unicode.T(optional=True) 

358 keywords = Unicode.T(optional=True) 

359 note = Unicode.T(optional=True) 

360 abstract = Unicode.T(optional=True) 

361 

362 @classmethod 

363 def from_bibtex(cls, filename=None, stream=None): 

364 

365 from pybtex.database.input import bibtex 

366 

367 parser = bibtex.Parser() 

368 

369 if filename is not None: 

370 bib_data = parser.parse_file(filename) 

371 elif stream is not None: 

372 bib_data = parser.parse_stream(stream) 

373 

374 references = [] 

375 

376 for id_, entry in bib_data.entries.items(): 

377 d = {} 

378 avail = entry.fields.keys() 

379 for prop in cls.T.properties: 

380 if prop.name == 'authors' or prop.name not in avail: 

381 continue 

382 

383 d[prop.name] = entry.fields[prop.name] 

384 

385 if 'author' in entry.persons: 

386 d['authors'] = [] 

387 for person in entry.persons['author']: 

388 d['authors'].append(new_str(person)) 

389 

390 c = Reference(id=id_, type=entry.type, **d) 

391 references.append(c) 

392 

393 return references 

394 

395 

396_fpat = r'[+-]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?' 

397_spat = StringID.pattern[1:-1] 

398_cake_pat = cake.PhaseDef.allowed_characters_pattern 

399_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic 

400 

401_ppat = r'(' \ 

402 r'cake:' + _cake_pat + \ 

403 r'|iaspei:' + _iaspei_pat + \ 

404 r'|vel_surface:' + _fpat + \ 

405 r'|vel:' + _fpat + \ 

406 r'|stored:' + _spat + \ 

407 r')' 

408 

409 

410timing_regex_old = re.compile( 

411 r'^((first|last)?\((' + _spat + r'(\|' + _spat + r')*)\)|(' + 

412 _spat + r'))?(' + _fpat + ')?$') 

413 

414 

415timing_regex = re.compile( 

416 r'^((first|last)?\{(' + _ppat + r'(\|' + _ppat + r')*)\}|(' + 

417 _ppat + r'))?(' + _fpat + '(S|%)?)?$') 

418 

419 

420class PhaseSelect(StringChoice): 

421 choices = ['', 'first', 'last'] 

422 

423 

424class InvalidTimingSpecification(ValidationError): 

425 pass 

426 

427 

428class TimingAttributeError(Exception): 

429 pass 

430 

431 

432class Timing(SObject): 

433 ''' 

434 Definition of a time instant relative to one or more named phase arrivals. 

435 

436 Instances of this class can be used e.g. in cutting and tapering 

437 operations. They can hold an absolute time or an offset to a named phase 

438 arrival or group of such arrivals. 

439 

440 Timings can be instantiated from a simple string defintion i.e. with 

441 ``Timing(str)`` where ``str`` is something like 

442 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``, 

443 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase 

444 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'`` 

445 is appended, it is interpreted as percent. If the an ``'S'`` is appended 

446 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km]. 

447 

448 Phase definitions can be specified in either of the following ways: 

449 

450 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table 

451 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake 

452 (see :py:class:`pyrocko.cake.PhaseDef`) 

453 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance / 

454 velocity [km/s] 

455 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s] 

456 

457 **Examples:** 

458 

459 * ``'100'`` : absolute time; 100 s 

460 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to 

461 stored travel time table named ``'P'`` 

462 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to 

463 stored travel time table named ``'P'`` 

464 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to 

465 stored travel time table named ``'P'`` 

466 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is 

467 undefined for a given geometry 

468 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A 

469 and B is chosen, if both phases are defined for a given geometry 

470 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen 

471 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of 

472 arrivals A, B, and C 

473 ''' 

474 

475 def __init__(self, s=None, **kwargs): 

476 

477 if s is not None: 

478 offset_is = None 

479 s = re.sub(r'\s+', '', s) 

480 try: 

481 offset = float(s.rstrip('S')) 

482 

483 if s.endswith('S'): 

484 offset_is = 'slowness' 

485 

486 phase_defs = [] 

487 select = '' 

488 

489 except ValueError: 

490 matched = False 

491 m = timing_regex.match(s) 

492 if m: 

493 if m.group(3): 

494 phase_defs = m.group(3).split('|') 

495 elif m.group(19): 

496 phase_defs = [m.group(19)] 

497 else: 

498 phase_defs = [] 

499 

500 select = m.group(2) or '' 

501 

502 offset = 0.0 

503 soff = m.group(27) 

504 if soff: 

505 offset = float(soff.rstrip('S%')) 

506 if soff.endswith('S'): 

507 offset_is = 'slowness' 

508 elif soff.endswith('%'): 

509 offset_is = 'percent' 

510 

511 matched = True 

512 

513 else: 

514 m = timing_regex_old.match(s) 

515 if m: 

516 if m.group(3): 

517 phase_defs = m.group(3).split('|') 

518 elif m.group(5): 

519 phase_defs = [m.group(5)] 

520 else: 

521 phase_defs = [] 

522 

523 select = m.group(2) or '' 

524 

525 offset = 0.0 

526 if m.group(6): 

527 offset = float(m.group(6)) 

528 

529 phase_defs = [ 

530 'stored:' + phase_def for phase_def in phase_defs] 

531 

532 matched = True 

533 

534 if not matched: 

535 raise InvalidTimingSpecification(s) 

536 

537 kwargs = dict( 

538 phase_defs=phase_defs, 

539 select=select, 

540 offset=offset, 

541 offset_is=offset_is) 

542 

543 SObject.__init__(self, **kwargs) 

544 

545 def __str__(self): 

546 s = [] 

547 if self.phase_defs: 

548 sphases = '|'.join(self.phase_defs) 

549 # if len(self.phase_defs) > 1 or self.select: 

550 sphases = '{%s}' % sphases 

551 

552 if self.select: 

553 sphases = self.select + sphases 

554 

555 s.append(sphases) 

556 

557 if self.offset != 0.0 or not self.phase_defs: 

558 s.append('%+g' % self.offset) 

559 if self.offset_is == 'slowness': 

560 s.append('S') 

561 elif self.offset_is == 'percent': 

562 s.append('%') 

563 

564 return ''.join(s) 

565 

566 def evaluate(self, get_phase, args, attributes=None): 

567 try: 

568 if self.offset_is == 'slowness' and self.offset != 0.0: 

569 phase_offset = get_phase( 

570 'vel_surface:%g' % (1.0/self.offset)) 

571 offset = phase_offset(args) 

572 else: 

573 offset = self.offset 

574 

575 if self.phase_defs: 

576 extra_args = () 

577 if attributes: 

578 extra_args = (attributes,) 

579 

580 phases = [ 

581 get_phase(phase_def, *extra_args) 

582 for phase_def in self.phase_defs] 

583 

584 results = [phase(args) for phase in phases] 

585 results = [ 

586 result if isinstance(result, tuple) else (result,) 

587 for result in results] 

588 

589 results = [ 

590 result for result in results 

591 if result[0] is not None] 

592 

593 if self.select == 'first': 

594 results.sort(key=lambda result: result[0]) 

595 elif self.select == 'last': 

596 results.sort(key=lambda result: -result[0]) 

597 

598 if not results: 

599 if attributes is not None: 

600 return (None,) * (len(attributes) + 1) 

601 else: 

602 return None 

603 

604 else: 

605 t, values = results[0][0], results[0][1:] 

606 if self.offset_is == 'percent': 

607 t = t*(1.+offset/100.) 

608 else: 

609 t = t + offset 

610 

611 if attributes is not None: 

612 return (t, ) + values 

613 else: 

614 return t 

615 

616 else: 

617 if attributes is not None: 

618 raise TimingAttributeError( 

619 'Cannot get phase attributes from offset-only ' 

620 'definition.') 

621 return offset 

622 

623 except spit.OutOfBounds: 

624 raise OutOfBounds(args) 

625 

626 phase_defs = List.T(String.T()) 

627 offset = Float.T(default=0.0) 

628 offset_is = String.T(optional=True) 

629 select = PhaseSelect.T( 

630 default='', 

631 help=("Can be either ``'%s'``, ``'%s'``, or ``'%s'``. " % 

632 tuple(PhaseSelect.choices))) 

633 

634 

635def mkdefs(s): 

636 defs = [] 

637 for x in s.split(','): 

638 try: 

639 defs.append(float(x)) 

640 except ValueError: 

641 if x.startswith('!'): 

642 defs.extend(cake.PhaseDef.classic(x[1:])) 

643 else: 

644 defs.append(cake.PhaseDef(x)) 

645 

646 return defs 

647 

648 

649class TPDef(Object): 

650 ''' 

651 Maps an arrival phase identifier to an arrival phase definition. 

652 ''' 

653 

654 id = StringID.T( 

655 help='name used to identify the phase') 

656 definition = String.T( 

657 help='definition of the phase in either cake syntax as defined in ' 

658 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an ' 

659 '``!``, as a *classic phase name*, or, if it is a simple ' 

660 'number, as a constant horizontal velocity.') 

661 

662 @property 

663 def phases(self): 

664 return [x for x in mkdefs(self.definition) 

665 if isinstance(x, cake.PhaseDef)] 

666 

667 @property 

668 def horizontal_velocities(self): 

669 return [x for x in mkdefs(self.definition) if isinstance(x, float)] 

670 

671 

672class OutOfBounds(Exception): 

673 def __init__(self, values=None, reason=None): 

674 Exception.__init__(self) 

675 self.values = values 

676 self.reason = reason 

677 self.context = None 

678 

679 def __str__(self): 

680 scontext = '' 

681 if self.context: 

682 scontext = '\n%s' % str(self.context) 

683 

684 if self.reason: 

685 scontext += '\n%s' % self.reason 

686 

687 if self.values: 

688 return 'out of bounds: (%s)%s' % ( 

689 ', '.join('%g' % x for x in self.values), scontext) 

690 else: 

691 return 'out of bounds%s' % scontext 

692 

693 

694class MultiLocation(Object): 

695 ''' 

696 Unstructured grid of point coordinates. 

697 ''' 

698 

699 lats = Array.T( 

700 optional=True, shape=(None,), dtype=float, 

701 help='Latitudes of targets.') 

702 lons = Array.T( 

703 optional=True, shape=(None,), dtype=float, 

704 help='Longitude of targets.') 

705 north_shifts = Array.T( 

706 optional=True, shape=(None,), dtype=float, 

707 help='North shifts of targets.') 

708 east_shifts = Array.T( 

709 optional=True, shape=(None,), dtype=float, 

710 help='East shifts of targets.') 

711 elevation = Array.T( 

712 optional=True, shape=(None,), dtype=float, 

713 help='Elevations of targets.') 

714 

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

716 self._coords5 = None 

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

718 

719 @property 

720 def coords5(self): 

721 if self._coords5 is not None: 

722 return self._coords5 

723 props = [self.lats, self.lons, self.north_shifts, self.east_shifts, 

724 self.elevation] 

725 sizes = [p.size for p in props if p is not None] 

726 if not sizes: 

727 raise AttributeError('Empty StaticTarget') 

728 if num.unique(sizes).size != 1: 

729 raise AttributeError('Inconsistent coordinate sizes.') 

730 

731 self._coords5 = num.zeros((sizes[0], 5)) 

732 for idx, p in enumerate(props): 

733 if p is not None: 

734 self._coords5[:, idx] = p.flatten() 

735 

736 return self._coords5 

737 

738 @property 

739 def ncoords(self): 

740 if self.coords5.shape[0] is None: 

741 return 0 

742 return int(self.coords5.shape[0]) 

743 

744 def get_latlon(self): 

745 ''' 

746 Get all coordinates as lat/lon. 

747 

748 :returns: Coordinates in Latitude, Longitude 

749 :rtype: :class:`numpy.ndarray`, (N, 2) 

750 ''' 

751 latlons = num.empty((self.ncoords, 2)) 

752 for i in range(self.ncoords): 

753 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4]) 

754 return latlons 

755 

756 def get_corner_coords(self): 

757 ''' 

758 Returns the corner coordinates of the multi-location object. 

759 

760 :returns: In lat/lon: lower left, upper left, upper right, lower right 

761 :rtype: tuple 

762 ''' 

763 latlon = self.get_latlon() 

764 ll = (latlon[:, 0].min(), latlon[:, 1].min()) 

765 ur = (latlon[:, 0].max(), latlon[:, 1].max()) 

766 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1])) 

767 

768 

769class Receiver(Location): 

770 codes = Tuple.T( 

771 3, String.T(), 

772 optional=True, 

773 help='network, station, and location codes') 

774 

775 def pyrocko_station(self): 

776 from pyrocko import model 

777 

778 lat, lon = self.effective_latlon 

779 

780 return model.Station( 

781 network=self.codes[0], 

782 station=self.codes[1], 

783 location=self.codes[2], 

784 lat=lat, 

785 lon=lon, 

786 depth=self.depth) 

787 

788 

789def g(x, d): 

790 if x is None: 

791 return d 

792 else: 

793 return x 

794 

795 

796class UnavailableScheme(Exception): 

797 ''' 

798 Raised when it is not possible to use the requested component scheme. 

799 ''' 

800 pass 

801 

802 

803class InvalidNComponents(Exception): 

804 ''' 

805 Raised when ``ncomponents`` is incompatible with ``component_scheme``. 

806 ''' 

807 pass 

808 

809 

810class DiscretizedSource(Object): 

811 ''' 

812 Base class for discretized sources. 

813 

814 To compute synthetic seismograms, the parameterized source models (see of 

815 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first 

816 discretized into a number of point sources. These spacio-temporal point 

817 source distributions are represented by subclasses of the 

818 :py:class:`DiscretizedSource`. For elastodynamic problems there is the 

819 :py:class:`DiscretizedMTSource` for moment tensor point source 

820 distributions and the :py:class:`DiscretizedExplosionSource` for pure 

821 explosion/implosion type source distributions. The geometry calculations 

822 are implemented in the base class. How Green's function components have to 

823 be weighted and sumed is defined in the derived classes. 

824 

825 Like in the :py:class:`pyrocko.model.location.Location` class, the 

826 positions of the point sources contained in the discretized source are 

827 defined by a common reference point (:py:attr:`lat`, :py:attr:`lon`) and 

828 cartesian offsets to this (:py:attr:`north_shifts`, :py:attr:`east_shifts`, 

829 :py:attr:`depths`). Alternatively latitude and longitude of each contained 

830 point source can be specified directly (:py:attr:`lats`, :py:attr:`lons`). 

831 ''' 

832 

833 times = Array.T(shape=(None,), dtype=float) 

834 lats = Array.T(shape=(None,), dtype=float, optional=True) 

835 lons = Array.T(shape=(None,), dtype=float, optional=True) 

836 lat = Float.T(optional=True) 

837 lon = Float.T(optional=True) 

838 north_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True) 

839 east_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True) 

840 depths = Array.T(shape=(None,), dtype=num.float64) 

841 dl = Float.T(optional=True) 

842 dw = Float.T(optional=True) 

843 nl = Float.T(optional=True) 

844 nw = Float.T(optional=True) 

845 

846 @classmethod 

847 def check_scheme(cls, scheme): 

848 ''' 

849 Check if given GF component scheme is supported by the class. 

850 

851 Raises :py:class:`UnavailableScheme` if the given scheme is not 

852 supported by this discretized source class. 

853 ''' 

854 

855 if scheme not in cls.provided_schemes: 

856 raise UnavailableScheme( 

857 'source type "%s" does not support GF component scheme "%s"' % 

858 (cls.__name__, scheme)) 

859 

860 def __init__(self, **kwargs): 

861 Object.__init__(self, **kwargs) 

862 self._latlons = None 

863 

864 def __setattr__(self, name, value): 

865 if name in ('lat', 'lon', 'north_shifts', 'east_shifts', 

866 'lats', 'lons'): 

867 self.__dict__['_latlons'] = None 

868 

869 Object.__setattr__(self, name, value) 

870 

871 def get_source_terms(self, scheme): 

872 raise NotImplementedError() 

873 

874 def make_weights(self, receiver, scheme): 

875 raise NotImplementedError() 

876 

877 @property 

878 def effective_latlons(self): 

879 ''' 

880 Property holding the offest-corrected lats and lons of all points. 

881 ''' 

882 

883 if self._latlons is None: 

884 if self.lats is not None and self.lons is not None: 

885 if (self.north_shifts is not None and 

886 self.east_shifts is not None): 

887 self._latlons = orthodrome.ne_to_latlon( 

888 self.lats, self.lons, 

889 self.north_shifts, self.east_shifts) 

890 else: 

891 self._latlons = self.lats, self.lons 

892 else: 

893 lat = g(self.lat, 0.0) 

894 lon = g(self.lon, 0.0) 

895 self._latlons = orthodrome.ne_to_latlon( 

896 lat, lon, self.north_shifts, self.east_shifts) 

897 

898 return self._latlons 

899 

900 @property 

901 def effective_north_shifts(self): 

902 if self.north_shifts is not None: 

903 return self.north_shifts 

904 else: 

905 return num.zeros(self.times.size) 

906 

907 @property 

908 def effective_east_shifts(self): 

909 if self.east_shifts is not None: 

910 return self.east_shifts 

911 else: 

912 return num.zeros(self.times.size) 

913 

914 def same_origin(self, receiver): 

915 ''' 

916 Check if receiver has same reference point. 

917 ''' 

918 

919 return (g(self.lat, 0.0) == receiver.lat and 

920 g(self.lon, 0.0) == receiver.lon and 

921 self.lats is None and self.lons is None) 

922 

923 def azibazis_to(self, receiver): 

924 ''' 

925 Compute azimuths and backaziumuths to/from receiver, for all contained 

926 points. 

927 ''' 

928 

929 if self.same_origin(receiver): 

930 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts, 

931 receiver.north_shift - self.north_shifts) 

932 bazis = azis + 180. 

933 else: 

934 slats, slons = self.effective_latlons 

935 rlat, rlon = receiver.effective_latlon 

936 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon) 

937 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons) 

938 

939 return azis, bazis 

940 

941 def distances_to(self, receiver): 

942 ''' 

943 Compute distances to receiver for all contained points. 

944 ''' 

945 if self.same_origin(receiver): 

946 return num.sqrt((self.north_shifts - receiver.north_shift)**2 + 

947 (self.east_shifts - receiver.east_shift)**2) 

948 

949 else: 

950 slats, slons = self.effective_latlons 

951 rlat, rlon = receiver.effective_latlon 

952 return orthodrome.distance_accurate50m_numpy(slats, slons, 

953 rlat, rlon) 

954 

955 def element_coords(self, i): 

956 if self.lats is not None and self.lons is not None: 

957 lat = float(self.lats[i]) 

958 lon = float(self.lons[i]) 

959 else: 

960 lat = self.lat 

961 lon = self.lon 

962 

963 if self.north_shifts is not None and self.east_shifts is not None: 

964 north_shift = float(self.north_shifts[i]) 

965 east_shift = float(self.east_shifts[i]) 

966 

967 else: 

968 north_shift = east_shift = 0.0 

969 

970 return lat, lon, north_shift, east_shift 

971 

972 def coords5(self): 

973 xs = num.zeros((self.nelements, 5)) 

974 

975 if self.lats is not None and self.lons is not None: 

976 xs[:, 0] = self.lats 

977 xs[:, 1] = self.lons 

978 else: 

979 xs[:, 0] = g(self.lat, 0.0) 

980 xs[:, 1] = g(self.lon, 0.0) 

981 

982 if self.north_shifts is not None and self.east_shifts is not None: 

983 xs[:, 2] = self.north_shifts 

984 xs[:, 3] = self.east_shifts 

985 

986 xs[:, 4] = self.depths 

987 

988 return xs 

989 

990 @property 

991 def nelements(self): 

992 return self.times.size 

993 

994 @classmethod 

995 def combine(cls, sources, **kwargs): 

996 ''' 

997 Combine several discretized source models. 

998 

999 Concatenenates all point sources in the given discretized ``sources``. 

1000 Care must be taken when using this function that the external amplitude 

1001 factors and reference times of the parameterized (undiscretized) 

1002 sources match or are accounted for. 

1003 ''' 

1004 

1005 first = sources[0] 

1006 

1007 if not all(type(s) is type(first) for s in sources): 

1008 raise Exception('DiscretizedSource.combine must be called with ' 

1009 'sources of same type.') 

1010 

1011 latlons = [] 

1012 for s in sources: 

1013 latlons.append(s.effective_latlons) 

1014 

1015 lats, lons = num.hstack(latlons) 

1016 

1017 if all((s.lats is None and s.lons is None) for s in sources): 

1018 rlats = num.array([s.lat for s in sources], dtype=float) 

1019 rlons = num.array([s.lon for s in sources], dtype=float) 

1020 same_ref = num.all( 

1021 rlats == rlats[0]) and num.all(rlons == rlons[0]) 

1022 else: 

1023 same_ref = False 

1024 

1025 cat = num.concatenate 

1026 times = cat([s.times for s in sources]) 

1027 print(times) 

1028 depths = cat([s.depths for s in sources]) 

1029 

1030 if same_ref: 

1031 lat = first.lat 

1032 lon = first.lon 

1033 north_shifts = cat([s.effective_north_shifts for s in sources]) 

1034 east_shifts = cat([s.effective_east_shifts for s in sources]) 

1035 lats = None 

1036 lons = None 

1037 else: 

1038 lat = None 

1039 lon = None 

1040 north_shifts = None 

1041 east_shifts = None 

1042 

1043 return cls( 

1044 times=times, lat=lat, lon=lon, lats=lats, lons=lons, 

1045 north_shifts=north_shifts, east_shifts=east_shifts, 

1046 depths=depths, **kwargs) 

1047 

1048 def centroid_position(self): 

1049 moments = self.moments() 

1050 norm = num.sum(moments) 

1051 if norm != 0.0: 

1052 w = moments / num.sum(moments) 

1053 else: 

1054 w = num.ones(moments.size) 

1055 

1056 if self.lats is not None and self.lons is not None: 

1057 lats, lons = self.effective_latlons 

1058 rlat, rlon = num.mean(lats), num.mean(lons) 

1059 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons) 

1060 else: 

1061 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0) 

1062 n, e = self.north_shifts, self.east_shifts 

1063 

1064 cn = num.sum(n*w) 

1065 ce = num.sum(e*w) 

1066 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce) 

1067 

1068 if self.lats is not None and self.lons is not None: 

1069 lat = clat 

1070 lon = clon 

1071 north_shift = 0. 

1072 east_shift = 0. 

1073 else: 

1074 lat = g(self.lat, 0.0) 

1075 lon = g(self.lon, 0.0) 

1076 north_shift = cn 

1077 east_shift = ce 

1078 

1079 depth = num.sum(self.depths*w) 

1080 time = num.sum(self.times*w) 

1081 return tuple(float(x) for x in 

1082 (time, lat, lon, north_shift, east_shift, depth)) 

1083 

1084 

1085class DiscretizedExplosionSource(DiscretizedSource): 

1086 m0s = Array.T(shape=(None,), dtype=float) 

1087 

1088 provided_schemes = ( 

1089 'elastic2', 

1090 'elastic8', 

1091 'elastic10', 

1092 'rotational8', 

1093 ) 

1094 

1095 def get_source_terms(self, scheme): 

1096 self.check_scheme(scheme) 

1097 

1098 if scheme == 'elastic2': 

1099 return self.m0s[:, num.newaxis].copy() 

1100 

1101 elif scheme in ('elastic8', 'elastic10', 'rotational8'): 

1102 m6s = num.zeros((self.m0s.size, 6)) 

1103 m6s[:, 0:3] = self.m0s[:, num.newaxis] 

1104 return m6s 

1105 else: 

1106 assert False 

1107 

1108 def make_weights(self, receiver, scheme): 

1109 self.check_scheme(scheme) 

1110 

1111 azis, bazis = self.azibazis_to(receiver) 

1112 

1113 sb = num.sin(bazis*d2r-num.pi) 

1114 cb = num.cos(bazis*d2r-num.pi) 

1115 

1116 m0s = self.m0s 

1117 n = azis.size 

1118 

1119 cat = num.concatenate 

1120 rep = num.repeat 

1121 

1122 if scheme == 'elastic2': 

1123 w_n = cb*m0s 

1124 g_n = filledi(0, n) 

1125 w_e = sb*m0s 

1126 g_e = filledi(0, n) 

1127 w_d = m0s 

1128 g_d = filledi(1, n) 

1129 

1130 elif scheme == 'elastic8': 

1131 w_n = cat((cb*m0s, cb*m0s)) 

1132 g_n = rep((0, 2), n) 

1133 w_e = cat((sb*m0s, sb*m0s)) 

1134 g_e = rep((0, 2), n) 

1135 w_d = cat((m0s, m0s)) 

1136 g_d = rep((5, 7), n) 

1137 

1138 elif scheme == 'elastic10': 

1139 w_n = cat((cb*m0s, cb*m0s, cb*m0s)) 

1140 g_n = rep((0, 2, 8), n) 

1141 w_e = cat((sb*m0s, sb*m0s, sb*m0s)) 

1142 g_e = rep((0, 2, 8), n) 

1143 w_d = cat((m0s, m0s, m0s)) 

1144 g_d = rep((5, 7, 9), n) 

1145 

1146 elif scheme == 'rotational8': 

1147 w_n = cat((sb*m0s, sb*m0s, sb*m0s)) 

1148 g_n = rep((2, 4, 7), n) 

1149 w_e = cat((sb*m0s, sb*m0s, sb*m0s)) 

1150 g_e = rep((2, 4, 7), n) 

1151 w_d = num.zeros(0) 

1152 g_d = num.zeros(0, dtype=int) 

1153 logger.warning('untested code executed, check for sign errors') 

1154 

1155 else: 

1156 assert False 

1157 

1158 return ( 

1159 ('displacement.n', w_n, g_n), 

1160 ('displacement.e', w_e, g_e), 

1161 ('displacement.d', w_d, g_d), 

1162 ) 

1163 

1164 def split(self): 

1165 from pyrocko.gf.seismosizer import ExplosionSource 

1166 sources = [] 

1167 for i in range(self.nelements): 

1168 lat, lon, north_shift, east_shift = self.element_coords(i) 

1169 sources.append(ExplosionSource( 

1170 time=float(self.times[i]), 

1171 lat=lat, 

1172 lon=lon, 

1173 north_shift=north_shift, 

1174 east_shift=east_shift, 

1175 depth=float(self.depths[i]), 

1176 moment=float(self.m0s[i]))) 

1177 

1178 return sources 

1179 

1180 def moments(self): 

1181 return self.m0s 

1182 

1183 def centroid(self): 

1184 from pyrocko.gf.seismosizer import ExplosionSource 

1185 time, lat, lon, north_shift, east_shift, depth = \ 

1186 self.centroid_position() 

1187 

1188 return ExplosionSource( 

1189 time=time, 

1190 lat=lat, 

1191 lon=lon, 

1192 north_shift=north_shift, 

1193 east_shift=east_shift, 

1194 depth=depth, 

1195 moment=float(num.sum(self.m0s))) 

1196 

1197 @classmethod 

1198 def combine(cls, sources, **kwargs): 

1199 ''' 

1200 Combine several discretized source models. 

1201 

1202 Concatenenates all point sources in the given discretized ``sources``. 

1203 Care must be taken when using this function that the external amplitude 

1204 factors and reference times of the parameterized (undiscretized) 

1205 sources match or are accounted for. 

1206 ''' 

1207 

1208 if 'm0s' not in kwargs: 

1209 kwargs['m0s'] = num.concatenate([s.m0s for s in sources]) 

1210 

1211 return super(DiscretizedExplosionSource, cls).combine(sources, 

1212 **kwargs) 

1213 

1214 

1215class DiscretizedSFSource(DiscretizedSource): 

1216 forces = Array.T(shape=(None, 3), dtype=float) 

1217 

1218 provided_schemes = ( 

1219 'elastic5', 

1220 ) 

1221 

1222 def get_source_terms(self, scheme): 

1223 self.check_scheme(scheme) 

1224 

1225 return self.forces 

1226 

1227 def make_weights(self, receiver, scheme): 

1228 self.check_scheme(scheme) 

1229 

1230 azis, bazis = self.azibazis_to(receiver) 

1231 

1232 sa = num.sin(azis*d2r) 

1233 ca = num.cos(azis*d2r) 

1234 sb = num.sin(bazis*d2r-num.pi) 

1235 cb = num.cos(bazis*d2r-num.pi) 

1236 

1237 forces = self.forces 

1238 fn = forces[:, 0] 

1239 fe = forces[:, 1] 

1240 fd = forces[:, 2] 

1241 

1242 f0 = fd 

1243 f1 = ca * fn + sa * fe 

1244 f2 = ca * fe - sa * fn 

1245 

1246 n = azis.size 

1247 

1248 if scheme == 'elastic5': 

1249 ioff = 0 

1250 

1251 cat = num.concatenate 

1252 rep = num.repeat 

1253 

1254 w_n = cat((cb*f0, cb*f1, -sb*f2)) 

1255 g_n = ioff + rep((0, 1, 2), n) 

1256 w_e = cat((sb*f0, sb*f1, cb*f2)) 

1257 g_e = ioff + rep((0, 1, 2), n) 

1258 w_d = cat((f0, f1)) 

1259 g_d = ioff + rep((3, 4), n) 

1260 

1261 return ( 

1262 ('displacement.n', w_n, g_n), 

1263 ('displacement.e', w_e, g_e), 

1264 ('displacement.d', w_d, g_d), 

1265 ) 

1266 

1267 @classmethod 

1268 def combine(cls, sources, **kwargs): 

1269 ''' 

1270 Combine several discretized source models. 

1271 

1272 Concatenenates all point sources in the given discretized ``sources``. 

1273 Care must be taken when using this function that the external amplitude 

1274 factors and reference times of the parameterized (undiscretized) 

1275 sources match or are accounted for. 

1276 ''' 

1277 

1278 if 'forces' not in kwargs: 

1279 kwargs['forces'] = num.vstack([s.forces for s in sources]) 

1280 

1281 return super(DiscretizedSFSource, cls).combine(sources, **kwargs) 

1282 

1283 def moments(self): 

1284 return num.sum(self.forces**2, axis=1) 

1285 

1286 def centroid(self): 

1287 from pyrocko.gf.seismosizer import SFSource 

1288 time, lat, lon, north_shift, east_shift, depth = \ 

1289 self.centroid_position() 

1290 

1291 fn, fe, fd = map(float, num.sum(self.forces, axis=0)) 

1292 return SFSource( 

1293 time=time, 

1294 lat=lat, 

1295 lon=lon, 

1296 north_shift=north_shift, 

1297 east_shift=east_shift, 

1298 depth=depth, 

1299 fn=fn, 

1300 fe=fe, 

1301 fd=fd) 

1302 

1303 

1304class DiscretizedMTSource(DiscretizedSource): 

1305 m6s = Array.T( 

1306 shape=(None, 6), dtype=float, 

1307 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)') 

1308 

1309 provided_schemes = ( 

1310 'elastic8', 

1311 'elastic10', 

1312 'elastic18', 

1313 'rotational8', 

1314 ) 

1315 

1316 def get_source_terms(self, scheme): 

1317 self.check_scheme(scheme) 

1318 return self.m6s 

1319 

1320 def make_weights(self, receiver, scheme): 

1321 self.check_scheme(scheme) 

1322 

1323 m6s = self.m6s 

1324 n = m6s.shape[0] 

1325 rep = num.repeat 

1326 

1327 if scheme == 'elastic18': 

1328 w_n = m6s.flatten() 

1329 g_n = num.tile(num.arange(0, 6), n) 

1330 w_e = m6s.flatten() 

1331 g_e = num.tile(num.arange(6, 12), n) 

1332 w_d = m6s.flatten() 

1333 g_d = num.tile(num.arange(12, 18), n) 

1334 

1335 else: 

1336 azis, bazis = self.azibazis_to(receiver) 

1337 

1338 sa = num.sin(azis*d2r) 

1339 ca = num.cos(azis*d2r) 

1340 s2a = num.sin(2.*azis*d2r) 

1341 c2a = num.cos(2.*azis*d2r) 

1342 sb = num.sin(bazis*d2r-num.pi) 

1343 cb = num.cos(bazis*d2r-num.pi) 

1344 

1345 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a 

1346 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa 

1347 f2 = m6s[:, 2] 

1348 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a 

1349 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa 

1350 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a 

1351 

1352 cat = num.concatenate 

1353 

1354 if scheme == 'elastic8': 

1355 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4)) 

1356 g_n = rep((0, 1, 2, 3, 4), n) 

1357 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4)) 

1358 g_e = rep((0, 1, 2, 3, 4), n) 

1359 w_d = cat((f0, f1, f2)) 

1360 g_d = rep((5, 6, 7), n) 

1361 

1362 elif scheme == 'elastic10': 

1363 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4)) 

1364 g_n = rep((0, 1, 2, 8, 3, 4), n) 

1365 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4)) 

1366 g_e = rep((0, 1, 2, 8, 3, 4), n) 

1367 w_d = cat((f0, f1, f2, f5)) 

1368 g_d = rep((5, 6, 7, 9), n) 

1369 

1370 elif scheme == 'rotational18': 

1371 w_n = cat((cb*f3, cb*f4, -sb*f0, -sb*f1, -sb*f2, -sb*f5)) 

1372 g_n = rep((0, 1, 2, 3, 4, 7), n) 

1373 w_n = cat((sb*f3, sb*f4, cb*f0, cb*f1, cb*f2, cb*f5)) 

1374 g_e = rep((0, 1, 2, 3, 4, 7), n) 

1375 w_d = cat((f3, f4)) 

1376 g_d = rep((5, 6), n) 

1377 

1378 else: 

1379 assert False 

1380 

1381 return ( 

1382 ('displacement.n', w_n, g_n), 

1383 ('displacement.e', w_e, g_e), 

1384 ('displacement.d', w_d, g_d), 

1385 ) 

1386 

1387 def split(self): 

1388 from pyrocko.gf.seismosizer import MTSource 

1389 sources = [] 

1390 for i in range(self.nelements): 

1391 lat, lon, north_shift, east_shift = self.element_coords(i) 

1392 sources.append(MTSource( 

1393 time=float(self.times[i]), 

1394 lat=lat, 

1395 lon=lon, 

1396 north_shift=north_shift, 

1397 east_shift=east_shift, 

1398 depth=float(self.depths[i]), 

1399 m6=self.m6s[i])) 

1400 

1401 return sources 

1402 

1403 def moments(self): 

1404 moments = num.array( 

1405 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6)) 

1406 for m6 in self.m6s]) 

1407 return num.linalg.norm(moments, axis=1) / num.sqrt(2.) 

1408 

1409 def get_moment_rate(self, deltat=None): 

1410 moments = self.moments() 

1411 times = self.times 

1412 times -= times.min() 

1413 

1414 t_max = times.max() 

1415 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat 

1416 mom_times[mom_times > t_max] = t_max 

1417 

1418 # Right open histrogram bins 

1419 mom, _ = num.histogram( 

1420 times, 

1421 bins=mom_times, 

1422 weights=moments) 

1423 

1424 deltat = num.diff(mom_times) 

1425 mom_rate = mom / deltat 

1426 

1427 return mom_rate, mom_times[1:] 

1428 

1429 def centroid(self): 

1430 from pyrocko.gf.seismosizer import MTSource 

1431 time, lat, lon, north_shift, east_shift, depth = \ 

1432 self.centroid_position() 

1433 

1434 return MTSource( 

1435 time=time, 

1436 lat=lat, 

1437 lon=lon, 

1438 north_shift=north_shift, 

1439 east_shift=east_shift, 

1440 depth=depth, 

1441 m6=num.sum(self.m6s, axis=0)) 

1442 

1443 @classmethod 

1444 def combine(cls, sources, **kwargs): 

1445 ''' 

1446 Combine several discretized source models. 

1447 

1448 Concatenenates all point sources in the given discretized ``sources``. 

1449 Care must be taken when using this function that the external amplitude 

1450 factors and reference times of the parameterized (undiscretized) 

1451 sources match or are accounted for. 

1452 ''' 

1453 

1454 if 'm6s' not in kwargs: 

1455 kwargs['m6s'] = num.vstack([s.m6s for s in sources]) 

1456 

1457 return super(DiscretizedMTSource, cls).combine(sources, **kwargs) 

1458 

1459 

1460class DiscretizedPorePressureSource(DiscretizedSource): 

1461 pp = Array.T(shape=(None,), dtype=float) 

1462 

1463 provided_schemes = ( 

1464 'poroelastic10', 

1465 ) 

1466 

1467 def get_source_terms(self, scheme): 

1468 self.check_scheme(scheme) 

1469 return self.pp[:, num.newaxis].copy() 

1470 

1471 def make_weights(self, receiver, scheme): 

1472 self.check_scheme(scheme) 

1473 

1474 azis, bazis = self.azibazis_to(receiver) 

1475 

1476 sb = num.sin(bazis*d2r-num.pi) 

1477 cb = num.cos(bazis*d2r-num.pi) 

1478 

1479 pp = self.pp 

1480 n = bazis.size 

1481 

1482 w_un = cb*pp 

1483 g_un = filledi(1, n) 

1484 w_ue = sb*pp 

1485 g_ue = filledi(1, n) 

1486 w_ud = pp 

1487 g_ud = filledi(0, n) 

1488 

1489 w_tn = cb*pp 

1490 g_tn = filledi(6, n) 

1491 w_te = sb*pp 

1492 g_te = filledi(6, n) 

1493 

1494 w_pp = pp 

1495 g_pp = filledi(7, n) 

1496 

1497 w_dvn = cb*pp 

1498 g_dvn = filledi(9, n) 

1499 w_dve = sb*pp 

1500 g_dve = filledi(9, n) 

1501 w_dvd = pp 

1502 g_dvd = filledi(8, n) 

1503 

1504 return ( 

1505 ('displacement.n', w_un, g_un), 

1506 ('displacement.e', w_ue, g_ue), 

1507 ('displacement.d', w_ud, g_ud), 

1508 ('vertical_tilt.n', w_tn, g_tn), 

1509 ('vertical_tilt.e', w_te, g_te), 

1510 ('pore_pressure', w_pp, g_pp), 

1511 ('darcy_velocity.n', w_dvn, g_dvn), 

1512 ('darcy_velocity.e', w_dve, g_dve), 

1513 ('darcy_velocity.d', w_dvd, g_dvd), 

1514 ) 

1515 

1516 def moments(self): 

1517 return self.pp 

1518 

1519 def centroid(self): 

1520 from pyrocko.gf.seismosizer import PorePressurePointSource 

1521 time, lat, lon, north_shift, east_shift, depth = \ 

1522 self.centroid_position() 

1523 

1524 return PorePressurePointSource( 

1525 time=time, 

1526 lat=lat, 

1527 lon=lon, 

1528 north_shift=north_shift, 

1529 east_shift=east_shift, 

1530 depth=depth, 

1531 pp=float(num.sum(self.pp))) 

1532 

1533 @classmethod 

1534 def combine(cls, sources, **kwargs): 

1535 ''' 

1536 Combine several discretized source models. 

1537 

1538 Concatenenates all point sources in the given discretized ``sources``. 

1539 Care must be taken when using this function that the external amplitude 

1540 factors and reference times of the parameterized (undiscretized) 

1541 sources match or are accounted for. 

1542 ''' 

1543 

1544 if 'pp' not in kwargs: 

1545 kwargs['pp'] = num.concatenate([s.pp for s in sources]) 

1546 

1547 return super(DiscretizedPorePressureSource, cls).combine(sources, 

1548 **kwargs) 

1549 

1550 

1551class Region(Object): 

1552 name = String.T(optional=True) 

1553 

1554 

1555class RectangularRegion(Region): 

1556 lat_min = Float.T() 

1557 lat_max = Float.T() 

1558 lon_min = Float.T() 

1559 lon_max = Float.T() 

1560 

1561 

1562class CircularRegion(Region): 

1563 lat = Float.T() 

1564 lon = Float.T() 

1565 radius = Float.T() 

1566 

1567 

1568class Config(Object): 

1569 ''' 

1570 Green's function store meta information. 

1571 

1572 Currently implemented :py:class:`~pyrocko.gf.store.Store` 

1573 configuration types are: 

1574 

1575 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or 

1576 spherical symmetry, 1D earth model, single receiver depth 

1577 

1578 * Problem is invariant to horizontal translations and rotations around 

1579 vertical axis. 

1580 * All receivers must be at the same depth (e.g. at the surface) 

1581 * High level index variables: ``(source_depth, receiver_distance, 

1582 component)`` 

1583 

1584 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or 

1585 spherical symmetry, 1D earth model, variable receiver depth 

1586 

1587 * Symmetries like in Type A but has additional index for receiver depth 

1588 * High level index variables: ``(source_depth, receiver_distance, 

1589 receiver_depth, component)`` 

1590 

1591 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical 

1592 constraints but fixed receiver positions 

1593 

1594 * Cartesian source volume around a reference point 

1595 * High level index variables: ``(ireceiver, source_depth, 

1596 source_east_shift, source_north_shift, component)`` 

1597 ''' 

1598 

1599 id = StringID.T( 

1600 help='Name of the store. May consist of upper and lower-case letters, ' 

1601 'digits, dots and underscores. The name must start with a ' 

1602 'letter.') 

1603 

1604 derived_from_id = StringID.T( 

1605 optional=True, 

1606 help='Name of the original store, if this store has been derived from ' 

1607 'another one (e.g. extracted subset).') 

1608 

1609 version = String.T( 

1610 default='1.0', 

1611 optional=True, 

1612 help='User-defined version string. Use <major>.<minor> format.') 

1613 

1614 modelling_code_id = StringID.T( 

1615 optional=True, 

1616 help='Identifier of the backend used to compute the store.') 

1617 

1618 author = Unicode.T( 

1619 optional=True, 

1620 help='Comma-separated list of author names.') 

1621 

1622 author_email = String.T( 

1623 optional=True, 

1624 help="Author's contact email address.") 

1625 

1626 created_time = Timestamp.T( 

1627 optional=True, 

1628 help='Time of creation of the store.') 

1629 

1630 regions = List.T( 

1631 Region.T(), 

1632 help='Geographical regions for which the store is representative.') 

1633 

1634 scope_type = ScopeType.T( 

1635 optional=True, 

1636 help='Distance range scope of the store ' 

1637 '(%s).' % fmt_choices(ScopeType)) 

1638 

1639 waveform_type = WaveType.T( 

1640 optional=True, 

1641 help='Wave type stored (%s).' % fmt_choices(WaveType)) 

1642 

1643 nearfield_terms = NearfieldTermsType.T( 

1644 optional=True, 

1645 help='Information about the inclusion of near-field terms in the ' 

1646 'modelling (%s).' % fmt_choices(NearfieldTermsType)) 

1647 

1648 description = String.T( 

1649 optional=True, 

1650 help='Free form textual description of the GF store.') 

1651 

1652 references = List.T( 

1653 Reference.T(), 

1654 help='Reference list to cite the modelling code, earth model or ' 

1655 'related work.') 

1656 

1657 earthmodel_1d = Earthmodel1D.T( 

1658 optional=True, 

1659 help='Layered earth model in ND (named discontinuity) format.') 

1660 

1661 earthmodel_receiver_1d = Earthmodel1D.T( 

1662 optional=True, 

1663 help='Receiver-side layered earth model in ND format.') 

1664 

1665 can_interpolate_source = Bool.T( 

1666 optional=True, 

1667 help='Hint to indicate if the spatial sampling of the store is dense ' 

1668 'enough for multi-linear interpolation at the source.') 

1669 

1670 can_interpolate_receiver = Bool.T( 

1671 optional=True, 

1672 help='Hint to indicate if the spatial sampling of the store is dense ' 

1673 'enough for multi-linear interpolation at the receiver.') 

1674 

1675 frequency_min = Float.T( 

1676 optional=True, 

1677 help='Hint to indicate the lower bound of valid frequencies [Hz].') 

1678 

1679 frequency_max = Float.T( 

1680 optional=True, 

1681 help='Hint to indicate the upper bound of valid frequencies [Hz].') 

1682 

1683 sample_rate = Float.T( 

1684 optional=True, 

1685 help='Sample rate of the GF store [Hz].') 

1686 

1687 factor = Float.T( 

1688 default=1.0, 

1689 help='Gain value, factored out of the stored GF samples. ' 

1690 '(may not work properly, keep at 1.0).', 

1691 optional=True) 

1692 

1693 component_scheme = ComponentScheme.T( 

1694 default='elastic10', 

1695 help='GF component scheme (%s).' % fmt_choices(ComponentScheme)) 

1696 

1697 stored_quantity = QuantityType.T( 

1698 optional=True, 

1699 help='Physical quantity of stored values (%s). If not given, a ' 

1700 'default is used based on the GF component scheme. The default ' 

1701 'for the ``"elastic*"`` family of component schemes is ' 

1702 '``"displacement"``.' % fmt_choices(QuantityType)) 

1703 

1704 tabulated_phases = List.T( 

1705 TPDef.T(), 

1706 help='Mapping of phase names to phase definitions, for which travel ' 

1707 'time tables are available in the GF store.') 

1708 

1709 ncomponents = Int.T( 

1710 optional=True, 

1711 help='Number of GF components. Use :gattr:`component_scheme` instead.') 

1712 

1713 uuid = String.T( 

1714 optional=True, 

1715 help='Heuristic hash value which can be used to uniquely identify the ' 

1716 'GF store for practical purposes.') 

1717 

1718 reference = String.T( 

1719 optional=True, 

1720 help="Store reference name composed of the store's :gattr:`id` and " 

1721 'the first six letters of its :gattr:`uuid`.') 

1722 

1723 def __init__(self, **kwargs): 

1724 self._do_auto_updates = False 

1725 Object.__init__(self, **kwargs) 

1726 self._index_function = None 

1727 self._indices_function = None 

1728 self._vicinity_function = None 

1729 self.validate(regularize=True, depth=1) 

1730 self._do_auto_updates = True 

1731 self.update() 

1732 

1733 def check_ncomponents(self): 

1734 ncomponents = component_scheme_to_description[ 

1735 self.component_scheme].ncomponents 

1736 

1737 if self.ncomponents is None: 

1738 self.ncomponents = ncomponents 

1739 elif ncomponents != self.ncomponents: 

1740 raise InvalidNComponents( 

1741 'ncomponents=%i incompatible with component_scheme="%s"' % ( 

1742 self.ncomponents, self.component_scheme)) 

1743 

1744 def __setattr__(self, name, value): 

1745 Object.__setattr__(self, name, value) 

1746 try: 

1747 self.T.get_property(name) 

1748 if self._do_auto_updates: 

1749 self.update() 

1750 

1751 except ValueError: 

1752 pass 

1753 

1754 def update(self): 

1755 self.check_ncomponents() 

1756 self._update() 

1757 self._make_index_functions() 

1758 

1759 def irecord(self, *args): 

1760 return self._index_function(*args) 

1761 

1762 def irecords(self, *args): 

1763 return self._indices_function(*args) 

1764 

1765 def vicinity(self, *args): 

1766 return self._vicinity_function(*args) 

1767 

1768 def vicinities(self, *args): 

1769 return self._vicinities_function(*args) 

1770 

1771 def grid_interpolation_coefficients(self, *args): 

1772 return self._grid_interpolation_coefficients(*args) 

1773 

1774 def nodes(self, level=None, minlevel=None): 

1775 return nodes(self.coords[minlevel:level]) 

1776 

1777 def iter_nodes(self, level=None, minlevel=None): 

1778 return nditer_outer(self.coords[minlevel:level]) 

1779 

1780 def iter_extraction(self, gdef, level=None): 

1781 i = 0 

1782 arrs = [] 

1783 ntotal = 1 

1784 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas): 

1785 if gdef and len(gdef) > i: 

1786 sssn = gdef[i] 

1787 else: 

1788 sssn = (None,)*4 

1789 

1790 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc)))) 

1791 ntotal *= len(arr) 

1792 

1793 arrs.append(arr) 

1794 i += 1 

1795 

1796 arrs.append(self.coords[-1]) 

1797 return nditer_outer(arrs[:level]) 

1798 

1799 def make_sum_params(self, source, receiver, implementation='c', 

1800 nthreads=0): 

1801 assert implementation in ['c', 'python'] 

1802 

1803 out = [] 

1804 delays = source.times 

1805 for comp, weights, icomponents in source.make_weights( 

1806 receiver, 

1807 self.component_scheme): 

1808 

1809 weights *= self.factor 

1810 

1811 args = self.make_indexing_args(source, receiver, icomponents) 

1812 delays_expanded = num.tile(delays, icomponents.size//delays.size) 

1813 out.append((comp, args, delays_expanded, weights)) 

1814 

1815 return out 

1816 

1817 def short_info(self): 

1818 raise NotImplementedError('should be implemented in subclass') 

1819 

1820 def get_shear_moduli(self, lat, lon, points, 

1821 interpolation=None): 

1822 ''' 

1823 Get shear moduli at given points from contained velocity model. 

1824 

1825 :param lat: surface origin for coordinate system of ``points`` 

1826 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1827 a point ``(north, east, depth)``, relative to origin at 

1828 ``(lat, lon)`` 

1829 :param interpolation: interpolation method. Choose from 

1830 ``('nearest_neighbor', 'multilinear')`` 

1831 :returns: NumPy array of length N with extracted shear moduli at each 

1832 point 

1833 

1834 The default implementation retrieves and interpolates the shear moduli 

1835 from the contained 1D velocity profile. 

1836 ''' 

1837 return self.get_material_property(lat, lon, points, 

1838 parameter='shear_moduli', 

1839 interpolation=interpolation) 

1840 

1841 def get_lambda_moduli(self, lat, lon, points, 

1842 interpolation=None): 

1843 ''' 

1844 Get lambda moduli at given points from contained velocity model. 

1845 

1846 :param lat: surface origin for coordinate system of ``points`` 

1847 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1848 a point ``(north, east, depth)``, relative to origin at 

1849 ``(lat, lon)`` 

1850 :param interpolation: interpolation method. Choose from 

1851 ``('nearest_neighbor', 'multilinear')`` 

1852 :returns: NumPy array of length N with extracted shear moduli at each 

1853 point 

1854 

1855 The default implementation retrieves and interpolates the lambda moduli 

1856 from the contained 1D velocity profile. 

1857 ''' 

1858 return self.get_material_property(lat, lon, points, 

1859 parameter='lambda_moduli', 

1860 interpolation=interpolation) 

1861 

1862 def get_bulk_moduli(self, lat, lon, points, 

1863 interpolation=None): 

1864 ''' 

1865 Get bulk moduli at given points from contained velocity model. 

1866 

1867 :param lat: surface origin for coordinate system of ``points`` 

1868 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1869 a point ``(north, east, depth)``, relative to origin at 

1870 ``(lat, lon)`` 

1871 :param interpolation: interpolation method. Choose from 

1872 ``('nearest_neighbor', 'multilinear')`` 

1873 :returns: NumPy array of length N with extracted shear moduli at each 

1874 point 

1875 

1876 The default implementation retrieves and interpolates the lambda moduli 

1877 from the contained 1D velocity profile. 

1878 ''' 

1879 lambda_moduli = self.get_material_property( 

1880 lat, lon, points, parameter='lambda_moduli', 

1881 interpolation=interpolation) 

1882 shear_moduli = self.get_material_property( 

1883 lat, lon, points, parameter='shear_moduli', 

1884 interpolation=interpolation) 

1885 return lambda_moduli + (2 / 3) * shear_moduli 

1886 

1887 def get_vs(self, lat, lon, points, interpolation=None): 

1888 ''' 

1889 Get Vs at given points from contained velocity model. 

1890 

1891 :param lat: surface origin for coordinate system of ``points`` 

1892 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1893 a point ``(north, east, depth)``, relative to origin at 

1894 ``(lat, lon)`` 

1895 :param interpolation: interpolation method. Choose from 

1896 ``('nearest_neighbor', 'multilinear')`` 

1897 :returns: NumPy array of length N with extracted shear moduli at each 

1898 point 

1899 

1900 The default implementation retrieves and interpolates Vs 

1901 from the contained 1D velocity profile. 

1902 ''' 

1903 return self.get_material_property(lat, lon, points, 

1904 parameter='vs', 

1905 interpolation=interpolation) 

1906 

1907 def get_vp(self, lat, lon, points, interpolation=None): 

1908 ''' 

1909 Get Vp at given points from contained velocity model. 

1910 

1911 :param lat: surface origin for coordinate system of ``points`` 

1912 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1913 a point ``(north, east, depth)``, relative to origin at 

1914 ``(lat, lon)`` 

1915 :param interpolation: interpolation method. Choose from 

1916 ``('nearest_neighbor', 'multilinear')`` 

1917 :returns: NumPy array of length N with extracted shear moduli at each 

1918 point 

1919 

1920 The default implementation retrieves and interpolates Vp 

1921 from the contained 1D velocity profile. 

1922 ''' 

1923 return self.get_material_property(lat, lon, points, 

1924 parameter='vp', 

1925 interpolation=interpolation) 

1926 

1927 def get_rho(self, lat, lon, points, interpolation=None): 

1928 ''' 

1929 Get rho at given points from contained velocity model. 

1930 

1931 :param lat: surface origin for coordinate system of ``points`` 

1932 :param points: NumPy array of shape ``(N, 3)``, where each row is 

1933 a point ``(north, east, depth)``, relative to origin at 

1934 ``(lat, lon)`` 

1935 :param interpolation: interpolation method. Choose from 

1936 ``('nearest_neighbor', 'multilinear')`` 

1937 :returns: NumPy array of length N with extracted shear moduli at each 

1938 point 

1939 

1940 The default implementation retrieves and interpolates rho 

1941 from the contained 1D velocity profile. 

1942 ''' 

1943 return self.get_material_property(lat, lon, points, 

1944 parameter='rho', 

1945 interpolation=interpolation) 

1946 

1947 def get_material_property(self, lat, lon, points, parameter='vs', 

1948 interpolation=None): 

1949 

1950 if interpolation is None: 

1951 raise TypeError('Interpolation method not defined! available: ' 

1952 'multilinear', 'nearest_neighbor') 

1953 

1954 earthmod = self.earthmodel_1d 

1955 store_depth_profile = self.get_source_depths() 

1956 z_profile = earthmod.profile('z') 

1957 

1958 if parameter == 'vs': 

1959 vs_profile = earthmod.profile('vs') 

1960 profile = num.interp( 

1961 store_depth_profile, z_profile, vs_profile) 

1962 

1963 elif parameter == 'vp': 

1964 vp_profile = earthmod.profile('vp') 

1965 profile = num.interp( 

1966 store_depth_profile, z_profile, vp_profile) 

1967 

1968 elif parameter == 'rho': 

1969 rho_profile = earthmod.profile('rho') 

1970 

1971 profile = num.interp( 

1972 store_depth_profile, z_profile, rho_profile) 

1973 

1974 elif parameter == 'shear_moduli': 

1975 vs_profile = earthmod.profile('vs') 

1976 rho_profile = earthmod.profile('rho') 

1977 

1978 store_vs_profile = num.interp( 

1979 store_depth_profile, z_profile, vs_profile) 

1980 store_rho_profile = num.interp( 

1981 store_depth_profile, z_profile, rho_profile) 

1982 

1983 profile = num.power(store_vs_profile, 2) * store_rho_profile 

1984 

1985 elif parameter == 'lambda_moduli': 

1986 vs_profile = earthmod.profile('vs') 

1987 vp_profile = earthmod.profile('vp') 

1988 rho_profile = earthmod.profile('rho') 

1989 

1990 store_vs_profile = num.interp( 

1991 store_depth_profile, z_profile, vs_profile) 

1992 store_vp_profile = num.interp( 

1993 store_depth_profile, z_profile, vp_profile) 

1994 store_rho_profile = num.interp( 

1995 store_depth_profile, z_profile, rho_profile) 

1996 

1997 profile = store_rho_profile * ( 

1998 num.power(store_vp_profile, 2) - 

1999 num.power(store_vs_profile, 2) * 2) 

2000 else: 

2001 raise TypeError( 

2002 'parameter %s not available' % parameter) 

2003 

2004 if interpolation == 'multilinear': 

2005 kind = 'linear' 

2006 elif interpolation == 'nearest_neighbor': 

2007 kind = 'nearest' 

2008 else: 

2009 raise TypeError( 

2010 'Interpolation method %s not available' % interpolation) 

2011 

2012 interpolator = interp1d(store_depth_profile, profile, kind=kind) 

2013 

2014 try: 

2015 return interpolator(points[:, 2]) 

2016 except ValueError: 

2017 raise OutOfBounds() 

2018 

2019 def is_static(self): 

2020 for code in ('psgrn_pscmp', 'poel'): 

2021 if self.modelling_code_id.startswith(code): 

2022 return True 

2023 return False 

2024 

2025 def is_dynamic(self): 

2026 return not self.is_static() 

2027 

2028 def get_source_depths(self): 

2029 raise NotImplementedError('must be implemented in subclass') 

2030 

2031 def get_tabulated_phase(self, phase_id): 

2032 ''' 

2033 Get tabulated phase definition. 

2034 ''' 

2035 

2036 for pdef in self.tabulated_phases: 

2037 if pdef.id == phase_id: 

2038 return pdef 

2039 

2040 raise StoreError('No such phase: %s' % phase_id) 

2041 

2042 def fix_ttt_holes(self, sptree, mode): 

2043 raise StoreError( 

2044 'Cannot fix travel time table holes in GF stores of type %s.' 

2045 % self.short_type) 

2046 

2047 @property 

2048 def effective_stored_quantity(self): 

2049 return self.stored_quantity if self.stored_quantity is not None else \ 

2050 component_scheme_to_description[self.component_scheme] \ 

2051 .default_stored_quantity 

2052 

2053 

2054class ConfigTypeA(Config): 

2055 ''' 

2056 Cylindrical symmetry, 1D earth model, single receiver depth 

2057 

2058 * Problem is invariant to horizontal translations and rotations around 

2059 vertical axis. 

2060 

2061 * All receivers must be at the same depth (e.g. at the surface) 

2062 High level index variables: ``(source_depth, distance, 

2063 component)`` 

2064 

2065 * The ``distance`` is the surface distance between source and receiver 

2066 points. 

2067 ''' 

2068 

2069 receiver_depth = Float.T( 

2070 default=0.0, 

2071 help='Fixed receiver depth [m].') 

2072 

2073 source_depth_min = Float.T( 

2074 help='Minimum source depth [m].') 

2075 

2076 source_depth_max = Float.T( 

2077 help='Maximum source depth [m].') 

2078 

2079 source_depth_delta = Float.T( 

2080 help='Grid spacing of source depths [m]') 

2081 

2082 distance_min = Float.T( 

2083 help='Minimum source-receiver surface distance [m].') 

2084 

2085 distance_max = Float.T( 

2086 help='Maximum source-receiver surface distance [m].') 

2087 

2088 distance_delta = Float.T( 

2089 help='Grid spacing of source-receiver surface distance [m].') 

2090 

2091 fd_distance_delta = Float.T( 

2092 optional=True, 

2093 help='Finite differences interval for FD stores [m].') 

2094 

2095 short_type = 'A' 

2096 

2097 provided_schemes = [ 

2098 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10', 

2099 'rotational8'] 

2100 

2101 def get_surface_distance(self, args): 

2102 return args[1] 

2103 

2104 def get_distance(self, args): 

2105 return math.sqrt(args[0]**2 + args[1]**2) 

2106 

2107 def get_source_depth(self, args): 

2108 return args[0] 

2109 

2110 def get_source_depths(self): 

2111 return self.coords[0] 

2112 

2113 def get_receiver_depth(self, args): 

2114 return self.receiver_depth 

2115 

2116 def _update(self): 

2117 self.mins = num.array( 

2118 [self.source_depth_min, self.distance_min], dtype=float) 

2119 self.maxs = num.array( 

2120 [self.source_depth_max, self.distance_max], dtype=float) 

2121 self.deltas = num.array( 

2122 [self.source_depth_delta, self.distance_delta], 

2123 dtype=float) 

2124 self.ns = num.floor((self.maxs - self.mins) / self.deltas + 

2125 vicinity_eps).astype(int) + 1 

2126 self.effective_maxs = self.mins + self.deltas * (self.ns - 1) 

2127 self.deltat = 1.0/self.sample_rate 

2128 self.nrecords = num.prod(self.ns) * self.ncomponents 

2129 self.coords = tuple(num.linspace(mi, ma, n) for 

2130 (mi, ma, n) in 

2131 zip(self.mins, self.effective_maxs, self.ns)) + \ 

2132 (num.arange(self.ncomponents),) 

2133 

2134 self.nsource_depths, self.ndistances = self.ns 

2135 

2136 def _make_index_functions(self): 

2137 

2138 amin, bmin = self.mins 

2139 da, db = self.deltas 

2140 na, nb = self.ns 

2141 

2142 ng = self.ncomponents 

2143 

2144 def index_function(a, b, ig): 

2145 ia = int(round((a - amin) / da)) 

2146 ib = int(round((b - bmin) / db)) 

2147 try: 

2148 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng)) 

2149 except ValueError: 

2150 raise OutOfBounds() 

2151 

2152 def indices_function(a, b, ig): 

2153 ia = num.round((a - amin) / da).astype(int) 

2154 ib = num.round((b - bmin) / db).astype(int) 

2155 try: 

2156 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng)) 

2157 except ValueError: 

2158 for ia_, ib_, ig_ in zip(ia, ib, ig): 

2159 try: 

2160 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng)) 

2161 except ValueError: 

2162 raise OutOfBounds() 

2163 

2164 def grid_interpolation_coefficients(a, b): 

2165 ias = indi12((a - amin) / da, na) 

2166 ibs = indi12((b - bmin) / db, nb) 

2167 return ias, ibs 

2168 

2169 def vicinity_function(a, b, ig): 

2170 ias, ibs = grid_interpolation_coefficients(a, b) 

2171 

2172 if not (0 <= ig < ng): 

2173 raise OutOfBounds() 

2174 

2175 indis = [] 

2176 weights = [] 

2177 for ia, va in ias: 

2178 iia = ia*nb*ng 

2179 for ib, vb in ibs: 

2180 indis.append(iia + ib*ng + ig) 

2181 weights.append(va*vb) 

2182 

2183 return num.array(indis), num.array(weights) 

2184 

2185 def vicinities_function(a, b, ig): 

2186 

2187 xa = (a - amin) / da 

2188 xb = (b - bmin) / db 

2189 

2190 xa_fl = num.floor(xa) 

2191 xa_ce = num.ceil(xa) 

2192 xb_fl = num.floor(xb) 

2193 xb_ce = num.ceil(xb) 

2194 va_fl = 1.0 - (xa - xa_fl) 

2195 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl) 

2196 vb_fl = 1.0 - (xb - xb_fl) 

2197 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl) 

2198 

2199 ia_fl = xa_fl.astype(int) 

2200 ia_ce = xa_ce.astype(int) 

2201 ib_fl = xb_fl.astype(int) 

2202 ib_ce = xb_ce.astype(int) 

2203 

2204 if num.any(ia_fl < 0) or num.any(ia_fl >= na): 

2205 raise OutOfBounds() 

2206 

2207 if num.any(ia_ce < 0) or num.any(ia_ce >= na): 

2208 raise OutOfBounds() 

2209 

2210 if num.any(ib_fl < 0) or num.any(ib_fl >= nb): 

2211 raise OutOfBounds() 

2212 

2213 if num.any(ib_ce < 0) or num.any(ib_ce >= nb): 

2214 raise OutOfBounds() 

2215 

2216 irecords = num.empty(a.size*4, dtype=int) 

2217 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig 

2218 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig 

2219 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig 

2220 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig 

2221 

2222 weights = num.empty(a.size*4, dtype=float) 

2223 weights[0::4] = va_fl * vb_fl 

2224 weights[1::4] = va_ce * vb_fl 

2225 weights[2::4] = va_fl * vb_ce 

2226 weights[3::4] = va_ce * vb_ce 

2227 

2228 return irecords, weights 

2229 

2230 self._index_function = index_function 

2231 self._indices_function = indices_function 

2232 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2233 self._vicinity_function = vicinity_function 

2234 self._vicinities_function = vicinities_function 

2235 

2236 def make_indexing_args(self, source, receiver, icomponents): 

2237 nc = icomponents.size 

2238 dists = source.distances_to(receiver) 

2239 n = dists.size 

2240 return (num.tile(source.depths, nc//n), 

2241 num.tile(dists, nc//n), 

2242 icomponents) 

2243 

2244 def make_indexing_args1(self, source, receiver): 

2245 return (source.depth, source.distance_to(receiver)) 

2246 

2247 @property 

2248 def short_extent(self): 

2249 return '%g:%g:%g x %g:%g:%g' % ( 

2250 self.source_depth_min/km, 

2251 self.source_depth_max/km, 

2252 self.source_depth_delta/km, 

2253 self.distance_min/km, 

2254 self.distance_max/km, 

2255 self.distance_delta/km) 

2256 

2257 def fix_ttt_holes(self, sptree, mode): 

2258 from pyrocko import eikonal_ext, spit 

2259 

2260 nodes = self.nodes(level=-1) 

2261 

2262 delta = self.deltas[-1] 

2263 assert num.all(delta == self.deltas) 

2264 

2265 nsources, ndistances = self.ns 

2266 

2267 points = num.zeros((nodes.shape[0], 3)) 

2268 points[:, 0] = nodes[:, 1] 

2269 points[:, 2] = nodes[:, 0] 

2270 

2271 speeds = self.get_material_property( 

2272 0., 0., points, 

2273 parameter='vp' if mode == cake.P else 'vs', 

2274 interpolation='multilinear') 

2275 

2276 speeds = speeds.reshape((nsources, ndistances)) 

2277 

2278 times = sptree.interpolate_many(nodes) 

2279 

2280 times[num.isnan(times)] = -1. 

2281 times = times.reshape(speeds.shape) 

2282 

2283 try: 

2284 eikonal_ext.eikonal_solver_fmm_cartesian( 

2285 speeds, times, delta) 

2286 except eikonal_ext.EikonalExtError as e: 

2287 if str(e).endswith('please check results'): 

2288 logger.debug( 

2289 'Got a warning from eikonal solver ' 

2290 '- may be ok...') 

2291 else: 

2292 raise 

2293 

2294 def func(x): 

2295 ibs, ics = \ 

2296 self.grid_interpolation_coefficients(*x) 

2297 

2298 t = 0 

2299 for ib, vb in ibs: 

2300 for ic, vc in ics: 

2301 t += times[ib, ic] * vb * vc 

2302 

2303 return t 

2304 

2305 return spit.SPTree( 

2306 f=func, 

2307 ftol=sptree.ftol, 

2308 xbounds=sptree.xbounds, 

2309 xtols=sptree.xtols) 

2310 

2311 

2312class ConfigTypeB(Config): 

2313 ''' 

2314 Cylindrical symmetry, 1D earth model, variable receiver depth 

2315 

2316 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for 

2317 receiver depth 

2318 

2319 * High level index variables: ``(receiver_depth, source_depth, 

2320 receiver_distance, component)`` 

2321 ''' 

2322 

2323 receiver_depth_min = Float.T( 

2324 help='Minimum receiver depth [m].') 

2325 

2326 receiver_depth_max = Float.T( 

2327 help='Maximum receiver depth [m].') 

2328 

2329 receiver_depth_delta = Float.T( 

2330 help='Grid spacing of receiver depths [m]') 

2331 

2332 source_depth_min = Float.T( 

2333 help='Minimum source depth [m].') 

2334 

2335 source_depth_max = Float.T( 

2336 help='Maximum source depth [m].') 

2337 

2338 source_depth_delta = Float.T( 

2339 help='Grid spacing of source depths [m]') 

2340 

2341 distance_min = Float.T( 

2342 help='Minimum source-receiver surface distance [m].') 

2343 

2344 distance_max = Float.T( 

2345 help='Maximum source-receiver surface distance [m].') 

2346 

2347 distance_delta = Float.T( 

2348 help='Grid spacing of source-receiver surface distances [m].') 

2349 

2350 fd_distance_delta = Float.T( 

2351 optional=True, 

2352 help='Finite differences interval for FD stores [m].') 

2353 

2354 fd_receiver_depth_delta = Float.T( 

2355 optional=True, 

2356 help='Finite differences interval for FD stores [m].') 

2357 

2358 short_type = 'B' 

2359 

2360 provided_schemes = [ 

2361 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10', 

2362 'rotational8'] 

2363 

2364 def get_distance(self, args): 

2365 return math.sqrt((args[1] - args[0])**2 + args[2]**2) 

2366 

2367 def get_surface_distance(self, args): 

2368 return args[2] 

2369 

2370 def get_source_depth(self, args): 

2371 return args[1] 

2372 

2373 def get_receiver_depth(self, args): 

2374 return args[0] 

2375 

2376 def get_source_depths(self): 

2377 return self.coords[1] 

2378 

2379 def _update(self): 

2380 self.mins = num.array([ 

2381 self.receiver_depth_min, 

2382 self.source_depth_min, 

2383 self.distance_min], 

2384 dtype=float) 

2385 

2386 self.maxs = num.array([ 

2387 self.receiver_depth_max, 

2388 self.source_depth_max, 

2389 self.distance_max], 

2390 dtype=float) 

2391 

2392 self.deltas = num.array([ 

2393 self.receiver_depth_delta, 

2394 self.source_depth_delta, 

2395 self.distance_delta], 

2396 dtype=float) 

2397 

2398 self.ns = num.floor((self.maxs - self.mins) / self.deltas + 

2399 vicinity_eps).astype(int) + 1 

2400 self.effective_maxs = self.mins + self.deltas * (self.ns - 1) 

2401 self.deltat = 1.0/self.sample_rate 

2402 self.nrecords = num.prod(self.ns) * self.ncomponents 

2403 self.coords = tuple(num.linspace(mi, ma, n) for 

2404 (mi, ma, n) in 

2405 zip(self.mins, self.effective_maxs, self.ns)) + \ 

2406 (num.arange(self.ncomponents),) 

2407 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns 

2408 

2409 def _make_index_functions(self): 

2410 

2411 amin, bmin, cmin = self.mins 

2412 da, db, dc = self.deltas 

2413 na, nb, nc = self.ns 

2414 ng = self.ncomponents 

2415 

2416 def index_function(a, b, c, ig): 

2417 ia = int(round((a - amin) / da)) 

2418 ib = int(round((b - bmin) / db)) 

2419 ic = int(round((c - cmin) / dc)) 

2420 try: 

2421 return num.ravel_multi_index((ia, ib, ic, ig), 

2422 (na, nb, nc, ng)) 

2423 except ValueError: 

2424 raise OutOfBounds() 

2425 

2426 def indices_function(a, b, c, ig): 

2427 ia = num.round((a - amin) / da).astype(int) 

2428 ib = num.round((b - bmin) / db).astype(int) 

2429 ic = num.round((c - cmin) / dc).astype(int) 

2430 try: 

2431 return num.ravel_multi_index((ia, ib, ic, ig), 

2432 (na, nb, nc, ng)) 

2433 except ValueError: 

2434 raise OutOfBounds() 

2435 

2436 def grid_interpolation_coefficients(a, b, c): 

2437 ias = indi12((a - amin) / da, na) 

2438 ibs = indi12((b - bmin) / db, nb) 

2439 ics = indi12((c - cmin) / dc, nc) 

2440 return ias, ibs, ics 

2441 

2442 def vicinity_function(a, b, c, ig): 

2443 ias, ibs, ics = grid_interpolation_coefficients(a, b, c) 

2444 

2445 if not (0 <= ig < ng): 

2446 raise OutOfBounds() 

2447 

2448 indis = [] 

2449 weights = [] 

2450 for ia, va in ias: 

2451 iia = ia*nb*nc*ng 

2452 for ib, vb in ibs: 

2453 iib = ib*nc*ng 

2454 for ic, vc in ics: 

2455 indis.append(iia + iib + ic*ng + ig) 

2456 weights.append(va*vb*vc) 

2457 

2458 return num.array(indis), num.array(weights) 

2459 

2460 def vicinities_function(a, b, c, ig): 

2461 

2462 xa = (a - amin) / da 

2463 xb = (b - bmin) / db 

2464 xc = (c - cmin) / dc 

2465 

2466 xa_fl = num.floor(xa) 

2467 xa_ce = num.ceil(xa) 

2468 xb_fl = num.floor(xb) 

2469 xb_ce = num.ceil(xb) 

2470 xc_fl = num.floor(xc) 

2471 xc_ce = num.ceil(xc) 

2472 va_fl = 1.0 - (xa - xa_fl) 

2473 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl) 

2474 vb_fl = 1.0 - (xb - xb_fl) 

2475 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl) 

2476 vc_fl = 1.0 - (xc - xc_fl) 

2477 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl) 

2478 

2479 ia_fl = xa_fl.astype(int) 

2480 ia_ce = xa_ce.astype(int) 

2481 ib_fl = xb_fl.astype(int) 

2482 ib_ce = xb_ce.astype(int) 

2483 ic_fl = xc_fl.astype(int) 

2484 ic_ce = xc_ce.astype(int) 

2485 

2486 if num.any(ia_fl < 0) or num.any(ia_fl >= na): 

2487 raise OutOfBounds() 

2488 

2489 if num.any(ia_ce < 0) or num.any(ia_ce >= na): 

2490 raise OutOfBounds() 

2491 

2492 if num.any(ib_fl < 0) or num.any(ib_fl >= nb): 

2493 raise OutOfBounds() 

2494 

2495 if num.any(ib_ce < 0) or num.any(ib_ce >= nb): 

2496 raise OutOfBounds() 

2497 

2498 if num.any(ic_fl < 0) or num.any(ic_fl >= nc): 

2499 raise OutOfBounds() 

2500 

2501 if num.any(ic_ce < 0) or num.any(ic_ce >= nc): 

2502 raise OutOfBounds() 

2503 

2504 irecords = num.empty(a.size*8, dtype=int) 

2505 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig 

2506 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig 

2507 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig 

2508 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig 

2509 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig 

2510 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig 

2511 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig 

2512 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig 

2513 

2514 weights = num.empty(a.size*8, dtype=float) 

2515 weights[0::8] = va_fl * vb_fl * vc_fl 

2516 weights[1::8] = va_ce * vb_fl * vc_fl 

2517 weights[2::8] = va_fl * vb_ce * vc_fl 

2518 weights[3::8] = va_ce * vb_ce * vc_fl 

2519 weights[4::8] = va_fl * vb_fl * vc_ce 

2520 weights[5::8] = va_ce * vb_fl * vc_ce 

2521 weights[6::8] = va_fl * vb_ce * vc_ce 

2522 weights[7::8] = va_ce * vb_ce * vc_ce 

2523 

2524 return irecords, weights 

2525 

2526 self._index_function = index_function 

2527 self._indices_function = indices_function 

2528 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2529 self._vicinity_function = vicinity_function 

2530 self._vicinities_function = vicinities_function 

2531 

2532 def make_indexing_args(self, source, receiver, icomponents): 

2533 nc = icomponents.size 

2534 dists = source.distances_to(receiver) 

2535 n = dists.size 

2536 receiver_depths = num.empty(nc) 

2537 receiver_depths.fill(receiver.depth) 

2538 return (receiver_depths, 

2539 num.tile(source.depths, nc//n), 

2540 num.tile(dists, nc//n), 

2541 icomponents) 

2542 

2543 def make_indexing_args1(self, source, receiver): 

2544 return (receiver.depth, 

2545 source.depth, 

2546 source.distance_to(receiver)) 

2547 

2548 @property 

2549 def short_extent(self): 

2550 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2551 self.receiver_depth_min/km, 

2552 self.receiver_depth_max/km, 

2553 self.receiver_depth_delta/km, 

2554 self.source_depth_min/km, 

2555 self.source_depth_max/km, 

2556 self.source_depth_delta/km, 

2557 self.distance_min/km, 

2558 self.distance_max/km, 

2559 self.distance_delta/km) 

2560 

2561 def fix_ttt_holes(self, sptree, mode): 

2562 from pyrocko import eikonal_ext, spit 

2563 

2564 nodes_sr = self.nodes(minlevel=1, level=-1) 

2565 

2566 delta = self.deltas[-1] 

2567 assert num.all(delta == self.deltas[1:]) 

2568 

2569 nreceivers, nsources, ndistances = self.ns 

2570 

2571 points = num.zeros((nodes_sr.shape[0], 3)) 

2572 points[:, 0] = nodes_sr[:, 1] 

2573 points[:, 2] = nodes_sr[:, 0] 

2574 

2575 speeds = self.get_material_property( 

2576 0., 0., points, 

2577 parameter='vp' if mode == cake.P else 'vs', 

2578 interpolation='multilinear') 

2579 

2580 speeds = speeds.reshape((nsources, ndistances)) 

2581 

2582 receiver_times = [] 

2583 for ireceiver in range(nreceivers): 

2584 nodes = num.hstack([ 

2585 num.full( 

2586 (nodes_sr.shape[0], 1), 

2587 self.coords[0][ireceiver]), 

2588 nodes_sr]) 

2589 

2590 times = sptree.interpolate_many(nodes) 

2591 

2592 times[num.isnan(times)] = -1. 

2593 

2594 times = times.reshape(speeds.shape) 

2595 

2596 try: 

2597 eikonal_ext.eikonal_solver_fmm_cartesian( 

2598 speeds, times, delta) 

2599 except eikonal_ext.EikonalExtError as e: 

2600 if str(e).endswith('please check results'): 

2601 logger.debug( 

2602 'Got a warning from eikonal solver ' 

2603 '- may be ok...') 

2604 else: 

2605 raise 

2606 

2607 receiver_times.append(times) 

2608 

2609 def func(x): 

2610 ias, ibs, ics = \ 

2611 self.grid_interpolation_coefficients(*x) 

2612 

2613 t = 0 

2614 for ia, va in ias: 

2615 times = receiver_times[ia] 

2616 for ib, vb in ibs: 

2617 for ic, vc in ics: 

2618 t += times[ib, ic] * va * vb * vc 

2619 

2620 return t 

2621 

2622 return spit.SPTree( 

2623 f=func, 

2624 ftol=sptree.ftol, 

2625 xbounds=sptree.xbounds, 

2626 xtols=sptree.xtols) 

2627 

2628 

2629class ConfigTypeC(Config): 

2630 ''' 

2631 No symmetrical constraints, one fixed receiver position. 

2632 

2633 * Cartesian 3D source volume around a reference point 

2634 

2635 * High level index variables: ``(source_depth, 

2636 source_east_shift, source_north_shift, component)`` 

2637 ''' 

2638 

2639 receiver = Receiver.T( 

2640 help='Receiver location') 

2641 

2642 source_origin = Location.T( 

2643 help='Origin of the source volume grid.') 

2644 

2645 source_depth_min = Float.T( 

2646 help='Minimum source depth [m].') 

2647 

2648 source_depth_max = Float.T( 

2649 help='Maximum source depth [m].') 

2650 

2651 source_depth_delta = Float.T( 

2652 help='Source depth grid spacing [m].') 

2653 

2654 source_east_shift_min = Float.T( 

2655 help='Minimum easting of source grid [m].') 

2656 

2657 source_east_shift_max = Float.T( 

2658 help='Maximum easting of source grid [m].') 

2659 

2660 source_east_shift_delta = Float.T( 

2661 help='Source volume grid spacing in east direction [m].') 

2662 

2663 source_north_shift_min = Float.T( 

2664 help='Minimum northing of source grid [m].') 

2665 

2666 source_north_shift_max = Float.T( 

2667 help='Maximum northing of source grid [m].') 

2668 

2669 source_north_shift_delta = Float.T( 

2670 help='Source volume grid spacing in north direction [m].') 

2671 

2672 short_type = 'C' 

2673 

2674 provided_schemes = ['elastic18'] 

2675 

2676 def get_surface_distance(self, args): 

2677 _, source_east_shift, source_north_shift, _ = args 

2678 sorig = self.source_origin 

2679 sloc = Location( 

2680 lat=sorig.lat, 

2681 lon=sorig.lon, 

2682 north_shift=sorig.north_shift + source_north_shift, 

2683 east_shift=sorig.east_shift + source_east_shift) 

2684 

2685 return self.receiver.distance_to(sloc) 

2686 

2687 def get_distance(self, args): 

2688 sdepth, source_east_shift, source_north_shift, _ = args 

2689 sorig = self.source_origin 

2690 sloc = Location( 

2691 lat=sorig.lat, 

2692 lon=sorig.lon, 

2693 north_shift=sorig.north_shift + source_north_shift, 

2694 east_shift=sorig.east_shift + source_east_shift) 

2695 

2696 return self.receiver.distance_3d_to(sloc) 

2697 

2698 def get_source_depth(self, args): 

2699 return args[0] 

2700 

2701 def get_receiver_depth(self, args): 

2702 return self.receiver.depth 

2703 

2704 def get_source_depths(self): 

2705 return self.coords[0] 

2706 

2707 def _update(self): 

2708 self.mins = num.array([ 

2709 self.source_depth_min, 

2710 self.source_east_shift_min, 

2711 self.source_north_shift_min], 

2712 dtype=float) 

2713 

2714 self.maxs = num.array([ 

2715 self.source_depth_max, 

2716 self.source_east_shift_max, 

2717 self.source_north_shift_max], 

2718 dtype=float) 

2719 

2720 self.deltas = num.array([ 

2721 self.source_depth_delta, 

2722 self.source_east_shift_delta, 

2723 self.source_north_shift_delta], 

2724 dtype=float) 

2725 

2726 self.ns = num.floor((self.maxs - self.mins) / self.deltas + 

2727 vicinity_eps).astype(int) + 1 

2728 self.effective_maxs = self.mins + self.deltas * (self.ns - 1) 

2729 self.deltat = 1.0/self.sample_rate 

2730 self.nrecords = num.prod(self.ns) * self.ncomponents 

2731 

2732 self.coords = tuple( 

2733 num.linspace(mi, ma, n) for (mi, ma, n) in 

2734 zip(self.mins, self.effective_maxs, self.ns)) + \ 

2735 (num.arange(self.ncomponents),) 

2736 

2737 def _make_index_functions(self): 

2738 

2739 amin, bmin, cmin = self.mins 

2740 da, db, dc = self.deltas 

2741 na, nb, nc = self.ns 

2742 ng = self.ncomponents 

2743 

2744 def index_function(a, b, c, ig): 

2745 ia = int(round((a - amin) / da)) 

2746 ib = int(round((b - bmin) / db)) 

2747 ic = int(round((c - cmin) / dc)) 

2748 try: 

2749 return num.ravel_multi_index((ia, ib, ic, ig), 

2750 (na, nb, nc, ng)) 

2751 except ValueError: 

2752 raise OutOfBounds(values=(a, b, c, ig)) 

2753 

2754 def indices_function(a, b, c, ig): 

2755 ia = num.round((a - amin) / da).astype(int) 

2756 ib = num.round((b - bmin) / db).astype(int) 

2757 ic = num.round((c - cmin) / dc).astype(int) 

2758 

2759 try: 

2760 return num.ravel_multi_index((ia, ib, ic, ig), 

2761 (na, nb, nc, ng)) 

2762 except ValueError: 

2763 raise OutOfBounds() 

2764 

2765 def vicinity_function(a, b, c, ig): 

2766 ias = indi12((a - amin) / da, na) 

2767 ibs = indi12((b - bmin) / db, nb) 

2768 ics = indi12((c - cmin) / dc, nc) 

2769 

2770 if not (0 <= ig < ng): 

2771 raise OutOfBounds() 

2772 

2773 indis = [] 

2774 weights = [] 

2775 for ia, va in ias: 

2776 iia = ia*nb*nc*ng 

2777 for ib, vb in ibs: 

2778 iib = ib*nc*ng 

2779 for ic, vc in ics: 

2780 indis.append(iia + iib + ic*ng + ig) 

2781 weights.append(va*vb*vc) 

2782 

2783 return num.array(indis), num.array(weights) 

2784 

2785 def vicinities_function(a, b, c, ig): 

2786 

2787 xa = (a-amin) / da 

2788 xb = (b-bmin) / db 

2789 xc = (c-cmin) / dc 

2790 

2791 xa_fl = num.floor(xa) 

2792 xa_ce = num.ceil(xa) 

2793 xb_fl = num.floor(xb) 

2794 xb_ce = num.ceil(xb) 

2795 xc_fl = num.floor(xc) 

2796 xc_ce = num.ceil(xc) 

2797 va_fl = 1.0 - (xa - xa_fl) 

2798 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl) 

2799 vb_fl = 1.0 - (xb - xb_fl) 

2800 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl) 

2801 vc_fl = 1.0 - (xc - xc_fl) 

2802 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl) 

2803 

2804 ia_fl = xa_fl.astype(int) 

2805 ia_ce = xa_ce.astype(int) 

2806 ib_fl = xb_fl.astype(int) 

2807 ib_ce = xb_ce.astype(int) 

2808 ic_fl = xc_fl.astype(int) 

2809 ic_ce = xc_ce.astype(int) 

2810 

2811 if num.any(ia_fl < 0) or num.any(ia_fl >= na): 

2812 raise OutOfBounds() 

2813 

2814 if num.any(ia_ce < 0) or num.any(ia_ce >= na): 

2815 raise OutOfBounds() 

2816 

2817 if num.any(ib_fl < 0) or num.any(ib_fl >= nb): 

2818 raise OutOfBounds() 

2819 

2820 if num.any(ib_ce < 0) or num.any(ib_ce >= nb): 

2821 raise OutOfBounds() 

2822 

2823 if num.any(ic_fl < 0) or num.any(ic_fl >= nc): 

2824 raise OutOfBounds() 

2825 

2826 if num.any(ic_ce < 0) or num.any(ic_ce >= nc): 

2827 raise OutOfBounds() 

2828 

2829 irecords = num.empty(a.size*8, dtype=int) 

2830 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig 

2831 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig 

2832 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig 

2833 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig 

2834 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig 

2835 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig 

2836 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig 

2837 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig 

2838 

2839 weights = num.empty(a.size*8, dtype=float) 

2840 weights[0::8] = va_fl * vb_fl * vc_fl 

2841 weights[1::8] = va_ce * vb_fl * vc_fl 

2842 weights[2::8] = va_fl * vb_ce * vc_fl 

2843 weights[3::8] = va_ce * vb_ce * vc_fl 

2844 weights[4::8] = va_fl * vb_fl * vc_ce 

2845 weights[5::8] = va_ce * vb_fl * vc_ce 

2846 weights[6::8] = va_fl * vb_ce * vc_ce 

2847 weights[7::8] = va_ce * vb_ce * vc_ce 

2848 

2849 return irecords, weights 

2850 

2851 self._index_function = index_function 

2852 self._indices_function = indices_function 

2853 self._vicinity_function = vicinity_function 

2854 self._vicinities_function = vicinities_function 

2855 

2856 def make_indexing_args(self, source, receiver, icomponents): 

2857 nc = icomponents.size 

2858 

2859 dists = source.distances_to(self.source_origin) 

2860 azis, _ = source.azibazis_to(self.source_origin) 

2861 

2862 source_north_shifts = - num.cos(d2r*azis) * dists 

2863 source_east_shifts = - num.sin(d2r*azis) * dists 

2864 source_depths = source.depths - self.source_origin.depth 

2865 

2866 n = dists.size 

2867 

2868 return (num.tile(source_depths, nc//n), 

2869 num.tile(source_east_shifts, nc//n), 

2870 num.tile(source_north_shifts, nc//n), 

2871 icomponents) 

2872 

2873 def make_indexing_args1(self, source, receiver): 

2874 dist = source.distance_to(self.source_origin) 

2875 azi, _ = source.azibazi_to(self.source_origin) 

2876 

2877 source_north_shift = - num.cos(d2r*azi) * dist 

2878 source_east_shift = - num.sin(d2r*azi) * dist 

2879 source_depth = source.depth - self.source_origin.depth 

2880 

2881 return (source_depth, 

2882 source_east_shift, 

2883 source_north_shift) 

2884 

2885 @property 

2886 def short_extent(self): 

2887 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2888 self.source_depth_min/km, 

2889 self.source_depth_max/km, 

2890 self.source_depth_delta/km, 

2891 self.source_east_shift_min/km, 

2892 self.source_east_shift_max/km, 

2893 self.source_east_shift_delta/km, 

2894 self.source_north_shift_min/km, 

2895 self.source_north_shift_max/km, 

2896 self.source_north_shift_delta/km) 

2897 

2898 

2899class Weighting(Object): 

2900 factor = Float.T(default=1.0) 

2901 

2902 

2903class Taper(Object): 

2904 tmin = Timing.T() 

2905 tmax = Timing.T() 

2906 tfade = Float.T(default=0.0) 

2907 shape = StringChoice.T( 

2908 choices=['cos', 'linear'], 

2909 default='cos', 

2910 optional=True) 

2911 

2912 

2913class SimplePattern(SObject): 

2914 

2915 _pool = {} 

2916 

2917 def __init__(self, pattern): 

2918 self._pattern = pattern 

2919 SObject.__init__(self) 

2920 

2921 def __str__(self): 

2922 return self._pattern 

2923 

2924 @property 

2925 def regex(self): 

2926 pool = SimplePattern._pool 

2927 if self.pattern not in pool: 

2928 rpat = '|'.join(fnmatch.translate(x) for 

2929 x in self.pattern.split('|')) 

2930 pool[self.pattern] = re.compile(rpat, re.I) 

2931 

2932 return pool[self.pattern] 

2933 

2934 def match(self, s): 

2935 return self.regex.match(s) 

2936 

2937 

2938class WaveformType(StringChoice): 

2939 choices = ['dis', 'vel', 'acc', 

2940 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2941 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2942 

2943 

2944class ChannelSelection(Object): 

2945 pattern = SimplePattern.T(optional=True) 

2946 min_sample_rate = Float.T(optional=True) 

2947 max_sample_rate = Float.T(optional=True) 

2948 

2949 

2950class StationSelection(Object): 

2951 includes = SimplePattern.T() 

2952 excludes = SimplePattern.T() 

2953 distance_min = Float.T(optional=True) 

2954 distance_max = Float.T(optional=True) 

2955 azimuth_min = Float.T(optional=True) 

2956 azimuth_max = Float.T(optional=True) 

2957 

2958 

2959class WaveformSelection(Object): 

2960 channel_selection = ChannelSelection.T(optional=True) 

2961 station_selection = StationSelection.T(optional=True) 

2962 taper = Taper.T() 

2963 # filter = FrequencyResponse.T() 

2964 waveform_type = WaveformType.T(default='dis') 

2965 weighting = Weighting.T(optional=True) 

2966 sample_rate = Float.T(optional=True) 

2967 gf_store_id = StringID.T(optional=True) 

2968 

2969 

2970def indi12(x, n): 

2971 ''' 

2972 Get linear interpolation index and weight. 

2973 ''' 

2974 

2975 r = round(x) 

2976 if abs(r - x) < vicinity_eps: 

2977 i = int(r) 

2978 if not (0 <= i < n): 

2979 raise OutOfBounds() 

2980 

2981 return ((int(r), 1.),) 

2982 else: 

2983 f = math.floor(x) 

2984 i = int(f) 

2985 if not (0 <= i < n-1): 

2986 raise OutOfBounds() 

2987 

2988 v = x-f 

2989 return ((i, 1.-v), (i + 1, v)) 

2990 

2991 

2992def float_or_none(s): 

2993 units = { 

2994 'k': 1e3, 

2995 'M': 1e6, 

2996 } 

2997 

2998 factor = 1.0 

2999 if s and s[-1] in units: 

3000 factor = units[s[-1]] 

3001 s = s[:-1] 

3002 if not s: 

3003 raise ValueError("unit without a number: '%s'" % s) 

3004 

3005 if s: 

3006 return float(s) * factor 

3007 else: 

3008 return None 

3009 

3010 

3011class GridSpecError(Exception): 

3012 def __init__(self, s): 

3013 Exception.__init__(self, 'invalid grid specification: %s' % s) 

3014 

3015 

3016def parse_grid_spec(spec): 

3017 try: 

3018 result = [] 

3019 for dspec in spec.split(','): 

3020 t = dspec.split('@') 

3021 num = start = stop = step = None 

3022 if len(t) == 2: 

3023 num = int(t[1]) 

3024 if num <= 0: 

3025 raise GridSpecError(spec) 

3026 

3027 elif len(t) > 2: 

3028 raise GridSpecError(spec) 

3029 

3030 s = t[0] 

3031 v = [float_or_none(x) for x in s.split(':')] 

3032 if len(v) == 1: 

3033 start = stop = v[0] 

3034 if len(v) >= 2: 

3035 start, stop = v[0:2] 

3036 if len(v) == 3: 

3037 step = v[2] 

3038 

3039 if len(v) > 3 or (len(v) > 2 and num is not None): 

3040 raise GridSpecError(spec) 

3041 

3042 if step == 0.0: 

3043 raise GridSpecError(spec) 

3044 

3045 result.append((start, stop, step, num)) 

3046 

3047 except ValueError: 

3048 raise GridSpecError(spec) 

3049 

3050 return result 

3051 

3052 

3053def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5): 

3054 swap = step is not None and step < 0. 

3055 if start is None: 

3056 start = ma if swap else mi 

3057 if stop is None: 

3058 stop = mi if swap else ma 

3059 if step is None: 

3060 step = -inc if ma < mi else inc 

3061 if num is None: 

3062 if (step < 0) != (stop-start < 0): 

3063 raise GridSpecError() 

3064 

3065 num = int(round((stop-start)/step))+1 

3066 stop2 = start + (num-1)*step 

3067 if abs(stop-stop2) > eps: 

3068 num = int(math.floor((stop-start)/step))+1 

3069 stop = start + (num-1)*step 

3070 else: 

3071 stop = stop2 

3072 

3073 if start == stop: 

3074 num = 1 

3075 

3076 return start, stop, num 

3077 

3078 

3079def nditer_outer(x): 

3080 return num.nditer( 

3081 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist()) 

3082 

3083 

3084def nodes(xs): 

3085 ns = [x.size for x in xs] 

3086 nnodes = num.prod(ns) 

3087 ndim = len(xs) 

3088 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype) 

3089 for idim in range(ndim-1, -1, -1): 

3090 x = xs[idim] 

3091 nrepeat = num.prod(ns[idim+1:], dtype=int) 

3092 ntile = num.prod(ns[:idim], dtype=int) 

3093 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat) 

3094 

3095 return nodes 

3096 

3097 

3098def filledi(x, n): 

3099 a = num.empty(n, dtype=int) 

3100 a.fill(x) 

3101 return a 

3102 

3103 

3104config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3105 

3106discretized_source_classes = [ 

3107 DiscretizedExplosionSource, 

3108 DiscretizedSFSource, 

3109 DiscretizedMTSource, 

3110 DiscretizedPorePressureSource] 

3111 

3112 

3113__all__ = ''' 

3114Earthmodel1D 

3115StringID 

3116ScopeType 

3117WaveformType 

3118QuantityType 

3119NearfieldTermsType 

3120Reference 

3121Region 

3122CircularRegion 

3123RectangularRegion 

3124PhaseSelect 

3125InvalidTimingSpecification 

3126Timing 

3127TPDef 

3128OutOfBounds 

3129Location 

3130Receiver 

3131'''.split() + [ 

3132 S.__name__ for S in discretized_source_classes + config_type_classes] + ''' 

3133ComponentScheme 

3134component_scheme_to_description 

3135component_schemes 

3136Config 

3137GridSpecError 

3138Weighting 

3139Taper 

3140SimplePattern 

3141WaveformType 

3142ChannelSelection 

3143StationSelection 

3144WaveformSelection 

3145nditer_outer 

3146dump 

3147load 

3148discretized_source_classes 

3149config_type_classes 

3150UnavailableScheme 

3151InterpolationMethod 

3152SeismosizerTrace 

3153SeismosizerResult 

3154Result 

3155StaticResult 

3156TimingAttributeError 

3157'''.split()