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

1425 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

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

8''' 

9 

10import math 

11import re 

12import fnmatch 

13import logging 

14 

15import numpy as num 

16from scipy.interpolate import interp1d 

17 

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

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

20 List, ValidationError, Timestamp, Tuple, Dict) 

21from pyrocko.guts import dump, load # noqa 

22from pyrocko.guts_array import literal, Array 

23from pyrocko.model import Location, gnss 

24 

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

26from pyrocko.util import num_full 

27 

28from .error import StoreError 

29 

30try: 

31 new_str = unicode 

32except NameError: 

33 new_str = str 

34 

35guts_prefix = 'pf' 

36 

37d2r = math.pi / 180. 

38r2d = 1.0 / d2r 

39km = 1000. 

40vicinity_eps = 1e-5 

41 

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

43 

44 

45def fmt_choices(cls): 

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

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

48 

49 

50class InterpolationMethod(StringChoice): 

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

52 

53 

54class SeismosizerTrace(Object): 

55 

56 codes = Tuple.T( 

57 4, String.T(), 

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

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

60 

61 data = Array.T( 

62 shape=(None,), 

63 dtype=num.float32, 

64 serialize_as='base64', 

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

66 help='numpy array with data samples') 

67 

68 deltat = Float.T( 

69 default=1.0, 

70 help='sampling interval [s]') 

71 

72 tmin = Timestamp.T( 

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

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

75 

76 def pyrocko_trace(self): 

77 c = self.codes 

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

79 ydata=self.data, 

80 deltat=self.deltat, 

81 tmin=self.tmin) 

82 

83 @classmethod 

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

85 d = dict( 

86 codes=tr.nslc_id, 

87 tmin=tr.tmin, 

88 deltat=tr.deltat, 

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

90 d.update(kwargs) 

91 return cls(**d) 

92 

93 

94class SeismosizerResult(Object): 

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

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

97 

98 

99class Result(SeismosizerResult): 

100 trace = SeismosizerTrace.T(optional=True) 

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

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

103 

104 

105class StaticResult(SeismosizerResult): 

106 result = Dict.T( 

107 String.T(), 

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

109 

110 

111class GNSSCampaignResult(StaticResult): 

112 campaign = gnss.GNSSCampaign.T( 

113 optional=True) 

114 

115 

116class SatelliteResult(StaticResult): 

117 

118 theta = Array.T( 

119 optional=True, 

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

121 

122 phi = Array.T( 

123 optional=True, 

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

125 

126 

127class KiteSceneResult(SatelliteResult): 

128 

129 shape = Tuple.T() 

130 

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

132 try: 

133 from kite import Scene 

134 except ImportError: 

135 raise ImportError('Kite not installed') 

136 

137 def reshape(arr): 

138 return arr.reshape(self.shape) 

139 

140 displacement = self.result[component] 

141 

142 displacement, theta, phi = map( 

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

144 

145 sc = Scene( 

146 displacement=displacement, 

147 theta=theta, phi=phi, 

148 config=self.config) 

149 

150 return sc 

151 

152 

153class ComponentSchemeDescription(Object): 

154 name = String.T() 

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

156 ncomponents = Int.T() 

157 default_stored_quantity = String.T(optional=True) 

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

159 

160 

161component_scheme_descriptions = [ 

162 ComponentSchemeDescription( 

163 name='elastic2', 

164 source_terms=['m0'], 

165 ncomponents=2, 

166 default_stored_quantity='displacement', 

167 provided_components=[ 

168 'n', 'e', 'd']), 

169 ComponentSchemeDescription( 

170 name='elastic8', 

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

172 ncomponents=8, 

173 default_stored_quantity='displacement', 

174 provided_components=[ 

175 'n', 'e', 'd']), 

176 ComponentSchemeDescription( 

177 name='elastic10', 

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

179 ncomponents=10, 

180 default_stored_quantity='displacement', 

181 provided_components=[ 

182 'n', 'e', 'd']), 

183 ComponentSchemeDescription( 

184 name='elastic18', 

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

186 ncomponents=18, 

187 default_stored_quantity='displacement', 

188 provided_components=[ 

189 'n', 'e', 'd']), 

190 ComponentSchemeDescription( 

191 name='elastic5', 

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

193 ncomponents=5, 

194 default_stored_quantity='displacement', 

195 provided_components=[ 

196 'n', 'e', 'd']), 

197 ComponentSchemeDescription( 

198 name='poroelastic10', 

199 source_terms=['pore_pressure'], 

200 ncomponents=10, 

201 default_stored_quantity=None, 

202 provided_components=[ 

203 'displacement.n', 'displacement.e', 'displacement.d', 

204 'vertical_tilt.n', 'vertical_tilt.e', 

205 'pore_pressure', 

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

207 

208 

209# new names? 

210# 'mt_to_displacement_1d' 

211# 'mt_to_displacement_1d_ff_only' 

212# 'mt_to_gravity_1d' 

213# 'mt_to_stress_1d' 

214# 'explosion_to_displacement_1d' 

215# 'sf_to_displacement_1d' 

216# 'mt_to_displacement_3d' 

217# 'mt_to_gravity_3d' 

218# 'mt_to_stress_3d' 

219# 'pore_pressure_to_displacement_1d' 

220# 'pore_pressure_to_vertical_tilt_1d' 

221# 'pore_pressure_to_pore_pressure_1d' 

222# 'pore_pressure_to_darcy_velocity_1d' 

223 

224 

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

226component_scheme_to_description = dict( 

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

228 

229 

230class ComponentScheme(StringChoice): 

231 ''' 

232 Different Green's Function component schemes are available: 

233 

234 ================= ========================================================= 

235 Name Description 

236 ================= ========================================================= 

237 ``elastic10`` Elastodynamic for 

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

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

240 sources only 

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

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

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

244 MT sources only 

245 ``elastic2`` Elastodynamic for 

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

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

248 isotropic sources only 

249 ``elastic5`` Elastodynamic for 

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

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

252 sources only 

253 ``elastic18`` Elastodynamic for 

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

255 sources only 

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

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

258 ================= ========================================================= 

259 ''' 

260 

261 choices = component_schemes 

262 

263 

264class Earthmodel1D(Object): 

265 dummy_for = cake.LayeredModel 

266 dummy_for_description = 'pyrocko.cake.LayeredModel' 

267 

268 class __T(TBase): 

269 def regularize_extra(self, val): 

270 if isinstance(val, str): 

271 val = cake.LayeredModel.from_scanlines( 

272 cake.read_nd_model_str(val)) 

273 

274 return val 

275 

276 def to_save(self, val): 

277 return literal(cake.write_nd_model_str(val)) 

278 

279 

280class StringID(StringPattern): 

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

282 

283 

284class ScopeType(StringChoice): 

285 choices = [ 

286 'global', 

287 'regional', 

288 'local', 

289 ] 

290 

291 

292class WaveType(StringChoice): 

293 choices = [ 

294 'full waveform', 

295 'bodywave', 

296 'P wave', 

297 'S wave', 

298 'surface wave', 

299 ] 

300 

301 

302class NearfieldTermsType(StringChoice): 

303 choices = [ 

304 'complete', 

305 'incomplete', 

306 'missing', 

307 ] 

308 

309 

310class QuantityType(StringChoice): 

311 choices = [ 

312 'displacement', 

313 'rotation', 

314 'velocity', 

315 'acceleration', 

316 'pressure', 

317 'tilt', 

318 'pore_pressure', 

319 'darcy_velocity', 

320 'vertical_tilt'] 

321 

322 

323class Reference(Object): 

324 id = StringID.T() 

325 type = String.T() 

326 title = Unicode.T() 

327 journal = Unicode.T(optional=True) 

328 volume = Unicode.T(optional=True) 

329 number = Unicode.T(optional=True) 

330 pages = Unicode.T(optional=True) 

331 year = String.T() 

332 note = Unicode.T(optional=True) 

333 issn = String.T(optional=True) 

334 doi = String.T(optional=True) 

335 url = String.T(optional=True) 

336 eprint = String.T(optional=True) 

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

338 publisher = Unicode.T(optional=True) 

339 keywords = Unicode.T(optional=True) 

340 note = Unicode.T(optional=True) 

341 abstract = Unicode.T(optional=True) 

342 

343 @classmethod 

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

345 

346 from pybtex.database.input import bibtex 

347 

348 parser = bibtex.Parser() 

349 

350 if filename is not None: 

351 bib_data = parser.parse_file(filename) 

352 elif stream is not None: 

353 bib_data = parser.parse_stream(stream) 

354 

355 references = [] 

356 

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

358 d = {} 

359 avail = entry.fields.keys() 

360 for prop in cls.T.properties: 

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

362 continue 

363 

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

365 

366 if 'author' in entry.persons: 

367 d['authors'] = [] 

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

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

370 

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

372 references.append(c) 

373 

374 return references 

375 

376 

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

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

379_cake_pat = cake.PhaseDef.allowed_characters_pattern 

380_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic 

381 

382_ppat = r'(' \ 

383 r'cake:' + _cake_pat + \ 

384 r'|iaspei:' + _iaspei_pat + \ 

385 r'|vel_surface:' + _fpat + \ 

386 r'|vel:' + _fpat + \ 

387 r'|stored:' + _spat + \ 

388 r')' 

389 

390 

391timing_regex_old = re.compile( 

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

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

394 

395 

396timing_regex = re.compile( 

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

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

399 

400 

401class PhaseSelect(StringChoice): 

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

403 

404 

405class InvalidTimingSpecification(ValidationError): 

406 pass 

407 

408 

409class TimingAttributeError(Exception): 

410 pass 

411 

412 

413class Timing(SObject): 

414 ''' 

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

