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

2470 statements  

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

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7High level synthetic seismogram synthesis. 

8 

9.. _coordinate-system-names: 

10 

11Coordinate systems 

12.................. 

13 

14Coordinate system names commonly used in source models. 

15 

16================= ============================================ 

17Name Description 

18================= ============================================ 

19``'xyz'`` northing, easting, depth in [m] 

20``'xy'`` northing, easting in [m] 

21``'latlon'`` latitude, longitude in [deg] 

22``'lonlat'`` longitude, latitude in [deg] 

23``'latlondepth'`` latitude, longitude in [deg], depth in [m] 

24================= ============================================ 

25''' 

26 

27from collections import defaultdict 

28from functools import cmp_to_key 

29import time 

30import math 

31import os 

32import re 

33import logging 

34try: 

35 import resource 

36except ImportError: 

37 resource = None 

38from hashlib import sha1 

39 

40import numpy as num 

41from scipy.interpolate import RegularGridInterpolator 

42 

43from pyrocko.guts import (Object, Float, String, StringChoice, List, 

44 Timestamp, Int, SObject, ArgumentError, Dict, 

45 ValidationError, Bool) 

46from pyrocko.guts_array import Array 

47 

48from pyrocko import moment_tensor as pmt 

49from pyrocko import trace, util, config, model, eikonal_ext 

50from pyrocko.orthodrome import ne_to_latlon 

51from pyrocko.model import Location 

52from pyrocko.modelling import OkadaSource, make_okada_coefficient_matrix, \ 

53 okada_ext, invert_fault_dislocations_bem 

54 

55from . import meta, store, ws 

56from .tractions import TractionField, DirectedTractions 

57from .targets import Target, StaticTarget, SatelliteTarget 

58 

59pjoin = os.path.join 

60 

61guts_prefix = 'pf' 

62 

63d2r = math.pi / 180. 

64r2d = 180. / math.pi 

65km = 1e3 

66 

67logger = logging.getLogger('pyrocko.gf.seismosizer') 

68 

69 

70def cmp_none_aware(a, b): 

71 if isinstance(a, tuple) and isinstance(b, tuple): 

72 for xa, xb in zip(a, b): 

73 rv = cmp_none_aware(xa, xb) 

74 if rv != 0: 

75 return rv 

76 

77 return 0 

78 

79 anone = a is None 

80 bnone = b is None 

81 

82 if anone and bnone: 

83 return 0 

84 

85 if anone: 

86 return -1 

87 

88 if bnone: 

89 return 1 

90 

91 return bool(a > b) - bool(a < b) 

92 

93 

94def xtime(): 

95 return time.time() 

96 

97 

98class SeismosizerError(Exception): 

99 pass 

100 

101 

102class BadRequest(SeismosizerError): 

103 pass 

104 

105 

106class DuplicateStoreId(Exception): 

107 pass 

108 

109 

110class NoDefaultStoreSet(Exception): 

111 ''' 

112 Raised, when a default store would be used but none is set. 

113 ''' 

114 pass 

115 

116 

117class ConversionError(Exception): 

118 pass 

119 

120 

121class NoSuchStore(BadRequest): 

122 

123 def __init__(self, store_id=None, dirs=None): 

124 BadRequest.__init__(self) 

125 self.store_id = store_id 

126 self.dirs = dirs 

127 

128 def __str__(self): 

129 if self.store_id is not None: 

130 rstr = 'no GF store with id "%s" found.' % self.store_id 

131 else: 

132 rstr = 'GF store not found.' 

133 

134 if self.dirs is not None: 

135 rstr += ' Searched folders:\n %s' % '\n '.join(sorted(self.dirs)) 

136 return rstr 

137 

138 

139def ufloat(s): 

140 units = { 

141 'k': 1e3, 

142 'M': 1e6, 

143 } 

144 

145 factor = 1.0 

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

147 factor = units[s[-1]] 

148 s = s[:-1] 

149 if not s: 

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

151 

152 return float(s) * factor 

153 

154 

155def ufloat_or_none(s): 

156 if s: 

157 return ufloat(s) 

158 else: 

159 return None 

160 

161 

162def int_or_none(s): 

163 if s: 

164 return int(s) 

165 else: 

166 return None 

167 

168 

169def nonzero(x, eps=1e-15): 

170 return abs(x) > eps 

171 

172 

173def permudef(ln, j=0): 

174 if j < len(ln): 

175 k, v = ln[j] 

176 for y in v: 

177 ln[j] = k, y 

178 for s in permudef(ln, j + 1): 

179 yield s 

180 

181 ln[j] = k, v 

182 return 

183 else: 

184 yield ln 

185 

186 

187def arr(x): 

188 return num.atleast_1d(num.asarray(x)) 

189 

190 

191def discretize_rect_source(deltas, deltat, time, north, east, depth, 

192 strike, dip, length, width, 

193 anchor, velocity=None, stf=None, 

194 nucleation_x=None, nucleation_y=None, 

195 decimation_factor=1, pointsonly=False, 

196 plane_coords=False, 

197 aggressive_oversampling=False): 

198 

199 if stf is None: 

200 stf = STF() 

201 

202 if not velocity and not pointsonly: 

203 raise AttributeError('velocity is required in time mode') 

204 

205 mindeltagf = float(num.min(deltas)) 

206 if velocity: 

207 mindeltagf = min(mindeltagf, deltat * velocity) 

208 

209 ln = length 

210 wd = width 

211 

212 if aggressive_oversampling: 

213 nl = int((2. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1 

214 nw = int((2. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1 

215 else: 

216 nl = int((1. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1 

217 nw = int((1. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1 

218 

219 n = int(nl * nw) 

220 

221 dl = ln / nl 

222 dw = wd / nw 

223 

224 xl = num.linspace(-0.5 * (ln - dl), 0.5 * (ln - dl), nl) 

225 xw = num.linspace(-0.5 * (wd - dw), 0.5 * (wd - dw), nw) 

226 

227 points = num.zeros((n, 3)) 

228 points[:, 0] = num.tile(xl, nw) 

229 points[:, 1] = num.repeat(xw, nl) 

230 

231 if nucleation_x is not None: 

232 dist_x = num.abs(nucleation_x - points[:, 0]) 

233 else: 

234 dist_x = num.zeros(n) 

235 

236 if nucleation_y is not None: 

237 dist_y = num.abs(nucleation_y - points[:, 1]) 

238 else: 

239 dist_y = num.zeros(n) 

240 

241 dist = num.sqrt(dist_x**2 + dist_y**2) 

242 times = dist / velocity 

243 

244 anch_x, anch_y = map_anchor[anchor] 

245 

246 points[:, 0] -= anch_x * 0.5 * length 

247 points[:, 1] -= anch_y * 0.5 * width 

248 

249 if plane_coords: 

250 return points, dl, dw, nl, nw 

251 

252 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0) 

253 points = num.dot(rotmat.T, points.T).T 

254 

255 points[:, 0] += north 

256 points[:, 1] += east 

257 points[:, 2] += depth 

258 

259 if pointsonly: 

260 return points, dl, dw, nl, nw 

261 

262 xtau, amplitudes = stf.discretize_t(deltat, time) 

263 nt = xtau.size 

264 

265 points2 = num.repeat(points, nt, axis=0) 

266 times2 = (times[:, num.newaxis] + xtau[num.newaxis, :]).ravel() 

267 amplitudes2 = num.tile(amplitudes, n) 

268 

269 return points2, times2, amplitudes2, dl, dw, nl, nw 

270 

271 

272def check_rect_source_discretisation(points2, nl, nw, store): 

273 # We assume a non-rotated fault plane 

274 N_CRITICAL = 8 

275 points = points2.T.reshape((3, nl, nw)) 

276 if points.size <= N_CRITICAL: 

277 logger.warning('RectangularSource is defined by only %d sub-sources!' 

278 % points.size) 

279 return True 

280 

281 distances = num.sqrt( 

282 (points[0, 0, :] - points[0, 1, :])**2 + 

283 (points[1, 0, :] - points[1, 1, :])**2 + 

284 (points[2, 0, :] - points[2, 1, :])**2) 

285 

286 depths = points[2, 0, :] 

287 vs_profile = store.config.get_vs( 

288 lat=0., lon=0., 

289 points=num.repeat(depths[:, num.newaxis], 3, axis=1), 

290 interpolation='multilinear') 

291 

292 min_wavelength = vs_profile * (store.config.deltat * 2) 

293 if not num.all(min_wavelength > distances / 2): 

294 return False 

295 return True 

296 

297 

298def outline_rect_source(strike, dip, length, width, anchor): 

299 ln = length 

300 wd = width 

301 points = num.array( 

302 [[-0.5 * ln, -0.5 * wd, 0.], 

303 [0.5 * ln, -0.5 * wd, 0.], 

304 [0.5 * ln, 0.5 * wd, 0.], 

305 [-0.5 * ln, 0.5 * wd, 0.], 

306 [-0.5 * ln, -0.5 * wd, 0.]]) 

307 

308 anch_x, anch_y = map_anchor[anchor] 

309 points[:, 0] -= anch_x * 0.5 * length 

310 points[:, 1] -= anch_y * 0.5 * width 

311 

312 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0) 

313 

314 return num.dot(rotmat.T, points.T).T 

315 

316 

317def from_plane_coords( 

318 strike, dip, length, width, depth, x_plane_coords, y_plane_coords, 

319 lat=0., lon=0., 

320 north_shift=0, east_shift=0, 

321 anchor='top', cs='xy'): 

322 

323 ln = length 

324 wd = width 

325 x_abs = [] 

326 y_abs = [] 

327 if not isinstance(x_plane_coords, list): 

328 x_plane_coords = [x_plane_coords] 

329 y_plane_coords = [y_plane_coords] 

330 

331 for x_plane, y_plane in zip(x_plane_coords, y_plane_coords): 

332 points = num.array( 

333 [[-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.], 

334 [0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.], 

335 [0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.], 

336 [-0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.], 

337 [-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.]]) 

338 

339 anch_x, anch_y = map_anchor[anchor] 

340 points[:, 0] -= anch_x * 0.5 * length 

341 points[:, 1] -= anch_y * 0.5 * width 

342 

343 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0) 

344 

345 points = num.dot(rotmat.T, points.T).T 

346 points[:, 0] += north_shift 

347 points[:, 1] += east_shift 

348 points[:, 2] += depth 

349 if cs in ('latlon', 'lonlat'): 

350 latlon = ne_to_latlon(lat, lon, 

351 points[:, 0], points[:, 1]) 

352 latlon = num.array(latlon).T 

353 x_abs.append(latlon[1:2, 1]) 

354 y_abs.append(latlon[2:3, 0]) 

355 if cs == 'xy': 

356 x_abs.append(points[1:2, 1]) 

357 y_abs.append(points[2:3, 0]) 

358 

359 if cs == 'lonlat': 

360 return y_abs, x_abs 

361 else: 

362 return x_abs, y_abs 

363 

364 

365def points_on_rect_source( 

366 strike, dip, length, width, anchor, 

367 discretized_basesource=None, points_x=None, points_y=None): 

368 

369 ln = length 

370 wd = width 

371 

372 if isinstance(points_x, list) or isinstance(points_x, float): 

373 points_x = num.array([points_x]) 

374 if isinstance(points_y, list) or isinstance(points_y, float): 

375 points_y = num.array([points_y]) 

376 

377 if discretized_basesource: 

378 ds = discretized_basesource 

379 

380 nl_patches = ds.nl + 1 

381 nw_patches = ds.nw + 1 

382 

383 npoints = nl_patches * nw_patches 

384 points = num.zeros((npoints, 3)) 

385 ln_patches = num.array([il for il in range(nl_patches)]) 

386 wd_patches = num.array([iw for iw in range(nw_patches)]) 

387 

388 points_ln =\ 

389 2 * ((ln_patches - num.min(ln_patches)) / num.ptp(ln_patches)) - 1 

390 points_wd =\ 

391 2 * ((wd_patches - num.min(wd_patches)) / num.ptp(wd_patches)) - 1 

392 

393 for il in range(nl_patches): 

394 for iw in range(nw_patches): 

395 points[il * nw_patches + iw, :] = num.array([ 

396 points_ln[il] * ln * 0.5, 

397 points_wd[iw] * wd * 0.5, 0.0]) 

398 

399 elif points_x.shape[0] > 0 and points_y.shape[0] > 0: 

400 points = num.zeros(shape=((len(points_x), 3))) 

401 for i, (x, y) in enumerate(zip(points_x, points_y)): 

402 points[i, :] = num.array( 

403 [x * 0.5 * ln, y * 0.5 * wd, 0.0]) 

404 

405 anch_x, anch_y = map_anchor[anchor] 

406 

407 points[:, 0] -= anch_x * 0.5 * ln 

408 points[:, 1] -= anch_y * 0.5 * wd 

409 

410 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0) 

411 

412 return num.dot(rotmat.T, points.T).T 

413 

414 

415class InvalidGridDef(Exception): 

416 pass 

417 

418 

419class Range(SObject): 

420 ''' 

421 Convenient range specification. 

422 

423 Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]:: 

424 

425 Range('0 .. 10k : 1k') 

426 Range(start=0., stop=10e3, step=1e3) 

427 Range(0, 10e3, 1e3) 

428 Range('0 .. 10k @ 11') 

429 Range(start=0., stop=10*km, n=11) 

430 

431 Range(0, 10e3, n=11) 

432 Range(values=[x*1e3 for x in range(11)]) 

433 

434 Depending on the use context, it can be possible to omit any part of the 

435 specification. E.g. in the context of extracting a subset of an already 

436 existing range, the existing range's specification values would be filled 

437 in where missing. 

438 

439 The values are distributed with equal spacing, unless the ``spacing`` 

440 argument is modified. The values can be created offset or relative to an 

441 external base value with the ``relative`` argument if the use context 

442 supports this. 

443 

444 The range specification can be expressed with a short string 

445 representation:: 

446 

447 'start .. stop @ num | spacing, relative' 

448 'start .. stop : step | spacing, relative' 

449 

450 most parts of the expression can be omitted if not needed. Whitespace is 

451 allowed for readability but can also be omitted. 

452 ''' 

453 

454 start = Float.T(optional=True) 

455 stop = Float.T(optional=True) 

456 step = Float.T(optional=True) 

457 n = Int.T(optional=True) 

458 values = Array.T(optional=True, dtype=float, shape=(None,)) 

459 

460 spacing = StringChoice.T( 

461 choices=['lin', 'log', 'symlog'], 

462 default='lin', 

463 optional=True) 

464 

465 relative = StringChoice.T( 

466 choices=['', 'add', 'mult'], 

467 default='', 

468 optional=True) 

469 

470 pattern = re.compile(r'^((?P<start>.*)\.\.(?P<stop>[^@|:]*))?' 

471 r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?' 

472 r'(\|(?P<stuff>.+))?$') 

473 

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

475 d = {} 

476 if len(args) == 1: 

477 d = self.parse(args[0]) 

478 elif len(args) in (2, 3): 

479 d['start'], d['stop'] = [float(x) for x in args[:2]] 

480 if len(args) == 3: 

481 d['step'] = float(args[2]) 

482 

483 for k, v in kwargs.items(): 

484 if k in d: 

485 raise ArgumentError('%s specified more than once' % k) 

486 

487 d[k] = v 

488 

489 SObject.__init__(self, **d) 

490 

491 def __str__(self): 

492 def sfloat(x): 

493 if x is not None: 

494 return '%g' % x 

495 else: 

496 return '' 

497 

498 if self.values: 

499 return ','.join('%g' % x for x in self.values) 

500 

501 if self.start is None and self.stop is None: 

502 s0 = '' 

503 else: 

504 s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop)) 

505 

506 s1 = '' 

507 if self.step is not None: 

508 s1 = [' : %g', ':%g'][s0 == ''] % self.step 

509 elif self.n is not None: 

510 s1 = [' @ %i', '@%i'][s0 == ''] % self.n 

511 

512 if self.spacing == 'lin' and self.relative == '': 

513 s2 = '' 

514 else: 

515 x = [] 

516 if self.spacing != 'lin': 

517 x.append(self.spacing) 

518 

519 if self.relative != '': 

520 x.append(self.relative) 

521 

522 s2 = ' | %s' % ','.join(x) 

523 

524 return s0 + s1 + s2 

525 

526 @classmethod 

527 def parse(cls, s): 

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

529 m = cls.pattern.match(s) 

530 if not m: 

531 try: 

532 vals = [ufloat(x) for x in s.split(',')] 

533 except Exception: 

534 raise InvalidGridDef( 

535 '"%s" is not a valid range specification' % s) 

536 

537 return dict(values=num.array(vals, dtype=float)) 

538 

539 d = m.groupdict() 

540 try: 

541 start = ufloat_or_none(d['start']) 

542 stop = ufloat_or_none(d['stop']) 

543 step = ufloat_or_none(d['step']) 

544 n = int_or_none(d['n']) 

545 except Exception: 

546 raise InvalidGridDef( 

547 '"%s" is not a valid range specification' % s) 

548 

549 spacing = 'lin' 

550 relative = '' 

551 

552 if d['stuff'] is not None: 

553 t = d['stuff'].split(',') 

554 for x in t: 

555 if x in cls.spacing.choices: 

556 spacing = x 

557 elif x and x in cls.relative.choices: 

558 relative = x 

559 else: 

560 raise InvalidGridDef( 

561 '"%s" is not a valid range specification' % s) 

562 

563 return dict(start=start, stop=stop, step=step, n=n, spacing=spacing, 

564 relative=relative) 

565 

566 def make(self, mi=None, ma=None, inc=None, base=None, eps=1e-5): 

567 if self.values: 

568 return self.values 

569 

570 start = self.start 

571 stop = self.stop 

572 step = self.step 

573 n = self.n 

574 

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

576 if start is None: 

577 start = [mi, ma][swap] 

578 if stop is None: 

579 stop = [ma, mi][swap] 

580 if step is None and inc is not None: 

581 step = [inc, -inc][ma < mi] 

582 

583 if start is None or stop is None: 

584 raise InvalidGridDef( 

585 'Cannot use range specification "%s" without start ' 

586 'and stop in this context' % self) 

587 

588 if step is None and n is None: 

589 step = stop - start 

590 

591 if n is None: 

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

593 raise InvalidGridDef( 

594 'Range specification "%s" has inconsistent ordering ' 

595 '(step < 0 => stop > start)' % self) 

596 

597 n = int(round((stop - start) / step)) + 1 

598 stop2 = start + (n - 1) * step 

599 if abs(stop - stop2) > eps: 

600 n = int(math.floor((stop - start) / step)) + 1 

601 stop = start + (n - 1) * step 

602 else: 

603 stop = stop2 

604 

605 if start == stop: 

606 n = 1 

607 

608 if self.spacing == 'lin': 

609 vals = num.linspace(start, stop, n) 

610 

611 elif self.spacing in ('log', 'symlog'): 

612 if start > 0. and stop > 0.: 

613 vals = num.exp(num.linspace(num.log(start), 

614 num.log(stop), n)) 

615 elif start < 0. and stop < 0.: 

616 vals = -num.exp(num.linspace(num.log(-start), 

617 num.log(-stop), n)) 

618 else: 

619 raise InvalidGridDef( 

620 'Log ranges should not include or cross zero ' 

621 '(in range specification "%s").' % self) 

622 

623 if self.spacing == 'symlog': 

624 nvals = - vals 

625 vals = num.concatenate((nvals[::-1], vals)) 

626 

627 if self.relative in ('add', 'mult') and base is None: 

628 raise InvalidGridDef( 

629 'Cannot use relative range specification in this context.') 

630 

631 vals = self.make_relative(base, vals) 

632 

633 return list(map(float, vals)) 

634 

635 def make_relative(self, base, vals): 

636 if self.relative == 'add': 

637 vals += base 

638 

639 if self.relative == 'mult': 

640 vals *= base 

641 

642 return vals 

643 

644 

645class GridDefElement(Object): 

646 

647 param = meta.StringID.T() 

648 rs = Range.T() 

649 

650 def __init__(self, shorthand=None, **kwargs): 

651 if shorthand is not None: 

652 t = shorthand.split('=') 

653 if len(t) != 2: 

654 raise InvalidGridDef( 

655 'Invalid grid specification element: %s' % shorthand) 

656 

657 sp, sr = t[0].strip(), t[1].strip() 

658 

659 kwargs['param'] = sp 

660 kwargs['rs'] = Range(sr) 

661 

662 Object.__init__(self, **kwargs) 

663 

664 def shorthand(self): 

665 return self.param + ' = ' + str(self.rs) 

666 

667 

668class GridDef(Object): 

669 

670 elements = List.T(GridDefElement.T()) 

671 

672 def __init__(self, shorthand=None, **kwargs): 

673 if shorthand is not None: 

674 t = shorthand.splitlines() 

675 tt = [] 

676 for x in t: 

677 x = x.strip() 

678 if x: 

679 tt.extend(x.split(';')) 

680 

681 elements = [] 

682 for se in tt: 

683 elements.append(GridDef(se)) 

684 

685 kwargs['elements'] = elements 

686 

687 Object.__init__(self, **kwargs) 

688 

689 def shorthand(self): 

690 return '; '.join(str(x) for x in self.elements) 

691 

692 

693class Cloneable(object): 

694 ''' 

695 Mix-in class for Guts objects, providing dict-like access and cloning. 

696 ''' 

697 

698 def __iter__(self): 

699 return iter(self.T.propnames) 

700 

701 def __getitem__(self, k): 

702 if k not in self.keys(): 

703 raise KeyError(k) 

704 

705 return getattr(self, k) 

706 

707 def __setitem__(self, k, v): 

708 if k not in self.keys(): 

709 raise KeyError(k) 

710 

711 return setattr(self, k, v) 

712 

713 def clone(self, **kwargs): 

714 ''' 

715 Make a copy of the object. 

716 

717 A new object of the same class is created and initialized with the 

718 parameters of the object on which this method is called on. If 

719 ``kwargs`` are given, these are used to override any of the 

720 initialization parameters. 

721 ''' 

722 

723 d = dict(self) 

724 for k in d: 

725 v = d[k] 

726 if isinstance(v, Cloneable): 

727 d[k] = v.clone() 

728 

729 d.update(kwargs) 

730 return self.__class__(**d) 

731 

732 @classmethod 

733 def keys(cls): 

734 ''' 

735 Get list of the source model's parameter names. 

736 ''' 

737 

738 return cls.T.propnames 

739 

740 

741class STF(Object, Cloneable): 

742 

743 ''' 

744 Base class for source time functions. 

745 ''' 

746 

747 def __init__(self, effective_duration=None, **kwargs): 

748 if effective_duration is not None: 

749 kwargs['duration'] = effective_duration / \ 

750 self.factor_duration_to_effective() 

751 

752 Object.__init__(self, **kwargs) 

753 

754 @classmethod 

755 def factor_duration_to_effective(cls): 

756 return 1.0 

757 

758 def centroid_time(self, tref): 

759 return tref 

760 

761 @property 

762 def effective_duration(self): 

763 return self.duration * self.factor_duration_to_effective() 

764 

765 def discretize_t(self, deltat, tref): 

766 tl = math.floor(tref / deltat) * deltat 

767 th = math.ceil(tref / deltat) * deltat 

768 if tl == th: 

769 return num.array([tl], dtype=float), num.ones(1) 

770 else: 

771 return ( 

772 num.array([tl, th], dtype=float), 

773 num.array([th - tref, tref - tl], dtype=float) / deltat) 

774 

775 def base_key(self): 

776 return (type(self).__name__,) 

777 

778 

779g_unit_pulse = STF() 

780 

781 

782def sshift(times, amplitudes, tshift, deltat): 

783 

784 t0 = math.floor(tshift / deltat) * deltat 

785 t1 = math.ceil(tshift / deltat) * deltat 

786 if t0 == t1: 

787 return times, amplitudes 

788 

789 amplitudes2 = num.zeros(amplitudes.size + 1, dtype=float) 

790 

791 amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes 

792 amplitudes2[1:] += (tshift - t0) / deltat * amplitudes 

793 

794 times2 = num.arange(times.size + 1, dtype=float) * \ 

795 deltat + times[0] + t0 

796 

797 return times2, amplitudes2 

798 

799 

800class BoxcarSTF(STF): 

801 

802 ''' 

803 Boxcar type source time function. 

804 

805 .. figure :: /static/stf-BoxcarSTF.svg 

806 :width: 40% 

807 :align: center 

808 :alt: boxcar source time function 

809 ''' 

810 

811 duration = Float.T( 

812 default=0.0, 

813 help='duration of the boxcar') 

814 

815 anchor = Float.T( 

816 default=0.0, 

817 help='anchor point with respect to source.time: (' 

818 '-1.0: left -> source duration [0, T] ~ hypocenter time, ' 

819 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' 

820 '+1.0: right -> source duration [-T, 0] ~ rupture end time)') 

821 

822 @classmethod 

823 def factor_duration_to_effective(cls): 

824 return 1.0 

825 

826 def centroid_time(self, tref): 

827 return tref - 0.5 * self.duration * self.anchor 

828 

829 def discretize_t(self, deltat, tref): 

830 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 

831 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 

832 tmin = round(tmin_stf / deltat) * deltat 

833 tmax = round(tmax_stf / deltat) * deltat 

834 nt = int(round((tmax - tmin) / deltat)) + 1 

835 times = num.linspace(tmin, tmax, nt) 

836 amplitudes = num.ones_like(times) 

837 if times.size > 1: 

838 t_edges = num.linspace( 

839 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1) 

840 t = tmin_stf + self.duration * num.array( 

841 [0.0, 0.0, 1.0, 1.0], dtype=float) 

842 f = num.array([0., 1., 1., 0.], dtype=float) 

843 amplitudes = util.plf_integrate_piecewise(t_edges, t, f) 

844 amplitudes /= num.sum(amplitudes) 

845 

846 tshift = (num.sum(amplitudes * times) - self.centroid_time(tref)) 

847 

848 return sshift(times, amplitudes, -tshift, deltat) 

849 

850 def base_key(self): 

851 return (type(self).__name__, self.duration, self.anchor) 

852 

853 

854class TriangularSTF(STF): 

855 

856 ''' 

857 Triangular type source time function. 

858 

859 .. figure :: /static/stf-TriangularSTF.svg 

860 :width: 40% 

861 :align: center 

862 :alt: triangular source time function 

