1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import math 

7import re 

8import fnmatch 

9import logging 

10 

11import numpy as num 

12from scipy.interpolate import interp1d 

13 

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

15 StringPattern, Unicode, Float, Bool, Int, TBase, 

16 List, ValidationError, Timestamp, Tuple, Dict) 

17from pyrocko.guts import dump, load # noqa 

18from pyrocko.guts_array import literal, Array 

19from pyrocko.model import Location, gnss 

20 

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

22from pyrocko.util import num_full 

23 

24from .error import StoreError 

25 

26try: 

27 new_str = unicode 

28except NameError: 

29 new_str = str 

30 

31guts_prefix = 'pf' 

32 

33d2r = math.pi / 180. 

34r2d = 1.0 / d2r 

35km = 1000. 

36vicinity_eps = 1e-5 

37 

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

39 

40 

41def fmt_choices(cls): 

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

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

44 

45 

46class InterpolationMethod(StringChoice): 

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

48 

49 

50class SeismosizerTrace(Object): 

51 

52 codes = Tuple.T( 

53 4, String.T(), 

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

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

56 

57 data = Array.T( 

58 shape=(None,), 

59 dtype=num.float32, 

60 serialize_as='base64', 

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

62 help='numpy array with data samples') 

63 

64 deltat = Float.T( 

65 default=1.0, 

66 help='sampling interval [s]') 

67 

68 tmin = Timestamp.T( 

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

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

71 

72 def pyrocko_trace(self): 

73 c = self.codes 

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

75 ydata=self.data, 

76 deltat=self.deltat, 

77 tmin=self.tmin) 

78 

79 @classmethod 

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

81 d = dict( 

82 codes=tr.nslc_id, 

83 tmin=tr.tmin, 

84 deltat=tr.deltat, 

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

86 d.update(kwargs) 

87 return cls(**d) 

88 

89 

90class SeismosizerResult(Object): 

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

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

93 

94 

95class Result(SeismosizerResult): 

96 trace = SeismosizerTrace.T(optional=True) 

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

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

99 

100 

101class StaticResult(SeismosizerResult): 

102 result = Dict.T( 

103 String.T(), 

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

105 

106 

107class GNSSCampaignResult(StaticResult): 

108 campaign = gnss.GNSSCampaign.T( 

109 optional=True) 

110 

111 

112class SatelliteResult(StaticResult): 

113 

114 theta = Array.T( 

115 optional=True, 

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

117 

118 phi = Array.T( 

119 optional=True, 

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

121 

122 

123class KiteSceneResult(SatelliteResult): 

124 

125 shape = Tuple.T() 

126 

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

128 try: 

129 from kite import Scene 

130 except ImportError: 

131 raise ImportError('Kite not installed') 

132 

133 def reshape(arr): 

134 return arr.reshape(self.shape) 

135 

136 displacement = self.result[component] 

137 

138 displacement, theta, phi = map( 

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

140 

141 sc = Scene( 

142 displacement=displacement, 

143 theta=theta, phi=phi, 

144 config=self.config) 

145 

146 return sc 

147 

148 

149class ComponentSchemeDescription(Object): 

150 name = String.T() 

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

152 ncomponents = Int.T() 

153 default_stored_quantity = String.T(optional=True) 

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

155 

156 

157component_scheme_descriptions = [ 

158 ComponentSchemeDescription( 

159 name='elastic2', 

160 source_terms=['m0'], 

161 ncomponents=2, 

162 default_stored_quantity='displacement', 

163 provided_components=[ 

164 'n', 'e', 'd']), 

165 ComponentSchemeDescription( 

166 name='elastic8', 

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

168 ncomponents=8, 

169 default_stored_quantity='displacement', 

170 provided_components=[ 

171 'n', 'e', 'd']), 

172 ComponentSchemeDescription( 

173 name='elastic10', 

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

175 ncomponents=10, 

176 default_stored_quantity='displacement', 

177 provided_components=[ 

178 'n', 'e', 'd']), 

179 ComponentSchemeDescription( 

180 name='elastic18', 

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

182 ncomponents=18, 

183 default_stored_quantity='displacement', 

184 provided_components=[ 

185 'n', 'e', 'd']), 

186 ComponentSchemeDescription( 

187 name='elastic5', 

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

189 ncomponents=5, 

190 default_stored_quantity='displacement', 

191 provided_components=[ 

192 'n', 'e', 'd']), 

193 ComponentSchemeDescription( 

194 name='poroelastic10', 

195 source_terms=['pore_pressure'], 

196 ncomponents=10, 

197 default_stored_quantity=None, 

198 provided_components=[ 

199 'displacement.n', 'displacement.e', 'displacement.d', 

200 'vertical_tilt.n', 'vertical_tilt.e', 

201 'pore_pressure', 

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

203 

204 

205# new names? 

206# 'mt_to_displacement_1d' 

207# 'mt_to_displacement_1d_ff_only' 

208# 'mt_to_gravity_1d' 

209# 'mt_to_stress_1d' 

210# 'explosion_to_displacement_1d' 

211# 'sf_to_displacement_1d' 

212# 'mt_to_displacement_3d' 

213# 'mt_to_gravity_3d' 

214# 'mt_to_stress_3d' 

215# 'pore_pressure_to_displacement_1d' 

216# 'pore_pressure_to_vertical_tilt_1d' 

217# 'pore_pressure_to_pore_pressure_1d' 

218# 'pore_pressure_to_darcy_velocity_1d' 

219 

220 

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

222component_scheme_to_description = dict( 

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

224 

225 

226class ComponentScheme(StringChoice): 

227 ''' 

228 Different Green's Function component schemes are available: 

229 

230 ================= ========================================================= 

231 Name Description 

232 ================= ========================================================= 

233 ``elastic10`` Elastodynamic for 

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

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

236 sources only 

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

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

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

240 MT sources only 

241 ``elastic2`` Elastodynamic for 

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

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

244 isotropic sources only 

245 ``elastic5`` Elastodynamic for 

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

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

248 sources only 

249 ``elastic18`` Elastodynamic for 

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

251 sources only 

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

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

254 ================= ========================================================= 

255 ''' 

256 

257 choices = component_schemes 

258 

259 

260class Earthmodel1D(Object): 

261 dummy_for = cake.LayeredModel 

262 

263 class __T(TBase): 

264 def regularize_extra(self, val): 

265 if isinstance(val, str): 

266 val = cake.LayeredModel.from_scanlines( 

267 cake.read_nd_model_str(val)) 

268 

269 return val 

270 

271 def to_save(self, val): 

272 return literal(cake.write_nd_model_str(val)) 

273 

274 

275class StringID(StringPattern): 

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

277 

278 

279class ScopeType(StringChoice): 

280 choices = [ 

281 'global', 

282 'regional', 

283 'local', 

284 ] 

285 

286 

287class WaveType(StringChoice): 

288 choices = [ 

289 'full waveform', 

290 'bodywave', 

291 'P wave', 

292 'S wave', 

293 'surface wave', 

294 ] 

295 

296 

297class NearfieldTermsType(StringChoice): 

298 choices = [ 

299 'complete', 

300 'incomplete', 

301 'missing', 

302 ] 

303 

304 

305class QuantityType(StringChoice): 

306 choices = [ 

307 'displacement', 

308 'rotation', 

309 'velocity', 

310 'acceleration', 

311 'pressure', 

312 'tilt', 

313 'pore_pressure', 

314 'darcy_velocity', 

315 'vertical_tilt'] 

316 

317 

318class Reference(Object): 

319 id = StringID.T() 

320 type = String.T() 

321 title = Unicode.T() 

322 journal = Unicode.T(optional=True) 

323 volume = Unicode.T(optional=True) 

324 number = Unicode.T(optional=True) 

325 pages = Unicode.T(optional=True) 

326 year = String.T() 

327 note = Unicode.T(optional=True) 

328 issn = String.T(optional=True) 

329 doi = String.T(optional=True) 

330 url = String.T(optional=True) 

331 eprint = String.T(optional=True) 

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

333 publisher = Unicode.T(optional=True) 

334 keywords = Unicode.T(optional=True) 

335 note = Unicode.T(optional=True) 

336 abstract = Unicode.T(optional=True) 

337 

338 @classmethod 

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

340 

341 from pybtex.database.input import bibtex 

342 

343 parser = bibtex.Parser() 

344 

345 if filename is not None: 

346 bib_data = parser.parse_file(filename) 

347 elif stream is not None: 

348 bib_data = parser.parse_stream(stream) 

349 

350 references = [] 

351 

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

353 d = {} 

354 avail = entry.fields.keys() 

355 for prop in cls.T.properties: 

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

357 continue 

358 

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

360 

361 if 'author' in entry.persons: 

362 d['authors'] = [] 

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

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

365 

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

367 references.append(c) 

368 

369 return references 

370 

371 

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

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

374_cake_pat = cake.PhaseDef.allowed_characters_pattern 

375_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic 

376 

377_ppat = r'(' \ 

378 r'cake:' + _cake_pat + \ 

379 r'|iaspei:' + _iaspei_pat + \ 

380 r'|vel_surface:' + _fpat + \ 

381 r'|vel:' + _fpat + \ 

382 r'|stored:' + _spat + \ 

383 r')' 

384 

385 

386timing_regex_old = re.compile( 

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

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

389 

390 

391timing_regex = re.compile( 

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

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

394 

395 

396class PhaseSelect(StringChoice): 

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

398 

399 

400class InvalidTimingSpecification(ValidationError): 

401 pass 

402 

403 

404class Timing(SObject): 

405 ''' 

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

407 

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

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

410 arrival or group of such arrivals. 

411 

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

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

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

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

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

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

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

419 

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

421 

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

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

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

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

426 velocity [km/s] 

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

428 

429 **Examples:** 

430 

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

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

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

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

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

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

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

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

439 undefined for a given geometry 

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

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

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

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

444 arrivals A, B, and C 

445 ''' 

446 

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

448 

449 if s is not None: 

450 offset_is = None 

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

452 try: 

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

454 

455 if s.endswith('S'): 

456 offset_is = 'slowness' 

457 

458 phase_defs = [] 

459 select = '' 

460 

461 except ValueError: 

462 matched = False 

463 m = timing_regex.match(s) 

464 if m: 

465 if m.group(3): 

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

467 elif m.group(19): 

468 phase_defs = [m.group(19)] 

469 else: 

470 phase_defs = [] 

471 

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

473 

474 offset = 0.0 

475 soff = m.group(27) 

476 if soff: 

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

478 if soff.endswith('S'): 

479 offset_is = 'slowness' 

480 elif soff.endswith('%'): 

481 offset_is = 'percent' 

482 

483 matched = True 

484 

485 else: 

486 m = timing_regex_old.match(s) 

487 if m: 

488 if m.group(3): 

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

490 elif m.group(5): 

491 phase_defs = [m.group(5)] 

492 else: 

493 phase_defs = [] 

494 

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

496 

497 offset = 0.0 

498 if m.group(6): 

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

500 

501 phase_defs = [ 

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

503 

504 matched = True 

505 

506 if not matched: 

507 raise InvalidTimingSpecification(s) 

508 

509 kwargs = dict( 

510 phase_defs=phase_defs, 

511 select=select, 

512 offset=offset, 

513 offset_is=offset_is) 

514 

515 SObject.__init__(self, **kwargs) 

516 

517 def __str__(self): 

518 s = [] 

519 if self.phase_defs: 

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

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

522 sphases = '{%s}' % sphases 

523 

524 if self.select: 

525 sphases = self.select + sphases 

526 

527 s.append(sphases) 

528 

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

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

531 if self.offset_is == 'slowness': 

532 s.append('S') 

533 elif self.offset_is == 'percent': 

534 s.append('%') 

535 

536 return ''.join(s) 

537 

538 def evaluate(self, get_phase, args): 

539 try: 

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

541 phase_offset = get_phase( 

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

543 offset = phase_offset(args) 

544 else: 

545 offset = self.offset 

546 

547 if self.phase_defs: 

548 phases = [ 

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

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

551 if self.offset_is == 'percent': 

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

553 for t in times if t is not None] 

554 else: 

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

556 

557 if not times: 

558 return None 

559 elif self.select == 'first': 

560 return min(times) 

561 elif self.select == 'last': 

562 return max(times) 

563 else: 

564 return times[0] 

565 else: 

566 return offset 

567 

568 except spit.OutOfBounds: 

569 raise OutOfBounds(args) 

570 

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

572 offset = Float.T(default=0.0) 

573 offset_is = String.T(optional=True) 

574 select = PhaseSelect.T( 

575 default='', 

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

577 tuple(PhaseSelect.choices))) 

578 

579 

580def mkdefs(s): 

581 defs = [] 

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

583 try: 

584 defs.append(float(x)) 

585 except ValueError: 

586 if x.startswith('!'): 

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

588 else: 

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

590 

591 return defs 

592 

593 

594class TPDef(Object): 

595 ''' 

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

597 ''' 

598 

599 id = StringID.T( 

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

601 definition = String.T( 

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

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

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

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

606 

607 @property 

608 def phases(self): 

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

610 if isinstance(x, cake.PhaseDef)] 

611 

612 @property 

613 def horizontal_velocities(self): 

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

615 

616 

617class OutOfBounds(Exception): 

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

619 Exception.__init__(self) 

620 self.values = values 

621 self.reason = reason 

622 self.context = None 

623 

624 def __str__(self): 

625 scontext = '' 

626 if self.context: 

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

628 

629 if self.reason: 

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

631 

632 if self.values: 

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

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

635 else: 

636 return 'out of bounds%s' % scontext 

637 

638 

639class MultiLocation(Object): 

640 

641 lats = Array.T( 

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

643 help='Latitudes of targets.') 

644 lons = Array.T( 

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

646 help='Longitude of targets.') 

647 north_shifts = Array.T( 

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

649 help='North shifts of targets.') 

650 east_shifts = Array.T( 

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

652 help='East shifts of targets.') 

653 elevation = Array.T( 

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

655 help='Elevations of targets.') 

656 

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

658 self._coords5 = None 

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

660 

661 @property 

662 def coords5(self): 

663 if self._coords5 is not None: 

664 return self._coords5 

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

666 self.elevation] 

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

668 if not sizes: 

669 raise AttributeError('Empty StaticTarget') 

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

671 raise AttributeError('Inconsistent coordinate sizes.') 

672 

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

674 for idx, p in enumerate(props): 

675 if p is not None: 

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

677 

678 return self._coords5 

679 

680 @property 

681 def ncoords(self): 

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

683 return 0 

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

685 

686 def get_latlon(self): 

687 ''' 

688 Get all coordinates as lat/lon. 

689 

690 :returns: Coordinates in Latitude, Longitude 

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

692 ''' 

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

694 for i in range(self.ncoords): 

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

696 return latlons 

697 

698 def get_corner_coords(self): 

699 ''' 

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

701 

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

703 :rtype: tuple 

704 ''' 

705 latlon = self.get_latlon() 

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

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

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

709 

710 

711class Receiver(Location): 

712 codes = Tuple.T( 

713 3, String.T(), 

714 optional=True, 

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

716 

717 def pyrocko_station(self): 

718 from pyrocko import model 

719 

720 lat, lon = self.effective_latlon 

721 

722 return model.Station( 

723 network=self.codes[0], 

724 station=self.codes[1], 

725 location=self.codes[2], 

726 lat=lat, 

727 lon=lon, 

728 depth=self.depth) 

729 

730 

731def g(x, d): 

732 if x is None: 

733 return d 

734 else: 

735 return x 

736 

737 

738class UnavailableScheme(Exception): 

739 pass 

740 

741 

742class InvalidNComponents(Exception): 

743 pass 

744 

745 

746class DiscretizedSource(Object): 

747 ''' 

748 Base class for discretized sources. 

749 

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

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

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

753 source distributions are represented by subclasses of the 

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

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

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

757 explosion/implosion type source distributions. The geometry calculations 

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

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

760 

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

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

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

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

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

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

767 ''' 

768 

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

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

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

772 lat = Float.T(optional=True) 

773 lon = Float.T(optional=True) 

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

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

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

777 dl = Float.T(optional=True) 

778 dw = Float.T(optional=True) 

779 nl = Float.T(optional=True) 

780 nw = Float.T(optional=True) 

781 

782 @classmethod 

783 def check_scheme(cls, scheme): 

784 ''' 

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

786 

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

788 supported by this discretized source class. 

789 ''' 

790 

791 if scheme not in cls.provided_schemes: 

792 raise UnavailableScheme( 

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

794 (cls.__name__, scheme)) 

795 

796 def __init__(self, **kwargs): 

797 Object.__init__(self, **kwargs) 

798 self._latlons = None 

799 

800 def __setattr__(self, name, value): 

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

802 'lats', 'lons'): 

803 self.__dict__['_latlons'] = None 

804 

805 Object.__setattr__(self, name, value) 

806 

807 def get_source_terms(self, scheme): 

808 raise NotImplementedError() 

809 

810 def make_weights(self, receiver, scheme): 

811 raise NotImplementedError() 

812 

813 @property 

814 def effective_latlons(self): 

815 ''' 

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

817 ''' 

818 

819 if self._latlons is None: 

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

821 if (self.north_shifts is not None and 

822 self.east_shifts is not None): 

823 self._latlons = orthodrome.ne_to_latlon( 

824 self.lats, self.lons, 

825 self.north_shifts, self.east_shifts) 

826 else: 

827 self._latlons = self.lats, self.lons 

828 else: 

829 lat = g(self.lat, 0.0) 

830 lon = g(self.lon, 0.0) 

831 self._latlons = orthodrome.ne_to_latlon( 

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

833 

834 return self._latlons 

835 

836 @property 

837 def effective_north_shifts(self): 

838 if self.north_shifts is not None: 

839 return self.north_shifts 

840 else: 

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

842 

843 @property 

844 def effective_east_shifts(self): 

845 if self.east_shifts is not None: 

846 return self.east_shifts 

847 else: 

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

849 

850 def same_origin(self, receiver): 

851 ''' 

852 Check if receiver has same reference point. 

853 ''' 

854 

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

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

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

858 

859 def azibazis_to(self, receiver): 

860 ''' 

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

862 points. 

863 ''' 

864 

865 if self.same_origin(receiver): 

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

867 receiver.north_shift - self.north_shifts) 

868 bazis = azis + 180. 

869 else: 

870 slats, slons = self.effective_latlons 

871 rlat, rlon = receiver.effective_latlon 

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

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

874 

875 return azis, bazis 

876 

877 def distances_to(self, receiver): 

878 ''' 

879 Compute distances to receiver for all contained points. 

880 ''' 

881 if self.same_origin(receiver): 

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

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

884 

885 else: 

886 slats, slons = self.effective_latlons 

887 rlat, rlon = receiver.effective_latlon 

888 return orthodrome.distance_accurate50m_numpy(slats, slons, 

889 rlat, rlon) 

890 

891 def element_coords(self, i): 

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

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

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

895 else: 

896 lat = self.lat 

897 lon = self.lon 

898 

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

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

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

902 

903 else: 

904 north_shift = east_shift = 0.0 

905 

906 return lat, lon, north_shift, east_shift 

907 

908 def coords5(self): 

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

910 

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

912 xs[:, 0] = self.lats 

913 xs[:, 1] = self.lons 

914 else: 

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

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

917 

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

919 xs[:, 2] = self.north_shifts 

920 xs[:, 3] = self.east_shifts 

921 

922 xs[:, 4] = self.depths 

923 

924 return xs 

925 

926 @property 

927 def nelements(self): 

928 return self.times.size 

929 

930 @classmethod 

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

932 ''' 

933 Combine several discretized source models. 

934 

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

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

937 factors and reference times of the parameterized (undiscretized) 

938 sources match or are accounted for. 

939 ''' 

940 

941 first = sources[0] 

942 

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

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

945 'sources of same type.') 

946 

947 latlons = [] 

948 for s in sources: 

949 latlons.append(s.effective_latlons) 

950 

951 lats, lons = num.hstack(latlons) 

952 

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

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

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

956 same_ref = num.all( 

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

958 else: 

959 same_ref = False 

960 

961 cat = num.concatenate 

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

963 print(times) 

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, one fixed receiver position. 

2529 

2530 * Cartesian 3D source volume around a reference point 

2531 

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

2533 source_east_shift, source_north_shift, component)`` 

2534 ''' 

2535 

2536 receiver = Receiver.T( 

2537 help='Receiver location') 

2538 

2539 source_origin = Location.T( 

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

2541 

2542 source_depth_min = Float.T( 

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

2544 

2545 source_depth_max = Float.T( 

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

2547 

2548 source_depth_delta = Float.T( 

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

2550 

2551 source_east_shift_min = Float.T( 

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

2553 

2554 source_east_shift_max = Float.T( 

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

2556 

2557 source_east_shift_delta = Float.T( 

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

2559 

2560 source_north_shift_min = Float.T( 

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

2562 

2563 source_north_shift_max = Float.T( 

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

2565 

2566 source_north_shift_delta = Float.T( 

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

2568 

2569 short_type = 'C' 

2570 

2571 provided_schemes = ['elastic18'] 

2572 

2573 def get_surface_distance(self, args): 

2574 _, source_east_shift, source_north_shift, _ = args 

2575 sorig = self.source_origin 

2576 sloc = Location( 

2577 lat=sorig.lat, 

2578 lon=sorig.lon, 

2579 north_shift=sorig.north_shift + source_north_shift, 

2580 east_shift=sorig.east_shift + source_east_shift) 

2581 

2582 return self.receiver.distance_to(sloc) 

2583 

2584 def get_distance(self, args): 

2585 sdepth, source_east_shift, source_north_shift, _ = args 

2586 sorig = self.source_origin 

2587 sloc = Location( 

2588 lat=sorig.lat, 

2589 lon=sorig.lon, 

2590 north_shift=sorig.north_shift + source_north_shift, 

2591 east_shift=sorig.east_shift + source_east_shift) 

2592 

2593 return self.receiver.distance_3d_to(sloc) 

2594 

2595 def get_source_depth(self, args): 

2596 return args[0] 

2597 

2598 def get_receiver_depth(self, args): 

2599 return self.receiver.depth 

2600 

2601 def get_source_depths(self): 

2602 return self.coords[0] 

2603 

2604 def _update(self): 

2605 self.mins = num.array([ 

2606 self.source_depth_min, 

2607 self.source_east_shift_min, 

2608 self.source_north_shift_min], 

2609 dtype=float) 

2610 

2611 self.maxs = num.array([ 

2612 self.source_depth_max, 

2613 self.source_east_shift_max, 

2614 self.source_north_shift_max], 

2615 dtype=float) 

2616 

2617 self.deltas = num.array([ 

2618 self.source_depth_delta, 

2619 self.source_east_shift_delta, 

2620 self.source_north_shift_delta], 

2621 dtype=float) 

2622 

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

2624 vicinity_eps).astype(int) + 1 

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

2626 self.deltat = 1.0/self.sample_rate 

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

2628 

2629 self.coords = tuple( 

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

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

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

2633 

2634 def _make_index_functions(self): 

2635 

2636 amin, bmin, cmin = self.mins 

2637 da, db, dc = self.deltas 

2638 na, nb, nc = self.ns 

2639 ng = self.ncomponents 

2640 

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

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

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

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

2645 try: 

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

2647 (na, nb, nc, ng)) 

2648 except ValueError: 

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

2650 

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

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

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

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

2655 

2656 try: 

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

2658 (na, nb, nc, ng)) 

2659 except ValueError: 

2660 raise OutOfBounds() 

2661 

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

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

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

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

2666 

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

2668 raise OutOfBounds() 

2669 

2670 indis = [] 

2671 weights = [] 

2672 for ia, va in ias: 

2673 iia = ia*nb*nc*ng 

2674 for ib, vb in ibs: 

2675 iib = ib*nc*ng 

2676 for ic, vc in ics: 

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

2678 weights.append(va*vb*vc) 

2679 

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

2681 

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

2683 

2684 xa = (a-amin) / da 

2685 xb = (b-bmin) / db 

2686 xc = (c-cmin) / dc 

2687 

2688 xa_fl = num.floor(xa) 

2689 xa_ce = num.ceil(xa) 

2690 xb_fl = num.floor(xb) 

2691 xb_ce = num.ceil(xb) 

2692 xc_fl = num.floor(xc) 

2693 xc_ce = num.ceil(xc) 

2694 va_fl = 1.0 - (xa - xa_fl) 

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

2696 vb_fl = 1.0 - (xb - xb_fl) 

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

2698 vc_fl = 1.0 - (xc - xc_fl) 

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

2700 

2701 ia_fl = xa_fl.astype(int) 

2702 ia_ce = xa_ce.astype(int) 

2703 ib_fl = xb_fl.astype(int) 

2704 ib_ce = xb_ce.astype(int) 

2705 ic_fl = xc_fl.astype(int) 

2706 ic_ce = xc_ce.astype(int) 

2707 

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

2709 raise OutOfBounds() 

2710 

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

2712 raise OutOfBounds() 

2713 

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

2715 raise OutOfBounds() 

2716 

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

2718 raise OutOfBounds() 

2719 

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

2721 raise OutOfBounds() 

2722 

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

2724 raise OutOfBounds() 

2725 

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

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

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

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

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

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

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

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

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

2735 

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

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

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

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

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

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

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

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

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

2745 

2746 return irecords, weights 

2747 

2748 self._index_function = index_function 

2749 self._indices_function = indices_function 

2750 self._vicinity_function = vicinity_function 

2751 self._vicinities_function = vicinities_function 

2752 

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

2754 nc = icomponents.size 

2755 

2756 dists = source.distances_to(self.source_origin) 

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

2758 

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

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

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

2762 

2763 n = dists.size 

2764 

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

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

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

2768 icomponents) 

2769 

2770 def make_indexing_args1(self, source, receiver): 

2771 dist = source.distance_to(self.source_origin) 

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

2773 

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

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

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

2777 

2778 return (source_depth, 

2779 source_east_shift, 

2780 source_north_shift) 

2781 

2782 @property 

2783 def short_extent(self): 

2784 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2785 self.source_depth_min/km, 

2786 self.source_depth_max/km, 

2787 self.source_depth_delta/km, 

2788 self.source_east_shift_min/km, 

2789 self.source_east_shift_max/km, 

2790 self.source_east_shift_delta/km, 

2791 self.source_north_shift_min/km, 

2792 self.source_north_shift_max/km, 

2793 self.source_north_shift_delta/km) 

2794 

2795 

2796class Weighting(Object): 

2797 factor = Float.T(default=1.0) 

2798 

2799 

2800class Taper(Object): 

2801 tmin = Timing.T() 

2802 tmax = Timing.T() 

2803 tfade = Float.T(default=0.0) 

2804 shape = StringChoice.T( 

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

2806 default='cos', 

2807 optional=True) 

2808 

2809 

2810class SimplePattern(SObject): 

2811 

2812 _pool = {} 

2813 

2814 def __init__(self, pattern): 

2815 self._pattern = pattern 

2816 SObject.__init__(self) 

2817 

2818 def __str__(self): 

2819 return self._pattern 

2820 

2821 @property 

2822 def regex(self): 

2823 pool = SimplePattern._pool 

2824 if self.pattern not in pool: 

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

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

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

2828 

2829 return pool[self.pattern] 

2830 

2831 def match(self, s): 

2832 return self.regex.match(s) 

2833 

2834 

2835class WaveformType(StringChoice): 

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

2837 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2838 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2839 

2840 

2841class ChannelSelection(Object): 

2842 pattern = SimplePattern.T(optional=True) 

2843 min_sample_rate = Float.T(optional=True) 

2844 max_sample_rate = Float.T(optional=True) 

2845 

2846 

2847class StationSelection(Object): 

2848 includes = SimplePattern.T() 

2849 excludes = SimplePattern.T() 

2850 distance_min = Float.T(optional=True) 

2851 distance_max = Float.T(optional=True) 

2852 azimuth_min = Float.T(optional=True) 

2853 azimuth_max = Float.T(optional=True) 

2854 

2855 

2856class WaveformSelection(Object): 

2857 channel_selection = ChannelSelection.T(optional=True) 

2858 station_selection = StationSelection.T(optional=True) 

2859 taper = Taper.T() 

2860 # filter = FrequencyResponse.T() 

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

2862 weighting = Weighting.T(optional=True) 

2863 sample_rate = Float.T(optional=True) 

2864 gf_store_id = StringID.T(optional=True) 

2865 

2866 

2867def indi12(x, n): 

2868 ''' 

2869 Get linear interpolation index and weight. 

2870 ''' 

2871 

2872 r = round(x) 

2873 if abs(r - x) < vicinity_eps: 

2874 i = int(r) 

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

2876 raise OutOfBounds() 

2877 

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

2879 else: 

2880 f = math.floor(x) 

2881 i = int(f) 

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

2883 raise OutOfBounds() 

2884 

2885 v = x-f 

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

2887 

2888 

2889def float_or_none(s): 

2890 units = { 

2891 'k': 1e3, 

2892 'M': 1e6, 

2893 } 

2894 

2895 factor = 1.0 

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

2897 factor = units[s[-1]] 

2898 s = s[:-1] 

2899 if not s: 

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

2901 

2902 if s: 

2903 return float(s) * factor 

2904 else: 

2905 return None 

2906 

2907 

2908class GridSpecError(Exception): 

2909 def __init__(self, s): 

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

2911 

2912 

2913def parse_grid_spec(spec): 

2914 try: 

2915 result = [] 

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

2917 t = dspec.split('@') 

2918 num = start = stop = step = None 

2919 if len(t) == 2: 

2920 num = int(t[1]) 

2921 if num <= 0: 

2922 raise GridSpecError(spec) 

2923 

2924 elif len(t) > 2: 

2925 raise GridSpecError(spec) 

2926 

2927 s = t[0] 

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

2929 if len(v) == 1: 

2930 start = stop = v[0] 

2931 if len(v) >= 2: 

2932 start, stop = v[0:2] 

2933 if len(v) == 3: 

2934 step = v[2] 

2935 

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

2937 raise GridSpecError(spec) 

2938 

2939 if step == 0.0: 

2940 raise GridSpecError(spec) 

2941 

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

2943 

2944 except ValueError: 

2945 raise GridSpecError(spec) 

2946 

2947 return result 

2948 

2949 

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

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

2952 if start is None: 

2953 start = ma if swap else mi 

2954 if stop is None: 

2955 stop = mi if swap else ma 

2956 if step is None: 

2957 step = -inc if ma < mi else inc 

2958 if num is None: 

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

2960 raise GridSpecError() 

2961 

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

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

2964 if abs(stop-stop2) > eps: 

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

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

2967 else: 

2968 stop = stop2 

2969 

2970 if start == stop: 

2971 num = 1 

2972 

2973 return start, stop, num 

2974 

2975 

2976def nditer_outer(x): 

2977 return num.nditer( 

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

2979 

2980 

2981def nodes(xs): 

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

2983 nnodes = num.prod(ns) 

2984 ndim = len(xs) 

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

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

2987 x = xs[idim] 

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

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

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

2991 

2992 return nodes 

2993 

2994 

2995def filledi(x, n): 

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

2997 a.fill(x) 

2998 return a 

2999 

3000 

3001config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3002 

3003discretized_source_classes = [ 

3004 DiscretizedExplosionSource, 

3005 DiscretizedSFSource, 

3006 DiscretizedMTSource, 

3007 DiscretizedPorePressureSource] 

3008 

3009 

3010__all__ = ''' 

3011Earthmodel1D 

3012StringID 

3013ScopeType 

3014WaveformType 

3015QuantityType 

3016NearfieldTermsType 

3017Reference 

3018Region 

3019CircularRegion 

3020RectangularRegion 

3021PhaseSelect 

3022InvalidTimingSpecification 

3023Timing 

3024TPDef 

3025OutOfBounds 

3026Location 

3027Receiver 

3028'''.split() + [ 

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

3030ComponentScheme 

3031component_scheme_to_description 

3032component_schemes 

3033Config 

3034GridSpecError 

3035Weighting 

3036Taper 

3037SimplePattern 

3038WaveformType 

3039ChannelSelection 

3040StationSelection 

3041WaveformSelection 

3042nditer_outer 

3043dump 

3044load 

3045discretized_source_classes 

3046config_type_classes 

3047UnavailableScheme 

3048InterpolationMethod 

3049SeismosizerTrace 

3050SeismosizerResult 

3051Result 

3052StaticResult 

3053'''.split()