416 

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

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

419 arrival or group of such arrivals. 

420 

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

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

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

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

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

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

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

428 

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

430 

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

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

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

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

435 velocity [km/s] 

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

437 

438 **Examples:** 

439 

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

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

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

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

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

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

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

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

448 undefined for a given geometry 

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

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

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

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

453 arrivals A, B, and C 

454 ''' 

455 

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

457 

458 if s is not None: 

459 offset_is = None 

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

461 try: 

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

463 

464 if s.endswith('S'): 

465 offset_is = 'slowness' 

466 

467 phase_defs = [] 

468 select = '' 

469 

470 except ValueError: 

471 matched = False 

472 m = timing_regex.match(s) 

473 if m: 

474 if m.group(3): 

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

476 elif m.group(19): 

477 phase_defs = [m.group(19)] 

478 else: 

479 phase_defs = [] 

480 

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

482 

483 offset = 0.0 

484 soff = m.group(27) 

485 if soff: 

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

487 if soff.endswith('S'): 

488 offset_is = 'slowness' 

489 elif soff.endswith('%'): 

490 offset_is = 'percent' 

491 

492 matched = True 

493 

494 else: 

495 m = timing_regex_old.match(s) 

496 if m: 

497 if m.group(3): 

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

499 elif m.group(5): 

500 phase_defs = [m.group(5)] 

501 else: 

502 phase_defs = [] 

503 

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

505 

506 offset = 0.0 

507 if m.group(6): 

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

509 

510 phase_defs = [ 

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

512 

513 matched = True 

514 

515 if not matched: 

516 raise InvalidTimingSpecification(s) 

517 

518 kwargs = dict( 

519 phase_defs=phase_defs, 

520 select=select, 

521 offset=offset, 

522 offset_is=offset_is) 

523 

524 SObject.__init__(self, **kwargs) 

525 

526 def __str__(self): 

527 s = [] 

528 if self.phase_defs: 

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

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

531 sphases = '{%s}' % sphases 

532 

533 if self.select: 

534 sphases = self.select + sphases 

535 

536 s.append(sphases) 

537 

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

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

540 if self.offset_is == 'slowness': 

541 s.append('S') 

542 elif self.offset_is == 'percent': 

543 s.append('%') 

544 

545 return ''.join(s) 

546 

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

548 try: 

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

550 phase_offset = get_phase( 

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

552 offset = phase_offset(args) 

553 else: 

554 offset = self.offset 

555 

556 if self.phase_defs: 

557 extra_args = () 

558 if attributes: 

559 extra_args = (attributes,) 

560 

561 phases = [ 

562 get_phase(phase_def, *extra_args) 

563 for phase_def in self.phase_defs] 

564 

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

566 results = [ 

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

568 for result in results] 

569 

570 results = [ 

571 result for result in results 

572 if result[0] is not None] 

573 

574 if self.select == 'first': 

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

576 elif self.select == 'last': 

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

578 

579 if not results: 

580 if attributes is not None: 

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

582 else: 

583 return None 

584 

585 else: 

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

587 if self.offset_is == 'percent': 

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

589 else: 

590 t = t + offset 

591 

592 if attributes is not None: 

593 return (t, ) + values 

594 else: 

595 return t 

596 

597 else: 

598 if attributes is not None: 

599 raise TimingAttributeError( 

600 'Cannot get phase attributes from offset-only ' 

601 'definition.') 

602 return offset 

603 

604 except spit.OutOfBounds: 

605 raise OutOfBounds(args) 

606 

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

608 offset = Float.T(default=0.0) 

609 offset_is = String.T(optional=True) 

610 select = PhaseSelect.T( 

611 default='', 

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

613 tuple(PhaseSelect.choices))) 

614 

615 

616def mkdefs(s): 

617 defs = [] 

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

619 try: 

620 defs.append(float(x)) 

621 except ValueError: 

622 if x.startswith('!'): 

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

624 else: 

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

626 

627 return defs 

628 

629 

630class TPDef(Object): 

631 ''' 

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

633 ''' 

634 

635 id = StringID.T( 

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

637 definition = String.T( 

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

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

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

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

642 

643 @property 

644 def phases(self): 

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

646 if isinstance(x, cake.PhaseDef)] 

647 

648 @property 

649 def horizontal_velocities(self): 

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

651 

652 

653class OutOfBounds(Exception): 

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

655 Exception.__init__(self) 

656 self.values = values 

657 self.reason = reason 

658 self.context = None 

659 

660 def __str__(self): 

661 scontext = '' 

662 if self.context: 

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

664 

665 if self.reason: 

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

667 

668 if self.values: 

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

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

671 else: 

672 return 'out of bounds%s' % scontext 

673 

674 

675class MultiLocation(Object): 

676 ''' 

677 Unstructured grid of point coordinates. 

678 ''' 

679 

680 lats = Array.T( 

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

682 help='Latitudes of targets.') 

683 lons = Array.T( 

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

685 help='Longitude of targets.') 

686 north_shifts = Array.T( 

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

688 help='North shifts of targets.') 

689 east_shifts = Array.T( 

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

691 help='East shifts of targets.') 

692 elevation = Array.T( 

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

694 help='Elevations of targets.') 

695 

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

697 self._coords5 = None 

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

699 

700 @property 

701 def coords5(self): 

702 if self._coords5 is not None: 

703 return self._coords5 

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

705 self.elevation] 

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

707 if not sizes: 

708 raise AttributeError('Empty StaticTarget') 

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

710 raise AttributeError('Inconsistent coordinate sizes.') 

711 

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

713 for idx, p in enumerate(props): 

714 if p is not None: 

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

716 

717 return self._coords5 

718 

719 @property 

720 def ncoords(self): 

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

722 return 0 

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

724 

725 def get_latlon(self): 

726 ''' 

727 Get all coordinates as lat/lon. 

728 

729 :returns: Coordinates in Latitude, Longitude 

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

731 ''' 

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

733 for i in range(self.ncoords): 

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

735 return latlons 

736 

737 def get_corner_coords(self): 

738 ''' 

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

740 

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

742 :rtype: tuple 

743 ''' 

744 latlon = self.get_latlon() 

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

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

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

748 

749 

750class Receiver(Location): 

751 codes = Tuple.T( 

752 3, String.T(), 

753 optional=True, 

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

755 

756 def pyrocko_station(self): 

757 from pyrocko import model 

758 

759 lat, lon = self.effective_latlon 

760 

761 return model.Station( 

762 network=self.codes[0], 

763 station=self.codes[1], 

764 location=self.codes[2], 

765 lat=lat, 

766 lon=lon, 

767 depth=self.depth) 

768 

769 

770def g(x, d): 

771 if x is None: 

772 return d 

773 else: 

774 return x 

775 

776 

777class UnavailableScheme(Exception): 

778 ''' 

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

780 ''' 

781 pass 

782 

783 

784class InvalidNComponents(Exception): 

785 ''' 

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

787 ''' 

788 pass 

789 

790 

791class DiscretizedSource(Object): 

792 ''' 

793 Base class for discretized sources. 

794 

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

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

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

798 source distributions are represented by subclasses of the 

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

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

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

802 explosion/implosion type source distributions. The geometry calculations 

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

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

805 

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

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

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

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

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

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

812 ''' 

813 

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

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

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

817 lat = Float.T(optional=True) 

818 lon = Float.T(optional=True) 

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

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

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

822 dl = Float.T(optional=True) 

823 dw = Float.T(optional=True) 

824 nl = Float.T(optional=True) 

825 nw = Float.T(optional=True) 

826 

827 @classmethod 

828 def check_scheme(cls, scheme): 

829 ''' 

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

831 

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

833 supported by this discretized source class. 