863 ''' 

864 

865 duration = Float.T( 

866 default=0.0, 

867 help='baseline of the triangle') 

868 

869 peak_ratio = Float.T( 

870 default=0.5, 

871 help='fraction of time compared to duration, ' 

872 'when the maximum amplitude is reached') 

873 

874 anchor = Float.T( 

875 default=0.0, 

876 help='anchor point with respect to source.time: (' 

877 '-1.0: left -> source duration [0, T] ~ hypocenter time, ' 

878 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' 

879 '+1.0: right -> source duration [-T, 0] ~ rupture end time)') 

880 

881 @classmethod 

882 def factor_duration_to_effective(cls, peak_ratio=None): 

883 if peak_ratio is None: 

884 peak_ratio = cls.peak_ratio.default() 

885 

886 return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0) 

887 

888 def __init__(self, effective_duration=None, **kwargs): 

889 if effective_duration is not None: 

890 kwargs['duration'] = effective_duration / \ 

891 self.factor_duration_to_effective( 

892 kwargs.get('peak_ratio', None)) 

893 

894 STF.__init__(self, **kwargs) 

895 

896 @property 

897 def centroid_ratio(self): 

898 ra = self.peak_ratio 

899 rb = 1.0 - ra 

900 return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb) 

901 

902 def centroid_time(self, tref): 

903 ca = self.centroid_ratio 

904 cb = 1.0 - ca 

905 if self.anchor <= 0.: 

906 return tref - ca * self.duration * self.anchor 

907 else: 

908 return tref - cb * self.duration * self.anchor 

909 

910 @property 

911 def effective_duration(self): 

912 return self.duration * self.factor_duration_to_effective( 

913 self.peak_ratio) 

914 

915 def tminmax_stf(self, tref): 

916 ca = self.centroid_ratio 

917 cb = 1.0 - ca 

918 if self.anchor <= 0.: 

919 tmin_stf = tref - ca * self.duration * (self.anchor + 1.) 

920 tmax_stf = tmin_stf + self.duration 

921 else: 

922 tmax_stf = tref + cb * self.duration * (1. - self.anchor) 

923 tmin_stf = tmax_stf - self.duration 

924 

925 return tmin_stf, tmax_stf 

926 

927 def discretize_t(self, deltat, tref): 

928 tmin_stf, tmax_stf = self.tminmax_stf(tref) 

929 

930 tmin = round(tmin_stf / deltat) * deltat 

931 tmax = round(tmax_stf / deltat) * deltat 

932 nt = int(round((tmax - tmin) / deltat)) + 1 

933 if nt > 1: 

934 t_edges = num.linspace( 

935 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1) 

936 t = tmin_stf + self.duration * num.array( 

937 [0.0, self.peak_ratio, 1.0], dtype=float) 

938 f = num.array([0., 1., 0.], dtype=float) 

939 amplitudes = util.plf_integrate_piecewise(t_edges, t, f) 

940 amplitudes /= num.sum(amplitudes) 

941 else: 

942 amplitudes = num.ones(1) 

943 

944 times = num.linspace(tmin, tmax, nt) 

945 return times, amplitudes 

946 

947 def base_key(self): 

948 return ( 

949 type(self).__name__, self.duration, self.peak_ratio, self.anchor) 

950 

951 

952class HalfSinusoidSTF(STF): 

953 

954 ''' 

955 Half sinusoid type source time function. 

956 

957 .. figure :: /static/stf-HalfSinusoidSTF.svg 

958 :width: 40% 

959 :align: center 

960 :alt: half-sinusouid source time function 

961 ''' 

962 

963 duration = Float.T( 

964 default=0.0, 

965 help='duration of the half-sinusoid (baseline)') 

966 

967 anchor = Float.T( 

968 default=0.0, 

969 help='anchor point with respect to source.time: (' 

970 '-1.0: left -> source duration [0, T] ~ hypocenter time, ' 

971 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' 

972 '+1.0: right -> source duration [-T, 0] ~ rupture end time)') 

973 

974 exponent = Int.T( 

975 default=1, 

976 help='set to 2 to use square of the half-period sinusoidal function.') 

977 

978 def __init__(self, effective_duration=None, **kwargs): 

979 if effective_duration is not None: 

980 kwargs['duration'] = effective_duration / \ 

981 self.factor_duration_to_effective( 

982 kwargs.get('exponent', 1)) 

983 

984 STF.__init__(self, **kwargs) 

985 

986 @classmethod 

987 def factor_duration_to_effective(cls, exponent): 

988 if exponent == 1: 

989 return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi 

990 elif exponent == 2: 

991 return math.sqrt(math.pi**2 - 6) / math.pi 

992 else: 

993 raise ValueError('Exponent for HalfSinusoidSTF must be 1 or 2.') 

994 

995 @property 

996 def effective_duration(self): 

997 return self.duration * self.factor_duration_to_effective(self.exponent) 

998 

999 def centroid_time(self, tref): 

1000 return tref - 0.5 * self.duration * self.anchor 

1001 

1002 def discretize_t(self, deltat, tref): 

1003 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 

1004 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 

1005 tmin = round(tmin_stf / deltat) * deltat 

1006 tmax = round(tmax_stf / deltat) * deltat 

1007 nt = int(round((tmax - tmin) / deltat)) + 1 

1008 if nt > 1: 

1009 t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace( 

1010 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1))) 

1011 

1012 if self.exponent == 1: 

1013 fint = -num.cos( 

1014 (t_edges - tmin_stf) * (math.pi / self.duration)) 

1015 

1016 elif self.exponent == 2: 

1017 fint = (t_edges - tmin_stf) / self.duration \ 

1018 - 1.0 / (2.0 * math.pi) * num.sin( 

1019 (t_edges - tmin_stf) * (2.0 * math.pi / self.duration)) 

1020 else: 

1021 raise ValueError( 

1022 'Exponent for HalfSinusoidSTF must be 1 or 2.') 

1023 

1024 amplitudes = fint[1:] - fint[:-1] 

1025 amplitudes /= num.sum(amplitudes) 

1026 else: 

1027 amplitudes = num.ones(1) 

1028 

1029 times = num.linspace(tmin, tmax, nt) 

1030 return times, amplitudes 

1031 

1032 def base_key(self): 

1033 return (type(self).__name__, self.duration, self.anchor) 

1034 

1035 

1036class SmoothRampSTF(STF): 

1037 ''' 

1038 Smooth-ramp type source time function for near-field displacement. 

1039 Based on moment function of double-couple point source proposed by [1]_. 

1040 

1041 .. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow 

1042 earthquakes from Love-wave modelling for regional distances, PEPI 32, 

1043 312-324. 

1044 

1045 .. figure :: /static/stf-SmoothRampSTF.svg 

1046 :width: 40% 

1047 :alt: smooth ramp source time function 

1048 ''' 

1049 duration = Float.T( 

1050 default=0.0, 

1051 help='duration of the ramp (baseline)') 

1052 

1053 rise_ratio = Float.T( 

1054 default=0.5, 

1055 help='fraction of time compared to duration, ' 

1056 'when the maximum amplitude is reached') 

1057 

1058 anchor = Float.T( 

1059 default=0.0, 

1060 help='anchor point with respect to source.time: (' 

1061 '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, ' 

1062 '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, ' 

1063 '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)') 

1064 

1065 def discretize_t(self, deltat, tref): 

1066 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 

1067 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 

1068 tmin = round(tmin_stf / deltat) * deltat 

1069 tmax = round(tmax_stf / deltat) * deltat 

1070 D = round((tmax - tmin) / deltat) * deltat 

1071 nt = int(round(D / deltat)) + 1 

1072 times = num.linspace(tmin, tmax, nt) 

1073 if nt > 1: 

1074 rise_time = self.rise_ratio * self.duration 

1075 amplitudes = num.ones_like(times) 

1076 tp = tmin + rise_time 

1077 ii = num.where(times <= tp) 

1078 t_inc = times[ii] 

1079 a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time) 

1080 b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0 

1081 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b) 

1082 

1083 amplitudes /= num.sum(amplitudes) 

1084 else: 

1085 amplitudes = num.ones(1) 

1086 

1087 return times, amplitudes 

1088 

1089 def base_key(self): 

1090 return (type(self).__name__, 

1091 self.duration, self.rise_ratio, self.anchor) 

1092 

1093 

1094class ResonatorSTF(STF): 

1095 ''' 

1096 Simple resonator like source time function. 

1097 

1098 .. math :: 

1099 

1100 f(t) = 0 for t < 0 

1101 f(t) = e^{-t/tau} * sin(2 * pi * f * t) 

1102 

1103 

1104 .. figure :: /static/stf-SmoothRampSTF.svg 

1105 :width: 40% 

1106 :alt: smooth ramp source time function 

1107 

1108 ''' 

1109 

1110 duration = Float.T( 

1111 default=0.0, 

1112 help='decay time') 

1113 

1114 frequency = Float.T( 

1115 default=1.0, 

1116 help='resonance frequency') 

1117 

1118 def discretize_t(self, deltat, tref): 

1119 tmin_stf = tref 

1120 tmax_stf = tref + self.duration * 3 

1121 tmin = math.floor(tmin_stf / deltat) * deltat 

1122 tmax = math.ceil(tmax_stf / deltat) * deltat 

1123 times = util.arange2(tmin, tmax, deltat) 

1124 amplitudes = num.exp(-(times - tref) / self.duration) \ 

1125 * num.sin(2.0 * num.pi * self.frequency * (times - tref)) 

1126 

1127 return times, amplitudes 

1128 

1129 def base_key(self): 

1130 return (type(self).__name__, 

1131 self.duration, self.frequency) 

1132 

1133 

1134class TremorSTF(STF): 

1135 ''' 

1136 Oszillating source time function. 

1137 

1138 .. math :: 

1139 

1140 f(t) = 0 for t < -tau/2 or t > tau/2 

1141 f(t) = cos(pi/tau*t) * sin(2 * pi * f * t) 

1142 

1143 ''' 

1144 

1145 duration = Float.T( 

1146 default=0.0, 

1147 help='Tremor duration [s]') 

1148 

1149 frequency = Float.T( 

1150 default=1.0, 

1151 help='Frequency [Hz]') 

1152 

1153 def discretize_t(self, deltat, tref): 

1154 tmin_stf = tref - 0.5 * self.duration 

1155 tmax_stf = tref + 0.5 * self.duration 

1156 tmin = math.floor(tmin_stf / deltat) * deltat 

1157 tmax = math.ceil(tmax_stf / deltat) * deltat 

1158 times = util.arange2(tmin, tmax, deltat) 

1159 mask = num.logical_and( 

1160 tref - 0.5 * self.duration < times, 

1161 times < tref + 0.5 * self.duration) 

1162 

1163 amplitudes = num.zeros_like(times) 

1164 amplitudes[mask] = num.cos(num.pi/self.duration*(times[mask] - tref)) \ 

1165 * num.sin(2.0 * num.pi * self.frequency * (times[mask] - tref)) 

1166 

1167 return times, amplitudes 

1168 

1169 def base_key(self): 

1170 return (type(self).__name__, 

1171 self.duration, self.frequency) 

1172 

1173 

1174class STFMode(StringChoice): 

1175 choices = ['pre', 'post'] 

1176 

1177 

1178class Source(Location, Cloneable): 

1179 ''' 

1180 Base class for all source models. 

1181 ''' 

1182 

1183 name = String.T(optional=True, default='') 

1184 

1185 time = Timestamp.T( 

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

1187 help='source origin time.') 

1188 

1189 stf = STF.T( 

1190 optional=True, 

1191 help='source time function.') 

1192 

1193 stf_mode = STFMode.T( 

1194 default='post', 

1195 help='whether to apply source time function in pre or ' 

1196 'post-processing.') 

1197 

1198 def __init__(self, **kwargs): 

1199 Location.__init__(self, **kwargs) 

1200 

1201 def update(self, **kwargs): 

1202 ''' 

1203 Change some of the source models parameters. 

1204 

1205 Example:: 

1206 

1207 >>> from pyrocko import gf 

1208 >>> s = gf.DCSource() 

1209 >>> s.update(strike=66., dip=33.) 

1210 >>> print(s) 

1211 --- !pf.DCSource 

1212 depth: 0.0 

1213 time: 1970-01-01 00:00:00 

1214 magnitude: 6.0 

1215 strike: 66.0 

1216 dip: 33.0 

1217 rake: 0.0 

1218 

1219 ''' 

1220 

1221 for (k, v) in kwargs.items(): 

1222 self[k] = v 

1223 

1224 def grid(self, **variables): 

1225 ''' 

1226 Create grid of source model variations. 

1227 

1228 :returns: :py:class:`SourceGrid` instance. 

1229 

1230 Example:: 

1231 

1232 >>> from pyrocko import gf 

1233 >>> base = DCSource() 

1234 >>> R = gf.Range 

