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 depths = cat([s.depths for s in sources]) 

964 

965 if same_ref: 

966 lat = first.lat 

967 lon = first.lon 

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

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

970 lats = None 

971 lons = None 

972 else: 

973 lat = None 

974 lon = None 

975 north_shifts = None 

976 east_shifts = None 

977 

978 return cls( 

979 times=times, lat=lat, lon=lon, lats=lats, lons=lons, 

980 north_shifts=north_shifts, east_shifts=east_shifts, 

981 depths=depths, **kwargs) 

982 

983 def centroid_position(self): 

984 moments = self.moments() 

985 norm = num.sum(moments) 

986 if norm != 0.0: 

987 w = moments / num.sum(moments) 

988 else: 

989 w = num.ones(moments.size) 

990 

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

992 lats, lons = self.effective_latlons 

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

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

995 else: 

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

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

998 

999 cn = num.sum(n*w) 

1000 ce = num.sum(e*w) 

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

1002 

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

1004 lat = clat 

1005 lon = clon 

1006 north_shift = 0. 

1007 east_shift = 0. 

1008 else: 

1009 lat = g(self.lat, 0.0) 

1010 lon = g(self.lon, 0.0) 

1011 north_shift = cn 

1012 east_shift = ce 

1013 

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

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

1016 return tuple(float(x) for x in 

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

1018 

1019 

1020class DiscretizedExplosionSource(DiscretizedSource): 

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

1022 

1023 provided_schemes = ( 

1024 'elastic2', 

1025 'elastic8', 

1026 'elastic10', 

1027 ) 

1028 

1029 def get_source_terms(self, scheme): 

1030 self.check_scheme(scheme) 

1031 

1032 if scheme == 'elastic2': 

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

1034 

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

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

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

1038 return m6s 

1039 else: 

1040 assert False 

1041 

1042 def make_weights(self, receiver, scheme): 

1043 self.check_scheme(scheme) 

1044 

1045 azis, bazis = self.azibazis_to(receiver) 

1046 

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

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

1049 

1050 m0s = self.m0s 

1051 n = azis.size 

1052 

1053 cat = num.concatenate 

1054 rep = num.repeat 

1055 

1056 if scheme == 'elastic2': 

1057 w_n = cb*m0s 

1058 g_n = filledi(0, n) 

1059 w_e = sb*m0s 

1060 g_e = filledi(0, n) 

1061 w_d = m0s 

1062 g_d = filledi(1, n) 

1063 

1064 elif scheme == 'elastic8': 

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

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

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

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

1069 w_d = cat((m0s, m0s)) 

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

1071 

1072 elif scheme == 'elastic10': 

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

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

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

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

1077 w_d = cat((m0s, m0s, m0s)) 

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

1079 

1080 else: 

1081 assert False 

1082 

1083 return ( 

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

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

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

1087 ) 

1088 

1089 def split(self): 

1090 from pyrocko.gf.seismosizer import ExplosionSource 

1091 sources = [] 

1092 for i in range(self.nelements): 

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

1094 sources.append(ExplosionSource( 

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

1096 lat=lat, 

1097 lon=lon, 

1098 north_shift=north_shift, 

1099 east_shift=east_shift, 

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

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

1102 

1103 return sources 

1104 

1105 def moments(self): 

1106 return self.m0s 

1107 

1108 def centroid(self): 

1109 from pyrocko.gf.seismosizer import ExplosionSource 

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

1111 self.centroid_position() 

1112 

1113 return ExplosionSource( 

1114 time=time, 

1115 lat=lat, 

1116 lon=lon, 

1117 north_shift=north_shift, 

1118 east_shift=east_shift, 

1119 depth=depth, 

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

1121 

1122 @classmethod 

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

1124 ''' 

1125 Combine several discretized source models. 

1126 

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

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

1129 factors and reference times of the parameterized (undiscretized) 

1130 sources match or are accounted for. 

1131 ''' 

1132 

1133 if 'm0s' not in kwargs: 

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

1135 

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

1137 **kwargs) 

1138 

1139 

1140class DiscretizedSFSource(DiscretizedSource): 

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

1142 

1143 provided_schemes = ( 

1144 'elastic5', 

1145 ) 

1146 

1147 def get_source_terms(self, scheme): 

1148 self.check_scheme(scheme) 

1149 

1150 return self.forces 

1151 

1152 def make_weights(self, receiver, scheme): 

1153 self.check_scheme(scheme) 

1154 

1155 azis, bazis = self.azibazis_to(receiver) 

1156 

1157 sa = num.sin(azis*d2r) 

1158 ca = num.cos(azis*d2r) 

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

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

1161 

1162 forces = self.forces 

1163 fn = forces[:, 0] 

1164 fe = forces[:, 1] 

1165 fd = forces[:, 2] 

1166 

1167 f0 = fd 

1168 f1 = ca * fn + sa * fe 

1169 f2 = ca * fe - sa * fn 

1170 

1171 n = azis.size 

1172 

1173 if scheme == 'elastic5': 

1174 ioff = 0 

1175 

1176 cat = num.concatenate 

1177 rep = num.repeat 

1178 

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

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

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

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

1183 w_d = cat((f0, f1)) 

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

1185 

1186 return ( 

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

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

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

1190 ) 

1191 

1192 @classmethod 

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

1194 ''' 

1195 Combine several discretized source models. 

1196 

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

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

1199 factors and reference times of the parameterized (undiscretized) 

1200 sources match or are accounted for. 

1201 ''' 

1202 

1203 if 'forces' not in kwargs: 

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

1205 

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

1207 

1208 def moments(self): 

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

1210 

1211 def centroid(self): 

1212 from pyrocko.gf.seismosizer import SFSource 

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

1214 self.centroid_position() 

1215 

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

1217 return SFSource( 

1218 time=time, 

1219 lat=lat, 

1220 lon=lon, 

1221 north_shift=north_shift, 

1222 east_shift=east_shift, 

1223 depth=depth, 

1224 fn=fn, 

1225 fe=fe, 

1226 fd=fd) 

1227 

1228 

1229class DiscretizedMTSource(DiscretizedSource): 

1230 m6s = Array.T( 

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

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

1233 

1234 provided_schemes = ( 

1235 'elastic8', 

1236 'elastic10', 

1237 'elastic18', 

1238 ) 

1239 

1240 def get_source_terms(self, scheme): 

1241 self.check_scheme(scheme) 

1242 return self.m6s 

1243 

1244 def make_weights(self, receiver, scheme): 

1245 self.check_scheme(scheme) 

1246 

1247 m6s = self.m6s 

1248 n = m6s.shape[0] 

1249 rep = num.repeat 

1250 

1251 if scheme == 'elastic18': 

1252 w_n = m6s.flatten() 

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

1254 w_e = m6s.flatten() 

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

1256 w_d = m6s.flatten() 

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

1258 

1259 else: 

1260 azis, bazis = self.azibazis_to(receiver) 

1261 

1262 sa = num.sin(azis*d2r) 

1263 ca = num.cos(azis*d2r) 

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

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

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

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

1268 

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

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

1271 f2 = m6s[:, 2] 

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

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

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

1275 

1276 cat = num.concatenate 

1277 

1278 if scheme == 'elastic8': 

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

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

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

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

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

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

1285 

1286 elif scheme == 'elastic10': 

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

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

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

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

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

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

1293 

1294 else: 

1295 assert False 

1296 

1297 return ( 

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

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

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

1301 ) 

1302 

1303 def split(self): 

1304 from pyrocko.gf.seismosizer import MTSource 

1305 sources = [] 

1306 for i in range(self.nelements): 

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

1308 sources.append(MTSource( 

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

1310 lat=lat, 

1311 lon=lon, 

1312 north_shift=north_shift, 

1313 east_shift=east_shift, 

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

1315 m6=self.m6s[i])) 

1316 

1317 return sources 

1318 

1319 def moments(self): 

1320 moments = num.array( 

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

1322 for m6 in self.m6s]) 

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

1324 

1325 def get_moment_rate(self, deltat=None): 

1326 moments = self.moments() 

1327 times = self.times 

1328 times -= times.min() 

1329 

1330 t_max = times.max() 

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

1332 mom_times[mom_times > t_max] = t_max 

1333 

1334 # Right open histrogram bins 

1335 mom, _ = num.histogram( 

1336 times, 

1337 bins=mom_times, 

1338 weights=moments) 

1339 

1340 deltat = num.diff(mom_times) 

1341 mom_rate = mom / deltat 

1342 

1343 return mom_rate, mom_times[1:] 

1344 

1345 def centroid(self): 

1346 from pyrocko.gf.seismosizer import MTSource 

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

1348 self.centroid_position() 

1349 

1350 return MTSource( 

1351 time=time, 

1352 lat=lat, 

1353 lon=lon, 

1354 north_shift=north_shift, 

1355 east_shift=east_shift, 

1356 depth=depth, 

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

1358 

1359 @classmethod 

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

1361 ''' 

1362 Combine several discretized source models. 

1363 

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

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

1366 factors and reference times of the parameterized (undiscretized) 

1367 sources match or are accounted for. 

1368 ''' 

1369 

1370 if 'm6s' not in kwargs: 

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

1372 

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

1374 

1375 

1376class DiscretizedPorePressureSource(DiscretizedSource): 

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

1378 

1379 provided_schemes = ( 

1380 'poroelastic10', 

1381 ) 

1382 

1383 def get_source_terms(self, scheme): 

1384 self.check_scheme(scheme) 

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

1386 

1387 def make_weights(self, receiver, scheme): 

1388 self.check_scheme(scheme) 

1389 

1390 azis, bazis = self.azibazis_to(receiver) 

1391 

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

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

1394 

1395 pp = self.pp 

1396 n = bazis.size 

1397 

1398 w_un = cb*pp 

1399 g_un = filledi(1, n) 

1400 w_ue = sb*pp 

1401 g_ue = filledi(1, n) 

1402 w_ud = pp 

1403 g_ud = filledi(0, n) 

1404 

1405 w_tn = cb*pp 

1406 g_tn = filledi(6, n) 

1407 w_te = sb*pp 

1408 g_te = filledi(6, n) 

1409 

1410 w_pp = pp 

1411 g_pp = filledi(7, n) 

1412 

1413 w_dvn = cb*pp 

1414 g_dvn = filledi(9, n) 

1415 w_dve = sb*pp 

1416 g_dve = filledi(9, n) 

1417 w_dvd = pp 

1418 g_dvd = filledi(8, n) 

1419 

1420 return ( 

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

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

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

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

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

1426 ('pore_pressure', w_pp, g_pp), 

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

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

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

1430 ) 

1431 

1432 def moments(self): 

1433 return self.pp 

1434 

1435 def centroid(self): 

1436 from pyrocko.gf.seismosizer import PorePressurePointSource 

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

1438 self.centroid_position() 

1439 

1440 return PorePressurePointSource( 

1441 time=time, 

1442 lat=lat, 

1443 lon=lon, 

1444 north_shift=north_shift, 

1445 east_shift=east_shift, 

1446 depth=depth, 

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

1448 

1449 @classmethod 

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

1451 ''' 

1452 Combine several discretized source models. 

1453 

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

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

1456 factors and reference times of the parameterized (undiscretized) 

1457 sources match or are accounted for. 

1458 ''' 

1459 

1460 if 'pp' not in kwargs: 

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

1462 

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

1464 **kwargs) 

1465 

1466 

1467class Region(Object): 

1468 name = String.T(optional=True) 

1469 

1470 

1471class RectangularRegion(Region): 

1472 lat_min = Float.T() 

1473 lat_max = Float.T() 

1474 lon_min = Float.T() 

1475 lon_max = Float.T() 

1476 

1477 

1478class CircularRegion(Region): 

1479 lat = Float.T() 

1480 lon = Float.T() 

1481 radius = Float.T() 

1482 

1483 

1484class Config(Object): 

1485 ''' 

1486 Green's function store meta information. 

1487 

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

1489 configuration types are: 

1490 

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

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

1493 

1494 * Problem is invariant to horizontal translations and rotations around 

1495 vertical axis. 

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

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

1498 component)`` 

1499 

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

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

1502 

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

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

1505 receiver_depth, component)`` 

1506 

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

1508 constraints but fixed receiver positions 

1509 

1510 * Cartesian source volume around a reference point 

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

1512 source_east_shift, source_north_shift, component)`` 

1513 ''' 

1514 

1515 id = StringID.T( 

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

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

1518 'letter.') 

1519 

1520 derived_from_id = StringID.T( 

1521 optional=True, 

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

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

1524 

1525 version = String.T( 

1526 default='1.0', 

1527 optional=True, 

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

1529 

1530 modelling_code_id = StringID.T( 

1531 optional=True, 

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

1533 

1534 author = Unicode.T( 

1535 optional=True, 

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

1537 

1538 author_email = String.T( 

1539 optional=True, 

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

1541 

1542 created_time = Timestamp.T( 

1543 optional=True, 

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

1545 

1546 regions = List.T( 

1547 Region.T(), 

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

1549 

1550 scope_type = ScopeType.T( 

1551 optional=True, 

1552 help='Distance range scope of the store ' 

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

1554 

1555 waveform_type = WaveType.T( 

1556 optional=True, 

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

1558 

1559 nearfield_terms = NearfieldTermsType.T( 

1560 optional=True, 

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

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

1563 

1564 description = String.T( 

1565 optional=True, 

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

1567 

1568 references = List.T( 

1569 Reference.T(), 

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

1571 'related work.') 

1572 

1573 earthmodel_1d = Earthmodel1D.T( 

1574 optional=True, 

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

1576 

1577 earthmodel_receiver_1d = Earthmodel1D.T( 

1578 optional=True, 

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

1580 

1581 can_interpolate_source = Bool.T( 

1582 optional=True, 

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

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

1585 

1586 can_interpolate_receiver = Bool.T( 

1587 optional=True, 

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

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

1590 

1591 frequency_min = Float.T( 

1592 optional=True, 

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

1594 

1595 frequency_max = Float.T( 

1596 optional=True, 

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

1598 

1599 sample_rate = Float.T( 

1600 optional=True, 

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

1602 

1603 factor = Float.T( 

1604 default=1.0, 

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

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

1607 optional=True) 

1608 

1609 component_scheme = ComponentScheme.T( 

1610 default='elastic10', 

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

1612 

1613 stored_quantity = QuantityType.T( 

1614 optional=True, 

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

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

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

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

1619 

1620 tabulated_phases = List.T( 

1621 TPDef.T(), 

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

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

1624 

1625 ncomponents = Int.T( 

1626 optional=True, 

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

1628 

1629 uuid = String.T( 

1630 optional=True, 

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

1632 'GF store for practical purposes.') 

1633 

1634 reference = String.T( 

1635 optional=True, 

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

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

1638 

1639 def __init__(self, **kwargs): 

1640 self._do_auto_updates = False 

1641 Object.__init__(self, **kwargs) 

1642 self._index_function = None 

1643 self._indices_function = None 

1644 self._vicinity_function = None 

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

1646 self._do_auto_updates = True 

1647 self.update() 

1648 

1649 def check_ncomponents(self): 

1650 ncomponents = component_scheme_to_description[ 

1651 self.component_scheme].ncomponents 

1652 

1653 if self.ncomponents is None: 

1654 self.ncomponents = ncomponents 

1655 elif ncomponents != self.ncomponents: 

1656 raise InvalidNComponents( 

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

1658 self.ncomponents, self.component_scheme)) 

1659 

1660 def __setattr__(self, name, value): 

1661 Object.__setattr__(self, name, value) 

1662 try: 

1663 self.T.get_property(name) 

1664 if self._do_auto_updates: 

1665 self.update() 

1666 

1667 except ValueError: 

1668 pass 

1669 

1670 def update(self): 

1671 self.check_ncomponents() 

1672 self._update() 

1673 self._make_index_functions() 

1674 

1675 def irecord(self, *args): 

1676 return self._index_function(*args) 

1677 

1678 def irecords(self, *args): 

1679 return self._indices_function(*args) 

1680 

1681 def vicinity(self, *args): 

1682 return self._vicinity_function(*args) 

1683 

1684 def vicinities(self, *args): 

1685 return self._vicinities_function(*args) 

1686 

1687 def grid_interpolation_coefficients(self, *args): 

1688 return self._grid_interpolation_coefficients(*args) 

1689 

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

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

1692 

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

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

1695 

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

1697 i = 0 

1698 arrs = [] 

1699 ntotal = 1 

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

1701 if gdef and len(gdef) > i: 

1702 sssn = gdef[i] 

1703 else: 

1704 sssn = (None,)*4 

1705 

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

1707 ntotal *= len(arr) 

1708 

1709 arrs.append(arr) 

1710 i += 1 

1711 

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

1713 return nditer_outer(arrs[:level]) 

1714 

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

1716 nthreads=0): 

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

1718 

1719 out = [] 

1720 delays = source.times 

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

1722 receiver, 

1723 self.component_scheme): 

1724 

1725 weights *= self.factor 

1726 

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

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

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

1730 

1731 return out 

1732 

1733 def short_info(self): 

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

1735 

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

1737 interpolation=None): 

1738 ''' 

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

1740 

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

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

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

1744 ``(lat, lon)`` 

1745 :param interpolation: interpolation method. Choose from 

1746 ``('nearest_neighbor', 'multilinear')`` 

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

1748 point 

1749 

1750 The default implementation retrieves and interpolates the shear moduli 

1751 from the contained 1D velocity profile. 

1752 ''' 

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

1754 parameter='shear_moduli', 

1755 interpolation=interpolation) 

1756 

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

1758 interpolation=None): 

1759 ''' 

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

1761 

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

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

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

1765 ``(lat, lon)`` 

1766 :param interpolation: interpolation method. Choose from 

1767 ``('nearest_neighbor', 'multilinear')`` 

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

1769 point 

1770 

1771 The default implementation retrieves and interpolates the lambda moduli 

1772 from the contained 1D velocity profile. 

1773 ''' 

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

1775 parameter='lambda_moduli', 

1776 interpolation=interpolation) 

1777 

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

1779 interpolation=None): 

1780 ''' 

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

1782 

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

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

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

1786 ``(lat, lon)`` 

1787 :param interpolation: interpolation method. Choose from 

1788 ``('nearest_neighbor', 'multilinear')`` 

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

1790 point 

1791 

1792 The default implementation retrieves and interpolates the lambda moduli 

1793 from the contained 1D velocity profile. 

1794 ''' 

1795 lambda_moduli = self.get_material_property( 

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

1797 interpolation=interpolation) 

1798 shear_moduli = self.get_material_property( 

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

1800 interpolation=interpolation) 

1801 return lambda_moduli + (2 / 3) * shear_moduli 

1802 

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

1804 ''' 

1805 Get Vs at given points from contained velocity model. 

1806 

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

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

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

1810 ``(lat, lon)`` 

1811 :param interpolation: interpolation method. Choose from 

1812 ``('nearest_neighbor', 'multilinear')`` 

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

1814 point 

1815 

1816 The default implementation retrieves and interpolates Vs 

1817 from the contained 1D velocity profile. 

1818 ''' 

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

1820 parameter='vs', 

1821 interpolation=interpolation) 

1822 

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

1824 ''' 

1825 Get Vp at given points from contained velocity model. 

1826 

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

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

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

1830 ``(lat, lon)`` 

1831 :param interpolation: interpolation method. Choose from 

1832 ``('nearest_neighbor', 'multilinear')`` 

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

1834 point 

1835 

1836 The default implementation retrieves and interpolates Vp 

1837 from the contained 1D velocity profile. 

1838 ''' 

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

1840 parameter='vp', 

1841 interpolation=interpolation) 

1842 

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

1844 ''' 

1845 Get rho at given points from contained velocity model. 

1846 

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

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

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

1850 ``(lat, lon)`` 

1851 :param interpolation: interpolation method. Choose from 

1852 ``('nearest_neighbor', 'multilinear')`` 

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

1854 point 

1855 

1856 The default implementation retrieves and interpolates rho 

1857 from the contained 1D velocity profile. 

1858 ''' 

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

1860 parameter='rho', 

1861 interpolation=interpolation) 

1862 

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

1864 interpolation=None): 

1865 

1866 if interpolation is None: 

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

1868 'multilinear', 'nearest_neighbor') 

1869 

1870 earthmod = self.earthmodel_1d 

1871 store_depth_profile = self.get_source_depths() 

1872 z_profile = earthmod.profile('z') 

1873 

1874 if parameter == 'vs': 

1875 vs_profile = earthmod.profile('vs') 

1876 profile = num.interp( 

1877 store_depth_profile, z_profile, vs_profile) 

1878 

1879 elif parameter == 'vp': 

1880 vp_profile = earthmod.profile('vp') 

1881 profile = num.interp( 

1882 store_depth_profile, z_profile, vp_profile) 

1883 

1884 elif parameter == 'rho': 

1885 rho_profile = earthmod.profile('rho') 

1886 

1887 profile = num.interp( 

1888 store_depth_profile, z_profile, rho_profile) 

1889 

1890 elif parameter == 'shear_moduli': 

1891 vs_profile = earthmod.profile('vs') 

1892 rho_profile = earthmod.profile('rho') 

1893 

1894 store_vs_profile = num.interp( 

1895 store_depth_profile, z_profile, vs_profile) 

1896 store_rho_profile = num.interp( 

1897 store_depth_profile, z_profile, rho_profile) 

1898 

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

1900 

1901 elif parameter == 'lambda_moduli': 

1902 vs_profile = earthmod.profile('vs') 

1903 vp_profile = earthmod.profile('vp') 

1904 rho_profile = earthmod.profile('rho') 

1905 

1906 store_vs_profile = num.interp( 

1907 store_depth_profile, z_profile, vs_profile) 

1908 store_vp_profile = num.interp( 

1909 store_depth_profile, z_profile, vp_profile) 

1910 store_rho_profile = num.interp( 

1911 store_depth_profile, z_profile, rho_profile) 

1912 

1913 profile = store_rho_profile * ( 

1914 num.power(store_vp_profile, 2) - 

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

1916 else: 

1917 raise TypeError( 

1918 'parameter %s not available' % parameter) 

1919 

1920 if interpolation == 'multilinear': 

1921 kind = 'linear' 

1922 elif interpolation == 'nearest_neighbor': 

1923 kind = 'nearest' 

1924 else: 

1925 raise TypeError( 

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

1927 

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

1929 

1930 try: 

1931 return interpolator(points[:, 2]) 

1932 except ValueError: 

1933 raise OutOfBounds() 

1934 

1935 def is_static(self): 

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

1937 if self.modelling_code_id.startswith(code): 

1938 return True 

1939 return False 

1940 

1941 def is_dynamic(self): 

1942 return not self.is_static() 

1943 

1944 def get_source_depths(self): 

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

1946 

1947 def get_tabulated_phase(self, phase_id): 

1948 ''' 

1949 Get tabulated phase definition. 

1950 ''' 

1951 

1952 for pdef in self.tabulated_phases: 

1953 if pdef.id == phase_id: 

1954 return pdef 

1955 

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

1957 

1958 def fix_ttt_holes(self, sptree, mode): 

1959 raise StoreError( 

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

1961 % self.short_type) 

1962 

1963 

1964class ConfigTypeA(Config): 

1965 ''' 

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

1967 

1968 * Problem is invariant to horizontal translations and rotations around 

1969 vertical axis. 

1970 

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

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

1973 component)`` 

1974 

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

1976 points. 

1977 ''' 

1978 

1979 receiver_depth = Float.T( 

1980 default=0.0, 

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

1982 

1983 source_depth_min = Float.T( 

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

1985 

1986 source_depth_max = Float.T( 

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

1988 

1989 source_depth_delta = Float.T( 

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

1991 

1992 distance_min = Float.T( 

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

1994 

1995 distance_max = Float.T( 

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

1997 

1998 distance_delta = Float.T( 

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

2000 

2001 short_type = 'A' 

2002 

2003 provided_schemes = [ 

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

2005 

2006 def get_surface_distance(self, args): 

2007 return args[1] 

2008 

2009 def get_distance(self, args): 

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

2011 

2012 def get_source_depth(self, args): 

2013 return args[0] 

2014 

2015 def get_source_depths(self): 

2016 return self.coords[0] 

2017 

2018 def get_receiver_depth(self, args): 

2019 return self.receiver_depth 

2020 

2021 def _update(self): 

2022 self.mins = num.array( 

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

2024 self.maxs = num.array( 

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

2026 self.deltas = num.array( 

2027 [self.source_depth_delta, self.distance_delta], 

2028 dtype=float) 

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

2030 vicinity_eps).astype(int) + 1 

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

2032 self.deltat = 1.0/self.sample_rate 

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

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

2035 (mi, ma, n) in 

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

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

2038 

2039 self.nsource_depths, self.ndistances = self.ns 

2040 

2041 def _make_index_functions(self): 

2042 

2043 amin, bmin = self.mins 

2044 da, db = self.deltas 

2045 na, nb = self.ns 

2046 

2047 ng = self.ncomponents 

2048 

2049 def index_function(a, b, ig): 

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

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

2052 try: 

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

2054 except ValueError: 

2055 raise OutOfBounds() 

2056 

2057 def indices_function(a, b, ig): 

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

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

2060 try: 

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

2062 except ValueError: 

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

2064 try: 

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

2066 except ValueError: 

2067 raise OutOfBounds() 

2068 

2069 def grid_interpolation_coefficients(a, b): 

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

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

2072 return ias, ibs 

2073 

2074 def vicinity_function(a, b, ig): 

2075 ias, ibs = grid_interpolation_coefficients(a, b) 

2076 

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

2078 raise OutOfBounds() 

2079 

2080 indis = [] 

2081 weights = [] 

2082 for ia, va in ias: 

2083 iia = ia*nb*ng 

2084 for ib, vb in ibs: 

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

2086 weights.append(va*vb) 

2087 

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

2089 

2090 def vicinities_function(a, b, ig): 

2091 

2092 xa = (a - amin) / da 

2093 xb = (b - bmin) / db 

2094 

2095 xa_fl = num.floor(xa) 

2096 xa_ce = num.ceil(xa) 

2097 xb_fl = num.floor(xb) 

2098 xb_ce = num.ceil(xb) 

2099 va_fl = 1.0 - (xa - xa_fl) 

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

2101 vb_fl = 1.0 - (xb - xb_fl) 

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

2103 

2104 ia_fl = xa_fl.astype(int) 

2105 ia_ce = xa_ce.astype(int) 

2106 ib_fl = xb_fl.astype(int) 

2107 ib_ce = xb_ce.astype(int) 

2108 

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

2110 raise OutOfBounds() 

2111 

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

2113 raise OutOfBounds() 

2114 

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

2116 raise OutOfBounds() 

2117 

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

2119 raise OutOfBounds() 

2120 

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

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

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

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

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

2126 

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

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

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

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

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

2132 

2133 return irecords, weights 

2134 

2135 self._index_function = index_function 

2136 self._indices_function = indices_function 

2137 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2138 self._vicinity_function = vicinity_function 

2139 self._vicinities_function = vicinities_function 

2140 

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

2142 nc = icomponents.size 

2143 dists = source.distances_to(receiver) 

2144 n = dists.size 

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

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

2147 icomponents) 

2148 

2149 def make_indexing_args1(self, source, receiver): 

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

2151 

2152 @property 

2153 def short_extent(self): 

2154 return '%g:%g:%g x %g:%g:%g' % ( 

2155 self.source_depth_min/km, 

2156 self.source_depth_max/km, 

2157 self.source_depth_delta/km, 

2158 self.distance_min/km, 

2159 self.distance_max/km, 

2160 self.distance_delta/km) 

2161 

2162 def fix_ttt_holes(self, sptree, mode): 

2163 from pyrocko import eikonal_ext, spit 

2164 

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

2166 

2167 delta = self.deltas[-1] 

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

2169 

2170 nsources, ndistances = self.ns 

2171 

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

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

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

2175 

2176 speeds = self.get_material_property( 

2177 0., 0., points, 

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

2179 interpolation='multilinear') 

2180 

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

2182 

2183 times = sptree.interpolate_many(nodes) 

2184 

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

2186 times = times.reshape(speeds.shape) 

2187 

2188 try: 

2189 eikonal_ext.eikonal_solver_fmm_cartesian( 

2190 speeds, times, delta) 

2191 except eikonal_ext.EikonalExtError as e: 

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

2193 logger.debug( 

2194 'Got a warning from eikonal solver ' 

2195 '- may be ok...') 

2196 else: 

2197 raise 

2198 

2199 def func(x): 

2200 ibs, ics = \ 

2201 self.grid_interpolation_coefficients(*x) 

2202 

2203 t = 0 

2204 for ib, vb in ibs: 

2205 for ic, vc in ics: 

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

2207 

2208 return t 

2209 

2210 return spit.SPTree( 

2211 f=func, 

2212 ftol=sptree.ftol, 

2213 xbounds=sptree.xbounds, 

2214 xtols=sptree.xtols) 

2215 

2216 

2217class ConfigTypeB(Config): 

2218 ''' 

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

2220 

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

2222 receiver depth 

2223 

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

2225 receiver_distance, component)`` 

2226 ''' 

2227 

2228 receiver_depth_min = Float.T( 

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

2230 

2231 receiver_depth_max = Float.T( 

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

2233 

2234 receiver_depth_delta = Float.T( 

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

2236 

2237 source_depth_min = Float.T( 

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

2239 

2240 source_depth_max = Float.T( 

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

2242 

2243 source_depth_delta = Float.T( 

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

2245 

2246 distance_min = Float.T( 

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

2248 

2249 distance_max = Float.T( 

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

2251 

2252 distance_delta = Float.T( 

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

2254 

2255 short_type = 'B' 

2256 

2257 provided_schemes = [ 

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

2259 

2260 def get_distance(self, args): 

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

2262 

2263 def get_surface_distance(self, args): 

2264 return args[2] 

2265 

2266 def get_source_depth(self, args): 

2267 return args[1] 

2268 

2269 def get_receiver_depth(self, args): 

2270 return args[0] 

2271 

2272 def get_source_depths(self): 

2273 return self.coords[1] 

2274 

2275 def _update(self): 

2276 self.mins = num.array([ 

2277 self.receiver_depth_min, 

2278 self.source_depth_min, 

2279 self.distance_min], 

2280 dtype=float) 

2281 

2282 self.maxs = num.array([ 

2283 self.receiver_depth_max, 

2284 self.source_depth_max, 

2285 self.distance_max], 

2286 dtype=float) 

2287 

2288 self.deltas = num.array([ 

2289 self.receiver_depth_delta, 

2290 self.source_depth_delta, 

2291 self.distance_delta], 

2292 dtype=float) 

2293 

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

2295 vicinity_eps).astype(int) + 1 

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

2297 self.deltat = 1.0/self.sample_rate 

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

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

2300 (mi, ma, n) in 

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

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

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

2304 

2305 def _make_index_functions(self): 

2306 

2307 amin, bmin, cmin = self.mins 

2308 da, db, dc = self.deltas 

2309 na, nb, nc = self.ns 

2310 ng = self.ncomponents 

2311 

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

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

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

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

2316 try: 

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

2318 (na, nb, nc, ng)) 

2319 except ValueError: 

2320 raise OutOfBounds() 

2321 

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

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

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

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

2326 try: 

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

2328 (na, nb, nc, ng)) 

2329 except ValueError: 

2330 raise OutOfBounds() 

2331 

2332 def grid_interpolation_coefficients(a, b, c): 

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

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

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

2336 return ias, ibs, ics 

2337 

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

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

2340 

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

2342 raise OutOfBounds() 

2343 

2344 indis = [] 

2345 weights = [] 

2346 for ia, va in ias: 

2347 iia = ia*nb*nc*ng 

2348 for ib, vb in ibs: 

2349 iib = ib*nc*ng 

2350 for ic, vc in ics: 

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

2352 weights.append(va*vb*vc) 

2353 

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

2355 

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

2357 

2358 xa = (a - amin) / da 

2359 xb = (b - bmin) / db 

2360 xc = (c - cmin) / dc 

2361 

2362 xa_fl = num.floor(xa) 

2363 xa_ce = num.ceil(xa) 

2364 xb_fl = num.floor(xb) 

2365 xb_ce = num.ceil(xb) 

2366 xc_fl = num.floor(xc) 

2367 xc_ce = num.ceil(xc) 

2368 va_fl = 1.0 - (xa - xa_fl) 

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

2370 vb_fl = 1.0 - (xb - xb_fl) 

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

2372 vc_fl = 1.0 - (xc - xc_fl) 

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

2374 

2375 ia_fl = xa_fl.astype(int) 

2376 ia_ce = xa_ce.astype(int) 

2377 ib_fl = xb_fl.astype(int) 

2378 ib_ce = xb_ce.astype(int) 

2379 ic_fl = xc_fl.astype(int) 

2380 ic_ce = xc_ce.astype(int) 

2381 

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

2383 raise OutOfBounds() 

2384 

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

2386 raise OutOfBounds() 

2387 

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

2389 raise OutOfBounds() 

2390 

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

2392 raise OutOfBounds() 

2393 

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

2395 raise OutOfBounds() 

2396 

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

2398 raise OutOfBounds() 

2399 

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

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

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

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

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

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

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

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

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

2409 

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

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

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

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

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

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

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

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

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

2419 

2420 return irecords, weights 

2421 

2422 self._index_function = index_function 

2423 self._indices_function = indices_function 

2424 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2425 self._vicinity_function = vicinity_function 

2426 self._vicinities_function = vicinities_function 

2427 

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

2429 nc = icomponents.size 

2430 dists = source.distances_to(receiver) 

2431 n = dists.size 

2432 receiver_depths = num.empty(nc) 

2433 receiver_depths.fill(receiver.depth) 

2434 return (receiver_depths, 

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

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

2437 icomponents) 

2438 

2439 def make_indexing_args1(self, source, receiver): 

2440 return (receiver.depth, 

2441 source.depth, 

2442 source.distance_to(receiver)) 

2443 

2444 @property 

2445 def short_extent(self): 

2446 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2447 self.receiver_depth_min/km, 

2448 self.receiver_depth_max/km, 

2449 self.receiver_depth_delta/km, 

2450 self.source_depth_min/km, 

2451 self.source_depth_max/km, 

2452 self.source_depth_delta/km, 

2453 self.distance_min/km, 

2454 self.distance_max/km, 

2455 self.distance_delta/km) 

2456 

2457 def fix_ttt_holes(self, sptree, mode): 

2458 from pyrocko import eikonal_ext, spit 

2459 

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

2461 

2462 delta = self.deltas[-1] 

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

2464 

2465 nreceivers, nsources, ndistances = self.ns 

2466 

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

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

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

2470 

2471 speeds = self.get_material_property( 

2472 0., 0., points, 

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

2474 interpolation='multilinear') 

2475 

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

2477 

2478 receiver_times = [] 

2479 for ireceiver in range(nreceivers): 

2480 nodes = num.hstack([ 

2481 num_full( 

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

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

2484 nodes_sr]) 

2485 

2486 times = sptree.interpolate_many(nodes) 

2487 

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

2489 

2490 times = times.reshape(speeds.shape) 

2491 

2492 try: 

2493 eikonal_ext.eikonal_solver_fmm_cartesian( 

2494 speeds, times, delta) 

2495 except eikonal_ext.EikonalExtError as e: 

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

2497 logger.debug( 

2498 'Got a warning from eikonal solver ' 

2499 '- may be ok...') 

2500 else: 

2501 raise 

2502 

2503 receiver_times.append(times) 

2504 

2505 def func(x): 

2506 ias, ibs, ics = \ 

2507 self.grid_interpolation_coefficients(*x) 

2508 

2509 t = 0 

2510 for ia, va in ias: 

2511 times = receiver_times[ia] 

2512 for ib, vb in ibs: 

2513 for ic, vc in ics: 

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

2515 

2516 return t 

2517 

2518 return spit.SPTree( 

2519 f=func, 

2520 ftol=sptree.ftol, 

2521 xbounds=sptree.xbounds, 

2522 xtols=sptree.xtols) 

2523 

2524 

2525class ConfigTypeC(Config): 

2526 ''' 

2527 No symmetrical constraints, one fixed receiver position. 

2528 

2529 * Cartesian 3D source volume around a reference point 

2530 

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

2532 source_east_shift, source_north_shift, component)`` 

2533 ''' 

2534 

2535 receiver = Receiver.T( 

2536 help='Receiver location') 

2537 

2538 source_origin = Location.T( 

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

2540 

2541 source_depth_min = Float.T( 

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

2543 

2544 source_depth_max = Float.T( 

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

2546 

2547 source_depth_delta = Float.T( 

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

2549 

2550 source_east_shift_min = Float.T( 

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

2552 

2553 source_east_shift_max = Float.T( 

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

2555 

2556 source_east_shift_delta = Float.T( 

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

2558 

2559 source_north_shift_min = Float.T( 

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

2561 

2562 source_north_shift_max = Float.T( 

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

2564 

2565 source_north_shift_delta = Float.T( 

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

2567 

2568 short_type = 'C' 

2569 

2570 provided_schemes = ['elastic18'] 

2571 

2572 def get_surface_distance(self, args): 

2573 _, source_east_shift, source_north_shift, _ = args 

2574 sorig = self.source_origin 

2575 sloc = Location( 

2576 lat=sorig.lat, 

2577 lon=sorig.lon, 

2578 north_shift=sorig.north_shift + source_north_shift, 

2579 east_shift=sorig.east_shift + source_east_shift) 

2580 

2581 return self.receiver.distance_to(sloc) 

2582 

2583 def get_distance(self, args): 

2584 sdepth, source_east_shift, source_north_shift, _ = args 

2585 sorig = self.source_origin 

2586 sloc = Location( 

2587 lat=sorig.lat, 

2588 lon=sorig.lon, 

2589 north_shift=sorig.north_shift + source_north_shift, 

2590 east_shift=sorig.east_shift + source_east_shift) 

2591 

2592 return self.receiver.distance_3d_to(sloc) 

2593 

2594 def get_source_depth(self, args): 

2595 return args[0] 

2596 

2597 def get_receiver_depth(self, args): 

2598 return self.receiver.depth 

2599 

2600 def get_source_depths(self): 

2601 return self.coords[0] 

2602 

2603 def _update(self): 

2604 self.mins = num.array([ 

2605 self.source_depth_min, 

2606 self.source_east_shift_min, 

2607 self.source_north_shift_min], 

2608 dtype=float) 

2609 

2610 self.maxs = num.array([ 

2611 self.source_depth_max, 

2612 self.source_east_shift_max, 

2613 self.source_north_shift_max], 

2614 dtype=float) 

2615 

2616 self.deltas = num.array([ 

2617 self.source_depth_delta, 

2618 self.source_east_shift_delta, 

2619 self.source_north_shift_delta], 

2620 dtype=float) 

2621 

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

2623 vicinity_eps).astype(int) + 1 

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

2625 self.deltat = 1.0/self.sample_rate 

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

2627 

2628 self.coords = tuple( 

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

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

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

2632 

2633 def _make_index_functions(self): 

2634 

2635 amin, bmin, cmin = self.mins 

2636 da, db, dc = self.deltas 

2637 na, nb, nc = self.ns 

2638 ng = self.ncomponents 

2639 

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

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

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

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

2644 try: 

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

2646 (na, nb, nc, ng)) 

2647 except ValueError: 

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

2649 

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

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

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

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

2654 

2655 try: 

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

2657 (na, nb, nc, ng)) 

2658 except ValueError: 

2659 raise OutOfBounds() 

2660 

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

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

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

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

2665 

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

2667 raise OutOfBounds() 

2668 

2669 indis = [] 

2670 weights = [] 

2671 for ia, va in ias: 

2672 iia = ia*nb*nc*ng 

2673 for ib, vb in ibs: 

2674 iib = ib*nc*ng 

2675 for ic, vc in ics: 

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

2677 weights.append(va*vb*vc) 

2678 

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

2680 

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

2682 

2683 xa = (a-amin) / da 

2684 xb = (b-bmin) / db 

2685 xc = (c-cmin) / dc 

2686 

2687 xa_fl = num.floor(xa) 

2688 xa_ce = num.ceil(xa) 

2689 xb_fl = num.floor(xb) 

2690 xb_ce = num.ceil(xb) 

2691 xc_fl = num.floor(xc) 

2692 xc_ce = num.ceil(xc) 

2693 va_fl = 1.0 - (xa - xa_fl) 

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

2695 vb_fl = 1.0 - (xb - xb_fl) 

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

2697 vc_fl = 1.0 - (xc - xc_fl) 

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

2699 

2700 ia_fl = xa_fl.astype(int) 

2701 ia_ce = xa_ce.astype(int) 

2702 ib_fl = xb_fl.astype(int) 

2703 ib_ce = xb_ce.astype(int) 

2704 ic_fl = xc_fl.astype(int) 

2705 ic_ce = xc_ce.astype(int) 

2706 

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

2708 raise OutOfBounds() 

2709 

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

2711 raise OutOfBounds() 

2712 

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

2714 raise OutOfBounds() 

2715 

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

2717 raise OutOfBounds() 

2718 

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

2720 raise OutOfBounds() 

2721 

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

2723 raise OutOfBounds() 

2724 

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

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

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

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

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

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

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

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

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

2734 

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

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

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

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

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

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

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

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

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

2744 

2745 return irecords, weights 

2746 

2747 self._index_function = index_function 

2748 self._indices_function = indices_function 

2749 self._vicinity_function = vicinity_function 

2750 self._vicinities_function = vicinities_function 

2751 

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

2753 nc = icomponents.size 

2754 

2755 dists = source.distances_to(self.source_origin) 

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

2757 

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

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

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

2761 

2762 n = dists.size 

2763 

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

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

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

2767 icomponents) 

2768 

2769 def make_indexing_args1(self, source, receiver): 

2770 dist = source.distance_to(self.source_origin) 

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

2772 

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

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

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

2776 

2777 return (source_depth, 

2778 source_east_shift, 

2779 source_north_shift) 

2780 

2781 @property 

2782 def short_extent(self): 

2783 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2784 self.source_depth_min/km, 

2785 self.source_depth_max/km, 

2786 self.source_depth_delta/km, 

2787 self.source_east_shift_min/km, 

2788 self.source_east_shift_max/km, 

2789 self.source_east_shift_delta/km, 

2790 self.source_north_shift_min/km, 

2791 self.source_north_shift_max/km, 

2792 self.source_north_shift_delta/km) 

2793 

2794 

2795class Weighting(Object): 

2796 factor = Float.T(default=1.0) 

2797 

2798 

2799class Taper(Object): 

2800 tmin = Timing.T() 

2801 tmax = Timing.T() 

2802 tfade = Float.T(default=0.0) 

2803 shape = StringChoice.T( 

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

2805 default='cos', 

2806 optional=True) 

2807 

2808 

2809class SimplePattern(SObject): 

2810 

2811 _pool = {} 

2812 

2813 def __init__(self, pattern): 

2814 self._pattern = pattern 

2815 SObject.__init__(self) 

2816 

2817 def __str__(self): 

2818 return self._pattern 

2819 

2820 @property 

2821 def regex(self): 

2822 pool = SimplePattern._pool 

2823 if self.pattern not in pool: 

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

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

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

2827 

2828 return pool[self.pattern] 

2829 

2830 def match(self, s): 

2831 return self.regex.match(s) 

2832 

2833 

2834class WaveformType(StringChoice): 

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

2836 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2837 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2838 

2839 

2840class ChannelSelection(Object): 

2841 pattern = SimplePattern.T(optional=True) 

2842 min_sample_rate = Float.T(optional=True) 

2843 max_sample_rate = Float.T(optional=True) 

2844 

2845 

2846class StationSelection(Object): 

2847 includes = SimplePattern.T() 

2848 excludes = SimplePattern.T() 

2849 distance_min = Float.T(optional=True) 

2850 distance_max = Float.T(optional=True) 

2851 azimuth_min = Float.T(optional=True) 

2852 azimuth_max = Float.T(optional=True) 

2853 

2854 

2855class WaveformSelection(Object): 

2856 channel_selection = ChannelSelection.T(optional=True) 

2857 station_selection = StationSelection.T(optional=True) 

2858 taper = Taper.T() 

2859 # filter = FrequencyResponse.T() 

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

2861 weighting = Weighting.T(optional=True) 

2862 sample_rate = Float.T(optional=True) 

2863 gf_store_id = StringID.T(optional=True) 

2864 

2865 

2866def indi12(x, n): 

2867 ''' 

2868 Get linear interpolation index and weight. 

2869 ''' 

2870 

2871 r = round(x) 

2872 if abs(r - x) < vicinity_eps: 

2873 i = int(r) 

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

2875 raise OutOfBounds() 

2876 

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

2878 else: 

2879 f = math.floor(x) 

2880 i = int(f) 

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

2882 raise OutOfBounds() 

2883 

2884 v = x-f 

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

2886 

2887 

2888def float_or_none(s): 

2889 units = { 

2890 'k': 1e3, 

2891 'M': 1e6, 

2892 } 

2893 

2894 factor = 1.0 

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

2896 factor = units[s[-1]] 

2897 s = s[:-1] 

2898 if not s: 

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

2900 

2901 if s: 

2902 return float(s) * factor 

2903 else: 

2904 return None 

2905 

2906 

2907class GridSpecError(Exception): 

2908 def __init__(self, s): 

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

2910 

2911 

2912def parse_grid_spec(spec): 

2913 try: 

2914 result = [] 

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

2916 t = dspec.split('@') 

2917 num = start = stop = step = None 

2918 if len(t) == 2: 

2919 num = int(t[1]) 

2920 if num <= 0: 

2921 raise GridSpecError(spec) 

2922 

2923 elif len(t) > 2: 

2924 raise GridSpecError(spec) 

2925 

2926 s = t[0] 

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

2928 if len(v) == 1: 

2929 start = stop = v[0] 

2930 if len(v) >= 2: 

2931 start, stop = v[0:2] 

2932 if len(v) == 3: 

2933 step = v[2] 

2934 

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

2936 raise GridSpecError(spec) 

2937 

2938 if step == 0.0: 

2939 raise GridSpecError(spec) 

2940 

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

2942 

2943 except ValueError: 

2944 raise GridSpecError(spec) 

2945 

2946 return result 

2947 

2948 

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

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

2951 if start is None: 

2952 start = ma if swap else mi 

2953 if stop is None: 

2954 stop = mi if swap else ma 

2955 if step is None: 

2956 step = -inc if ma < mi else inc 

2957 if num is None: 

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

2959 raise GridSpecError() 

2960 

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

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

2963 if abs(stop-stop2) > eps: 

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

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

2966 else: 

2967 stop = stop2 

2968 

2969 if start == stop: 

2970 num = 1 

2971 

2972 return start, stop, num 

2973 

2974 

2975def nditer_outer(x): 

2976 return num.nditer( 

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

2978 

2979 

2980def nodes(xs): 

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

2982 nnodes = num.prod(ns) 

2983 ndim = len(xs) 

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

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

2986 x = xs[idim] 

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

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

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

2990 

2991 return nodes 

2992 

2993 

2994def filledi(x, n): 

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

2996 a.fill(x) 

2997 return a 

2998 

2999 

3000config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3001 

3002discretized_source_classes = [ 

3003 DiscretizedExplosionSource, 

3004 DiscretizedSFSource, 

3005 DiscretizedMTSource, 

3006 DiscretizedPorePressureSource] 

3007 

3008 

3009__all__ = ''' 

3010Earthmodel1D 

3011StringID 

3012ScopeType 

3013WaveformType 

3014QuantityType 

3015NearfieldTermsType 

3016Reference 

3017Region 

3018CircularRegion 

3019RectangularRegion 

3020PhaseSelect 

3021InvalidTimingSpecification 

3022Timing 

3023TPDef 

3024OutOfBounds 

3025Location 

3026Receiver 

3027'''.split() + [ 

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

3029ComponentScheme 

3030component_scheme_to_description 

3031component_schemes 

3032Config 

3033GridSpecError 

3034Weighting 

3035Taper 

3036SimplePattern 

3037WaveformType 

3038ChannelSelection 

3039StationSelection 

3040WaveformSelection 

3041nditer_outer 

3042dump 

3043load 

3044discretized_source_classes 

3045config_type_classes 

3046UnavailableScheme 

3047InterpolationMethod 

3048SeismosizerTrace 

3049SeismosizerResult 

3050Result 

3051StaticResult 

3052'''.split()