834 ''' 

835 

836 if scheme not in cls.provided_schemes: 

837 raise UnavailableScheme( 

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

839 (cls.__name__, scheme)) 

840 

841 def __init__(self, **kwargs): 

842 Object.__init__(self, **kwargs) 

843 self._latlons = None 

844 

845 def __setattr__(self, name, value): 

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

847 'lats', 'lons'): 

848 self.__dict__['_latlons'] = None 

849 

850 Object.__setattr__(self, name, value) 

851 

852 def get_source_terms(self, scheme): 

853 raise NotImplementedError() 

854 

855 def make_weights(self, receiver, scheme): 

856 raise NotImplementedError() 

857 

858 @property 

859 def effective_latlons(self): 

860 ''' 

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

862 ''' 

863 

864 if self._latlons is None: 

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

866 if (self.north_shifts is not None and 

867 self.east_shifts is not None): 

868 self._latlons = orthodrome.ne_to_latlon( 

869 self.lats, self.lons, 

870 self.north_shifts, self.east_shifts) 

871 else: 

872 self._latlons = self.lats, self.lons 

873 else: 

874 lat = g(self.lat, 0.0) 

875 lon = g(self.lon, 0.0) 

876 self._latlons = orthodrome.ne_to_latlon( 

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

878 

879 return self._latlons 

880 

881 @property 

882 def effective_north_shifts(self): 

883 if self.north_shifts is not None: 

884 return self.north_shifts 

885 else: 

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

887 

888 @property 

889 def effective_east_shifts(self): 

890 if self.east_shifts is not None: 

891 return self.east_shifts 

892 else: 

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

894 

895 def same_origin(self, receiver): 

896 ''' 

897 Check if receiver has same reference point. 

898 ''' 

899 

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

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

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

903 

904 def azibazis_to(self, receiver): 

905 ''' 

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

907 points. 

908 ''' 

909 

910 if self.same_origin(receiver): 

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

912 receiver.north_shift - self.north_shifts) 

913 bazis = azis + 180. 

914 else: 

915 slats, slons = self.effective_latlons 

916 rlat, rlon = receiver.effective_latlon 

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

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

919 

920 return azis, bazis 

921 

922 def distances_to(self, receiver): 

923 ''' 

924 Compute distances to receiver for all contained points. 

925 ''' 

926 if self.same_origin(receiver): 

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

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

929 

930 else: 

931 slats, slons = self.effective_latlons 

932 rlat, rlon = receiver.effective_latlon 

933 return orthodrome.distance_accurate50m_numpy(slats, slons, 

934 rlat, rlon) 

935 

936 def element_coords(self, i): 

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

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

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

940 else: 

941 lat = self.lat 

942 lon = self.lon 

943 

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

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

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

947 

948 else: 

949 north_shift = east_shift = 0.0 

950 

951 return lat, lon, north_shift, east_shift 

952 

953 def coords5(self): 

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

955 

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

957 xs[:, 0] = self.lats 

958 xs[:, 1] = self.lons 

959 else: 

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

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

962 

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

964 xs[:, 2] = self.north_shifts 

965 xs[:, 3] = self.east_shifts 

966 

967 xs[:, 4] = self.depths 

968 

969 return xs 

970 

971 @property 

972 def nelements(self): 

973 return self.times.size 

974 

975 @classmethod 

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

977 ''' 

978 Combine several discretized source models. 

979 

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

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

982 factors and reference times of the parameterized (undiscretized) 

983 sources match or are accounted for. 

984 ''' 

985 

986 first = sources[0] 

987 

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

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

990 'sources of same type.') 

991 

992 latlons = [] 

993 for s in sources: 

994 latlons.append(s.effective_latlons) 

995 

996 lats, lons = num.hstack(latlons) 

997 

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

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

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

1001 same_ref = num.all( 

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

1003 else: 

1004 same_ref = False 

1005 

1006 cat = num.concatenate 

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

1008 print(times) 

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

1010 

1011 if same_ref: 

1012 lat = first.lat 

1013 lon = first.lon 

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

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

1016 lats = None 

1017 lons = None 

1018 else: 

1019 lat = None 

1020 lon = None 

1021 north_shifts = None 

1022 east_shifts = None 

1023 

1024 return cls( 

1025 times=times, lat=lat, lon=lon, lats=lats, lons=lons, 

1026 north_shifts=north_shifts, east_shifts=east_shifts, 

1027 depths=depths, **kwargs) 

1028 

1029 def centroid_position(self): 

1030 moments = self.moments() 

1031 norm = num.sum(moments) 

1032 if norm != 0.0: 

1033 w = moments / num.sum(moments) 

1034 else: 

1035 w = num.ones(moments.size) 

1036 

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

1038 lats, lons = self.effective_latlons 

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

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

1041 else: 

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

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

1044 

1045 cn = num.sum(n*w) 

1046 ce = num.sum(e*w) 

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

1048 

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

1050 lat = clat 

1051 lon = clon 

1052 north_shift = 0. 

1053 east_shift = 0. 

1054 else: 

1055 lat = g(self.lat, 0.0) 

1056 lon = g(self.lon, 0.0) 

1057 north_shift = cn 

1058 east_shift = ce 

1059 

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

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

1062 return tuple(float(x) for x in 

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

1064 

1065 

1066class DiscretizedExplosionSource(DiscretizedSource): 

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

1068 

1069 provided_schemes = ( 

1070 'elastic2', 

1071 'elastic8', 

1072 'elastic10', 

1073 ) 

1074 

1075 def get_source_terms(self, scheme): 

1076 self.check_scheme(scheme) 

1077 

1078 if scheme == 'elastic2': 

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

1080 

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

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

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

1084 return m6s 

1085 else: 

1086 assert False 

1087 

1088 def make_weights(self, receiver, scheme): 

1089 self.check_scheme(scheme) 

1090 

1091 azis, bazis = self.azibazis_to(receiver) 

1092 

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

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

1095 

1096 m0s = self.m0s 

1097 n = azis.size 

1098 

1099 cat = num.concatenate 

1100 rep = num.repeat 

1101 

1102 if scheme == 'elastic2': 

1103 w_n = cb*m0s 

1104 g_n = filledi(0, n) 

1105 w_e = sb*m0s 

1106 g_e = filledi(0, n) 

1107 w_d = m0s 

1108 g_d = filledi(1, n) 

1109 

1110 elif scheme == 'elastic8': 

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

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

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

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

1115 w_d = cat((m0s, m0s)) 

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

1117 

1118 elif scheme == 'elastic10': 

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

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

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

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

1123 w_d = cat((m0s, m0s, m0s)) 

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

1125 

1126 else: 

1127 assert False 

1128 

1129 return ( 

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

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

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

1133 ) 

1134 

1135 def split(self): 

1136 from pyrocko.gf.seismosizer import ExplosionSource 

1137 sources = [] 

1138 for i in range(self.nelements): 

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

1140 sources.append(ExplosionSource( 

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

1142 lat=lat, 

1143 lon=lon, 

1144 north_shift=north_shift, 

1145 east_shift=east_shift, 

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

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

1148 

1149 return sources 

1150 

1151 def moments(self): 

1152 return self.m0s 

1153 

1154 def centroid(self): 

1155 from pyrocko.gf.seismosizer import ExplosionSource 

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

1157 self.centroid_position() 

1158 

1159 return ExplosionSource( 

1160 time=time, 

1161 lat=lat, 

1162 lon=lon, 

1163 north_shift=north_shift, 

1164 east_shift=east_shift, 

1165 depth=depth, 

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

1167 

1168 @classmethod 

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

1170 ''' 

1171 Combine several discretized source models. 

1172 

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

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

1175 factors and reference times of the parameterized (undiscretized) 

1176 sources match or are accounted for. 

1177 ''' 

1178 

1179 if 'm0s' not in kwargs: 

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

1181 

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

1183 **kwargs) 

1184 

1185 

1186class DiscretizedSFSource(DiscretizedSource): 

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

1188 

1189 provided_schemes = ( 

1190 'elastic5', 

1191 ) 

1192 

1193 def get_source_terms(self, scheme): 

1194 self.check_scheme(scheme) 

1195 

1196 return self.forces 

1197 

1198 def make_weights(self, receiver, scheme): 

1199 self.check_scheme(scheme) 

1200 

1201 azis, bazis = self.azibazis_to(receiver) 

1202 

1203 sa = num.sin(azis*d2r) 

1204 ca = num.cos(azis*d2r) 

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

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

1207 

1208 forces = self.forces 

1209 fn = forces[:, 0] 

1210 fe = forces[:, 1] 

1211 fd = forces[:, 2] 

1212 

1213 f0 = fd 

1214 f1 = ca * fn + sa * fe 

1215 f2 = ca * fe - sa * fn 

1216 

1217 n = azis.size 

1218 

1219 if scheme == 'elastic5': 

1220 ioff = 0 

1221 

1222 cat = num.concatenate 

1223 rep = num.repeat 

1224 

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

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

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

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

1229 w_d = cat((f0, f1)) 

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

1231 

1232 return ( 

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

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

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

1236 ) 

1237 

1238 @classmethod 

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

1240 ''' 

1241 Combine several discretized source models. 