1235 >>> for s in base.grid(R(' 

1236 

1237 ''' 

1238 return SourceGrid(base=self, variables=variables) 

1239 

1240 def base_key(self): 

1241 ''' 

1242 Get key to decide about source discretization / GF stack sharing. 

1243 

1244 When two source models differ only in amplitude and origin time, the 

1245 discretization and the GF stacking can be done only once for a unit 

1246 amplitude and a zero origin time and the amplitude and origin times of 

1247 the seismograms can be applied during post-processing of the synthetic 

1248 seismogram. 

1249 

1250 For any derived parameterized source model, this method is called to 

1251 decide if discretization and stacking of the source should be shared. 

1252 When two source models return an equal vector of values discretization 

1253 is shared. 

1254 ''' 

1255 return (self.depth, self.lat, self.north_shift, 

1256 self.lon, self.east_shift, self.time, type(self).__name__) + \ 

1257 self.effective_stf_pre().base_key() 

1258 

1259 def get_factor(self): 

1260 ''' 

1261 Get the scaling factor to be applied during post-processing. 

1262 

1263 Discretization of the base seismogram is usually done for a unit 

1264 amplitude, because a common factor can be efficiently multiplied to 

1265 final seismograms. This eliminates to do repeat the stacking when 

1266 creating seismograms for a series of source models only differing in 

1267 amplitude. 

1268 

1269 This method should return the scaling factor to apply in the 

1270 post-processing (often this is simply the scalar moment of the source). 

1271 ''' 

1272 

1273 return 1.0 

1274 

1275 def effective_stf_pre(self): 

1276 ''' 

1277 Return the STF applied before stacking of the Green's functions. 

1278 

1279 This STF is used during discretization of the parameterized source 

1280 models, i.e. to produce a temporal distribution of point sources. 

1281 

1282 Handling of the STF before stacking of the GFs is less efficient but 

1283 allows to use different source time functions for different parts of 

1284 the source. 

1285 ''' 

1286 

1287 if self.stf is not None and self.stf_mode == 'pre': 

1288 return self.stf 

1289 else: 

1290 return g_unit_pulse 

1291 

1292 def effective_stf_post(self): 

1293 ''' 

1294 Return the STF applied after stacking of the Green's fuctions. 

1295 

1296 This STF is used in the post-processing of the synthetic seismograms. 

1297 

1298 Handling of the STF after stacking of the GFs is usually more efficient 

1299 but is only possible when a common STF is used for all subsources. 

1300 ''' 

1301 

1302 if self.stf is not None and self.stf_mode == 'post': 

1303 return self.stf 

1304 else: 

1305 return g_unit_pulse 

1306 

1307 def _dparams_base(self): 

1308 return dict(times=arr(self.time), 

1309 lat=self.lat, lon=self.lon, 

1310 north_shifts=arr(self.north_shift), 

1311 east_shifts=arr(self.east_shift), 

1312 depths=arr(self.depth)) 

1313 

1314 def _hash(self): 

1315 sha = sha1() 

1316 for k in self.base_key(): 

1317 sha.update(str(k).encode()) 

1318 return sha.hexdigest() 

1319 

1320 def _dparams_base_repeated(self, times): 

1321 if times is None: 

1322 return self._dparams_base() 

1323 

1324 nt = times.size 

1325 north_shifts = num.repeat(self.north_shift, nt) 

1326 east_shifts = num.repeat(self.east_shift, nt) 

1327 depths = num.repeat(self.depth, nt) 

1328 return dict(times=times, 

1329 lat=self.lat, lon=self.lon, 

1330 north_shifts=north_shifts, 

1331 east_shifts=east_shifts, 

1332 depths=depths) 

1333 

1334 def pyrocko_event(self, store=None, target=None, **kwargs): 

1335 duration = None 

1336 if self.stf: 

1337 duration = self.stf.effective_duration 

1338 

1339 return model.Event( 

1340 lat=self.lat, 

1341 lon=self.lon, 

1342 north_shift=self.north_shift, 

1343 east_shift=self.east_shift, 

1344 time=self.time, 

1345 name=self.name, 

1346 depth=self.depth, 

1347 duration=duration, 

1348 **kwargs) 

1349 

1350 def geometry(self, **kwargs): 

1351 raise NotImplementedError 

1352 

1353 def outline(self, cs='xyz'): 

1354 points = num.atleast_2d(num.zeros([1, 3])) 

1355 

1356 points[:, 0] += self.north_shift 

1357 points[:, 1] += self.east_shift 

1358 points[:, 2] += self.depth 

1359 if cs == 'xyz': 

1360 return points 

1361 elif cs == 'xy': 

1362 return points[:, :2] 

1363 elif cs in ('latlon', 'lonlat'): 

1364 latlon = ne_to_latlon( 

1365 self.lat, self.lon, points[:, 0], points[:, 1]) 

1366 

1367 latlon = num.array(latlon).T 

1368 if cs == 'latlon': 

1369 return latlon 

1370 else: 

1371 return latlon[:, ::-1] 

1372 

1373 @classmethod 

1374 def from_pyrocko_event(cls, ev, **kwargs): 

1375 if ev.depth is None: 

1376 raise ConversionError( 

1377 'Cannot convert event object to source object: ' 

1378 'no depth information available') 

1379 

1380 stf = None 

1381 if ev.duration is not None: 

1382 stf = HalfSinusoidSTF(effective_duration=ev.duration) 

1383 

1384 d = dict( 

1385 name=ev.name, 

1386 time=ev.time, 

1387 lat=ev.lat, 

1388 lon=ev.lon, 

1389 north_shift=ev.north_shift, 

1390 east_shift=ev.east_shift, 

1391 depth=ev.depth, 

1392 stf=stf) 

1393 d.update(kwargs) 

1394 return cls(**d) 

1395 

1396 def get_magnitude(self): 

1397 raise NotImplementedError( 

1398 '%s does not implement get_magnitude()' 

1399 % self.__class__.__name__) 

1400 

1401 

1402class SourceWithMagnitude(Source): 

1403 ''' 

1404 Base class for sources containing a moment magnitude. 

1405 ''' 

1406 

1407 magnitude = Float.T( 

1408 default=6.0, 

1409 help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]') 

1410 

1411 def __init__(self, **kwargs): 

1412 if 'moment' in kwargs: 

1413 mom = kwargs.pop('moment') 

1414 if 'magnitude' not in kwargs: 

1415 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom)) 

1416 

1417 Source.__init__(self, **kwargs) 

1418 

1419 @property 

1420 def moment(self): 

1421 return float(pmt.magnitude_to_moment(self.magnitude)) 

1422 

1423 @moment.setter 

1424 def moment(self, value): 

1425 self.magnitude = float(pmt.moment_to_magnitude(value)) 

1426 

1427 def pyrocko_event(self, store=None, target=None, **kwargs): 

1428 return Source.pyrocko_event( 

1429 self, store, target, 

1430 magnitude=self.magnitude, 

1431 **kwargs) 

1432 

1433 @classmethod 

1434 def from_pyrocko_event(cls, ev, **kwargs): 

1435 d = {} 

1436 if ev.magnitude: 

1437 d.update(magnitude=ev.magnitude) 

1438 

1439 d.update(kwargs) 

1440 return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d) 

1441 

1442 def get_magnitude(self): 

1443 return self.magnitude 

1444 

1445 

1446class DerivedMagnitudeError(ValidationError): 

1447 ''' 

1448 Raised when conversion between magnitude, moment, volume change or 

1449 displacement failed. 

1450 ''' 

1451 pass 

1452 

1453 

1454class SourceWithDerivedMagnitude(Source): 

1455 

1456 class __T(Source.T): 

1457 

1458 def validate_extra(self, val): 

1459 Source.T.validate_extra(self, val) 

1460 val.check_conflicts() 

1461 

1462 def check_conflicts(self): 

1463 ''' 

1464 Check for parameter conflicts. 

1465 

1466 To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError` 

1467 on conflicts. 

1468 ''' 

1469 pass 

1470 

1471 def get_magnitude(self, store=None, target=None): 

1472 raise DerivedMagnitudeError('No magnitude set.') 

1473 

1474 def get_moment(self, store=None, target=None): 

1475 return float(pmt.magnitude_to_moment( 

1476 self.get_magnitude(store, target))) 

1477 

1478 def pyrocko_moment_tensor(self, store=None, target=None): 

1479 raise NotImplementedError( 

1480 '%s does not implement pyrocko_moment_tensor()' 

1481 % self.__class__.__name__) 

1482 

1483 def pyrocko_event(self, store=None, target=None, **kwargs): 

1484 try: 

1485 mt = self.pyrocko_moment_tensor(store, target) 

1486 magnitude = self.get_magnitude() 

1487 except (DerivedMagnitudeError, NotImplementedError): 

1488 mt = None 

1489 magnitude = None 

1490 

1491 return Source.pyrocko_event( 

1492 self, store, target, 

1493 moment_tensor=mt, 

1494 magnitude=magnitude, 

1495 **kwargs) 

1496 

1497 

1498class ExplosionSource(SourceWithDerivedMagnitude): 

1499 ''' 

1500 An isotropic explosion point source. 

1501 ''' 

1502 

1503 magnitude = Float.T( 

1504 optional=True, 

1505 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]') 

1506 

1507 volume_change = Float.T( 

1508 optional=True, 

1509 help='volume change of the explosion/implosion or ' 

1510 'the contracting/extending magmatic source. [m^3]') 

1511 

1512 discretized_source_class = meta.DiscretizedExplosionSource 

1513 

1514 def __init__(self, **kwargs): 

1515 if 'moment' in kwargs: 

1516 mom = kwargs.pop('moment') 

1517 if 'magnitude' not in kwargs: 

1518 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom)) 

1519 

1520 SourceWithDerivedMagnitude.__init__(self, **kwargs) 

1521 

1522 def base_key(self): 

1523 return SourceWithDerivedMagnitude.base_key(self) + \ 

1524 (self.volume_change,) 

1525 

1526 def check_conflicts(self): 

1527 if self.magnitude is not None and self.volume_change is not None: 

1528 raise DerivedMagnitudeError( 

1529 'Magnitude and volume_change are both defined.') 

1530 

1531 def get_magnitude(self, store=None, target=None): 

1532 self.check_conflicts() 

1533 

1534 if self.magnitude is not None: 

1535 return self.magnitude 

1536 

1537 elif self.volume_change is not None: 

1538 moment = self.volume_change * \ 

1539 self.get_moment_to_volume_change_ratio(store, target) 

1540 

1541 return float(pmt.moment_to_magnitude(abs(moment))) 

1542 else: 

1543 return float(pmt.moment_to_magnitude(1.0)) 

1544 

1545 def get_volume_change(self, store=None, target=None): 

1546 self.check_conflicts() 

1547 

1548 if self.volume_change is not None: 

1549 return self.volume_change 

1550 

1551 elif self.magnitude is not None: 

1552 moment = float(pmt.magnitude_to_moment(self.magnitude)) 

1553 return moment / self.get_moment_to_volume_change_ratio( 

1554 store, target) 

1555 

1556 else: 

1557 return 1.0 / self.get_moment_to_volume_change_ratio(store) 

1558 

1559 def get_moment_to_volume_change_ratio(self, store, target=None): 

1560 if store is None: 

1561 raise DerivedMagnitudeError( 

1562 'Need earth model to convert between volume change and ' 

1563 'magnitude.') 

1564 

1565 points = num.array( 

1566 [[self.north_shift, self.east_shift, self.depth]], dtype=float) 

1567 

1568 interpolation = target.interpolation if target else 'multilinear' 

1569 try: 

1570 shear_moduli = store.config.get_shear_moduli( 

1571 self.lat, self.lon, 

1572 points=points, 

1573 interpolation=interpolation)[0] 

1574 

1575 bulk_moduli = store.config.get_bulk_moduli( 

1576 self.lat, self.lon, 

1577 points=points, 

1578 interpolation=interpolation)[0] 

1579 except meta.OutOfBounds: 

1580 raise DerivedMagnitudeError( 

1581 'Could not get shear modulus at source position.') 

1582 

1583 return float(2. * shear_moduli + bulk_moduli) 

1584 

1585 def get_factor(self): 

1586 return 1.0 

1587 

1588 def discretize_basesource(self, store, target=None): 

1589 times, amplitudes = self.effective_stf_pre().discretize_t( 

1590 store.config.deltat, self.time) 

1591 

1592 amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.) 

1593 

1594 if self.volume_change is not None: 

1595 if self.volume_change < 0.: 

1596 amplitudes *= -1 

1597 

1598 return meta.DiscretizedExplosionSource( 

1599 m0s=amplitudes, 

1600 **self._dparams_base_repeated(times)) 

1601 

1602 def pyrocko_moment_tensor(self, store=None, target=None): 

1603 a = self.get_moment(store, target) * math.sqrt(2. / 3.) 

1604 return pmt.MomentTensor(m=pmt.symmat6(a, a, a, 0., 0., 0.)) 

1605 

1606 

1607class RectangularExplosionSource(ExplosionSource): 

1608 ''' 

1609 Rectangular or line explosion source. 

1610 ''' 

1611 

1612 discretized_source_class = meta.DiscretizedExplosionSource 

1613 

1614 strike = Float.T( 

1615 default=0.0, 

1616 help='strike direction in [deg], measured clockwise from north') 

1617 

1618 dip = Float.T( 

1619 default=90.0, 

1620 help='dip angle in [deg], measured downward from horizontal') 

1621 

1622 length = Float.T( 

1623 default=0., 

1624 help='length of rectangular source area [m]') 

1625 

1626 width = Float.T( 

1627 default=0., 

1628 help='width of rectangular source area [m]') 

1629 

1630 anchor = StringChoice.T( 

1631 choices=['top', 'top_left', 'top_right', 'center', 'bottom', 

1632 'bottom_left', 'bottom_right'], 

1633 default='center', 

1634 optional=True, 

1635 help='Anchor point for positioning the plane, can be: top, center or' 

1636 'bottom and also top_left, top_right,bottom_left,' 

1637 'bottom_right, center_left and center right') 

1638 

1639 nucleation_x = Float.T( 

1640 optional=True, 

1641 help='horizontal position of rupture nucleation in normalized fault ' 

1642 'plane coordinates (-1 = left edge, +1 = right edge)') 

1643 

1644 nucleation_y = Float.T( 

1645 optional=True, 

1646 help='down-dip position of rupture nucleation in normalized fault ' 

1647 'plane coordinates (-1 = upper edge, +1 = lower edge)') 

1648 

1649 velocity = Float.T( 

1650 default=3500., 

1651 help='speed of explosion front [m/s]') 

1652 

1653 aggressive_oversampling = Bool.T( 

1654 default=False, 

1655 help='Aggressive oversampling for basesource discretization. ' 

1656 "When using 'multilinear' interpolation oversampling has" 

1657 ' practically no effect.') 

1658 

1659 def base_key(self): 

1660 return Source.base_key(self) + (self.strike, self.dip, self.length, 

1661 self.width, self.nucleation_x, 

1662 self.nucleation_y, self.velocity, 

1663 self.anchor) 

1664 

1665 def discretize_basesource(self, store, target=None): 

1666 

1667 if self.nucleation_x is not None: 

1668 nucx = self.nucleation_x * 0.5 * self.length 

1669 else: 

1670 nucx = None 

1671 

1672 if self.nucleation_y is not None: 

1673 nucy = self.nucleation_y * 0.5 * self.width 

1674 else: 

1675 nucy = None 

1676 

1677 stf = self.effective_stf_pre() 

1678 

1679 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source( 

1680 store.config.deltas, store.config.deltat, 

1681 self.time, self.north_shift, self.east_shift, self.depth, 

1682 self.strike, self.dip, self.length, self.width, self.anchor, 

1683 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy) 

1684 

1685 amplitudes /= num.sum(amplitudes) 

1686 amplitudes *= self.get_moment(store, target) 

1687 

1688 return meta.DiscretizedExplosionSource( 

1689 lat=self.lat, 

1690 lon=self.lon, 

1691 times=times, 

1692 north_shifts=points[:, 0], 

1693 east_shifts=points[:, 1], 

1694 depths=points[:, 2], 

1695 m0s=amplitudes) 

1696 

1697 def outline(self, cs='xyz'): 

1698 points = outline_rect_source(self.strike, self.dip, self.length, 

1699 self.width, self.anchor) 

1700 

1701 points[:, 0] += self.north_shift 

1702 points[:, 1] += self.east_shift 

1703 points[:, 2] += self.depth 

1704 if cs == 'xyz': 

1705 return points 

1706 elif cs == 'xy': 

1707 return points[:, :2] 

1708 elif cs in ('latlon', 'lonlat'): 

1709 latlon = ne_to_latlon( 

1710 self.lat, self.lon, points[:, 0], points[:, 1]) 

1711 

1712 latlon = num.array(latlon).T 

1713 if cs == 'latlon': 

1714 return latlon 

1715 else: 

1716 return latlon[:, ::-1] 

1717 

1718 def get_nucleation_abs_coord(self, cs='xy'): 

1719 

1720 if self.nucleation_x is None: 

1721 return None, None 

1722 

1723 coords = from_plane_coords(self.strike, self.dip, self.length, 

1724 self.width, self.depth, self.nucleation_x, 

1725 self.nucleation_y, lat=self.lat, 

1726 lon=self.lon, north_shift=self.north_shift, 

1727 east_shift=self.east_shift, cs=cs) 

1728 return coords 

1729 

1730 

1731class DCSource(SourceWithMagnitude): 

1732 ''' 

1733 A double-couple point source. 

1734 ''' 

1735 

1736 strike = Float.T( 

1737 default=0.0, 

1738 help='strike direction in [deg], measured clockwise from north') 

1739 

1740 dip = Float.T( 

1741 default=90.0, 

1742 help='dip angle in [deg], measured downward from horizontal') 

1743 

1744 rake = Float.T( 

1745 default=0.0, 

1746 help='rake angle in [deg], ' 

1747 'measured counter-clockwise from right-horizontal ' 

1748 'in on-plane view') 

1749 

1750 discretized_source_class = meta.DiscretizedMTSource 

1751 

1752 def base_key(self): 

1753 return Source.base_key(self) + (self.strike, self.dip, self.rake) 

1754 

1755 def get_factor(self): 

1756 return float(pmt.magnitude_to_moment(self.magnitude)) 

1757 

1758 def discretize_basesource(self, store, target=None): 

1759 mot = pmt.MomentTensor( 

1760 strike=self.strike, dip=self.dip, rake=self.rake) 

1761 

1762 times, amplitudes = self.effective_stf_pre().discretize_t( 

1763 store.config.deltat, self.time) 

1764 return meta.DiscretizedMTSource( 

1765 m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis], 

1766 **self._dparams_base_repeated(times)) 

1767 

1768 def pyrocko_moment_tensor(self, store=None, target=None): 

1769 return pmt.MomentTensor( 

1770 strike=self.strike, 

1771 dip=self.dip, 

1772 rake=self.rake, 

1773 scalar_moment=self.moment) 

1774 

1775 def pyrocko_event(self, store=None, target=None, **kwargs): 

1776 return SourceWithMagnitude.pyrocko_event( 

1777 self, store, target, 

1778 moment_tensor=self.pyrocko_moment_tensor(store, target), 

1779 **kwargs) 

1780 

1781 @classmethod 

1782 def from_pyrocko_event(cls, ev, **kwargs): 

1783 d = {} 

1784 mt = ev.moment_tensor 

1785 if mt: 

1786 (strike, dip, rake), _ = mt.both_strike_dip_rake() 

1787 d.update( 

1788 strike=float(strike), 

1789 dip=float(dip), 

1790 rake=float(rake), 

1791 magnitude=float(mt.moment_magnitude())) 

1792 

1793 d.update(kwargs) 

1794 return super(DCSource, cls).from_pyrocko_event(ev, **d) 

1795 

1796 

1797class CLVDSource(SourceWithMagnitude): 

1798 ''' 

1799 A pure CLVD point source. 

1800 ''' 

1801 

1802 discretized_source_class = meta.DiscretizedMTSource 

1803 

1804 azimuth = Float.T( 

1805 default=0.0, 

1806 help='azimuth direction of largest dipole, clockwise from north [deg]') 

1807 

1808 dip = Float.T( 

1809 default=90., 

1810 help='dip direction of largest dipole, downward from horizontal [deg]') 

1811 

1812 def base_key(self): 

1813 return Source.base_key(self) + (self.azimuth, self.dip) 

1814 

1815 def get_factor(self): 

1816 return float(pmt.magnitude_to_moment(self.magnitude)) 

1817 

1818 @property 

1819 def m6(self): 

1820 a = math.sqrt(4. / 3.) * self.get_factor() 

1821 m = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.) 

1822 rotmat1 = pmt.euler_to_matrix( 

1823 d2r * (self.dip - 90.), 

1824 d2r * (self.azimuth - 90.), 

1825 0.) 

1826 m = num.dot(rotmat1.T, num.dot(m, rotmat1)) 

1827 return pmt.to6(m) 

1828 

1829 @property 

1830 def m6_astuple(self): 

1831 return tuple(self.m6.tolist()) 

1832 

1833 def discretize_basesource(self, store, target=None): 

1834 factor = self.get_factor() 

1835 times, amplitudes = self.effective_stf_pre().discretize_t( 

1836 store.config.deltat, self.time) 

1837 return meta.DiscretizedMTSource( 

1838 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor, 

1839 **self._dparams_base_repeated(times)) 

1840 

1841 def pyrocko_moment_tensor(self, store=None, target=None): 

1842 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple)) 

1843 

1844 def pyrocko_event(self, store=None, target=None, **kwargs): 

1845 mt = self.pyrocko_moment_tensor(store, target) 

1846 return Source.pyrocko_event( 

1847 self, store, target, 

1848 moment_tensor=self.pyrocko_moment_tensor(store, target), 

1849 magnitude=float(mt.moment_magnitude()), 

1850 **kwargs) 

1851 

1852 

1853class VLVDSource(SourceWithDerivedMagnitude): 

1854 ''' 

1855 Volumetric linear vector dipole source. 

1856 

1857 This source is a parameterization for a restricted moment tensor point 

1858 source, useful to represent dyke or sill like inflation or deflation 

1859 sources. The restriction is such that the moment tensor is rotational 

1860 symmetric. It can be represented by a superposition of a linear vector 

1861 dipole (here we use a CLVD for convenience) and an isotropic component. The 

1862 restricted moment tensor has 4 degrees of freedom: 2 independent 

1863 eigenvalues and 2 rotation angles orienting the the symmetry axis. 

1864 

1865 In this parameterization, the isotropic component is controlled by 

1866 ``volume_change``. To define the moment tensor, it must be converted to the 

1867 scalar moment of the the MT's isotropic component. For the conversion, the 

1868 shear modulus at the source's position must be known. This value is 

1869 extracted from the earth model defined in the GF store in use. 

1870 

1871 The CLVD part by controlled by its scalar moment :math:`M_0`: 

1872 ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a 

1873 positiv or negativ CLVD (the sign of the largest eigenvalue). 

1874 ''' 

1875 

1876 discretized_source_class = meta.DiscretizedMTSource 

1877 

1878 azimuth = Float.T( 

1879 default=0.0, 

1880 help='azimuth direction of symmetry axis, clockwise from north [deg].') 

1881 

1882 dip = Float.T( 

1883 default=90., 

1884 help='dip direction of symmetry axis, downward from horizontal [deg].') 

1885 

1886 volume_change = Float.T( 

1887 default=0., 

1888 help='volume change of the inflation/deflation [m^3].') 

1889 

1890 clvd_moment = Float.T( 

1891 default=0., 

1892 help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign ' 

1893 'controls the sign of the CLVD (the sign of its largest ' 

1894 'eigenvalue).') 

1895 

1896 def get_moment_to_volume_change_ratio(self, store, target): 

1897 if store is None or target is None: 

1898 raise DerivedMagnitudeError( 

1899 'Need earth model to convert between volume change and ' 

1900 'magnitude.') 

1901 

1902 points = num.array( 

1903 [[self.north_shift, self.east_shift, self.depth]], dtype=float) 

1904 

1905 try: 

1906 shear_moduli = store.config.get_shear_moduli( 

1907 self.lat, self.lon, 

1908 points=points, 

1909 interpolation=target.interpolation)[0] 

1910 

1911 bulk_moduli = store.config.get_bulk_moduli( 

1912 self.lat, self.lon, 

1913 points=points, 

1914 interpolation=target.interpolation)[0] 

1915 except meta.OutOfBounds: 

1916 raise DerivedMagnitudeError( 

1917 'Could not get shear modulus at source position.') 

1918 

1919 return float(2. * shear_moduli + bulk_moduli) 

1920 

1921 def base_key(self): 

1922 return Source.base_key(self) + \ 

1923 (self.azimuth, self.dip, self.volume_change, self.clvd_moment) 

1924 

1925 def get_magnitude(self, store=None, target=None): 

1926 mt = self.pyrocko_moment_tensor(store, target) 

1927 return float(pmt.moment_to_magnitude(mt.moment)) 

1928 

1929 def get_m6(self, store, target): 

1930 a = math.sqrt(4. / 3.) * self.clvd_moment 

1931 m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.) 

1932 

1933 rotmat1 = pmt.euler_to_matrix( 

1934 d2r * (self.dip - 90.), 

1935 d2r * (self.azimuth - 90.), 

1936 0.) 

1937 m_clvd = num.dot(rotmat1.T, num.dot(m_clvd, rotmat1)) 

1938 

1939 m_iso = self.volume_change * \ 

1940 self.get_moment_to_volume_change_ratio(store, target) 

1941 

1942 m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0., 

1943 0., 0.,) * math.sqrt(2. / 3) 

1944 

1945 m = pmt.to6(m_clvd) + pmt.to6(m_iso) 

1946 return m 

1947 

1948 def get_moment(self, store=None, target=None): 

1949 return float(pmt.magnitude_to_moment( 

1950 self.get_magnitude(store, target))) 

1951 

1952 def get_m6_astuple(self, store, target): 

1953 m6 = self.get_m6(store, target) 

1954 return tuple(m6.tolist()) 

1955 

1956 def discretize_basesource(self, store, target=None): 

1957 times, amplitudes = self.effective_stf_pre().discretize_t( 

1958 store.config.deltat, self.time) 

1959 

1960 m6 = self.get_m6(store, target) 

1961 m6 *= amplitudes / self.get_factor() 

1962 

1963 return meta.DiscretizedMTSource( 

1964 m6s=m6[num.newaxis, :], 

1965 **self._dparams_base_repeated(times)) 

1966 

1967 def pyrocko_moment_tensor(self, store=None, target=None): 

1968 m6_astuple = self.get_m6_astuple(store, target) 

1969 return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple)) 

1970 

1971 

1972class MTSource(Source): 

1973 ''' 

1974 A moment tensor point source. 

1975 ''' 

1976 

1977 discretized_source_class = meta.DiscretizedMTSource 

1978 

1979 mnn = Float.T( 

1980 default=1., 

1981 help='north-north component of moment tensor in [Nm]') 

1982 

1983 mee = Float.T( 

1984 default=1., 

1985 help='east-east component of moment tensor in [Nm]') 

1986 

1987 mdd = Float.T( 

1988 default=1., 

1989 help='down-down component of moment tensor in [Nm]') 

1990 

1991 mne = Float.T( 

1992 default=0., 

1993 help='north-east component of moment tensor in [Nm]') 

1994 

1995 mnd = Float.T( 

1996 default=0., 

1997 help='north-down component of moment tensor in [Nm]') 

1998 

1999 med = Float.T( 

2000 default=0., 

2001 help='east-down component of moment tensor in [Nm]') 

2002 

2003 def __init__(self, **kwargs): 

2004 if 'm6' in kwargs: 

2005 for (k, v) in zip('mnn mee mdd mne mnd med'.split(), 

2006 kwargs.pop('m6')): 

2007 kwargs[k] = float(v) 

2008 

2009 Source.__init__(self, **kwargs) 

2010 

2011 @property 

2012 def m6(self): 

2013 return num.array(self.m6_astuple) 

2014 

2015 @property 

2016 def m6_astuple(self): 

2017 return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med) 

2018 

2019 @m6.setter 

2020 def m6(self, value): 

2021 self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value 

2022 

2023 def base_key(self): 

2024 return Source.base_key(self) + self.m6_astuple 

2025 

2026 def discretize_basesource(self, store, target=None): 

2027 times, amplitudes = self.effective_stf_pre().discretize_t( 

2028 store.config.deltat, self.time) 

2029 return meta.DiscretizedMTSource( 

2030 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis], 

2031 **self._dparams_base_repeated(times)) 

2032 

2033 def get_magnitude(self, store=None, target=None): 

2034 m6 = self.m6 

2035 return pmt.moment_to_magnitude( 

2036 math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) / 

2037 math.sqrt(2.)) 

2038 

2039 def pyrocko_moment_tensor(self, store=None, target=None): 

2040 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple)) 

2041 

2042 def pyrocko_event(self, store=None, target=None, **kwargs): 

2043 mt = self.pyrocko_moment_tensor(store, target) 

2044 return Source.pyrocko_event( 

2045 self, store, target, 

2046 moment_tensor=self.pyrocko_moment_tensor(store, target), 

2047 magnitude=float(mt.moment_magnitude()), 

2048 **kwargs) 

2049 

2050 @classmethod 

2051 def from_pyrocko_event(cls, ev, **kwargs): 

2052 d = {} 

2053 mt = ev.moment_tensor 

2054 if mt: 

2055 d.update(m6=tuple(map(float, mt.m6()))) 

2056 else: 

2057 if ev.magnitude is not None: 

2058 mom = pmt.magnitude_to_moment(ev.magnitude) 

2059 v = math.sqrt(2. / 3.) * mom 

2060 d.update(m6=(v, v, v, 0., 0., 0.)) 

2061 

2062 d.update(kwargs) 

2063 return super(MTSource, cls).from_pyrocko_event(ev, **d) 

2064 

2065 

2066map_anchor = { 

2067 'center': (0.0, 0.0), 

2068 'center_left': (-1.0, 0.0), 

2069 'center_right': (1.0, 0.0), 

2070 'top': (0.0, -1.0), 

2071 'top_left': (-1.0, -1.0), 

2072 'top_right': (1.0, -1.0), 

2073 'bottom': (0.0, 1.0), 

2074 'bottom_left': (-1.0, 1.0), 

2075 'bottom_right': (1.0, 1.0)} 

2076 

2077 

2078class RectangularSource(SourceWithDerivedMagnitude): 

2079 ''' 

2080 Classical Haskell source model modified for bilateral rupture. 

2081 ''' 

2082 

2083 discretized_source_class = meta.DiscretizedMTSource 

2084 

2085 magnitude = Float.T( 

2086 optional=True, 

2087 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]') 

2088 

2089 strike = Float.T( 

2090 default=0.0, 

2091 help='strike direction in [deg], measured clockwise from north') 

2092 

2093 dip = Float.T( 

2094 default=90.0, 

2095 help='dip angle in [deg], measured downward from horizontal') 

2096 

2097 rake = Float.T( 

2098 default=0.0, 

2099 help='rake angle in [deg], ' 

2100 'measured counter-clockwise from right-horizontal ' 

2101 'in on-plane view') 

2102 

2103 length = Float.T( 

2104 default=0., 

2105 help='length of rectangular source area [m]') 

2106 

2107 width = Float.T( 

2108 default=0., 

2109 help='width of rectangular source area [m]') 

2110 

2111 anchor = StringChoice.T( 

2112 choices=['top', 'top_left', 'top_right', 'center', 'bottom', 

2113 'bottom_left', 'bottom_right'], 

2114 default='center', 

2115 optional=True, 

2116 help='Anchor point for positioning the plane, can be: ``top, center ' 

2117 'bottom, top_left, top_right,bottom_left,' 

2118 'bottom_right, center_left, center right``.') 

2119 

2120 nucleation_x = Float.T( 

2121 optional=True, 

2122 help='horizontal position of rupture nucleation in normalized fault ' 

2123 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge)') 

2124 

2125 nucleation_y = Float.T( 

2126 optional=True, 

2127 help='down-dip position of rupture nucleation in normalized fault ' 

2128 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge)') 

2129 

2130 velocity = Float.T( 

2131 default=3500., 

2132 help='speed of rupture front [m/s]') 

2133 

2134 slip = Float.T( 

2135 optional=True, 

2136 help='Slip on the rectangular source area [m]') 

2137 

2138 opening_fraction = Float.T( 

2139 default=0., 

2140 help='Determines fraction of slip related to opening. ' 

2141 '(``-1``: pure tensile closing, ' 

2142 '``0``: pure shear, ' 

2143 '``1``: pure tensile opening)') 

2144 

2145 decimation_factor = Int.T( 

2146 optional=True, 

2147 default=1, 

2148 help='Sub-source decimation factor, a larger decimation will' 

2149 ' make the result inaccurate but shorten the necessary' 

2150 ' computation time (use for testing puposes only).') 

2151 

2152 aggressive_oversampling = Bool.T( 

2153 default=False, 

2154 help='Aggressive oversampling for basesource discretization. ' 

2155 "When using 'multilinear' interpolation oversampling has" 

2156 ' practically no effect.') 

2157 

2158 def __init__(self, **kwargs): 

2159 if 'moment' in kwargs: 

2160 mom = kwargs.pop('moment') 

2161 if 'magnitude' not in kwargs: 

2162 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom)) 

2163 

2164 SourceWithDerivedMagnitude.__init__(self, **kwargs) 

2165 

2166 def base_key(self): 

2167 return SourceWithDerivedMagnitude.base_key(self) + ( 

2168 self.magnitude, 

2169 self.slip, 

2170 self.strike, 

2171 self.dip, 

2172 self.rake, 

2173 self.length, 

2174 self.width, 

2175 self.nucleation_x, 

2176 self.nucleation_y, 

2177 self.velocity, 

2178 self.decimation_factor, 

2179 self.anchor) 

2180 

2181 def check_conflicts(self): 

2182 if self.magnitude is not None and self.slip is not None: 

2183 raise DerivedMagnitudeError( 

2184 'Magnitude and slip are both defined.') 

2185 

2186 def get_magnitude(self, store=None, target=None): 

2187 self.check_conflicts() 

2188 if self.magnitude is not None: 

2189 return self.magnitude 

2190 

2191 elif self.slip is not None: 

2192 if None in (store, target): 

2193 raise DerivedMagnitudeError( 

2194 'Magnitude for a rectangular source with slip defined ' 

2195 'can only be derived when earth model and target ' 

2196 'interpolation method are available.') 

2197 

2198 amplitudes = self._discretize(store, target)[2] 

2199 if amplitudes.ndim == 2: 

2200 # CLVD component has no net moment, leave out 

2201 return float(pmt.moment_to_magnitude( 

2202 num.sum(num.abs(amplitudes[0:2, :]).sum()))) 

2203 else: 

2204 return float(pmt.moment_to_magnitude(num.sum(amplitudes))) 

2205 

2206 else: 

2207 return float(pmt.moment_to_magnitude(1.0)) 

2208 

2209 def get_factor(self): 

2210 return 1.0 

2211 

2212 def get_slip_tensile(self): 

2213 return self.slip * self.opening_fraction 

2214 

2215 def get_slip_shear(self): 

2216 return self.slip - abs(self.get_slip_tensile) 

2217 

2218 def _discretize(self, store, target): 

2219 if self.nucleation_x is not None: 

2220 nucx = self.nucleation_x * 0.5 * self.length 

2221 else: 

2222 nucx = None 

2223 

2224 if self.nucleation_y is not None: 

2225 nucy = self.nucleation_y * 0.5 * self.width 

2226 else: 

2227 nucy = None 

2228 

2229 stf = self.effective_stf_pre() 

2230 

2231 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source( 

2232 store.config.deltas, store.config.deltat, 

2233 self.time, self.north_shift, self.east_shift, self.depth, 

2234 self.strike, self.dip, self.length, self.width, self.anchor, 

2235 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy, 

2236 decimation_factor=self.decimation_factor, 

2237 aggressive_oversampling=self.aggressive_oversampling) 

2238 

2239 if self.slip is not None: 

2240 if target is not None: 

2241 interpolation = target.interpolation 

2242 else: 

2243 interpolation = 'nearest_neighbor' 

2244 logger.warning( 

2245 'no target information available, will use ' 

2246 '"nearest_neighbor" interpolation when extracting shear ' 

2247 'modulus from earth model') 

2248 

2249 shear_moduli = store.config.get_shear_moduli( 

2250 self.lat, self.lon, 

2251 points=points, 

2252 interpolation=interpolation) 

2253 

2254 tensile_slip = self.get_slip_tensile() 

2255 shear_slip = self.slip - abs(tensile_slip) 

2256 

2257 amplitudes_total = [shear_moduli * shear_slip] 

2258 if tensile_slip != 0: 

2259 bulk_moduli = store.config.get_bulk_moduli( 

2260 self.lat, self.lon, 

2261 points=points, 

2262 interpolation=interpolation) 

2263 

2264 tensile_iso = bulk_moduli * tensile_slip 

2265 tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip 

2266 

2267 amplitudes_total.extend([tensile_iso, tensile_clvd]) 

2268 

2269 amplitudes_total = num.vstack(amplitudes_total).squeeze() * \ 

2270 amplitudes * dl * dw 

2271 

2272 else: 

2273 # normalization to retain total moment 

2274 amplitudes_norm = amplitudes / num.sum(amplitudes) 

2275 moment = self.get_moment(store, target) 

2276 

2277 amplitudes_total = [ 

2278 amplitudes_norm * moment * (1 - abs(self.opening_fraction))] 

2279 if self.opening_fraction != 0.: 

2280 amplitudes_total.append( 

2281 amplitudes_norm * self.opening_fraction * moment) 

2282 

2283 amplitudes_total = num.vstack(amplitudes_total).squeeze() 

2284 

2285 return points, times, num.atleast_1d(amplitudes_total), dl, dw, nl, nw 

2286 

2287 def discretize_basesource(self, store, target=None): 

2288 

2289 points, times, amplitudes, dl, dw, nl, nw = self._discretize( 

2290 store, target) 

2291 

2292 mot = pmt.MomentTensor( 

2293 strike=self.strike, dip=self.dip, rake=self.rake) 

2294 m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0) 

2295 

2296 if amplitudes.ndim == 1: 

2297 m6s[:, :] *= amplitudes[:, num.newaxis] 

2298 elif amplitudes.ndim == 2: 

2299 # shear MT components 

2300 rotmat1 = pmt.euler_to_matrix( 

2301 d2r * self.dip, d2r * self.strike, d2r * -self.rake) 

2302 m6s[:, :] *= amplitudes[0, :][:, num.newaxis] 

2303 

2304 if amplitudes.shape[0] == 2: 

2305 # tensile MT components - moment/magnitude input 

2306 tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.) 

2307 rot_tensile = pmt.to6( 

2308 num.dot(rotmat1.T, num.dot(tensile, rotmat1))) 

2309 

2310 m6s_tensile = rot_tensile[ 

2311 num.newaxis, :] * amplitudes[1, :][:, num.newaxis] 

2312 m6s += m6s_tensile 

2313 

2314 elif amplitudes.shape[0] == 3: 

2315 # tensile MT components - slip input 

2316 iso = pmt.symmat6(1., 1., 1., 0., 0., 0.) 

2317 clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.) 

2318 

2319 rot_iso = pmt.to6( 

2320 num.dot(rotmat1.T, num.dot(iso, rotmat1))) 

2321 rot_clvd = pmt.to6( 

2322 num.dot(rotmat1.T, num.dot(clvd, rotmat1))) 

2323 

2324 m6s_iso = rot_iso[ 

2325 num.newaxis, :] * amplitudes[1, :][:, num.newaxis] 

2326 m6s_clvd = rot_clvd[ 

2327 num.newaxis, :] * amplitudes[2, :][:, num.newaxis] 

2328 m6s += m6s_iso + m6s_clvd 

2329 else: 

2330 raise ValueError('Unknwown amplitudes shape!') 

2331 else: 

2332 raise ValueError( 

2333 'Unexpected dimension of {}'.format(amplitudes.ndim)) 

2334 

2335 ds = meta.DiscretizedMTSource( 

2336 lat=self.lat, 

2337 lon=self.lon, 

2338 times=times, 

2339 north_shifts=points[:, 0], 

2340 east_shifts=points[:, 1], 

2341 depths=points[:, 2], 

2342 m6s=m6s, 

2343 dl=dl, 

2344 dw=dw, 

2345 nl=nl, 

2346 nw=nw) 

2347 

2348 return ds 

2349 

2350 def xy_to_coord(self, x, y, cs='xyz'): 

2351 ln, wd = self.length, self.width 

2352 strike, dip = self.strike, self.dip 

2353 

2354 def array_check(variable): 

2355 if not isinstance(variable, num.ndarray): 

2356 return num.array(variable) 

2357 else: 

2358 return variable 

2359 

2360 x, y = array_check(x), array_check(y) 

2361 

2362 if x.shape[0] != y.shape[0]: 

2363 raise ValueError('Shapes of x and y mismatch') 

2364 

2365 x, y = x * 0.5 * ln, y * 0.5 * wd 

2366 

2367 points = num.hstack(( 

2368 x.reshape(-1, 1), y.reshape(-1, 1), num.zeros((x.shape[0], 1)))) 

2369 

2370 anch_x, anch_y = map_anchor[self.anchor] 

2371 points[:, 0] -= anch_x * 0.5 * ln 

2372 points[:, 1] -= anch_y * 0.5 * wd 

2373 

2374 rotmat = num.asarray( 

2375 pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)) 

2376 

2377 points_rot = num.dot(rotmat.T, points.T).T 

2378 

2379 points_rot[:, 0] += self.north_shift 

2380 points_rot[:, 1] += self.east_shift 

2381 points_rot[:, 2] += self.depth 

2382 

2383 if cs == 'xyz': 

2384 return points_rot 

2385 elif cs == 'xy': 

2386 return points_rot[:, :2] 

2387 elif cs in ('latlon', 'lonlat', 'latlondepth'): 

2388 latlon = ne_to_latlon( 

2389 self.lat, self.lon, points_rot[:, 0], points_rot[:, 1]) 

2390 latlon = num.array(latlon).T 

2391 if cs == 'latlon': 

2392 return latlon 

2393 elif cs == 'lonlat': 

2394 return latlon[:, ::-1] 

2395 else: 

2396 return num.concatenate( 

2397 (latlon, points_rot[:, 2].reshape((len(points_rot), 1))), 

2398 axis=1) 

2399 

2400 def outline(self, cs='xyz'): 

2401 x = num.array([-1., 1., 1., -1., -1.]) 

2402 y = num.array([-1., -1., 1., 1., -1.]) 

2403 

2404 return self.xy_to_coord(x, y, cs=cs) 

2405 

2406 def points_on_source(self, cs='xyz', **kwargs): 

2407 

2408 points = points_on_rect_source( 

2409 self.strike, self.dip, self.length, self.width, 

2410 self.anchor, **kwargs) 

2411 

2412 points[:, 0] += self.north_shift 

2413 points[:, 1] += self.east_shift 

2414 points[:, 2] += self.depth 

2415 if cs == 'xyz': 

2416 return points 

2417 elif cs == 'xy': 

2418 return points[:, :2] 

2419 elif cs in ('latlon', 'lonlat', 'latlondepth'): 

2420 latlon = ne_to_latlon( 

2421 self.lat, self.lon, points[:, 0], points[:, 1]) 

2422 

2423 latlon = num.array(latlon).T 

2424 if cs == 'latlon': 

2425 return latlon 

2426 elif cs == 'lonlat': 

2427 return latlon[:, ::-1] 

2428 else: 

2429 return num.concatenate( 

2430 (latlon, points[:, 2].reshape((len(points), 1))), 

2431 axis=1) 

2432 

2433 def geometry(self, *args, **kwargs): 

2434 from pyrocko.model import Geometry 

2435 

2436 ds = self.discretize_basesource(*args, **kwargs) 

2437 nx, ny = ds.nl, ds.nw 

2438 

2439 def patch_outlines_xy(nx, ny): 

2440 points = num.zeros((nx * ny, 2)) 

2441 points[:, 0] = num.tile(num.linspace(-1., 1., nx), ny) 

2442 points[:, 1] = num.repeat(num.linspace(-1., 1., ny), nx) 

2443 

2444 return points 

2445 

2446 points_ds = patch_outlines_xy(nx + 1, ny + 1) 

2447 npoints = (nx + 1) * (ny + 1) 

2448 

2449 vertices = num.hstack(( 

2450 num.ones((npoints, 1)) * self.lat, 

2451 num.ones((npoints, 1)) * self.lon, 

2452 self.xy_to_coord(points_ds[:, 0], points_ds[:, 1], cs='xyz'))) 

2453 

2454 faces = num.array([[ 

2455 iy * (nx + 1) + ix, 

2456 iy * (nx + 1) + ix + 1, 

2457 (iy + 1) * (nx + 1) + ix + 1, 

2458 (iy + 1) * (nx + 1) + ix, 

2459 iy * (nx + 1) + ix] 

2460 for iy in range(ny) for ix in range(nx)]) 

2461 

2462 xyz = self.outline('xyz') 

2463 latlon = num.ones((5, 2)) * num.array([self.lat, self.lon]) 

2464 patchverts = num.hstack((latlon, xyz)) 

2465 

2466 geom = Geometry() 

2467 geom.setup(vertices, faces) 

2468 geom.set_outlines([patchverts]) 

2469 

2470 if self.stf: 

2471 geom.times = num.unique(ds.times) 

2472 

2473 if self.nucleation_x is not None and self.nucleation_y is not None: 

2474 geom.add_property('t_arrival', ds.times) 

2475 

2476 geom.add_property( 

2477 'moment', ds.moments().reshape(ds.nl*ds.nw, -1)) 

2478 

2479 geom.add_property( 

2480 'slip', num.ones_like(ds.times) * self.slip) 

2481 

2482 return geom 

2483 

2484 def get_nucleation_abs_coord(self, cs='xy'): 

2485 

2486 if self.nucleation_x is None: 

2487 return None, None 

2488 

2489 coords = from_plane_coords(self.strike, self.dip, self.length, 

2490 self.width, self.depth, self.nucleation_x, 

2491 self.nucleation_y, lat=self.lat, 

2492 lon=self.lon, north_shift=self.north_shift, 

2493 east_shift=self.east_shift, cs=cs) 

2494 return coords 

2495 

2496 def pyrocko_moment_tensor(self, store=None, target=None): 

2497 return pmt.MomentTensor( 

2498 strike=self.strike, 

2499 dip=self.dip, 

2500 rake=self.rake, 

2501 scalar_moment=self.get_moment(store, target)) 

2502 

2503 def pyrocko_event(self, store=None, target=None, **kwargs): 

2504 return SourceWithDerivedMagnitude.pyrocko_event( 

2505 self, store, target, 

2506 **kwargs) 

2507 

2508 @classmethod 

2509 def from_pyrocko_event(cls, ev, **kwargs): 

2510 d = {} 

2511 mt = ev.moment_tensor 

2512 if mt: 

2513 (strike, dip, rake), _ = mt.both_strike_dip_rake() 

2514 d.update( 

2515 strike=float(strike), 

2516 dip=float(dip), 

2517 rake=float(rake), 

2518 magnitude=float(mt.moment_magnitude())) 

2519 

2520 d.update(kwargs) 

2521 return super(RectangularSource, cls).from_pyrocko_event(ev, **d) 

2522 

2523 

2524class PseudoDynamicRupture(SourceWithDerivedMagnitude): 

2525 ''' 

2526 Combined Eikonal and Okada quasi-dynamic rupture model. 

2527 

2528 Details are described in :doc:`/topics/pseudo-dynamic-rupture`. 

2529 Note: attribute `stf` is not used so far, but kept for future applications. 

2530 ''' 

2531 

2532 discretized_source_class = meta.DiscretizedMTSource 

2533 

2534 strike = Float.T( 

2535 default=0.0, 

2536 help='Strike direction in [deg], measured clockwise from north.') 

2537 

2538 dip = Float.T( 

2539 default=0.0, 

2540 help='Dip angle in [deg], measured downward from horizontal.') 

2541 

2542 length = Float.T( 

2543 default=10. * km, 

2544 help='Length of rectangular source area in [m].') 

2545 

2546 width = Float.T( 

2547 default=5. * km, 

2548 help='Width of rectangular source area in [m].') 

2549 

2550 anchor = StringChoice.T( 

2551 choices=['top', 'top_left', 'top_right', 'center', 'bottom', 

2552 'bottom_left', 'bottom_right'], 

2553 default='center', 

2554 optional=True, 

2555 help='Anchor point for positioning the plane, can be: ``top, center, ' 

2556 'bottom, top_left, top_right, bottom_left, ' 

2557 'bottom_right, center_left, center_right``.') 

2558 

2559 nucleation_x__ = Array.T( 

2560 default=num.array([0.]), 

2561 dtype=num.float64, 

2562 serialize_as='list', 

2563 help='Horizontal position of rupture nucleation in normalized fault ' 

2564 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge).') 

2565 

2566 nucleation_y__ = Array.T( 

2567 default=num.array([0.]), 

2568 dtype=num.float64, 

2569 serialize_as='list', 

2570 help='Down-dip position of rupture nucleation in normalized fault ' 

2571 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge).') 

2572 

2573 nucleation_time__ = Array.T( 

2574 optional=True, 

2575 help='Time in [s] after origin, when nucleation points defined by ' 

2576 '``nucleation_x`` and ``nucleation_y`` rupture.', 

2577 dtype=num.float64, 

2578 serialize_as='list') 

2579 

2580 gamma = Float.T( 

2581 default=0.8, 

2582 help='Scaling factor between rupture velocity and S-wave velocity: ' 

2583 r':math:`v_r = \gamma * v_s`.') 

2584 

2585 nx = Int.T( 

2586 default=2, 

2587 help='Number of discrete source patches in x direction (along ' 

2588 'strike).') 

2589 

2590 ny = Int.T( 

2591 default=2, 

2592 help='Number of discrete source patches in y direction (down dip).') 

2593 

2594 slip = Float.T( 

2595 optional=True, 

2596 help='Maximum slip of the rectangular source [m]. ' 

2597 'Setting the slip the tractions/stress field ' 

2598 'will be normalized to accomodate the desired maximum slip.') 

2599 

2600 rake = Float.T( 

2601 optional=True, 

2602 help='Rake angle in [deg], ' 

2603 'measured counter-clockwise from right-horizontal ' 

2604 'in on-plane view. Rake is translated into homogenous tractions ' 

2605 'in strike and up-dip direction. ``rake`` is mutually exclusive ' 

2606 'with tractions parameter.') 

2607 

2608 patches = List.T( 

2609 OkadaSource.T(), 

2610 optional=True, 

2611 help='List of all boundary elements/sub faults/fault patches.') 

2612 

2613 patch_mask__ = Array.T( 

2614 dtype=bool, 

2615 serialize_as='list', 

2616 shape=(None,), 

2617 optional=True, 

2618 help='Mask for all boundary elements/sub faults/fault patches. True ' 

2619 'leaves the patch in the calculation, False excludes the patch.') 

2620 

2621 tractions = TractionField.T( 

2622 optional=True, 

2623 help='Traction field the rupture plane is exposed to. See the ' 

2624 ':py:mod:`pyrocko.gf.tractions` module for more details. ' 

2625 'If ``tractions=None`` and ``rake`` is given' 

2626 ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will' 

2627 ' be used.') 

2628 

2629 coef_mat = Array.T( 

2630 optional=True, 

2631 help='Coefficient matrix linking traction and dislocation field.', 

2632 dtype=num.float64, 

2633 shape=(None, None)) 

2634 

2635 eikonal_decimation = Int.T( 

2636 optional=True, 

2637 default=1, 

2638 help='Sub-source eikonal factor, a smaller eikonal factor will' 

2639 ' increase the accuracy of rupture front calculation but' 

2640 ' increases also the computation time.') 

2641 

2642 decimation_factor = Int.T( 

2643 optional=True, 

2644 default=1, 

2645 help='Sub-source decimation factor, a larger decimation will' 

2646 ' make the result inaccurate but shorten the necessary' 

2647 ' computation time (use for testing puposes only).') 

2648 

2649 nthreads = Int.T( 

2650 optional=True, 

2651 default=1, 

2652 help='Number of threads for Okada forward modelling, ' 

2653 'matrix inversion and calculation of point subsources. ' 

2654 'Note: for small/medium matrices 1 thread is most efficient.') 

2655 

2656 pure_shear = Bool.T( 

2657 optional=True, 

2658 default=False, 

2659 help='Calculate only shear tractions and omit tensile tractions.') 

2660 

2661 smooth_rupture = Bool.T( 

2662 default=True, 

2663 help='Smooth the tractions by weighting partially ruptured' 

2664 ' fault patches.') 

2665 

2666 aggressive_oversampling = Bool.T( 

2667 default=False, 

2668 help='Aggressive oversampling for basesource discretization. ' 

2669 "When using 'multilinear' interpolation oversampling has" 

2670 ' practically no effect.') 

2671 

2672 def __init__(self, **kwargs): 

2673 SourceWithDerivedMagnitude.__init__(self, **kwargs) 

2674 self._interpolators = {} 

2675 self.check_conflicts() 

2676 

2677 @property 

2678 def nucleation_x(self): 

2679 return self.nucleation_x__ 

2680 

2681 @nucleation_x.setter 

2682 def nucleation_x(self, nucleation_x): 

2683 if isinstance(nucleation_x, list): 

2684 nucleation_x = num.array(nucleation_x) 

2685 

2686 elif not isinstance( 

2687 nucleation_x, num.ndarray) and nucleation_x is not None: 

2688 

2689 nucleation_x = num.array([nucleation_x]) 

2690 self.nucleation_x__ = nucleation_x 

2691 

2692 @property 

2693 def nucleation_y(self): 

2694 return self.nucleation_y__ 

2695 

2696 @nucleation_y.setter 

2697 def nucleation_y(self, nucleation_y): 

2698 if isinstance(nucleation_y, list): 

2699 nucleation_y = num.array(nucleation_y) 

2700 

2701 elif not isinstance(nucleation_y, num.ndarray) \ 

2702 and nucleation_y is not None: 

2703 nucleation_y = num.array([nucleation_y]) 

2704 

2705 self.nucleation_y__ = nucleation_y 

2706 

2707 @property 

2708 def nucleation(self): 

2709 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y 

2710 

2711 if (nucl_x is None) or (nucl_y is None): 

2712 return None 

2713 

2714 assert nucl_x.shape[0] == nucl_y.shape[0] 

2715 

2716 return num.concatenate( 

2717 (nucl_x[:, num.newaxis], nucl_y[:, num.newaxis]), axis=1) 

2718 

2719 @nucleation.setter 

2720 def nucleation(self, nucleation): 

2721 if isinstance(nucleation, list): 

2722 nucleation = num.array(nucleation) 

2723 

2724 assert nucleation.shape[1] == 2 

2725 

2726 self.nucleation_x = nucleation[:, 0] 

2727 self.nucleation_y = nucleation[:, 1] 

2728 

2729 @property 

2730 def nucleation_time(self): 

2731 return self.nucleation_time__ 

2732 

2733 @nucleation_time.setter 

2734 def nucleation_time(self, nucleation_time): 

2735 if not isinstance(nucleation_time, num.ndarray) \ 

2736 and nucleation_time is not None: 

2737 nucleation_time = num.array([nucleation_time]) 

2738 

2739 self.nucleation_time__ = nucleation_time 

2740 

2741 @property 

2742 def patch_mask(self): 

2743 if (self.patch_mask__ is not None and 

2744 self.patch_mask__.shape == (self.nx * self.ny,)): 

2745 

2746 return self.patch_mask__ 

2747 else: 

2748 return num.ones(self.nx * self.ny, dtype=bool) 

2749 

2750 @patch_mask.setter 

2751 def patch_mask(self, patch_mask): 

2752 if isinstance(patch_mask, list): 

2753 patch_mask = num.array(patch_mask) 

2754 

2755 self.patch_mask__ = patch_mask 

2756 

2757 def get_tractions(self): 

2758 ''' 

2759 Get source traction vectors. 

2760 

2761 If :py:attr:`rake` is given, unit length directed traction vectors 

2762 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are returned, 

2763 else the given :py:attr:`tractions` are used. 

2764 

2765 :returns: 

2766 Traction vectors per patch. 

2767 :rtype: 

2768 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``. 

2769 ''' 

2770 

2771 if self.rake is not None: 

2772 if num.isnan(self.rake): 

2773 raise ValueError('Rake must be a real number, not NaN.') 

2774 

2775 logger.warning( 

2776 'Tractions are derived based on the given source rake.') 

2777 tractions = DirectedTractions(rake=self.rake) 

2778 else: 

2779 tractions = self.tractions 

2780 return tractions.get_tractions(self.nx, self.ny, self.patches) 

2781 

2782 def get_scaled_tractions(self, store): 

2783 ''' 

2784 Get traction vectors rescaled to given slip. 

2785 

2786 Opposing to :py:meth:`get_tractions` traction vectors 

2787 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are rescaled to 

2788 the given :py:attr:`slip` before returning. If no :py:attr:`slip` and 

2789 :py:attr:`rake` are provided, the given :py:attr:`tractions` are 

2790 returned without scaling. 

2791 

2792 :param store: 

2793 Green's function database (needs to cover whole region of the 

2794 source). 

2795 :type store: 

2796 :py:class:`~pyrocko.gf.store.Store` 

2797 

2798 :returns: 

2799 Traction vectors per patch. 

2800 :rtype: 

2801 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``. 

2802 ''' 

2803 tractions = self.tractions 

2804 factor = 1. 

2805 

2806 if self.rake is not None and self.slip is not None: 

2807 if num.isnan(self.rake): 

2808 raise ValueError('Rake must be a real number, not NaN.') 

2809 

2810 self.discretize_patches(store) 

2811 slip_0t = max(num.linalg.norm( 

2812 self.get_slip(scale_slip=False), 

2813 axis=1)) 

2814 

2815 factor = self.slip / slip_0t 

2816 tractions = DirectedTractions(rake=self.rake) 

2817 

2818 return tractions.get_tractions(self.nx, self.ny, self.patches) * factor 

2819 

2820 def base_key(self): 

2821 return SourceWithDerivedMagnitude.base_key(self) + ( 

2822 self.slip, 

2823 self.strike, 

2824 self.dip, 

2825 self.rake, 

2826 self.length, 

2827 self.width, 

2828 float(self.nucleation_x.mean()), 

2829 float(self.nucleation_y.mean()), 

2830 self.decimation_factor, 

2831 self.anchor, 

2832 self.pure_shear, 

2833 self.gamma, 

2834 tuple(self.patch_mask)) 

2835 

2836 def check_conflicts(self): 

2837 if self.tractions and self.rake: 

2838 raise AttributeError( 

2839 'Tractions and rake are mutually exclusive.') 

2840 if self.tractions is None and self.rake is None: 

2841 self.rake = 0. 

2842 

2843 def get_magnitude(self, store=None, target=None): 

2844 ''' 

2845 Get total seismic moment magnitude Mw. 

2846 

2847 :param store: 

2848 GF store to guide the discretization and providing the earthmodel 

2849 which is needed to calculate moment from slip. 

2850 :type store: 

2851 :py:class:`~pyrocko.gf.store.Store` 

2852 

2853 :param target: 

2854 Target, used to get GF interpolation settings. 

2855 :type target: 

2856 :py:class:`pyrocko.gf.targets.Target` 

2857 

2858 :returns: 

2859 Moment magnitude 

2860 :rtype: 

2861 float 

2862 ''' 

2863 self.check_conflicts() 

2864 if self.slip is not None or self.tractions is not None: 

2865 if store is None: 

2866 raise DerivedMagnitudeError( 

2867 'Magnitude for a rectangular source with slip or ' 

2868 'tractions defined can only be derived when earth model ' 

2869 'is set.') 

2870 

2871 moment_rate, calc_times = self.discretize_basesource( 

2872 store, target=target).get_moment_rate(store.config.deltat) 

2873 

2874 deltat = num.concatenate(( 

2875 (num.diff(calc_times)[0],), 

2876 num.diff(calc_times))) 

2877 

2878 return float(pmt.moment_to_magnitude( 

2879 num.sum(moment_rate * deltat))) 

2880 

2881 else: 

2882 return float(pmt.moment_to_magnitude(1.0)) 

2883 

2884 def get_factor(self): 

2885 return 1.0 

2886 

2887 def outline(self, cs='xyz'): 

2888 ''' 

2889 Get source outline corner coordinates. 

2890 

2891 :param cs: 

2892 :ref:`Output coordinate system <coordinate-system-names>`. 

2893 :type cs: 

2894 str 

2895 

2896 :returns: 

2897 Corner points in desired coordinate system. 

2898 :rtype: 

2899 :py:class:`~numpy.ndarray`: ``(5, [2, 3])``. 

2900 ''' 

2901 points = outline_rect_source(self.strike, self.dip, self.length, 

2902 self.width, self.anchor) 

2903 

2904 points[:, 0] += self.north_shift 

2905 points[:, 1] += self.east_shift 

2906 points[:, 2] += self.depth 

2907 if cs == 'xyz': 

2908 return points 

2909 elif cs == 'xy': 

2910 return points[:, :2] 

2911 elif cs in ('latlon', 'lonlat', 'latlondepth'): 

2912 latlon = ne_to_latlon( 

2913 self.lat, self.lon, points[:, 0], points[:, 1]) 

2914 

2915 latlon = num.array(latlon).T 

2916 if cs == 'latlon': 

2917 return latlon 

2918 elif cs == 'lonlat': 

2919 return latlon[:, ::-1] 

2920 else: 

2921 return num.concatenate( 

2922 (latlon, points[:, 2].reshape((len(points), 1))), 

2923 axis=1) 

2924 

2925 def points_on_source(self, cs='xyz', **kwargs): 

2926 ''' 

2927 Convert relative plane coordinates to geographical coordinates. 

2928 

2929 Given x and y coordinates (relative source coordinates between -1. 

2930 and 1.) are converted to desired geographical coordinates. Coordinates 

2931 need to be given as :py:class:`~numpy.ndarray` arguments ``points_x`` 

2932 and ``points_y``. 

2933 

2934 :param cs: 

2935 :ref:`Output coordinate system <coordinate-system-names>`. 

2936 :type cs: 

2937 str 

2938 

2939 :returns: 

2940 Point coordinates in desired coordinate system. 

2941 :rtype: 

2942 :py:class:`~numpy.ndarray`: ``(n_points, [2, 3])``. 

2943 ''' 

2944 points = points_on_rect_source( 

2945 self.strike, self.dip, self.length, self.width, 

2946 self.anchor, **kwargs) 

2947 

2948 points[:, 0] += self.north_shift 

2949 points[:, 1] += self.east_shift 

2950 points[:, 2] += self.depth 

2951 if cs == 'xyz': 

2952 return points 

2953 elif cs == 'xy': 

2954 return points[:, :2] 

2955 elif cs in ('latlon', 'lonlat', 'latlondepth'): 

2956 latlon = ne_to_latlon( 

2957 self.lat, self.lon, points[:, 0], points[:, 1]) 

2958 

2959 latlon = num.array(latlon).T 

2960 if cs == 'latlon': 

2961 return latlon 

2962 elif cs == 'lonlat': 

2963 return latlon[:, ::-1] 

2964 else: 

2965 return num.concatenate( 

2966 (latlon, points[:, 2].reshape((len(points), 1))), 

2967 axis=1) 

2968 

2969 def pyrocko_moment_tensor(self, store=None, target=None): 

2970 ''' 

2971 Get overall moment tensor of the rupture. 

2972 

2973 :param store: 

2974 GF store to guide the discretization and providing the earthmodel 

2975 which is needed to calculate moment from slip. 

2976 :type store: 

2977 :py:class:`~pyrocko.gf.store.Store` 

2978 

2979 :param target: 

2980 Target, used to get GF interpolation settings. 

2981 :type target: 

2982 :py:class:`pyrocko.gf.targets.Target` 

2983 

2984 :returns: 

2985 Moment tensor. 

2986 :rtype: 

2987 :py:class:`~pyrocko.moment_tensor.MomentTensor` 

2988 ''' 

2989 if store is not None: 

2990 if not self.patches: 

2991 self.discretize_patches(store) 

2992 

2993 data = self.get_slip() 

2994 else: 

2995 data = self.get_tractions() 

2996 

2997 weights = num.linalg.norm(data, axis=1) 

2998 weights /= weights.sum() 

2999 

3000 rakes = num.arctan2(data[:, 1], data[:, 0]) * r2d 

3001 rake = num.average(rakes, weights=weights) 

3002 

3003 return pmt.MomentTensor( 

3004 strike=self.strike, 

3005 dip=self.dip, 

3006 rake=rake, 

3007 scalar_moment=self.get_moment(store, target)) 

3008 

3009 def pyrocko_event(self, store=None, target=None, **kwargs): 

3010 return SourceWithDerivedMagnitude.pyrocko_event( 

3011 self, store, target, 

3012 **kwargs) 

3013 

3014 @classmethod 

3015 def from_pyrocko_event(cls, ev, **kwargs): 

3016 d = {} 

3017 mt = ev.moment_tensor 

3018 if mt: 

3019 (strike, dip, rake), _ = mt.both_strike_dip_rake() 

3020 d.update( 

3021 strike=float(strike), 

3022 dip=float(dip), 

3023 rake=float(rake)) 

3024 

3025 d.update(kwargs) 

3026 return super(PseudoDynamicRupture, cls).from_pyrocko_event(ev, **d) 

3027 

3028 def _discretize_points(self, store, *args, **kwargs): 

3029 ''' 

3030 Discretize source plane with equal vertical and horizontal spacing. 

3031 

3032 Additional ``*args`` and ``**kwargs`` are passed to 

3033 :py:meth:`points_on_source`. 

3034 

3035 :param store: 

3036 Green's function database (needs to cover whole region of the 

3037 source). 

3038 :type store: 

3039 :py:class:`~pyrocko.gf.store.Store` 

3040 

3041 :returns: 

3042 Number of points in strike and dip direction, distance 

3043 between adjacent points, coordinates (latlondepth) and coordinates 

3044 (xy on fault) for discrete points. 

3045 :rtype: 

3046 (int, int, float, :py:class:`~numpy.ndarray`, 

3047 :py:class:`~numpy.ndarray`). 

3048 ''' 

3049 anch_x, anch_y = map_anchor[self.anchor] 

3050 

3051 npoints = int(self.width // km) + 1 

3052 points = num.zeros((npoints, 3)) 

3053 points[:, 1] = num.linspace(-1., 1., npoints) 

3054 points[:, 1] = (points[:, 1] - anch_y) * self.width/2 

3055 

3056 rotmat = pmt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.0) 

3057 points = num.dot(rotmat.T, points.T).T 

3058 points[:, 2] += self.depth 

3059 

3060 vs_min = store.config.get_vs( 

3061 self.lat, self.lon, points, 

3062 interpolation='nearest_neighbor') 

3063 vr_min = max(vs_min.min(), .5*km) * self.gamma 

3064 

3065 oversampling = 10. 

3066 delta_l = self.length / (self.nx * oversampling) 

3067 delta_w = self.width / (self.ny * oversampling) 

3068 

3069 delta = self.eikonal_decimation * num.min([ 

3070 store.config.deltat * vr_min / oversampling, 

3071 delta_l, delta_w] + [ 

3072 deltas for deltas in store.config.deltas]) 

3073 

3074 delta = delta_w / num.ceil(delta_w / delta) 

3075 

3076 nx = int(num.ceil(self.length / delta)) + 1 

3077 ny = int(num.ceil(self.width / delta)) + 1 

3078 

3079 rem_l = (nx-1)*delta - self.length 

3080 lim_x = rem_l / self.length 

3081 

3082 points_xy = num.zeros((nx * ny, 2)) 

3083 points_xy[:, 0] = num.repeat( 

3084 num.linspace(-1.-lim_x, 1.+lim_x, nx), ny) 

3085 points_xy[:, 1] = num.tile( 

3086 num.linspace(-1., 1., ny), nx) 

3087 

3088 points = self.points_on_source( 

3089 points_x=points_xy[:, 0], 

3090 points_y=points_xy[:, 1], 

3091 **kwargs) 

3092 

3093 return nx, ny, delta, points, points_xy 

3094 

3095 def _discretize_rupture_v(self, store, interpolation='nearest_neighbor', 

3096 points=None): 

3097 ''' 

3098 Get rupture velocity for discrete points on source plane. 

3099 

3100 :param store: 

3101 Green's function database (needs to cover the whole region of the 

3102 source) 

3103 :type store: 

3104 :py:class:`~pyrocko.gf.store.Store` 

3105 

3106 :param interpolation: 

3107 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3108 and ``'multilinear'``). 

3109 :type interpolation: 

3110 str 

3111 

3112 :param points: 

3113 Coordinates on fault (-1.:1.) of discrete points. 

3114 :type points: 

3115 :py:class:`~numpy.ndarray`: ``(n_points, 2)`` 

3116 

3117 :returns: 

3118 Rupture velocity assumed as :math:`v_s * \\gamma` for discrete 

3119 points. 

3120 :rtype: 

3121 :py:class:`~numpy.ndarray`: ``(n_points, )``. 

3122 ''' 

3123 

3124 if points is None: 

3125 _, _, _, points, _ = self._discretize_points(store, cs='xyz') 

3126 

3127 return store.config.get_vs( 

3128 self.lat, self.lon, 

3129 points=points, 

3130 interpolation=interpolation) * self.gamma 

3131 

3132 def discretize_time( 

3133 self, store, interpolation='nearest_neighbor', 

3134 vr=None, times=None, *args, **kwargs): 

3135 ''' 

3136 Get rupture start time for discrete points on source plane. 

3137 

3138 :param store: 

3139 Green's function database (needs to cover whole region of the 

3140 source) 

3141 :type store: 

3142 :py:class:`~pyrocko.gf.store.Store` 

3143 

3144 :param interpolation: 

3145 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3146 and ``'multilinear'``). 

3147 :type interpolation: 

3148 str 

3149 

3150 :param vr: 

3151 Array, containing rupture user defined rupture velocity values. 

3152 :type vr: 

3153 :py:class:`~numpy.ndarray` 

3154 

3155 :param times: 

3156 Array, containing zeros, where rupture is starting, real positive 

3157 numbers at later secondary nucleation points and -1, where time 

3158 will be calculated. If not given, rupture starts at nucleation_x, 

3159 nucleation_y. Times are given for discrete points with equal 

3160 horizontal and vertical spacing. 

3161 :type times: 

3162 :py:class:`~numpy.ndarray` 

3163 

3164 :returns: 

3165 Coordinates (latlondepth), coordinates (xy), rupture velocity, 

3166 rupture propagation time of discrete points. 

3167 :rtype: 

3168 :py:class:`~numpy.ndarray`: ``(n_points, 3)``, 

3169 :py:class:`~numpy.ndarray`: ``(n_points, 2)``, 

3170 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``, 

3171 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``. 

3172 ''' 

3173 nx, ny, delta, points, points_xy = self._discretize_points( 

3174 store, cs='xyz') 

3175 

3176 if vr is None or vr.shape != tuple((nx, ny)): 

3177 if vr: 

3178 logger.warning( 

3179 'Given rupture velocities are not in right shape: ' 

3180 '(%i, %i), but needed is (%i, %i).', *vr.shape + (nx, ny)) 

3181 vr = self._discretize_rupture_v(store, interpolation, points)\ 

3182 .reshape(nx, ny) 

3183 

3184 if vr.shape != tuple((nx, ny)): 

3185 logger.warning( 

3186 'Given rupture velocities are not in right shape. Therefore' 

3187 ' standard rupture velocity array is used.') 

3188 

3189 def initialize_times(): 

3190 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y 

3191 

3192 if nucl_x.shape != nucl_y.shape: 

3193 raise ValueError( 

3194 'Nucleation coordinates have different shape.') 

3195 

3196 dist_points = num.array([ 

3197 num.linalg.norm(points_xy - num.array([x, y]).ravel(), axis=1) 

3198 for x, y in zip(nucl_x, nucl_y)]) 

3199 nucl_indices = num.argmin(dist_points, axis=1) 

3200 

3201 if self.nucleation_time is None: 

3202 nucl_times = num.zeros_like(nucl_indices) 

3203 else: 

3204 if self.nucleation_time.shape == nucl_x.shape: 

3205 nucl_times = self.nucleation_time 

3206 else: 

3207 raise ValueError( 

3208 'Nucleation coordinates and times have different ' 

3209 'shapes') 

3210 

3211 t = num.full(nx * ny, -1.) 

3212 t[nucl_indices] = nucl_times 

3213 return t.reshape(nx, ny) 

3214 

3215 if times is None: 

3216 times = initialize_times() 

3217 elif times.shape != tuple((nx, ny)): 

3218 times = initialize_times() 

3219 logger.warning( 

3220 'Given times are not in right shape. Therefore standard time' 

3221 ' array is used.') 

3222 

3223 eikonal_ext.eikonal_solver_fmm_cartesian( 

3224 speeds=vr, times=times, delta=delta) 

3225 

3226 return points, points_xy, vr, times 

3227 

3228 def get_vr_time_interpolators( 

3229 self, store, interpolation='nearest_neighbor', force=False, 

3230 **kwargs): 

3231 ''' 

3232 Get interpolators for rupture velocity and rupture time. 

3233 

3234 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`. 

3235 

3236 :param store: 

3237 Green's function database (needs to cover whole region of the 

3238 source). 

3239 :type store: 

3240 :py:class:`~pyrocko.gf.store.Store` 

3241 

3242 :param interpolation: 

3243 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3244 and ``'multilinear'``). 

