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 TimingAttributeError(Exception): 

405 pass 

406 

407 

408class Timing(SObject): 

409 ''' 

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

411 

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

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

414 arrival or group of such arrivals. 

415 

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

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

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

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

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

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

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

423 

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

425 

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

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

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

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

430 velocity [km/s] 

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

432 

433 **Examples:** 

434 

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

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

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

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

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

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

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

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

443 undefined for a given geometry 

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

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

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

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

448 arrivals A, B, and C 

449 ''' 

450 

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

452 

453 if s is not None: 

454 offset_is = None 

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

456 try: 

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

458 

459 if s.endswith('S'): 

460 offset_is = 'slowness' 

461 

462 phase_defs = [] 

463 select = '' 

464 

465 except ValueError: 

466 matched = False 

467 m = timing_regex.match(s) 

468 if m: 

469 if m.group(3): 

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

471 elif m.group(19): 

472 phase_defs = [m.group(19)] 

473 else: 

474 phase_defs = [] 

475 

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

477 

478 offset = 0.0 

479 soff = m.group(27) 

480 if soff: 

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

482 if soff.endswith('S'): 

483 offset_is = 'slowness' 

484 elif soff.endswith('%'): 

485 offset_is = 'percent' 

486 

487 matched = True 

488 

489 else: 

490 m = timing_regex_old.match(s) 

491 if m: 

492 if m.group(3): 

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

494 elif m.group(5): 

495 phase_defs = [m.group(5)] 

496 else: 

497 phase_defs = [] 

498 

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

500 

501 offset = 0.0 

502 if m.group(6): 

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

504 

505 phase_defs = [ 

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

507 

508 matched = True 

509 

510 if not matched: 

511 raise InvalidTimingSpecification(s) 

512 

513 kwargs = dict( 

514 phase_defs=phase_defs, 

515 select=select, 

516 offset=offset, 

517 offset_is=offset_is) 

518 

519 SObject.__init__(self, **kwargs) 

520 

521 def __str__(self): 

522 s = [] 

523 if self.phase_defs: 

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

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

526 sphases = '{%s}' % sphases 

527 

528 if self.select: 

529 sphases = self.select + sphases 

530 

531 s.append(sphases) 

532 

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

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

535 if self.offset_is == 'slowness': 

536 s.append('S') 

537 elif self.offset_is == 'percent': 

538 s.append('%') 

539 

540 return ''.join(s) 

541 

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

543 try: 

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

545 phase_offset = get_phase( 

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

547 offset = phase_offset(args) 

548 else: 

549 offset = self.offset 

550 

551 if self.phase_defs: 

552 extra_args = () 

553 if attributes: 

554 extra_args = (attributes,) 

555 

556 phases = [ 

557 get_phase(phase_def, *extra_args) 

558 for phase_def in self.phase_defs] 

559 

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

561 results = [ 

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

563 for result in results] 

564 

565 results = [ 

566 result for result in results 

567 if result[0] is not None] 

568 

569 if self.select == 'first': 

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

571 elif self.select == 'last': 

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

573 

574 if not results: 

575 if attributes is not None: 

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

577 else: 

578 return None 

579 

580 else: 

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

582 if self.offset_is == 'percent': 

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

584 else: 

585 t = t + offset 

586 

587 if attributes is not None: 

588 return (t, ) + values 

589 else: 

590 return t 

591 

592 else: 

593 if attributes is not None: 

594 raise TimingAttributeError( 

595 'Cannot get phase attributes from offset-only ' 

596 'definition.') 

597 return offset 

598 

599 except spit.OutOfBounds: 

600 raise OutOfBounds(args) 

601 

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

603 offset = Float.T(default=0.0) 

604 offset_is = String.T(optional=True) 

605 select = PhaseSelect.T( 

606 default='', 

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

608 tuple(PhaseSelect.choices))) 

609 

610 

611def mkdefs(s): 

612 defs = [] 

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

614 try: 

615 defs.append(float(x)) 

616 except ValueError: 

617 if x.startswith('!'): 

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

619 else: 

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

621 

622 return defs 

623 

624 

625class TPDef(Object): 

626 ''' 

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

628 ''' 

629 

630 id = StringID.T( 

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

632 definition = String.T( 

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

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

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

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

637 

638 @property 

639 def phases(self): 

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

641 if isinstance(x, cake.PhaseDef)] 

642 

643 @property 

644 def horizontal_velocities(self): 

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

646 

647 

648class OutOfBounds(Exception): 

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

650 Exception.__init__(self) 

651 self.values = values 

652 self.reason = reason 

653 self.context = None 

654 

655 def __str__(self): 

656 scontext = '' 

657 if self.context: 

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

659 

660 if self.reason: 

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

662 

663 if self.values: 

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

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

666 else: 

667 return 'out of bounds%s' % scontext 

668 

669 

670class MultiLocation(Object): 

671 

672 lats = Array.T( 

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

674 help='Latitudes of targets.') 

675 lons = Array.T( 

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

677 help='Longitude of targets.') 

678 north_shifts = Array.T( 

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

680 help='North shifts of targets.') 

681 east_shifts = Array.T( 

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

683 help='East shifts of targets.') 

684 elevation = Array.T( 

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

686 help='Elevations of targets.') 

687 

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

689 self._coords5 = None 

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

691 

692 @property 

693 def coords5(self): 

694 if self._coords5 is not None: 

695 return self._coords5 

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

697 self.elevation] 

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

699 if not sizes: 

700 raise AttributeError('Empty StaticTarget') 

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

702 raise AttributeError('Inconsistent coordinate sizes.') 

703 

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

705 for idx, p in enumerate(props): 

706 if p is not None: 

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

708 

709 return self._coords5 

710 

711 @property 

712 def ncoords(self): 

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

714 return 0 

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

716 

717 def get_latlon(self): 

718 ''' 

719 Get all coordinates as lat/lon. 

720 

721 :returns: Coordinates in Latitude, Longitude 

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

723 ''' 

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

725 for i in range(self.ncoords): 

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

727 return latlons 

728 

729 def get_corner_coords(self): 

730 ''' 

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

732 

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

734 :rtype: tuple 

735 ''' 

736 latlon = self.get_latlon() 

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

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

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

740 

741 

742class Receiver(Location): 

743 codes = Tuple.T( 

744 3, String.T(), 

745 optional=True, 

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

747 

748 def pyrocko_station(self): 

749 from pyrocko import model 

750 

751 lat, lon = self.effective_latlon 

752 

753 return model.Station( 

754 network=self.codes[0], 

755 station=self.codes[1], 

756 location=self.codes[2], 

757 lat=lat, 

758 lon=lon, 

759 depth=self.depth) 

760 

761 

762def g(x, d): 

763 if x is None: 

764 return d 

765 else: 

766 return x 

767 

768 

769class UnavailableScheme(Exception): 

770 pass 

771 

772 

773class InvalidNComponents(Exception): 

774 pass 

775 

776 

777class DiscretizedSource(Object): 

778 ''' 

779 Base class for discretized sources. 

780 

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

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

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

784 source distributions are represented by subclasses of the 

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

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

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

788 explosion/implosion type source distributions. The geometry calculations 

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

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

791 

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

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

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

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

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

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

798 ''' 

799 

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

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

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

803 lat = Float.T(optional=True) 

804 lon = Float.T(optional=True) 

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

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

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

808 dl = Float.T(optional=True) 

809 dw = Float.T(optional=True) 

810 nl = Float.T(optional=True) 

811 nw = Float.T(optional=True) 

812 

813 @classmethod 

814 def check_scheme(cls, scheme): 

815 ''' 

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

817 

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

819 supported by this discretized source class. 

820 ''' 

821 

822 if scheme not in cls.provided_schemes: 

823 raise UnavailableScheme( 

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

825 (cls.__name__, scheme)) 

826 

827 def __init__(self, **kwargs): 

828 Object.__init__(self, **kwargs) 

829 self._latlons = None 

830 

831 def __setattr__(self, name, value): 

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

833 'lats', 'lons'): 

834 self.__dict__['_latlons'] = None 

835 

836 Object.__setattr__(self, name, value) 

837 

838 def get_source_terms(self, scheme): 

839 raise NotImplementedError() 

840 

841 def make_weights(self, receiver, scheme): 

842 raise NotImplementedError() 

843 

844 @property 

845 def effective_latlons(self): 

846 ''' 

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

848 ''' 

849 

850 if self._latlons is None: 

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

852 if (self.north_shifts is not None and 

853 self.east_shifts is not None): 

854 self._latlons = orthodrome.ne_to_latlon( 

855 self.lats, self.lons, 

856 self.north_shifts, self.east_shifts) 

857 else: 

858 self._latlons = self.lats, self.lons 

859 else: 

860 lat = g(self.lat, 0.0) 

861 lon = g(self.lon, 0.0) 

862 self._latlons = orthodrome.ne_to_latlon( 

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

864 

865 return self._latlons 

866 

867 @property 

868 def effective_north_shifts(self): 

869 if self.north_shifts is not None: 

870 return self.north_shifts 

871 else: 

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

873 

874 @property 

875 def effective_east_shifts(self): 

876 if self.east_shifts is not None: 

877 return self.east_shifts 

878 else: 

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

880 

881 def same_origin(self, receiver): 

882 ''' 

883 Check if receiver has same reference point. 

884 ''' 

885 

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

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

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

889 

890 def azibazis_to(self, receiver): 

891 ''' 

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

893 points. 

894 ''' 

895 

896 if self.same_origin(receiver): 

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

898 receiver.north_shift - self.north_shifts) 

899 bazis = azis + 180. 

900 else: 

901 slats, slons = self.effective_latlons 

902 rlat, rlon = receiver.effective_latlon 

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

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

905 

906 return azis, bazis 

907 

908 def distances_to(self, receiver): 

909 ''' 

910 Compute distances to receiver for all contained points. 

911 ''' 

912 if self.same_origin(receiver): 

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

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

915 

916 else: 

917 slats, slons = self.effective_latlons 

918 rlat, rlon = receiver.effective_latlon 

919 return orthodrome.distance_accurate50m_numpy(slats, slons, 

920 rlat, rlon) 

921 

922 def element_coords(self, i): 

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

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

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

926 else: 

927 lat = self.lat 

928 lon = self.lon 

929 

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

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

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

933 

934 else: 

935 north_shift = east_shift = 0.0 

936 

937 return lat, lon, north_shift, east_shift 

938 

939 def coords5(self): 

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

941 

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

943 xs[:, 0] = self.lats 

944 xs[:, 1] = self.lons 

945 else: 

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

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

948 

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

950 xs[:, 2] = self.north_shifts 

951 xs[:, 3] = self.east_shifts 

952 

953 xs[:, 4] = self.depths 

954 

955 return xs 

956 

957 @property 

958 def nelements(self): 

959 return self.times.size 

960 

961 @classmethod 

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

963 ''' 

964 Combine several discretized source models. 

965 

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

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

968 factors and reference times of the parameterized (undiscretized) 

969 sources match or are accounted for. 

970 ''' 

971 

972 first = sources[0] 

973 

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

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

976 'sources of same type.') 