1242 

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

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

1245 factors and reference times of the parameterized (undiscretized) 

1246 sources match or are accounted for. 

1247 ''' 

1248 

1249 if 'forces' not in kwargs: 

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

1251 

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

1253 

1254 def moments(self): 

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

1256 

1257 def centroid(self): 

1258 from pyrocko.gf.seismosizer import SFSource 

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

1260 self.centroid_position() 

1261 

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

1263 return SFSource( 

1264 time=time, 

1265 lat=lat, 

1266 lon=lon, 

1267 north_shift=north_shift, 

1268 east_shift=east_shift, 

1269 depth=depth, 

1270 fn=fn, 

1271 fe=fe, 

1272 fd=fd) 

1273 

1274 

1275class DiscretizedMTSource(DiscretizedSource): 

1276 m6s = Array.T( 

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

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

1279 

1280 provided_schemes = ( 

1281 'elastic8', 

1282 'elastic10', 

1283 'elastic18', 

1284 ) 

1285 

1286 def get_source_terms(self, scheme): 

1287 self.check_scheme(scheme) 

1288 return self.m6s 

1289 

1290 def make_weights(self, receiver, scheme): 

1291 self.check_scheme(scheme) 

1292 

1293 m6s = self.m6s 

1294 n = m6s.shape[0] 

1295 rep = num.repeat 

1296 

1297 if scheme == 'elastic18': 

1298 w_n = m6s.flatten() 

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

1300 w_e = m6s.flatten() 

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

1302 w_d = m6s.flatten() 

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

1304 

1305 else: 

1306 azis, bazis = self.azibazis_to(receiver) 

1307 

1308 sa = num.sin(azis*d2r) 

1309 ca = num.cos(azis*d2r) 

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

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

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

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

1314 

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

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

1317 f2 = m6s[:, 2] 

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

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

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

1321 

1322 cat = num.concatenate 

1323 

1324 if scheme == 'elastic8': 

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

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

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

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

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

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

1331 

1332 elif scheme == 'elastic10': 

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

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

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

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

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

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

1339 

1340 else: 

1341 assert False 

1342 

1343 return ( 

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

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

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

1347 ) 

1348 

1349 def split(self): 

1350 from pyrocko.gf.seismosizer import MTSource 

1351 sources = [] 

1352 for i in range(self.nelements): 

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

1354 sources.append(MTSource( 

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

1356 lat=lat, 

1357 lon=lon, 

1358 north_shift=north_shift, 

1359 east_shift=east_shift, 

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

1361 m6=self.m6s[i])) 

1362 

1363 return sources 

1364 

1365 def moments(self): 

1366 moments = num.array( 

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

1368 for m6 in self.m6s]) 

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

1370 

1371 def get_moment_rate(self, deltat=None): 

1372 moments = self.moments() 

1373 times = self.times 

1374 times -= times.min() 

1375 

1376 t_max = times.max() 

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

1378 mom_times[mom_times > t_max] = t_max 

1379 

1380 # Right open histrogram bins 

1381 mom, _ = num.histogram( 

1382 times, 

1383 bins=mom_times, 

1384 weights=moments) 

1385 

1386 deltat = num.diff(mom_times) 

1387 mom_rate = mom / deltat 

1388 

1389 return mom_rate, mom_times[1:] 

1390 

1391 def centroid(self): 

1392 from pyrocko.gf.seismosizer import MTSource 

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

1394 self.centroid_position() 

1395 

1396 return MTSource( 

1397 time=time, 

1398 lat=lat, 

1399 lon=lon, 

1400 north_shift=north_shift, 

1401 east_shift=east_shift, 

1402 depth=depth, 

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

1404 

1405 @classmethod 

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

1407 ''' 

1408 Combine several discretized source models. 

1409 

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

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

1412 factors and reference times of the parameterized (undiscretized) 

1413 sources match or are accounted for. 

1414 ''' 

1415 

1416 if 'm6s' not in kwargs: 

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

1418 

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

1420 

1421 

1422class DiscretizedPorePressureSource(DiscretizedSource): 

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

1424 

1425 provided_schemes = ( 

1426 'poroelastic10', 

1427 ) 

1428 

1429 def get_source_terms(self, scheme): 

1430 self.check_scheme(scheme) 

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

1432 

1433 def make_weights(self, receiver, scheme): 

1434 self.check_scheme(scheme) 

1435 

1436 azis, bazis = self.azibazis_to(receiver) 

1437 

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

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

1440 

1441 pp = self.pp 

1442 n = bazis.size 

1443 

1444 w_un = cb*pp 

1445 g_un = filledi(1, n) 

1446 w_ue = sb*pp 

1447 g_ue = filledi(1, n) 

1448 w_ud = pp 

1449 g_ud = filledi(0, n) 

1450 

1451 w_tn = cb*pp 

1452 g_tn = filledi(6, n) 

1453 w_te = sb*pp 

1454 g_te = filledi(6, n) 

1455 

1456 w_pp = pp 

1457 g_pp = filledi(7, n) 

1458 

1459 w_dvn = cb*pp 

1460 g_dvn = filledi(9, n) 

1461 w_dve = sb*pp 

1462 g_dve = filledi(9, n) 

1463 w_dvd = pp 

1464 g_dvd = filledi(8, n) 

1465 

1466 return ( 

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

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

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

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

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

1472 ('pore_pressure', w_pp, g_pp), 

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

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

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

1476 ) 

1477 

1478 def moments(self): 

1479 return self.pp 

1480 

1481 def centroid(self): 

1482 from pyrocko.gf.seismosizer import PorePressurePointSource 

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

1484 self.centroid_position() 

1485 

1486 return PorePressurePointSource( 

1487 time=time, 

1488 lat=lat, 

1489 lon=lon, 

1490 north_shift=north_shift, 

1491 east_shift=east_shift, 

1492 depth=depth, 

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

1494 

1495 @classmethod 

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

1497 ''' 

1498 Combine several discretized source models. 

1499 

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

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

1502 factors and reference times of the parameterized (undiscretized) 

1503 sources match or are accounted for. 

1504 ''' 

1505 

1506 if 'pp' not in kwargs: 

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

1508 

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

1510 **kwargs) 

1511 

1512 

1513class Region(Object): 

1514 name = String.T(optional=True) 

1515 

1516 

1517class RectangularRegion(Region): 

1518 lat_min = Float.T() 

1519 lat_max = Float.T() 

1520 lon_min = Float.T() 

1521 lon_max = Float.T() 

1522 

1523 

1524class CircularRegion(Region): 

1525 lat = Float.T() 

1526 lon = Float.T() 

1527 radius = Float.T() 

1528 

1529 

1530class Config(Object): 

1531 ''' 

1532 Green's function store meta information. 

1533 

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

1535 configuration types are: 

1536 

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

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

1539 

1540 * Problem is invariant to horizontal translations and rotations around 

1541 vertical axis. 

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

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

1544 component)`` 

1545 

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

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

1548 

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

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

1551 receiver_depth, component)`` 

1552 

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

1554 constraints but fixed receiver positions 

1555 

1556 * Cartesian source volume around a reference point 

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

1558 source_east_shift, source_north_shift, component)`` 