3245 :type interpolation: 

3246 str 

3247 

3248 :param force: 

3249 Force recalculation of the interpolators (e.g. after change of 

3250 nucleation point locations/times). Default is ``False``. 

3251 :type force: 

3252 bool 

3253 ''' 

3254 interp_map = {'multilinear': 'linear', 'nearest_neighbor': 'nearest'} 

3255 if interpolation not in interp_map: 

3256 raise TypeError( 

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

3258 

3259 if not self._interpolators.get(interpolation, False) or force: 

3260 _, points_xy, vr, times = self.discretize_time( 

3261 store, **kwargs) 

3262 

3263 if self.length <= 0.: 

3264 raise ValueError( 

3265 'length must be larger then 0. not %g' % self.length) 

3266 

3267 if self.width <= 0.: 

3268 raise ValueError( 

3269 'width must be larger then 0. not %g' % self.width) 

3270 

3271 nx, ny = times.shape 

3272 anch_x, anch_y = map_anchor[self.anchor] 

3273 

3274 points_xy[:, 0] = (points_xy[:, 0] - anch_x) * self.length / 2. 

3275 points_xy[:, 1] = (points_xy[:, 1] - anch_y) * self.width / 2. 

3276 

3277 ascont = num.ascontiguousarray 

3278 

3279 self._interpolators[interpolation] = ( 

3280 nx, ny, times, vr, 

3281 RegularGridInterpolator( 

3282 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])), 

3283 times, 

3284 method=interp_map[interpolation]), 

3285 RegularGridInterpolator( 

3286 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])), 

3287 vr, 

3288 method=interp_map[interpolation])) 

3289 

3290 return self._interpolators[interpolation] 

3291 

3292 def discretize_patches( 

3293 self, store, interpolation='nearest_neighbor', force=False, 

3294 grid_shape=(), 

3295 **kwargs): 

3296 ''' 

