1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import, division 

6 

7import math 

8import re 

9import fnmatch 

10import logging 

11 

12import numpy as num 

13from scipy.interpolate import interp1d 

14 

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

16 StringPattern, Unicode, Float, Bool, Int, TBase, 

17 List, ValidationError, Timestamp, Tuple, Dict) 

18from pyrocko.guts import dump, load # noqa 

19from pyrocko.guts_array import literal, Array 

20from pyrocko.model import Location, gnss 

21 

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

23from pyrocko.util import num_full 

24 

25from .error import StoreError 

26 

27try: 

28 new_str = unicode 

29except NameError: 

30 new_str = str 

31 

32guts_prefix = 'pf' 

33 

34d2r = math.pi / 180. 

35r2d = 1.0 / d2r 

36km = 1000. 

37vicinity_eps = 1e-5 

38 

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

40 

41 

42def fmt_choices(cls): 

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

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

45 

46 

47class InterpolationMethod(StringChoice): 

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

49 

50 

51class SeismosizerTrace(Object): 

52 

53 codes = Tuple.T( 

54 4, String.T(), 

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

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

57 

58 data = Array.T( 

59 shape=(None,), 

60 dtype=num.float32, 

61 serialize_as='base64', 

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

63 help='numpy array with data samples') 

64 

65 deltat = Float.T( 

66 default=1.0, 

67 help='sampling interval [s]') 

68 

69 tmin = Timestamp.T( 

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

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

72 

73 def pyrocko_trace(self): 

74 c = self.codes 

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

76 ydata=self.data, 

77 deltat=self.deltat, 

78 tmin=self.tmin) 

79 

80 @classmethod 

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

82 d = dict( 

83 codes=tr.nslc_id, 

84 tmin=tr.tmin, 

85 deltat=tr.deltat, 

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

87 d.update(kwargs) 

88 return cls(**d) 

89 

90 

91class SeismosizerResult(Object): 

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

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

94 

95 

96class Result(SeismosizerResult): 

97 trace = SeismosizerTrace.T(optional=True) 

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

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

100 

101 

102class StaticResult(SeismosizerResult): 

103 result = Dict.T( 

104 String.T(), 

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

106 

107 

108class GNSSCampaignResult(StaticResult): 

109 campaign = gnss.GNSSCampaign.T( 

110 optional=True) 

111 

112 

113class SatelliteResult(StaticResult): 

114 

115 theta = Array.T( 

116 optional=True, 

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

118 

119 phi = Array.T( 

120 optional=True, 

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

122 

123 

124class KiteSceneResult(SatelliteResult): 

125 

126 shape = Tuple.T() 

127 

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

129 try: 

130 from kite import Scene 

131 except ImportError: 

132 raise ImportError('Kite not installed') 

133 

134 def reshape(arr): 

135 return arr.reshape(self.shape) 

136 

137 displacement = self.result[component] 

138 

139 displacement, theta, phi = map( 

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

141 

142 sc = Scene( 

143 displacement=displacement, 

144 theta=theta, phi=phi, 

145 config=self.config) 

146 

147 return sc 

148 

149 

150class ComponentSchemeDescription(Object): 

151 name = String.T() 

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

153 ncomponents = Int.T() 

154 default_stored_quantity = String.T(optional=True) 

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

156 

157 

158component_scheme_descriptions = [ 

159 ComponentSchemeDescription( 

160 name='elastic2', 

161 source_terms=['m0'], 

162 ncomponents=2, 

163 default_stored_quantity='displacement', 

164 provided_components=[ 

165 'n', 'e', 'd']), 

166 ComponentSchemeDescription( 

167 name='elastic8', 

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

169 ncomponents=8, 

170 default_stored_quantity='displacement', 

171 provided_components=[ 

172 'n', 'e', 'd']), 

173 ComponentSchemeDescription( 

174 name='elastic10', 

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

176 ncomponents=10, 

177 default_stored_quantity='displacement', 

178 provided_components=[ 

179 'n', 'e', 'd']), 

180 ComponentSchemeDescription( 

181 name='elastic18', 

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

183 ncomponents=18, 

184 default_stored_quantity='displacement', 

185 provided_components=[ 

186 'n', 'e', 'd']), 

187 ComponentSchemeDescription( 

188 name='elastic5', 

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

190 ncomponents=5, 

191 default_stored_quantity='displacement', 

192 provided_components=[ 

193 'n', 'e', 'd']), 

194 ComponentSchemeDescription( 

195 name='poroelastic10', 

196 source_terms=['pore_pressure'], 

197 ncomponents=10, 

198 default_stored_quantity=None, 

199 provided_components=[ 

200 'displacement.n', 'displacement.e', 'displacement.d', 

201 'vertical_tilt.n', 'vertical_tilt.e', 

202 'pore_pressure', 

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

204 

205 

206# new names? 

207# 'mt_to_displacement_1d' 

208# 'mt_to_displacement_1d_ff_only' 

209# 'mt_to_gravity_1d' 

210# 'mt_to_stress_1d' 

211# 'explosion_to_displacement_1d' 

212# 'sf_to_displacement_1d' 

213# 'mt_to_displacement_3d' 

214# 'mt_to_gravity_3d' 

215# 'mt_to_stress_3d' 

216# 'pore_pressure_to_displacement_1d' 

217# 'pore_pressure_to_vertical_tilt_1d' 

218# 'pore_pressure_to_pore_pressure_1d' 

219# 'pore_pressure_to_darcy_velocity_1d' 

220 

221 

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

223component_scheme_to_description = dict( 

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

225 

226 

227class ComponentScheme(StringChoice): 

228 ''' 

229 Different Green's Function component schemes are available: 

230 

231 ================= ========================================================= 

232 Name Description 

233 ================= ========================================================= 

234 ``elastic10`` Elastodynamic for 

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

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

237 sources only 

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

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

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

241 MT sources only 

242 ``elastic2`` Elastodynamic for 

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

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

245 isotropic sources only 

246 ``elastic5`` Elastodynamic for 

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

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

249 sources only 

250 ``elastic18`` Elastodynamic for 

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

252 sources only 

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

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

255 ================= ========================================================= 

256 ''' 

257 

258 choices = component_schemes 

259 

260 

261class Earthmodel1D(Object): 

262 dummy_for = cake.LayeredModel 

263 

264 class __T(TBase): 

265 def regularize_extra(self, val): 

266 if isinstance(val, str): 

267 val = cake.LayeredModel.from_scanlines( 

268 cake.read_nd_model_str(val)) 

269 

270 return val 

271 

272 def to_save(self, val): 

273 return literal(cake.write_nd_model_str(val)) 

274 

275 

276class StringID(StringPattern): 

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

278 

279 

280class ScopeType(StringChoice): 

281 choices = [ 

282 'global', 

283 'regional', 

284 'local', 

285 ] 

286 

287 

288class WaveType(StringChoice): 

289 choices = [ 

290 'full waveform', 

291 'bodywave', 

292 'P wave', 

293 'S wave', 

294 'surface wave', 

295 ] 

296 

297 

298class NearfieldTermsType(StringChoice): 

299 choices = [ 

300 'complete', 

301 'incomplete', 

302 'missing', 

303 ] 

304 

305 

306class QuantityType(StringChoice): 

307 choices = [ 

308 'displacement', 

309 'rotation', 

310 'velocity', 

311 'acceleration', 

312 'pressure', 

313 'tilt', 

314 'pore_pressure', 

315 'darcy_velocity', 

316 'vertical_tilt'] 

317 

318 

319class Reference(Object): 

320 id = StringID.T() 

321 type = String.T() 

322 title = Unicode.T() 

323 journal = Unicode.T(optional=True) 

324 volume = Unicode.T(optional=True) 

325 number = Unicode.T(optional=True) 

326 pages = Unicode.T(optional=True) 

327 year = String.T() 

328 note = Unicode.T(optional=True) 

329 issn = String.T(optional=True) 

330 doi = String.T(optional=True) 

331 url = String.T(optional=True) 

332 eprint = String.T(optional=True) 

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

334 publisher = Unicode.T(optional=True) 

335 keywords = Unicode.T(optional=True) 

336 note = Unicode.T(optional=True) 

337 abstract = Unicode.T(optional=True) 

338 

339 @classmethod 

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

341 

342 from pybtex.database.input import bibtex 

343 

344 parser = bibtex.Parser() 

345 

346 if filename is not None: 

347 bib_data = parser.parse_file(filename) 

348 elif stream is not None: 

349 bib_data = parser.parse_stream(stream) 

350 

351 references = [] 

352 

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

354 d = {} 

355 avail = entry.fields.keys() 

356 for prop in cls.T.properties: 

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

358 continue 

359 

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

361 

362 if 'author' in entry.persons: 

363 d['authors'] = [] 

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

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

366 

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

368 references.append(c) 

369 

370 return references 

371 

372 

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

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

375_cake_pat = cake.PhaseDef.allowed_characters_pattern 

376_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic 

377 

378_ppat = r'(' \ 

379 r'cake:' + _cake_pat + \ 

380 r'|iaspei:' + _iaspei_pat + \ 

381 r'|vel_surface:' + _fpat + \ 

382 r'|vel:' + _fpat + \ 

383 r'|stored:' + _spat + \ 

384 r')' 

385 

386 

387timing_regex_old = re.compile( 

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

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

390 

391 

392timing_regex = re.compile( 

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

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

395 

396 

397class PhaseSelect(StringChoice): 

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

399 

400 

401class InvalidTimingSpecification(ValidationError): 

402 pass 

403 

404 

405class Timing(SObject): 

406 ''' 

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

408 

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

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

411 arrival or group of such arrivals. 

412 

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

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

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

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

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

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

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

420 

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

422 

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

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

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

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

427 velocity [km/s] 

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

429 

430 **Examples:** 

431 

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

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

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

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

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

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

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

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

440 undefined for a given geometry 

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

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

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

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

445 arrivals A, B, and C 

446 ''' 

447 

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

449 

450 if s is not None: 

451 offset_is = None 

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

453 try: 

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

455 

456 if s.endswith('S'): 

457 offset_is = 'slowness' 

458 

459 phase_defs = [] 

460 select = '' 

461 

462 except ValueError: 

463 matched = False 

464 m = timing_regex.match(s) 

465 if m: 

466 if m.group(3): 

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

468 elif m.group(19): 

469 phase_defs = [m.group(19)] 

470 else: 

471 phase_defs = [] 

472 

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

474 

475 offset = 0.0 

476 soff = m.group(27) 

477 if soff: 

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

479 if soff.endswith('S'): 

480 offset_is = 'slowness' 

481 elif soff.endswith('%'): 

482 offset_is = 'percent' 

483 

484 matched = True 

485 

486 else: 

487 m = timing_regex_old.match(s) 

488 if m: 

489 if m.group(3): 

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

491 elif m.group(5): 

492 phase_defs = [m.group(5)] 

493 else: 

494 phase_defs = [] 

495 

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

497 

498 offset = 0.0 

499 if m.group(6): 

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

501 

502 phase_defs = [ 

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

504 

505 matched = True 

506 

507 if not matched: 

508 raise InvalidTimingSpecification(s) 

509 

510 kwargs = dict( 

511 phase_defs=phase_defs, 

512 select=select, 

513 offset=offset, 

514 offset_is=offset_is) 

515 

516 SObject.__init__(self, **kwargs) 

517 

518 def __str__(self): 

519 s = [] 

520 if self.phase_defs: 

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

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

523 sphases = '{%s}' % sphases 

524 

525 if self.select: 

526 sphases = self.select + sphases 

527 

528 s.append(sphases) 

529 

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

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

532 if self.offset_is == 'slowness': 

533 s.append('S') 

534 elif self.offset_is == 'percent': 

535 s.append('%') 

536 

537 return ''.join(s) 

538 

539 def evaluate(self, get_phase, args): 

540 try: 

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

542 phase_offset = get_phase( 

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

544 offset = phase_offset(args) 

545 else: 

546 offset = self.offset 

547 

548 if self.phase_defs: 

549 phases = [ 

550 get_phase(phase_def) for phase_def in self.phase_defs] 

551 times = [phase(args) for phase in phases] 

552 if self.offset_is == 'percent': 

553 times = [t*(1.+offset/100.) 

554 for t in times if t is not None] 

555 else: 

556 times = [t+offset for t in times if t is not None] 

557 

558 if not times: 

559 return None 

560 elif self.select == 'first': 

561 return min(times) 

562 elif self.select == 'last': 

563 return max(times) 

564 else: 

565 return times[0] 

566 else: 

567 return offset 

568 

569 except spit.OutOfBounds: 

570 raise OutOfBounds(args) 

571 

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

573 offset = Float.T(default=0.0) 

574 offset_is = String.T(optional=True) 

575 select = PhaseSelect.T( 

576 default='', 

577 help=('Can be either ``\'%s\'``, ``\'%s\'``, or ``\'%s\'``. ' % 

578 tuple(PhaseSelect.choices))) 

579 

580 

581def mkdefs(s): 

582 defs = [] 

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

584 try: 

585 defs.append(float(x)) 

586 except ValueError: 

587 if x.startswith('!'): 

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

589 else: 

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

591 

592 return defs 

593 

594 

595class TPDef(Object): 

596 ''' 

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

598 ''' 

599 

600 id = StringID.T( 

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

602 definition = String.T( 

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

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

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

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

607 

608 @property 

609 def phases(self): 

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

611 if isinstance(x, cake.PhaseDef)] 

612 

613 @property 

614 def horizontal_velocities(self): 

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

616 

617 

618class OutOfBounds(Exception): 

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

620 Exception.__init__(self) 

621 self.values = values 

622 self.reason = reason 

623 self.context = None 

624 

625 def __str__(self): 

626 scontext = '' 

627 if self.context: 

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

629 

630 if self.reason: 

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

632 

633 if self.values: 

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

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

636 else: 

637 return 'out of bounds%s' % scontext 

638 

639 

640class MultiLocation(Object): 

641 

642 lats = Array.T( 

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

644 help='Latitudes of targets.') 

645 lons = Array.T( 

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

647 help='Longitude of targets.') 

648 north_shifts = Array.T( 

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

650 help='North shifts of targets.') 

651 east_shifts = Array.T( 

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

653 help='East shifts of targets.') 

654 elevation = Array.T( 

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

656 help='Elevations of targets.') 

657 

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

659 self._coords5 = None 

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

661 

662 @property 

663 def coords5(self): 

664 if self._coords5 is not None: 

665 return self._coords5 

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

667 self.elevation] 

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

669 if not sizes: 

670 raise AttributeError('Empty StaticTarget') 

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

672 raise AttributeError('Inconsistent coordinate sizes.') 

673 

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

675 for idx, p in enumerate(props): 

676 if p is not None: 

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

678 

679 return self._coords5 

680 

681 @property 

682 def ncoords(self): 

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

684 return 0 

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

686 

687 def get_latlon(self): 

688 ''' 

689 Get all coordinates as lat/lon. 

690 

691 :returns: Coordinates in Latitude, Longitude 

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

693 ''' 

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

695 for i in range(self.ncoords): 

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

697 return latlons 

698 

699 def get_corner_coords(self): 

700 ''' 

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

702 

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

704 :rtype: tuple 

705 ''' 

706 latlon = self.get_latlon() 

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

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

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

710 

711 

712class Receiver(Location): 

713 codes = Tuple.T( 

714 3, String.T(), 

715 optional=True, 

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

717 

718 def pyrocko_station(self): 

719 from pyrocko import model 

720 

721 lat, lon = self.effective_latlon 

722 

723 return model.Station( 

724 network=self.codes[0], 

725 station=self.codes[1], 

726 location=self.codes[2], 

727 lat=lat, 

728 lon=lon, 

729 depth=self.depth) 

730 

731 

732def g(x, d): 

733 if x is None: 

734 return d 

735 else: 

736 return x 

737 

738 

739class UnavailableScheme(Exception): 

740 pass 

741 

742 

743class InvalidNComponents(Exception): 

744 pass 

745 

746 

747class DiscretizedSource(Object): 

748 ''' 

749 Base class for discretized sources. 

750 

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

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

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

754 source distributions are represented by subclasses of the 

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

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

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

758 explosion/implosion type source distributions. The geometry calculations 

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

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

761 

762 Like in the :py:class:`Location` class, the positions of the point sources 

763 contained in the discretized source are defined by a common reference point 

764 (:py:attr:`lat`, :py:attr:`lon`) and cartesian offsets to this 

765 (:py:attr:`north_shifts`, :py:attr:`east_shifts`, :py:attr:`depths`). 

766 Alternatively latitude and longitude of each contained point source can be 

767 specified directly (:py:attr:`lats`, :py:attr:`lons`). 

768 ''' 

769 

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

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

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

773 lat = Float.T(optional=True) 

774 lon = Float.T(optional=True) 

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

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

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

778 dl = Float.T(optional=True) 

779 dw = Float.T(optional=True) 

780 nl = Float.T(optional=True) 

781 nw = Float.T(optional=True) 

782 

783 @classmethod 

784 def check_scheme(cls, scheme): 

785 ''' 

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

787 

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

789 supported by this discretized source class. 

790 ''' 

791 

792 if scheme not in cls.provided_schemes: 

793 raise UnavailableScheme( 

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

795 (cls.__name__, scheme)) 

796 

797 def __init__(self, **kwargs): 

798 Object.__init__(self, **kwargs) 

799 self._latlons = None 

800 

801 def __setattr__(self, name, value): 

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

803 'lats', 'lons'): 

804 self.__dict__['_latlons'] = None 

805 

806 Object.__setattr__(self, name, value) 

807 

808 def get_source_terms(self, scheme): 

809 raise NotImplementedError() 

810 

811 def make_weights(self, receiver, scheme): 

812 raise NotImplementedError() 

813 

814 @property 

815 def effective_latlons(self): 

816 ''' 

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

818 ''' 

819 

820 if self._latlons is None: 

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

822 if (self.north_shifts is not None and 

823 self.east_shifts is not None): 

824 self._latlons = orthodrome.ne_to_latlon( 

825 self.lats, self.lons, 

826 self.north_shifts, self.east_shifts) 

827 else: 

828 self._latlons = self.lats, self.lons 

829 else: 

830 lat = g(self.lat, 0.0) 

831 lon = g(self.lon, 0.0) 

832 self._latlons = orthodrome.ne_to_latlon( 

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

834 

835 return self._latlons 

836 

837 @property 

838 def effective_north_shifts(self): 

839 if self.north_shifts is not None: 

840 return self.north_shifts 

841 else: 

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

843 

844 @property 

845 def effective_east_shifts(self): 

846 if self.east_shifts is not None: 

847 return self.east_shifts 

848 else: 

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

850 

851 def same_origin(self, receiver): 

852 ''' 

853 Check if receiver has same reference point. 

854 ''' 

855 

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

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

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

859 

860 def azibazis_to(self, receiver): 

861 ''' 

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

863 points. 

864 ''' 

865 

866 if self.same_origin(receiver): 

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

868 receiver.north_shift - self.north_shifts) 

869 bazis = azis + 180. 

870 else: 

871 slats, slons = self.effective_latlons 

872 rlat, rlon = receiver.effective_latlon 

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

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

875 

876 return azis, bazis 

877 

878 def distances_to(self, receiver): 

879 ''' 

880 Compute distances to receiver for all contained points. 

881 ''' 

882 if self.same_origin(receiver): 

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

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

885 

886 else: 

887 slats, slons = self.effective_latlons 

888 rlat, rlon = receiver.effective_latlon 

889 return orthodrome.distance_accurate50m_numpy(slats, slons, 

890 rlat, rlon) 

891 

892 def element_coords(self, i): 

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

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

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

896 else: 

897 lat = self.lat 

898 lon = self.lon 

899 

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

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

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

903 

904 else: 

905 north_shift = east_shift = 0.0 

906 

907 return lat, lon, north_shift, east_shift 

908 

909 def coords5(self): 

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

911 

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

913 xs[:, 0] = self.lats 

914 xs[:, 1] = self.lons 

915 else: 

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

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

918 

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

920 xs[:, 2] = self.north_shifts 

921 xs[:, 3] = self.east_shifts 

922 

923 xs[:, 4] = self.depths 

924 

925 return xs 

926 

927 @property 

928 def nelements(self): 

929 return self.times.size 

930 

931 @classmethod 

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

933 ''' 

934 Combine several discretized source models. 

935 

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

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

938 factors and reference times of the parameterized (undiscretized) 

939 sources match or are accounted for. 

940 ''' 

941 

942 first = sources[0] 

943 

944 if not all(type(s) == type(first) for s in sources): 

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

946 'sources of same type.') 

947 

948 latlons = [] 

949 for s in sources: 

950 latlons.append(s.effective_latlons) 

951 

952 lats, lons = num.hstack(latlons) 

953 

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

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

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

957 same_ref = num.all( 

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

959 else: 

960 same_ref = False 

961 

962 cat = num.concatenate 

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

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

965 

966 if same_ref: 

967 lat = first.lat 

968 lon = first.lon 

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

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

971 lats = None 

972 lons = None 

973 else: 

974 lat = None 

975 lon = None 

976 north_shifts = None 

977 east_shifts = None 

978 

979 return cls( 

980 times=times, lat=lat, lon=lon, lats=lats, lons=lons, 

981 north_shifts=north_shifts, east_shifts=east_shifts, 

982 depths=depths, **kwargs) 

983 

984 def centroid_position(self): 

985 moments = self.moments() 

986 norm = num.sum(moments) 

987 if norm != 0.0: 

988 w = moments / num.sum(moments) 

989 else: 

990 w = num.ones(moments.size) 

991 

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

993 lats, lons = self.effective_latlons 

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

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

996 else: 

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

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

999 

1000 cn = num.sum(n*w) 

1001 ce = num.sum(e*w) 

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

1003 

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

1005 lat = clat 

1006 lon = clon 

1007 north_shift = 0. 

1008 east_shift = 0. 

1009 else: 

1010 lat = g(self.lat, 0.0) 

1011 lon = g(self.lon, 0.0) 

1012 north_shift = cn 

1013 east_shift = ce 

1014 

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

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

1017 return tuple(float(x) for x in 

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

1019 

1020 

1021class DiscretizedExplosionSource(DiscretizedSource): 

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

1023 

1024 provided_schemes = ( 

1025 'elastic2', 

1026 'elastic8', 

1027 'elastic10', 

1028 ) 

1029 

1030 def get_source_terms(self, scheme): 

1031 self.check_scheme(scheme) 

1032 

1033 if scheme == 'elastic2': 

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

1035 

1036 elif scheme in ('elastic8', 'elastic10'): 

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

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

1039 return m6s 

1040 else: 

1041 assert False 

1042 

1043 def make_weights(self, receiver, scheme): 

1044 self.check_scheme(scheme) 

1045 

1046 azis, bazis = self.azibazis_to(receiver) 

1047 

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

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

1050 

1051 m0s = self.m0s 

1052 n = azis.size 

1053 

1054 cat = num.concatenate 

1055 rep = num.repeat 

1056 

1057 if scheme == 'elastic2': 

1058 w_n = cb*m0s 

1059 g_n = filledi(0, n) 

1060 w_e = sb*m0s 

1061 g_e = filledi(0, n) 

1062 w_d = m0s 

1063 g_d = filledi(1, n) 

1064 

1065 elif scheme == 'elastic8': 

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

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

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

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

1070 w_d = cat((m0s, m0s)) 

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

1072 

1073 elif scheme == 'elastic10': 

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

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

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

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

1078 w_d = cat((m0s, m0s, m0s)) 

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

1080 

1081 else: 

1082 assert False 

1083 

1084 return ( 

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

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

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

1088 ) 

1089 

1090 def split(self): 

1091 from pyrocko.gf.seismosizer import ExplosionSource 

1092 sources = [] 

1093 for i in range(self.nelements): 

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

1095 sources.append(ExplosionSource( 

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

1097 lat=lat, 

1098 lon=lon, 

1099 north_shift=north_shift, 

1100 east_shift=east_shift, 

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

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

1103 

1104 return sources 

1105 

1106 def moments(self): 

1107 return self.m0s 

1108 

1109 def centroid(self): 

1110 from pyrocko.gf.seismosizer import ExplosionSource 

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

1112 self.centroid_position() 

1113 

1114 return ExplosionSource( 

1115 time=time, 

1116 lat=lat, 

1117 lon=lon, 

1118 north_shift=north_shift, 

1119 east_shift=east_shift, 

1120 depth=depth, 

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

1122 

1123 @classmethod 

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

1125 ''' 

1126 Combine several discretized source models. 

1127 

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

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

1130 factors and reference times of the parameterized (undiscretized) 

1131 sources match or are accounted for. 

1132 ''' 

1133 

1134 if 'm0s' not in kwargs: 

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

1136 

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

1138 **kwargs) 

1139 

1140 

1141class DiscretizedSFSource(DiscretizedSource): 

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

1143 

1144 provided_schemes = ( 

1145 'elastic5', 

1146 ) 

1147 

1148 def get_source_terms(self, scheme): 

1149 self.check_scheme(scheme) 

1150 

1151 return self.forces 

1152 

1153 def make_weights(self, receiver, scheme): 

1154 self.check_scheme(scheme) 

1155 

1156 azis, bazis = self.azibazis_to(receiver) 

1157 

1158 sa = num.sin(azis*d2r) 

1159 ca = num.cos(azis*d2r) 

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

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

1162 

1163 forces = self.forces 

1164 fn = forces[:, 0] 

1165 fe = forces[:, 1] 

1166 fd = forces[:, 2] 

1167 

1168 f0 = fd 

1169 f1 = ca * fn + sa * fe 

1170 f2 = ca * fe - sa * fn 

1171 

1172 n = azis.size 

1173 

1174 if scheme == 'elastic5': 

1175 ioff = 0 

1176 

1177 cat = num.concatenate 

1178 rep = num.repeat 

1179 

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

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

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

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

1184 w_d = cat((f0, f1)) 

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

1186 

1187 return ( 

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

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

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

1191 ) 

1192 

1193 @classmethod 

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

1195 ''' 

1196 Combine several discretized source models. 

1197 

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

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

1200 factors and reference times of the parameterized (undiscretized) 

1201 sources match or are accounted for. 

1202 ''' 

1203 

1204 if 'forces' not in kwargs: 

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

1206 

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

1208 

1209 def moments(self): 

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

1211 

1212 def centroid(self): 

1213 from pyrocko.gf.seismosizer import SFSource 

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

1215 self.centroid_position() 

1216 

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

1218 return SFSource( 

1219 time=time, 

1220 lat=lat, 

1221 lon=lon, 

1222 north_shift=north_shift, 

1223 east_shift=east_shift, 

1224 depth=depth, 

1225 fn=fn, 

1226 fe=fe, 

1227 fd=fd) 

1228 

1229 

1230class DiscretizedMTSource(DiscretizedSource): 

1231 m6s = Array.T( 

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

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

1234 

1235 provided_schemes = ( 

1236 'elastic8', 

1237 'elastic10', 

1238 'elastic18', 

1239 ) 

1240 

1241 def get_source_terms(self, scheme): 

1242 self.check_scheme(scheme) 

1243 return self.m6s 

1244 

1245 def make_weights(self, receiver, scheme): 

1246 self.check_scheme(scheme) 

1247 

1248 m6s = self.m6s 

1249 n = m6s.shape[0] 

1250 rep = num.repeat 

1251 

1252 if scheme == 'elastic18': 

1253 w_n = m6s.flatten() 

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

1255 w_e = m6s.flatten() 

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

1257 w_d = m6s.flatten() 

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

1259 

1260 else: 

1261 azis, bazis = self.azibazis_to(receiver) 

1262 

1263 sa = num.sin(azis*d2r) 

1264 ca = num.cos(azis*d2r) 

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

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

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

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

1269 

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

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

1272 f2 = m6s[:, 2] 

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

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

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

1276 

1277 cat = num.concatenate 

1278 

1279 if scheme == 'elastic8': 

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

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

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

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

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

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

1286 

1287 elif scheme == 'elastic10': 

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

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

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

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

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

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

1294 

1295 else: 

1296 assert False 

1297 

1298 return ( 

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

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

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

1302 ) 

1303 

1304 def split(self): 

1305 from pyrocko.gf.seismosizer import MTSource 

1306 sources = [] 

1307 for i in range(self.nelements): 

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

1309 sources.append(MTSource( 

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

1311 lat=lat, 

1312 lon=lon, 

1313 north_shift=north_shift, 

1314 east_shift=east_shift, 

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

1316 m6=self.m6s[i])) 

1317 

1318 return sources 

1319 

1320 def moments(self): 

1321 moments = num.array( 

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

1323 for m6 in self.m6s]) 

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

1325 

1326 def get_moment_rate(self, deltat=None): 

1327 moments = self.moments() 

1328 times = self.times 

1329 times -= times.min() 

1330 

1331 t_max = times.max() 

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

1333 mom_times[mom_times > t_max] = t_max 

1334 

1335 # Right open histrogram bins 

1336 mom, _ = num.histogram( 

1337 times, 

1338 bins=mom_times, 

1339 weights=moments) 

1340 

1341 deltat = num.diff(mom_times) 

1342 mom_rate = mom / deltat 

1343 

1344 return mom_rate, mom_times[1:] 

1345 

1346 def centroid(self): 

1347 from pyrocko.gf.seismosizer import MTSource 

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

1349 self.centroid_position() 

1350 

1351 return MTSource( 

1352 time=time, 

1353 lat=lat, 

1354 lon=lon, 

1355 north_shift=north_shift, 

1356 east_shift=east_shift, 

1357 depth=depth, 

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

1359 

1360 @classmethod 

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

1362 ''' 

1363 Combine several discretized source models. 

1364 

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

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

1367 factors and reference times of the parameterized (undiscretized) 

1368 sources match or are accounted for. 

1369 ''' 

1370 

1371 if 'm6s' not in kwargs: 

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

1373 

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

1375 

1376 

1377class DiscretizedPorePressureSource(DiscretizedSource): 

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

1379 

1380 provided_schemes = ( 

1381 'poroelastic10', 

1382 ) 

1383 

1384 def get_source_terms(self, scheme): 

1385 self.check_scheme(scheme) 

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

1387 

1388 def make_weights(self, receiver, scheme): 

1389 self.check_scheme(scheme) 

1390 

1391 azis, bazis = self.azibazis_to(receiver) 

1392 

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

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

1395 

1396 pp = self.pp 

1397 n = bazis.size 

1398 

1399 w_un = cb*pp 

1400 g_un = filledi(1, n) 

1401 w_ue = sb*pp 

1402 g_ue = filledi(1, n) 

1403 w_ud = pp 

1404 g_ud = filledi(0, n) 

1405 

1406 w_tn = cb*pp 

1407 g_tn = filledi(6, n) 

1408 w_te = sb*pp 

1409 g_te = filledi(6, n) 

1410 

1411 w_pp = pp 

1412 g_pp = filledi(7, n) 

1413 

1414 w_dvn = cb*pp 

1415 g_dvn = filledi(9, n) 

1416 w_dve = sb*pp 

1417 g_dve = filledi(9, n) 

1418 w_dvd = pp 

1419 g_dvd = filledi(8, n) 

1420 

1421 return ( 

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

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

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

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

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

1427 ('pore_pressure', w_pp, g_pp), 

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

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

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

1431 ) 

1432 

1433 def moments(self): 

1434 return self.pp 

1435 

1436 def centroid(self): 

1437 from pyrocko.gf.seismosizer import PorePressurePointSource 

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

1439 self.centroid_position() 

1440 

1441 return PorePressurePointSource( 

1442 time=time, 

1443 lat=lat, 

1444 lon=lon, 

1445 north_shift=north_shift, 

1446 east_shift=east_shift, 

1447 depth=depth, 

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

1449 

1450 @classmethod 

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

1452 ''' 

1453 Combine several discretized source models. 

1454 

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

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

1457 factors and reference times of the parameterized (undiscretized) 

1458 sources match or are accounted for. 

1459 ''' 

1460 

1461 if 'pp' not in kwargs: 

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

1463 

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

1465 **kwargs) 

1466 

1467 

1468class Region(Object): 

1469 name = String.T(optional=True) 

1470 

1471 

1472class RectangularRegion(Region): 

1473 lat_min = Float.T() 

1474 lat_max = Float.T() 

1475 lon_min = Float.T() 

1476 lon_max = Float.T() 

1477 

1478 

1479class CircularRegion(Region): 

1480 lat = Float.T() 

1481 lon = Float.T() 

1482 radius = Float.T() 

1483 

1484 

1485class Config(Object): 

1486 ''' 

1487 Green's function store meta information. 

1488 

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

1490 configuration types are: 

1491 

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

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

1494 

1495 * Problem is invariant to horizontal translations and rotations around 

1496 vertical axis. 

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

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

1499 component)`` 

1500 

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

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

1503 

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

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

1506 receiver_depth, component)`` 

1507 

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

1509 constraints but fixed receiver positions 

1510 

1511 * Cartesian source volume around a reference point 

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

1513 source_east_shift, source_north_shift, component)`` 

1514 ''' 

1515 

1516 id = StringID.T( 

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

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

1519 'letter.') 

1520 

1521 derived_from_id = StringID.T( 

1522 optional=True, 

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

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

1525 

1526 version = String.T( 

1527 default='1.0', 

1528 optional=True, 

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

1530 

1531 modelling_code_id = StringID.T( 

1532 optional=True, 

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

1534 

1535 author = Unicode.T( 

1536 optional=True, 

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

1538 

1539 author_email = String.T( 

1540 optional=True, 

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

1542 

1543 created_time = Timestamp.T( 

1544 optional=True, 

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

1546 

1547 regions = List.T( 

1548 Region.T(), 

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

1550 

1551 scope_type = ScopeType.T( 

1552 optional=True, 

1553 help='Distance range scope of the store ' 

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

1555 

1556 waveform_type = WaveType.T( 

1557 optional=True, 

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

1559 

1560 nearfield_terms = NearfieldTermsType.T( 

1561 optional=True, 

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

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

1564 

1565 description = String.T( 

1566 optional=True, 

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

1568 

1569 references = List.T( 

1570 Reference.T(), 

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

1572 'related work.') 

1573 

1574 earthmodel_1d = Earthmodel1D.T( 

1575 optional=True, 

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

1577 

1578 earthmodel_receiver_1d = Earthmodel1D.T( 

1579 optional=True, 

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

1581 

1582 can_interpolate_source = Bool.T( 

1583 optional=True, 

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

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

1586 

1587 can_interpolate_receiver = Bool.T( 

1588 optional=True, 

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

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

1591 

1592 frequency_min = Float.T( 

1593 optional=True, 

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

1595 

1596 frequency_max = Float.T( 

1597 optional=True, 

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

1599 

1600 sample_rate = Float.T( 

1601 optional=True, 

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

1603 

1604 factor = Float.T( 

1605 default=1.0, 

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

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

1608 optional=True) 

1609 

1610 component_scheme = ComponentScheme.T( 

1611 default='elastic10', 

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

1613 

1614 stored_quantity = QuantityType.T( 

1615 optional=True, 

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

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

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

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

1620 

1621 tabulated_phases = List.T( 

1622 TPDef.T(), 

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

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

1625 

1626 ncomponents = Int.T( 

1627 optional=True, 

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

1629 

1630 uuid = String.T( 

1631 optional=True, 

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

1633 'GF store for practical purposes.') 

1634 

1635 reference = String.T( 

1636 optional=True, 

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

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

1639 

1640 def __init__(self, **kwargs): 

1641 self._do_auto_updates = False 

1642 Object.__init__(self, **kwargs) 

1643 self._index_function = None 

1644 self._indices_function = None 

1645 self._vicinity_function = None 

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

1647 self._do_auto_updates = True 

1648 self.update() 

1649 

1650 def check_ncomponents(self): 

1651 ncomponents = component_scheme_to_description[ 

1652 self.component_scheme].ncomponents 

1653 

1654 if self.ncomponents is None: 

1655 self.ncomponents = ncomponents 

1656 elif ncomponents != self.ncomponents: 

1657 raise InvalidNComponents( 

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

1659 self.ncomponents, self.component_scheme)) 

1660 

1661 def __setattr__(self, name, value): 

1662 Object.__setattr__(self, name, value) 

1663 try: 

1664 self.T.get_property(name) 

1665 if self._do_auto_updates: 

1666 self.update() 

1667 

1668 except ValueError: 

1669 pass 

1670 

1671 def update(self): 

1672 self.check_ncomponents() 

1673 self._update() 

1674 self._make_index_functions() 

1675 

1676 def irecord(self, *args): 

1677 return self._index_function(*args) 

1678 

1679 def irecords(self, *args): 

1680 return self._indices_function(*args) 

1681 

1682 def vicinity(self, *args): 

1683 return self._vicinity_function(*args) 

1684 

1685 def vicinities(self, *args): 

1686 return self._vicinities_function(*args) 

1687 

1688 def grid_interpolation_coefficients(self, *args): 

1689 return self._grid_interpolation_coefficients(*args) 

1690 

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

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

1693 

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

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

1696 

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

1698 i = 0 

1699 arrs = [] 

1700 ntotal = 1 

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

1702 if gdef and len(gdef) > i: 

1703 sssn = gdef[i] 

1704 else: 

1705 sssn = (None,)*4 

1706 

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

1708 ntotal *= len(arr) 

1709 

1710 arrs.append(arr) 

1711 i += 1 

1712 

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

1714 return nditer_outer(arrs[:level]) 

1715 

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

1717 nthreads=0): 

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

1719 

1720 out = [] 

1721 delays = source.times 

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

1723 receiver, 

1724 self.component_scheme): 

1725 

1726 weights *= self.factor 

1727 

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

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

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

1731 

1732 return out 

1733 

1734 def short_info(self): 

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

1736 

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

1738 interpolation=None): 

1739 ''' 

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

1741 

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

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

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

1745 ``(lat, lon)`` 

1746 :param interpolation: interpolation method. Choose from 

1747 ``('nearest_neighbor', 'multilinear')`` 

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

1749 point 

1750 

1751 The default implementation retrieves and interpolates the shear moduli 

1752 from the contained 1D velocity profile. 

1753 ''' 

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

1755 parameter='shear_moduli', 

1756 interpolation=interpolation) 

1757 

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

1759 interpolation=None): 

1760 ''' 

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

1762 

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

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

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

1766 ``(lat, lon)`` 

1767 :param interpolation: interpolation method. Choose from 

1768 ``('nearest_neighbor', 'multilinear')`` 

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

1770 point 

1771 

1772 The default implementation retrieves and interpolates the lambda moduli 

1773 from the contained 1D velocity profile. 

1774 ''' 

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

1776 parameter='lambda_moduli', 

1777 interpolation=interpolation) 

1778 

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

1780 interpolation=None): 

1781 ''' 

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

1783 

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

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

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

1787 ``(lat, lon)`` 

1788 :param interpolation: interpolation method. Choose from 

1789 ``('nearest_neighbor', 'multilinear')`` 

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

1791 point 

1792 

1793 The default implementation retrieves and interpolates the lambda moduli 

1794 from the contained 1D velocity profile. 

1795 ''' 

1796 lambda_moduli = self.get_material_property( 

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

1798 interpolation=interpolation) 

1799 shear_moduli = self.get_material_property( 

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

1801 interpolation=interpolation) 

1802 return lambda_moduli + (2 / 3) * shear_moduli 

1803 

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

1805 ''' 

1806 Get Vs at given points from contained velocity model. 

1807 

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

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

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

1811 ``(lat, lon)`` 

1812 :param interpolation: interpolation method. Choose from 

1813 ``('nearest_neighbor', 'multilinear')`` 

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

1815 point 

1816 

1817 The default implementation retrieves and interpolates Vs 

1818 from the contained 1D velocity profile. 

1819 ''' 

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

1821 parameter='vs', 

1822 interpolation=interpolation) 

1823 

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

1825 ''' 

1826 Get Vp at given points from contained velocity model. 

1827 

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

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

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

1831 ``(lat, lon)`` 

1832 :param interpolation: interpolation method. Choose from 

1833 ``('nearest_neighbor', 'multilinear')`` 

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

1835 point 

1836 

1837 The default implementation retrieves and interpolates Vp 

1838 from the contained 1D velocity profile. 

1839 ''' 

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

1841 parameter='vp', 

1842 interpolation=interpolation) 

1843 

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

1845 ''' 

1846 Get rho at given points from contained velocity model. 

1847 

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

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

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

1851 ``(lat, lon)`` 

1852 :param interpolation: interpolation method. Choose from 

1853 ``('nearest_neighbor', 'multilinear')`` 

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

1855 point 

1856 

1857 The default implementation retrieves and interpolates rho 

1858 from the contained 1D velocity profile. 

1859 ''' 

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

1861 parameter='rho', 

1862 interpolation=interpolation) 

1863 

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

1865 interpolation=None): 

1866 

1867 if interpolation is None: 

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

1869 "multilinear", "nearest_neighbor") 

1870 

1871 earthmod = self.earthmodel_1d 

1872 store_depth_profile = self.get_source_depths() 

1873 z_profile = earthmod.profile('z') 

1874 

1875 if parameter == 'vs': 

1876 vs_profile = earthmod.profile('vs') 

1877 profile = num.interp( 

1878 store_depth_profile, z_profile, vs_profile) 

1879 

1880 elif parameter == 'vp': 

1881 vp_profile = earthmod.profile('vp') 

1882 profile = num.interp( 

1883 store_depth_profile, z_profile, vp_profile) 

1884 

1885 elif parameter == 'rho': 

1886 rho_profile = earthmod.profile('rho') 

1887 

1888 profile = num.interp( 

1889 store_depth_profile, z_profile, rho_profile) 

1890 

1891 elif parameter == 'shear_moduli': 

1892 vs_profile = earthmod.profile('vs') 

1893 rho_profile = earthmod.profile('rho') 

1894 

1895 store_vs_profile = num.interp( 

1896 store_depth_profile, z_profile, vs_profile) 

1897 store_rho_profile = num.interp( 

1898 store_depth_profile, z_profile, rho_profile) 

1899 

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

1901 

1902 elif parameter == 'lambda_moduli': 

1903 vs_profile = earthmod.profile('vs') 

1904 vp_profile = earthmod.profile('vp') 

1905 rho_profile = earthmod.profile('rho') 

1906 

1907 store_vs_profile = num.interp( 

1908 store_depth_profile, z_profile, vs_profile) 

1909 store_vp_profile = num.interp( 

1910 store_depth_profile, z_profile, vp_profile) 

1911 store_rho_profile = num.interp( 

1912 store_depth_profile, z_profile, rho_profile) 

1913 

1914 profile = store_rho_profile * ( 

1915 num.power(store_vp_profile, 2) - 

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

1917 else: 

1918 raise TypeError( 

1919 'parameter %s not available' % parameter) 

1920 

1921 if interpolation == 'multilinear': 

1922 kind = 'linear' 

1923 elif interpolation == 'nearest_neighbor': 

1924 kind = 'nearest' 

1925 else: 

1926 raise TypeError( 

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

1928 

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

1930 

1931 try: 

1932 return interpolator(points[:, 2]) 

1933 except ValueError: 

1934 raise OutOfBounds() 

1935 

1936 def is_static(self): 

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

1938 if self.modelling_code_id.startswith(code): 

1939 return True 

1940 return False 

1941 

1942 def is_dynamic(self): 

1943 return not self.is_static() 

1944 

1945 def get_source_depths(self): 

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

1947 

1948 def get_tabulated_phase(self, phase_id): 

1949 ''' 

1950 Get tabulated phase definition. 

1951 ''' 

1952 

1953 for pdef in self.tabulated_phases: 

1954 if pdef.id == phase_id: 

1955 return pdef 

1956 

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

1958 

1959 def fix_ttt_holes(self, sptree, mode): 

1960 raise StoreError( 

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

1962 % self.short_type) 

1963 

1964 

1965class ConfigTypeA(Config): 

1966 ''' 

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

1968 

1969 * Problem is invariant to horizontal translations and rotations around 

1970 vertical axis. 

1971 

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

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

1974 component)`` 

1975 

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

1977 points. 

1978 ''' 

1979 

1980 receiver_depth = Float.T( 

1981 default=0.0, 

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

1983 

1984 source_depth_min = Float.T( 

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

1986 

1987 source_depth_max = Float.T( 

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

1989 

1990 source_depth_delta = Float.T( 

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

1992 

1993 distance_min = Float.T( 

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

1995 

1996 distance_max = Float.T( 

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

1998 

1999 distance_delta = Float.T( 

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

2001 

2002 short_type = 'A' 

2003 

2004 provided_schemes = [ 

2005 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10'] 

2006 

2007 def get_surface_distance(self, args): 

2008 return args[1] 

2009 

2010 def get_distance(self, args): 

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

2012 

2013 def get_source_depth(self, args): 

2014 return args[0] 

2015 

2016 def get_source_depths(self): 

2017 return self.coords[0] 

2018 

2019 def get_receiver_depth(self, args): 

2020 return self.receiver_depth 

2021 

2022 def _update(self): 

2023 self.mins = num.array( 

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

2025 self.maxs = num.array( 

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

2027 self.deltas = num.array( 

2028 [self.source_depth_delta, self.distance_delta], 

2029 dtype=float) 

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

2031 vicinity_eps).astype(int) + 1 

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

2033 self.deltat = 1.0/self.sample_rate 

2034 self.nrecords = num.product(self.ns) * self.ncomponents 

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

2036 (mi, ma, n) in 

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

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

2039 

2040 self.nsource_depths, self.ndistances = self.ns 

2041 

2042 def _make_index_functions(self): 

2043 

2044 amin, bmin = self.mins 

2045 da, db = self.deltas 

2046 na, nb = self.ns 

2047 

2048 ng = self.ncomponents 

2049 

2050 def index_function(a, b, ig): 

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

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

2053 try: 

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

2055 except ValueError: 

2056 raise OutOfBounds() 

2057 

2058 def indices_function(a, b, ig): 

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

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

2061 try: 

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

2063 except ValueError: 

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

2065 try: 

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

2067 except ValueError: 

2068 raise OutOfBounds() 

2069 

2070 def grid_interpolation_coefficients(a, b): 

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

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

2073 return ias, ibs 

2074 

2075 def vicinity_function(a, b, ig): 

2076 ias, ibs = grid_interpolation_coefficients(a, b) 

2077 

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

2079 raise OutOfBounds() 

2080 

2081 indis = [] 

2082 weights = [] 

2083 for ia, va in ias: 

2084 iia = ia*nb*ng 

2085 for ib, vb in ibs: 

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

2087 weights.append(va*vb) 

2088 

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

2090 

2091 def vicinities_function(a, b, ig): 

2092 

2093 xa = (a - amin) / da 

2094 xb = (b - bmin) / db 

2095 

2096 xa_fl = num.floor(xa) 

2097 xa_ce = num.ceil(xa) 

2098 xb_fl = num.floor(xb) 

2099 xb_ce = num.ceil(xb) 

2100 va_fl = 1.0 - (xa - xa_fl) 

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

2102 vb_fl = 1.0 - (xb - xb_fl) 

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

2104 

2105 ia_fl = xa_fl.astype(int) 

2106 ia_ce = xa_ce.astype(int) 

2107 ib_fl = xb_fl.astype(int) 

2108 ib_ce = xb_ce.astype(int) 

2109 

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

2111 raise OutOfBounds() 

2112 

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

2114 raise OutOfBounds() 

2115 

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

2117 raise OutOfBounds() 

2118 

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

2120 raise OutOfBounds() 

2121 

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

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

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

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

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

2127 

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

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

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

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

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

2133 

2134 return irecords, weights 

2135 

2136 self._index_function = index_function 

2137 self._indices_function = indices_function 

2138 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2139 self._vicinity_function = vicinity_function 

2140 self._vicinities_function = vicinities_function 

2141 

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

2143 nc = icomponents.size 

2144 dists = source.distances_to(receiver) 

2145 n = dists.size 

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

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

2148 icomponents) 

2149 

2150 def make_indexing_args1(self, source, receiver): 

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

2152 

2153 @property 

2154 def short_extent(self): 

2155 return '%g:%g:%g x %g:%g:%g' % ( 

2156 self.source_depth_min/km, 

2157 self.source_depth_max/km, 

2158 self.source_depth_delta/km, 

2159 self.distance_min/km, 

2160 self.distance_max/km, 

2161 self.distance_delta/km) 

2162 

2163 def fix_ttt_holes(self, sptree, mode): 

2164 from pyrocko import eikonal_ext, spit 

2165 

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

2167 

2168 delta = self.deltas[-1] 

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

2170 

2171 nsources, ndistances = self.ns 

2172 

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

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

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

2176 

2177 speeds = self.get_material_property( 

2178 0., 0., points, 

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

2180 interpolation='multilinear') 

2181 

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

2183 

2184 times = sptree.interpolate_many(nodes) 

2185 

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

2187 times = times.reshape(speeds.shape) 

2188 

2189 try: 

2190 eikonal_ext.eikonal_solver_fmm_cartesian( 

2191 speeds, times, delta) 

2192 except eikonal_ext.EikonalExtError as e: 

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

2194 logger.debug( 

2195 'Got a warning from eikonal solver ' 

2196 '- may be ok...') 

2197 else: 

2198 raise 

2199 

2200 def func(x): 

2201 ibs, ics = \ 

2202 self.grid_interpolation_coefficients(*x) 

2203 

2204 t = 0 

2205 for ib, vb in ibs: 

2206 for ic, vc in ics: 

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

2208 

2209 return t 

2210 

2211 return spit.SPTree( 

2212 f=func, 

2213 ftol=sptree.ftol, 

2214 xbounds=sptree.xbounds, 

2215 xtols=sptree.xtols) 

2216 

2217 

2218class ConfigTypeB(Config): 

2219 ''' 

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

2221 

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

2223 receiver depth 

2224 

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

2226 receiver_distance, component)`` 

2227 ''' 

2228 

2229 receiver_depth_min = Float.T( 

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

2231 

2232 receiver_depth_max = Float.T( 

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

2234 

2235 receiver_depth_delta = Float.T( 

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

2237 

2238 source_depth_min = Float.T( 

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

2240 

2241 source_depth_max = Float.T( 

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

2243 

2244 source_depth_delta = Float.T( 

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

2246 

2247 distance_min = Float.T( 

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

2249 

2250 distance_max = Float.T( 

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

2252 

2253 distance_delta = Float.T( 

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

2255 

2256 short_type = 'B' 

2257 

2258 provided_schemes = [ 

2259 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10'] 

2260 

2261 def get_distance(self, args): 

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

2263 

2264 def get_surface_distance(self, args): 

2265 return args[2] 

2266 

2267 def get_source_depth(self, args): 

2268 return args[1] 

2269 

2270 def get_receiver_depth(self, args): 

2271 return args[0] 

2272 

2273 def get_source_depths(self): 

2274 return self.coords[1] 

2275 

2276 def _update(self): 

2277 self.mins = num.array([ 

2278 self.receiver_depth_min, 

2279 self.source_depth_min, 

2280 self.distance_min], 

2281 dtype=float) 

2282 

2283 self.maxs = num.array([ 

2284 self.receiver_depth_max, 

2285 self.source_depth_max, 

2286 self.distance_max], 

2287 dtype=float) 

2288 

2289 self.deltas = num.array([ 

2290 self.receiver_depth_delta, 

2291 self.source_depth_delta, 

2292 self.distance_delta], 

2293 dtype=float) 

2294 

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

2296 vicinity_eps).astype(int) + 1 

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

2298 self.deltat = 1.0/self.sample_rate 

2299 self.nrecords = num.product(self.ns) * self.ncomponents 

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

2301 (mi, ma, n) in 

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

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

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

2305 

2306 def _make_index_functions(self): 

2307 

2308 amin, bmin, cmin = self.mins 

2309 da, db, dc = self.deltas 

2310 na, nb, nc = self.ns 

2311 ng = self.ncomponents 

2312 

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

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

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

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

2317 try: 

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

2319 (na, nb, nc, ng)) 

2320 except ValueError: 

2321 raise OutOfBounds() 

2322 

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

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

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

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

2327 try: 

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

2329 (na, nb, nc, ng)) 

2330 except ValueError: 

2331 raise OutOfBounds() 

2332 

2333 def grid_interpolation_coefficients(a, b, c): 

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

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

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

2337 return ias, ibs, ics 

2338 

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

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

2341 

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

2343 raise OutOfBounds() 

2344 

2345 indis = [] 

2346 weights = [] 

2347 for ia, va in ias: 

2348 iia = ia*nb*nc*ng 

2349 for ib, vb in ibs: 

2350 iib = ib*nc*ng 

2351 for ic, vc in ics: 

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

2353 weights.append(va*vb*vc) 

2354 

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

2356 

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

2358 

2359 xa = (a - amin) / da 

2360 xb = (b - bmin) / db 

2361 xc = (c - cmin) / dc 

2362 

2363 xa_fl = num.floor(xa) 

2364 xa_ce = num.ceil(xa) 

2365 xb_fl = num.floor(xb) 

2366 xb_ce = num.ceil(xb) 

2367 xc_fl = num.floor(xc) 

2368 xc_ce = num.ceil(xc) 

2369 va_fl = 1.0 - (xa - xa_fl) 

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

2371 vb_fl = 1.0 - (xb - xb_fl) 

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

2373 vc_fl = 1.0 - (xc - xc_fl) 

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

2375 

2376 ia_fl = xa_fl.astype(int) 

2377 ia_ce = xa_ce.astype(int) 

2378 ib_fl = xb_fl.astype(int) 

2379 ib_ce = xb_ce.astype(int) 

2380 ic_fl = xc_fl.astype(int) 

2381 ic_ce = xc_ce.astype(int) 

2382 

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

2384 raise OutOfBounds() 

2385 

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

2387 raise OutOfBounds() 

2388 

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

2390 raise OutOfBounds() 

2391 

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

2393 raise OutOfBounds() 

2394 

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

2396 raise OutOfBounds() 

2397 

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

2399 raise OutOfBounds() 

2400 

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

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

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

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

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

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

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

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

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

2410 

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

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

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

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

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

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

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

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

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

2420 

2421 return irecords, weights 

2422 

2423 self._index_function = index_function 

2424 self._indices_function = indices_function 

2425 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2426 self._vicinity_function = vicinity_function 

2427 self._vicinities_function = vicinities_function 

2428 

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

2430 nc = icomponents.size 

2431 dists = source.distances_to(receiver) 

2432 n = dists.size 

2433 receiver_depths = num.empty(nc) 

2434 receiver_depths.fill(receiver.depth) 

2435 return (receiver_depths, 

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

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

2438 icomponents) 

2439 

2440 def make_indexing_args1(self, source, receiver): 

2441 return (receiver.depth, 

2442 source.depth, 

2443 source.distance_to(receiver)) 

2444 

2445 @property 

2446 def short_extent(self): 

2447 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2448 self.receiver_depth_min/km, 

2449 self.receiver_depth_max/km, 

2450 self.receiver_depth_delta/km, 

2451 self.source_depth_min/km, 

2452 self.source_depth_max/km, 

2453 self.source_depth_delta/km, 

2454 self.distance_min/km, 

2455 self.distance_max/km, 

2456 self.distance_delta/km) 

2457 

2458 def fix_ttt_holes(self, sptree, mode): 

2459 from pyrocko import eikonal_ext, spit 

2460 

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

2462 

2463 delta = self.deltas[-1] 

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

2465 

2466 nreceivers, nsources, ndistances = self.ns 

2467 

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

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

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

2471 

2472 speeds = self.get_material_property( 

2473 0., 0., points, 

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

2475 interpolation='multilinear') 

2476 

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

2478 

2479 receiver_times = [] 

2480 for ireceiver in range(nreceivers): 

2481 nodes = num.hstack([ 

2482 num_full( 

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

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

2485 nodes_sr]) 

2486 

2487 times = sptree.interpolate_many(nodes) 

2488 

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

2490 

2491 times = times.reshape(speeds.shape) 

2492 

2493 try: 

2494 eikonal_ext.eikonal_solver_fmm_cartesian( 

2495 speeds, times, delta) 

2496 except eikonal_ext.EikonalExtError as e: 

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

2498 logger.debug( 

2499 'Got a warning from eikonal solver ' 

2500 '- may be ok...') 

2501 else: 

2502 raise 

2503 

2504 receiver_times.append(times) 

2505 

2506 def func(x): 

2507 ias, ibs, ics = \ 

2508 self.grid_interpolation_coefficients(*x) 

2509 

2510 t = 0 

2511 for ia, va in ias: 

2512 times = receiver_times[ia] 

2513 for ib, vb in ibs: 

2514 for ic, vc in ics: 

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

2516 

2517 return t 

2518 

2519 return spit.SPTree( 

2520 f=func, 

2521 ftol=sptree.ftol, 

2522 xbounds=sptree.xbounds, 

2523 xtols=sptree.xtols) 

2524 

2525 

2526class ConfigTypeC(Config): 

2527 ''' 

2528 No symmetrical constraints but fixed receiver positions. 

2529 

2530 * Cartesian 3D source volume around a reference point 

2531 

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

2533 source_east_shift, source_north_shift, component)`` 

2534 ''' 

2535 

2536 receivers = List.T( 

2537 Receiver.T(), 

2538 help='List of fixed receivers.') 

2539 

2540 source_origin = Location.T( 

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

2542 

2543 source_depth_min = Float.T( 

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

2545 

2546 source_depth_max = Float.T( 

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

2548 

2549 source_depth_delta = Float.T( 

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

2551 

2552 source_east_shift_min = Float.T( 

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

2554 

2555 source_east_shift_max = Float.T( 

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

2557 

2558 source_east_shift_delta = Float.T( 

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

2560 

2561 source_north_shift_min = Float.T( 

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

2563 

2564 source_north_shift_max = Float.T( 

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

2566 

2567 source_north_shift_delta = Float.T( 

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

2569 

2570 short_type = 'C' 

2571 

2572 provided_schemes = ['elastic18'] 

2573 

2574 def get_surface_distance(self, args): 

2575 ireceiver, _, source_east_shift, source_north_shift, _ = args 

2576 sorig = self.source_origin 

2577 sloc = Location( 

2578 lat=sorig.lat, 

2579 lon=sorig.lon, 

2580 north_shift=sorig.north_shift + source_north_shift, 

2581 east_shift=sorig.east_shift + source_east_shift) 

2582 

2583 return self.receivers[args[0]].distance_to(sloc) 

2584 

2585 def get_distance(self, args): 

2586 # to be improved... 

2587 ireceiver, sdepth, source_east_shift, source_north_shift, _ = args 

2588 sorig = self.source_origin 

2589 sloc = Location( 

2590 lat=sorig.lat, 

2591 lon=sorig.lon, 

2592 north_shift=sorig.north_shift + source_north_shift, 

2593 east_shift=sorig.east_shift + source_east_shift) 

2594 

2595 return math.sqrt( 

2596 self.receivers[args[0]].distance_to(sloc)**2 + sdepth**2) 

2597 

2598 def get_source_depth(self, args): 

2599 return args[1] 

2600 

2601 def get_receiver_depth(self, args): 

2602 return self.receivers[args[0]].depth 

2603 

2604 def get_source_depths(self): 

2605 return self.coords[0] 

2606 

2607 def _update(self): 

2608 self.mins = num.array([ 

2609 self.source_depth_min, 

2610 self.source_east_shift_min, 

2611 self.source_north_shift_min], 

2612 dtype=float) 

2613 

2614 self.maxs = num.array([ 

2615 self.source_depth_max, 

2616 self.source_east_shift_max, 

2617 self.source_north_shift_max], 

2618 dtype=float) 

2619 

2620 self.deltas = num.array([ 

2621 self.source_depth_delta, 

2622 self.source_east_shift_delta, 

2623 self.source_north_shift_delta], 

2624 dtype=float) 

2625 

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

2627 vicinity_eps).astype(int) + 1 

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

2629 self.deltat = 1.0/self.sample_rate 

2630 self.nreceivers = len(self.receivers) 

2631 self.nrecords = \ 

2632 self.nreceivers * num.product(self.ns) * self.ncomponents 

2633 

2634 self.coords = (num.arange(self.nreceivers),) + \ 

2635 tuple(num.linspace(mi, ma, n) for (mi, ma, n) in 

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

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

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

2639 

2640 self._distances_cache = {} 

2641 

2642 def _make_index_functions(self): 

2643 

2644 amin, bmin, cmin = self.mins 

2645 da, db, dc = self.deltas 

2646 na, nb, nc = self.ns 

2647 ng = self.ncomponents 

2648 nr = self.nreceivers 

2649 

2650 def index_function(ir, a, b, c, ig): 

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

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

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

2654 try: 

2655 return num.ravel_multi_index((ir, ia, ib, ic, ig), 

2656 (nr, na, nb, nc, ng)) 

2657 except ValueError: 

2658 raise OutOfBounds() 

2659 

2660 def indices_function(ir, a, b, c, ig): 

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

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

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

2664 

2665 try: 

2666 return num.ravel_multi_index((ir, ia, ib, ic, ig), 

2667 (nr, na, nb, nc, ng)) 

2668 except ValueError: 

2669 raise OutOfBounds() 

2670 

2671 def vicinity_function(ir, a, b, c, ig): 

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

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

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

2675 

2676 if not (0 <= ir < nr): 

2677 raise OutOfBounds() 

2678 

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

2680 raise OutOfBounds() 

2681 

2682 indis = [] 

2683 weights = [] 

2684 iir = ir*na*nb*nc*ng 

2685 for ia, va in ias: 

2686 iia = ia*nb*nc*ng 

2687 for ib, vb in ibs: 

2688 iib = ib*nc*ng 

2689 for ic, vc in ics: 

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

2691 weights.append(va*vb*vc) 

2692 

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

2694 

2695 def vicinities_function(ir, a, b, c, ig): 

2696 

2697 xa = (a-amin) / da 

2698 xb = (b-bmin) / db 

2699 xc = (c-cmin) / dc 

2700 

2701 xa_fl = num.floor(xa) 

2702 xa_ce = num.ceil(xa) 

2703 xb_fl = num.floor(xb) 

2704 xb_ce = num.ceil(xb) 

2705 xc_fl = num.floor(xc) 

2706 xc_ce = num.ceil(xc) 

2707 va_fl = 1.0 - (xa - xa_fl) 

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

2709 vb_fl = 1.0 - (xb - xb_fl) 

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

2711 vc_fl = 1.0 - (xc - xc_fl) 

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

2713 

2714 ia_fl = xa_fl.astype(int) 

2715 ia_ce = xa_ce.astype(int) 

2716 ib_fl = xb_fl.astype(int) 

2717 ib_ce = xb_ce.astype(int) 

2718 ic_fl = xc_fl.astype(int) 

2719 ic_ce = xc_ce.astype(int) 

2720 

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

2722 raise OutOfBounds() 

2723 

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

2725 raise OutOfBounds() 

2726 

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

2728 raise OutOfBounds() 

2729 

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

2731 raise OutOfBounds() 

2732 

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

2734 raise OutOfBounds() 

2735 

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

2737 raise OutOfBounds() 

2738 

2739 irig = ir*na*nb*nc*ng + ig 

2740 

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

2742 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + irig 

2743 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + irig 

2744 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + irig 

2745 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + irig 

2746 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + irig 

2747 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + irig 

2748 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + irig 

2749 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + irig 

2750 

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

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

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

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

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

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

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

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

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

2760 

2761 return irecords, weights 

2762 

2763 self._index_function = index_function 

2764 self._indices_function = indices_function 

2765 self._vicinity_function = vicinity_function 

2766 self._vicinities_function = vicinities_function 

2767 

2768 def lookup_ireceiver(self, receiver): 

2769 k = (receiver.lat, receiver.lon, 

2770 receiver.north_shift, receiver.east_shift) 

2771 dh = min(self.source_north_shift_delta, self.source_east_shift_delta) 

2772 dv = self.source_depth_delta 

2773 

2774 for irec, rec in enumerate(self.receivers): 

2775 if (k, irec) not in self._distances_cache: 

2776 self._distances_cache[k, irec] = math.sqrt( 

2777 (receiver.distance_to(rec)/dh)**2 + 

2778 ((rec.depth - receiver.depth)/dv)**2) 

2779 

2780 if self._distances_cache[k, irec] < 0.1: 

2781 return irec 

2782 

2783 raise OutOfBounds( 

2784 reason='No GFs available for receiver at (%g, %g).' % 

2785 receiver.effective_latlon) 

2786 

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

2788 nc = icomponents.size 

2789 

2790 dists = source.distances_to(self.source_origin) 

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

2792 

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

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

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

2796 

2797 n = dists.size 

2798 ireceivers = num.empty(nc, dtype=int) 

2799 ireceivers.fill(self.lookup_ireceiver(receiver)) 

2800 

2801 return (ireceivers, 

2802 num.tile(source_depths, nc//n), 

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

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

2805 icomponents) 

2806 

2807 def make_indexing_args1(self, source, receiver): 

2808 dist = source.distance_to(self.source_origin) 

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

2810 

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

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

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

2814 

2815 return (self.lookup_ireceiver(receiver), 

2816 source_depth, 

2817 source_east_shift, 

2818 source_north_shift) 

2819 

2820 @property 

2821 def short_extent(self): 

2822 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2823 self.source_depth_min/km, 

2824 self.source_depth_max/km, 

2825 self.source_depth_delta/km, 

2826 self.source_east_shift_min/km, 

2827 self.source_east_shift_max/km, 

2828 self.source_east_shift_delta/km, 

2829 self.source_north_shift_min/km, 

2830 self.source_north_shift_max/km, 

2831 self.source_north_shift_delta/km) 

2832 

2833 

2834class Weighting(Object): 

2835 factor = Float.T(default=1.0) 

2836 

2837 

2838class Taper(Object): 

2839 tmin = Timing.T() 

2840 tmax = Timing.T() 

2841 tfade = Float.T(default=0.0) 

2842 shape = StringChoice.T( 

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

2844 default='cos', 

2845 optional=True) 

2846 

2847 

2848class SimplePattern(SObject): 

2849 

2850 _pool = {} 

2851 

2852 def __init__(self, pattern): 

2853 self._pattern = pattern 

2854 SObject.__init__(self) 

2855 

2856 def __str__(self): 

2857 return self._pattern 

2858 

2859 @property 

2860 def regex(self): 

2861 pool = SimplePattern._pool 

2862 if self.pattern not in pool: 

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

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

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

2866 

2867 return pool[self.pattern] 

2868 

2869 def match(self, s): 

2870 return self.regex.match(s) 

2871 

2872 

2873class WaveformType(StringChoice): 

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

2875 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2876 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2877 

2878 

2879class ChannelSelection(Object): 

2880 pattern = SimplePattern.T(optional=True) 

2881 min_sample_rate = Float.T(optional=True) 

2882 max_sample_rate = Float.T(optional=True) 

2883 

2884 

2885class StationSelection(Object): 

2886 includes = SimplePattern.T() 

2887 excludes = SimplePattern.T() 

2888 distance_min = Float.T(optional=True) 

2889 distance_max = Float.T(optional=True) 

2890 azimuth_min = Float.T(optional=True) 

2891 azimuth_max = Float.T(optional=True) 

2892 

2893 

2894class WaveformSelection(Object): 

2895 channel_selection = ChannelSelection.T(optional=True) 

2896 station_selection = StationSelection.T(optional=True) 

2897 taper = Taper.T() 

2898 # filter = FrequencyResponse.T() 

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

2900 weighting = Weighting.T(optional=True) 

2901 sample_rate = Float.T(optional=True) 

2902 gf_store_id = StringID.T(optional=True) 

2903 

2904 

2905def indi12(x, n): 

2906 ''' 

2907 Get linear interpolation index and weight. 

2908 ''' 

2909 

2910 r = round(x) 

2911 if abs(r - x) < vicinity_eps: 

2912 i = int(r) 

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

2914 raise OutOfBounds() 

2915 

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

2917 else: 

2918 f = math.floor(x) 

2919 i = int(f) 

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

2921 raise OutOfBounds() 

2922 

2923 v = x-f 

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

2925 

2926 

2927def float_or_none(s): 

2928 units = { 

2929 'k': 1e3, 

2930 'M': 1e6, 

2931 } 

2932 

2933 factor = 1.0 

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

2935 factor = units[s[-1]] 

2936 s = s[:-1] 

2937 if not s: 

2938 raise ValueError('unit without a number: \'%s\'' % s) 

2939 

2940 if s: 

2941 return float(s) * factor 

2942 else: 

2943 return None 

2944 

2945 

2946class GridSpecError(Exception): 

2947 def __init__(self, s): 

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

2949 

2950 

2951def parse_grid_spec(spec): 

2952 try: 

2953 result = [] 

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

2955 t = dspec.split('@') 

2956 num = start = stop = step = None 

2957 if len(t) == 2: 

2958 num = int(t[1]) 

2959 if num <= 0: 

2960 raise GridSpecError(spec) 

2961 

2962 elif len(t) > 2: 

2963 raise GridSpecError(spec) 

2964 

2965 s = t[0] 

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

2967 if len(v) == 1: 

2968 start = stop = v[0] 

2969 if len(v) >= 2: 

2970 start, stop = v[0:2] 

2971 if len(v) == 3: 

2972 step = v[2] 

2973 

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

2975 raise GridSpecError(spec) 

2976 

2977 if step == 0.0: 

2978 raise GridSpecError(spec) 

2979 

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

2981 

2982 except ValueError: 

2983 raise GridSpecError(spec) 

2984 

2985 return result 

2986 

2987 

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

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

2990 if start is None: 

2991 start = ma if swap else mi 

2992 if stop is None: 

2993 stop = mi if swap else ma 

2994 if step is None: 

2995 step = -inc if ma < mi else inc 

2996 if num is None: 

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

2998 raise GridSpecError() 

2999 

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

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

3002 if abs(stop-stop2) > eps: 

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

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

3005 else: 

3006 stop = stop2 

3007 

3008 if start == stop: 

3009 num = 1 

3010 

3011 return start, stop, num 

3012 

3013 

3014def nditer_outer(x): 

3015 return num.nditer( 

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

3017 

3018 

3019def nodes(xs): 

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

3021 nnodes = num.prod(ns) 

3022 ndim = len(xs) 

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

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

3025 x = xs[idim] 

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

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

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

3029 

3030 return nodes 

3031 

3032 

3033def filledi(x, n): 

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

3035 a.fill(x) 

3036 return a 

3037 

3038 

3039config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3040 

3041discretized_source_classes = [ 

3042 DiscretizedExplosionSource, 

3043 DiscretizedSFSource, 

3044 DiscretizedMTSource, 

3045 DiscretizedPorePressureSource] 

3046 

3047 

3048__all__ = ''' 

3049Earthmodel1D 

3050StringID 

3051ScopeType 

3052WaveformType 

3053QuantityType 

3054NearfieldTermsType 

3055Reference 

3056Region 

3057CircularRegion 

3058RectangularRegion 

3059PhaseSelect 

3060InvalidTimingSpecification 

3061Timing 

3062TPDef 

3063OutOfBounds 

3064Location 

3065Receiver 

3066'''.split() + [ 

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

3068ComponentScheme 

3069component_scheme_to_description 

3070component_schemes 

3071Config 

3072GridSpecError 

3073Weighting 

3074Taper 

3075SimplePattern 

3076WaveformType 

3077ChannelSelection 

3078StationSelection 

3079WaveformSelection 

3080nditer_outer 

3081dump 

3082load 

3083discretized_source_classes 

3084config_type_classes 

3085UnavailableScheme 

3086InterpolationMethod 

3087SeismosizerTrace 

3088SeismosizerResult 

3089Result 

3090StaticResult 

3091'''.split()