1559 ''' 

1560 

1561 id = StringID.T( 

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

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

1564 'letter.') 

1565 

1566 derived_from_id = StringID.T( 

1567 optional=True, 

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

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

1570 

1571 version = String.T( 

1572 default='1.0', 

1573 optional=True, 

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

1575 

1576 modelling_code_id = StringID.T( 

1577 optional=True, 

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

1579 

1580 author = Unicode.T( 

1581 optional=True, 

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

1583 

1584 author_email = String.T( 

1585 optional=True, 

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

1587 

1588 created_time = Timestamp.T( 

1589 optional=True, 

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

1591 

1592 regions = List.T( 

1593 Region.T(), 

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

1595 

1596 scope_type = ScopeType.T( 

1597 optional=True, 

1598 help='Distance range scope of the store ' 

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

1600 

1601 waveform_type = WaveType.T( 

1602 optional=True, 

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

1604 

1605 nearfield_terms = NearfieldTermsType.T( 

1606 optional=True, 

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

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

1609 

1610 description = String.T( 

1611 optional=True, 

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

1613 

1614 references = List.T( 

1615 Reference.T(), 

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

1617 'related work.') 

1618 

1619 earthmodel_1d = Earthmodel1D.T( 

1620 optional=True, 

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

1622 

1623 earthmodel_receiver_1d = Earthmodel1D.T( 

1624 optional=True, 

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

1626 

1627 can_interpolate_source = Bool.T( 

1628 optional=True, 

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

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

1631 

1632 can_interpolate_receiver = Bool.T( 

1633 optional=True, 

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

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

1636 

1637 frequency_min = Float.T( 

1638 optional=True, 

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

1640 

1641 frequency_max = Float.T( 

1642 optional=True, 

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

1644 

1645 sample_rate = Float.T( 

1646 optional=True, 

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

1648 

1649 factor = Float.T( 

1650 default=1.0, 

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

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

1653 optional=True) 

1654 

1655 component_scheme = ComponentScheme.T( 

1656 default='elastic10', 

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

1658 

1659 stored_quantity = QuantityType.T( 

1660 optional=True, 

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

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

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

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

1665 

1666 tabulated_phases = List.T( 

1667 TPDef.T(), 

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

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

1670 

1671 ncomponents = Int.T( 

1672 optional=True, 

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

1674 

1675 uuid = String.T( 

1676 optional=True, 

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

1678 'GF store for practical purposes.') 

1679 

1680 reference = String.T( 

1681 optional=True, 

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

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

1684 

1685 def __init__(self, **kwargs): 

1686 self._do_auto_updates = False 

1687 Object.__init__(self, **kwargs) 

1688 self._index_function = None 

1689 self._indices_function = None 

1690 self._vicinity_function = None 

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

1692 self._do_auto_updates = True 

1693 self.update() 

1694 

1695 def check_ncomponents(self): 

1696 ncomponents = component_scheme_to_description[ 

1697 self.component_scheme].ncomponents 

1698 

1699 if self.ncomponents is None: 

1700 self.ncomponents = ncomponents 

1701 elif ncomponents != self.ncomponents: 

1702 raise InvalidNComponents( 

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

1704 self.ncomponents, self.component_scheme)) 

1705 

1706 def __setattr__(self, name, value): 

1707 Object.__setattr__(self, name, value) 

1708 try: 

1709 self.T.get_property(name) 

1710 if self._do_auto_updates: 

1711 self.update() 

1712 

1713 except ValueError: 

1714 pass 

1715 

1716 def update(self): 

1717 self.check_ncomponents() 

1718 self._update() 

1719 self._make_index_functions() 

1720 

1721 def irecord(self, *args): 

1722 return self._index_function(*args) 

1723 

1724 def irecords(self, *args): 

1725 return self._indices_function(*args) 

1726 

1727 def vicinity(self, *args): 

1728 return self._vicinity_function(*args) 

1729 

1730 def vicinities(self, *args): 

1731 return self._vicinities_function(*args) 

1732 

1733 def grid_interpolation_coefficients(self, *args): 

1734 return self._grid_interpolation_coefficients(*args) 

1735 

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

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

1738 

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

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

1741 

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

1743 i = 0 

1744 arrs = [] 

1745 ntotal = 1 

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

1747 if gdef and len(gdef) > i: 

1748 sssn = gdef[i] 

1749 else: 

1750 sssn = (None,)*4 

1751 

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

1753 ntotal *= len(arr) 

1754 

1755 arrs.append(arr) 

1756 i += 1 

1757 

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

1759 return nditer_outer(arrs[:level]) 

1760 

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

1762 nthreads=0): 

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

1764 

1765 out = [] 

1766 delays = source.times 

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

1768 receiver, 

1769 self.component_scheme): 

1770 

1771 weights *= self.factor 

1772 

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

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

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

1776 

1777 return out 

1778 

1779 def short_info(self): 

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

1781 

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

1783 interpolation=None): 

1784 ''' 

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

1786 

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

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

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

1790 ``(lat, lon)`` 

1791 :param interpolation: interpolation method. Choose from 

1792 ``('nearest_neighbor', 'multilinear')`` 

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

1794 point 

1795 

1796 The default implementation retrieves and interpolates the shear moduli 

1797 from the contained 1D velocity profile. 

1798 ''' 

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

1800 parameter='shear_moduli', 

1801 interpolation=interpolation) 

1802 

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

1804 interpolation=None): 

1805 ''' 

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

1807 

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

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

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

1811 ``(lat, lon)`` 

1812 :param interpolation: interpolation method. Choose from 

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

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

1815 point 

1816 

1817 The default implementation retrieves and interpolates the lambda moduli 

1818 from the contained 1D velocity profile. 

1819 ''' 

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

1821 parameter='lambda_moduli', 

1822 interpolation=interpolation) 

1823 

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

1825 interpolation=None): 

1826 ''' 

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

1828 

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

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

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

1832 ``(lat, lon)`` 

1833 :param interpolation: interpolation method. Choose from 

1834 ``('nearest_neighbor', 'multilinear')`` 

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

1836 point 

1837 

1838 The default implementation retrieves and interpolates the lambda moduli 

1839 from the contained 1D velocity profile. 

1840 ''' 

1841 lambda_moduli = self.get_material_property( 

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

1843 interpolation=interpolation) 

1844 shear_moduli = self.get_material_property( 

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

1846 interpolation=interpolation) 

1847 return lambda_moduli + (2 / 3) * shear_moduli 

1848 

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

1850 ''' 

1851 Get Vs at given points from contained velocity model. 

1852 

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

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

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

1856 ``(lat, lon)`` 

1857 :param interpolation: interpolation method. Choose from 

1858 ``('nearest_neighbor', 'multilinear')`` 

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

1860 point 

1861 

1862 The default implementation retrieves and interpolates Vs 

1863 from the contained 1D velocity profile. 

1864 ''' 

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

1866 parameter='vs', 

1867 interpolation=interpolation) 

1868 

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

1870 ''' 

1871 Get Vp at given points from contained velocity model. 

1872 

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

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

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

1876 ``(lat, lon)`` 

1877 :param interpolation: interpolation method. Choose from 

1878 ``('nearest_neighbor', 'multilinear')`` 

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

1880 point 

1881 

1882 The default implementation retrieves and interpolates Vp 

1883 from the contained 1D velocity profile. 

1884 ''' 

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

1886 parameter='vp', 

1887 interpolation=interpolation) 

1888 

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

1890 ''' 

1891 Get rho at given points from contained velocity model. 

1892 

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

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

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

1896 ``(lat, lon)`` 

1897 :param interpolation: interpolation method. Choose from 

1898 ``('nearest_neighbor', 'multilinear')`` 

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

1900 point 

1901 

1902 The default implementation retrieves and interpolates rho 

1903 from the contained 1D velocity profile. 

1904 ''' 

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

1906 parameter='rho', 

1907 interpolation=interpolation) 

1908 

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

1910 interpolation=None): 

1911 

1912 if interpolation is None: 

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

1914 'multilinear', 'nearest_neighbor') 

1915 

1916 earthmod = self.earthmodel_1d 

1917 store_depth_profile = self.get_source_depths() 

1918 z_profile = earthmod.profile('z') 

1919 

1920 if parameter == 'vs': 

1921 vs_profile = earthmod.profile('vs') 

1922 profile = num.interp( 

1923 store_depth_profile, z_profile, vs_profile) 

1924 

1925 elif parameter == 'vp': 

1926 vp_profile = earthmod.profile('vp') 

1927 profile = num.interp( 

1928 store_depth_profile, z_profile, vp_profile) 

1929 

1930 elif parameter == 'rho': 

1931 rho_profile = earthmod.profile('rho') 

1932 

1933 profile = num.interp( 

1934 store_depth_profile, z_profile, rho_profile) 

1935 

1936 elif parameter == 'shear_moduli': 

1937 vs_profile = earthmod.profile('vs') 

1938 rho_profile = earthmod.profile('rho') 

1939 

1940 store_vs_profile = num.interp( 

1941 store_depth_profile, z_profile, vs_profile) 

1942 store_rho_profile = num.interp( 

1943 store_depth_profile, z_profile, rho_profile) 

1944 

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

1946 

1947 elif parameter == 'lambda_moduli': 

1948 vs_profile = earthmod.profile('vs') 

1949 vp_profile = earthmod.profile('vp') 

1950 rho_profile = earthmod.profile('rho') 

1951 

1952 store_vs_profile = num.interp( 

1953 store_depth_profile, z_profile, vs_profile) 

1954 store_vp_profile = num.interp( 

1955 store_depth_profile, z_profile, vp_profile) 

1956 store_rho_profile = num.interp( 

1957 store_depth_profile, z_profile, rho_profile) 

1958 

1959 profile = store_rho_profile * ( 

1960 num.power(store_vp_profile, 2) - 

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

1962 else: 

1963 raise TypeError( 

1964 'parameter %s not available' % parameter) 

1965 

1966 if interpolation == 'multilinear': 

1967 kind = 'linear' 

1968 elif interpolation == 'nearest_neighbor': 

1969 kind = 'nearest' 

1970 else: 

1971 raise TypeError( 

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

1973 

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

1975 

1976 try: 

1977 return interpolator(points[:, 2]) 

1978 except ValueError: 

1979 raise OutOfBounds() 

1980 

1981 def is_static(self): 

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

1983 if self.modelling_code_id.startswith(code): 

1984 return True 

1985 return False 

1986 

1987 def is_dynamic(self): 

1988 return not self.is_static() 

1989 

1990 def get_source_depths(self): 

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

1992 

1993 def get_tabulated_phase(self, phase_id): 

1994 ''' 