3297 Get rupture start time and OkadaSource elements for points on rupture. 

3298 

3299 All source elements and their corresponding center points are 

3300 calculated and stored in the :py:attr:`patches` attribute. 

3301 

3302 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`. 

3303 

3304 :param store: 

3305 Green's function database (needs to cover whole region of the 

3306 source). 

3307 :type store: 

3308 :py:class:`~pyrocko.gf.store.Store` 

3309 

3310 :param interpolation: 

3311 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3312 and ``'multilinear'``). 

3313 :type interpolation: 

3314 str 

3315 

3316 :param force: 

3317 Force recalculation of the vr and time interpolators ( e.g. after 

3318 change of nucleation point locations/times). Default is ``False``. 

3319 :type force: 

3320 bool 

3321 

3322 :param grid_shape: 

3323 Desired sub fault patch grid size (nlength, nwidth). Either factor 

3324 or grid_shape should be set. 

3325 :type grid_shape: 

3326 :py:class:`tuple` of :py:class:`int` 

3327 ''' 

3328 nx, ny, times, vr, time_interpolator, vr_interpolator = \ 

3329 self.get_vr_time_interpolators( 

3330 store, 

3331 interpolation=interpolation, force=force, **kwargs) 

3332 anch_x, anch_y = map_anchor[self.anchor] 

3333 

3334 al = self.length / 2. 

3335 aw = self.width / 2. 

3336 al1 = -(al + anch_x * al) 

3337 al2 = al - anch_x * al 

3338 aw1 = -aw + anch_y * aw 

3339 aw2 = aw + anch_y * aw 

3340 assert num.abs([al1, al2]).sum() == self.length 

3341 assert num.abs([aw1, aw2]).sum() == self.width 

3342 

3343 def get_lame(*a, **kw): 

3344 shear_mod = store.config.get_shear_moduli(*a, **kw) 

3345 lamb = store.config.get_vp(*a, **kw)**2 \ 

3346 * store.config.get_rho(*a, **kw) - 2. * shear_mod 

3347 return shear_mod, lamb / (2. * (lamb + shear_mod)) 

3348 

3349 shear_mod, poisson = get_lame( 

3350 self.lat, self.lon, 

3351 num.array([[self.north_shift, self.east_shift, self.depth]]), 

3352 interpolation=interpolation) 

3353 

3354 okada_src = OkadaSource( 

3355 lat=self.lat, lon=self.lon, 

3356 strike=self.strike, dip=self.dip, 

3357 north_shift=self.north_shift, east_shift=self.east_shift, 

3358 depth=self.depth, 

3359 al1=al1, al2=al2, aw1=aw1, aw2=aw2, 

3360 poisson=poisson.mean(), 

3361 shearmod=shear_mod.mean(), 

3362 opening=kwargs.get('opening', 0.)) 

3363 

3364 if not (self.nx and self.ny): 

3365 if grid_shape: 

3366 self.nx, self.ny = grid_shape 

3367 else: 

3368 self.nx = nx 

3369 self.ny = ny 

3370 

3371 source_disc, source_points = okada_src.discretize(self.nx, self.ny) 

3372 

3373 shear_mod, poisson = get_lame( 

3374 self.lat, self.lon, 

3375 num.array([src.source_patch()[:3] for src in source_disc]), 

3376 interpolation=interpolation) 

3377 

3378 if (self.nx, self.ny) != (nx, ny): 

3379 times_interp = time_interpolator( 

3380 num.ascontiguousarray(source_points[:, :2])) 

3381 vr_interp = vr_interpolator( 

3382 num.ascontiguousarray(source_points[:, :2])) 

3383 else: 

3384 times_interp = times.T.ravel() 

3385 vr_interp = vr.T.ravel() 

3386 

3387 for isrc, src in enumerate(source_disc): 

3388 src.vr = vr_interp[isrc] 

3389 src.time = times_interp[isrc] + self.time 

3390 

3391 self.patches = source_disc 

3392 

3393 def discretize_basesource(self, store, target=None): 

3394 ''' 

3395 Prepare source for synthetic waveform calculation. 

3396 

3397 :param store: 

3398 Green's function database (needs to cover whole region of the 

3399 source). 

3400 :type store: 

3401 :py:class:`~pyrocko.gf.store.Store` 

3402 

3403 :param target: 

3404 Target information. 

3405 :type target: 

3406 :py:class:`~pyrocko.gf.targets.Target` 

3407 

3408 :returns: 

3409 Source discretized by a set of moment tensors and times. 

3410 :rtype: 

3411 :py:class:`~pyrocko.gf.meta.DiscretizedMTSource` 

3412 ''' 

3413 if not target: 

3414 interpolation = 'nearest_neighbor' 

3415 else: 

3416 interpolation = target.interpolation 

3417 

3418 if not self.patches: 

3419 self.discretize_patches(store, interpolation) 

3420 

3421 if self.coef_mat is None: 

3422 self.calc_coef_mat() 

3423 

3424 delta_slip, slip_times = self.get_delta_slip(store) 

3425 npatches = self.nx * self.ny 

3426 ntimes = slip_times.size 

3427 

3428 anch_x, anch_y = map_anchor[self.anchor] 

3429 

3430 pln = self.length / self.nx 

3431 pwd = self.width / self.ny 

3432 

3433 patch_coords = num.array([ 

3434 (p.ix, p.iy) 

3435 for p in self.patches]).reshape(self.nx, self.ny, 2) 

3436 

3437 # boundary condition is zero-slip 

3438 # is not valid to avoid unwished interpolation effects 

3439 slip_grid = num.zeros((self.nx + 2, self.ny + 2, ntimes, 3)) 

3440 slip_grid[1:-1, 1:-1, :, :] = \ 

3441 delta_slip.reshape(self.nx, self.ny, ntimes, 3) 

3442 

3443 slip_grid[0, 0, :, :] = slip_grid[1, 1, :, :] 

3444 slip_grid[0, -1, :, :] = slip_grid[1, -2, :, :] 

3445 slip_grid[-1, 0, :, :] = slip_grid[-2, 1, :, :] 

3446 slip_grid[-1, -1, :, :] = slip_grid[-2, -2, :, :] 

3447 

3448 slip_grid[1:-1, 0, :, :] = slip_grid[1:-1, 1, :, :] 

3449 slip_grid[1:-1, -1, :, :] = slip_grid[1:-1, -2, :, :] 

3450 slip_grid[0, 1:-1, :, :] = slip_grid[1, 1:-1, :, :] 

3451 slip_grid[-1, 1:-1, :, :] = slip_grid[-2, 1:-1, :, :] 

3452 

3453 def make_grid(patch_parameter): 

3454 grid = num.zeros((self.nx + 2, self.ny + 2)) 

3455 grid[1:-1, 1:-1] = patch_parameter.reshape(self.nx, self.ny) 

3456 

3457 grid[0, 0] = grid[1, 1] 

3458 grid[0, -1] = grid[1, -2] 

3459 grid[-1, 0] = grid[-2, 1] 

3460 grid[-1, -1] = grid[-2, -2] 

3461 

3462 grid[1:-1, 0] = grid[1:-1, 1] 

3463 grid[1:-1, -1] = grid[1:-1, -2] 

3464 grid[0, 1:-1] = grid[1, 1:-1] 

3465 grid[-1, 1:-1] = grid[-2, 1:-1] 

3466 

3467 return grid 

3468 

3469 lamb = self.get_patch_attribute('lamb') 

3470 mu = self.get_patch_attribute('shearmod') 

3471 

3472 lamb_grid = make_grid(lamb) 

3473 mu_grid = make_grid(mu) 

3474 

3475 coords_x = num.zeros(self.nx + 2) 

3476 coords_x[1:-1] = patch_coords[:, 0, 0] 

3477 coords_x[0] = coords_x[1] - pln / 2 

3478 coords_x[-1] = coords_x[-2] + pln / 2 

3479 

3480 coords_y = num.zeros(self.ny + 2) 

3481 coords_y[1:-1] = patch_coords[0, :, 1] 

3482 coords_y[0] = coords_y[1] - pwd / 2 

3483 coords_y[-1] = coords_y[-2] + pwd / 2 

3484 

3485 slip_interp = RegularGridInterpolator( 

3486 (coords_x, coords_y, slip_times), 

3487 slip_grid, method='nearest') 

3488 

3489 lamb_interp = RegularGridInterpolator( 

3490 (coords_x, coords_y), 

3491 lamb_grid, method='nearest') 

3492 

3493 mu_interp = RegularGridInterpolator( 

3494 (coords_x, coords_y), 

3495 mu_grid, method='nearest') 

3496 

3497 # discretize basesources 

3498 mindeltagf = min(tuple( 

3499 (self.length / self.nx, self.width / self.ny) + 

3500 tuple(store.config.deltas))) 

3501 

3502 nl = int((1. / self.decimation_factor) * 

3503 num.ceil(pln / mindeltagf)) + 1 

3504 nw = int((1. / self.decimation_factor) * 

3505 num.ceil(pwd / mindeltagf)) + 1 

3506 nsrc_patch = int(nl * nw) 

3507 dl = pln / nl 

3508 dw = pwd / nw 

3509 

3510 patch_area = dl * dw 

3511 

3512 xl = num.linspace(-0.5 * (pln - dl), 0.5 * (pln - dl), nl) 

3513 xw = num.linspace(-0.5 * (pwd - dw), 0.5 * (pwd - dw), nw) 

3514 

3515 base_coords = num.zeros((nsrc_patch, 3)) 

3516 base_coords[:, 0] = num.tile(xl, nw) 

3517 base_coords[:, 1] = num.repeat(xw, nl) 

3518 base_coords = num.tile(base_coords, (npatches, 1)) 

3519 

3520 center_coords = num.zeros((npatches, 3)) 

3521 center_coords[:, 0] = num.repeat( 

3522 num.arange(self.nx) * pln + pln / 2, self.ny) - self.length / 2 

3523 center_coords[:, 1] = num.tile( 

3524 num.arange(self.ny) * pwd + pwd / 2, self.nx) - self.width / 2 

3525 

3526 base_coords -= center_coords.repeat(nsrc_patch, axis=0) 

3527 nbaselocs = base_coords.shape[0] 

3528 

3529 base_interp = base_coords.repeat(ntimes, axis=0) 

3530 

3531 base_times = num.tile(slip_times, nbaselocs) 

3532 base_interp[:, 0] -= anch_x * self.length / 2 

3533 base_interp[:, 1] -= anch_y * self.width / 2 

3534 base_interp[:, 2] = base_times 

3535 

3536 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators( 

3537 store, interpolation=interpolation) 

3538 

3539 time_eikonal_max = time_interpolator.values.max() 

3540 

3541 nbasesrcs = base_interp.shape[0] 

3542 delta_slip = slip_interp(base_interp).reshape(nbaselocs, ntimes, 3) 

3543 lamb = lamb_interp(base_interp[:, :2]).ravel() 

3544 mu = mu_interp(base_interp[:, :2]).ravel() 

3545 

3546 if False: 

3547 try: 

3548 import matplotlib.pyplot as plt 

3549 coords = base_coords.copy() 

3550 norm = num.sum(num.linalg.norm(delta_slip, axis=2), axis=1) 

3551 plt.scatter(coords[:, 0], coords[:, 1], c=norm) 

3552 plt.show() 

3553 except AttributeError: 

3554 pass 

3555 

3556 base_interp[:, 2] = 0. 

3557 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0) 

3558 base_interp = num.dot(rotmat.T, base_interp.T).T 

3559 base_interp[:, 0] += self.north_shift 

3560 base_interp[:, 1] += self.east_shift 

3561 base_interp[:, 2] += self.depth 

3562 

3563 slip_strike = delta_slip[:, :, 0].ravel() 

3564 slip_dip = delta_slip[:, :, 1].ravel() 

3565 slip_norm = delta_slip[:, :, 2].ravel() 

3566 

3567 slip_shear = num.linalg.norm([slip_strike, slip_dip], axis=0) 

3568 slip_rake = r2d * num.arctan2(slip_dip, slip_strike) 

3569 

3570 m6s = okada_ext.patch2m6( 

3571 strikes=num.full(nbasesrcs, self.strike, dtype=float), 

3572 dips=num.full(nbasesrcs, self.dip, dtype=float), 

3573 rakes=slip_rake, 

3574 disl_shear=slip_shear, 

3575 disl_norm=slip_norm, 

3576 lamb=lamb, 

3577 mu=mu, 

3578 nthreads=self.nthreads) 

3579 

3580 m6s *= patch_area 

3581 

3582 dl = -self.patches[0].al1 + self.patches[0].al2 

3583 dw = -self.patches[0].aw1 + self.patches[0].aw2 

3584 

3585 base_times[base_times > time_eikonal_max] = time_eikonal_max 

3586 

3587 ds = meta.DiscretizedMTSource( 

3588 lat=self.lat, 

3589 lon=self.lon, 

3590 times=base_times + self.time, 

3591 north_shifts=base_interp[:, 0], 

3592 east_shifts=base_interp[:, 1], 

3593 depths=base_interp[:, 2], 

3594 m6s=m6s, 

3595 dl=dl, 

3596 dw=dw, 

3597 nl=self.nx, 

3598 nw=self.ny) 

3599 

3600 return ds 

3601 

3602 def calc_coef_mat(self): 

3603 ''' 

3604 Calculate coefficients connecting tractions and dislocations. 

3605 ''' 

3606 if not self.patches: 

3607 raise ValueError( 

3608 'Patches are needed. Please calculate them first.') 

3609 

3610 self.coef_mat = make_okada_coefficient_matrix( 

3611 self.patches, nthreads=self.nthreads, pure_shear=self.pure_shear) 

3612 

3613 def get_patch_attribute(self, attr): 

3614 ''' 

3615 Get patch attributes. 

3616 

3617 :param attr: 

3618 Name of selected attribute (see 

3619 :py:class`pyrocko.modelling.okada.OkadaSource`). 

3620 :type attr: 

3621 str 

3622 

3623 :returns: 

3624 Array with attribute value for each fault patch. 

3625 :rtype: 

3626 :py:class:`~numpy.ndarray` 

3627 

3628 ''' 

3629 if not self.patches: 

3630 raise ValueError( 

3631 'Patches are needed. Please calculate them first.') 

3632 return num.array([getattr(p, attr) for p in self.patches]) 

3633 

3634 def get_slip( 

3635 self, 

3636 time=None, 

3637 scale_slip=True, 

3638 interpolation='nearest_neighbor', 

3639 **kwargs): 

3640 ''' 

3641 Get slip per subfault patch for given time after rupture start. 

3642 

3643 :param time: 

3644 Time after origin [s], for which slip is computed. If not 

3645 given, final static slip is returned. 

3646 :type time: 

3647 float 

3648 

3649 :param scale_slip: 

3650 If ``True`` and :py:attr:`slip` given, all slip values are scaled 

3651 to fit the given maximum slip. 

3652 :type scale_slip: 

3653 bool 

3654 

3655 :param interpolation: 

3656 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3657 and ``'multilinear'``). 

3658 :type interpolation: 

3659 str 

3660 

3661 :returns: 

3662 Inverted dislocations (:math:`u_{strike}, u_{dip}, u_{tensile}`) 

3663 for each source patch. 

3664 :rtype: 

3665 :py:class:`~numpy.ndarray`: ``(n_sources, 3)`` 

3666 ''' 

3667 

3668 if self.patches is None: 

3669 raise ValueError( 

3670 'Please discretize the source first (discretize_patches())') 

3671 npatches = len(self.patches) 

3672 tractions = self.get_tractions() 

3673 time_patch_max = self.get_patch_attribute('time').max() - self.time 

3674 

3675 time_patch = time 

3676 if time is None: 

3677 time_patch = time_patch_max 

3678 

3679 if self.coef_mat is None: 

3680 self.calc_coef_mat() 

3681 

3682 if tractions.shape != (npatches, 3): 

3683 raise AttributeError( 

3684 'The traction vector is of invalid shape.' 

3685 ' Required shape is (npatches, 3)') 

3686 

3687 patch_mask = num.ones(npatches, dtype=bool) 

3688 if self.patch_mask is not None: 

3689 patch_mask = self.patch_mask 

3690 

3691 times = self.get_patch_attribute('time') - self.time 

3692 times[~patch_mask] = time_patch + 1. # exlcude unmasked patches 

3693 relevant_sources = num.nonzero(times <= time_patch)[0] 

3694 disloc_est = num.zeros_like(tractions) 

3695 

3696 if self.smooth_rupture: 

3697 patch_activation = num.zeros(npatches) 

3698 

3699 nx, ny, times, vr, time_interpolator, vr_interpolator = \ 

3700 self.get_vr_time_interpolators( 

3701 store, interpolation=interpolation) 

3702 

3703 # Getting the native Eikonal grid, bit hackish 

3704 points_x = num.round(time_interpolator.grid[0], decimals=2) 

3705 points_y = num.round(time_interpolator.grid[1], decimals=2) 

3706 times_eikonal = time_interpolator.values 

3707 

3708 time_max = time 

3709 if time is None: 

3710 time_max = times_eikonal.max() 

3711 

3712 for ip, p in enumerate(self.patches): 

3713 ul = num.round((p.ix + p.al1, p.iy + p.aw1), decimals=2) 

3714 lr = num.round((p.ix + p.al2, p.iy + p.aw2), decimals=2) 

3715 

3716 idx_length = num.logical_and( 

3717 points_x >= ul[0], points_x <= lr[0]) 

3718 idx_width = num.logical_and( 

3719 points_y >= ul[1], points_y <= lr[1]) 

3720 

3721 times_patch = times_eikonal[num.ix_(idx_length, idx_width)] 

3722 if times_patch.size == 0: 

3723 raise AttributeError('could not use smooth_rupture') 

3724 

3725 patch_activation[ip] = \ 

3726 (times_patch <= time_max).sum() / times_patch.size 

3727 

3728 if time_patch == 0 and time_patch != time_patch_max: 

3729 patch_activation[ip] = 0. 

3730 

3731 patch_activation[~patch_mask] = 0. # exlcude unmasked patches 

3732 

3733 relevant_sources = num.nonzero(patch_activation > 0.)[0] 

3734 

3735 if relevant_sources.size == 0: 

3736 return disloc_est 

3737 

3738 indices_disl = num.repeat(relevant_sources * 3, 3) 

3739 indices_disl[1::3] += 1 

3740 indices_disl[2::3] += 2 

3741 

3742 disloc_est[relevant_sources] = invert_fault_dislocations_bem( 

3743 stress_field=tractions[relevant_sources, :].ravel(), 

3744 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl], 

3745 pure_shear=self.pure_shear, nthreads=self.nthreads, 

3746 epsilon=None, 

3747 **kwargs) 

3748 

3749 if self.smooth_rupture: 

3750 disloc_est *= patch_activation[:, num.newaxis] 

3751 

3752 if scale_slip and self.slip is not None: 

3753 disloc_tmax = num.zeros(npatches) 

3754 

3755 indices_disl = num.repeat(num.nonzero(patch_mask)[0] * 3, 3) 

3756 indices_disl[1::3] += 1 

3757 indices_disl[2::3] += 2 

3758 

3759 disloc_tmax[patch_mask] = num.linalg.norm( 

3760 invert_fault_dislocations_bem( 

3761 stress_field=tractions[patch_mask, :].ravel(), 

3762 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl], 

3763 pure_shear=self.pure_shear, nthreads=self.nthreads, 

3764 epsilon=None, 

3765 **kwargs), axis=1) 

3766 

3767 disloc_tmax_max = disloc_tmax.max() 

3768 if disloc_tmax_max == 0.: 

3769 logger.warning( 

3770 'slip scaling not performed. Maximum slip is 0.') 

3771 

3772 disloc_est *= self.slip / disloc_tmax_max 

3773 

3774 return disloc_est 

3775 

3776 def get_delta_slip( 

3777 self, 

3778 store=None, 

3779 deltat=None, 

3780 delta=True, 

3781 interpolation='nearest_neighbor', 

3782 **kwargs): 

3783 ''' 

3784 Get slip change snapshots. 

3785 

3786 The time interval, within which the slip changes are computed is 

3787 determined by the sampling rate of the Green's function database or 

3788 ``deltat``. Additional ``**kwargs`` are passed to :py:meth:`get_slip`. 

3789 

3790 :param store: 

3791 Green's function database (needs to cover whole region of of the 

3792 source). Its sampling interval is used as time increment for slip 

3793 difference calculation. Either ``deltat`` or ``store`` should be 

3794 given. 

3795 :type store: 

3796 :py:class:`~pyrocko.gf.store.Store` 

3797 

3798 :param deltat: 

3799 Time interval for slip difference calculation [s]. Either 

3800 ``deltat`` or ``store`` should be given. 

3801 :type deltat: 

3802 float 

3803 

3804 :param delta: 

3805 If ``True``, slip differences between two time steps are given. If 

3806 ``False``, cumulative slip for all time steps. 

3807 :type delta: 

3808 bool 

3809 

3810 :param interpolation: 

3811 Interpolation method to use (choose between ``'nearest_neighbor'`` 

3812 and ``'multilinear'``). 

3813 :type interpolation: 

3814 str 

3815 

3816 :returns: 

3817 Displacement changes(:math:`\\Delta u_{strike}, 