977 

978 latlons = [] 

979 for s in sources: 

980 latlons.append(s.effective_latlons) 

981 

982 lats, lons = num.hstack(latlons) 

983 

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

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

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

987 same_ref = num.all( 

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

989 else: 

990 same_ref = False 

991 

992 cat = num.concatenate 

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

994 print(times) 

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

996 

997 if same_ref: 

998 lat = first.lat 

999 lon = first.lon 

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

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

1002 lats = None 

1003 lons = None 

1004 else: 

1005 lat = None 

1006 lon = None 

1007 north_shifts = None 

1008 east_shifts = None 

1009 

1010 return cls( 

1011 times=times, lat=lat, lon=lon, lats=lats, lons=lons, 

1012 north_shifts=north_shifts, east_shifts=east_shifts, 

1013 depths=depths, **kwargs) 

1014 

1015 def centroid_position(self): 

1016 moments = self.moments() 

1017 norm = num.sum(moments) 

1018 if norm != 0.0: 

1019 w = moments / num.sum(moments) 

1020 else: 

1021 w = num.ones(moments.size) 

1022 

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

1024 lats, lons = self.effective_latlons 

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

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

1027 else: 

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

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

1030 

1031 cn = num.sum(n*w) 

1032 ce = num.sum(e*w) 

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

1034 

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

1036 lat = clat 

1037 lon = clon 

1038 north_shift = 0. 

1039 east_shift = 0. 

1040 else: 

1041 lat = g(self.lat, 0.0) 

1042 lon = g(self.lon, 0.0) 

1043 north_shift = cn 

1044 east_shift = ce 

1045 

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

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

1048 return tuple(float(x) for x in 

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

1050 

1051 

1052class DiscretizedExplosionSource(DiscretizedSource): 

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

1054 

1055 provided_schemes = ( 

1056 'elastic2', 

1057 'elastic8', 

1058 'elastic10', 

1059 ) 

1060 

1061 def get_source_terms(self, scheme): 

1062 self.check_scheme(scheme) 

1063 

1064 if scheme == 'elastic2': 

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

1066 

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

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

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

1070 return m6s 

1071 else: 

1072 assert False 

1073 

1074 def make_weights(self, receiver, scheme): 

1075 self.check_scheme(scheme) 

1076 

1077 azis, bazis = self.azibazis_to(receiver) 

1078 

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

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

1081 

1082 m0s = self.m0s 

1083 n = azis.size 

1084 

1085 cat = num.concatenate 

1086 rep = num.repeat 

1087 

1088 if scheme == 'elastic2': 

1089 w_n = cb*m0s 

1090 g_n = filledi(0, n) 

1091 w_e = sb*m0s 

1092 g_e = filledi(0, n) 

1093 w_d = m0s 

1094 g_d = filledi(1, n) 

1095 

1096 elif scheme == 'elastic8': 

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

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

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

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

1101 w_d = cat((m0s, m0s)) 

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

1103 

1104 elif scheme == 'elastic10': 

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

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

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

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

1109 w_d = cat((m0s, m0s, m0s)) 

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

1111 

1112 else: 

1113 assert False 

1114 

1115 return ( 

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

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

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

1119 ) 

1120 

1121 def split(self): 

1122 from pyrocko.gf.seismosizer import ExplosionSource 

1123 sources = [] 

1124 for i in range(self.nelements): 

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

1126 sources.append(ExplosionSource( 

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

1128 lat=lat, 

1129 lon=lon, 

1130 north_shift=north_shift, 

1131 east_shift=east_shift, 

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

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

1134 

1135 return sources 

1136 

1137 def moments(self): 

1138 return self.m0s 

1139 

1140 def centroid(self): 

1141 from pyrocko.gf.seismosizer import ExplosionSource 

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

1143 self.centroid_position() 

1144 

1145 return ExplosionSource( 

1146 time=time, 

1147 lat=lat, 

1148 lon=lon, 

1149 north_shift=north_shift, 

1150 east_shift=east_shift, 

1151 depth=depth, 

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

1153 

1154 @classmethod 

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

1156 ''' 

1157 Combine several discretized source models. 

1158 

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

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

1161 factors and reference times of the parameterized (undiscretized) 

1162 sources match or are accounted for. 

1163 ''' 

1164 

1165 if 'm0s' not in kwargs: 

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

1167 

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

1169 **kwargs) 

1170 

1171 

1172class DiscretizedSFSource(DiscretizedSource): 

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

1174 

1175 provided_schemes = ( 

1176 'elastic5', 

1177 ) 

1178 

1179 def get_source_terms(self, scheme): 

1180 self.check_scheme(scheme) 

1181 

1182 return self.forces 

1183 

1184 def make_weights(self, receiver, scheme): 

1185 self.check_scheme(scheme) 

1186 

1187 azis, bazis = self.azibazis_to(receiver) 

1188 

1189 sa = num.sin(azis*d2r) 

1190 ca = num.cos(azis*d2r) 

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

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

1193 

1194 forces = self.forces 

1195 fn = forces[:, 0] 

1196 fe = forces[:, 1] 

1197 fd = forces[:, 2] 

1198 

1199 f0 = fd 

1200 f1 = ca * fn + sa * fe 

1201 f2 = ca * fe - sa * fn 

1202 

1203 n = azis.size 

1204 

1205 if scheme == 'elastic5': 

1206 ioff = 0 

1207 

1208 cat = num.concatenate 

1209 rep = num.repeat 

1210 

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

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

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

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

1215 w_d = cat((f0, f1)) 

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

1217 

1218 return ( 

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

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

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

1222 ) 

1223 

1224 @classmethod 

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

1226 ''' 

1227 Combine several discretized source models. 

1228 

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

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

1231 factors and reference times of the parameterized (undiscretized) 

1232 sources match or are accounted for. 

1233 ''' 

1234 

1235 if 'forces' not in kwargs: 

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

1237 

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

1239 

1240 def moments(self): 

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

1242 

1243 def centroid(self): 

1244 from pyrocko.gf.seismosizer import SFSource 

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

1246 self.centroid_position() 

1247 

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

1249 return SFSource( 

1250 time=time, 

1251 lat=lat, 

1252 lon=lon, 

1253 north_shift=north_shift, 

1254 east_shift=east_shift, 

1255 depth=depth, 

1256 fn=fn, 

1257 fe=fe, 

1258 fd=fd) 

1259 

1260 

1261class DiscretizedMTSource(DiscretizedSource): 

1262 m6s = Array.T( 

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

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

1265 

1266 provided_schemes = ( 

1267 'elastic8', 

1268 'elastic10', 

1269 'elastic18', 

1270 ) 

1271 

1272 def get_source_terms(self, scheme): 

1273 self.check_scheme(scheme) 

1274 return self.m6s 

1275 

1276 def make_weights(self, receiver, scheme): 

1277 self.check_scheme(scheme) 

1278 

1279 m6s = self.m6s 

1280 n = m6s.shape[0] 

1281 rep = num.repeat 

1282 

1283 if scheme == 'elastic18': 

1284 w_n = m6s.flatten() 

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

1286 w_e = m6s.flatten() 

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

1288 w_d = m6s.flatten() 

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

1290 

1291 else: 

1292 azis, bazis = self.azibazis_to(receiver) 

1293 

1294 sa = num.sin(azis*d2r) 

1295 ca = num.cos(azis*d2r) 

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

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

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

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

1300 

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

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

1303 f2 = m6s[:, 2] 

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

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

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

1307 

1308 cat = num.concatenate 

1309 

1310 if scheme == 'elastic8': 

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

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

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

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

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

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

1317 

1318 elif scheme == 'elastic10': 

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

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

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

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

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

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

1325 

1326 else: 

1327 assert False 

1328 

1329 return ( 

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

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

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

1333 ) 

1334 

1335 def split(self): 

1336 from pyrocko.gf.seismosizer import MTSource 

1337 sources = [] 

1338 for i in range(self.nelements): 

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

1340 sources.append(MTSource( 

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

1342 lat=lat, 

1343 lon=lon, 

1344 north_shift=north_shift, 

1345 east_shift=east_shift, 

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

1347 m6=self.m6s[i])) 

1348 

1349 return sources 

1350 

1351 def moments(self): 

1352 moments = num.array( 

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

1354 for m6 in self.m6s]) 

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

1356 

1357 def get_moment_rate(self, deltat=None): 

1358 moments = self.moments() 

1359 times = self.times 

1360 times -= times.min() 

1361 

1362 t_max = times.max() 

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

1364 mom_times[mom_times > t_max] = t_max 

1365 

1366 # Right open histrogram bins 

1367 mom, _ = num.histogram( 

1368 times, 

1369 bins=mom_times, 

1370 weights=moments) 

1371 

1372 deltat = num.diff(mom_times) 

1373 mom_rate = mom / deltat 

1374 

1375 return mom_rate, mom_times[1:] 

1376 

1377 def centroid(self): 

1378 from pyrocko.gf.seismosizer import MTSource 

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

1380 self.centroid_position() 

1381 

1382 return MTSource( 

1383 time=time, 

1384 lat=lat, 

1385 lon=lon, 

1386 north_shift=north_shift, 

1387 east_shift=east_shift, 

1388 depth=depth, 

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

1390 

1391 @classmethod 

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

1393 ''' 

1394 Combine several discretized source models. 

1395 

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

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

1398 factors and reference times of the parameterized (undiscretized) 

1399 sources match or are accounted for. 

1400 ''' 

1401 

1402 if 'm6s' not in kwargs: 

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

1404 

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

1406 

1407 

1408class DiscretizedPorePressureSource(DiscretizedSource): 

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

1410 

1411 provided_schemes = ( 

1412 'poroelastic10', 

1413 ) 

1414 

1415 def get_source_terms(self, scheme): 

1416 self.check_scheme(scheme) 

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

1418 

1419 def make_weights(self, receiver, scheme): 

1420 self.check_scheme(scheme) 

1421 

1422 azis, bazis = self.azibazis_to(receiver) 

1423 

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

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

1426 

1427 pp = self.pp 

1428 n = bazis.size 

1429 

1430 w_un = cb*pp 

1431 g_un = filledi(1, n) 

1432 w_ue = sb*pp 

1433 g_ue = filledi(1, n) 

1434 w_ud = pp 

1435 g_ud = filledi(0, n) 

1436 

1437 w_tn = cb*pp 

1438 g_tn = filledi(6, n) 

1439 w_te = sb*pp 

1440 g_te = filledi(6, n) 

1441 

1442 w_pp = pp 

1443 g_pp = filledi(7, n) 

1444 

1445 w_dvn = cb*pp 

1446 g_dvn = filledi(9, n) 

1447 w_dve = sb*pp 

1448 g_dve = filledi(9, n) 

1449 w_dvd = pp 

1450 g_dvd = filledi(8, n) 

1451 

1452 return ( 

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

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

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

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

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

1458 ('pore_pressure', w_pp, g_pp), 

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

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

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

1462 ) 

1463 

1464 def moments(self): 

1465 return self.pp 

1466 

1467 def centroid(self): 

1468 from pyrocko.gf.seismosizer import PorePressurePointSource 

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

1470 self.centroid_position() 

1471 

1472 return PorePressurePointSource( 

1473 time=time, 

1474 lat=lat, 

1475 lon=lon, 

1476 north_shift=north_shift, 

1477 east_shift=east_shift, 

1478 depth=depth, 

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

1480 

1481 @classmethod 

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

1483 ''' 

1484 Combine several discretized source models. 

1485 

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

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

1488 factors and reference times of the parameterized (undiscretized) 

1489 sources match or are accounted for. 

1490 ''' 

1491 

1492 if 'pp' not in kwargs: 

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

1494 

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

1496 **kwargs) 

1497 

1498 

1499class Region(Object): 

1500 name = String.T(optional=True) 

1501 

1502 

1503class RectangularRegion(Region): 

1504 lat_min = Float.T() 

1505 lat_max = Float.T() 

1506 lon_min = Float.T() 

1507 lon_max = Float.T() 

1508 

1509 

1510class CircularRegion(Region): 

1511 lat = Float.T() 

1512 lon = Float.T() 

1513 radius = Float.T() 

1514 

1515 

1516class Config(Object): 

1517 ''' 

1518 Green's function store meta information. 

1519 

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

1521 configuration types are: 

1522 

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

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

1525 

1526 * Problem is invariant to horizontal translations and rotations around 

1527 vertical axis. 

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

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

1530 component)`` 

1531 

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

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

1534 

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

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

1537 receiver_depth, component)`` 

1538 

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

1540 constraints but fixed receiver positions 

1541 

1542 * Cartesian source volume around a reference point 

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

1544 source_east_shift, source_north_shift, component)`` 

1545 ''' 

1546 

1547 id = StringID.T( 

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

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

1550 'letter.') 

1551 

1552 derived_from_id = StringID.T( 

1553 optional=True, 

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

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

1556 

1557 version = String.T( 

1558 default='1.0', 

1559 optional=True, 

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

1561 

1562 modelling_code_id = StringID.T( 

1563 optional=True, 

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

1565 

1566 author = Unicode.T( 

1567 optional=True, 

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

1569 

1570 author_email = String.T( 

1571 optional=True, 

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

1573 

1574 created_time = Timestamp.T( 

1575 optional=True, 

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

1577 

1578 regions = List.T( 

1579 Region.T(), 

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

1581 

1582 scope_type = ScopeType.T( 

1583 optional=True, 

1584 help='Distance range scope of the store ' 

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

1586 

1587 waveform_type = WaveType.T( 

1588 optional=True, 

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

1590 

1591 nearfield_terms = NearfieldTermsType.T( 

1592 optional=True, 

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

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

1595 

1596 description = String.T( 

1597 optional=True, 

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

1599 

1600 references = List.T( 

1601 Reference.T(), 

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

1603 'related work.') 

1604 

1605 earthmodel_1d = Earthmodel1D.T( 

1606 optional=True, 

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

1608 

1609 earthmodel_receiver_1d = Earthmodel1D.T( 

1610 optional=True, 

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

1612 

1613 can_interpolate_source = Bool.T( 

1614 optional=True, 

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

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

1617 

1618 can_interpolate_receiver = Bool.T( 

1619 optional=True, 

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

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

1622 

1623 frequency_min = Float.T( 

1624 optional=True, 

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

1626 

1627 frequency_max = Float.T( 

1628 optional=True, 

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

1630 

1631 sample_rate = Float.T( 

1632 optional=True, 

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

1634 

1635 factor = Float.T( 

1636 default=1.0, 

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

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

1639 optional=True) 

1640 

1641 component_scheme = ComponentScheme.T( 

1642 default='elastic10', 

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

1644 

1645 stored_quantity = QuantityType.T( 

1646 optional=True, 

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

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

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

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

1651 

1652 tabulated_phases = List.T( 

1653 TPDef.T(), 

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

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

1656 

1657 ncomponents = Int.T( 

1658 optional=True, 

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

1660 

1661 uuid = String.T( 

1662 optional=True, 

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

1664 'GF store for practical purposes.') 

1665 

1666 reference = String.T( 

1667 optional=True, 

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

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

1670 

1671 def __init__(self, **kwargs): 

1672 self._do_auto_updates = False 

1673 Object.__init__(self, **kwargs) 

1674 self._index_function = None 

1675 self._indices_function = None 

1676 self._vicinity_function = None 

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

1678 self._do_auto_updates = True 

1679 self.update() 

1680 

1681 def check_ncomponents(self): 

1682 ncomponents = component_scheme_to_description[ 

1683 self.component_scheme].ncomponents 

1684 

1685 if self.ncomponents is None: 

1686 self.ncomponents = ncomponents 

1687 elif ncomponents != self.ncomponents: 

1688 raise InvalidNComponents( 

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

1690 self.ncomponents, self.component_scheme)) 

1691 

1692 def __setattr__(self, name, value): 

1693 Object.__setattr__(self, name, value) 

1694 try: 

1695 self.T.get_property(name) 

1696 if self._do_auto_updates: 

1697 self.update() 

1698 

1699 except ValueError: 

1700 pass 

1701 

1702 def update(self): 

1703 self.check_ncomponents() 

1704 self._update() 

1705 self._make_index_functions() 

1706 

1707 def irecord(self, *args): 

1708 return self._index_function(*args) 

1709 

1710 def irecords(self, *args): 

1711 return self._indices_function(*args) 

1712 

1713 def vicinity(self, *args): 

1714 return self._vicinity_function(*args) 

1715 

1716 def vicinities(self, *args): 

1717 return self._vicinities_function(*args) 

1718 

1719 def grid_interpolation_coefficients(self, *args): 

1720 return self._grid_interpolation_coefficients(*args) 

1721 

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

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

1724 

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

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

1727 

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

1729 i = 0 

1730 arrs = [] 

1731 ntotal = 1 

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

1733 if gdef and len(gdef) > i: 

1734 sssn = gdef[i] 

1735 else: 

1736 sssn = (None,)*4 

1737 

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

1739 ntotal *= len(arr) 

1740 

1741 arrs.append(arr) 

1742 i += 1 

1743 

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

1745 return nditer_outer(arrs[:level]) 

1746 

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

1748 nthreads=0): 

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

1750 

1751 out = [] 

1752 delays = source.times 

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

1754 receiver, 

1755 self.component_scheme): 

1756 

1757 weights *= self.factor 

1758 

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

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

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

1762 

1763 return out 

1764 

1765 def short_info(self): 

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

1767 

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

1769 interpolation=None): 

1770 ''' 

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

1772 

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

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

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

1776 ``(lat, lon)`` 

1777 :param interpolation: interpolation method. Choose from 

1778 ``('nearest_neighbor', 'multilinear')`` 

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

1780 point 

1781 

1782 The default implementation retrieves and interpolates the shear moduli 

1783 from the contained 1D velocity profile. 

1784 ''' 

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

1786 parameter='shear_moduli', 

1787 interpolation=interpolation) 

1788 

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

1790 interpolation=None): 

1791 ''' 

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

1793 

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

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

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

1797 ``(lat, lon)`` 

1798 :param interpolation: interpolation method. Choose from 

1799 ``('nearest_neighbor', 'multilinear')`` 

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

1801 point 

1802 

1803 The default implementation retrieves and interpolates the lambda moduli 

1804 from the contained 1D velocity profile. 

1805 ''' 

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

1807 parameter='lambda_moduli', 

1808 interpolation=interpolation) 

1809 

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

1811 interpolation=None): 

1812 ''' 

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

1814 

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

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

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

1818 ``(lat, lon)`` 

1819 :param interpolation: interpolation method. Choose from 

1820 ``('nearest_neighbor', 'multilinear')`` 

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

1822 point 

1823 

1824 The default implementation retrieves and interpolates the lambda moduli 

1825 from the contained 1D velocity profile. 

1826 ''' 

1827 lambda_moduli = self.get_material_property( 

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

1829 interpolation=interpolation) 

1830 shear_moduli = self.get_material_property( 

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

1832 interpolation=interpolation) 

1833 return lambda_moduli + (2 / 3) * shear_moduli 

1834 

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

1836 ''' 

1837 Get Vs at given points from contained velocity model. 

1838 

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

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

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

1842 ``(lat, lon)`` 

1843 :param interpolation: interpolation method. Choose from 

1844 ``('nearest_neighbor', 'multilinear')`` 

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

1846 point 

1847 

1848 The default implementation retrieves and interpolates Vs 

1849 from the contained 1D velocity profile. 

1850 ''' 

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

1852 parameter='vs', 

1853 interpolation=interpolation) 

1854 

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

1856 ''' 

1857 Get Vp at given points from contained velocity model. 

1858 

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

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

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

1862 ``(lat, lon)`` 

1863 :param interpolation: interpolation method. Choose from 

1864 ``('nearest_neighbor', 'multilinear')`` 

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

1866 point 

1867 

1868 The default implementation retrieves and interpolates Vp 

1869 from the contained 1D velocity profile. 

1870 ''' 

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

1872 parameter='vp', 

1873 interpolation=interpolation) 

1874 

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

1876 ''' 

1877 Get rho at given points from contained velocity model. 

1878 

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

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

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

1882 ``(lat, lon)`` 

1883 :param interpolation: interpolation method. Choose from 

1884 ``('nearest_neighbor', 'multilinear')`` 

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

1886 point 

1887 

1888 The default implementation retrieves and interpolates rho 

1889 from the contained 1D velocity profile. 

1890 ''' 

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

1892 parameter='rho', 

1893 interpolation=interpolation) 

1894 

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

1896 interpolation=None): 

1897 

1898 if interpolation is None: 

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

1900 'multilinear', 'nearest_neighbor') 

1901 

1902 earthmod = self.earthmodel_1d 

1903 store_depth_profile = self.get_source_depths() 

1904 z_profile = earthmod.profile('z') 

1905 

1906 if parameter == 'vs': 

1907 vs_profile = earthmod.profile('vs') 

1908 profile = num.interp( 

1909 store_depth_profile, z_profile, vs_profile) 

1910 

1911 elif parameter == 'vp': 

1912 vp_profile = earthmod.profile('vp') 

1913 profile = num.interp( 

1914 store_depth_profile, z_profile, vp_profile) 

1915 

1916 elif parameter == 'rho': 

1917 rho_profile = earthmod.profile('rho') 

1918 

1919 profile = num.interp( 

1920 store_depth_profile, z_profile, rho_profile) 

1921 

1922 elif parameter == 'shear_moduli': 

1923 vs_profile = earthmod.profile('vs') 

1924 rho_profile = earthmod.profile('rho') 

1925 

1926 store_vs_profile = num.interp( 

1927 store_depth_profile, z_profile, vs_profile) 

1928 store_rho_profile = num.interp( 

1929 store_depth_profile, z_profile, rho_profile) 

1930 

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

1932 

1933 elif parameter == 'lambda_moduli': 

1934 vs_profile = earthmod.profile('vs') 

1935 vp_profile = earthmod.profile('vp') 

1936 rho_profile = earthmod.profile('rho') 

1937 

1938 store_vs_profile = num.interp( 

1939 store_depth_profile, z_profile, vs_profile) 

1940 store_vp_profile = num.interp( 

1941 store_depth_profile, z_profile, vp_profile) 

1942 store_rho_profile = num.interp( 

1943 store_depth_profile, z_profile, rho_profile) 

1944 

1945 profile = store_rho_profile * ( 

1946 num.power(store_vp_profile, 2) - 

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

1948 else: 

1949 raise TypeError( 

1950 'parameter %s not available' % parameter) 

1951 

1952 if interpolation == 'multilinear': 

1953 kind = 'linear' 

1954 elif interpolation == 'nearest_neighbor': 

1955 kind = 'nearest' 

1956 else: 

1957 raise TypeError( 

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

1959 

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

1961 

1962 try: 

1963 return interpolator(points[:, 2]) 

1964 except ValueError: 

1965 raise OutOfBounds() 

1966 

1967 def is_static(self): 

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

1969 if self.modelling_code_id.startswith(code): 

1970 return True 

1971 return False 

1972 

1973 def is_dynamic(self): 

1974 return not self.is_static() 

1975 

1976 def get_source_depths(self): 

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

1978 

1979 def get_tabulated_phase(self, phase_id): 

1980 ''' 

1981 Get tabulated phase definition. 

1982 ''' 

1983 

1984 for pdef in self.tabulated_phases: 

1985 if pdef.id == phase_id: 

1986 return pdef 

1987 

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

1989 

1990 def fix_ttt_holes(self, sptree, mode): 

1991 raise StoreError( 

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

1993 % self.short_type) 

1994 

1995 

1996class ConfigTypeA(Config): 

1997 ''' 

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

1999 

2000 * Problem is invariant to horizontal translations and rotations around 

2001 vertical axis. 

2002 

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

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

2005 component)`` 

2006 

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

2008 points. 

2009 ''' 

2010 

2011 receiver_depth = Float.T( 

2012 default=0.0, 

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

2014 

2015 source_depth_min = Float.T( 

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

2017 

2018 source_depth_max = Float.T( 

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

2020 

2021 source_depth_delta = Float.T( 

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

2023 

2024 distance_min = Float.T( 

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

2026 

2027 distance_max = Float.T( 

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

2029 

2030 distance_delta = Float.T( 

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

2032 

2033 short_type = 'A' 

2034 

2035 provided_schemes = [ 

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

2037 

2038 def get_surface_distance(self, args): 

2039 return args[1] 

2040 

2041 def get_distance(self, args): 

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

2043 

2044 def get_source_depth(self, args): 

2045 return args[0] 

2046 

2047 def get_source_depths(self): 

2048 return self.coords[0] 

2049 

2050 def get_receiver_depth(self, args): 

2051 return self.receiver_depth 

2052 

2053 def _update(self): 

2054 self.mins = num.array( 

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

2056 self.maxs = num.array( 

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

2058 self.deltas = num.array( 

2059 [self.source_depth_delta, self.distance_delta], 

2060 dtype=float) 

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

2062 vicinity_eps).astype(int) + 1 

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

2064 self.deltat = 1.0/self.sample_rate 

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

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

2067 (mi, ma, n) in 

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

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

2070 

2071 self.nsource_depths, self.ndistances = self.ns 

2072 

2073 def _make_index_functions(self): 

2074 

2075 amin, bmin = self.mins 

2076 da, db = self.deltas 

2077 na, nb = self.ns 

2078 

2079 ng = self.ncomponents 

2080 

2081 def index_function(a, b, ig): 

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

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

2084 try: 

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

2086 except ValueError: 

2087 raise OutOfBounds() 

2088 

2089 def indices_function(a, b, ig): 

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

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

2092 try: 

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

2094 except ValueError: 

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

2096 try: 

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

2098 except ValueError: 

2099 raise OutOfBounds() 

2100 

2101 def grid_interpolation_coefficients(a, b): 

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

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

2104 return ias, ibs 

2105 

2106 def vicinity_function(a, b, ig): 

2107 ias, ibs = grid_interpolation_coefficients(a, b) 

2108 

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

2110 raise OutOfBounds() 

2111 

2112 indis = [] 

2113 weights = [] 

2114 for ia, va in ias: 

2115 iia = ia*nb*ng 

2116 for ib, vb in ibs: 

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

2118 weights.append(va*vb) 

2119 

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

2121 

2122 def vicinities_function(a, b, ig): 

2123 

2124 xa = (a - amin) / da 

2125 xb = (b - bmin) / db 

2126 

2127 xa_fl = num.floor(xa) 

2128 xa_ce = num.ceil(xa) 

2129 xb_fl = num.floor(xb) 

2130 xb_ce = num.ceil(xb) 

2131 va_fl = 1.0 - (xa - xa_fl) 

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

2133 vb_fl = 1.0 - (xb - xb_fl) 

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

2135 

2136 ia_fl = xa_fl.astype(int) 

2137 ia_ce = xa_ce.astype(int) 

2138 ib_fl = xb_fl.astype(int) 

2139 ib_ce = xb_ce.astype(int) 

2140 

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

2142 raise OutOfBounds() 

2143 

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

2145 raise OutOfBounds() 

2146 

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

2148 raise OutOfBounds() 

2149 

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

2151 raise OutOfBounds() 

2152 

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

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

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

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

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

2158 

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

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

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

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

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

2164 

2165 return irecords, weights 

2166 

2167 self._index_function = index_function 

2168 self._indices_function = indices_function 

2169 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2170 self._vicinity_function = vicinity_function 

2171 self._vicinities_function = vicinities_function 

2172 

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

2174 nc = icomponents.size 

2175 dists = source.distances_to(receiver) 

2176 n = dists.size 

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

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

2179 icomponents) 

2180 

2181 def make_indexing_args1(self, source, receiver): 

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

2183 

2184 @property 

2185 def short_extent(self): 

2186 return '%g:%g:%g x %g:%g:%g' % ( 

2187 self.source_depth_min/km, 

2188 self.source_depth_max/km, 

2189 self.source_depth_delta/km, 

2190 self.distance_min/km, 

2191 self.distance_max/km, 

2192 self.distance_delta/km) 

2193 

2194 def fix_ttt_holes(self, sptree, mode): 

2195 from pyrocko import eikonal_ext, spit 

2196 

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

2198 

2199 delta = self.deltas[-1] 

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

2201 

2202 nsources, ndistances = self.ns 

2203 

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

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

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

2207 

2208 speeds = self.get_material_property( 

2209 0., 0., points, 

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

2211 interpolation='multilinear') 

2212 

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

2214 

2215 times = sptree.interpolate_many(nodes) 

2216 

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

2218 times = times.reshape(speeds.shape) 

2219 

2220 try: 

2221 eikonal_ext.eikonal_solver_fmm_cartesian( 

2222 speeds, times, delta) 

2223 except eikonal_ext.EikonalExtError as e: 

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

2225 logger.debug( 

2226 'Got a warning from eikonal solver ' 

2227 '- may be ok...') 

2228 else: 

2229 raise 

2230 

2231 def func(x): 

2232 ibs, ics = \ 

2233 self.grid_interpolation_coefficients(*x) 

2234 

2235 t = 0 

2236 for ib, vb in ibs: 

2237 for ic, vc in ics: 

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

2239 

2240 return t 

2241 

2242 return spit.SPTree( 

2243 f=func, 

2244 ftol=sptree.ftol, 

2245 xbounds=sptree.xbounds, 

2246 xtols=sptree.xtols) 

2247 

2248 

2249class ConfigTypeB(Config): 

2250 ''' 

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

2252 

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

2254 receiver depth 

2255 

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

2257 receiver_distance, component)`` 

2258 ''' 

2259 

2260 receiver_depth_min = Float.T( 

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

2262 

2263 receiver_depth_max = Float.T( 

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

2265 

2266 receiver_depth_delta = Float.T( 

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

2268 

2269 source_depth_min = Float.T( 

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

2271 

2272 source_depth_max = Float.T( 

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

2274 

2275 source_depth_delta = Float.T( 

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

2277 

2278 distance_min = Float.T( 

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

2280 

2281 distance_max = Float.T( 

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

2283 

2284 distance_delta = Float.T( 

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

2286 

2287 short_type = 'B' 

2288 

2289 provided_schemes = [ 

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

2291 

2292 def get_distance(self, args): 

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

2294 

2295 def get_surface_distance(self, args): 

2296 return args[2] 

2297 

2298 def get_source_depth(self, args): 

2299 return args[1] 

2300 

2301 def get_receiver_depth(self, args): 

2302 return args[0] 

2303 

2304 def get_source_depths(self): 

2305 return self.coords[1] 

2306 

2307 def _update(self): 

2308 self.mins = num.array([ 

2309 self.receiver_depth_min, 

2310 self.source_depth_min, 

2311 self.distance_min], 

2312 dtype=float) 

2313 

2314 self.maxs = num.array([ 

2315 self.receiver_depth_max, 

2316 self.source_depth_max, 

2317 self.distance_max], 

2318 dtype=float) 

2319 

2320 self.deltas = num.array([ 

2321 self.receiver_depth_delta, 

2322 self.source_depth_delta, 

2323 self.distance_delta], 

2324 dtype=float) 

2325 

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

2327 vicinity_eps).astype(int) + 1 

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

2329 self.deltat = 1.0/self.sample_rate 

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

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

2332 (mi, ma, n) in 

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

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

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

2336 

2337 def _make_index_functions(self): 

2338 

2339 amin, bmin, cmin = self.mins 

2340 da, db, dc = self.deltas 

2341 na, nb, nc = self.ns 

2342 ng = self.ncomponents 

2343 

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

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

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

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

2348 try: 

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

2350 (na, nb, nc, ng)) 

2351 except ValueError: 

2352 raise OutOfBounds() 

2353 

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

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

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

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

2358 try: 

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

2360 (na, nb, nc, ng)) 

2361 except ValueError: 

2362 raise OutOfBounds() 

2363 

2364 def grid_interpolation_coefficients(a, b, c): 

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

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

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

2368 return ias, ibs, ics 

2369 

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

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

2372 

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

2374 raise OutOfBounds() 

2375 

2376 indis = [] 

2377 weights = [] 

2378 for ia, va in ias: 

2379 iia = ia*nb*nc*ng 

2380 for ib, vb in ibs: 

2381 iib = ib*nc*ng 

2382 for ic, vc in ics: 

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

2384 weights.append(va*vb*vc) 

2385 

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

2387 

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

2389 

2390 xa = (a - amin) / da 

2391 xb = (b - bmin) / db 

2392 xc = (c - cmin) / dc 

2393 

2394 xa_fl = num.floor(xa) 

2395 xa_ce = num.ceil(xa) 

2396 xb_fl = num.floor(xb) 

2397 xb_ce = num.ceil(xb) 

2398 xc_fl = num.floor(xc) 

2399 xc_ce = num.ceil(xc) 

2400 va_fl = 1.0 - (xa - xa_fl) 

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

2402 vb_fl = 1.0 - (xb - xb_fl) 

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

2404 vc_fl = 1.0 - (xc - xc_fl) 

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

2406 

2407 ia_fl = xa_fl.astype(int) 

2408 ia_ce = xa_ce.astype(int) 

2409 ib_fl = xb_fl.astype(int) 

2410 ib_ce = xb_ce.astype(int) 

2411 ic_fl = xc_fl.astype(int) 

2412 ic_ce = xc_ce.astype(int) 

2413 

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

2415 raise OutOfBounds() 

2416 

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

2418 raise OutOfBounds() 

2419 

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

2421 raise OutOfBounds() 

2422 

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

2424 raise OutOfBounds() 

2425 

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

2427 raise OutOfBounds() 

2428 

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

2430 raise OutOfBounds() 

2431 

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

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

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

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

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

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

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

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

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

2441 

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

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

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

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

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

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

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

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

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

2451 

2452 return irecords, weights 

2453 

2454 self._index_function = index_function 

2455 self._indices_function = indices_function 

2456 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2457 self._vicinity_function = vicinity_function 

2458 self._vicinities_function = vicinities_function 

2459 

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

2461 nc = icomponents.size 

2462 dists = source.distances_to(receiver) 

2463 n = dists.size 

2464 receiver_depths = num.empty(nc) 

2465 receiver_depths.fill(receiver.depth) 

2466 return (receiver_depths, 

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

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

2469 icomponents) 

2470 

2471 def make_indexing_args1(self, source, receiver): 

2472 return (receiver.depth, 

2473 source.depth, 

2474 source.distance_to(receiver)) 

2475 

2476 @property 

2477 def short_extent(self): 

2478 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2479 self.receiver_depth_min/km, 

2480 self.receiver_depth_max/km, 

2481 self.receiver_depth_delta/km, 

2482 self.source_depth_min/km, 

2483 self.source_depth_max/km, 

2484 self.source_depth_delta/km, 

2485 self.distance_min/km, 

2486 self.distance_max/km, 

2487 self.distance_delta/km) 

2488 

2489 def fix_ttt_holes(self, sptree, mode): 

2490 from pyrocko import eikonal_ext, spit 

2491 

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

2493 

2494 delta = self.deltas[-1] 

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

2496 

2497 nreceivers, nsources, ndistances = self.ns 

2498 

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

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

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

2502 

2503 speeds = self.get_material_property( 

2504 0., 0., points, 

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

2506 interpolation='multilinear') 

2507 

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

2509 

2510 receiver_times = [] 

2511 for ireceiver in range(nreceivers): 

2512 nodes = num.hstack([ 

2513 num_full( 

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

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

2516 nodes_sr]) 

2517 

2518 times = sptree.interpolate_many(nodes) 

2519 

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

2521 

2522 times = times.reshape(speeds.shape) 

2523 

2524 try: 

2525 eikonal_ext.eikonal_solver_fmm_cartesian( 

2526 speeds, times, delta) 

2527 except eikonal_ext.EikonalExtError as e: 

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

2529 logger.debug( 

2530 'Got a warning from eikonal solver ' 

2531 '- may be ok...') 

2532 else: 

2533 raise 

2534 

2535 receiver_times.append(times) 

2536 

2537 def func(x): 

2538 ias, ibs, ics = \ 

2539 self.grid_interpolation_coefficients(*x) 

2540 

2541 t = 0 

2542 for ia, va in ias: 

2543 times = receiver_times[ia] 

2544 for ib, vb in ibs: 

2545 for ic, vc in ics: 

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

2547 

2548 return t 

2549 

2550 return spit.SPTree( 

2551 f=func, 

2552 ftol=sptree.ftol, 

2553 xbounds=sptree.xbounds, 

2554 xtols=sptree.xtols) 

2555 

2556 

2557class ConfigTypeC(Config): 

2558 ''' 

2559 No symmetrical constraints, one fixed receiver position. 

2560 

2561 * Cartesian 3D source volume around a reference point 

2562 

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

2564 source_east_shift, source_north_shift, component)`` 

2565 ''' 

2566 

2567 receiver = Receiver.T( 

2568 help='Receiver location') 

2569 

2570 source_origin = Location.T( 

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

2572 

2573 source_depth_min = Float.T( 

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

2575 

2576 source_depth_max = Float.T( 

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

2578 

2579 source_depth_delta = Float.T( 

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

2581 

2582 source_east_shift_min = Float.T( 

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

2584 

2585 source_east_shift_max = Float.T( 

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

2587 

2588 source_east_shift_delta = Float.T( 

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

2590 

2591 source_north_shift_min = Float.T( 

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

2593 

2594 source_north_shift_max = Float.T( 

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

2596 

2597 source_north_shift_delta = Float.T( 

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

2599 

2600 short_type = 'C' 

2601 

2602 provided_schemes = ['elastic18'] 

2603 

2604 def get_surface_distance(self, args): 

2605 _, source_east_shift, source_north_shift, _ = args 

2606 sorig = self.source_origin 

2607 sloc = Location( 

2608 lat=sorig.lat, 

2609 lon=sorig.lon, 

2610 north_shift=sorig.north_shift + source_north_shift, 

2611 east_shift=sorig.east_shift + source_east_shift) 

2612 

2613 return self.receiver.distance_to(sloc) 

2614 

2615 def get_distance(self, args): 

2616 sdepth, source_east_shift, source_north_shift, _ = args 

2617 sorig = self.source_origin 

2618 sloc = Location( 

2619 lat=sorig.lat, 

2620 lon=sorig.lon, 

2621 north_shift=sorig.north_shift + source_north_shift, 

2622 east_shift=sorig.east_shift + source_east_shift) 

2623 

2624 return self.receiver.distance_3d_to(sloc) 

2625 

2626 def get_source_depth(self, args): 

2627 return args[0] 

2628 

2629 def get_receiver_depth(self, args): 

2630 return self.receiver.depth 

2631 

2632 def get_source_depths(self): 

2633 return self.coords[0] 

2634 

2635 def _update(self): 

2636 self.mins = num.array([ 

2637 self.source_depth_min, 

2638 self.source_east_shift_min, 

2639 self.source_north_shift_min], 

2640 dtype=float) 

2641 

2642 self.maxs = num.array([ 

2643 self.source_depth_max, 

2644 self.source_east_shift_max, 

2645 self.source_north_shift_max], 

2646 dtype=float) 

2647 

2648 self.deltas = num.array([ 

2649 self.source_depth_delta, 

2650 self.source_east_shift_delta, 

2651 self.source_north_shift_delta], 

2652 dtype=float) 

2653 

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

2655 vicinity_eps).astype(int) + 1 

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

2657 self.deltat = 1.0/self.sample_rate 

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

2659 

2660 self.coords = tuple( 

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

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

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

2664 

2665 def _make_index_functions(self): 

2666 

2667 amin, bmin, cmin = self.mins 

2668 da, db, dc = self.deltas 

2669 na, nb, nc = self.ns 

2670 ng = self.ncomponents 

2671 

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

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

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

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

2676 try: 

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

2678 (na, nb, nc, ng)) 

2679 except ValueError: 

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

2681 

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

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

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

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

2686 

2687 try: 

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

2689 (na, nb, nc, ng)) 

2690 except ValueError: 

2691 raise OutOfBounds() 

2692 

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

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

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

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

2697 

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

2699 raise OutOfBounds() 

2700 

2701 indis = [] 

2702 weights = [] 

2703 for ia, va in ias: 

2704 iia = ia*nb*nc*ng 

2705 for ib, vb in ibs: 

2706 iib = ib*nc*ng 

2707 for ic, vc in ics: 

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

2709 weights.append(va*vb*vc) 

2710 

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

2712 

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

2714 

2715 xa = (a-amin) / da 

2716 xb = (b-bmin) / db 

2717 xc = (c-cmin) / dc 

2718 

2719 xa_fl = num.floor(xa) 

2720 xa_ce = num.ceil(xa) 

2721 xb_fl = num.floor(xb) 

2722 xb_ce = num.ceil(xb) 

2723 xc_fl = num.floor(xc) 

2724 xc_ce = num.ceil(xc) 

2725 va_fl = 1.0 - (xa - xa_fl) 

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

2727 vb_fl = 1.0 - (xb - xb_fl) 

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

2729 vc_fl = 1.0 - (xc - xc_fl) 

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

2731 

2732 ia_fl = xa_fl.astype(int) 

2733 ia_ce = xa_ce.astype(int) 

2734 ib_fl = xb_fl.astype(int) 

2735 ib_ce = xb_ce.astype(int) 

2736 ic_fl = xc_fl.astype(int) 

2737 ic_ce = xc_ce.astype(int) 

2738 

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

2740 raise OutOfBounds() 

2741 

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

2743 raise OutOfBounds() 

2744 

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

2746 raise OutOfBounds() 

2747 

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

2749 raise OutOfBounds() 

2750 

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

2752 raise OutOfBounds() 

2753 

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

2755 raise OutOfBounds() 

2756 

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

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

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

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

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

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

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

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

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

2766 

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

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

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

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

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

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

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

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

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

2776 

2777 return irecords, weights 

2778 

2779 self._index_function = index_function 

2780 self._indices_function = indices_function 

2781 self._vicinity_function = vicinity_function 

2782 self._vicinities_function = vicinities_function 

2783 

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

2785 nc = icomponents.size 

2786 

2787 dists = source.distances_to(self.source_origin) 

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

2789 

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

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

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

2793 

2794 n = dists.size 

2795 

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

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

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

2799 icomponents) 

2800 

2801 def make_indexing_args1(self, source, receiver): 

2802 dist = source.distance_to(self.source_origin) 

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

2804 

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

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

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

2808 

2809 return (source_depth, 

2810 source_east_shift, 

2811 source_north_shift) 

2812 

2813 @property 

2814 def short_extent(self): 

2815 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2816 self.source_depth_min/km, 

2817 self.source_depth_max/km, 

2818 self.source_depth_delta/km, 

2819 self.source_east_shift_min/km, 

2820 self.source_east_shift_max/km, 

2821 self.source_east_shift_delta/km, 

2822 self.source_north_shift_min/km, 

2823 self.source_north_shift_max/km, 

2824 self.source_north_shift_delta/km) 

2825 

2826 

2827class Weighting(Object): 

2828 factor = Float.T(default=1.0) 

2829 

2830 

2831class Taper(Object): 

2832 tmin = Timing.T() 

2833 tmax = Timing.T() 

2834 tfade = Float.T(default=0.0) 

2835 shape = StringChoice.T( 

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

2837 default='cos', 

2838 optional=True) 

2839 

2840 

2841class SimplePattern(SObject): 

2842 

2843 _pool = {} 

2844 

2845 def __init__(self, pattern): 

2846 self._pattern = pattern 

2847 SObject.__init__(self) 

2848 

2849 def __str__(self): 

2850 return self._pattern 

2851 

2852 @property 

2853 def regex(self): 

2854 pool = SimplePattern._pool 

2855 if self.pattern not in pool: 

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

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

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

2859 

2860 return pool[self.pattern] 

2861 

2862 def match(self, s): 

2863 return self.regex.match(s) 

2864 

2865 

2866class WaveformType(StringChoice): 

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

2868 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2869 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2870 

2871 

2872class ChannelSelection(Object): 

2873 pattern = SimplePattern.T(optional=True) 

2874 min_sample_rate = Float.T(optional=True) 

2875 max_sample_rate = Float.T(optional=True) 

2876 

2877 

2878class StationSelection(Object): 

2879 includes = SimplePattern.T() 

2880 excludes = SimplePattern.T() 

2881 distance_min = Float.T(optional=True) 

2882 distance_max = Float.T(optional=True) 

2883 azimuth_min = Float.T(optional=True) 

2884 azimuth_max = Float.T(optional=True) 

2885 

2886 

2887class WaveformSelection(Object): 

2888 channel_selection = ChannelSelection.T(optional=True) 

2889 station_selection = StationSelection.T(optional=True) 

2890 taper = Taper.T() 

2891 # filter = FrequencyResponse.T() 

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

2893 weighting = Weighting.T(optional=True) 

2894 sample_rate = Float.T(optional=True) 

2895 gf_store_id = StringID.T(optional=True) 

2896 

2897 

2898def indi12(x, n): 

2899 ''' 

2900 Get linear interpolation index and weight. 

2901 ''' 

2902 

2903 r = round(x) 

2904 if abs(r - x) < vicinity_eps: 

2905 i = int(r) 

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

2907 raise OutOfBounds() 

2908 

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

2910 else: 

2911 f = math.floor(x) 

2912 i = int(f) 

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

2914 raise OutOfBounds() 

2915 

2916 v = x-f 

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

2918 

2919 

2920def float_or_none(s): 

2921 units = { 

2922 'k': 1e3, 

2923 'M': 1e6, 

2924 } 

2925 

2926 factor = 1.0 

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

2928 factor = units[s[-1]] 

2929 s = s[:-1] 

2930 if not s: 

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

2932 

2933 if s: 

2934 return float(s) * factor 

2935 else: 

2936 return None 

2937 

2938 

2939class GridSpecError(Exception): 

2940 def __init__(self, s): 

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

2942 

2943 

2944def parse_grid_spec(spec): 

2945 try: 

2946 result = [] 

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

2948 t = dspec.split('@') 

2949 num = start = stop = step = None 

2950 if len(t) == 2: 

2951 num = int(t[1]) 

2952 if num <= 0: 

2953 raise GridSpecError(spec) 

2954 

2955 elif len(t) > 2: 

2956 raise GridSpecError(spec) 

2957 

2958 s = t[0] 

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

2960 if len(v) == 1: 

2961 start = stop = v[0] 

2962 if len(v) >= 2: 

2963 start, stop = v[0:2] 

2964 if len(v) == 3: 

2965 step = v[2] 

2966 

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

2968 raise GridSpecError(spec) 

2969 

2970 if step == 0.0: 

2971 raise GridSpecError(spec) 

2972 

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

2974 

2975 except ValueError: 

2976 raise GridSpecError(spec) 

2977 

2978 return result 

2979 

2980 

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

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

2983 if start is None: 

2984 start = ma if swap else mi 

2985 if stop is None: 

2986 stop = mi if swap else ma 

2987 if step is None: 

2988 step = -inc if ma < mi else inc 

2989 if num is None: 

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

2991 raise GridSpecError() 

2992 

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

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

2995 if abs(stop-stop2) > eps: 

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

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

2998 else: 

2999 stop = stop2 

3000 

3001 if start == stop: 

3002 num = 1 

3003 

3004 return start, stop, num 

3005 

3006 

3007def nditer_outer(x): 

3008 return num.nditer( 

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

3010 

3011 

3012def nodes(xs): 

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

3014 nnodes = num.prod(ns) 

3015 ndim = len(xs) 

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

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

3018 x = xs[idim] 

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

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

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

3022 

3023 return nodes 

3024 

3025 

3026def filledi(x, n): 

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

3028 a.fill(x) 

3029 return a 

3030 

3031 

3032config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3033 

3034discretized_source_classes = [ 

3035 DiscretizedExplosionSource, 

3036 DiscretizedSFSource, 

3037 DiscretizedMTSource, 

3038 DiscretizedPorePressureSource] 

3039 

3040 

3041__all__ = ''' 

3042Earthmodel1D 

3043StringID 

3044ScopeType 

3045WaveformType 

3046QuantityType 

3047NearfieldTermsType 

3048Reference 

3049Region 

3050CircularRegion 

3051RectangularRegion 

3052PhaseSelect 

3053InvalidTimingSpecification 

3054Timing 

3055TPDef 

3056OutOfBounds 

3057Location 

3058Receiver 

3059'''.split() + [ 

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

3061ComponentScheme 

3062component_scheme_to_description 

3063component_schemes 

3064Config 

3065GridSpecError 

3066Weighting 

3067Taper 

3068SimplePattern 

3069WaveformType 

3070ChannelSelection 

3071StationSelection 

3072WaveformSelection 

3073nditer_outer 

3074dump 

3075load 

3076discretized_source_classes 

3077config_type_classes 

3078UnavailableScheme 

3079InterpolationMethod 

3080SeismosizerTrace 

3081SeismosizerResult 

3082Result 

3083StaticResult 

3084TimingAttributeError 

3085'''.split()