1995 Get tabulated phase definition. 

1996 ''' 

1997 

1998 for pdef in self.tabulated_phases: 

1999 if pdef.id == phase_id: 

2000 return pdef 

2001 

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

2003 

2004 def fix_ttt_holes(self, sptree, mode): 

2005 raise StoreError( 

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

2007 % self.short_type) 

2008 

2009 

2010class ConfigTypeA(Config): 

2011 ''' 

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

2013 

2014 * Problem is invariant to horizontal translations and rotations around 

2015 vertical axis. 

2016 

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

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

2019 component)`` 

2020 

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

2022 points. 

2023 ''' 

2024 

2025 receiver_depth = Float.T( 

2026 default=0.0, 

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

2028 

2029 source_depth_min = Float.T( 

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

2031 

2032 source_depth_max = Float.T( 

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

2034 

2035 source_depth_delta = Float.T( 

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

2037 

2038 distance_min = Float.T( 

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

2040 

2041 distance_max = Float.T( 

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

2043 

2044 distance_delta = Float.T( 

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

2046 

2047 short_type = 'A' 

2048 

2049 provided_schemes = [ 

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

2051 

2052 def get_surface_distance(self, args): 

2053 return args[1] 

2054 

2055 def get_distance(self, args): 

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

2057 

2058 def get_source_depth(self, args): 

2059 return args[0] 

2060 

2061 def get_source_depths(self): 

2062 return self.coords[0] 

2063 

2064 def get_receiver_depth(self, args): 

2065 return self.receiver_depth 

2066 

2067 def _update(self): 

2068 self.mins = num.array( 

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

2070 self.maxs = num.array( 

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

2072 self.deltas = num.array( 

2073 [self.source_depth_delta, self.distance_delta], 

2074 dtype=float) 

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

2076 vicinity_eps).astype(int) + 1 

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

2078 self.deltat = 1.0/self.sample_rate 

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

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

2081 (mi, ma, n) in 

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

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

2084 

2085 self.nsource_depths, self.ndistances = self.ns 

2086 

2087 def _make_index_functions(self): 

2088 

2089 amin, bmin = self.mins 

2090 da, db = self.deltas 

2091 na, nb = self.ns 

2092 

2093 ng = self.ncomponents 

2094 

2095 def index_function(a, b, ig): 

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

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

2098 try: 

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

2100 except ValueError: 

2101 raise OutOfBounds() 

2102 

2103 def indices_function(a, b, ig): 

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

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

2106 try: 

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

2108 except ValueError: 

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

2110 try: 

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

2112 except ValueError: 

2113 raise OutOfBounds() 

2114 

2115 def grid_interpolation_coefficients(a, b): 

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

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

2118 return ias, ibs 

2119 

2120 def vicinity_function(a, b, ig): 

2121 ias, ibs = grid_interpolation_coefficients(a, b) 

2122 

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

2124 raise OutOfBounds() 

2125 

2126 indis = [] 

2127 weights = [] 

2128 for ia, va in ias: 

2129 iia = ia*nb*ng 

2130 for ib, vb in ibs: 

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

2132 weights.append(va*vb) 

2133 

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

2135 

2136 def vicinities_function(a, b, ig): 

2137 

2138 xa = (a - amin) / da 

2139 xb = (b - bmin) / db 

2140 

2141 xa_fl = num.floor(xa) 

2142 xa_ce = num.ceil(xa) 

2143 xb_fl = num.floor(xb) 

2144 xb_ce = num.ceil(xb) 

2145 va_fl = 1.0 - (xa - xa_fl) 

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

2147 vb_fl = 1.0 - (xb - xb_fl) 

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

2149 

2150 ia_fl = xa_fl.astype(int) 

2151 ia_ce = xa_ce.astype(int) 

2152 ib_fl = xb_fl.astype(int) 

2153 ib_ce = xb_ce.astype(int) 

2154 

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

2156 raise OutOfBounds() 

2157 

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

2159 raise OutOfBounds() 

2160 

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

2162 raise OutOfBounds() 

2163 

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

2165 raise OutOfBounds() 

2166 

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

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

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

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

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

2172 

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

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

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

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

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

2178 

2179 return irecords, weights 

2180 

2181 self._index_function = index_function 

2182 self._indices_function = indices_function 

2183 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2184 self._vicinity_function = vicinity_function 

2185 self._vicinities_function = vicinities_function 

2186 

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

2188 nc = icomponents.size 

2189 dists = source.distances_to(receiver) 

2190 n = dists.size 

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

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

2193 icomponents) 

2194 

2195 def make_indexing_args1(self, source, receiver): 

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

2197 

2198 @property 

2199 def short_extent(self): 

2200 return '%g:%g:%g x %g:%g:%g' % ( 

2201 self.source_depth_min/km, 

2202 self.source_depth_max/km, 

2203 self.source_depth_delta/km, 

2204 self.distance_min/km, 

2205 self.distance_max/km, 

2206 self.distance_delta/km) 

2207 

2208 def fix_ttt_holes(self, sptree, mode): 

2209 from pyrocko import eikonal_ext, spit 

2210 

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

2212 

2213 delta = self.deltas[-1] 

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

2215 

2216 nsources, ndistances = self.ns 

2217 

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

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

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

2221 

2222 speeds = self.get_material_property( 

2223 0., 0., points, 

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

2225 interpolation='multilinear') 

2226 

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

2228 

2229 times = sptree.interpolate_many(nodes) 

2230 

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

2232 times = times.reshape(speeds.shape) 

2233 

2234 try: 

2235 eikonal_ext.eikonal_solver_fmm_cartesian( 

2236 speeds, times, delta) 

2237 except eikonal_ext.EikonalExtError as e: 

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

2239 logger.debug( 

2240 'Got a warning from eikonal solver ' 

2241 '- may be ok...') 

2242 else: 

2243 raise 

2244 

2245 def func(x): 

2246 ibs, ics = \ 

2247 self.grid_interpolation_coefficients(*x) 

2248 

2249 t = 0 

2250 for ib, vb in ibs: 

2251 for ic, vc in ics: 

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

2253 

2254 return t 

2255 

2256 return spit.SPTree( 

2257 f=func, 

2258 ftol=sptree.ftol, 

2259 xbounds=sptree.xbounds, 

2260 xtols=sptree.xtols) 

2261 

2262 

2263class ConfigTypeB(Config): 

2264 ''' 

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

2266 

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

2268 receiver depth 

2269 

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

2271 receiver_distance, component)`` 