3818 \\Delta u_{dip} , \\Delta u_{tensile}`) for each source patch and 

3819 time; corner times, for which delta slip is computed. The order of 

3820 displacement changes array is: 

3821 

3822 .. math:: 

3823 

3824 &[[\\\\ 

3825 &[\\Delta u_{strike, patch1, t1}, 

3826 \\Delta u_{dip, patch1, t1}, 

3827 \\Delta u_{tensile, patch1, t1}],\\\\ 

3828 &[\\Delta u_{strike, patch1, t2}, 

3829 \\Delta u_{dip, patch1, t2}, 

3830 \\Delta u_{tensile, patch1, t2}]\\\\ 

3831 &], [\\\\ 

3832 &[\\Delta u_{strike, patch2, t1}, ...],\\\\ 

3833 &[\\Delta u_{strike, patch2, t2}, ...]]]\\\\ 

3834 

3835 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``, 

3836 :py:class:`~numpy.ndarray`: ``(n_times, )`` 

3837 ''' 

3838 if store and deltat: 

3839 raise AttributeError( 

3840 'Argument collision. ' 

3841 'Please define only the store or the deltat argument.') 

3842 

3843 if store: 

3844 deltat = store.config.deltat 

3845 

3846 if not deltat: 

3847 raise AttributeError('Please give a GF store or set deltat.') 

3848 

3849 npatches = len(self.patches) 

3850 

3851 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators( 

3852 store, interpolation=interpolation) 

3853 tmax = time_interpolator.values.max() 

3854 

3855 calc_times = num.arange(0., tmax + deltat, deltat) 

3856 calc_times[calc_times > tmax] = tmax 

3857 

3858 disloc_est = num.zeros((npatches, calc_times.size, 3)) 

3859 

3860 for itime, t in enumerate(calc_times): 

3861 disloc_est[:, itime, :] = self.get_slip( 

3862 time=t, scale_slip=False, **kwargs) 

3863 

3864 if self.slip: 

3865 disloc_tmax = num.linalg.norm( 

3866 self.get_slip(scale_slip=False, time=tmax), 

3867 axis=1) 

3868 

3869 disloc_tmax_max = disloc_tmax.max() 

3870 if disloc_tmax_max == 0.: 

3871 logger.warning( 

3872 'Slip scaling not performed. Maximum slip is 0.') 

3873 else: 

3874 disloc_est *= self.slip / disloc_tmax_max 

3875 

3876 if not delta: 

3877 return disloc_est, calc_times 

3878 

3879 # if we have only one timestep there is no gradient 

3880 if calc_times.size > 1: 

3881 disloc_init = disloc_est[:, 0, :] 

3882 disloc_est = num.diff(disloc_est, axis=1) 

3883 disloc_est = num.concatenate(( 

3884 disloc_init[:, num.newaxis, :], disloc_est), axis=1) 

3885 

3886 calc_times = calc_times 

3887 

3888 return disloc_est, calc_times 

3889 

3890 def get_slip_rate(self, *args, **kwargs): 

3891 ''' 

3892 Get slip rate inverted from patches. 

3893 

3894 The time interval, within which the slip rates are computed is 

3895 determined by the sampling rate of the Green's function database or 

3896 ``deltat``. Additional ``*args`` and ``**kwargs`` are passed to 

3897 :py:meth:`get_delta_slip`. 

3898 

3899 :returns: 

3900 Slip rates (:math:`\\Delta u_{strike}/\\Delta t`, 

3901 :math:`\\Delta u_{dip}/\\Delta t, \\Delta u_{tensile}/\\Delta t`) 

3902 for each source patch and time; corner times, for which slip rate 

3903 is computed. The order of sliprate array is: 

3904 

3905 .. math:: 

3906 

3907 &[[\\\\ 

3908 &[\\Delta u_{strike, patch1, t1}/\\Delta t, 

3909 \\Delta u_{dip, patch1, t1}/\\Delta t, 

3910 \\Delta u_{tensile, patch1, t1}/\\Delta t],\\\\ 

3911 &[\\Delta u_{strike, patch1, t2}/\\Delta t, 

3912 \\Delta u_{dip, patch1, t2}/\\Delta t, 

3913 \\Delta u_{tensile, patch1, t2}/\\Delta t]], [\\\\ 

3914 &[\\Delta u_{strike, patch2, t1}/\\Delta t, ...],\\\\ 

3915 &[\\Delta u_{strike, patch2, t2}/\\Delta t, ...]]]\\\\ 

3916 

3917 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``, 

3918 :py:class:`~numpy.ndarray`: ``(n_times, )`` 

3919 ''' 

3920 ddisloc_est, calc_times = self.get_delta_slip( 

3921 *args, delta=True, **kwargs) 

3922 

3923 dt = num.concatenate( 

3924 [(num.diff(calc_times)[0], ), num.diff(calc_times)]) 

3925 slip_rate = num.linalg.norm(ddisloc_est, axis=2) / dt 

3926 

3927 return slip_rate, calc_times 

3928 

3929 def get_moment_rate_patches(self, *args, **kwargs): 

3930 ''' 

3931 Get scalar seismic moment rate for each patch individually. 

3932 

3933 Additional ``*args`` and ``**kwargs`` are passed to 

3934 :py:meth:`get_slip_rate`. 

3935 

3936 :returns: 

3937 Seismic moment rate for each source patch and time; corner times, 

3938 for which patch moment rate is computed based on slip rate. The 

3939 order of the moment rate array is: 

3940 

3941 .. math:: 

3942 

3943 &[\\\\ 

3944 &[(\\Delta M / \\Delta t)_{patch1, t1}, 

3945 (\\Delta M / \\Delta t)_{patch1, t2}, ...],\\\\ 

3946 &[(\\Delta M / \\Delta t)_{patch2, t1}, 

3947 (\\Delta M / \\Delta t)_{patch, t2}, ...],\\\\ 

3948 &[...]]\\\\ 

3949 

3950 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times)``, 

3951 :py:class:`~numpy.ndarray`: ``(n_times, )`` 

3952 ''' 

3953 slip_rate, calc_times = self.get_slip_rate(*args, **kwargs) 

3954 

3955 shear_mod = self.get_patch_attribute('shearmod') 

3956 p_length = self.get_patch_attribute('length') 

3957 p_width = self.get_patch_attribute('width') 

3958 

3959 dA = p_length * p_width 

3960 

3961 mom_rate = shear_mod[:, num.newaxis] * slip_rate * dA[:, num.newaxis] 

3962 

3963 return mom_rate, calc_times 

3964 

3965 def get_moment_rate(self, store, target=None, deltat=None): 

3966 ''' 

3967 Get seismic source moment rate for the total source (STF). 

3968 

3969 :param store: 

3970 Green's function database (needs to cover whole region of of the 

3971 source). Its ``deltat`` [s] is used as time increment for slip 

3972 difference calculation. Either ``deltat`` or ``store`` should be 

3973 given. 

3974 :type store: 

3975 :py:class:`~pyrocko.gf.store.Store` 

3976 

3977 :param target: 

3978 Target information, needed for interpolation method. 

3979 :type target: 

3980 :py:class:`~pyrocko.gf.targets.Target` 

3981 

3982 :param deltat: 

3983 Time increment for slip difference calculation [s]. If not given 

3984 ``store.deltat`` is used. 

3985 :type deltat: 

3986 float 

3987 

3988 :return: 

3989 Seismic moment rate [Nm/s] for each time; corner times, for which 

3990 moment rate is computed. The order of the moment rate array is: 

3991 

3992 .. math:: 

3993 

3994 &[\\\\ 

3995 &(\\Delta M / \\Delta t)_{t1},\\\\ 

3996 &(\\Delta M / \\Delta t)_{t2},\\\\ 

3997 &...]\\\\ 

3998 

3999 :rtype: 

4000 :py:class:`~numpy.ndarray`: ``(n_times, )``, 

4001 :py:class:`~numpy.ndarray`: ``(n_times, )`` 

4002 ''' 

4003 if not deltat: 

4004 deltat = store.config.deltat 

4005 return self.discretize_basesource( 

4006 store, target=target).get_moment_rate(deltat) 

4007 

4008 def get_moment(self, *args, **kwargs): 

4009 ''' 

4010 Get cumulative seismic moment. 

4011 

4012 Additional ``*args`` and ``**kwargs`` are passed to 

4013 :py:meth:`get_magnitude`. 

4014 

4015 :returns: 

4016 Cumulative seismic moment in [Nm]. 

4017 :rtype: 

4018 float 

4019 ''' 

4020 return float(pmt.magnitude_to_moment(self.get_magnitude( 

4021 *args, **kwargs))) 

4022 

4023 def rescale_slip(self, magnitude=None, moment=None, **kwargs): 

4024 ''' 

4025 Rescale source slip based on given target magnitude or seismic moment. 

4026 

4027 Rescale the maximum source slip to fit the source moment magnitude or 

4028 seismic moment to the given target values. Either ``magnitude`` or 

4029 ``moment`` need to be given. Additional ``**kwargs`` are passed to 

4030 :py:meth:`get_moment`. 

4031 

4032 :param magnitude: 

4033 Target moment magnitude :math:`M_\\mathrm{w}` as in 

4034 [Hanks and Kanamori, 1979] 

4035 :type magnitude: 

4036 float 

4037 

4038 :param moment: 

4039 Target seismic moment :math:`M_0` [Nm]. 

4040 :type moment: 

4041 float 

4042 ''' 

4043 if self.slip is None: 

4044 self.slip = 1. 

4045 logger.warning('No slip found for rescaling. ' 

4046 'An initial slip of 1 m is assumed.') 

4047 

4048 if magnitude is None and moment is None: 

4049 raise ValueError( 

4050 'Either target magnitude or moment need to be given.') 

4051 

4052 moment_init = self.get_moment(**kwargs) 

4053 

4054 if magnitude is not None: 

4055 moment = pmt.magnitude_to_moment(magnitude) 

4056 

4057 self.slip *= moment / moment_init 

4058 

4059 def get_centroid(self, store, *args, **kwargs): 

4060 ''' 

4061 Centroid of the pseudo dynamic rupture model. 

4062 

4063 The centroid location and time are derived from the locations and times 

4064 of the individual patches weighted with their moment contribution. 

4065 Additional ``**kwargs`` are passed to :py:meth:`pyrocko_moment_tensor`. 

4066 

4067 :param store: 

4068 Green's function database (needs to cover whole region of of the 

4069 source). Its ``deltat`` [s] is used as time increment for slip 

4070 difference calculation. Either ``deltat`` or ``store`` should be 

4071 given. 

4072 :type store: 

4073 :py:class:`~pyrocko.gf.store.Store` 

4074 

4075 :returns: 

4076 The centroid location and associated moment tensor. 

4077 :rtype: 

4078 :py:class:`pyrocko.model.event.Event` 

4079 ''' 

4080 _, _, _, _, time, _ = self.get_vr_time_interpolators(store) 

4081 t_max = time.values.max() 

4082 

4083 moment_rate, times = self.get_moment_rate_patches(deltat=t_max) 

4084 

4085 moment = num.sum(moment_rate * times, axis=1) 

4086 weights = moment / moment.sum() 

4087 

4088 norths = self.get_patch_attribute('north_shift') 

4089 easts = self.get_patch_attribute('east_shift') 

4090 depths = self.get_patch_attribute('depth') 

4091 

4092 centroid_n = num.sum(weights * norths) 

4093 centroid_e = num.sum(weights * easts) 

4094 centroid_d = num.sum(weights * depths) 

4095 

4096 centroid_lat, centroid_lon = ne_to_latlon( 

4097 self.lat, self.lon, centroid_n, centroid_e) 

4098 

4099 moment_rate_, times = self.get_moment_rate(store) 

4100 delta_times = num.concatenate(( 

4101 [times[1] - times[0]], 

4102 num.diff(times))) 

4103 moment_src = delta_times * moment_rate 

4104 

4105 centroid_t = num.sum( 

4106 moment_src / num.sum(moment_src) * times) + self.time 

4107 

4108 mt = self.pyrocko_moment_tensor(store, *args, **kwargs) 

4109 

4110 return model.Event( 

4111 lat=centroid_lat, 

4112 lon=centroid_lon, 

4113 depth=centroid_d, 

4114 time=centroid_t, 

4115 moment_tensor=mt, 

4116 magnitude=mt.magnitude, 

4117 duration=t_max) 

4118 

4119 def get_coulomb_failure_stress( 

4120 self, 

4121 receiver_points, 

4122 friction, 

4123 pressure, 

4124 strike, 

4125 dip, 

4126 rake, 

4127 time=None, 

4128 *args, 

4129 **kwargs): 

4130 ''' 

4131 Calculate Coulomb failure stress change CFS. 

4132 

4133 The function obtains the Coulomb failure stress change :math:`\\Delta 

4134 \\sigma_C` at arbitrary receiver points with a commonly oriented 

4135 receiver plane assuming: 

4136 

4137 .. math:: 

4138 

4139 \\Delta \\sigma_C = \\sigma_S - \\mu (\\sigma_N - \\Delta p) 

4140 

4141 with the shear stress :math:`\\sigma_S`, the coefficient of friction 

4142 :math:`\\mu`, the normal stress :math:`\\sigma_N`, and the pore fluid 

4143 pressure change :math:`\\Delta p`. Each receiver point is characterized 

4144 by its geographical coordinates, and depth. The required receiver plane 

4145 orientation is defined by ``strike``, ``dip``, and ``rake``. The 

4146 Coulomb failure stress change is calculated for a given time after 

4147 rupture origin time. 

4148 

4149 :param receiver_points: 

4150 Location of the receiver points in Northing, Easting, and depth in 

4151 [m]. 

4152 :type receiver_points: 

4153 :py:class:`~numpy.ndarray`: ``(n_receiver, 3)`` 

4154 

4155 :param friction: 

4156 Coefficient of friction. 

4157 :type friction: 

4158 float 

4159 

4160 :param pressure: 

4161 Pore pressure change in [Pa]. 

4162 :type pressure: 

4163 float 

4164 

4165 :param strike: 

4166 Strike of the receiver plane in [deg]. 

4167 :type strike: 

4168 float 

4169 

4170 :param dip: 

4171 Dip of the receiver plane in [deg]. 

4172 :type dip: 

4173 float 

4174 

4175 :param rake: 

4176 Rake of the receiver plane in [deg]. 

4177 :type rake: 

4178 float 

4179 

4180 :param time: 

4181 Time after origin [s], for which the resulting :math:`\\Delta 