2272 ''' 

2273 

2274 receiver_depth_min = Float.T( 

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

2276 

2277 receiver_depth_max = Float.T( 

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

2279 

2280 receiver_depth_delta = Float.T( 

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

2282 

2283 source_depth_min = Float.T( 

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

2285 

2286 source_depth_max = Float.T( 

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

2288 

2289 source_depth_delta = Float.T( 

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

2291 

2292 distance_min = Float.T( 

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

2294 

2295 distance_max = Float.T( 

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

2297 

2298 distance_delta = Float.T( 

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

2300 

2301 short_type = 'B' 

2302 

2303 provided_schemes = [ 

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

2305 

2306 def get_distance(self, args): 

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

2308 

2309 def get_surface_distance(self, args): 

2310 return args[2] 

2311 

2312 def get_source_depth(self, args): 

2313 return args[1] 

2314 

2315 def get_receiver_depth(self, args): 

2316 return args[0] 

2317 

2318 def get_source_depths(self): 

2319 return self.coords[1] 

2320 

2321 def _update(self): 

2322 self.mins = num.array([ 

2323 self.receiver_depth_min, 

2324 self.source_depth_min, 

2325 self.distance_min], 

2326 dtype=float) 

2327 

2328 self.maxs = num.array([ 

2329 self.receiver_depth_max, 

2330 self.source_depth_max, 

2331 self.distance_max], 

2332 dtype=float) 

2333 

2334 self.deltas = num.array([ 

2335 self.receiver_depth_delta, 

2336 self.source_depth_delta, 

2337 self.distance_delta], 

2338 dtype=float) 

2339 

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

2341 vicinity_eps).astype(int) + 1 

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

2343 self.deltat = 1.0/self.sample_rate 

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

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

2346 (mi, ma, n) in 

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

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

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

2350 

2351 def _make_index_functions(self): 

2352 

2353 amin, bmin, cmin = self.mins 

2354 da, db, dc = self.deltas 

2355 na, nb, nc = self.ns 

2356 ng = self.ncomponents 

2357 

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

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

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

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

2362 try: 

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

2364 (na, nb, nc, ng)) 

2365 except ValueError: 

2366 raise OutOfBounds() 

2367 

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

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

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

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

2372 try: 

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

2374 (na, nb, nc, ng)) 

2375 except ValueError: 

2376 raise OutOfBounds() 

2377 

2378 def grid_interpolation_coefficients(a, b, c): 

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

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

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

2382 return ias, ibs, ics 

2383 

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

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

2386 

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

2388 raise OutOfBounds() 

2389 

2390 indis = [] 

2391 weights = [] 

2392 for ia, va in ias: 

2393 iia = ia*nb*nc*ng 

2394 for ib, vb in ibs: 

2395 iib = ib*nc*ng 

2396 for ic, vc in ics: 

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

2398 weights.append(va*vb*vc) 

2399 

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

2401 

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

2403 

2404 xa = (a - amin) / da 

2405 xb = (b - bmin) / db 

2406 xc = (c - cmin) / dc 

2407 

2408 xa_fl = num.floor(xa) 

2409 xa_ce = num.ceil(xa) 

2410 xb_fl = num.floor(xb) 

2411 xb_ce = num.ceil(xb) 

2412 xc_fl = num.floor(xc) 

2413 xc_ce = num.ceil(xc) 

2414 va_fl = 1.0 - (xa - xa_fl) 

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

2416 vb_fl = 1.0 - (xb - xb_fl) 

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

2418 vc_fl = 1.0 - (xc - xc_fl) 

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

2420 

2421 ia_fl = xa_fl.astype(int) 

2422 ia_ce = xa_ce.astype(int) 

2423 ib_fl = xb_fl.astype(int) 

2424 ib_ce = xb_ce.astype(int) 

2425 ic_fl = xc_fl.astype(int) 

2426 ic_ce = xc_ce.astype(int) 

2427 

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

2429 raise OutOfBounds() 

2430 

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

2432 raise OutOfBounds() 

2433 

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

2435 raise OutOfBounds() 

2436 

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

2438 raise OutOfBounds() 

2439 

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

2441 raise OutOfBounds() 

2442 

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

2444 raise OutOfBounds() 

2445 

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

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

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

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

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

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

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

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

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

2455 

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

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

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

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

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

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

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

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

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

2465 

2466 return irecords, weights 

2467 

2468 self._index_function = index_function 

2469 self._indices_function = indices_function 

2470 self._grid_interpolation_coefficients = grid_interpolation_coefficients 

2471 self._vicinity_function = vicinity_function 

2472 self._vicinities_function = vicinities_function 

2473 

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

2475 nc = icomponents.size 

2476 dists = source.distances_to(receiver) 

2477 n = dists.size 

2478 receiver_depths = num.empty(nc) 

2479 receiver_depths.fill(receiver.depth) 

2480 return (receiver_depths, 

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

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

2483 icomponents) 

2484 

2485 def make_indexing_args1(self, source, receiver): 

2486 return (receiver.depth, 

2487 source.depth, 

2488 source.distance_to(receiver)) 

2489 

2490 @property 

2491 def short_extent(self): 

2492 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2493 self.receiver_depth_min/km, 

2494 self.receiver_depth_max/km, 

2495 self.receiver_depth_delta/km, 

2496 self.source_depth_min/km, 

2497 self.source_depth_max/km, 

2498 self.source_depth_delta/km, 

2499 self.distance_min/km, 

2500 self.distance_max/km, 

2501 self.distance_delta/km) 

2502 

2503 def fix_ttt_holes(self, sptree, mode): 

2504 from pyrocko import eikonal_ext, spit 

2505 

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

2507 

2508 delta = self.deltas[-1] 

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

2510 

2511 nreceivers, nsources, ndistances = self.ns 

2512 

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

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

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

2516 

2517 speeds = self.get_material_property( 

2518 0., 0., points, 

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

2520 interpolation='multilinear') 

2521 

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

2523 

2524 receiver_times = [] 

2525 for ireceiver in range(nreceivers): 

2526 nodes = num.hstack([ 

2527 num_full( 

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

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

2530 nodes_sr]) 

2531 

2532 times = sptree.interpolate_many(nodes) 

2533 

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

2535 

2536 times = times.reshape(speeds.shape) 

2537 

2538 try: 

2539 eikonal_ext.eikonal_solver_fmm_cartesian( 

2540 speeds, times, delta) 

2541 except eikonal_ext.EikonalExtError as e: 

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

2543 logger.debug( 

2544 'Got a warning from eikonal solver ' 

2545 '- may be ok...') 

2546 else: 

2547 raise 

2548 

2549 receiver_times.append(times) 

2550 

2551 def func(x): 

2552 ias, ibs, ics = \ 

2553 self.grid_interpolation_coefficients(*x) 

2554 

2555 t = 0 

2556 for ia, va in ias: 

2557 times = receiver_times[ia] 

2558 for ib, vb in ibs: 

2559 for ic, vc in ics: 

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

2561 

2562 return t 

2563 

2564 return spit.SPTree( 

2565 f=func, 

2566 ftol=sptree.ftol, 

2567 xbounds=sptree.xbounds, 

2568 xtols=sptree.xtols) 

2569 

2570 

2571class ConfigTypeC(Config): 

2572 ''' 

2573 No symmetrical constraints, one fixed receiver position. 

2574 

2575 * Cartesian 3D source volume around a reference point 

2576 

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

2578 source_east_shift, source_north_shift, component)`` 