4182 \\Sigma_c` is computed. If not given, :math:`\\Delta \\Sigma_c` is 

4183 derived based on the final static slip. 

4184 :type time: 

4185 float 

4186 

4187 :returns: 

4188 The Coulomb failure stress change :math:`\\Delta \\Sigma_c` at each 

4189 receiver point in [Pa]. 

4190 :rtype: 

4191 :py:class:`~numpy.ndarray`: ``(n_receiver,)`` 

4192 ''' 

4193 # dislocation at given time 

4194 source_slip = self.get_slip(time=time, scale_slip=True) 

4195 

4196 # source planes 

4197 source_patches = num.array([ 

4198 src.source_patch() for src in self.patches]) 

4199 

4200 # earth model 

4201 lambda_mean = num.mean([src.lamb for src in self.patches]) 

4202 mu_mean = num.mean([src.shearmod for src in self.patches]) 

4203 

4204 # Dislocation and spatial derivatives from okada in NED 

4205 results = okada_ext.okada( 

4206 source_patches, 

4207 source_slip, 

4208 receiver_points, 

4209 lambda_mean, 

4210 mu_mean, 

4211 rotate_sdn=False, # TODO Check 

4212 stack_sources=0, # TODO Check 

4213 *args, **kwargs) 

4214 

4215 # resolve stress tensor (sum!) 

4216 diag_ind = [0, 4, 8] 

4217 kron = num.zeros(9) 

4218 kron[diag_ind] = 1. 

4219 kron = kron[num.newaxis, num.newaxis, :] 

4220 

4221 eps = 0.5 * ( 

4222 results[:, :, 3:] + 

4223 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)]) 

4224 

4225 dilatation \ 

4226 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis] 

4227 

4228 stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps 

4229 

4230 # superposed stress of all sources at receiver locations 

4231 stress_sum = num.sum(stress, axis=0) 

4232 

4233 # get shear and normal stress from stress tensor 

4234 strike_rad = d2r * strike 

4235 dip_rad = d2r * dip 

4236 rake_rad = d2r * rake 

4237 

4238 n_rec = receiver_points.shape[0] 

4239 stress_normal = num.zeros(n_rec) 

4240 tau = num.zeros(n_rec) 

4241 

4242 # Get vectors in receiver fault normal (ns), strike (rst) and 

4243 # dip (rdi) direction 

4244 ns = num.zeros(3) 

4245 rst = num.zeros(3) 

4246 rdi = num.zeros(3) 

4247 

4248 ns[0] = num.sin(dip_rad) * num.cos(strike_rad + 0.5 * num.pi) 

4249 ns[1] = num.sin(dip_rad) * num.sin(strike_rad + 0.5 * num.pi) 

4250 ns[2] = -num.cos(dip_rad) 

4251 

4252 rst[0] = num.cos(strike_rad) 

4253 rst[1] = num.sin(strike_rad) 

4254 rst[2] = 0.0 

4255 

4256 rdi[0] = num.cos(dip_rad) * num.cos(strike_rad + 0.5 * num.pi) 

4257 rdi[1] = num.cos(dip_rad) * num.sin(strike_rad + 0.5 * num.pi) 

4258 rdi[2] = num.sin(dip_rad) 

4259 

4260 ts = rst * num.cos(rake_rad) - rdi * num.sin(rake_rad) 

4261 

4262 stress_normal = num.sum( 

4263 num.tile(ns, 3) * stress_sum * num.repeat(ns, 3), axis=1) 

4264 

4265 tau = num.sum( 

4266 num.tile(ts, 3) * stress_sum * num.repeat(ns, 3), axis=1) 

4267 

4268 # calculate cfs using formula above and return 

4269 return tau + friction * (stress_normal + pressure) 

4270 

4271 

4272class DoubleDCSource(SourceWithMagnitude): 

4273 ''' 

4274 Two double-couple point sources separated in space and time. 

4275 Moment share between the sub-sources is controlled by the 

4276 parameter mix. 

4277 The position of the subsources is dependent on the moment 

4278 distribution between the two sources. Depth, east and north 

4279 shift are given for the centroid between the two double-couples. 

4280 The subsources will positioned according to their moment shares 

4281 around this centroid position. 

4282 This is done according to their delta parameters, which are 

4283 therefore in relation to that centroid. 

4284 Note that depth of the subsources therefore can be 

4285 depth+/-delta_depth. For shallow earthquakes therefore 

4286 the depth has to be chosen deeper to avoid sampling 

4287 above surface. 

4288 ''' 

4289 

4290 strike1 = Float.T( 

4291 default=0.0, 

4292 help='strike direction in [deg], measured clockwise from north') 

4293 

4294 dip1 = Float.T( 

4295 default=90.0, 

4296 help='dip angle in [deg], measured downward from horizontal') 

4297 

4298 azimuth = Float.T( 

4299 default=0.0, 

4300 help='azimuth to second double-couple [deg], ' 

4301 'measured at first, clockwise from north') 

4302 

4303 rake1 = Float.T( 

4304 default=0.0, 

4305 help='rake angle in [deg], ' 

4306 'measured counter-clockwise from right-horizontal ' 

4307 'in on-plane view') 

4308 

4309 strike2 = Float.T( 

4310 default=0.0, 

4311 help='strike direction in [deg], measured clockwise from north') 

4312 

4313 dip2 = Float.T( 

4314 default=90.0, 

4315 help='dip angle in [deg], measured downward from horizontal') 

4316 

4317 rake2 = Float.T( 

4318 default=0.0, 

4319 help='rake angle in [deg], ' 

4320 'measured counter-clockwise from right-horizontal ' 

4321 'in on-plane view') 

4322 

4323 delta_time = Float.T( 

4324 default=0.0, 

4325 help='separation of double-couples in time (t2-t1) [s]') 

4326 

4327 delta_depth = Float.T( 

4328 default=0.0, 

4329 help='difference in depth (z2-z1) [m]') 

4330 

4331 distance = Float.T( 

4332 default=0.0, 

4333 help='distance between the two double-couples [m]') 

4334 

4335 mix = Float.T( 

4336 default=0.5, 

4337 help='how to distribute the moment to the two doublecouples ' 

4338 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1') 

4339 

4340 stf1 = STF.T( 

4341 optional=True, 

4342 help='Source time function of subsource 1 ' 

4343 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)') 

4344 

4345 stf2 = STF.T( 

4346 optional=True, 

4347 help='Source time function of subsource 2 ' 

4348 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)') 

4349 

4350 discretized_source_class = meta.DiscretizedMTSource 

4351 

4352 def base_key(self): 

4353 return ( 

4354 self.time, self.depth, self.lat, self.north_shift, 

4355 self.lon, self.east_shift, type(self).__name__) + \ 

4356 self.effective_stf1_pre().base_key() + \ 

4357 self.effective_stf2_pre().base_key() + ( 

4358 self.strike1, self.dip1, self.rake1, 

4359 self.strike2, self.dip2, self.rake2, 

4360 self.delta_time, self.delta_depth, 

4361 self.azimuth, self.distance, self.mix) 

4362 

4363 def get_factor(self): 

4364 return self.moment 

4365 

4366 def effective_stf1_pre(self): 

4367 return self.stf1 or self.stf or g_unit_pulse 

4368 

4369 def effective_stf2_pre(self): 

4370 return self.stf2 or self.stf or g_unit_pulse 

4371 

4372 def effective_stf_post(self): 

4373 return g_unit_pulse 

4374 

4375 def split(self): 

4376 a1 = 1.0 - self.mix 

4377 a2 = self.mix 

4378 delta_north = math.cos(self.azimuth * d2r) * self.distance 

4379 delta_east = math.sin(self.azimuth * d2r) * self.distance 

4380 

4381 dc1 = DCSource( 

4382 lat=self.lat, 

4383 lon=self.lon, 

4384 time=self.time - self.delta_time * a2, 

4385 north_shift=self.north_shift - delta_north * a2, 

4386 east_shift=self.east_shift - delta_east * a2, 

4387 depth=self.depth - self.delta_depth * a2, 

4388 moment=self.moment * a1, 

4389 strike=self.strike1, 

4390 dip=self.dip1, 

4391 rake=self.rake1, 

4392 stf=self.stf1 or self.stf) 

4393 

4394 dc2 = DCSource( 

4395 lat=self.lat, 

4396 lon=self.lon, 

4397 time=self.time + self.delta_time * a1, 

4398 north_shift=self.north_shift + delta_north * a1, 

4399 east_shift=self.east_shift + delta_east * a1, 

4400 depth=self.depth + self.delta_depth * a1, 

4401 moment=self.moment * a2, 

4402 strike=self.strike2, 

4403 dip=self.dip2, 

4404 rake=self.rake2, 

4405 stf=self.stf2 or self.stf) 

4406 

4407 return [dc1, dc2] 

4408 

4409 def discretize_basesource(self, store, target=None): 

4410 a1 = 1.0 - self.mix 

4411 a2 = self.mix 

4412 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1, 

4413 rake=self.rake1, scalar_moment=a1) 

4414 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2, 

4415 rake=self.rake2, scalar_moment=a2) 

4416 

4417 delta_north = math.cos(self.azimuth * d2r) * self.distance 

4418 delta_east = math.sin(self.azimuth * d2r) * self.distance 

4419 

4420 times1, amplitudes1 = self.effective_stf1_pre().discretize_t( 

4421 store.config.deltat, self.time - self.delta_time * a2) 

4422 

4423 times2, amplitudes2 = self.effective_stf2_pre().discretize_t( 

4424 store.config.deltat, self.time + self.delta_time * a1) 

4425 

4426 nt1 = times1.size 

4427 nt2 = times2.size 

4428 

4429 ds = meta.DiscretizedMTSource( 

4430 lat=self.lat, 

4431 lon=self.lon, 

4432 times=num.concatenate((times1, times2)), 

4433 north_shifts=num.concatenate(( 

4434 num.repeat(self.north_shift - delta_north * a2, nt1), 

4435 num.repeat(self.north_shift + delta_north * a1, nt2))), 

4436 east_shifts=num.concatenate(( 

4437 num.repeat(self.east_shift - delta_east * a2, nt1), 

4438 num.repeat(self.east_shift + delta_east * a1, nt2))), 

4439 depths=num.concatenate(( 

4440 num.repeat(self.depth - self.delta_depth * a2, nt1), 

4441 num.repeat(self.depth + self.delta_depth * a1, nt2))), 

4442 m6s=num.vstack(( 

4443 mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis], 

4444 mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis]))) 

4445 

4446 return ds 

4447 

4448 def pyrocko_moment_tensor(self, store=None, target=None): 

4449 a1 = 1.0 - self.mix 

4450 a2 = self.mix 

4451 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1, 

4452 rake=self.rake1, 

4453 scalar_moment=a1 * self.moment) 

4454 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2, 

4455 rake=self.rake2, 

4456 scalar_moment=a2 * self.moment) 

4457 return pmt.MomentTensor(m=mot1.m() + mot2.m()) 

4458 

4459 def pyrocko_event(self, store=None, target=None, **kwargs): 

4460 return SourceWithMagnitude.pyrocko_event( 

4461 self, store, target, 

4462 moment_tensor=self.pyrocko_moment_tensor(store, target), 

4463 **kwargs) 

4464 

4465 @classmethod 

4466 def from_pyrocko_event(cls, ev, **kwargs): 

4467 d = {} 

4468 mt = ev.moment_tensor 

4469 if mt: 

4470 (strike, dip, rake), _ = mt.both_strike_dip_rake() 

4471 d.update( 

4472 strike1=float(strike), 

4473 dip1=float(dip), 

4474 rake1=float(rake), 

4475 strike2=float(strike), 

4476 dip2=float(dip), 

4477 rake2=float(rake), 

4478 mix=0.0, 

4479 magnitude=float(mt.moment_magnitude())) 

4480 

4481 d.update(kwargs) 

4482 source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d) 

4483 source.stf1 = source.stf 

4484 source.stf2 = HalfSinusoidSTF(effective_duration=0.) 

4485 source.stf = None 

4486 return source 

4487 

4488 

4489class RingfaultSource(SourceWithMagnitude): 

4490 ''' 

4491 A ring fault with vertical doublecouples. 

4492 ''' 

4493 

4494 diameter = Float.T( 

4495 default=1.0, 

4496 help='diameter of the ring in [m]') 

4497 

4498 sign = Float.T( 

4499 default=1.0, 

4500 help='inside of the ring moves up (+1) or down (-1)') 

4501 

4502 strike = Float.T( 

4503 default=0.0, 

4504 help='strike direction of the ring plane, clockwise from north,' 

4505 ' in [deg]') 

4506 

4507 dip = Float.T( 

4508 default=0.0, 

4509 help='dip angle of the ring plane from horizontal in [deg]') 

4510 

4511 npointsources = Int.T( 

4512 default=360, 

4513 help='number of point sources to use') 

4514 

4515 discretized_source_class = meta.DiscretizedMTSource 

4516 

4517 def base_key(self): 

4518 return Source.base_key(self) + ( 

4519 self.strike, self.dip, self.diameter, self.npointsources) 

4520 

4521 def get_factor(self): 

4522 return self.sign * self.moment 

4523 

4524 def discretize_basesource(self, store=None, target=None): 

4525 n = self.npointsources 

4526 phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False) 

4527 

4528 points = num.zeros((n, 3)) 

4529 points[:, 0] = num.cos(phi) * 0.5 * self.diameter 

4530 points[:, 1] = num.sin(phi) * 0.5 * self.diameter 

4531 

4532 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0) 

4533 points = num.dot(rotmat.T, points.T).T # !!! ? 

4534 

4535 points[:, 0] += self.north_shift 

4536 points[:, 1] += self.east_shift 

4537 points[:, 2] += self.depth 

4538 

4539 m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90., 

4540 scalar_moment=1.0 / n).m()) 

4541 

4542 rotmats = num.transpose( 

4543 [[num.cos(phi), num.sin(phi), num.zeros(n)], 

4544 [-num.sin(phi), num.cos(phi), num.zeros(n)], 

4545 [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1)) 

4546 

4547 ms = num.zeros((n, 3, 3)) 

4548 for i in range(n): 

4549 mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i])) 

4550 ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat)) 

4551 

4552 m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2], 

4553 ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T 

4554 

4555 times, amplitudes = self.effective_stf_pre().discretize_t( 

4556 store.config.deltat, self.time) 

4557 

4558 nt = times.size 

4559 

4560 return meta.DiscretizedMTSource( 

4561 times=num.tile(times, n), 

4562 lat=self.lat, 

4563 lon=self.lon, 

4564 north_shifts=num.repeat(points[:, 0], nt), 

4565 east_shifts=num.repeat(points[:, 1], nt), 

4566 depths=num.repeat(points[:, 2], nt), 

4567 m6s=num.repeat(m6s, nt, axis=0) * num.tile( 

4568 amplitudes, n)[:, num.newaxis]) 

4569 

4570 

4571class CombiSource(Source): 

4572 ''' 

4573 Composite source model. 

4574 ''' 

4575 

4576 discretized_source_class = meta.DiscretizedMTSource 

4577 

4578 subsources = List.T(Source.T()) 

4579 

4580 def __init__(self, subsources=[], **kwargs): 

4581 if not subsources: 

4582 raise BadRequest( 

4583 'Need at least one sub-source to create a CombiSource object.') 

4584 

4585 lats = num.array( 

4586 [subsource.lat for subsource in subsources], dtype=float) 

4587 lons = num.array( 

4588 [subsource.lon for subsource in subsources], dtype=float) 

4589 

4590 lat, lon = lats[0], lons[0] 

4591 if not num.all(lats == lat) and num.all(lons == lon): 

4592 subsources = [s.clone() for s in subsources] 

4593 for subsource in subsources[1:]: 

4594 subsource.set_origin(lat, lon) 

4595 

4596 depth = float(num.mean([p.depth for p in subsources])) 

4597 time = float(num.mean([p.time for p in subsources])) 

4598 north_shift = float(num.mean([p.north_shift for p in subsources])) 

4599 east_shift = float(num.mean([p.east_shift for p in subsources])) 

4600 kwargs.update( 

4601 time=time, 

4602 lat=float(lat), 

4603 lon=float(lon), 

4604 north_shift=north_shift, 

4605 east_shift=east_shift, 

4606 depth=depth) 

4607 

4608 Source.__init__(self, subsources=subsources, **kwargs) 

4609 

4610 def get_factor(self): 

4611 return 1.0 

4612 

4613 def discretize_basesource(self, store, target=None): 

4614 dsources = [] 

4615 for sf in self.subsources: 

4616 ds = sf.discretize_basesource(store, target) 

4617 ds.m6s *= sf.get_factor() 

4618 dsources.append(ds) 

4619 

4620 return meta.DiscretizedMTSource.combine(dsources) 

4621 

4622 

4623class SFSource(Source): 

4624 ''' 

4625 A single force point source. 

4626 

4627 Supported GF schemes: `'elastic5'`. 

4628 ''' 

4629 

4630 discretized_source_class = meta.DiscretizedSFSource 

4631 

4632 fn = Float.T( 

4633 default=0., 

4634 help='northward component of single force [N]') 

4635 

4636 fe = Float.T( 

4637 default=0., 

4638 help='eastward component of single force [N]') 

4639 

4640 fd = Float.T( 

4641 default=0., 

4642 help='downward component of single force [N]') 

4643 

4644 def __init__(self, **kwargs): 

4645 Source.__init__(self, **kwargs) 

4646 

4647 def base_key(self): 

4648 return Source.base_key(self) + (self.fn, self.fe, self.fd) 

4649 

4650 def get_factor(self): 

4651 return 1.0 

4652 

4653 def discretize_basesource(self, store, target=None): 

4654 times, amplitudes = self.effective_stf_pre().discretize_t( 

4655 store.config.deltat, self.time) 

4656 forces = amplitudes[:, num.newaxis] * num.array( 

4657 [[self.fn, self.fe, self.fd]], dtype=float) 

4658 

4659 return meta.DiscretizedSFSource(forces=forces, 

4660 **self._dparams_base_repeated(times)) 

4661 

4662 def pyrocko_event(self, store=None, target=None, **kwargs): 

4663 return Source.pyrocko_event( 

4664 self, store, target, 

4665 **kwargs) 

4666 

4667 @classmethod 

4668 def from_pyrocko_event(cls, ev, **kwargs): 

4669 d = {} 

4670 d.update(kwargs) 

4671 return super(SFSource, cls).from_pyrocko_event(ev, **d) 

4672 

4673 

4674class PorePressurePointSource(Source): 

4675 ''' 

4676 Excess pore pressure point source. 

4677 

4678 For poro-elastic initial value problem where an excess pore pressure is 

4679 brought into a small source volume. 

4680 ''' 

4681 

4682 discretized_source_class = meta.DiscretizedPorePressureSource 

4683 

4684 pp = Float.T( 

4685 default=1.0, 

4686 help='initial excess pore pressure in [Pa]') 

4687 

4688 def base_key(self): 

4689 return Source.base_key(self) 

4690 

4691 def get_factor(self): 

4692 return self.pp 

4693 

4694 def discretize_basesource(self, store, target=None): 

4695 return meta.DiscretizedPorePressureSource(pp=arr(1.0), 

4696 **self._dparams_base()) 

4697 

4698 

4699class PorePressureLineSource(Source): 

4700 ''' 

4701 Excess pore pressure line source. 

4702 

4703 The line source is centered at (north_shift, east_shift, depth). 

4704 ''' 

4705 

4706 discretized_source_class = meta.DiscretizedPorePressureSource 

4707 

4708 pp = Float.T( 

4709 default=1.0, 

4710 help='initial excess pore pressure in [Pa]') 

4711 

4712 length = Float.T( 

4713 default=0.0, 

4714 help='length of the line source [m]') 

4715 

4716 azimuth = Float.T( 

4717 default=0.0, 

4718 help='azimuth direction, clockwise from north [deg]') 

4719 

4720 dip = Float.T( 

4721 default=90., 

4722 help='dip direction, downward from horizontal [deg]') 

4723 

4724 def base_key(self): 

4725 return Source.base_key(self) + (self.azimuth, self.dip, self.length) 

4726 

4727 def get_factor(self): 

4728 return self.pp 

4729 

4730 def discretize_basesource(self, store, target=None): 

4731 

4732 n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1 

4733 

4734 a = num.linspace(-0.5 * self.length, 0.5 * self.length, n) 

4735 

4736 sa = math.sin(self.azimuth * d2r) 

4737 ca = math.cos(self.azimuth * d2r) 

4738 sd = math.sin(self.dip * d2r) 

4739 cd = math.cos(self.dip * d2r) 

4740 

4741 points = num.zeros((n, 3)) 

4742 points[:, 0] = self.north_shift + a * ca * cd 

4743 points[:, 1] = self.east_shift + a * sa * cd 

4744 points[:, 2] = self.depth + a * sd 

4745 

4746 return meta.DiscretizedPorePressureSource( 

4747 times=util.num_full(n, self.time), 

4748 lat=self.lat, 

4749 lon=self.lon, 

4750 north_shifts=points[:, 0], 

4751 east_shifts=points[:, 1], 

4752 depths=points[:, 2], 

4753 pp=num.ones(n) / n) 

4754 

4755 

4756class Request(Object): 

4757 ''' 

4758 Synthetic seismogram computation request. 

4759 

4760 :: 

4761 

4762 Request(**kwargs) 

4763 Request(sources, targets, **kwargs) 

4764 ''' 

4765 

4766 sources = List.T( 

4767 Source.T(), 

4768 help='list of sources for which to produce synthetics.') 

4769 

4770 targets = List.T( 

4771 Target.T(), 

4772 help='list of targets for which to produce synthetics.') 

4773 

4774 @classmethod 

4775 def args2kwargs(cls, args): 

4776 if len(args) not in (0, 2, 3): 

4777 raise BadRequest('Invalid arguments.') 

4778 

4779 if len(args) == 2: 

4780 return dict(sources=args[0], targets=args[1]) 

4781 else: 

4782 return {} 

4783 

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

4785 kwargs.update(self.args2kwargs(args)) 

4786 sources = kwargs.pop('sources', []) 

4787 targets = kwargs.pop('targets', []) 

4788 

4789 if isinstance(sources, Source): 

4790 sources = [sources] 

4791 

4792 if isinstance(targets, Target) or isinstance(targets, StaticTarget): 

4793 targets = [targets] 

4794 

4795 Object.__init__(self, sources=sources, targets=targets, **kwargs) 

4796 

4797 @property 

4798 def targets_dynamic(self): 

4799 return [t for t in self.targets if isinstance(t, Target)] 

4800 

4801 @property 

4802 def targets_static(self): 

4803 return [t for t in self.targets if isinstance(t, StaticTarget)] 

4804 

4805 @property 

4806 def has_dynamic(self): 

4807 return True if len(self.targets_dynamic) > 0 else False 

4808 

4809 @property 

4810 def has_statics(self): 

4811 return True if len(self.targets_static) > 0 else False 

4812 

4813 def subsources_map(self): 

4814 m = defaultdict(list) 

4815 for source in self.sources: 

4816 m[source.base_key()].append(source) 

4817 

4818 return m 

4819 

4820 def subtargets_map(self): 

4821 m = defaultdict(list) 

4822 for target in self.targets: 

4823 m[target.base_key()].append(target) 

4824 

4825 return m 

4826 

4827 def subrequest_map(self): 

4828 ms = self.subsources_map() 

4829 mt = self.subtargets_map() 

4830 m = {} 

4831 for (ks, ls) in ms.items(): 

4832 for (kt, lt) in mt.items(): 

4833 m[ks, kt] = (ls, lt) 

4834 

4835 return m 

4836 

4837 

4838class ProcessingStats(Object): 

4839 t_perc_get_store_and_receiver = Float.T(default=0.) 

4840 t_perc_discretize_source = Float.T(default=0.) 

4841 t_perc_make_base_seismogram = Float.T(default=0.) 

4842 t_perc_make_same_span = Float.T(default=0.) 

4843 t_perc_post_process = Float.T(default=0.) 

4844 t_perc_optimize = Float.T(default=0.) 

4845 t_perc_stack = Float.T(default=0.) 

4846 t_perc_static_get_store = Float.T(default=0.) 

4847 t_perc_static_discretize_basesource = Float.T(default=0.) 

4848 t_perc_static_sum_statics = Float.T(default=0.) 

4849 t_perc_static_post_process = Float.T(default=0.) 

4850 t_wallclock = Float.T(default=0.) 

4851 t_cpu = Float.T(default=0.) 

4852 n_read_blocks = Int.T(default=0) 

4853 n_results = Int.T(default=0) 

4854 n_subrequests = Int.T(default=0) 

4855 n_stores = Int.T(default=0) 

4856 n_records_stacked = Int.T(default=0) 

4857 

4858 

4859class Response(Object): 

4860 ''' 

4861 Resonse object to a synthetic seismogram computation request. 

4862 ''' 

4863 

4864 request = Request.T() 

4865 results_list = List.T(List.T(meta.SeismosizerResult.T())) 

4866 stats = ProcessingStats.T() 

4867 

4868 def pyrocko_traces(self): 

4869 ''' 

4870 Return a list of requested 

4871 :class:`~pyrocko.trace.Trace` instances. 

4872 ''' 

4873 

4874 traces = [] 

4875 for results in self.results_list: 

4876 for result in results: 

4877 if not isinstance(result, meta.Result): 

4878 continue 

4879 traces.append(result.trace.pyrocko_trace()) 

4880 

4881 return traces 

4882 

4883 def kite_scenes(self): 

4884 ''' 

4885 Return a list of requested 

4886 :class:`kite.Scene` instances. 

4887 ''' 

4888 kite_scenes = [] 

4889 for results in self.results_list: 

4890 for result in results: 

4891 if isinstance(result, meta.KiteSceneResult): 

4892 sc = result.get_scene() 

4893 kite_scenes.append(sc) 

4894 

4895 return kite_scenes 

4896 

4897 def static_results(self): 

4898 ''' 

4899 Return a list of requested 

4900 :class:`~pyrocko.gf.meta.StaticResult` instances. 

4901 ''' 

4902 statics = [] 

4903 for results in self.results_list: 

4904 for result in results: 

4905 if not isinstance(result, meta.StaticResult): 

4906 continue 

4907 statics.append(result) 

4908 

4909 return statics 

4910 

4911 def iter_results(self, get='pyrocko_traces'): 

4912 ''' 

4913 Generator function to iterate over results of request. 

4914 

4915 Yields associated :py:class:`Source`, 

4916 :class:`~pyrocko.gf.targets.Target`, 

4917 :class:`~pyrocko.trace.Trace` instances in each iteration. 

4918 ''' 

4919 

4920 for isource, source in enumerate(self.request.sources): 

4921 for itarget, target in enumerate(self.request.targets): 

4922 result = self.results_list[isource][itarget] 

4923 if get == 'pyrocko_traces': 

4924 yield source, target, result.trace.pyrocko_trace() 

4925 elif get == 'results': 

4926 yield source, target, result 

4927 

4928 def snuffle(self, **kwargs): 

4929 ''' 

4930 Open *snuffler* with requested traces. 

4931 ''' 

4932 

4933 trace.snuffle(self.pyrocko_traces(), **kwargs) 

4934 

4935 

4936class Engine(Object): 

4937 ''' 

4938 Base class for synthetic seismogram calculators. 

4939 ''' 

4940 

4941 def get_store_ids(self): 

4942 ''' 

4943 Get list of available GF store IDs 

4944 ''' 

4945 

4946 return [] 

4947 

4948 

4949class Rule(object): 

4950 pass 

4951 

4952 

4953class VectorRule(Rule): 

4954 

4955 def __init__(self, quantity, differentiate=0, integrate=0): 

4956 self.components = [quantity + '.' + c for c in 'ned'] 

4957 self.differentiate = differentiate 

4958 self.integrate = integrate 

4959 

4960 def required_components(self, target): 

4961 n, e, d = self.components 

4962 sa, ca, sd, cd = target.get_sin_cos_factors() 

4963 

4964 comps = [] 

4965 if nonzero(ca * cd): 

4966 comps.append(n) 

4967 

4968 if nonzero(sa * cd): 

4969 comps.append(e) 

4970 

4971 if nonzero(sd): 

4972 comps.append(d) 

4973 

4974 return tuple(comps) 

4975 

4976 def apply_(self, target, base_seismogram): 

4977 n, e, d = self.components 

4978 sa, ca, sd, cd = target.get_sin_cos_factors() 

4979 

4980 if nonzero(ca * cd): 

4981 data = base_seismogram[n].data * (ca * cd) 

4982 deltat = base_seismogram[n].deltat 

4983 else: 

4984 data = 0.0 

4985 

4986 if nonzero(sa * cd): 

4987 data = data + base_seismogram[e].data * (sa * cd) 

4988 deltat = base_seismogram[e].deltat 

4989 

4990 if nonzero(sd): 

4991 data = data + base_seismogram[d].data * sd 

4992 deltat = base_seismogram[d].deltat 

4993 

4994 if self.differentiate: 

4995 data = util.diff_fd(self.differentiate, 4, deltat, data) 

4996 

4997 if self.integrate: 

4998 raise NotImplementedError('Integration is not implemented yet.') 

4999 

5000 return data 

5001 

5002 

5003class HorizontalVectorRule(Rule): 

5004 

5005 def __init__(self, quantity, differentiate=0, integrate=0): 

5006 self.components = [quantity + '.' + c for c in 'ne'] 

5007 self.differentiate = differentiate 

5008 self.integrate = integrate 

5009 

5010 def required_components(self, target): 

5011 n, e = self.components 

5012 sa, ca, _, _ = target.get_sin_cos_factors() 

5013 

5014 comps = [] 

5015 if nonzero(ca): 

5016 comps.append(n) 

5017 

5018 if nonzero(sa): 

5019 comps.append(e) 

5020 

5021 return tuple(comps) 

5022 

5023 def apply_(self, target, base_seismogram): 

5024 n, e = self.components 

5025 sa, ca, _, _ = target.get_sin_cos_factors() 

5026 

5027 if nonzero(ca): 

5028 data = base_seismogram[n].data * ca 

5029 else: 

5030 data = 0.0 

5031 

5032 if nonzero(sa): 

5033 data = data + base_seismogram[e].data * sa 

5034 

5035 if self.differentiate: 

5036 deltat = base_seismogram[e].deltat 

5037 data = util.diff_fd(self.differentiate, 4, deltat, data) 

5038 

5039 if self.integrate: 

5040 raise NotImplementedError('Integration is not implemented yet.') 

5041 

5042 return data 

5043 

5044 

5045class ScalarRule(Rule): 

5046 

5047 def __init__(self, quantity, differentiate=0): 

5048 self.c = quantity 

5049 

5050 def required_components(self, target): 

5051 return (self.c, ) 

5052 

5053 def apply_(self, target, base_seismogram): 

5054 data = base_seismogram[self.c].data.copy() 

5055 deltat = base_seismogram[self.c].deltat 

5056 if self.differentiate: 

5057 data = util.diff_fd(self.differentiate, 4, deltat, data) 

5058 

5059 return data 

5060 

5061 

5062class StaticDisplacement(Rule): 

5063 

5064 def required_components(self, target): 

5065 return tuple(['displacement.%s' % c for c in list('ned')]) 

5066 

5067 def apply_(self, target, base_statics): 

5068 if isinstance(target, SatelliteTarget): 

5069 los_fac = target.get_los_factors() 

5070 base_statics['displacement.los'] =\ 

5071 (los_fac[:, 0] * -base_statics['displacement.d'] + 

5072 los_fac[:, 1] * base_statics['displacement.e'] + 

5073 los_fac[:, 2] * base_statics['displacement.n']) 

5074 return base_statics 

5075 

5076 

5077channel_rules = { 

5078 'displacement': [VectorRule('displacement')], 

5079 'rotation': [VectorRule('rotation')], 

5080 'velocity': [ 

5081 VectorRule('velocity'), 

5082 VectorRule('displacement', differentiate=1)], 

5083 'acceleration': [ 

5084 VectorRule('acceleration'), 

5085 VectorRule('velocity', differentiate=1), 

5086 VectorRule('displacement', differentiate=2)], 

5087 'pore_pressure': [ScalarRule('pore_pressure')], 

5088 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')], 

5089 'darcy_velocity': [VectorRule('darcy_velocity')], 

5090} 

5091 

5092static_rules = { 

5093 'displacement': [StaticDisplacement()] 

5094} 

5095 

5096 

5097class OutOfBoundsContext(Object): 

5098 source = Source.T() 

5099 target = Target.T() 

5100 distance = Float.T() 

5101 components = List.T(String.T()) 

5102 

5103 

5104def process_dynamic_timeseries(work, psources, ptargets, engine, nthreads=0): 

5105 dsource_cache = {} 

5106 tcounters = list(range(6)) 

5107 

5108 store_ids = set() 

5109 sources = set() 

5110 targets = set() 

5111 

5112 for itarget, target in enumerate(ptargets): 

5113 target._id = itarget 

5114 

5115 for w in work: 

5116 _, _, isources, itargets = w 

5117 

5118 sources.update([psources[isource] for isource in isources]) 

5119 targets.update([ptargets[itarget] for itarget in itargets]) 

5120 

5121 store_ids = set([t.store_id for t in targets]) 

5122 

5123 for isource, source in enumerate(psources): 

5124 

5125 components = set() 

5126 for itarget, target in enumerate(targets): 

5127 rule = engine.get_rule(source, target) 

5128 components.update(rule.required_components(target)) 

5129 

5130 for store_id in store_ids: 

5131 store_targets = [t for t in targets if t.store_id == store_id] 

5132 

5133 sample_rates = set([t.sample_rate for t in store_targets]) 

5134 interpolations = set([t.interpolation for t in store_targets]) 

5135 

5136 base_seismograms = [] 

5137 store_targets_out = [] 

5138 

5139 for samp_rate in sample_rates: 

5140 for interp in interpolations: 

5141 engine_targets = [ 

5142 t for t in store_targets if t.sample_rate == samp_rate 

5143 and t.interpolation == interp] 

5144 

5145 if not engine_targets: 

5146 continue 

5147 

5148 store_targets_out += engine_targets 

5149 

5150 base_seismograms += engine.base_seismograms( 

5151 source, 

5152 engine_targets, 

5153 components, 

5154 dsource_cache, 

5155 nthreads) 

5156 

5157 for iseis, seismogram in enumerate(base_seismograms): 

5158 for tr in seismogram.values(): 

5159 if tr.err != store.SeismosizerErrorEnum.SUCCESS: 

5160 e = SeismosizerError( 

5161 'Seismosizer failed with return code %i\n%s' % ( 

5162 tr.err, str( 

5163 OutOfBoundsContext( 

5164 source=source, 

5165 target=store_targets[iseis], 

5166 distance=source.distance_to( 

5167 store_targets[iseis]), 

5168 components=components)))) 

5169 raise e 

5170 

5171 for seismogram, target in zip(base_seismograms, store_targets_out): 

5172 

5173 try: 

5174 result = engine._post_process_dynamic( 

5175 seismogram, source, target) 

5176 except SeismosizerError as e: 

5177 result = e 

5178 

5179 yield (isource, target._id, result), tcounters 

5180 

5181 

5182def process_dynamic(work, psources, ptargets, engine, nthreads=0): 

5183 dsource_cache = {} 

5184 

5185 for w in work: 

5186 _, _, isources, itargets = w 

5187 

5188 sources = [psources[isource] for isource in isources] 

5189 targets = [ptargets[itarget] for itarget in itargets] 

5190 

5191 components = set() 

5192 for target in targets: 

5193 rule = engine.get_rule(sources[0], target) 

5194 components.update(rule.required_components(target)) 

5195 

5196 for isource, source in zip(isources, sources): 

5197 for itarget, target in zip(itargets, targets): 

5198 

5199 try: 

5200 base_seismogram, tcounters = engine.base_seismogram( 

5201 source, target, components, dsource_cache, nthreads) 

5202 except meta.OutOfBounds as e: 

5203 e.context = OutOfBoundsContext( 

5204 source=sources[0], 

5205 target=targets[0], 

5206 distance=sources[0].distance_to(targets[0]), 

5207 components=components) 

5208 raise 

5209 

5210 n_records_stacked = 0 

5211 t_optimize = 0.0 

5212 t_stack = 0.0 

5213 

5214 for _, tr in base_seismogram.items(): 

5215 n_records_stacked += tr.n_records_stacked 

5216 t_optimize += tr.t_optimize 

5217 t_stack += tr.t_stack 

5218 

5219 try: 

5220 result = engine._post_process_dynamic( 

5221 base_seismogram, source, target) 

5222 result.n_records_stacked = n_records_stacked 

5223 result.n_shared_stacking = len(sources) *\ 

5224 len(targets) 

5225 result.t_optimize = t_optimize 

5226 result.t_stack = t_stack 

5227 except SeismosizerError as e: 

5228 result = e 

5229 

5230 tcounters.append(xtime()) 

5231 yield (isource, itarget, result), tcounters 

5232 

5233 

5234def process_static(work, psources, ptargets, engine, nthreads=0): 

5235 for w in work: 

5236 _, _, isources, itargets = w 

5237 

5238 sources = [psources[isource] for isource in isources] 

5239 targets = [ptargets[itarget] for itarget in itargets] 

5240 

5241 for isource, source in zip(isources, sources): 

5242 for itarget, target in zip(itargets, targets): 

5243 components = engine.get_rule(source, target)\ 

5244 .required_components(target) 

5245 

5246 try: 

5247 base_statics, tcounters = engine.base_statics( 

5248 source, target, components, nthreads) 

5249 except meta.OutOfBounds as e: 

5250 e.context = OutOfBoundsContext( 

5251 source=sources[0], 

5252 target=targets[0], 

5253 distance=float('nan'), 

5254 components=components) 

5255 raise 

5256 result = engine._post_process_statics( 

5257 base_statics, source, target) 

5258 tcounters.append(xtime()) 

5259 

5260 yield (isource, itarget, result), tcounters 

5261 

5262 

5263class LocalEngine(Engine): 

5264 ''' 

5265 Offline synthetic seismogram calculator. 

5266 

5267 :param use_env: if ``True``, fill :py:attr:`store_superdirs` and 

5268 :py:attr:`store_dirs` with paths set in environment variables 

5269 GF_STORE_SUPERDIRS and GF_STORE_DIRS. 

5270 :param use_config: if ``True``, fill :py:attr:`store_superdirs` and 

5271 :py:attr:`store_dirs` with paths set in the user's config file. 

5272 

5273 The config file can be found at :file:`~/.pyrocko/config.pf` 

5274 

5275 .. code-block :: python 

5276 

5277 gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/'] 

5278 gf_store_superdirs: ['/home/pyrocko/gf_stores/'] 

5279 ''' 

5280 

5281 store_superdirs = List.T( 

5282 String.T(), 

5283 help="directories which are searched for Green's function stores") 

5284 

5285 store_dirs = List.T( 

5286 String.T(), 

5287 help="additional individual Green's function store directories") 

5288 

5289 default_store_id = String.T( 

5290 optional=True, 

5291 help='default store ID to be used when a request does not provide ' 

5292 'one') 

5293 

5294 def __init__(self, **kwargs): 

5295 use_env = kwargs.pop('use_env', False) 

5296 use_config = kwargs.pop('use_config', False) 

5297 Engine.__init__(self, **kwargs) 

5298 if use_env: 

5299 env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '') 

5300 env_store_dirs = os.environ.get('GF_STORE_DIRS', '') 

5301 if env_store_superdirs: 

5302 self.store_superdirs.extend(env_store_superdirs.split(':')) 

5303 

5304 if env_store_dirs: 

5305 self.store_dirs.extend(env_store_dirs.split(':')) 

5306 

5307 if use_config: 

5308 c = config.config() 

5309 self.store_superdirs.extend(c.gf_store_superdirs) 

5310 self.store_dirs.extend(c.gf_store_dirs) 

5311 

5312 self._check_store_dirs_type() 

5313 self._id_to_store_dir = {} 

5314 self._open_stores = {} 

5315 self._effective_default_store_id = None 

5316 

5317 def _check_store_dirs_type(self): 

5318 for sdir in ['store_dirs', 'store_superdirs']: 

5319 if not isinstance(self.__getattribute__(sdir), list): 

5320 raise TypeError('{} of {} is not of type list'.format( 

5321 sdir, self.__class__.__name__)) 

5322 

5323 def _get_store_id(self, store_dir): 

5324 store_ = store.Store(store_dir) 

5325 store_id = store_.config.id 

5326 store_.close() 

5327 return store_id 

5328 

5329 def _looks_like_store_dir(self, store_dir): 

5330 return os.path.isdir(store_dir) and \ 

5331 all(os.path.isfile(pjoin(store_dir, x)) for x in 

5332 ('index', 'traces', 'config')) 

5333 

5334 def iter_store_dirs(self): 

5335 store_dirs = set() 

5336 for d in self.store_superdirs: 

5337 if not os.path.exists(d): 

5338 logger.warning('store_superdir not available: %s' % d) 

5339 continue 

5340 

5341 for entry in os.listdir(d): 

5342 store_dir = os.path.realpath(pjoin(d, entry)) 

5343 if self._looks_like_store_dir(store_dir): 

5344 store_dirs.add(store_dir) 

5345 

5346 for store_dir in self.store_dirs: 

5347 store_dirs.add(os.path.realpath(store_dir)) 

5348 

5349 return store_dirs 

5350 

5351 def _scan_stores(self): 

5352 for store_dir in self.iter_store_dirs(): 

5353 store_id = self._get_store_id(store_dir) 

5354 if store_id not in self._id_to_store_dir: 

5355 self._id_to_store_dir[store_id] = store_dir 

5356 else: 

5357 if store_dir != self._id_to_store_dir[store_id]: 

5358 raise DuplicateStoreId( 

5359 'GF store ID %s is used in (at least) two ' 

5360 'different stores. Locations are: %s and %s' % 

5361 (store_id, self._id_to_store_dir[store_id], store_dir)) 

5362 

5363 def get_store_dir(self, store_id): 

5364 ''' 

5365 Lookup directory given a GF store ID. 

5366 ''' 

5367 

5368 if store_id not in self._id_to_store_dir: 

5369 self._scan_stores() 

5370 

5371 if store_id not in self._id_to_store_dir: 

5372 raise NoSuchStore(store_id, self.iter_store_dirs()) 

5373 

5374 return self._id_to_store_dir[store_id] 

5375 

5376 def get_store_ids(self): 

5377 ''' 

5378 Get list of available store IDs. 

5379 ''' 

5380 

5381 self._scan_stores() 

5382 return sorted(self._id_to_store_dir.keys()) 

5383 

5384 def effective_default_store_id(self): 

5385 if self._effective_default_store_id is None: 

5386 if self.default_store_id is None: 

5387 store_ids = self.get_store_ids() 

5388 if len(store_ids) == 1: 

5389 self._effective_default_store_id = self.get_store_ids()[0] 

5390 else: 

5391 raise NoDefaultStoreSet() 

5392 else: 

5393 self._effective_default_store_id = self.default_store_id 

5394 

5395 return self._effective_default_store_id 

5396 

5397 def get_store(self, store_id=None): 

5398 ''' 

5399 Get a store from the engine. 

5400 

5401 :param store_id: identifier of the store (optional) 

5402 :returns: :py:class:`~pyrocko.gf.store.Store` object 

5403 

5404 If no ``store_id`` is provided the store 

5405 associated with the :py:gattr:`default_store_id` is returned. 

5406 Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is 

5407 undefined. 

5408 ''' 

5409 

5410 if store_id is None: 

5411 store_id = self.effective_default_store_id() 

5412 

5413 if store_id not in self._open_stores: 

5414 store_dir = self.get_store_dir(store_id) 

5415 self._open_stores[store_id] = store.Store(store_dir) 

5416 

5417 return self._open_stores[store_id] 

5418 

5419 def get_store_config(self, store_id): 

5420 store = self.get_store(store_id) 

5421 return store.config 

5422 

5423 def get_store_extra(self, store_id, key): 

5424 store = self.get_store(store_id) 

5425 return store.get_extra(key) 

5426 

5427 def close_cashed_stores(self): 

5428 ''' 

5429 Close and remove ids from cashed stores. 

5430 ''' 

5431 store_ids = [] 

5432 for store_id, store_ in self._open_stores.items(): 

5433 store_.close() 

5434 store_ids.append(store_id) 

5435 

5436 for store_id in store_ids: 

5437 self._open_stores.pop(store_id) 

5438 

5439 def get_rule(self, source, target): 

5440 cprovided = self.get_store(target.store_id).get_provided_components() 

5441 

5442 if isinstance(target, StaticTarget): 

5443 quantity = target.quantity 

5444 available_rules = static_rules 

5445 elif isinstance(target, Target): 

5446 quantity = target.effective_quantity() 

5447 available_rules = channel_rules 

5448 

5449 try: 

5450 for rule in available_rules[quantity]: 

5451 cneeded = rule.required_components(target) 

5452 if all(c in cprovided for c in cneeded): 

5453 return rule 

5454 

5455 except KeyError: 

5456 pass 

5457 

5458 raise BadRequest( 

5459 'No rule to calculate "%s" with GFs from store "%s" ' 

5460 'for source model "%s".' % ( 

5461 target.effective_quantity(), 

5462 target.store_id, 

5463 source.__class__.__name__)) 

5464 

5465 def _cached_discretize_basesource(self, source, store, cache, target): 

5466 if (source, store) not in cache: 

5467 cache[source, store] = source.discretize_basesource(store, target) 

5468 

5469 return cache[source, store] 

5470 

5471 def base_seismograms(self, source, targets, components, dsource_cache, 

5472 nthreads=0): 

5473 

5474 target = targets[0] 

5475 

5476 interp = set([t.interpolation for t in targets]) 

5477 if len(interp) > 1: 

5478 raise BadRequest('Targets have different interpolation schemes.') 

5479 

5480 rates = set([t.sample_rate for t in targets]) 

5481 if len(rates) > 1: 

5482 raise BadRequest('Targets have different sample rates.') 

5483 

5484 store_ = self.get_store(target.store_id) 

5485 receivers = [t.receiver(store_) for t in targets] 

5486 

5487 if target.sample_rate is not None: 

5488 deltat = 1. / target.sample_rate 

5489 rate = target.sample_rate 

5490 else: 

5491 deltat = None 

5492 rate = store_.config.sample_rate 

5493 

5494 tmin = num.fromiter( 

5495 (t.tmin for t in targets), dtype=float, count=len(targets)) 

5496 tmax = num.fromiter( 

5497 (t.tmax for t in targets), dtype=float, count=len(targets)) 

5498 

5499 mask = num.logical_and(num.isfinite(tmin), num.isfinite(tmax)) 

5500 

5501 itmin = num.zeros_like(tmin, dtype=num.int64) 

5502 itmax = num.zeros_like(tmin, dtype=num.int64) 

5503 nsamples = num.full_like(tmin, -1, dtype=num.int64) 

5504 

5505 itmin[mask] = num.floor(tmin[mask] * rate).astype(num.int64) 

5506 itmax[mask] = num.ceil(tmax[mask] * rate).astype(num.int64) 

5507 nsamples = itmax - itmin + 1 

5508 nsamples[num.logical_not(mask)] = -1 

5509 

5510 base_source = self._cached_discretize_basesource( 

5511 source, store_, dsource_cache, target) 

5512 

5513 base_seismograms = store_.calc_seismograms( 

5514 base_source, receivers, components, 

5515 deltat=deltat, 

5516 itmin=itmin, nsamples=nsamples, 

5517 interpolation=target.interpolation, 

5518 optimization=target.optimization, 

5519 nthreads=nthreads) 

5520 

5521 for i, base_seismogram in enumerate(base_seismograms): 

5522 base_seismograms[i] = store.make_same_span(base_seismogram) 

5523 

5524 return base_seismograms 

5525 

5526 def base_seismogram(self, source, target, components, dsource_cache, 

5527 nthreads): 

5528 

5529 tcounters = [xtime()] 

5530 

5531 store_ = self.get_store(target.store_id) 

5532 receiver = target.receiver(store_) 

5533 

5534 if target.tmin and target.tmax is not None: 

5535 rate = store_.config.sample_rate 

5536 itmin = int(num.floor(target.tmin * rate)) 

5537 itmax = int(num.ceil(target.tmax * rate)) 

5538 nsamples = itmax - itmin + 1 

5539 else: 

5540 itmin = None 

5541 nsamples = None 

5542 

5543 tcounters.append(xtime()) 

5544 base_source = self._cached_discretize_basesource( 

5545 source, store_, dsource_cache, target) 

5546 

5547 tcounters.append(xtime()) 

5548 

5549 if target.sample_rate is not None: 

5550 deltat = 1. / target.sample_rate 

5551 else: 

5552 deltat = None 

5553 

5554 base_seismogram = store_.seismogram( 

5555 base_source, receiver, components, 

5556 deltat=deltat, 

5557 itmin=itmin, nsamples=nsamples, 

5558 interpolation=target.interpolation, 

5559 optimization=target.optimization, 

5560 nthreads=nthreads) 

5561 

5562 tcounters.append(xtime()) 

5563 

5564 base_seismogram = store.make_same_span(base_seismogram) 

5565 

5566 tcounters.append(xtime()) 

5567 

5568 return base_seismogram, tcounters 

5569 

5570 def base_statics(self, source, target, components, nthreads): 

5571 tcounters = [xtime()] 

5572 store_ = self.get_store(target.store_id) 

5573 

5574 if target.tsnapshot is not None: 

5575 rate = store_.config.sample_rate 

5576 itsnapshot = int(num.floor(target.tsnapshot * rate)) 

5577 else: 

5578 itsnapshot = None 

5579 tcounters.append(xtime()) 

5580 

5581 base_source = source.discretize_basesource(store_, target=target) 

5582 

5583 tcounters.append(xtime()) 

5584 

5585 base_statics = store_.statics( 

5586 base_source, 

5587 target, 

5588 itsnapshot, 

5589 components, 

5590 target.interpolation, 

5591 nthreads) 

5592 

5593 tcounters.append(xtime()) 

5594 

5595 return base_statics, tcounters 

5596 

5597 def _post_process_dynamic(self, base_seismogram, source, target): 

5598 base_any = next(iter(base_seismogram.values())) 

5599 deltat = base_any.deltat 

5600 itmin = base_any.itmin 

5601 

5602 rule = self.get_rule(source, target) 

5603 data = rule.apply_(target, base_seismogram) 

5604 

5605 factor = source.get_factor() * target.get_factor() 

5606 if factor != 1.0: 

5607 data = data * factor 

5608 

5609 stf = source.effective_stf_post() 

5610 

5611 times, amplitudes = stf.discretize_t( 

5612 deltat, 0.0) 

5613 

5614 # repeat end point to prevent boundary effects 

5615 padded_data = num.empty(data.size + amplitudes.size, dtype=float) 

5616 padded_data[:data.size] = data 

5617 padded_data[data.size:] = data[-1] 

5618 data = num.convolve(amplitudes, padded_data) 

5619 

5620 tmin = itmin * deltat + times[0] 

5621 

5622 tr = meta.SeismosizerTrace( 

5623 codes=target.codes, 

5624 data=data[:-amplitudes.size], 

5625 deltat=deltat, 

5626 tmin=tmin) 

5627 

5628 return target.post_process(self, source, tr) 

5629 

5630 def _post_process_statics(self, base_statics, source, starget): 

5631 rule = self.get_rule(source, starget) 

5632 data = rule.apply_(starget, base_statics) 

5633 

5634 factor = source.get_factor() 

5635 if factor != 1.0: 

5636 for v in data.values(): 

5637 v *= factor 

5638 

5639 return starget.post_process(self, source, base_statics) 

5640 

5641 def process(self, *args, **kwargs): 

5642 ''' 

5643 Process a request. 

5644 

5645 :: 

5646 

5647 process(**kwargs) 

5648 process(request, **kwargs) 

5649 process(sources, targets, **kwargs) 

5650 

5651 The request can be given a a :py:class:`Request` object, or such an 

5652 object is created using ``Request(**kwargs)`` for convenience. 

5653 

5654 :returns: :py:class:`Response` object 

5655 ''' 

5656 

5657 if len(args) not in (0, 1, 2): 

5658 raise BadRequest('Invalid arguments.') 

5659 

5660 if len(args) == 1: 

5661 kwargs['request'] = args[0] 

5662 

5663 elif len(args) == 2: 

5664 kwargs.update(Request.args2kwargs(args)) 

5665 

5666 request = kwargs.pop('request', None) 

5667 status_callback = kwargs.pop('status_callback', None) 

5668 calc_timeseries = kwargs.pop('calc_timeseries', True) 

5669 

5670 nprocs = kwargs.pop('nprocs', None) 

5671 nthreads = kwargs.pop('nthreads', 1) 

5672 if nprocs is not None: 

5673 nthreads = nprocs 

5674 

5675 if request is None: 

5676 request = Request(**kwargs) 

5677 

5678 if resource: 

5679 rs0 = resource.getrusage(resource.RUSAGE_SELF) 

5680 rc0 = resource.getrusage(resource.RUSAGE_CHILDREN) 

5681 tt0 = xtime() 

5682 

5683 # make sure stores are open before fork() 

5684 store_ids = set(target.store_id for target in request.targets) 

5685 for store_id in store_ids: 

5686 self.get_store(store_id) 

5687 

5688 source_index = dict((x, i) for (i, x) in 

5689 enumerate(request.sources)) 

5690 target_index = dict((x, i) for (i, x) in 

5691 enumerate(request.targets)) 

5692 

5693 m = request.subrequest_map() 

5694 

5695 skeys = sorted(m.keys(), key=cmp_to_key(cmp_none_aware)) 

5696 results_list = [] 

5697 

5698 for i in range(len(request.sources)): 

5699 results_list.append([None] * len(request.targets)) 

5700 

5701 tcounters_dyn_list = [] 

5702 tcounters_static_list = [] 

5703 nsub = len(skeys) 

5704 isub = 0 

5705 

5706 # Processing dynamic targets through 

5707 # parimap(process_subrequest_dynamic) 

5708 

5709 if calc_timeseries: 

5710 _process_dynamic = process_dynamic_timeseries 

5711 else: 

5712 _process_dynamic = process_dynamic 

5713 

5714 if request.has_dynamic: 

5715 work_dynamic = [ 

5716 (i, nsub, 

5717 [source_index[source] for source in m[k][0]], 

5718 [target_index[target] for target in m[k][1] 

5719 if not isinstance(target, StaticTarget)]) 

5720 for (i, k) in enumerate(skeys)] 

5721 

5722 for ii_results, tcounters_dyn in _process_dynamic( 

5723 work_dynamic, request.sources, request.targets, self, 

5724 nthreads): 

5725 

5726 tcounters_dyn_list.append(num.diff(tcounters_dyn)) 

5727 isource, itarget, result = ii_results 

5728 results_list[isource][itarget] = result 

5729 

5730 if status_callback: 

5731 status_callback(isub, nsub) 

5732 

5733 isub += 1 

5734 

5735 # Processing static targets through process_static 

5736 if request.has_statics: 

5737 work_static = [ 

5738 (i, nsub, 

5739 [source_index[source] for source in m[k][0]], 

5740 [target_index[target] for target in m[k][1] 

5741 if isinstance(target, StaticTarget)]) 

5742 for (i, k) in enumerate(skeys)] 

5743 

5744 for ii_results, tcounters_static in process_static( 

5745 work_static, request.sources, request.targets, self, 

5746 nthreads=nthreads): 

5747 

5748 tcounters_static_list.append(num.diff(tcounters_static)) 

5749 isource, itarget, result = ii_results 

5750 results_list[isource][itarget] = result 

5751 

5752 if status_callback: 

5753 status_callback(isub, nsub) 

5754 

5755 isub += 1 

5756 

5757 if status_callback: 

5758 status_callback(nsub, nsub) 

5759 

5760 tt1 = time.time() 

5761 if resource: 

5762 rs1 = resource.getrusage(resource.RUSAGE_SELF) 

5763 rc1 = resource.getrusage(resource.RUSAGE_CHILDREN) 

5764 

5765 s = ProcessingStats() 

5766 

5767 if request.has_dynamic: 

5768 tcumu_dyn = num.sum(num.vstack(tcounters_dyn_list), axis=0) 

5769 t_dyn = float(num.sum(tcumu_dyn)) 

5770 perc_dyn = map(float, tcumu_dyn / t_dyn * 100.) 

5771 (s.t_perc_get_store_and_receiver, 

5772 s.t_perc_discretize_source, 

5773 s.t_perc_make_base_seismogram, 

5774 s.t_perc_make_same_span, 

5775 s.t_perc_post_process) = perc_dyn 

5776 else: 

5777 t_dyn = 0. 

5778 

5779 if request.has_statics: 

5780 tcumu_static = num.sum(num.vstack(tcounters_static_list), axis=0) 

5781 t_static = num.sum(tcumu_static) 

5782 perc_static = map(float, tcumu_static / t_static * 100.) 

5783 (s.t_perc_static_get_store, 

5784 s.t_perc_static_discretize_basesource, 

5785 s.t_perc_static_sum_statics, 

5786 s.t_perc_static_post_process) = perc_static 

5787 

5788 s.t_wallclock = tt1 - tt0 

5789 if resource: 

5790 s.t_cpu = ( 

5791 (rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) - 

5792 (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime)) 

5793 s.n_read_blocks = ( 

5794 (rs1.ru_inblock + rc1.ru_inblock) - 

5795 (rs0.ru_inblock + rc0.ru_inblock)) 

5796 

5797 n_records_stacked = 0. 

5798 for results in results_list: 

5799 for result in results: 

5800 if not isinstance(result, meta.Result): 

5801 continue 

5802 shr = float(result.n_shared_stacking) 

5803 n_records_stacked += result.n_records_stacked / shr 

5804 s.t_perc_optimize += result.t_optimize / shr 

5805 s.t_perc_stack += result.t_stack / shr 

5806 s.n_records_stacked = int(n_records_stacked) 

5807 if t_dyn != 0.: 

5808 s.t_perc_optimize /= t_dyn * 100 

5809 s.t_perc_stack /= t_dyn * 100 

5810 

5811 return Response( 

5812 request=request, 

5813 results_list=results_list, 

5814 stats=s) 

5815 

5816 

5817class RemoteEngine(Engine): 

5818 ''' 

5819 Client for remote synthetic seismogram calculator. 

5820 ''' 

5821 

5822 site = String.T(default=ws.g_default_site, optional=True) 

5823 url = String.T(default=ws.g_url, optional=True) 

5824 

5825 def process(self, request=None, status_callback=None, **kwargs): 

5826 

5827 if request is None: 

5828 request = Request(**kwargs) 

5829 

5830 return ws.seismosizer(url=self.url, site=self.site, request=request) 

5831 

5832 

5833g_engine = None 

5834 

5835 

5836def get_engine(store_superdirs=[]): 

5837 global g_engine 

5838 if g_engine is None: 

5839 g_engine = LocalEngine(use_env=True, use_config=True) 

5840 

5841 for d in store_superdirs: 

5842 if d not in g_engine.store_superdirs: 

5843 g_engine.store_superdirs.append(d) 

5844 

5845 return g_engine 

5846 

5847 

5848class SourceGroup(Object): 

5849 

5850 def __getattr__(self, k): 

5851 return num.fromiter((getattr(s, k) for s in self), 

5852 dtype=float) 

5853 

5854 def __iter__(self): 

5855 raise NotImplementedError( 

5856 'This method should be implemented in subclass.') 

5857 

5858 def __len__(self): 

5859 raise NotImplementedError( 

5860 'This method should be implemented in subclass.') 

5861 

5862 

5863class SourceList(SourceGroup): 

5864 sources = List.T(Source.T()) 

5865 

5866 def append(self, s): 

5867 self.sources.append(s) 

5868 

5869 def __iter__(self): 

5870 return iter(self.sources) 

5871 

5872 def __len__(self): 

5873 return len(self.sources) 

5874 

5875 

5876class SourceGrid(SourceGroup): 

5877 

5878 base = Source.T() 

5879 variables = Dict.T(String.T(), Range.T()) 

5880 order = List.T(String.T()) 

5881 

5882 def __len__(self): 

5883 n = 1 

5884 for (k, v) in self.make_coords(self.base): 

5885 n *= len(list(v)) 

5886 

5887 return n 

5888 

5889 def __iter__(self): 

5890 for items in permudef(self.make_coords(self.base)): 

5891 s = self.base.clone(**{k: v for (k, v) in items}) 

5892 s.regularize() 

5893 yield s 

5894 

5895 def ordered_params(self): 

5896 ks = list(self.variables.keys()) 

5897 for k in self.order + list(self.base.keys()): 

5898 if k in ks: 

5899 yield k 

5900 ks.remove(k) 

5901 if ks: 

5902 raise Exception('Invalid parameter "%s" for source type "%s".' % 

5903 (ks[0], self.base.__class__.__name__)) 

5904 

5905 def make_coords(self, base): 

5906 return [(param, self.variables[param].make(base=base[param])) 

5907 for param in self.ordered_params()] 

5908 

5909 

5910source_classes = [ 

5911 Source, 

5912 SourceWithMagnitude, 

5913 SourceWithDerivedMagnitude, 

5914 ExplosionSource, 

5915 RectangularExplosionSource, 

5916 DCSource, 

5917 CLVDSource, 

5918 VLVDSource, 

5919 MTSource, 

5920 RectangularSource, 

5921 PseudoDynamicRupture, 

5922 DoubleDCSource, 

5923 RingfaultSource, 

5924 CombiSource, 

5925 SFSource, 

5926 PorePressurePointSource, 

5927 PorePressureLineSource, 

5928] 

5929 

5930stf_classes = [ 

5931 STF, 

5932 BoxcarSTF, 

5933 TriangularSTF, 

5934 HalfSinusoidSTF, 

5935 ResonatorSTF, 

5936 TremorSTF, 

5937] 

5938 

5939__all__ = ''' 

5940Cloneable 

5941NoDefaultStoreSet 

5942SeismosizerError 

5943BadRequest 

5944NoSuchStore 

5945DerivedMagnitudeError 

5946STFMode 

5947'''.split() + [S.__name__ for S in source_classes + stf_classes] + ''' 

5948Request 

5949ProcessingStats 

5950Response 

5951Engine 

5952LocalEngine 

5953RemoteEngine 

5954source_classes 

5955get_engine 

5956Range 

5957SourceGroup 

5958SourceList 

5959SourceGrid 

5960map_anchor 

5961'''.split()