2579 ''' 

2580 

2581 receiver = Receiver.T( 

2582 help='Receiver location') 

2583 

2584 source_origin = Location.T( 

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

2586 

2587 source_depth_min = Float.T( 

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

2589 

2590 source_depth_max = Float.T( 

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

2592 

2593 source_depth_delta = Float.T( 

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

2595 

2596 source_east_shift_min = Float.T( 

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

2598 

2599 source_east_shift_max = Float.T( 

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

2601 

2602 source_east_shift_delta = Float.T( 

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

2604 

2605 source_north_shift_min = Float.T( 

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

2607 

2608 source_north_shift_max = Float.T( 

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

2610 

2611 source_north_shift_delta = Float.T( 

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

2613 

2614 short_type = 'C' 

2615 

2616 provided_schemes = ['elastic18'] 

2617 

2618 def get_surface_distance(self, args): 

2619 _, source_east_shift, source_north_shift, _ = args 

2620 sorig = self.source_origin 

2621 sloc = Location( 

2622 lat=sorig.lat, 

2623 lon=sorig.lon, 

2624 north_shift=sorig.north_shift + source_north_shift, 

2625 east_shift=sorig.east_shift + source_east_shift) 

2626 

2627 return self.receiver.distance_to(sloc) 

2628 

2629 def get_distance(self, args): 

2630 sdepth, source_east_shift, source_north_shift, _ = args 

2631 sorig = self.source_origin 

2632 sloc = Location( 

2633 lat=sorig.lat, 

2634 lon=sorig.lon, 

2635 north_shift=sorig.north_shift + source_north_shift, 

2636 east_shift=sorig.east_shift + source_east_shift) 

2637 

2638 return self.receiver.distance_3d_to(sloc) 

2639 

2640 def get_source_depth(self, args): 

2641 return args[0] 

2642 

2643 def get_receiver_depth(self, args): 

2644 return self.receiver.depth 

2645 

2646 def get_source_depths(self): 

2647 return self.coords[0] 

2648 

2649 def _update(self): 

2650 self.mins = num.array([ 

2651 self.source_depth_min, 

2652 self.source_east_shift_min, 

2653 self.source_north_shift_min], 

2654 dtype=float) 

2655 

2656 self.maxs = num.array([ 

2657 self.source_depth_max, 

2658 self.source_east_shift_max, 

2659 self.source_north_shift_max], 

2660 dtype=float) 

2661 

2662 self.deltas = num.array([ 

2663 self.source_depth_delta, 

2664 self.source_east_shift_delta, 

2665 self.source_north_shift_delta], 

2666 dtype=float) 

2667 

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

2669 vicinity_eps).astype(int) + 1 

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

2671 self.deltat = 1.0/self.sample_rate 

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

2673 

2674 self.coords = tuple( 

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

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

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

2678 

2679 def _make_index_functions(self): 

2680 

2681 amin, bmin, cmin = self.mins 

2682 da, db, dc = self.deltas 

2683 na, nb, nc = self.ns 

2684 ng = self.ncomponents 

2685 

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

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

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

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

2690 try: 

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

2692 (na, nb, nc, ng)) 

2693 except ValueError: 

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

2695 

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

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

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

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

2700 

2701 try: 

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

2703 (na, nb, nc, ng)) 

2704 except ValueError: 

2705 raise OutOfBounds() 

2706 

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

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

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

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

2711 

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

2713 raise OutOfBounds() 

2714 

2715 indis = [] 

2716 weights = [] 

2717 for ia, va in ias: 

2718 iia = ia*nb*nc*ng 

2719 for ib, vb in ibs: 

2720 iib = ib*nc*ng 

2721 for ic, vc in ics: 

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

2723 weights.append(va*vb*vc) 

2724 

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

2726 

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

2728 

2729 xa = (a-amin) / da 

2730 xb = (b-bmin) / db 

2731 xc = (c-cmin) / dc 

2732 

2733 xa_fl = num.floor(xa) 

2734 xa_ce = num.ceil(xa) 

2735 xb_fl = num.floor(xb) 

2736 xb_ce = num.ceil(xb) 

2737 xc_fl = num.floor(xc) 

2738 xc_ce = num.ceil(xc) 

2739 va_fl = 1.0 - (xa - xa_fl) 

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

2741 vb_fl = 1.0 - (xb - xb_fl) 

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

2743 vc_fl = 1.0 - (xc - xc_fl) 

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

2745 

2746 ia_fl = xa_fl.astype(int) 

2747 ia_ce = xa_ce.astype(int) 

2748 ib_fl = xb_fl.astype(int) 

2749 ib_ce = xb_ce.astype(int) 

2750 ic_fl = xc_fl.astype(int) 

2751 ic_ce = xc_ce.astype(int) 

2752 

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

2754 raise OutOfBounds() 

2755 

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

2757 raise OutOfBounds() 

2758 

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

2760 raise OutOfBounds() 

2761 

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

2763 raise OutOfBounds() 

2764 

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

2766 raise OutOfBounds() 

2767 

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

2769 raise OutOfBounds() 

2770 

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

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

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

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

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

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

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

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

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

2780 

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

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

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

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

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

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

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

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

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

2790 

2791 return irecords, weights 

2792 

2793 self._index_function = index_function 

2794 self._indices_function = indices_function 

2795 self._vicinity_function = vicinity_function 

2796 self._vicinities_function = vicinities_function 

2797 

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

2799 nc = icomponents.size 

2800 

2801 dists = source.distances_to(self.source_origin) 

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

2803 

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

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

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

2807 

2808 n = dists.size 

2809 

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

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

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

2813 icomponents) 

2814 

2815 def make_indexing_args1(self, source, receiver): 

2816 dist = source.distance_to(self.source_origin) 

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

2818 

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

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

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

2822 

2823 return (source_depth, 

2824 source_east_shift, 

2825 source_north_shift) 

2826 

2827 @property 

2828 def short_extent(self): 

2829 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % ( 

2830 self.source_depth_min/km, 

2831 self.source_depth_max/km, 

2832 self.source_depth_delta/km, 

2833 self.source_east_shift_min/km, 

2834 self.source_east_shift_max/km, 

2835 self.source_east_shift_delta/km, 

2836 self.source_north_shift_min/km, 

2837 self.source_north_shift_max/km, 

2838 self.source_north_shift_delta/km) 

2839 

2840 

2841class Weighting(Object): 

2842 factor = Float.T(default=1.0) 

2843 

2844 

2845class Taper(Object): 

2846 tmin = Timing.T() 

2847 tmax = Timing.T() 

2848 tfade = Float.T(default=0.0) 

2849 shape = StringChoice.T( 

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

2851 default='cos', 

2852 optional=True) 

2853 

2854 

2855class SimplePattern(SObject): 

2856 

2857 _pool = {} 

2858 

2859 def __init__(self, pattern): 

2860 self._pattern = pattern 

2861 SObject.__init__(self) 

2862 

2863 def __str__(self): 

2864 return self._pattern 

2865 

2866 @property 

2867 def regex(self): 

2868 pool = SimplePattern._pool 

2869 if self.pattern not in pool: 

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

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

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

2873 

2874 return pool[self.pattern] 

2875 

2876 def match(self, s): 

2877 return self.regex.match(s) 

2878 

2879 

2880class WaveformType(StringChoice): 

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

2882 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc', 

2883 'envelope_dis', 'envelope_vel', 'envelope_acc'] 

2884 

2885 

2886class ChannelSelection(Object): 

2887 pattern = SimplePattern.T(optional=True) 

2888 min_sample_rate = Float.T(optional=True) 

2889 max_sample_rate = Float.T(optional=True) 

2890 

2891 

2892class StationSelection(Object): 

2893 includes = SimplePattern.T() 

2894 excludes = SimplePattern.T() 

2895 distance_min = Float.T(optional=True) 

2896 distance_max = Float.T(optional=True) 

2897 azimuth_min = Float.T(optional=True) 

2898 azimuth_max = Float.T(optional=True) 

2899 

2900 

2901class WaveformSelection(Object): 

2902 channel_selection = ChannelSelection.T(optional=True) 

2903 station_selection = StationSelection.T(optional=True) 

2904 taper = Taper.T() 

2905 # filter = FrequencyResponse.T() 

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

2907 weighting = Weighting.T(optional=True) 

2908 sample_rate = Float.T(optional=True) 

2909 gf_store_id = StringID.T(optional=True) 

2910 

2911 

2912def indi12(x, n): 

2913 ''' 

2914 Get linear interpolation index and weight. 

2915 ''' 

2916 

2917 r = round(x) 

2918 if abs(r - x) < vicinity_eps: 

2919 i = int(r) 

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

2921 raise OutOfBounds() 

2922 

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

2924 else: 

2925 f = math.floor(x) 

2926 i = int(f) 

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

2928 raise OutOfBounds() 

2929 

2930 v = x-f 

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

2932 

2933 

2934def float_or_none(s): 

2935 units = { 

2936 'k': 1e3, 

2937 'M': 1e6, 

2938 } 

2939 

2940 factor = 1.0 

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

2942 factor = units[s[-1]] 

2943 s = s[:-1] 

2944 if not s: 

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

2946 

2947 if s: 

2948 return float(s) * factor 

2949 else: 

2950 return None 

2951 

2952 

2953class GridSpecError(Exception): 

2954 def __init__(self, s): 

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

2956 

2957 

2958def parse_grid_spec(spec): 

2959 try: 

2960 result = [] 

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

2962 t = dspec.split('@') 

2963 num = start = stop = step = None 

2964 if len(t) == 2: 

2965 num = int(t[1]) 

2966 if num <= 0: 

2967 raise GridSpecError(spec) 

2968 

2969 elif len(t) > 2: 

2970 raise GridSpecError(spec) 

2971 

2972 s = t[0] 

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

2974 if len(v) == 1: 

2975 start = stop = v[0] 

2976 if len(v) >= 2: 

2977 start, stop = v[0:2] 

2978 if len(v) == 3: 

2979 step = v[2] 

2980 

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

2982 raise GridSpecError(spec) 

2983 

2984 if step == 0.0: 

2985 raise GridSpecError(spec) 

2986 

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

2988 

2989 except ValueError: 

2990 raise GridSpecError(spec) 

2991 

2992 return result 

2993 

2994 

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

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

2997 if start is None: 

2998 start = ma if swap else mi 

2999 if stop is None: 

3000 stop = mi if swap else ma 

3001 if step is None: 

3002 step = -inc if ma < mi else inc 

3003 if num is None: 

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

3005 raise GridSpecError() 

3006 

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

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

3009 if abs(stop-stop2) > eps: 

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

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

3012 else: 

3013 stop = stop2 

3014 

3015 if start == stop: 

3016 num = 1 

3017 

3018 return start, stop, num 

3019 

3020 

3021def nditer_outer(x): 

3022 return num.nditer( 

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

3024 

3025 

3026def nodes(xs): 

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

3028 nnodes = num.prod(ns) 

3029 ndim = len(xs) 

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

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

3032 x = xs[idim] 

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

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

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

3036 

3037 return nodes 

3038 

3039 

3040def filledi(x, n): 

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

3042 a.fill(x) 

3043 return a 

3044 

3045 

3046config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC] 

3047 

3048discretized_source_classes = [ 

3049 DiscretizedExplosionSource, 

3050 DiscretizedSFSource, 

3051 DiscretizedMTSource, 

3052 DiscretizedPorePressureSource] 

3053 

3054 

3055__all__ = ''' 

3056Earthmodel1D 

3057StringID 

3058ScopeType 

3059WaveformType 

3060QuantityType 

3061NearfieldTermsType 

3062Reference 

3063Region 

3064CircularRegion 

3065RectangularRegion 

3066PhaseSelect 

3067InvalidTimingSpecification 

3068Timing 

3069TPDef 

3070OutOfBounds 

3071Location 

3072Receiver 

3073'''.split() + [ 

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

3075ComponentScheme 

3076component_scheme_to_description 

3077component_schemes 

3078Config 

3079GridSpecError 

3080Weighting 

3081Taper 

3082SimplePattern 

3083WaveformType 

3084ChannelSelection 

3085StationSelection 

3086WaveformSelection 

3087nditer_outer 

3088dump 

3089load 

3090discretized_source_classes 

3091config_type_classes 

3092UnavailableScheme 

3093InterpolationMethod 

3094SeismosizerTrace 

3095SeismosizerResult 

3096Result 

3097StaticResult 

3098TimingAttributeError 

3099'''.split()