Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/gf/seismosizer.py: 83%
2581 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7High level synthetic seismogram synthesis.
9.. _coordinate-system-names:
11Coordinate systems
12..................
14Coordinate system names commonly used in source models.
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'''
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
40import numpy as num
41from scipy.interpolate import RegularGridInterpolator
43from pyrocko.guts import (Object, Float, String, StringChoice, List,
44 Timestamp, Int, SObject, ArgumentError, Dict,
45 ValidationError, Bool)
46from pyrocko.guts_array import Array
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
55from . import meta, store, ws
56from .tractions import TractionField, DirectedTractions
57from .targets import Target, StaticTarget, SatelliteTarget
59pjoin = os.path.join
61guts_prefix = 'pf'
63d2r = math.pi / 180.
64r2d = 180. / math.pi
65km = 1e3
67logger = logging.getLogger('pyrocko.gf.seismosizer')
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
77 return 0
79 anone = a is None
80 bnone = b is None
82 if anone and bnone:
83 return 0
85 if anone:
86 return -1
88 if bnone:
89 return 1
91 return bool(a > b) - bool(a < b)
94def xtime():
95 return time.time()
98class SeismosizerError(Exception):
99 pass
102class BadRequest(SeismosizerError):
103 pass
106class DuplicateStoreId(Exception):
107 pass
110class NoDefaultStoreSet(Exception):
111 '''
112 Raised, when a default store would be used but none is set.
113 '''
114 pass
117class ConversionError(Exception):
118 pass
121class STFError(SeismosizerError):
122 pass
125class NoSuchStore(BadRequest):
127 def __init__(self, store_id, store_dirs, store_superdirs, store_ids_avail):
128 BadRequest.__init__(self)
129 self.store_id = store_id
130 self.store_dirs = store_dirs
131 self.store_superdirs = store_superdirs
132 self.store_ids_avail = store_ids_avail
134 def __str__(self):
135 lines = ['No GF store with id "%s" found.' % self.store_id]
137 def add_entries(lines, name, entries):
138 lines.append(' %s:' % name)
139 if not entries:
140 lines.append(' *empty*')
141 else:
142 for entry in entries:
143 lines.append(' - %s' % entry)
145 add_entries(lines, 'store_superdirs searched', self.store_superdirs)
146 add_entries(lines, 'store_dirs searched', self.store_dirs)
147 add_entries(lines, 'store IDs available', self.store_ids_avail)
148 return '\n'.join(lines)
151def ufloat(s):
152 units = {
153 'k': 1e3,
154 'M': 1e6,
155 }
157 factor = 1.0
158 if s and s[-1] in units:
159 factor = units[s[-1]]
160 s = s[:-1]
161 if not s:
162 raise ValueError("unit without a number: '%s'" % s)
164 return float(s) * factor
167def ufloat_or_none(s):
168 if s:
169 return ufloat(s)
170 else:
171 return None
174def int_or_none(s):
175 if s:
176 return int(s)
177 else:
178 return None
181def nonzero(x, eps=1e-15):
182 return abs(x) > eps
185def permudef(ln, j=0):
186 if j < len(ln):
187 k, v = ln[j]
188 for y in v:
189 ln[j] = k, y
190 for s in permudef(ln, j + 1):
191 yield s
193 ln[j] = k, v
194 return
195 else:
196 yield ln
199def arr(x):
200 return num.atleast_1d(num.asarray(x))
203def discretize_rect_source(deltas, deltat, time, north, east, depth,
204 strike, dip, length, width,
205 anchor, velocity=None, stf=None,
206 nucleation_x=None, nucleation_y=None,
207 decimation_factor=1, pointsonly=False,
208 plane_coords=False,
209 aggressive_oversampling=False):
211 if stf is None:
212 stf = STF()
214 if not velocity and not pointsonly:
215 raise AttributeError('velocity is required in time mode')
217 mindeltagf = float(num.min(deltas))
218 if velocity:
219 mindeltagf = min(mindeltagf, deltat * velocity)
221 ln = length
222 wd = width
224 if aggressive_oversampling:
225 nl = int((2. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
226 nw = int((2. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
227 else:
228 nl = int((1. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
229 nw = int((1. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
231 n = int(nl * nw)
233 dl = ln / nl
234 dw = wd / nw
236 xl = num.linspace(-0.5 * (ln - dl), 0.5 * (ln - dl), nl)
237 xw = num.linspace(-0.5 * (wd - dw), 0.5 * (wd - dw), nw)
239 points = num.zeros((n, 3))
240 points[:, 0] = num.tile(xl, nw)
241 points[:, 1] = num.repeat(xw, nl)
243 if nucleation_x is not None:
244 dist_x = num.abs(nucleation_x - points[:, 0])
245 else:
246 dist_x = num.zeros(n)
248 if nucleation_y is not None:
249 dist_y = num.abs(nucleation_y - points[:, 1])
250 else:
251 dist_y = num.zeros(n)
253 dist = num.sqrt(dist_x**2 + dist_y**2)
254 times = dist / velocity
256 anch_x, anch_y = map_anchor[anchor]
258 points[:, 0] -= anch_x * 0.5 * length
259 points[:, 1] -= anch_y * 0.5 * width
261 if plane_coords:
262 return points, dl, dw, nl, nw
264 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
265 points = num.dot(rotmat.T, points.T).T
267 points[:, 0] += north
268 points[:, 1] += east
269 points[:, 2] += depth
271 if pointsonly:
272 return points, dl, dw, nl, nw
274 xtau, amplitudes = stf.discretize_t(deltat, time)
275 nt = xtau.size
277 points2 = num.repeat(points, nt, axis=0)
278 times2 = (times[:, num.newaxis] + xtau[num.newaxis, :]).ravel()
279 amplitudes2 = num.tile(amplitudes, n)
281 return points2, times2, amplitudes2, dl, dw, nl, nw
284def check_rect_source_discretisation(points2, nl, nw, store):
285 # We assume a non-rotated fault plane
286 N_CRITICAL = 8
287 points = points2.T.reshape((3, nl, nw))
288 if points.size <= N_CRITICAL:
289 logger.warning('RectangularSource is defined by only %d sub-sources!'
290 % points.size)
291 return True
293 distances = num.sqrt(
294 (points[0, 0, :] - points[0, 1, :])**2 +
295 (points[1, 0, :] - points[1, 1, :])**2 +
296 (points[2, 0, :] - points[2, 1, :])**2)
298 depths = points[2, 0, :]
299 vs_profile = store.config.get_vs(
300 lat=0., lon=0.,
301 points=num.repeat(depths[:, num.newaxis], 3, axis=1),
302 interpolation='multilinear')
304 min_wavelength = vs_profile * (store.config.deltat * 2)
305 if not num.all(min_wavelength > distances / 2):
306 return False
307 return True
310def outline_rect_source(strike, dip, length, width, anchor):
311 ln = length
312 wd = width
313 points = num.array(
314 [[-0.5 * ln, -0.5 * wd, 0.],
315 [0.5 * ln, -0.5 * wd, 0.],
316 [0.5 * ln, 0.5 * wd, 0.],
317 [-0.5 * ln, 0.5 * wd, 0.],
318 [-0.5 * ln, -0.5 * wd, 0.]])
320 anch_x, anch_y = map_anchor[anchor]
321 points[:, 0] -= anch_x * 0.5 * length
322 points[:, 1] -= anch_y * 0.5 * width
324 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
326 return num.dot(rotmat.T, points.T).T
329def from_plane_coords(
330 strike, dip, length, width, depth, x_plane_coords, y_plane_coords,
331 lat=0., lon=0.,
332 north_shift=0, east_shift=0,
333 anchor='top', cs='xy'):
335 ln = length
336 wd = width
337 x_abs = []
338 y_abs = []
339 if not isinstance(x_plane_coords, list):
340 x_plane_coords = [x_plane_coords]
341 y_plane_coords = [y_plane_coords]
343 for x_plane, y_plane in zip(x_plane_coords, y_plane_coords):
344 points = num.array(
345 [[-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
346 [0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
347 [0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
348 [-0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
349 [-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.]])
351 anch_x, anch_y = map_anchor[anchor]
352 points[:, 0] -= anch_x * 0.5 * length
353 points[:, 1] -= anch_y * 0.5 * width
355 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
357 points = num.dot(rotmat.T, points.T).T
358 points[:, 0] += north_shift
359 points[:, 1] += east_shift
360 points[:, 2] += depth
361 if cs in ('latlon', 'lonlat'):
362 latlon = ne_to_latlon(lat, lon,
363 points[:, 0], points[:, 1])
364 latlon = num.array(latlon).T
365 x_abs.append(latlon[1:2, 1])
366 y_abs.append(latlon[2:3, 0])
367 if cs == 'xy':
368 x_abs.append(points[1:2, 1])
369 y_abs.append(points[2:3, 0])
371 if cs == 'lonlat':
372 return y_abs, x_abs
373 else:
374 return x_abs, y_abs
377def points_on_rect_source(
378 strike, dip, length, width, anchor,
379 discretized_basesource=None, points_x=None, points_y=None):
381 ln = length
382 wd = width
384 if isinstance(points_x, list) or isinstance(points_x, float):
385 points_x = num.array([points_x])
386 if isinstance(points_y, list) or isinstance(points_y, float):
387 points_y = num.array([points_y])
389 if discretized_basesource:
390 ds = discretized_basesource
392 nl_patches = ds.nl + 1
393 nw_patches = ds.nw + 1
395 npoints = nl_patches * nw_patches
396 points = num.zeros((npoints, 3))
397 ln_patches = num.array([il for il in range(nl_patches)])
398 wd_patches = num.array([iw for iw in range(nw_patches)])
400 points_ln =\
401 2 * ((ln_patches - num.min(ln_patches)) / num.ptp(ln_patches)) - 1
402 points_wd =\
403 2 * ((wd_patches - num.min(wd_patches)) / num.ptp(wd_patches)) - 1
405 for il in range(nl_patches):
406 for iw in range(nw_patches):
407 points[il * nw_patches + iw, :] = num.array([
408 points_ln[il] * ln * 0.5,
409 points_wd[iw] * wd * 0.5, 0.0])
411 elif points_x.shape[0] > 0 and points_y.shape[0] > 0:
412 points = num.zeros(shape=((len(points_x), 3)))
413 for i, (x, y) in enumerate(zip(points_x, points_y)):
414 points[i, :] = num.array(
415 [x * 0.5 * ln, y * 0.5 * wd, 0.0])
417 anch_x, anch_y = map_anchor[anchor]
419 points[:, 0] -= anch_x * 0.5 * ln
420 points[:, 1] -= anch_y * 0.5 * wd
422 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
424 return num.dot(rotmat.T, points.T).T
427class InvalidGridDef(Exception):
428 pass
431class Range(SObject):
432 '''
433 Convenient range specification.
435 Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]::
437 Range('0 .. 10k : 1k')
438 Range(start=0., stop=10e3, step=1e3)
439 Range(0, 10e3, 1e3)
440 Range('0 .. 10k @ 11')
441 Range(start=0., stop=10*km, n=11)
443 Range(0, 10e3, n=11)
444 Range(values=[x*1e3 for x in range(11)])
446 Depending on the use context, it can be possible to omit any part of the
447 specification. E.g. in the context of extracting a subset of an already
448 existing range, the existing range's specification values would be filled
449 in where missing.
451 The values are distributed with equal spacing, unless the ``spacing``
452 argument is modified. The values can be created offset or relative to an
453 external base value with the ``relative`` argument if the use context
454 supports this.
456 The range specification can be expressed with a short string
457 representation::
459 'start .. stop @ num | spacing, relative'
460 'start .. stop : step | spacing, relative'
462 most parts of the expression can be omitted if not needed. Whitespace is
463 allowed for readability but can also be omitted.
464 '''
466 start = Float.T(optional=True)
467 stop = Float.T(optional=True)
468 step = Float.T(optional=True)
469 n = Int.T(optional=True)
470 values = Array.T(optional=True, dtype=float, shape=(None,))
472 spacing = StringChoice.T(
473 choices=['lin', 'log', 'symlog'],
474 default='lin',
475 optional=True)
477 relative = StringChoice.T(
478 choices=['', 'add', 'mult'],
479 default='',
480 optional=True)
482 pattern = re.compile(r'^((?P<start>.*)\.\.(?P<stop>[^@|:]*))?'
483 r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?'
484 r'(\|(?P<stuff>.+))?$')
486 def __init__(self, *args, **kwargs):
487 d = {}
488 if len(args) == 1:
489 d = self.parse(args[0])
490 elif len(args) in (2, 3):
491 d['start'], d['stop'] = [float(x) for x in args[:2]]
492 if len(args) == 3:
493 d['step'] = float(args[2])
495 for k, v in kwargs.items():
496 if k in d:
497 raise ArgumentError('%s specified more than once' % k)
499 d[k] = v
501 SObject.__init__(self, **d)
503 def __str__(self):
504 def sfloat(x):
505 if x is not None:
506 return '%g' % x
507 else:
508 return ''
510 if self.values:
511 return ','.join('%g' % x for x in self.values)
513 if self.start is None and self.stop is None:
514 s0 = ''
515 else:
516 s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop))
518 s1 = ''
519 if self.step is not None:
520 s1 = [' : %g', ':%g'][s0 == ''] % self.step
521 elif self.n is not None:
522 s1 = [' @ %i', '@%i'][s0 == ''] % self.n
524 if self.spacing == 'lin' and self.relative == '':
525 s2 = ''
526 else:
527 x = []
528 if self.spacing != 'lin':
529 x.append(self.spacing)
531 if self.relative != '':
532 x.append(self.relative)
534 s2 = ' | %s' % ','.join(x)
536 return s0 + s1 + s2
538 @classmethod
539 def parse(cls, s):
540 s = re.sub(r'\s+', '', s)
541 m = cls.pattern.match(s)
542 if not m:
543 try:
544 vals = [ufloat(x) for x in s.split(',')]
545 except Exception:
546 raise InvalidGridDef(
547 '"%s" is not a valid range specification' % s)
549 return dict(values=num.array(vals, dtype=float))
551 d = m.groupdict()
552 try:
553 start = ufloat_or_none(d['start'])
554 stop = ufloat_or_none(d['stop'])
555 step = ufloat_or_none(d['step'])
556 n = int_or_none(d['n'])
557 except Exception:
558 raise InvalidGridDef(
559 '"%s" is not a valid range specification' % s)
561 spacing = 'lin'
562 relative = ''
564 if d['stuff'] is not None:
565 t = d['stuff'].split(',')
566 for x in t:
567 if x in cls.spacing.choices:
568 spacing = x
569 elif x and x in cls.relative.choices:
570 relative = x
571 else:
572 raise InvalidGridDef(
573 '"%s" is not a valid range specification' % s)
575 return dict(start=start, stop=stop, step=step, n=n, spacing=spacing,
576 relative=relative)
578 def make(self, mi=None, ma=None, inc=None, base=None, eps=1e-5):
579 if self.values:
580 return self.values
582 start = self.start
583 stop = self.stop
584 step = self.step
585 n = self.n
587 swap = step is not None and step < 0.
588 if start is None:
589 start = [mi, ma][swap]
590 if stop is None:
591 stop = [ma, mi][swap]
592 if step is None and inc is not None:
593 step = [inc, -inc][ma < mi]
595 if start is None or stop is None:
596 raise InvalidGridDef(
597 'Cannot use range specification "%s" without start '
598 'and stop in this context' % self)
600 if step is None and n is None:
601 step = stop - start
603 if n is None:
604 if (step < 0) != (stop - start < 0):
605 raise InvalidGridDef(
606 'Range specification "%s" has inconsistent ordering '
607 '(step < 0 => stop > start)' % self)
609 n = int(round((stop - start) / step)) + 1
610 stop2 = start + (n - 1) * step
611 if abs(stop - stop2) > eps:
612 n = int(math.floor((stop - start) / step)) + 1
613 stop = start + (n - 1) * step
614 else:
615 stop = stop2
617 if start == stop:
618 n = 1
620 if self.spacing == 'lin':
621 vals = num.linspace(start, stop, n)
623 elif self.spacing in ('log', 'symlog'):
624 if start > 0. and stop > 0.:
625 vals = num.exp(num.linspace(num.log(start),
626 num.log(stop), n))
627 elif start < 0. and stop < 0.:
628 vals = -num.exp(num.linspace(num.log(-start),
629 num.log(-stop), n))
630 else:
631 raise InvalidGridDef(
632 'Log ranges should not include or cross zero '
633 '(in range specification "%s").' % self)
635 if self.spacing == 'symlog':
636 nvals = - vals
637 vals = num.concatenate((nvals[::-1], vals))
639 if self.relative in ('add', 'mult') and base is None:
640 raise InvalidGridDef(
641 'Cannot use relative range specification in this context.')
643 vals = self.make_relative(base, vals)
645 return list(map(float, vals))
647 def make_relative(self, base, vals):
648 if self.relative == 'add':
649 vals += base
651 if self.relative == 'mult':
652 vals *= base
654 return vals
657class GridDefElement(Object):
659 param = meta.StringID.T()
660 rs = Range.T()
662 def __init__(self, shorthand=None, **kwargs):
663 if shorthand is not None:
664 t = shorthand.split('=')
665 if len(t) != 2:
666 raise InvalidGridDef(
667 'Invalid grid specification element: %s' % shorthand)
669 sp, sr = t[0].strip(), t[1].strip()
671 kwargs['param'] = sp
672 kwargs['rs'] = Range(sr)
674 Object.__init__(self, **kwargs)
676 def shorthand(self):
677 return self.param + ' = ' + str(self.rs)
680class GridDef(Object):
682 elements = List.T(GridDefElement.T())
684 def __init__(self, shorthand=None, **kwargs):
685 if shorthand is not None:
686 t = shorthand.splitlines()
687 tt = []
688 for x in t:
689 x = x.strip()
690 if x:
691 tt.extend(x.split(';'))
693 elements = []
694 for se in tt:
695 elements.append(GridDef(se))
697 kwargs['elements'] = elements
699 Object.__init__(self, **kwargs)
701 def shorthand(self):
702 return '; '.join(str(x) for x in self.elements)
705class Cloneable(object):
706 '''
707 Mix-in class for Guts objects, providing dict-like access and cloning.
708 '''
710 def __iter__(self):
711 return iter(self.T.propnames)
713 def __getitem__(self, k):
714 if k not in self.keys():
715 raise KeyError(k)
717 return getattr(self, k)
719 def __setitem__(self, k, v):
720 if k not in self.keys():
721 raise KeyError(k)
723 return setattr(self, k, v)
725 def clone(self, **kwargs):
726 '''
727 Make a copy of the object.
729 A new object of the same class is created and initialized with the
730 parameters of the object on which this method is called on. If
731 ``kwargs`` are given, these are used to override any of the
732 initialization parameters.
733 '''
735 d = dict(self)
736 for k in d:
737 v = d[k]
738 if isinstance(v, Cloneable):
739 d[k] = v.clone()
741 d.update(kwargs)
742 return self.__class__(**d)
744 @classmethod
745 def keys(cls):
746 '''
747 Get list of the source model's parameter names.
748 '''
750 return cls.T.propnames
753class STF(Object, Cloneable):
755 '''
756 Base class for source time functions.
757 '''
759 def __init__(self, effective_duration=None, **kwargs):
760 if effective_duration is not None:
761 kwargs['duration'] = effective_duration / \
762 self.factor_duration_to_effective()
764 Object.__init__(self, **kwargs)
766 @classmethod
767 def factor_duration_to_effective(cls):
768 return 1.0
770 def centroid_time(self, tref):
771 return tref
773 @property
774 def effective_duration(self):
775 return self.duration * self.factor_duration_to_effective()
777 def discretize_t(self, deltat, tref):
778 tl = math.floor(tref / deltat) * deltat
779 th = math.ceil(tref / deltat) * deltat
780 if tl == th:
781 return num.array([tl], dtype=float), num.ones(1)
782 else:
783 return (
784 num.array([tl, th], dtype=float),
785 num.array([th - tref, tref - tl], dtype=float) / deltat)
787 def base_key(self):
788 return (type(self).__name__,)
791g_unit_pulse = STF()
794def sshift(times, amplitudes, tshift, deltat):
796 t0 = math.floor(tshift / deltat) * deltat
797 t1 = math.ceil(tshift / deltat) * deltat
798 if t0 == t1:
799 return times, amplitudes
801 amplitudes2 = num.zeros(amplitudes.size + 1, dtype=float)
803 amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes
804 amplitudes2[1:] += (tshift - t0) / deltat * amplitudes
806 times2 = num.arange(times.size + 1, dtype=float) * \
807 deltat + times[0] + t0
809 return times2, amplitudes2
812class BoxcarSTF(STF):
814 '''
815 Boxcar type source time function.
817 .. figure :: /static/stf-BoxcarSTF.svg
818 :width: 40%
819 :align: center
820 :alt: boxcar source time function
821 '''
823 duration = Float.T(
824 default=0.0,
825 help='duration of the boxcar')
827 anchor = Float.T(
828 default=0.0,
829 help='anchor point with respect to source.time: ('
830 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
831 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
832 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
834 @classmethod
835 def factor_duration_to_effective(cls):
836 return 1.0
838 def centroid_time(self, tref):
839 return tref - 0.5 * self.duration * self.anchor
841 def discretize_t(self, deltat, tref):
842 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
843 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
844 tmin = round(tmin_stf / deltat) * deltat
845 tmax = round(tmax_stf / deltat) * deltat
846 nt = int(round((tmax - tmin) / deltat)) + 1
847 times = num.linspace(tmin, tmax, nt)
848 amplitudes = num.ones_like(times)
849 if times.size > 1:
850 t_edges = num.linspace(
851 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
852 t = tmin_stf + self.duration * num.array(
853 [0.0, 0.0, 1.0, 1.0], dtype=float)
854 f = num.array([0., 1., 1., 0.], dtype=float)
855 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
856 amplitudes /= num.sum(amplitudes)
858 tshift = (num.sum(amplitudes * times) - self.centroid_time(tref))
860 return sshift(times, amplitudes, -tshift, deltat)
862 def base_key(self):
863 return (type(self).__name__, self.duration, self.anchor)
866class TriangularSTF(STF):
868 '''
869 Triangular type source time function.
871 .. figure :: /static/stf-TriangularSTF.svg
872 :width: 40%
873 :align: center
874 :alt: triangular source time function
875 '''
877 duration = Float.T(
878 default=0.0,
879 help='baseline of the triangle')
881 peak_ratio = Float.T(
882 default=0.5,
883 help='fraction of time compared to duration, '
884 'when the maximum amplitude is reached')
886 anchor = Float.T(
887 default=0.0,
888 help='anchor point with respect to source.time: ('
889 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
890 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
891 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
893 @classmethod
894 def factor_duration_to_effective(cls, peak_ratio=None):
895 if peak_ratio is None:
896 peak_ratio = cls.peak_ratio.default()
898 return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0)
900 def __init__(self, effective_duration=None, **kwargs):
901 if effective_duration is not None:
902 kwargs['duration'] = effective_duration / \
903 self.factor_duration_to_effective(
904 kwargs.get('peak_ratio', None))
906 STF.__init__(self, **kwargs)
908 @property
909 def centroid_ratio(self):
910 ra = self.peak_ratio
911 rb = 1.0 - ra
912 return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb)
914 def centroid_time(self, tref):
915 ca = self.centroid_ratio
916 cb = 1.0 - ca
917 if self.anchor <= 0.:
918 return tref - ca * self.duration * self.anchor
919 else:
920 return tref - cb * self.duration * self.anchor
922 @property
923 def effective_duration(self):
924 return self.duration * self.factor_duration_to_effective(
925 self.peak_ratio)
927 def tminmax_stf(self, tref):
928 ca = self.centroid_ratio
929 cb = 1.0 - ca
930 if self.anchor <= 0.:
931 tmin_stf = tref - ca * self.duration * (self.anchor + 1.)
932 tmax_stf = tmin_stf + self.duration
933 else:
934 tmax_stf = tref + cb * self.duration * (1. - self.anchor)
935 tmin_stf = tmax_stf - self.duration
937 return tmin_stf, tmax_stf
939 def discretize_t(self, deltat, tref):
940 tmin_stf, tmax_stf = self.tminmax_stf(tref)
942 tmin = round(tmin_stf / deltat) * deltat
943 tmax = round(tmax_stf / deltat) * deltat
944 nt = int(round((tmax - tmin) / deltat)) + 1
945 if nt > 1:
946 t_edges = num.linspace(
947 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
948 t = tmin_stf + self.duration * num.array(
949 [0.0, self.peak_ratio, 1.0], dtype=float)
950 f = num.array([0., 1., 0.], dtype=float)
951 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
952 amplitudes /= num.sum(amplitudes)
953 else:
954 amplitudes = num.ones(1)
956 times = num.linspace(tmin, tmax, nt)
957 return times, amplitudes
959 def base_key(self):
960 return (
961 type(self).__name__, self.duration, self.peak_ratio, self.anchor)
964class HalfSinusoidSTF(STF):
966 '''
967 Half sinusoid type source time function.
969 .. figure :: /static/stf-HalfSinusoidSTF.svg
970 :width: 40%
971 :align: center
972 :alt: half-sinusouid source time function
973 '''
975 duration = Float.T(
976 default=0.0,
977 help='duration of the half-sinusoid (baseline)')
979 anchor = Float.T(
980 default=0.0,
981 help='anchor point with respect to source.time: ('
982 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
983 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
984 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
986 exponent = Int.T(
987 default=1,
988 help='set to 2 to use square of the half-period sinusoidal function.')
990 def __init__(self, effective_duration=None, **kwargs):
991 if effective_duration is not None:
992 kwargs['duration'] = effective_duration / \
993 self.factor_duration_to_effective(
994 kwargs.get('exponent', 1))
996 STF.__init__(self, **kwargs)
998 @classmethod
999 def factor_duration_to_effective(cls, exponent):
1000 if exponent == 1:
1001 return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi
1002 elif exponent == 2:
1003 return math.sqrt(math.pi**2 - 6) / math.pi
1004 else:
1005 raise ValueError('Exponent for HalfSinusoidSTF must be 1 or 2.')
1007 @property
1008 def effective_duration(self):
1009 return self.duration * self.factor_duration_to_effective(self.exponent)
1011 def centroid_time(self, tref):
1012 return tref - 0.5 * self.duration * self.anchor
1014 def discretize_t(self, deltat, tref):
1015 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1016 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1017 tmin = round(tmin_stf / deltat) * deltat
1018 tmax = round(tmax_stf / deltat) * deltat
1019 nt = int(round((tmax - tmin) / deltat)) + 1
1020 if nt > 1:
1021 t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace(
1022 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)))
1024 if self.exponent == 1:
1025 fint = -num.cos(
1026 (t_edges - tmin_stf) * (math.pi / self.duration))
1028 elif self.exponent == 2:
1029 fint = (t_edges - tmin_stf) / self.duration \
1030 - 1.0 / (2.0 * math.pi) * num.sin(
1031 (t_edges - tmin_stf) * (2.0 * math.pi / self.duration))
1032 else:
1033 raise ValueError(
1034 'Exponent for HalfSinusoidSTF must be 1 or 2.')
1036 amplitudes = fint[1:] - fint[:-1]
1037 amplitudes /= num.sum(amplitudes)
1038 else:
1039 amplitudes = num.ones(1)
1041 times = num.linspace(tmin, tmax, nt)
1042 return times, amplitudes
1044 def base_key(self):
1045 return (type(self).__name__, self.duration, self.anchor)
1048class SmoothRampSTF(STF):
1049 '''
1050 Smooth-ramp type source time function for near-field displacement.
1051 Based on moment function of double-couple point source proposed by [1]_.
1053 .. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow
1054 earthquakes from Love-wave modelling for regional distances, PEPI 32,
1055 312-324.
1057 .. figure :: /static/stf-SmoothRampSTF.svg
1058 :width: 40%
1059 :alt: smooth ramp source time function
1060 '''
1061 duration = Float.T(
1062 default=0.0,
1063 help='duration of the ramp (baseline)')
1065 rise_ratio = Float.T(
1066 default=0.5,
1067 help='fraction of time compared to duration, '
1068 'when the maximum amplitude is reached')
1070 anchor = Float.T(
1071 default=0.0,
1072 help='anchor point with respect to source.time: ('
1073 '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, '
1074 '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, '
1075 '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)')
1077 def discretize_t(self, deltat, tref):
1078 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1079 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1080 tmin = round(tmin_stf / deltat) * deltat
1081 tmax = round(tmax_stf / deltat) * deltat
1082 D = round((tmax - tmin) / deltat) * deltat
1083 nt = int(round(D / deltat)) + 1
1084 times = num.linspace(tmin, tmax, nt)
1085 if nt > 1:
1086 rise_time = self.rise_ratio * self.duration
1087 amplitudes = num.ones_like(times)
1088 tp = tmin + rise_time
1089 ii = num.where(times <= tp)
1090 t_inc = times[ii]
1091 a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time)
1092 b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0
1093 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b)
1095 amplitudes /= num.sum(amplitudes)
1096 else:
1097 amplitudes = num.ones(1)
1099 return times, amplitudes
1101 def base_key(self):
1102 return (type(self).__name__,
1103 self.duration, self.rise_ratio, self.anchor)
1106class ResonatorSTF(STF):
1107 '''
1108 Simple resonator like source time function.
1110 .. math ::
1112 f(t) = 0 for t < 0
1113 f(t) = e^{-t/tau} * sin(2 * pi * f * t)
1116 .. figure :: /static/stf-SmoothRampSTF.svg
1117 :width: 40%
1118 :alt: smooth ramp source time function
1120 '''
1122 duration = Float.T(
1123 default=0.0,
1124 help='decay time')
1126 frequency = Float.T(
1127 default=1.0,
1128 help='resonance frequency')
1130 def discretize_t(self, deltat, tref):
1131 tmin_stf = tref
1132 tmax_stf = tref + self.duration * 3
1133 tmin = math.floor(tmin_stf / deltat) * deltat
1134 tmax = math.ceil(tmax_stf / deltat) * deltat
1135 times = util.arange2(tmin, tmax, deltat)
1136 amplitudes = num.exp(-(times - tref) / self.duration) \
1137 * num.sin(2.0 * num.pi * self.frequency * (times - tref))
1139 return times, amplitudes
1141 def base_key(self):
1142 return (type(self).__name__,
1143 self.duration, self.frequency)
1146class TremorSTF(STF):
1147 '''
1148 Oszillating source time function.
1150 .. math ::
1152 f(t) = 0 for t < -tau/2 or t > tau/2
1153 f(t) = cos(pi/tau*t) * sin(2 * pi * f * t)
1155 '''
1157 duration = Float.T(
1158 default=0.0,
1159 help='Tremor duration [s]')
1161 frequency = Float.T(
1162 default=1.0,
1163 help='Frequency [Hz]')
1165 def discretize_t(self, deltat, tref):
1166 tmin_stf = tref - 0.5 * self.duration
1167 tmax_stf = tref + 0.5 * self.duration
1168 tmin = math.floor(tmin_stf / deltat) * deltat
1169 tmax = math.ceil(tmax_stf / deltat) * deltat
1170 times = util.arange2(tmin, tmax, deltat)
1171 mask = num.logical_and(
1172 tref - 0.5 * self.duration < times,
1173 times < tref + 0.5 * self.duration)
1175 amplitudes = num.zeros_like(times)
1176 amplitudes[mask] = num.cos(num.pi/self.duration*(times[mask] - tref)) \
1177 * num.sin(2.0 * num.pi * self.frequency * (times[mask] - tref))
1178 amplitudes[mask] *= deltat
1180 return times, amplitudes
1182 def base_key(self):
1183 return (type(self).__name__,
1184 self.duration, self.frequency)
1187class SimpleLandslideSTF(STF):
1189 '''
1190 Doublepulse land-slide STF which respects conservation of momentum.
1191 '''
1193 duration_acceleration = Float.T(
1194 default=1.0,
1195 help='Duratian of the acceleration phase [s].')
1197 duration_deceleration = Float.T(
1198 default=1.0,
1199 help='Duration of the deceleration phase [s].')
1201 mute_acceleration = Bool.T(
1202 default=False,
1203 help='set acceleration to zero (for testing)')
1205 mute_deceleration = Bool.T(
1206 default=False,
1207 help='set acceleration to zero (for testing)')
1209 def discretize_t(self, deltat, tref):
1211 d_acc = self.duration_acceleration
1212 d_dec = self.duration_deceleration
1214 tmin_stf = tref
1215 tmax_stf = tref + d_acc + d_dec
1216 tmin = math.floor(tmin_stf / deltat) * deltat
1217 tmax = math.ceil(tmax_stf / deltat) * deltat
1218 times = util.arange2(tmin, tmax, deltat, epsilon=1e-3)
1220 mask_acc = num.logical_and(
1221 tref <= times,
1222 times < tref + d_acc)
1224 mask_dec = num.logical_and(
1225 tref + d_acc <= times,
1226 times < tref + d_acc + d_dec)
1228 n_acc = num.sum(mask_acc)
1229 if n_acc < 1:
1230 raise STFError(
1231 'SimpleLandslideSTF: `duration_acceleration` must be longer '
1232 'than sampling interval.')
1234 n_dec = num.sum(mask_dec)
1235 if n_dec < 1:
1236 raise STFError(
1237 'SimpleLandslideSTF: `duration_deceleration` must be longer '
1238 'than sampling interval.')
1240 amplitudes = num.zeros_like(times)
1242 amplitudes[mask_acc] = - num.sin(
1243 (times[mask_acc] - tref) / d_acc * num.pi)
1245 amplitudes[mask_dec] = num.sin(
1246 (times[mask_dec] - (tref + d_acc)) / d_dec * num.pi)
1248 sum_acc = num.abs(num.sum(amplitudes[mask_acc]))
1249 if sum_acc != 0.0:
1250 amplitudes[mask_acc] /= sum_acc
1251 else:
1252 amplitudes[mask_acc] = 1.0 / n_acc
1254 sum_dec = num.abs(num.sum(amplitudes[mask_dec]))
1255 if sum_dec != 0.0:
1256 amplitudes[mask_dec] /= sum_dec
1257 else:
1258 amplitudes[mask_dec] = 1.0 / n_dec
1260 if self.mute_acceleration:
1261 amplitudes[mask_acc] = 0.0
1263 if self.mute_deceleration:
1264 amplitudes[mask_dec] = 0.0
1266 return times, amplitudes
1269class STFMode(StringChoice):
1270 choices = ['pre', 'post']
1273class Source(Location, Cloneable):
1274 '''
1275 Base class for all source models.
1276 '''
1278 name = String.T(optional=True, default='')
1280 time = Timestamp.T(
1281 default=Timestamp.D('1970-01-01 00:00:00'),
1282 help='source origin time.')
1284 stf = STF.T(
1285 optional=True,
1286 help='source time function.')
1288 stf_mode = STFMode.T(
1289 default='post',
1290 help='whether to apply source time function in pre or '
1291 'post-processing.')
1293 def __init__(self, **kwargs):
1294 Location.__init__(self, **kwargs)
1296 def update(self, **kwargs):
1297 '''
1298 Change some of the source models parameters.
1300 Example::
1302 >>> from pyrocko import gf
1303 >>> s = gf.DCSource()
1304 >>> s.update(strike=66., dip=33.)
1305 >>> print(s)
1306 --- !pf.DCSource
1307 depth: 0.0
1308 time: 1970-01-01 00:00:00
1309 magnitude: 6.0
1310 strike: 66.0
1311 dip: 33.0
1312 rake: 0.0
1314 '''
1316 for (k, v) in kwargs.items():
1317 self[k] = v
1319 def grid(self, **variables):
1320 '''
1321 Create grid of source model variations.
1323 :returns: :py:class:`SourceGrid` instance.
1325 Example::
1327 >>> from pyrocko import gf
1328 >>> base = DCSource()
1329 >>> R = gf.Range
1330 >>> for s in base.grid(R('
1332 '''
1333 return SourceGrid(base=self, variables=variables)
1335 def base_key(self):
1336 '''
1337 Get key to decide about source discretization / GF stack sharing.
1339 When two source models differ only in amplitude and origin time, the
1340 discretization and the GF stacking can be done only once for a unit
1341 amplitude and a zero origin time and the amplitude and origin times of
1342 the seismograms can be applied during post-processing of the synthetic
1343 seismogram.
1345 For any derived parameterized source model, this method is called to
1346 decide if discretization and stacking of the source should be shared.
1347 When two source models return an equal vector of values discretization
1348 is shared.
1349 '''
1350 return (self.depth, self.lat, self.north_shift,
1351 self.lon, self.east_shift, self.time, type(self).__name__) + \
1352 self.effective_stf_pre().base_key()
1354 def get_factor(self):
1355 '''
1356 Get the scaling factor to be applied during post-processing.
1358 Discretization of the base seismogram is usually done for a unit
1359 amplitude, because a common factor can be efficiently multiplied to
1360 final seismograms. This eliminates to do repeat the stacking when
1361 creating seismograms for a series of source models only differing in
1362 amplitude.
1364 This method should return the scaling factor to apply in the
1365 post-processing (often this is simply the scalar moment of the source).
1366 '''
1368 return 1.0
1370 def effective_stf_pre(self):
1371 '''
1372 Return the STF applied before stacking of the Green's functions.
1374 This STF is used during discretization of the parameterized source
1375 models, i.e. to produce a temporal distribution of point sources.
1377 Handling of the STF before stacking of the GFs is less efficient but
1378 allows to use different source time functions for different parts of
1379 the source.
1380 '''
1382 if self.stf is not None and self.stf_mode == 'pre':
1383 return self.stf
1384 else:
1385 return g_unit_pulse
1387 def effective_stf_post(self):
1388 '''
1389 Return the STF applied after stacking of the Green's fuctions.
1391 This STF is used in the post-processing of the synthetic seismograms.
1393 Handling of the STF after stacking of the GFs is usually more efficient
1394 but is only possible when a common STF is used for all subsources.
1395 '''
1397 if self.stf is not None and self.stf_mode == 'post':
1398 return self.stf
1399 else:
1400 return g_unit_pulse
1402 def _dparams_base(self):
1403 return dict(times=arr(self.time),
1404 lat=self.lat, lon=self.lon,
1405 north_shifts=arr(self.north_shift),
1406 east_shifts=arr(self.east_shift),
1407 depths=arr(self.depth))
1409 def _hash(self):
1410 sha = sha1()
1411 for k in self.base_key():
1412 sha.update(str(k).encode())
1413 return sha.hexdigest()
1415 def _dparams_base_repeated(self, times):
1416 if times is None:
1417 return self._dparams_base()
1419 nt = times.size
1420 north_shifts = num.repeat(self.north_shift, nt)
1421 east_shifts = num.repeat(self.east_shift, nt)
1422 depths = num.repeat(self.depth, nt)
1423 return dict(times=times,
1424 lat=self.lat, lon=self.lon,
1425 north_shifts=north_shifts,
1426 east_shifts=east_shifts,
1427 depths=depths)
1429 def pyrocko_event(self, store=None, target=None, **kwargs):
1430 duration = None
1431 if self.stf:
1432 duration = self.stf.effective_duration
1434 return model.Event(
1435 lat=self.lat,
1436 lon=self.lon,
1437 north_shift=self.north_shift,
1438 east_shift=self.east_shift,
1439 time=self.time,
1440 name=self.name,
1441 depth=self.depth,
1442 duration=duration,
1443 **kwargs)
1445 def geometry(self, **kwargs):
1446 raise NotImplementedError
1448 def outline(self, cs='xyz'):
1449 points = num.atleast_2d(num.zeros([1, 3]))
1451 points[:, 0] += self.north_shift
1452 points[:, 1] += self.east_shift
1453 points[:, 2] += self.depth
1454 if cs == 'xyz':
1455 return points
1456 elif cs == 'xy':
1457 return points[:, :2]
1458 elif cs in ('latlon', 'lonlat'):
1459 latlon = ne_to_latlon(
1460 self.lat, self.lon, points[:, 0], points[:, 1])
1462 latlon = num.array(latlon).T
1463 if cs == 'latlon':
1464 return latlon
1465 else:
1466 return latlon[:, ::-1]
1468 @classmethod
1469 def from_pyrocko_event(cls, ev, **kwargs):
1470 if ev.depth is None:
1471 raise ConversionError(
1472 'Cannot convert event object to source object: '
1473 'no depth information available')
1475 stf = None
1476 if ev.duration is not None:
1477 stf = HalfSinusoidSTF(effective_duration=ev.duration)
1479 d = dict(
1480 name=ev.name,
1481 time=ev.time,
1482 lat=ev.lat,
1483 lon=ev.lon,
1484 north_shift=ev.north_shift,
1485 east_shift=ev.east_shift,
1486 depth=ev.depth,
1487 stf=stf)
1488 d.update(kwargs)
1489 return cls(**d)
1491 def get_magnitude(self):
1492 raise NotImplementedError(
1493 '%s does not implement get_magnitude()'
1494 % self.__class__.__name__)
1497class SourceWithMagnitude(Source):
1498 '''
1499 Base class for sources containing a moment magnitude.
1500 '''
1502 magnitude = Float.T(
1503 default=6.0,
1504 help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1506 def __init__(self, **kwargs):
1507 if 'moment' in kwargs:
1508 mom = kwargs.pop('moment')
1509 if 'magnitude' not in kwargs:
1510 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1512 Source.__init__(self, **kwargs)
1514 @property
1515 def moment(self):
1516 return float(pmt.magnitude_to_moment(self.magnitude))
1518 @moment.setter
1519 def moment(self, value):
1520 self.magnitude = float(pmt.moment_to_magnitude(value))
1522 def pyrocko_event(self, store=None, target=None, **kwargs):
1523 return Source.pyrocko_event(
1524 self, store, target,
1525 magnitude=self.magnitude,
1526 **kwargs)
1528 @classmethod
1529 def from_pyrocko_event(cls, ev, **kwargs):
1530 d = {}
1531 if ev.magnitude:
1532 d.update(magnitude=ev.magnitude)
1534 d.update(kwargs)
1535 return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
1537 def get_magnitude(self):
1538 return self.magnitude
1541class DerivedMagnitudeError(ValidationError):
1542 '''
1543 Raised when conversion between magnitude, moment, volume change or
1544 displacement failed.
1545 '''
1546 pass
1549class SourceWithDerivedMagnitude(Source):
1551 class __T(Source.T):
1553 def validate_extra(self, val):
1554 Source.T.validate_extra(self, val)
1555 val.check_conflicts()
1557 def check_conflicts(self):
1558 '''
1559 Check for parameter conflicts.
1561 To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError`
1562 on conflicts.
1563 '''
1564 pass
1566 def get_magnitude(self, store=None, target=None):
1567 raise DerivedMagnitudeError('No magnitude set.')
1569 def get_moment(self, store=None, target=None):
1570 return float(pmt.magnitude_to_moment(
1571 self.get_magnitude(store, target)))
1573 def pyrocko_moment_tensor(self, store=None, target=None):
1574 raise NotImplementedError(
1575 '%s does not implement pyrocko_moment_tensor()'
1576 % self.__class__.__name__)
1578 def pyrocko_event(self, store=None, target=None, **kwargs):
1579 try:
1580 mt = self.pyrocko_moment_tensor(store, target)
1581 magnitude = self.get_magnitude()
1582 except (DerivedMagnitudeError, NotImplementedError):
1583 mt = None
1584 magnitude = None
1586 return Source.pyrocko_event(
1587 self, store, target,
1588 moment_tensor=mt,
1589 magnitude=magnitude,
1590 **kwargs)
1593class ExplosionSource(SourceWithDerivedMagnitude):
1594 '''
1595 An isotropic explosion point source.
1596 '''
1598 magnitude = Float.T(
1599 optional=True,
1600 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1602 volume_change = Float.T(
1603 optional=True,
1604 help='volume change of the explosion/implosion or '
1605 'the contracting/extending magmatic source. [m^3]')
1607 discretized_source_class = meta.DiscretizedExplosionSource
1609 def __init__(self, **kwargs):
1610 if 'moment' in kwargs:
1611 mom = kwargs.pop('moment')
1612 if 'magnitude' not in kwargs:
1613 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1615 SourceWithDerivedMagnitude.__init__(self, **kwargs)
1617 def base_key(self):
1618 return SourceWithDerivedMagnitude.base_key(self) + \
1619 (self.volume_change,)
1621 def check_conflicts(self):
1622 if self.magnitude is not None and self.volume_change is not None:
1623 raise DerivedMagnitudeError(
1624 'Magnitude and volume_change are both defined.')
1626 def get_magnitude(self, store=None, target=None):
1627 self.check_conflicts()
1629 if self.magnitude is not None:
1630 return self.magnitude
1632 elif self.volume_change is not None:
1633 moment = self.volume_change * \
1634 self.get_moment_to_volume_change_ratio(store, target)
1636 return float(pmt.moment_to_magnitude(abs(moment)))
1637 else:
1638 return float(pmt.moment_to_magnitude(1.0))
1640 def get_volume_change(self, store=None, target=None):
1641 self.check_conflicts()
1643 if self.volume_change is not None:
1644 return self.volume_change
1646 elif self.magnitude is not None:
1647 moment = float(pmt.magnitude_to_moment(self.magnitude))
1648 return moment / self.get_moment_to_volume_change_ratio(
1649 store, target)
1651 else:
1652 return 1.0 / self.get_moment_to_volume_change_ratio(store)
1654 def get_moment_to_volume_change_ratio(self, store, target=None):
1655 if store is None:
1656 raise DerivedMagnitudeError(
1657 'Need earth model to convert between volume change and '
1658 'magnitude.')
1660 points = num.array(
1661 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1663 interpolation = target.interpolation if target else 'multilinear'
1664 try:
1665 shear_moduli = store.config.get_shear_moduli(
1666 self.lat, self.lon,
1667 points=points,
1668 interpolation=interpolation)[0]
1670 bulk_moduli = store.config.get_bulk_moduli(
1671 self.lat, self.lon,
1672 points=points,
1673 interpolation=interpolation)[0]
1674 except meta.OutOfBounds:
1675 raise DerivedMagnitudeError(
1676 'Could not get shear modulus at source position.')
1678 return float(2. * shear_moduli + bulk_moduli)
1680 def get_factor(self):
1681 return 1.0
1683 def discretize_basesource(self, store, target=None):
1684 times, amplitudes = self.effective_stf_pre().discretize_t(
1685 store.config.deltat, self.time)
1687 amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.)
1689 if self.volume_change is not None:
1690 if self.volume_change < 0.:
1691 amplitudes *= -1
1693 return meta.DiscretizedExplosionSource(
1694 m0s=amplitudes,
1695 **self._dparams_base_repeated(times))
1697 def pyrocko_moment_tensor(self, store=None, target=None):
1698 a = self.get_moment(store, target) * math.sqrt(2. / 3.)
1699 return pmt.MomentTensor(m=pmt.symmat6(a, a, a, 0., 0., 0.))
1702class RectangularExplosionSource(ExplosionSource):
1703 '''
1704 Rectangular or line explosion source.
1705 '''
1707 discretized_source_class = meta.DiscretizedExplosionSource
1709 strike = Float.T(
1710 default=0.0,
1711 help='strike direction in [deg], measured clockwise from north')
1713 dip = Float.T(
1714 default=90.0,
1715 help='dip angle in [deg], measured downward from horizontal')
1717 length = Float.T(
1718 default=0.,
1719 help='length of rectangular source area [m]')
1721 width = Float.T(
1722 default=0.,
1723 help='width of rectangular source area [m]')
1725 anchor = StringChoice.T(
1726 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
1727 'bottom_left', 'bottom_right'],
1728 default='center',
1729 optional=True,
1730 help='Anchor point for positioning the plane, can be: top, center or'
1731 'bottom and also top_left, top_right,bottom_left,'
1732 'bottom_right, center_left and center right')
1734 nucleation_x = Float.T(
1735 optional=True,
1736 help='horizontal position of rupture nucleation in normalized fault '
1737 'plane coordinates (-1 = left edge, +1 = right edge)')
1739 nucleation_y = Float.T(
1740 optional=True,
1741 help='down-dip position of rupture nucleation in normalized fault '
1742 'plane coordinates (-1 = upper edge, +1 = lower edge)')
1744 velocity = Float.T(
1745 default=3500.,
1746 help='speed of explosion front [m/s]')
1748 aggressive_oversampling = Bool.T(
1749 default=False,
1750 help='Aggressive oversampling for basesource discretization. '
1751 "When using 'multilinear' interpolation oversampling has"
1752 ' practically no effect.')
1754 def base_key(self):
1755 return Source.base_key(self) + (self.strike, self.dip, self.length,
1756 self.width, self.nucleation_x,
1757 self.nucleation_y, self.velocity,
1758 self.anchor)
1760 def discretize_basesource(self, store, target=None):
1762 if self.nucleation_x is not None:
1763 nucx = self.nucleation_x * 0.5 * self.length
1764 else:
1765 nucx = None
1767 if self.nucleation_y is not None:
1768 nucy = self.nucleation_y * 0.5 * self.width
1769 else:
1770 nucy = None
1772 stf = self.effective_stf_pre()
1774 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
1775 store.config.deltas, store.config.deltat,
1776 self.time, self.north_shift, self.east_shift, self.depth,
1777 self.strike, self.dip, self.length, self.width, self.anchor,
1778 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy)
1780 amplitudes /= num.sum(amplitudes)
1781 amplitudes *= self.get_moment(store, target)
1783 return meta.DiscretizedExplosionSource(
1784 lat=self.lat,
1785 lon=self.lon,
1786 times=times,
1787 north_shifts=points[:, 0],
1788 east_shifts=points[:, 1],
1789 depths=points[:, 2],
1790 m0s=amplitudes)
1792 def outline(self, cs='xyz'):
1793 points = outline_rect_source(self.strike, self.dip, self.length,
1794 self.width, self.anchor)
1796 points[:, 0] += self.north_shift
1797 points[:, 1] += self.east_shift
1798 points[:, 2] += self.depth
1799 if cs == 'xyz':
1800 return points
1801 elif cs == 'xy':
1802 return points[:, :2]
1803 elif cs in ('latlon', 'lonlat'):
1804 latlon = ne_to_latlon(
1805 self.lat, self.lon, points[:, 0], points[:, 1])
1807 latlon = num.array(latlon).T
1808 if cs == 'latlon':
1809 return latlon
1810 else:
1811 return latlon[:, ::-1]
1813 def get_nucleation_abs_coord(self, cs='xy'):
1815 if self.nucleation_x is None:
1816 return None, None
1818 coords = from_plane_coords(self.strike, self.dip, self.length,
1819 self.width, self.depth, self.nucleation_x,
1820 self.nucleation_y, lat=self.lat,
1821 lon=self.lon, north_shift=self.north_shift,
1822 east_shift=self.east_shift, cs=cs)
1823 return coords
1826class DCSource(SourceWithMagnitude):
1827 '''
1828 A double-couple point source.
1829 '''
1831 strike = Float.T(
1832 default=0.0,
1833 help='strike direction in [deg], measured clockwise from north')
1835 dip = Float.T(
1836 default=90.0,
1837 help='dip angle in [deg], measured downward from horizontal')
1839 rake = Float.T(
1840 default=0.0,
1841 help='rake angle in [deg], '
1842 'measured counter-clockwise from right-horizontal '
1843 'in on-plane view')
1845 discretized_source_class = meta.DiscretizedMTSource
1847 def base_key(self):
1848 return Source.base_key(self) + (self.strike, self.dip, self.rake)
1850 def get_factor(self):
1851 return float(pmt.magnitude_to_moment(self.magnitude))
1853 def discretize_basesource(self, store, target=None):
1854 mot = pmt.MomentTensor(
1855 strike=self.strike, dip=self.dip, rake=self.rake)
1857 times, amplitudes = self.effective_stf_pre().discretize_t(
1858 store.config.deltat, self.time)
1859 return meta.DiscretizedMTSource(
1860 m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis],
1861 **self._dparams_base_repeated(times))
1863 def pyrocko_moment_tensor(self, store=None, target=None):
1864 return pmt.MomentTensor(
1865 strike=self.strike,
1866 dip=self.dip,
1867 rake=self.rake,
1868 scalar_moment=self.moment)
1870 def pyrocko_event(self, store=None, target=None, **kwargs):
1871 return SourceWithMagnitude.pyrocko_event(
1872 self, store, target,
1873 moment_tensor=self.pyrocko_moment_tensor(store, target),
1874 **kwargs)
1876 @classmethod
1877 def from_pyrocko_event(cls, ev, **kwargs):
1878 d = {}
1879 mt = ev.moment_tensor
1880 if mt:
1881 (strike, dip, rake), _ = mt.both_strike_dip_rake()
1882 d.update(
1883 strike=float(strike),
1884 dip=float(dip),
1885 rake=float(rake),
1886 magnitude=float(mt.moment_magnitude()))
1888 d.update(kwargs)
1889 return super(DCSource, cls).from_pyrocko_event(ev, **d)
1892class CLVDSource(SourceWithMagnitude):
1893 '''
1894 A pure CLVD point source.
1895 '''
1897 discretized_source_class = meta.DiscretizedMTSource
1899 azimuth = Float.T(
1900 default=0.0,
1901 help='azimuth direction of largest dipole, clockwise from north [deg]')
1903 dip = Float.T(
1904 default=90.,
1905 help='dip direction of largest dipole, downward from horizontal [deg]')
1907 def base_key(self):
1908 return Source.base_key(self) + (self.azimuth, self.dip)
1910 def get_factor(self):
1911 return float(pmt.magnitude_to_moment(self.magnitude))
1913 @property
1914 def m6(self):
1915 a = math.sqrt(4. / 3.) * self.get_factor()
1916 m = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1917 rotmat1 = pmt.euler_to_matrix(
1918 d2r * (self.dip - 90.),
1919 d2r * (self.azimuth - 90.),
1920 0.)
1921 m = num.dot(rotmat1.T, num.dot(m, rotmat1))
1922 return pmt.to6(m)
1924 @property
1925 def m6_astuple(self):
1926 return tuple(self.m6.tolist())
1928 def discretize_basesource(self, store, target=None):
1929 factor = self.get_factor()
1930 times, amplitudes = self.effective_stf_pre().discretize_t(
1931 store.config.deltat, self.time)
1932 return meta.DiscretizedMTSource(
1933 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor,
1934 **self._dparams_base_repeated(times))
1936 def pyrocko_moment_tensor(self, store=None, target=None):
1937 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
1939 def pyrocko_event(self, store=None, target=None, **kwargs):
1940 mt = self.pyrocko_moment_tensor(store, target)
1941 return Source.pyrocko_event(
1942 self, store, target,
1943 moment_tensor=self.pyrocko_moment_tensor(store, target),
1944 magnitude=float(mt.moment_magnitude()),
1945 **kwargs)
1948class VLVDSource(SourceWithDerivedMagnitude):
1949 '''
1950 Volumetric linear vector dipole source.
1952 This source is a parameterization for a restricted moment tensor point
1953 source, useful to represent dyke or sill like inflation or deflation
1954 sources. The restriction is such that the moment tensor is rotational
1955 symmetric. It can be represented by a superposition of a linear vector
1956 dipole (here we use a CLVD for convenience) and an isotropic component. The
1957 restricted moment tensor has 4 degrees of freedom: 2 independent
1958 eigenvalues and 2 rotation angles orienting the the symmetry axis.
1960 In this parameterization, the isotropic component is controlled by
1961 ``volume_change``. To define the moment tensor, it must be converted to the
1962 scalar moment of the the MT's isotropic component. For the conversion, the
1963 shear modulus at the source's position must be known. This value is
1964 extracted from the earth model defined in the GF store in use.
1966 The CLVD part by controlled by its scalar moment :math:`M_0`:
1967 ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a
1968 positiv or negativ CLVD (the sign of the largest eigenvalue).
1969 '''
1971 discretized_source_class = meta.DiscretizedMTSource
1973 azimuth = Float.T(
1974 default=0.0,
1975 help='azimuth direction of symmetry axis, clockwise from north [deg].')
1977 dip = Float.T(
1978 default=90.,
1979 help='dip direction of symmetry axis, downward from horizontal [deg].')
1981 volume_change = Float.T(
1982 default=0.,
1983 help='volume change of the inflation/deflation [m^3].')
1985 clvd_moment = Float.T(
1986 default=0.,
1987 help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign '
1988 'controls the sign of the CLVD (the sign of its largest '
1989 'eigenvalue).')
1991 def get_moment_to_volume_change_ratio(self, store, target):
1992 if store is None or target is None:
1993 raise DerivedMagnitudeError(
1994 'Need earth model to convert between volume change and '
1995 'magnitude.')
1997 points = num.array(
1998 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
2000 try:
2001 shear_moduli = store.config.get_shear_moduli(
2002 self.lat, self.lon,
2003 points=points,
2004 interpolation=target.interpolation)[0]
2006 bulk_moduli = store.config.get_bulk_moduli(
2007 self.lat, self.lon,
2008 points=points,
2009 interpolation=target.interpolation)[0]
2010 except meta.OutOfBounds:
2011 raise DerivedMagnitudeError(
2012 'Could not get shear modulus at source position.')
2014 return float(2. * shear_moduli + bulk_moduli)
2016 def base_key(self):
2017 return Source.base_key(self) + \
2018 (self.azimuth, self.dip, self.volume_change, self.clvd_moment)
2020 def get_magnitude(self, store=None, target=None):
2021 mt = self.pyrocko_moment_tensor(store, target)
2022 return float(pmt.moment_to_magnitude(mt.moment))
2024 def get_m6(self, store, target):
2025 a = math.sqrt(4. / 3.) * self.clvd_moment
2026 m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
2028 rotmat1 = pmt.euler_to_matrix(
2029 d2r * (self.dip - 90.),
2030 d2r * (self.azimuth - 90.),
2031 0.)
2032 m_clvd = num.dot(rotmat1.T, num.dot(m_clvd, rotmat1))
2034 m_iso = self.volume_change * \
2035 self.get_moment_to_volume_change_ratio(store, target)
2037 m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0.,
2038 0., 0.,) * math.sqrt(2. / 3)
2040 m = pmt.to6(m_clvd) + pmt.to6(m_iso)
2041 return m
2043 def get_moment(self, store=None, target=None):
2044 return float(pmt.magnitude_to_moment(
2045 self.get_magnitude(store, target)))
2047 def get_m6_astuple(self, store, target):
2048 m6 = self.get_m6(store, target)
2049 return tuple(m6.tolist())
2051 def discretize_basesource(self, store, target=None):
2052 times, amplitudes = self.effective_stf_pre().discretize_t(
2053 store.config.deltat, self.time)
2055 m6 = self.get_m6(store, target)
2056 m6 *= amplitudes / self.get_factor()
2058 return meta.DiscretizedMTSource(
2059 m6s=m6[num.newaxis, :],
2060 **self._dparams_base_repeated(times))
2062 def pyrocko_moment_tensor(self, store=None, target=None):
2063 m6_astuple = self.get_m6_astuple(store, target)
2064 return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple))
2067class MTSource(Source):
2068 '''
2069 A moment tensor point source.
2070 '''
2072 discretized_source_class = meta.DiscretizedMTSource
2074 mnn = Float.T(
2075 default=1.,
2076 help='north-north component of moment tensor in [Nm]')
2078 mee = Float.T(
2079 default=1.,
2080 help='east-east component of moment tensor in [Nm]')
2082 mdd = Float.T(
2083 default=1.,
2084 help='down-down component of moment tensor in [Nm]')
2086 mne = Float.T(
2087 default=0.,
2088 help='north-east component of moment tensor in [Nm]')
2090 mnd = Float.T(
2091 default=0.,
2092 help='north-down component of moment tensor in [Nm]')
2094 med = Float.T(
2095 default=0.,
2096 help='east-down component of moment tensor in [Nm]')
2098 def __init__(self, **kwargs):
2099 if 'm6' in kwargs:
2100 for (k, v) in zip('mnn mee mdd mne mnd med'.split(),
2101 kwargs.pop('m6')):
2102 kwargs[k] = float(v)
2104 Source.__init__(self, **kwargs)
2106 @property
2107 def m6(self):
2108 return num.array(self.m6_astuple)
2110 @property
2111 def m6_astuple(self):
2112 return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)
2114 @m6.setter
2115 def m6(self, value):
2116 self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value
2118 def base_key(self):
2119 return Source.base_key(self) + self.m6_astuple
2121 def discretize_basesource(self, store, target=None):
2122 times, amplitudes = self.effective_stf_pre().discretize_t(
2123 store.config.deltat, self.time)
2124 return meta.DiscretizedMTSource(
2125 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis],
2126 **self._dparams_base_repeated(times))
2128 def get_magnitude(self, store=None, target=None):
2129 m6 = self.m6
2130 return pmt.moment_to_magnitude(
2131 math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) /
2132 math.sqrt(2.))
2134 def pyrocko_moment_tensor(self, store=None, target=None):
2135 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
2137 def pyrocko_event(self, store=None, target=None, **kwargs):
2138 mt = self.pyrocko_moment_tensor(store, target)
2139 return Source.pyrocko_event(
2140 self, store, target,
2141 moment_tensor=self.pyrocko_moment_tensor(store, target),
2142 magnitude=float(mt.moment_magnitude()),
2143 **kwargs)
2145 @classmethod
2146 def from_pyrocko_event(cls, ev, **kwargs):
2147 d = {}
2148 mt = ev.moment_tensor
2149 if mt:
2150 d.update(m6=tuple(map(float, mt.m6())))
2151 else:
2152 if ev.magnitude is not None:
2153 mom = pmt.magnitude_to_moment(ev.magnitude)
2154 v = math.sqrt(2. / 3.) * mom
2155 d.update(m6=(v, v, v, 0., 0., 0.))
2157 d.update(kwargs)
2158 return super(MTSource, cls).from_pyrocko_event(ev, **d)
2161map_anchor = {
2162 'center': (0.0, 0.0),
2163 'center_left': (-1.0, 0.0),
2164 'center_right': (1.0, 0.0),
2165 'top': (0.0, -1.0),
2166 'top_left': (-1.0, -1.0),
2167 'top_right': (1.0, -1.0),
2168 'bottom': (0.0, 1.0),
2169 'bottom_left': (-1.0, 1.0),
2170 'bottom_right': (1.0, 1.0)}
2173class RectangularSource(SourceWithDerivedMagnitude):
2174 '''
2175 Classical Haskell source model modified for bilateral rupture.
2176 '''
2178 discretized_source_class = meta.DiscretizedMTSource
2180 magnitude = Float.T(
2181 optional=True,
2182 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
2184 strike = Float.T(
2185 default=0.0,
2186 help='strike direction in [deg], measured clockwise from north')
2188 dip = Float.T(
2189 default=90.0,
2190 help='dip angle in [deg], measured downward from horizontal')
2192 rake = Float.T(
2193 default=0.0,
2194 help='rake angle in [deg], '
2195 'measured counter-clockwise from right-horizontal '
2196 'in on-plane view')
2198 length = Float.T(
2199 default=0.,
2200 help='length of rectangular source area [m]')
2202 width = Float.T(
2203 default=0.,
2204 help='width of rectangular source area [m]')
2206 anchor = StringChoice.T(
2207 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2208 'bottom_left', 'bottom_right'],
2209 default='center',
2210 optional=True,
2211 help='Anchor point for positioning the plane, can be: ``top, center '
2212 'bottom, top_left, top_right,bottom_left,'
2213 'bottom_right, center_left, center right``.')
2215 nucleation_x = Float.T(
2216 optional=True,
2217 help='horizontal position of rupture nucleation in normalized fault '
2218 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge)')
2220 nucleation_y = Float.T(
2221 optional=True,
2222 help='down-dip position of rupture nucleation in normalized fault '
2223 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge)')
2225 velocity = Float.T(
2226 default=3500.,
2227 help='speed of rupture front [m/s]')
2229 slip = Float.T(
2230 optional=True,
2231 help='Slip on the rectangular source area [m]')
2233 opening_fraction = Float.T(
2234 default=0.,
2235 help='Determines fraction of slip related to opening. '
2236 '(``-1``: pure tensile closing, '
2237 '``0``: pure shear, '
2238 '``1``: pure tensile opening)')
2240 decimation_factor = Int.T(
2241 optional=True,
2242 default=1,
2243 help='Sub-source decimation factor, a larger decimation will'
2244 ' make the result inaccurate but shorten the necessary'
2245 ' computation time (use for testing puposes only).')
2247 aggressive_oversampling = Bool.T(
2248 default=False,
2249 help='Aggressive oversampling for basesource discretization. '
2250 "When using 'multilinear' interpolation oversampling has"
2251 ' practically no effect.')
2253 def __init__(self, **kwargs):
2254 if 'moment' in kwargs:
2255 mom = kwargs.pop('moment')
2256 if 'magnitude' not in kwargs:
2257 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
2259 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2261 def base_key(self):
2262 return SourceWithDerivedMagnitude.base_key(self) + (
2263 self.magnitude,
2264 self.slip,
2265 self.strike,
2266 self.dip,
2267 self.rake,
2268 self.length,
2269 self.width,
2270 self.nucleation_x,
2271 self.nucleation_y,
2272 self.velocity,
2273 self.decimation_factor,
2274 self.anchor)
2276 def check_conflicts(self):
2277 if self.magnitude is not None and self.slip is not None:
2278 raise DerivedMagnitudeError(
2279 'Magnitude and slip are both defined.')
2281 def get_magnitude(self, store=None, target=None):
2282 self.check_conflicts()
2283 if self.magnitude is not None:
2284 return self.magnitude
2286 elif self.slip is not None:
2287 if None in (store, target):
2288 raise DerivedMagnitudeError(
2289 'Magnitude for a rectangular source with slip defined '
2290 'can only be derived when earth model and target '
2291 'interpolation method are available.')
2293 amplitudes = self._discretize(store, target)[2]
2294 if amplitudes.ndim == 2:
2295 # CLVD component has no net moment, leave out
2296 return float(pmt.moment_to_magnitude(
2297 num.sum(num.abs(amplitudes[0:2, :]).sum())))
2298 else:
2299 return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
2301 else:
2302 return float(pmt.moment_to_magnitude(1.0))
2304 def get_factor(self):
2305 return 1.0
2307 def get_slip_tensile(self):
2308 return self.slip * self.opening_fraction
2310 def get_slip_shear(self):
2311 return self.slip - abs(self.get_slip_tensile)
2313 def _discretize(self, store, target):
2314 if self.nucleation_x is not None:
2315 nucx = self.nucleation_x * 0.5 * self.length
2316 else:
2317 nucx = None
2319 if self.nucleation_y is not None:
2320 nucy = self.nucleation_y * 0.5 * self.width
2321 else:
2322 nucy = None
2324 stf = self.effective_stf_pre()
2326 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
2327 store.config.deltas, store.config.deltat,
2328 self.time, self.north_shift, self.east_shift, self.depth,
2329 self.strike, self.dip, self.length, self.width, self.anchor,
2330 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy,
2331 decimation_factor=self.decimation_factor,
2332 aggressive_oversampling=self.aggressive_oversampling)
2334 if self.slip is not None:
2335 if target is not None:
2336 interpolation = target.interpolation
2337 else:
2338 interpolation = 'nearest_neighbor'
2339 logger.warning(
2340 'no target information available, will use '
2341 '"nearest_neighbor" interpolation when extracting shear '
2342 'modulus from earth model')
2344 shear_moduli = store.config.get_shear_moduli(
2345 self.lat, self.lon,
2346 points=points,
2347 interpolation=interpolation)
2349 tensile_slip = self.get_slip_tensile()
2350 shear_slip = self.slip - abs(tensile_slip)
2352 amplitudes_total = [shear_moduli * shear_slip]
2353 if tensile_slip != 0:
2354 bulk_moduli = store.config.get_bulk_moduli(
2355 self.lat, self.lon,
2356 points=points,
2357 interpolation=interpolation)
2359 tensile_iso = bulk_moduli * tensile_slip
2360 tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip
2362 amplitudes_total.extend([tensile_iso, tensile_clvd])
2364 amplitudes_total = num.vstack(amplitudes_total).squeeze() * \
2365 amplitudes * dl * dw
2367 else:
2368 # normalization to retain total moment
2369 amplitudes_norm = amplitudes / num.sum(amplitudes)
2370 moment = self.get_moment(store, target)
2372 amplitudes_total = [
2373 amplitudes_norm * moment * (1 - abs(self.opening_fraction))]
2374 if self.opening_fraction != 0.:
2375 amplitudes_total.append(
2376 amplitudes_norm * self.opening_fraction * moment)
2378 amplitudes_total = num.vstack(amplitudes_total).squeeze()
2380 return points, times, num.atleast_1d(amplitudes_total), dl, dw, nl, nw
2382 def discretize_basesource(self, store, target=None):
2384 points, times, amplitudes, dl, dw, nl, nw = self._discretize(
2385 store, target)
2387 mot = pmt.MomentTensor(
2388 strike=self.strike, dip=self.dip, rake=self.rake)
2389 m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0)
2391 if amplitudes.ndim == 1:
2392 m6s[:, :] *= amplitudes[:, num.newaxis]
2393 elif amplitudes.ndim == 2:
2394 # shear MT components
2395 rotmat1 = pmt.euler_to_matrix(
2396 d2r * self.dip, d2r * self.strike, d2r * -self.rake)
2397 m6s[:, :] *= amplitudes[0, :][:, num.newaxis]
2399 if amplitudes.shape[0] == 2:
2400 # tensile MT components - moment/magnitude input
2401 tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.)
2402 rot_tensile = pmt.to6(
2403 num.dot(rotmat1.T, num.dot(tensile, rotmat1)))
2405 m6s_tensile = rot_tensile[
2406 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2407 m6s += m6s_tensile
2409 elif amplitudes.shape[0] == 3:
2410 # tensile MT components - slip input
2411 iso = pmt.symmat6(1., 1., 1., 0., 0., 0.)
2412 clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.)
2414 rot_iso = pmt.to6(
2415 num.dot(rotmat1.T, num.dot(iso, rotmat1)))
2416 rot_clvd = pmt.to6(
2417 num.dot(rotmat1.T, num.dot(clvd, rotmat1)))
2419 m6s_iso = rot_iso[
2420 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2421 m6s_clvd = rot_clvd[
2422 num.newaxis, :] * amplitudes[2, :][:, num.newaxis]
2423 m6s += m6s_iso + m6s_clvd
2424 else:
2425 raise ValueError('Unknwown amplitudes shape!')
2426 else:
2427 raise ValueError(
2428 'Unexpected dimension of {}'.format(amplitudes.ndim))
2430 ds = meta.DiscretizedMTSource(
2431 lat=self.lat,
2432 lon=self.lon,
2433 times=times,
2434 north_shifts=points[:, 0],
2435 east_shifts=points[:, 1],
2436 depths=points[:, 2],
2437 m6s=m6s,
2438 dl=dl,
2439 dw=dw,
2440 nl=nl,
2441 nw=nw)
2443 return ds
2445 def xy_to_coord(self, x, y, cs='xyz'):
2446 ln, wd = self.length, self.width
2447 strike, dip = self.strike, self.dip
2449 def array_check(variable):
2450 if not isinstance(variable, num.ndarray):
2451 return num.array(variable)
2452 else:
2453 return variable
2455 x, y = array_check(x), array_check(y)
2457 if x.shape[0] != y.shape[0]:
2458 raise ValueError('Shapes of x and y mismatch')
2460 x, y = x * 0.5 * ln, y * 0.5 * wd
2462 points = num.hstack((
2463 x.reshape(-1, 1), y.reshape(-1, 1), num.zeros((x.shape[0], 1))))
2465 anch_x, anch_y = map_anchor[self.anchor]
2466 points[:, 0] -= anch_x * 0.5 * ln
2467 points[:, 1] -= anch_y * 0.5 * wd
2469 rotmat = num.asarray(
2470 pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
2472 points_rot = num.dot(rotmat.T, points.T).T
2474 points_rot[:, 0] += self.north_shift
2475 points_rot[:, 1] += self.east_shift
2476 points_rot[:, 2] += self.depth
2478 if cs == 'xyz':
2479 return points_rot
2480 elif cs == 'xy':
2481 return points_rot[:, :2]
2482 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2483 latlon = ne_to_latlon(
2484 self.lat, self.lon, points_rot[:, 0], points_rot[:, 1])
2485 latlon = num.array(latlon).T
2486 if cs == 'latlon':
2487 return latlon
2488 elif cs == 'lonlat':
2489 return latlon[:, ::-1]
2490 else:
2491 return num.concatenate(
2492 (latlon, points_rot[:, 2].reshape((len(points_rot), 1))),
2493 axis=1)
2495 def outline(self, cs='xyz'):
2496 x = num.array([-1., 1., 1., -1., -1.])
2497 y = num.array([-1., -1., 1., 1., -1.])
2499 return self.xy_to_coord(x, y, cs=cs)
2501 def points_on_source(self, cs='xyz', **kwargs):
2503 points = points_on_rect_source(
2504 self.strike, self.dip, self.length, self.width,
2505 self.anchor, **kwargs)
2507 points[:, 0] += self.north_shift
2508 points[:, 1] += self.east_shift
2509 points[:, 2] += self.depth
2510 if cs == 'xyz':
2511 return points
2512 elif cs == 'xy':
2513 return points[:, :2]
2514 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2515 latlon = ne_to_latlon(
2516 self.lat, self.lon, points[:, 0], points[:, 1])
2518 latlon = num.array(latlon).T
2519 if cs == 'latlon':
2520 return latlon
2521 elif cs == 'lonlat':
2522 return latlon[:, ::-1]
2523 else:
2524 return num.concatenate(
2525 (latlon, points[:, 2].reshape((len(points), 1))),
2526 axis=1)
2528 def geometry(self, *args, **kwargs):
2529 from pyrocko.model import Geometry
2531 ds = self.discretize_basesource(*args, **kwargs)
2532 nx, ny = ds.nl, ds.nw
2534 def patch_outlines_xy(nx, ny):
2535 points = num.zeros((nx * ny, 2))
2536 points[:, 0] = num.tile(num.linspace(-1., 1., nx), ny)
2537 points[:, 1] = num.repeat(num.linspace(-1., 1., ny), nx)
2539 return points
2541 points_ds = patch_outlines_xy(nx + 1, ny + 1)
2542 npoints = (nx + 1) * (ny + 1)
2544 vertices = num.hstack((
2545 num.ones((npoints, 1)) * self.lat,
2546 num.ones((npoints, 1)) * self.lon,
2547 self.xy_to_coord(points_ds[:, 0], points_ds[:, 1], cs='xyz')))
2549 faces = num.array([[
2550 iy * (nx + 1) + ix,
2551 iy * (nx + 1) + ix + 1,
2552 (iy + 1) * (nx + 1) + ix + 1,
2553 (iy + 1) * (nx + 1) + ix,
2554 iy * (nx + 1) + ix]
2555 for iy in range(ny) for ix in range(nx)])
2557 xyz = self.outline('xyz')
2558 latlon = num.ones((5, 2)) * num.array([self.lat, self.lon])
2559 patchverts = num.hstack((latlon, xyz))
2561 geom = Geometry()
2562 geom.setup(vertices, faces)
2563 geom.set_outlines([patchverts])
2565 if self.stf:
2566 geom.times = num.unique(ds.times)
2568 if self.nucleation_x is not None and self.nucleation_y is not None:
2569 geom.add_property('t_arrival', ds.times)
2571 geom.add_property(
2572 'moment', ds.moments().reshape(ds.nl*ds.nw, -1))
2574 geom.add_property(
2575 'slip', num.ones_like(ds.times) * self.slip)
2577 return geom
2579 def get_nucleation_abs_coord(self, cs='xy'):
2581 if self.nucleation_x is None:
2582 return None, None
2584 coords = from_plane_coords(self.strike, self.dip, self.length,
2585 self.width, self.depth, self.nucleation_x,
2586 self.nucleation_y, lat=self.lat,
2587 lon=self.lon, north_shift=self.north_shift,
2588 east_shift=self.east_shift, cs=cs)
2589 return coords
2591 def pyrocko_moment_tensor(self, store=None, target=None):
2592 return pmt.MomentTensor(
2593 strike=self.strike,
2594 dip=self.dip,
2595 rake=self.rake,
2596 scalar_moment=self.get_moment(store, target))
2598 def pyrocko_event(self, store=None, target=None, **kwargs):
2599 return SourceWithDerivedMagnitude.pyrocko_event(
2600 self, store, target,
2601 **kwargs)
2603 @classmethod
2604 def from_pyrocko_event(cls, ev, **kwargs):
2605 d = {}
2606 mt = ev.moment_tensor
2607 if mt:
2608 (strike, dip, rake), _ = mt.both_strike_dip_rake()
2609 d.update(
2610 strike=float(strike),
2611 dip=float(dip),
2612 rake=float(rake),
2613 magnitude=float(mt.moment_magnitude()))
2615 d.update(kwargs)
2616 return super(RectangularSource, cls).from_pyrocko_event(ev, **d)
2619class PseudoDynamicRupture(SourceWithDerivedMagnitude):
2620 '''
2621 Combined Eikonal and Okada quasi-dynamic rupture model.
2623 Details are described in :doc:`/topics/pseudo-dynamic-rupture`.
2624 Note: attribute `stf` is not used so far, but kept for future applications.
2625 '''
2627 discretized_source_class = meta.DiscretizedMTSource
2629 strike = Float.T(
2630 default=0.0,
2631 help='Strike direction in [deg], measured clockwise from north.')
2633 dip = Float.T(
2634 default=0.0,
2635 help='Dip angle in [deg], measured downward from horizontal.')
2637 length = Float.T(
2638 default=10. * km,
2639 help='Length of rectangular source area in [m].')
2641 width = Float.T(
2642 default=5. * km,
2643 help='Width of rectangular source area in [m].')
2645 anchor = StringChoice.T(
2646 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2647 'bottom_left', 'bottom_right'],
2648 default='center',
2649 optional=True,
2650 help='Anchor point for positioning the plane, can be: ``top, center, '
2651 'bottom, top_left, top_right, bottom_left, '
2652 'bottom_right, center_left, center_right``.')
2654 nucleation_x__ = Array.T(
2655 default=num.array([0.]),
2656 dtype=num.float64,
2657 serialize_as='list',
2658 help='Horizontal position of rupture nucleation in normalized fault '
2659 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge).')
2661 nucleation_y__ = Array.T(
2662 default=num.array([0.]),
2663 dtype=num.float64,
2664 serialize_as='list',
2665 help='Down-dip position of rupture nucleation in normalized fault '
2666 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge).')
2668 nucleation_time__ = Array.T(
2669 optional=True,
2670 help='Time in [s] after origin, when nucleation points defined by '
2671 '``nucleation_x`` and ``nucleation_y`` rupture.',
2672 dtype=num.float64,
2673 serialize_as='list')
2675 gamma = Float.T(
2676 default=0.8,
2677 help='Scaling factor between rupture velocity and S-wave velocity: '
2678 r':math:`v_r = \gamma * v_s`.')
2680 nx = Int.T(
2681 default=2,
2682 help='Number of discrete source patches in x direction (along '
2683 'strike).')
2685 ny = Int.T(
2686 default=2,
2687 help='Number of discrete source patches in y direction (down dip).')
2689 slip = Float.T(
2690 optional=True,
2691 help='Maximum slip of the rectangular source [m]. '
2692 'Setting the slip the tractions/stress field '
2693 'will be normalized to accomodate the desired maximum slip.')
2695 rake = Float.T(
2696 optional=True,
2697 help='Rake angle in [deg], '
2698 'measured counter-clockwise from right-horizontal '
2699 'in on-plane view. Rake is translated into homogenous tractions '
2700 'in strike and up-dip direction. ``rake`` is mutually exclusive '
2701 'with tractions parameter.')
2703 patches = List.T(
2704 OkadaSource.T(),
2705 optional=True,
2706 help='List of all boundary elements/sub faults/fault patches.')
2708 patch_mask__ = Array.T(
2709 dtype=bool,
2710 serialize_as='list',
2711 shape=(None,),
2712 optional=True,
2713 help='Mask for all boundary elements/sub faults/fault patches. True '
2714 'leaves the patch in the calculation, False excludes the patch.')
2716 tractions = TractionField.T(
2717 optional=True,
2718 help='Traction field the rupture plane is exposed to. See the '
2719 ':py:mod:`pyrocko.gf.tractions` module for more details. '
2720 'If ``tractions=None`` and ``rake`` is given'
2721 ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will'
2722 ' be used.')
2724 coef_mat = Array.T(
2725 optional=True,
2726 help='Coefficient matrix linking traction and dislocation field.',
2727 dtype=num.float64,
2728 shape=(None, None))
2730 eikonal_decimation = Int.T(
2731 optional=True,
2732 default=1,
2733 help='Sub-source eikonal factor, a smaller eikonal factor will'
2734 ' increase the accuracy of rupture front calculation but'
2735 ' increases also the computation time.')
2737 decimation_factor = Int.T(
2738 optional=True,
2739 default=1,
2740 help='Sub-source decimation factor, a larger decimation will'
2741 ' make the result inaccurate but shorten the necessary'
2742 ' computation time (use for testing puposes only).')
2744 nthreads = Int.T(
2745 optional=True,
2746 default=1,
2747 help='Number of threads for Okada forward modelling, '
2748 'matrix inversion and calculation of point subsources. '
2749 'Note: for small/medium matrices 1 thread is most efficient.')
2751 pure_shear = Bool.T(
2752 optional=True,
2753 default=False,
2754 help='Calculate only shear tractions and omit tensile tractions.')
2756 smooth_rupture = Bool.T(
2757 default=True,
2758 help='Smooth the tractions by weighting partially ruptured'
2759 ' fault patches.')
2761 aggressive_oversampling = Bool.T(
2762 default=False,
2763 help='Aggressive oversampling for basesource discretization. '
2764 "When using 'multilinear' interpolation oversampling has"
2765 ' practically no effect.')
2767 def __init__(self, **kwargs):
2768 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2769 self._interpolators = {}
2770 self.check_conflicts()
2772 @property
2773 def nucleation_x(self):
2774 return self.nucleation_x__
2776 @nucleation_x.setter
2777 def nucleation_x(self, nucleation_x):
2778 if isinstance(nucleation_x, list):
2779 nucleation_x = num.array(nucleation_x)
2781 elif not isinstance(
2782 nucleation_x, num.ndarray) and nucleation_x is not None:
2784 nucleation_x = num.array([nucleation_x])
2785 self.nucleation_x__ = nucleation_x
2787 @property
2788 def nucleation_y(self):
2789 return self.nucleation_y__
2791 @nucleation_y.setter
2792 def nucleation_y(self, nucleation_y):
2793 if isinstance(nucleation_y, list):
2794 nucleation_y = num.array(nucleation_y)
2796 elif not isinstance(nucleation_y, num.ndarray) \
2797 and nucleation_y is not None:
2798 nucleation_y = num.array([nucleation_y])
2800 self.nucleation_y__ = nucleation_y
2802 @property
2803 def nucleation(self):
2804 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
2806 if (nucl_x is None) or (nucl_y is None):
2807 return None
2809 assert nucl_x.shape[0] == nucl_y.shape[0]
2811 return num.concatenate(
2812 (nucl_x[:, num.newaxis], nucl_y[:, num.newaxis]), axis=1)
2814 @nucleation.setter
2815 def nucleation(self, nucleation):
2816 if isinstance(nucleation, list):
2817 nucleation = num.array(nucleation)
2819 assert nucleation.shape[1] == 2
2821 self.nucleation_x = nucleation[:, 0]
2822 self.nucleation_y = nucleation[:, 1]
2824 @property
2825 def nucleation_time(self):
2826 return self.nucleation_time__
2828 @nucleation_time.setter
2829 def nucleation_time(self, nucleation_time):
2830 if not isinstance(nucleation_time, num.ndarray) \
2831 and nucleation_time is not None:
2832 nucleation_time = num.array([nucleation_time])
2834 self.nucleation_time__ = nucleation_time
2836 @property
2837 def patch_mask(self):
2838 if (self.patch_mask__ is not None and
2839 self.patch_mask__.shape == (self.nx * self.ny,)):
2841 return self.patch_mask__
2842 else:
2843 return num.ones(self.nx * self.ny, dtype=bool)
2845 @patch_mask.setter
2846 def patch_mask(self, patch_mask):
2847 if isinstance(patch_mask, list):
2848 patch_mask = num.array(patch_mask)
2850 self.patch_mask__ = patch_mask
2852 def get_tractions(self):
2853 '''
2854 Get source traction vectors.
2856 If :py:attr:`rake` is given, unit length directed traction vectors
2857 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are returned,
2858 else the given :py:attr:`tractions` are used.
2860 :returns:
2861 Traction vectors per patch.
2862 :rtype:
2863 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2864 '''
2866 if self.rake is not None:
2867 if num.isnan(self.rake):
2868 raise ValueError('Rake must be a real number, not NaN.')
2870 logger.warning(
2871 'Tractions are derived based on the given source rake.')
2872 tractions = DirectedTractions(rake=self.rake)
2873 else:
2874 tractions = self.tractions
2875 return tractions.get_tractions(self.nx, self.ny, self.patches)
2877 def get_scaled_tractions(self, store):
2878 '''
2879 Get traction vectors rescaled to given slip.
2881 Opposing to :py:meth:`get_tractions` traction vectors
2882 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are rescaled to
2883 the given :py:attr:`slip` before returning. If no :py:attr:`slip` and
2884 :py:attr:`rake` are provided, the given :py:attr:`tractions` are
2885 returned without scaling.
2887 :param store:
2888 Green's function database (needs to cover whole region of the
2889 source).
2890 :type store:
2891 :py:class:`~pyrocko.gf.store.Store`
2893 :returns:
2894 Traction vectors per patch.
2895 :rtype:
2896 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2897 '''
2898 tractions = self.tractions
2899 factor = 1.
2901 if self.rake is not None and self.slip is not None:
2902 if num.isnan(self.rake):
2903 raise ValueError('Rake must be a real number, not NaN.')
2905 self.discretize_patches(store)
2906 slip_0t = max(num.linalg.norm(
2907 self.get_slip(scale_slip=False),
2908 axis=1))
2910 factor = self.slip / slip_0t
2911 tractions = DirectedTractions(rake=self.rake)
2913 return tractions.get_tractions(self.nx, self.ny, self.patches) * factor
2915 def base_key(self):
2916 return SourceWithDerivedMagnitude.base_key(self) + (
2917 self.slip,
2918 self.strike,
2919 self.dip,
2920 self.rake,
2921 self.length,
2922 self.width,
2923 float(self.nucleation_x.mean()),
2924 float(self.nucleation_y.mean()),
2925 self.decimation_factor,
2926 self.anchor,
2927 self.pure_shear,
2928 self.gamma,
2929 tuple(self.patch_mask))
2931 def check_conflicts(self):
2932 if self.tractions and self.rake:
2933 raise AttributeError(
2934 'Tractions and rake are mutually exclusive.')
2935 if self.tractions is None and self.rake is None:
2936 self.rake = 0.
2938 def get_magnitude(self, store=None, target=None):
2939 '''
2940 Get total seismic moment magnitude Mw.
2942 :param store:
2943 GF store to guide the discretization and providing the earthmodel
2944 which is needed to calculate moment from slip.
2945 :type store:
2946 :py:class:`~pyrocko.gf.store.Store`
2948 :param target:
2949 Target, used to get GF interpolation settings.
2950 :type target:
2951 :py:class:`pyrocko.gf.targets.Target`
2953 :returns:
2954 Moment magnitude
2955 :rtype:
2956 float
2957 '''
2958 self.check_conflicts()
2959 if self.slip is not None or self.tractions is not None:
2960 if store is None:
2961 raise DerivedMagnitudeError(
2962 'Magnitude for a rectangular source with slip or '
2963 'tractions defined can only be derived when earth model '
2964 'is set.')
2966 moment_rate, calc_times = self.discretize_basesource(
2967 store, target=target).get_moment_rate(store.config.deltat)
2969 deltat = num.concatenate((
2970 (num.diff(calc_times)[0],),
2971 num.diff(calc_times)))
2973 return float(pmt.moment_to_magnitude(
2974 num.sum(moment_rate * deltat)))
2976 else:
2977 return float(pmt.moment_to_magnitude(1.0))
2979 def get_factor(self):
2980 return 1.0
2982 def outline(self, cs='xyz'):
2983 '''
2984 Get source outline corner coordinates.
2986 :param cs:
2987 :ref:`Output coordinate system <coordinate-system-names>`.
2988 :type cs:
2989 str
2991 :returns:
2992 Corner points in desired coordinate system.
2993 :rtype:
2994 :py:class:`~numpy.ndarray`: ``(5, [2, 3])``.
2995 '''
2996 points = outline_rect_source(self.strike, self.dip, self.length,
2997 self.width, self.anchor)
2999 points[:, 0] += self.north_shift
3000 points[:, 1] += self.east_shift
3001 points[:, 2] += self.depth
3002 if cs == 'xyz':
3003 return points
3004 elif cs == 'xy':
3005 return points[:, :2]
3006 elif cs in ('latlon', 'lonlat', 'latlondepth'):
3007 latlon = ne_to_latlon(
3008 self.lat, self.lon, points[:, 0], points[:, 1])
3010 latlon = num.array(latlon).T
3011 if cs == 'latlon':
3012 return latlon
3013 elif cs == 'lonlat':
3014 return latlon[:, ::-1]
3015 else:
3016 return num.concatenate(
3017 (latlon, points[:, 2].reshape((len(points), 1))),
3018 axis=1)
3020 def points_on_source(self, cs='xyz', **kwargs):
3021 '''
3022 Convert relative plane coordinates to geographical coordinates.
3024 Given x and y coordinates (relative source coordinates between -1.
3025 and 1.) are converted to desired geographical coordinates. Coordinates
3026 need to be given as :py:class:`~numpy.ndarray` arguments ``points_x``
3027 and ``points_y``.
3029 :param cs:
3030 :ref:`Output coordinate system <coordinate-system-names>`.
3031 :type cs:
3032 str
3034 :returns:
3035 Point coordinates in desired coordinate system.
3036 :rtype:
3037 :py:class:`~numpy.ndarray`: ``(n_points, [2, 3])``.
3038 '''
3039 points = points_on_rect_source(
3040 self.strike, self.dip, self.length, self.width,
3041 self.anchor, **kwargs)
3043 points[:, 0] += self.north_shift
3044 points[:, 1] += self.east_shift
3045 points[:, 2] += self.depth
3046 if cs == 'xyz':
3047 return points
3048 elif cs == 'xy':
3049 return points[:, :2]
3050 elif cs in ('latlon', 'lonlat', 'latlondepth'):
3051 latlon = ne_to_latlon(
3052 self.lat, self.lon, points[:, 0], points[:, 1])
3054 latlon = num.array(latlon).T
3055 if cs == 'latlon':
3056 return latlon
3057 elif cs == 'lonlat':
3058 return latlon[:, ::-1]
3059 else:
3060 return num.concatenate(
3061 (latlon, points[:, 2].reshape((len(points), 1))),
3062 axis=1)
3064 def pyrocko_moment_tensor(self, store=None, target=None):
3065 '''
3066 Get overall moment tensor of the rupture.
3068 :param store:
3069 GF store to guide the discretization and providing the earthmodel
3070 which is needed to calculate moment from slip.
3071 :type store:
3072 :py:class:`~pyrocko.gf.store.Store`
3074 :param target:
3075 Target, used to get GF interpolation settings.
3076 :type target:
3077 :py:class:`pyrocko.gf.targets.Target`
3079 :returns:
3080 Moment tensor.
3081 :rtype:
3082 :py:class:`~pyrocko.moment_tensor.MomentTensor`
3083 '''
3084 if store is not None:
3085 if not self.patches:
3086 self.discretize_patches(store)
3088 data = self.get_slip()
3089 else:
3090 data = self.get_tractions()
3092 weights = num.linalg.norm(data, axis=1)
3093 weights /= weights.sum()
3095 rakes = num.arctan2(data[:, 1], data[:, 0]) * r2d
3096 rake = num.average(rakes, weights=weights)
3098 return pmt.MomentTensor(
3099 strike=self.strike,
3100 dip=self.dip,
3101 rake=rake,
3102 scalar_moment=self.get_moment(store, target))
3104 def pyrocko_event(self, store=None, target=None, **kwargs):
3105 return SourceWithDerivedMagnitude.pyrocko_event(
3106 self, store, target,
3107 **kwargs)
3109 @classmethod
3110 def from_pyrocko_event(cls, ev, **kwargs):
3111 d = {}
3112 mt = ev.moment_tensor
3113 if mt:
3114 (strike, dip, rake), _ = mt.both_strike_dip_rake()
3115 d.update(
3116 strike=float(strike),
3117 dip=float(dip),
3118 rake=float(rake))
3120 d.update(kwargs)
3121 return super(PseudoDynamicRupture, cls).from_pyrocko_event(ev, **d)
3123 def _discretize_points(self, store, *args, **kwargs):
3124 '''
3125 Discretize source plane with equal vertical and horizontal spacing.
3127 Additional ``*args`` and ``**kwargs`` are passed to
3128 :py:meth:`points_on_source`.
3130 :param store:
3131 Green's function database (needs to cover whole region of the
3132 source).
3133 :type store:
3134 :py:class:`~pyrocko.gf.store.Store`
3136 :returns:
3137 Number of points in strike and dip direction, distance
3138 between adjacent points, coordinates (latlondepth) and coordinates
3139 (xy on fault) for discrete points.
3140 :rtype:
3141 (int, int, float, :py:class:`~numpy.ndarray`,
3142 :py:class:`~numpy.ndarray`).
3143 '''
3144 anch_x, anch_y = map_anchor[self.anchor]
3146 npoints = int(self.width // km) + 1
3147 points = num.zeros((npoints, 3))
3148 points[:, 1] = num.linspace(-1., 1., npoints)
3149 points[:, 1] = (points[:, 1] - anch_y) * self.width/2
3151 rotmat = pmt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.0)
3152 points = num.dot(rotmat.T, points.T).T
3153 points[:, 2] += self.depth
3155 vs_min = store.config.get_vs(
3156 self.lat, self.lon, points,
3157 interpolation='nearest_neighbor')
3158 vr_min = max(vs_min.min(), .5*km) * self.gamma
3160 oversampling = 10.
3161 delta_l = self.length / (self.nx * oversampling)
3162 delta_w = self.width / (self.ny * oversampling)
3164 delta = self.eikonal_decimation * num.min([
3165 store.config.deltat * vr_min / oversampling,
3166 delta_l, delta_w] + [
3167 deltas for deltas in store.config.deltas])
3169 delta = delta_w / num.ceil(delta_w / delta)
3171 nx = int(num.ceil(self.length / delta)) + 1
3172 ny = int(num.ceil(self.width / delta)) + 1
3174 rem_l = (nx-1)*delta - self.length
3175 lim_x = rem_l / self.length
3177 points_xy = num.zeros((nx * ny, 2))
3178 points_xy[:, 0] = num.repeat(
3179 num.linspace(-1.-lim_x, 1.+lim_x, nx), ny)
3180 points_xy[:, 1] = num.tile(
3181 num.linspace(-1., 1., ny), nx)
3183 points = self.points_on_source(
3184 points_x=points_xy[:, 0],
3185 points_y=points_xy[:, 1],
3186 **kwargs)
3188 return nx, ny, delta, points, points_xy
3190 def _discretize_rupture_v(self, store, interpolation='nearest_neighbor',
3191 points=None):
3192 '''
3193 Get rupture velocity for discrete points on source plane.
3195 :param store:
3196 Green's function database (needs to cover the whole region of the
3197 source)
3198 :type store:
3199 :py:class:`~pyrocko.gf.store.Store`
3201 :param interpolation:
3202 Interpolation method to use (choose between ``'nearest_neighbor'``
3203 and ``'multilinear'``).
3204 :type interpolation:
3205 str
3207 :param points:
3208 Coordinates on fault (-1.:1.) of discrete points.
3209 :type points:
3210 :py:class:`~numpy.ndarray`: ``(n_points, 2)``
3212 :returns:
3213 Rupture velocity assumed as :math:`v_s * \\gamma` for discrete
3214 points.
3215 :rtype:
3216 :py:class:`~numpy.ndarray`: ``(n_points, )``.
3217 '''
3219 if points is None:
3220 _, _, _, points, _ = self._discretize_points(store, cs='xyz')
3222 return store.config.get_vs(
3223 self.lat, self.lon,
3224 points=points,
3225 interpolation=interpolation) * self.gamma
3227 def discretize_time(
3228 self, store, interpolation='nearest_neighbor',
3229 vr=None, times=None, *args, **kwargs):
3230 '''
3231 Get rupture start time for discrete points on source plane.
3233 :param store:
3234 Green's function database (needs to cover whole region of the
3235 source)
3236 :type store:
3237 :py:class:`~pyrocko.gf.store.Store`
3239 :param interpolation:
3240 Interpolation method to use (choose between ``'nearest_neighbor'``
3241 and ``'multilinear'``).
3242 :type interpolation:
3243 str
3245 :param vr:
3246 Array, containing rupture user defined rupture velocity values.
3247 :type vr:
3248 :py:class:`~numpy.ndarray`
3250 :param times:
3251 Array, containing zeros, where rupture is starting, real positive
3252 numbers at later secondary nucleation points and -1, where time
3253 will be calculated. If not given, rupture starts at nucleation_x,
3254 nucleation_y. Times are given for discrete points with equal
3255 horizontal and vertical spacing.
3256 :type times:
3257 :py:class:`~numpy.ndarray`
3259 :returns:
3260 Coordinates (latlondepth), coordinates (xy), rupture velocity,
3261 rupture propagation time of discrete points.
3262 :rtype:
3263 :py:class:`~numpy.ndarray`: ``(n_points, 3)``,
3264 :py:class:`~numpy.ndarray`: ``(n_points, 2)``,
3265 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``,
3266 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``.
3267 '''
3268 nx, ny, delta, points, points_xy = self._discretize_points(
3269 store, cs='xyz')
3271 if vr is None or vr.shape != tuple((nx, ny)):
3272 if vr:
3273 logger.warning(
3274 'Given rupture velocities are not in right shape: '
3275 '(%i, %i), but needed is (%i, %i).', *vr.shape + (nx, ny))
3276 vr = self._discretize_rupture_v(store, interpolation, points)\
3277 .reshape(nx, ny)
3279 if vr.shape != tuple((nx, ny)):
3280 logger.warning(
3281 'Given rupture velocities are not in right shape. Therefore'
3282 ' standard rupture velocity array is used.')
3284 def initialize_times():
3285 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
3287 if nucl_x.shape != nucl_y.shape:
3288 raise ValueError(
3289 'Nucleation coordinates have different shape.')
3291 dist_points = num.array([
3292 num.linalg.norm(points_xy - num.array([x, y]).ravel(), axis=1)
3293 for x, y in zip(nucl_x, nucl_y)])
3294 nucl_indices = num.argmin(dist_points, axis=1)
3296 if self.nucleation_time is None:
3297 nucl_times = num.zeros_like(nucl_indices)
3298 else:
3299 if self.nucleation_time.shape == nucl_x.shape:
3300 nucl_times = self.nucleation_time
3301 else:
3302 raise ValueError(
3303 'Nucleation coordinates and times have different '
3304 'shapes')
3306 t = num.full(nx * ny, -1.)
3307 t[nucl_indices] = nucl_times
3308 return t.reshape(nx, ny)
3310 if times is None:
3311 times = initialize_times()
3312 elif times.shape != tuple((nx, ny)):
3313 times = initialize_times()
3314 logger.warning(
3315 'Given times are not in right shape. Therefore standard time'
3316 ' array is used.')
3318 eikonal_ext.eikonal_solver_fmm_cartesian(
3319 speeds=vr, times=times, delta=delta)
3321 return points, points_xy, vr, times
3323 def get_vr_time_interpolators(
3324 self, store, interpolation='nearest_neighbor', force=False,
3325 **kwargs):
3326 '''
3327 Get interpolators for rupture velocity and rupture time.
3329 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3331 :param store:
3332 Green's function database (needs to cover whole region of the
3333 source).
3334 :type store:
3335 :py:class:`~pyrocko.gf.store.Store`
3337 :param interpolation:
3338 Interpolation method to use (choose between ``'nearest_neighbor'``
3339 and ``'multilinear'``).
3340 :type interpolation:
3341 str
3343 :param force:
3344 Force recalculation of the interpolators (e.g. after change of
3345 nucleation point locations/times). Default is ``False``.
3346 :type force:
3347 bool
3348 '''
3349 interp_map = {'multilinear': 'linear', 'nearest_neighbor': 'nearest'}
3350 if interpolation not in interp_map:
3351 raise TypeError(
3352 'Interpolation method %s not available' % interpolation)
3354 if not self._interpolators.get(interpolation, False) or force:
3355 _, points_xy, vr, times = self.discretize_time(
3356 store, **kwargs)
3358 if self.length <= 0.:
3359 raise ValueError(
3360 'length must be larger then 0. not %g' % self.length)
3362 if self.width <= 0.:
3363 raise ValueError(
3364 'width must be larger then 0. not %g' % self.width)
3366 nx, ny = times.shape
3367 anch_x, anch_y = map_anchor[self.anchor]
3369 points_xy[:, 0] = (points_xy[:, 0] - anch_x) * self.length / 2.
3370 points_xy[:, 1] = (points_xy[:, 1] - anch_y) * self.width / 2.
3372 ascont = num.ascontiguousarray
3374 self._interpolators[interpolation] = (
3375 nx, ny, times, vr,
3376 RegularGridInterpolator(
3377 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3378 times,
3379 method=interp_map[interpolation]),
3380 RegularGridInterpolator(
3381 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3382 vr,
3383 method=interp_map[interpolation]))
3385 return self._interpolators[interpolation]
3387 def discretize_patches(
3388 self, store, interpolation='nearest_neighbor', force=False,
3389 grid_shape=(),
3390 **kwargs):
3391 '''
3392 Get rupture start time and OkadaSource elements for points on rupture.
3394 All source elements and their corresponding center points are
3395 calculated and stored in the :py:attr:`patches` attribute.
3397 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3399 :param store:
3400 Green's function database (needs to cover whole region of the
3401 source).
3402 :type store:
3403 :py:class:`~pyrocko.gf.store.Store`
3405 :param interpolation:
3406 Interpolation method to use (choose between ``'nearest_neighbor'``
3407 and ``'multilinear'``).
3408 :type interpolation:
3409 str
3411 :param force:
3412 Force recalculation of the vr and time interpolators ( e.g. after
3413 change of nucleation point locations/times). Default is ``False``.
3414 :type force:
3415 bool
3417 :param grid_shape:
3418 Desired sub fault patch grid size (nlength, nwidth). Either factor
3419 or grid_shape should be set.
3420 :type grid_shape:
3421 :py:class:`tuple` of :py:class:`int`
3422 '''
3423 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3424 self.get_vr_time_interpolators(
3425 store,
3426 interpolation=interpolation, force=force, **kwargs)
3427 anch_x, anch_y = map_anchor[self.anchor]
3429 al = self.length / 2.
3430 aw = self.width / 2.
3431 al1 = -(al + anch_x * al)
3432 al2 = al - anch_x * al
3433 aw1 = -aw + anch_y * aw
3434 aw2 = aw + anch_y * aw
3435 assert num.abs([al1, al2]).sum() == self.length
3436 assert num.abs([aw1, aw2]).sum() == self.width
3438 def get_lame(*a, **kw):
3439 shear_mod = store.config.get_shear_moduli(*a, **kw)
3440 lamb = store.config.get_vp(*a, **kw)**2 \
3441 * store.config.get_rho(*a, **kw) - 2. * shear_mod
3442 return shear_mod, lamb / (2. * (lamb + shear_mod))
3444 shear_mod, poisson = get_lame(
3445 self.lat, self.lon,
3446 num.array([[self.north_shift, self.east_shift, self.depth]]),
3447 interpolation=interpolation)
3449 okada_src = OkadaSource(
3450 lat=self.lat, lon=self.lon,
3451 strike=self.strike, dip=self.dip,
3452 north_shift=self.north_shift, east_shift=self.east_shift,
3453 depth=self.depth,
3454 al1=al1, al2=al2, aw1=aw1, aw2=aw2,
3455 poisson=poisson.mean(),
3456 shearmod=shear_mod.mean(),
3457 opening=kwargs.get('opening', 0.))
3459 if not (self.nx and self.ny):
3460 if grid_shape:
3461 self.nx, self.ny = grid_shape
3462 else:
3463 self.nx = nx
3464 self.ny = ny
3466 source_disc, source_points = okada_src.discretize(self.nx, self.ny)
3468 shear_mod, poisson = get_lame(
3469 self.lat, self.lon,
3470 num.array([src.source_patch()[:3] for src in source_disc]),
3471 interpolation=interpolation)
3473 if (self.nx, self.ny) != (nx, ny):
3474 times_interp = time_interpolator(
3475 num.ascontiguousarray(source_points[:, :2]))
3476 vr_interp = vr_interpolator(
3477 num.ascontiguousarray(source_points[:, :2]))
3478 else:
3479 times_interp = times.T.ravel()
3480 vr_interp = vr.T.ravel()
3482 for isrc, src in enumerate(source_disc):
3483 src.vr = vr_interp[isrc]
3484 src.time = times_interp[isrc] + self.time
3486 self.patches = source_disc
3488 def discretize_basesource(self, store, target=None):
3489 '''
3490 Prepare source for synthetic waveform calculation.
3492 :param store:
3493 Green's function database (needs to cover whole region of the
3494 source).
3495 :type store:
3496 :py:class:`~pyrocko.gf.store.Store`
3498 :param target:
3499 Target information.
3500 :type target:
3501 :py:class:`~pyrocko.gf.targets.Target`
3503 :returns:
3504 Source discretized by a set of moment tensors and times.
3505 :rtype:
3506 :py:class:`~pyrocko.gf.meta.DiscretizedMTSource`
3507 '''
3508 if not target:
3509 interpolation = 'nearest_neighbor'
3510 else:
3511 interpolation = target.interpolation
3513 if not self.patches:
3514 self.discretize_patches(store, interpolation)
3516 if self.coef_mat is None:
3517 self.calc_coef_mat()
3519 delta_slip, slip_times = self.get_delta_slip(store)
3520 npatches = self.nx * self.ny
3521 ntimes = slip_times.size
3523 anch_x, anch_y = map_anchor[self.anchor]
3525 pln = self.length / self.nx
3526 pwd = self.width / self.ny
3528 patch_coords = num.array([
3529 (p.ix, p.iy)
3530 for p in self.patches]).reshape(self.nx, self.ny, 2)
3532 # boundary condition is zero-slip
3533 # is not valid to avoid unwished interpolation effects
3534 slip_grid = num.zeros((self.nx + 2, self.ny + 2, ntimes, 3))
3535 slip_grid[1:-1, 1:-1, :, :] = \
3536 delta_slip.reshape(self.nx, self.ny, ntimes, 3)
3538 slip_grid[0, 0, :, :] = slip_grid[1, 1, :, :]
3539 slip_grid[0, -1, :, :] = slip_grid[1, -2, :, :]
3540 slip_grid[-1, 0, :, :] = slip_grid[-2, 1, :, :]
3541 slip_grid[-1, -1, :, :] = slip_grid[-2, -2, :, :]
3543 slip_grid[1:-1, 0, :, :] = slip_grid[1:-1, 1, :, :]
3544 slip_grid[1:-1, -1, :, :] = slip_grid[1:-1, -2, :, :]
3545 slip_grid[0, 1:-1, :, :] = slip_grid[1, 1:-1, :, :]
3546 slip_grid[-1, 1:-1, :, :] = slip_grid[-2, 1:-1, :, :]
3548 def make_grid(patch_parameter):
3549 grid = num.zeros((self.nx + 2, self.ny + 2))
3550 grid[1:-1, 1:-1] = patch_parameter.reshape(self.nx, self.ny)
3552 grid[0, 0] = grid[1, 1]
3553 grid[0, -1] = grid[1, -2]
3554 grid[-1, 0] = grid[-2, 1]
3555 grid[-1, -1] = grid[-2, -2]
3557 grid[1:-1, 0] = grid[1:-1, 1]
3558 grid[1:-1, -1] = grid[1:-1, -2]
3559 grid[0, 1:-1] = grid[1, 1:-1]
3560 grid[-1, 1:-1] = grid[-2, 1:-1]
3562 return grid
3564 lamb = self.get_patch_attribute('lamb')
3565 mu = self.get_patch_attribute('shearmod')
3567 lamb_grid = make_grid(lamb)
3568 mu_grid = make_grid(mu)
3570 coords_x = num.zeros(self.nx + 2)
3571 coords_x[1:-1] = patch_coords[:, 0, 0]
3572 coords_x[0] = coords_x[1] - pln / 2
3573 coords_x[-1] = coords_x[-2] + pln / 2
3575 coords_y = num.zeros(self.ny + 2)
3576 coords_y[1:-1] = patch_coords[0, :, 1]
3577 coords_y[0] = coords_y[1] - pwd / 2
3578 coords_y[-1] = coords_y[-2] + pwd / 2
3580 slip_interp = RegularGridInterpolator(
3581 (coords_x, coords_y, slip_times),
3582 slip_grid, method='nearest')
3584 lamb_interp = RegularGridInterpolator(
3585 (coords_x, coords_y),
3586 lamb_grid, method='nearest')
3588 mu_interp = RegularGridInterpolator(
3589 (coords_x, coords_y),
3590 mu_grid, method='nearest')
3592 # discretize basesources
3593 mindeltagf = min(tuple(
3594 (self.length / self.nx, self.width / self.ny) +
3595 tuple(store.config.deltas)))
3597 nl = int((1. / self.decimation_factor) *
3598 num.ceil(pln / mindeltagf)) + 1
3599 nw = int((1. / self.decimation_factor) *
3600 num.ceil(pwd / mindeltagf)) + 1
3601 nsrc_patch = int(nl * nw)
3602 dl = pln / nl
3603 dw = pwd / nw
3605 patch_area = dl * dw
3607 xl = num.linspace(-0.5 * (pln - dl), 0.5 * (pln - dl), nl)
3608 xw = num.linspace(-0.5 * (pwd - dw), 0.5 * (pwd - dw), nw)
3610 base_coords = num.zeros((nsrc_patch, 3))
3611 base_coords[:, 0] = num.tile(xl, nw)
3612 base_coords[:, 1] = num.repeat(xw, nl)
3613 base_coords = num.tile(base_coords, (npatches, 1))
3615 center_coords = num.zeros((npatches, 3))
3616 center_coords[:, 0] = num.repeat(
3617 num.arange(self.nx) * pln + pln / 2, self.ny) - self.length / 2
3618 center_coords[:, 1] = num.tile(
3619 num.arange(self.ny) * pwd + pwd / 2, self.nx) - self.width / 2
3621 base_coords -= center_coords.repeat(nsrc_patch, axis=0)
3622 nbaselocs = base_coords.shape[0]
3624 base_interp = base_coords.repeat(ntimes, axis=0)
3626 base_times = num.tile(slip_times, nbaselocs)
3627 base_interp[:, 0] -= anch_x * self.length / 2
3628 base_interp[:, 1] -= anch_y * self.width / 2
3629 base_interp[:, 2] = base_times
3631 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3632 store, interpolation=interpolation)
3634 time_eikonal_max = time_interpolator.values.max()
3636 nbasesrcs = base_interp.shape[0]
3637 delta_slip = slip_interp(base_interp).reshape(nbaselocs, ntimes, 3)
3638 lamb = lamb_interp(base_interp[:, :2]).ravel()
3639 mu = mu_interp(base_interp[:, :2]).ravel()
3641 if False:
3642 try:
3643 import matplotlib.pyplot as plt
3644 coords = base_coords.copy()
3645 norm = num.sum(num.linalg.norm(delta_slip, axis=2), axis=1)
3646 plt.scatter(coords[:, 0], coords[:, 1], c=norm)
3647 plt.show()
3648 except AttributeError:
3649 pass
3651 base_interp[:, 2] = 0.
3652 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
3653 base_interp = num.dot(rotmat.T, base_interp.T).T
3654 base_interp[:, 0] += self.north_shift
3655 base_interp[:, 1] += self.east_shift
3656 base_interp[:, 2] += self.depth
3658 slip_strike = delta_slip[:, :, 0].ravel()
3659 slip_dip = delta_slip[:, :, 1].ravel()
3660 slip_norm = delta_slip[:, :, 2].ravel()
3662 slip_shear = num.linalg.norm([slip_strike, slip_dip], axis=0)
3663 slip_rake = r2d * num.arctan2(slip_dip, slip_strike)
3665 m6s = okada_ext.patch2m6(
3666 strikes=num.full(nbasesrcs, self.strike, dtype=float),
3667 dips=num.full(nbasesrcs, self.dip, dtype=float),
3668 rakes=slip_rake,
3669 disl_shear=slip_shear,
3670 disl_norm=slip_norm,
3671 lamb=lamb,
3672 mu=mu,
3673 nthreads=self.nthreads)
3675 m6s *= patch_area
3677 dl = -self.patches[0].al1 + self.patches[0].al2
3678 dw = -self.patches[0].aw1 + self.patches[0].aw2
3680 base_times[base_times > time_eikonal_max] = time_eikonal_max
3682 ds = meta.DiscretizedMTSource(
3683 lat=self.lat,
3684 lon=self.lon,
3685 times=base_times + self.time,
3686 north_shifts=base_interp[:, 0],
3687 east_shifts=base_interp[:, 1],
3688 depths=base_interp[:, 2],
3689 m6s=m6s,
3690 dl=dl,
3691 dw=dw,
3692 nl=self.nx,
3693 nw=self.ny)
3695 return ds
3697 def calc_coef_mat(self):
3698 '''
3699 Calculate coefficients connecting tractions and dislocations.
3700 '''
3701 if not self.patches:
3702 raise ValueError(
3703 'Patches are needed. Please calculate them first.')
3705 self.coef_mat = make_okada_coefficient_matrix(
3706 self.patches, nthreads=self.nthreads, pure_shear=self.pure_shear)
3708 def get_patch_attribute(self, attr):
3709 '''
3710 Get patch attributes.
3712 :param attr:
3713 Name of selected attribute (see
3714 :py:class`pyrocko.modelling.okada.OkadaSource`).
3715 :type attr:
3716 str
3718 :returns:
3719 Array with attribute value for each fault patch.
3720 :rtype:
3721 :py:class:`~numpy.ndarray`
3723 '''
3724 if not self.patches:
3725 raise ValueError(
3726 'Patches are needed. Please calculate them first.')
3727 return num.array([getattr(p, attr) for p in self.patches])
3729 def get_slip(
3730 self,
3731 time=None,
3732 scale_slip=True,
3733 interpolation='nearest_neighbor',
3734 **kwargs):
3735 '''
3736 Get slip per subfault patch for given time after rupture start.
3738 :param time:
3739 Time after origin [s], for which slip is computed. If not
3740 given, final static slip is returned.
3741 :type time:
3742 float
3744 :param scale_slip:
3745 If ``True`` and :py:attr:`slip` given, all slip values are scaled
3746 to fit the given maximum slip.
3747 :type scale_slip:
3748 bool
3750 :param interpolation:
3751 Interpolation method to use (choose between ``'nearest_neighbor'``
3752 and ``'multilinear'``).
3753 :type interpolation:
3754 str
3756 :returns:
3757 Inverted dislocations (:math:`u_{strike}, u_{dip}, u_{tensile}`)
3758 for each source patch.
3759 :rtype:
3760 :py:class:`~numpy.ndarray`: ``(n_sources, 3)``
3761 '''
3763 if self.patches is None:
3764 raise ValueError(
3765 'Please discretize the source first (discretize_patches())')
3766 npatches = len(self.patches)
3767 tractions = self.get_tractions()
3768 time_patch_max = self.get_patch_attribute('time').max() - self.time
3770 time_patch = time
3771 if time is None:
3772 time_patch = time_patch_max
3774 if self.coef_mat is None:
3775 self.calc_coef_mat()
3777 if tractions.shape != (npatches, 3):
3778 raise AttributeError(
3779 'The traction vector is of invalid shape.'
3780 ' Required shape is (npatches, 3)')
3782 patch_mask = num.ones(npatches, dtype=bool)
3783 if self.patch_mask is not None:
3784 patch_mask = self.patch_mask
3786 times = self.get_patch_attribute('time') - self.time
3787 times[~patch_mask] = time_patch + 1. # exlcude unmasked patches
3788 relevant_sources = num.nonzero(times <= time_patch)[0]
3789 disloc_est = num.zeros_like(tractions)
3791 if self.smooth_rupture:
3792 patch_activation = num.zeros(npatches)
3794 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3795 self.get_vr_time_interpolators(
3796 store, interpolation=interpolation)
3798 # Getting the native Eikonal grid, bit hackish
3799 points_x = num.round(time_interpolator.grid[0], decimals=2)
3800 points_y = num.round(time_interpolator.grid[1], decimals=2)
3801 times_eikonal = time_interpolator.values
3803 time_max = time
3804 if time is None:
3805 time_max = times_eikonal.max()
3807 for ip, p in enumerate(self.patches):
3808 ul = num.round((p.ix + p.al1, p.iy + p.aw1), decimals=2)
3809 lr = num.round((p.ix + p.al2, p.iy + p.aw2), decimals=2)
3811 idx_length = num.logical_and(
3812 points_x >= ul[0], points_x <= lr[0])
3813 idx_width = num.logical_and(
3814 points_y >= ul[1], points_y <= lr[1])
3816 times_patch = times_eikonal[num.ix_(idx_length, idx_width)]
3817 if times_patch.size == 0:
3818 raise AttributeError('could not use smooth_rupture')
3820 patch_activation[ip] = \
3821 (times_patch <= time_max).sum() / times_patch.size
3823 if time_patch == 0 and time_patch != time_patch_max:
3824 patch_activation[ip] = 0.
3826 patch_activation[~patch_mask] = 0. # exlcude unmasked patches
3828 relevant_sources = num.nonzero(patch_activation > 0.)[0]
3830 if relevant_sources.size == 0:
3831 return disloc_est
3833 indices_disl = num.repeat(relevant_sources * 3, 3)
3834 indices_disl[1::3] += 1
3835 indices_disl[2::3] += 2
3837 disloc_est[relevant_sources] = invert_fault_dislocations_bem(
3838 stress_field=tractions[relevant_sources, :].ravel(),
3839 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3840 pure_shear=self.pure_shear, nthreads=self.nthreads,
3841 epsilon=None,
3842 **kwargs)
3844 if self.smooth_rupture:
3845 disloc_est *= patch_activation[:, num.newaxis]
3847 if scale_slip and self.slip is not None:
3848 disloc_tmax = num.zeros(npatches)
3850 indices_disl = num.repeat(num.nonzero(patch_mask)[0] * 3, 3)
3851 indices_disl[1::3] += 1
3852 indices_disl[2::3] += 2
3854 disloc_tmax[patch_mask] = num.linalg.norm(
3855 invert_fault_dislocations_bem(
3856 stress_field=tractions[patch_mask, :].ravel(),
3857 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3858 pure_shear=self.pure_shear, nthreads=self.nthreads,
3859 epsilon=None,
3860 **kwargs), axis=1)
3862 disloc_tmax_max = disloc_tmax.max()
3863 if disloc_tmax_max == 0.:
3864 logger.warning(
3865 'slip scaling not performed. Maximum slip is 0.')
3867 disloc_est *= self.slip / disloc_tmax_max
3869 return disloc_est
3871 def get_delta_slip(
3872 self,
3873 store=None,
3874 deltat=None,
3875 delta=True,
3876 interpolation='nearest_neighbor',
3877 **kwargs):
3878 '''
3879 Get slip change snapshots.
3881 The time interval, within which the slip changes are computed is
3882 determined by the sampling rate of the Green's function database or
3883 ``deltat``. Additional ``**kwargs`` are passed to :py:meth:`get_slip`.
3885 :param store:
3886 Green's function database (needs to cover whole region of of the
3887 source). Its sampling interval is used as time increment for slip
3888 difference calculation. Either ``deltat`` or ``store`` should be
3889 given.
3890 :type store:
3891 :py:class:`~pyrocko.gf.store.Store`
3893 :param deltat:
3894 Time interval for slip difference calculation [s]. Either
3895 ``deltat`` or ``store`` should be given.
3896 :type deltat:
3897 float
3899 :param delta:
3900 If ``True``, slip differences between two time steps are given. If
3901 ``False``, cumulative slip for all time steps.
3902 :type delta:
3903 bool
3905 :param interpolation:
3906 Interpolation method to use (choose between ``'nearest_neighbor'``
3907 and ``'multilinear'``).
3908 :type interpolation:
3909 str
3911 :returns:
3912 Displacement changes(:math:`\\Delta u_{strike},
3913 \\Delta u_{dip} , \\Delta u_{tensile}`) for each source patch and
3914 time; corner times, for which delta slip is computed. The order of
3915 displacement changes array is:
3917 .. math::
3919 &[[\\\\
3920 &[\\Delta u_{strike, patch1, t1},
3921 \\Delta u_{dip, patch1, t1},
3922 \\Delta u_{tensile, patch1, t1}],\\\\
3923 &[\\Delta u_{strike, patch1, t2},
3924 \\Delta u_{dip, patch1, t2},
3925 \\Delta u_{tensile, patch1, t2}]\\\\
3926 &], [\\\\
3927 &[\\Delta u_{strike, patch2, t1}, ...],\\\\
3928 &[\\Delta u_{strike, patch2, t2}, ...]]]\\\\
3930 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3931 :py:class:`~numpy.ndarray`: ``(n_times, )``
3932 '''
3933 if store and deltat:
3934 raise AttributeError(
3935 'Argument collision. '
3936 'Please define only the store or the deltat argument.')
3938 if store:
3939 deltat = store.config.deltat
3941 if not deltat:
3942 raise AttributeError('Please give a GF store or set deltat.')
3944 npatches = len(self.patches)
3946 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3947 store, interpolation=interpolation)
3948 tmax = time_interpolator.values.max()
3950 calc_times = num.arange(0., tmax + deltat, deltat)
3951 calc_times[calc_times > tmax] = tmax
3953 disloc_est = num.zeros((npatches, calc_times.size, 3))
3955 for itime, t in enumerate(calc_times):
3956 disloc_est[:, itime, :] = self.get_slip(
3957 time=t, scale_slip=False, **kwargs)
3959 if self.slip:
3960 disloc_tmax = num.linalg.norm(
3961 self.get_slip(scale_slip=False, time=tmax),
3962 axis=1)
3964 disloc_tmax_max = disloc_tmax.max()
3965 if disloc_tmax_max == 0.:
3966 logger.warning(
3967 'Slip scaling not performed. Maximum slip is 0.')
3968 else:
3969 disloc_est *= self.slip / disloc_tmax_max
3971 if not delta:
3972 return disloc_est, calc_times
3974 # if we have only one timestep there is no gradient
3975 if calc_times.size > 1:
3976 disloc_init = disloc_est[:, 0, :]
3977 disloc_est = num.diff(disloc_est, axis=1)
3978 disloc_est = num.concatenate((
3979 disloc_init[:, num.newaxis, :], disloc_est), axis=1)
3981 calc_times = calc_times
3983 return disloc_est, calc_times
3985 def get_slip_rate(self, *args, **kwargs):
3986 '''
3987 Get slip rate inverted from patches.
3989 The time interval, within which the slip rates are computed is
3990 determined by the sampling rate of the Green's function database or
3991 ``deltat``. Additional ``*args`` and ``**kwargs`` are passed to
3992 :py:meth:`get_delta_slip`.
3994 :returns:
3995 Slip rates (:math:`\\Delta u_{strike}/\\Delta t`,
3996 :math:`\\Delta u_{dip}/\\Delta t, \\Delta u_{tensile}/\\Delta t`)
3997 for each source patch and time; corner times, for which slip rate
3998 is computed. The order of sliprate array is:
4000 .. math::
4002 &[[\\\\
4003 &[\\Delta u_{strike, patch1, t1}/\\Delta t,
4004 \\Delta u_{dip, patch1, t1}/\\Delta t,
4005 \\Delta u_{tensile, patch1, t1}/\\Delta t],\\\\
4006 &[\\Delta u_{strike, patch1, t2}/\\Delta t,
4007 \\Delta u_{dip, patch1, t2}/\\Delta t,
4008 \\Delta u_{tensile, patch1, t2}/\\Delta t]], [\\\\
4009 &[\\Delta u_{strike, patch2, t1}/\\Delta t, ...],\\\\
4010 &[\\Delta u_{strike, patch2, t2}/\\Delta t, ...]]]\\\\
4012 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
4013 :py:class:`~numpy.ndarray`: ``(n_times, )``
4014 '''
4015 ddisloc_est, calc_times = self.get_delta_slip(
4016 *args, delta=True, **kwargs)
4018 dt = num.concatenate(
4019 [(num.diff(calc_times)[0], ), num.diff(calc_times)])
4020 slip_rate = num.linalg.norm(ddisloc_est, axis=2) / dt
4022 return slip_rate, calc_times
4024 def get_moment_rate_patches(self, *args, **kwargs):
4025 '''
4026 Get scalar seismic moment rate for each patch individually.
4028 Additional ``*args`` and ``**kwargs`` are passed to
4029 :py:meth:`get_slip_rate`.
4031 :returns:
4032 Seismic moment rate for each source patch and time; corner times,
4033 for which patch moment rate is computed based on slip rate. The
4034 order of the moment rate array is:
4036 .. math::
4038 &[\\\\
4039 &[(\\Delta M / \\Delta t)_{patch1, t1},
4040 (\\Delta M / \\Delta t)_{patch1, t2}, ...],\\\\
4041 &[(\\Delta M / \\Delta t)_{patch2, t1},
4042 (\\Delta M / \\Delta t)_{patch, t2}, ...],\\\\
4043 &[...]]\\\\
4045 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times)``,
4046 :py:class:`~numpy.ndarray`: ``(n_times, )``
4047 '''
4048 slip_rate, calc_times = self.get_slip_rate(*args, **kwargs)
4050 shear_mod = self.get_patch_attribute('shearmod')
4051 p_length = self.get_patch_attribute('length')
4052 p_width = self.get_patch_attribute('width')
4054 dA = p_length * p_width
4056 mom_rate = shear_mod[:, num.newaxis] * slip_rate * dA[:, num.newaxis]
4058 return mom_rate, calc_times
4060 def get_moment_rate(self, store, target=None, deltat=None):
4061 '''
4062 Get seismic source moment rate for the total source (STF).
4064 :param store:
4065 Green's function database (needs to cover whole region of of the
4066 source). Its ``deltat`` [s] is used as time increment for slip
4067 difference calculation. Either ``deltat`` or ``store`` should be
4068 given.
4069 :type store:
4070 :py:class:`~pyrocko.gf.store.Store`
4072 :param target:
4073 Target information, needed for interpolation method.
4074 :type target:
4075 :py:class:`~pyrocko.gf.targets.Target`
4077 :param deltat:
4078 Time increment for slip difference calculation [s]. If not given
4079 ``store.deltat`` is used.
4080 :type deltat:
4081 float
4083 :return:
4084 Seismic moment rate [Nm/s] for each time; corner times, for which
4085 moment rate is computed. The order of the moment rate array is:
4087 .. math::
4089 &[\\\\
4090 &(\\Delta M / \\Delta t)_{t1},\\\\
4091 &(\\Delta M / \\Delta t)_{t2},\\\\
4092 &...]\\\\
4094 :rtype:
4095 :py:class:`~numpy.ndarray`: ``(n_times, )``,
4096 :py:class:`~numpy.ndarray`: ``(n_times, )``
4097 '''
4098 if not deltat:
4099 deltat = store.config.deltat
4100 return self.discretize_basesource(
4101 store, target=target).get_moment_rate(deltat)
4103 def get_moment(self, *args, **kwargs):
4104 '''
4105 Get cumulative seismic moment.
4107 Additional ``*args`` and ``**kwargs`` are passed to
4108 :py:meth:`get_magnitude`.
4110 :returns:
4111 Cumulative seismic moment in [Nm].
4112 :rtype:
4113 float
4114 '''
4115 return float(pmt.magnitude_to_moment(self.get_magnitude(
4116 *args, **kwargs)))
4118 def rescale_slip(self, magnitude=None, moment=None, **kwargs):
4119 '''
4120 Rescale source slip based on given target magnitude or seismic moment.
4122 Rescale the maximum source slip to fit the source moment magnitude or
4123 seismic moment to the given target values. Either ``magnitude`` or
4124 ``moment`` need to be given. Additional ``**kwargs`` are passed to
4125 :py:meth:`get_moment`.
4127 :param magnitude:
4128 Target moment magnitude :math:`M_\\mathrm{w}` as in
4129 [Hanks and Kanamori, 1979]
4130 :type magnitude:
4131 float
4133 :param moment:
4134 Target seismic moment :math:`M_0` [Nm].
4135 :type moment:
4136 float
4137 '''
4138 if self.slip is None:
4139 self.slip = 1.
4140 logger.warning('No slip found for rescaling. '
4141 'An initial slip of 1 m is assumed.')
4143 if magnitude is None and moment is None:
4144 raise ValueError(
4145 'Either target magnitude or moment need to be given.')
4147 moment_init = self.get_moment(**kwargs)
4149 if magnitude is not None:
4150 moment = pmt.magnitude_to_moment(magnitude)
4152 self.slip *= moment / moment_init
4154 def get_centroid(self, store, *args, **kwargs):
4155 '''
4156 Centroid of the pseudo dynamic rupture model.
4158 The centroid location and time are derived from the locations and times
4159 of the individual patches weighted with their moment contribution.
4160 Additional ``**kwargs`` are passed to :py:meth:`pyrocko_moment_tensor`.
4162 :param store:
4163 Green's function database (needs to cover whole region of of the
4164 source). Its ``deltat`` [s] is used as time increment for slip
4165 difference calculation. Either ``deltat`` or ``store`` should be
4166 given.
4167 :type store:
4168 :py:class:`~pyrocko.gf.store.Store`
4170 :returns:
4171 The centroid location and associated moment tensor.
4172 :rtype:
4173 :py:class:`pyrocko.model.event.Event`
4174 '''
4175 _, _, _, _, time, _ = self.get_vr_time_interpolators(store)
4176 t_max = time.values.max()
4178 moment_rate, times = self.get_moment_rate_patches(deltat=t_max)
4180 moment = num.sum(moment_rate * times, axis=1)
4181 weights = moment / moment.sum()
4183 norths = self.get_patch_attribute('north_shift')
4184 easts = self.get_patch_attribute('east_shift')
4185 depths = self.get_patch_attribute('depth')
4187 centroid_n = num.sum(weights * norths)
4188 centroid_e = num.sum(weights * easts)
4189 centroid_d = num.sum(weights * depths)
4191 centroid_lat, centroid_lon = ne_to_latlon(
4192 self.lat, self.lon, centroid_n, centroid_e)
4194 moment_rate_, times = self.get_moment_rate(store)
4195 delta_times = num.concatenate((
4196 [times[1] - times[0]],
4197 num.diff(times)))
4198 moment_src = delta_times * moment_rate
4200 centroid_t = num.sum(
4201 moment_src / num.sum(moment_src) * times) + self.time
4203 mt = self.pyrocko_moment_tensor(store, *args, **kwargs)
4205 return model.Event(
4206 lat=centroid_lat,
4207 lon=centroid_lon,
4208 depth=centroid_d,
4209 time=centroid_t,
4210 moment_tensor=mt,
4211 magnitude=mt.magnitude,
4212 duration=t_max)
4214 def get_coulomb_failure_stress(
4215 self,
4216 receiver_points,
4217 friction,
4218 pressure,
4219 strike,
4220 dip,
4221 rake,
4222 time=None,
4223 *args,
4224 **kwargs):
4225 '''
4226 Calculate Coulomb failure stress change CFS.
4228 The function obtains the Coulomb failure stress change :math:`\\Delta
4229 \\sigma_C` at arbitrary receiver points with a commonly oriented
4230 receiver plane assuming:
4232 .. math::
4234 \\Delta \\sigma_C = \\sigma_S - \\mu (\\sigma_N - \\Delta p)
4236 with the shear stress :math:`\\sigma_S`, the coefficient of friction
4237 :math:`\\mu`, the normal stress :math:`\\sigma_N`, and the pore fluid
4238 pressure change :math:`\\Delta p`. Each receiver point is characterized
4239 by its geographical coordinates, and depth. The required receiver plane
4240 orientation is defined by ``strike``, ``dip``, and ``rake``. The
4241 Coulomb failure stress change is calculated for a given time after
4242 rupture origin time.
4244 :param receiver_points:
4245 Location of the receiver points in Northing, Easting, and depth in
4246 [m].
4247 :type receiver_points:
4248 :py:class:`~numpy.ndarray`: ``(n_receiver, 3)``
4250 :param friction:
4251 Coefficient of friction.
4252 :type friction:
4253 float
4255 :param pressure:
4256 Pore pressure change in [Pa].
4257 :type pressure:
4258 float
4260 :param strike:
4261 Strike of the receiver plane in [deg].
4262 :type strike:
4263 float
4265 :param dip:
4266 Dip of the receiver plane in [deg].
4267 :type dip:
4268 float
4270 :param rake:
4271 Rake of the receiver plane in [deg].
4272 :type rake:
4273 float
4275 :param time:
4276 Time after origin [s], for which the resulting :math:`\\Delta
4277 \\Sigma_c` is computed. If not given, :math:`\\Delta \\Sigma_c` is
4278 derived based on the final static slip.
4279 :type time:
4280 float
4282 :returns:
4283 The Coulomb failure stress change :math:`\\Delta \\Sigma_c` at each
4284 receiver point in [Pa].
4285 :rtype:
4286 :py:class:`~numpy.ndarray`: ``(n_receiver,)``
4287 '''
4288 # dislocation at given time
4289 source_slip = self.get_slip(time=time, scale_slip=True)
4291 # source planes
4292 source_patches = num.array([
4293 src.source_patch() for src in self.patches])
4295 # earth model
4296 lambda_mean = num.mean([src.lamb for src in self.patches])
4297 mu_mean = num.mean([src.shearmod for src in self.patches])
4299 # Dislocation and spatial derivatives from okada in NED
4300 results = okada_ext.okada(
4301 source_patches,
4302 source_slip,
4303 receiver_points,
4304 lambda_mean,
4305 mu_mean,
4306 rotate_sdn=False, # TODO Check
4307 stack_sources=0, # TODO Check
4308 *args, **kwargs)
4310 # resolve stress tensor (sum!)
4311 diag_ind = [0, 4, 8]
4312 kron = num.zeros(9)
4313 kron[diag_ind] = 1.
4314 kron = kron[num.newaxis, num.newaxis, :]
4316 eps = 0.5 * (
4317 results[:, :, 3:] +
4318 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
4320 dilatation \
4321 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
4323 stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps
4325 # superposed stress of all sources at receiver locations
4326 stress_sum = num.sum(stress, axis=0)
4328 # get shear and normal stress from stress tensor
4329 strike_rad = d2r * strike
4330 dip_rad = d2r * dip
4331 rake_rad = d2r * rake
4333 n_rec = receiver_points.shape[0]
4334 stress_normal = num.zeros(n_rec)
4335 tau = num.zeros(n_rec)
4337 # Get vectors in receiver fault normal (ns), strike (rst) and
4338 # dip (rdi) direction
4339 ns = num.zeros(3)
4340 rst = num.zeros(3)
4341 rdi = num.zeros(3)
4343 ns[0] = num.sin(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4344 ns[1] = num.sin(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4345 ns[2] = -num.cos(dip_rad)
4347 rst[0] = num.cos(strike_rad)
4348 rst[1] = num.sin(strike_rad)
4349 rst[2] = 0.0
4351 rdi[0] = num.cos(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4352 rdi[1] = num.cos(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4353 rdi[2] = num.sin(dip_rad)
4355 ts = rst * num.cos(rake_rad) - rdi * num.sin(rake_rad)
4357 stress_normal = num.sum(
4358 num.tile(ns, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4360 tau = num.sum(
4361 num.tile(ts, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4363 # calculate cfs using formula above and return
4364 return tau + friction * (stress_normal + pressure)
4367class DoubleDCSource(SourceWithMagnitude):
4368 '''
4369 Two double-couple point sources separated in space and time.
4370 Moment share between the sub-sources is controlled by the
4371 parameter mix.
4372 The position of the subsources is dependent on the moment
4373 distribution between the two sources. Depth, east and north
4374 shift are given for the centroid between the two double-couples.
4375 The subsources will positioned according to their moment shares
4376 around this centroid position.
4377 This is done according to their delta parameters, which are
4378 therefore in relation to that centroid.
4379 Note that depth of the subsources therefore can be
4380 depth+/-delta_depth. For shallow earthquakes therefore
4381 the depth has to be chosen deeper to avoid sampling
4382 above surface.
4383 '''
4385 strike1 = Float.T(
4386 default=0.0,
4387 help='strike direction in [deg], measured clockwise from north')
4389 dip1 = Float.T(
4390 default=90.0,
4391 help='dip angle in [deg], measured downward from horizontal')
4393 azimuth = Float.T(
4394 default=0.0,
4395 help='azimuth to second double-couple [deg], '
4396 'measured at first, clockwise from north')
4398 rake1 = Float.T(
4399 default=0.0,
4400 help='rake angle in [deg], '
4401 'measured counter-clockwise from right-horizontal '
4402 'in on-plane view')
4404 strike2 = Float.T(
4405 default=0.0,
4406 help='strike direction in [deg], measured clockwise from north')
4408 dip2 = Float.T(
4409 default=90.0,
4410 help='dip angle in [deg], measured downward from horizontal')
4412 rake2 = Float.T(
4413 default=0.0,
4414 help='rake angle in [deg], '
4415 'measured counter-clockwise from right-horizontal '
4416 'in on-plane view')
4418 delta_time = Float.T(
4419 default=0.0,
4420 help='separation of double-couples in time (t2-t1) [s]')
4422 delta_depth = Float.T(
4423 default=0.0,
4424 help='difference in depth (z2-z1) [m]')
4426 distance = Float.T(
4427 default=0.0,
4428 help='distance between the two double-couples [m]')
4430 mix = Float.T(
4431 default=0.5,
4432 help='how to distribute the moment to the two doublecouples '
4433 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1')
4435 stf1 = STF.T(
4436 optional=True,
4437 help='Source time function of subsource 1 '
4438 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4440 stf2 = STF.T(
4441 optional=True,
4442 help='Source time function of subsource 2 '
4443 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4445 discretized_source_class = meta.DiscretizedMTSource
4447 def base_key(self):
4448 return (
4449 self.time, self.depth, self.lat, self.north_shift,
4450 self.lon, self.east_shift, type(self).__name__) + \
4451 self.effective_stf1_pre().base_key() + \
4452 self.effective_stf2_pre().base_key() + (
4453 self.strike1, self.dip1, self.rake1,
4454 self.strike2, self.dip2, self.rake2,
4455 self.delta_time, self.delta_depth,
4456 self.azimuth, self.distance, self.mix)
4458 def get_factor(self):
4459 return self.moment
4461 def effective_stf1_pre(self):
4462 return self.stf1 or self.stf or g_unit_pulse
4464 def effective_stf2_pre(self):
4465 return self.stf2 or self.stf or g_unit_pulse
4467 def effective_stf_post(self):
4468 return g_unit_pulse
4470 def split(self):
4471 a1 = 1.0 - self.mix
4472 a2 = self.mix
4473 delta_north = math.cos(self.azimuth * d2r) * self.distance
4474 delta_east = math.sin(self.azimuth * d2r) * self.distance
4476 dc1 = DCSource(
4477 lat=self.lat,
4478 lon=self.lon,
4479 time=self.time - self.delta_time * a2,
4480 north_shift=self.north_shift - delta_north * a2,
4481 east_shift=self.east_shift - delta_east * a2,
4482 depth=self.depth - self.delta_depth * a2,
4483 moment=self.moment * a1,
4484 strike=self.strike1,
4485 dip=self.dip1,
4486 rake=self.rake1,
4487 stf=self.stf1 or self.stf)
4489 dc2 = DCSource(
4490 lat=self.lat,
4491 lon=self.lon,
4492 time=self.time + self.delta_time * a1,
4493 north_shift=self.north_shift + delta_north * a1,
4494 east_shift=self.east_shift + delta_east * a1,
4495 depth=self.depth + self.delta_depth * a1,
4496 moment=self.moment * a2,
4497 strike=self.strike2,
4498 dip=self.dip2,
4499 rake=self.rake2,
4500 stf=self.stf2 or self.stf)
4502 return [dc1, dc2]
4504 def discretize_basesource(self, store, target=None):
4505 a1 = 1.0 - self.mix
4506 a2 = self.mix
4507 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4508 rake=self.rake1, scalar_moment=a1)
4509 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4510 rake=self.rake2, scalar_moment=a2)
4512 delta_north = math.cos(self.azimuth * d2r) * self.distance
4513 delta_east = math.sin(self.azimuth * d2r) * self.distance
4515 times1, amplitudes1 = self.effective_stf1_pre().discretize_t(
4516 store.config.deltat, self.time - self.delta_time * a2)
4518 times2, amplitudes2 = self.effective_stf2_pre().discretize_t(
4519 store.config.deltat, self.time + self.delta_time * a1)
4521 nt1 = times1.size
4522 nt2 = times2.size
4524 ds = meta.DiscretizedMTSource(
4525 lat=self.lat,
4526 lon=self.lon,
4527 times=num.concatenate((times1, times2)),
4528 north_shifts=num.concatenate((
4529 num.repeat(self.north_shift - delta_north * a2, nt1),
4530 num.repeat(self.north_shift + delta_north * a1, nt2))),
4531 east_shifts=num.concatenate((
4532 num.repeat(self.east_shift - delta_east * a2, nt1),
4533 num.repeat(self.east_shift + delta_east * a1, nt2))),
4534 depths=num.concatenate((
4535 num.repeat(self.depth - self.delta_depth * a2, nt1),
4536 num.repeat(self.depth + self.delta_depth * a1, nt2))),
4537 m6s=num.vstack((
4538 mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis],
4539 mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis])))
4541 return ds
4543 def pyrocko_moment_tensor(self, store=None, target=None):
4544 a1 = 1.0 - self.mix
4545 a2 = self.mix
4546 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4547 rake=self.rake1,
4548 scalar_moment=a1 * self.moment)
4549 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4550 rake=self.rake2,
4551 scalar_moment=a2 * self.moment)
4552 return pmt.MomentTensor(m=mot1.m() + mot2.m())
4554 def pyrocko_event(self, store=None, target=None, **kwargs):
4555 return SourceWithMagnitude.pyrocko_event(
4556 self, store, target,
4557 moment_tensor=self.pyrocko_moment_tensor(store, target),
4558 **kwargs)
4560 @classmethod
4561 def from_pyrocko_event(cls, ev, **kwargs):
4562 d = {}
4563 mt = ev.moment_tensor
4564 if mt:
4565 (strike, dip, rake), _ = mt.both_strike_dip_rake()
4566 d.update(
4567 strike1=float(strike),
4568 dip1=float(dip),
4569 rake1=float(rake),
4570 strike2=float(strike),
4571 dip2=float(dip),
4572 rake2=float(rake),
4573 mix=0.0,
4574 magnitude=float(mt.moment_magnitude()))
4576 d.update(kwargs)
4577 source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d)
4578 source.stf1 = source.stf
4579 source.stf2 = HalfSinusoidSTF(effective_duration=0.)
4580 source.stf = None
4581 return source
4584class RingfaultSource(SourceWithMagnitude):
4585 '''
4586 A ring fault with vertical doublecouples.
4587 '''
4589 diameter = Float.T(
4590 default=1.0,
4591 help='diameter of the ring in [m]')
4593 sign = Float.T(
4594 default=1.0,
4595 help='inside of the ring moves up (+1) or down (-1)')
4597 strike = Float.T(
4598 default=0.0,
4599 help='strike direction of the ring plane, clockwise from north,'
4600 ' in [deg]')
4602 dip = Float.T(
4603 default=0.0,
4604 help='dip angle of the ring plane from horizontal in [deg]')
4606 npointsources = Int.T(
4607 default=360,
4608 help='number of point sources to use')
4610 discretized_source_class = meta.DiscretizedMTSource
4612 def base_key(self):
4613 return Source.base_key(self) + (
4614 self.strike, self.dip, self.diameter, self.npointsources)
4616 def get_factor(self):
4617 return self.sign * self.moment
4619 def discretize_basesource(self, store=None, target=None):
4620 n = self.npointsources
4621 phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False)
4623 points = num.zeros((n, 3))
4624 points[:, 0] = num.cos(phi) * 0.5 * self.diameter
4625 points[:, 1] = num.sin(phi) * 0.5 * self.diameter
4627 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
4628 points = num.dot(rotmat.T, points.T).T # !!! ?
4630 points[:, 0] += self.north_shift
4631 points[:, 1] += self.east_shift
4632 points[:, 2] += self.depth
4634 m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90.,
4635 scalar_moment=1.0 / n).m())
4637 rotmats = num.transpose(
4638 [[num.cos(phi), num.sin(phi), num.zeros(n)],
4639 [-num.sin(phi), num.cos(phi), num.zeros(n)],
4640 [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1))
4642 ms = num.zeros((n, 3, 3))
4643 for i in range(n):
4644 mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i]))
4645 ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat))
4647 m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2],
4648 ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T
4650 times, amplitudes = self.effective_stf_pre().discretize_t(
4651 store.config.deltat, self.time)
4653 nt = times.size
4655 return meta.DiscretizedMTSource(
4656 times=num.tile(times, n),
4657 lat=self.lat,
4658 lon=self.lon,
4659 north_shifts=num.repeat(points[:, 0], nt),
4660 east_shifts=num.repeat(points[:, 1], nt),
4661 depths=num.repeat(points[:, 2], nt),
4662 m6s=num.repeat(m6s, nt, axis=0) * num.tile(
4663 amplitudes, n)[:, num.newaxis])
4666class CombiSource(Source):
4667 '''
4668 Composite source model.
4669 '''
4671 discretized_source_class = meta.DiscretizedMTSource
4673 subsources = List.T(Source.T())
4675 def __init__(self, subsources=[], **kwargs):
4676 if not subsources:
4677 raise BadRequest(
4678 'Need at least one sub-source to create a CombiSource object.')
4680 lats = num.array(
4681 [subsource.lat for subsource in subsources], dtype=float)
4682 lons = num.array(
4683 [subsource.lon for subsource in subsources], dtype=float)
4685 lat, lon = lats[0], lons[0]
4686 if not num.all(lats == lat) and num.all(lons == lon):
4687 subsources = [s.clone() for s in subsources]
4688 for subsource in subsources[1:]:
4689 subsource.set_origin(lat, lon)
4691 depth = float(num.mean([p.depth for p in subsources]))
4692 time = float(num.mean([p.time for p in subsources]))
4693 north_shift = float(num.mean([p.north_shift for p in subsources]))
4694 east_shift = float(num.mean([p.east_shift for p in subsources]))
4695 kwargs.update(
4696 time=time,
4697 lat=float(lat),
4698 lon=float(lon),
4699 north_shift=north_shift,
4700 east_shift=east_shift,
4701 depth=depth)
4703 Source.__init__(self, subsources=subsources, **kwargs)
4705 def get_factor(self):
4706 return 1.0
4708 def discretize_basesource(self, store, target=None):
4709 dsources = []
4710 for sf in self.subsources:
4711 ds = sf.discretize_basesource(store, target)
4712 ds.m6s *= sf.get_factor()
4713 dsources.append(ds)
4715 return meta.DiscretizedMTSource.combine(dsources)
4718class CombiSFSource(Source):
4719 '''
4720 Composite source model.
4721 '''
4723 discretized_source_class = meta.DiscretizedSFSource
4725 subsources = List.T(Source.T())
4727 def __init__(self, subsources=[], **kwargs):
4728 if not subsources:
4729 raise BadRequest(
4730 'Need at least one sub-source to create a CombiSFSource '
4731 'object.')
4733 lats = num.array(
4734 [subsource.lat for subsource in subsources], dtype=float)
4735 lons = num.array(
4736 [subsource.lon for subsource in subsources], dtype=float)
4738 lat, lon = lats[0], lons[0]
4739 if not num.all(lats == lat) and num.all(lons == lon):
4740 subsources = [s.clone() for s in subsources]
4741 for subsource in subsources[1:]:
4742 subsource.set_origin(lat, lon)
4744 depth = float(num.mean([p.depth for p in subsources]))
4745 time = float(num.mean([p.time for p in subsources]))
4746 north_shift = float(num.mean([p.north_shift for p in subsources]))
4747 east_shift = float(num.mean([p.east_shift for p in subsources]))
4748 kwargs.update(
4749 time=time,
4750 lat=float(lat),
4751 lon=float(lon),
4752 north_shift=north_shift,
4753 east_shift=east_shift,
4754 depth=depth)
4756 Source.__init__(self, subsources=subsources, **kwargs)
4758 def get_factor(self):
4759 return 1.0
4761 def discretize_basesource(self, store, target=None):
4762 dsources = []
4763 for sf in self.subsources:
4764 ds = sf.discretize_basesource(store, target)
4765 ds.forces *= sf.get_factor()
4766 dsources.append(ds)
4768 return meta.DiscretizedSFSource.combine(dsources)
4771class SFSource(Source):
4772 '''
4773 A single force point source.
4775 Supported GF schemes: `'elastic5'`.
4776 '''
4778 discretized_source_class = meta.DiscretizedSFSource
4780 fn = Float.T(
4781 default=0.,
4782 help='northward component of single force [N]')
4784 fe = Float.T(
4785 default=0.,
4786 help='eastward component of single force [N]')
4788 fd = Float.T(
4789 default=0.,
4790 help='downward component of single force [N]')
4792 def __init__(self, **kwargs):
4793 Source.__init__(self, **kwargs)
4795 def base_key(self):
4796 return Source.base_key(self) + (self.fn, self.fe, self.fd)
4798 def get_factor(self):
4799 return 1.0
4801 @property
4802 def force(self):
4803 return math.sqrt(self.fn**2 + self.fe**2 + self.fd**2)
4805 def discretize_basesource(self, store, target=None):
4806 times, amplitudes = self.effective_stf_pre().discretize_t(
4807 store.config.deltat, self.time)
4808 forces = amplitudes[:, num.newaxis] * num.array(
4809 [[self.fn, self.fe, self.fd]], dtype=float)
4811 return meta.DiscretizedSFSource(forces=forces,
4812 **self._dparams_base_repeated(times))
4814 def pyrocko_event(self, store=None, target=None, **kwargs):
4815 return Source.pyrocko_event(
4816 self, store, target,
4817 **kwargs)
4819 @classmethod
4820 def from_pyrocko_event(cls, ev, **kwargs):
4821 d = {}
4822 d.update(kwargs)
4823 return super(SFSource, cls).from_pyrocko_event(ev, **d)
4826class SimpleLandslideSource(Source):
4827 '''
4828 A single force landslide source respecting conservation of momentum.
4830 The landslide is modelled point-like in space but with individual source
4831 time functions for each force component. The source time functions
4832 :py:class:`SimpleLandslideSTF` impose the constraint that everything is at
4833 rest before and after the event but are allowed to have different
4834 acceleration and deceleration durations. It should thus be suitable as a
4835 first order approximation of a rock fall or landslide.
4836 For realistic landslides, the horizontal accelerations and decelerations
4837 can be delayed with respect to the vertical ones but cannot precede them.
4839 Supported GF schemes: `'elastic5'`.
4840 '''
4842 discretized_source_class = meta.DiscretizedSFSource
4844 stf_mode = STFMode.T(
4845 default='pre',
4846 help='SimpleLandslideSource only works with `stf_mode == "pre"`.')
4848 impulse_n = Float.T(
4849 default=0.,
4850 help='northward component of impulse [Ns]')
4852 impulse_e = Float.T(
4853 default=0.,
4854 help='eastward component of impulse [Ns]')
4856 impulse_d = Float.T(
4857 default=0.,
4858 help='downward component of impulse [Ns]')
4860 azimuth = Float.T(
4861 default=0.,
4862 help='azimuth direction of the mass movement [deg]')
4864 stf_v = SimpleLandslideSTF.T(
4865 default=SimpleLandslideSTF.D(),
4866 help='source time function for vertical force component')
4868 stf_h = SimpleLandslideSTF.T(
4869 default=SimpleLandslideSTF.D(),
4870 help='source time function for horizontal force component')
4872 anchor_stf = StringChoice.T(
4873 choices=['onset', 'centroid'],
4874 default='onset',
4875 help='``"onset"``: STFs start at origin time ``"centroid"``: STFs all '
4876 'start at the same time but so that the centroid is at the given '
4877 'origin time.')
4879 def __init__(self, **kwargs):
4880 Source.__init__(self, **kwargs)
4882 def base_key(self):
4883 return Source.base_key(self) + (
4884 self.impulse_n, self.impulse_e, self.impulse_d) \
4885 + self.stf_v.base_key() + self.stf_h.base_key()
4887 def get_factor(self):
4888 return 1.0
4890 def discretize_basesource(self, store, target=None):
4891 if self.stf_mode != 'pre':
4892 raise Exception(
4893 'SimpleLandslideSource: '
4894 'Only works with `stf_mode == "pre"`.')
4896 if self.stf is not None:
4897 raise Exception(
4898 'SimpleLandslideSource: '
4899 'Setting `stf` is not supported: use `stf_v` and `stf_h`.')
4901 if self.anchor_stf == 'centroid':
4902 duration_acc = num.array([
4903 self.stf_h.duration_acceleration,
4904 self.stf_h.duration_acceleration,
4905 self.stf_v.duration_acceleration], dtype=float)
4907 impulse = num.array([
4908 self.impulse_n,
4909 self.impulse_e,
4910 self.impulse_d], dtype=float)
4912 tshift_centroid = \
4913 - num.sum(duration_acc * impulse**2) \
4914 / num.sum(impulse**2)
4916 elif self.anchor_stf == 'onset':
4917 tshift_centroid = 0.0
4919 times, amplitudes = self.stf_v.discretize_t(
4920 store.config.deltat,
4921 self.time + tshift_centroid)
4923 forces_v = num.zeros((times.size, 3))
4924 forces_v[:, 2] = amplitudes * self.impulse_d
4926 dsource_v = meta.DiscretizedSFSource(
4927 forces=forces_v,
4928 **self._dparams_base_repeated(times))
4930 times, amplitudes = self.stf_h.discretize_t(
4931 store.config.deltat,
4932 self.time + tshift_centroid)
4934 forces_h = num.zeros((times.size, 3))
4935 forces_h[:, 0] = \
4936 amplitudes * self.impulse_n
4937 forces_h[:, 1] = \
4938 amplitudes * self.impulse_e
4940 dsource_h = meta.DiscretizedSFSource(
4941 forces=forces_h,
4942 **self._dparams_base_repeated(times))
4944 return meta.DiscretizedSFSource.combine([dsource_v, dsource_h])
4946 def pyrocko_event(self, store=None, target=None, **kwargs):
4947 return Source.pyrocko_event(
4948 self, store, target,
4949 **kwargs)
4951 @classmethod
4952 def from_pyrocko_event(cls, ev, **kwargs):
4953 d = {}
4954 d.update(kwargs)
4955 return super(SimpleLandslideSource, cls).from_pyrocko_event(ev, **d)
4958class PorePressurePointSource(Source):
4959 '''
4960 Excess pore pressure point source.
4962 For poro-elastic initial value problem where an excess pore pressure is
4963 brought into a small source volume.
4964 '''
4966 discretized_source_class = meta.DiscretizedPorePressureSource
4968 pp = Float.T(
4969 default=1.0,
4970 help='initial excess pore pressure in [Pa]')
4972 def base_key(self):
4973 return Source.base_key(self)
4975 def get_factor(self):
4976 return self.pp
4978 def discretize_basesource(self, store, target=None):
4979 return meta.DiscretizedPorePressureSource(pp=arr(1.0),
4980 **self._dparams_base())
4983class PorePressureLineSource(Source):
4984 '''
4985 Excess pore pressure line source.
4987 The line source is centered at (north_shift, east_shift, depth).
4988 '''
4990 discretized_source_class = meta.DiscretizedPorePressureSource
4992 pp = Float.T(
4993 default=1.0,
4994 help='initial excess pore pressure in [Pa]')
4996 length = Float.T(
4997 default=0.0,
4998 help='length of the line source [m]')
5000 azimuth = Float.T(
5001 default=0.0,
5002 help='azimuth direction, clockwise from north [deg]')
5004 dip = Float.T(
5005 default=90.,
5006 help='dip direction, downward from horizontal [deg]')
5008 def base_key(self):
5009 return Source.base_key(self) + (self.azimuth, self.dip, self.length)
5011 def get_factor(self):
5012 return self.pp
5014 def discretize_basesource(self, store, target=None):
5016 n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1
5018 a = num.linspace(-0.5 * self.length, 0.5 * self.length, n)
5020 sa = math.sin(self.azimuth * d2r)
5021 ca = math.cos(self.azimuth * d2r)
5022 sd = math.sin(self.dip * d2r)
5023 cd = math.cos(self.dip * d2r)
5025 points = num.zeros((n, 3))
5026 points[:, 0] = self.north_shift + a * ca * cd
5027 points[:, 1] = self.east_shift + a * sa * cd
5028 points[:, 2] = self.depth + a * sd
5030 return meta.DiscretizedPorePressureSource(
5031 times=num.full(n, self.time),
5032 lat=self.lat,
5033 lon=self.lon,
5034 north_shifts=points[:, 0],
5035 east_shifts=points[:, 1],
5036 depths=points[:, 2],
5037 pp=num.ones(n) / n)
5040class Request(Object):
5041 '''
5042 Synthetic seismogram computation request.
5044 ::
5046 Request(**kwargs)
5047 Request(sources, targets, **kwargs)
5048 '''
5050 sources = List.T(
5051 Source.T(),
5052 help='list of sources for which to produce synthetics.')
5054 targets = List.T(
5055 Target.T(),
5056 help='list of targets for which to produce synthetics.')
5058 @classmethod
5059 def args2kwargs(cls, args):
5060 if len(args) not in (0, 2, 3):
5061 raise BadRequest('Invalid arguments.')
5063 if len(args) == 2:
5064 return dict(sources=args[0], targets=args[1])
5065 else:
5066 return {}
5068 def __init__(self, *args, **kwargs):
5069 kwargs.update(self.args2kwargs(args))
5070 sources = kwargs.pop('sources', [])
5071 targets = kwargs.pop('targets', [])
5073 if isinstance(sources, Source):
5074 sources = [sources]
5076 if isinstance(targets, Target) or isinstance(targets, StaticTarget):
5077 targets = [targets]
5079 Object.__init__(self, sources=sources, targets=targets, **kwargs)
5081 @property
5082 def targets_dynamic(self):
5083 return [t for t in self.targets if isinstance(t, Target)]
5085 @property
5086 def targets_static(self):
5087 return [t for t in self.targets if isinstance(t, StaticTarget)]
5089 @property
5090 def has_dynamic(self):
5091 return True if len(self.targets_dynamic) > 0 else False
5093 @property
5094 def has_statics(self):
5095 return True if len(self.targets_static) > 0 else False
5097 def subsources_map(self):
5098 m = defaultdict(list)
5099 for source in self.sources:
5100 m[source.base_key()].append(source)
5102 return m
5104 def subtargets_map(self):
5105 m = defaultdict(list)
5106 for target in self.targets:
5107 m[target.base_key()].append(target)
5109 return m
5111 def subrequest_map(self):
5112 ms = self.subsources_map()
5113 mt = self.subtargets_map()
5114 m = {}
5115 for (ks, ls) in ms.items():
5116 for (kt, lt) in mt.items():
5117 m[ks, kt] = (ls, lt)
5119 return m
5122class ProcessingStats(Object):
5123 t_perc_get_store_and_receiver = Float.T(default=0.)
5124 t_perc_discretize_source = Float.T(default=0.)
5125 t_perc_make_base_seismogram = Float.T(default=0.)
5126 t_perc_make_same_span = Float.T(default=0.)
5127 t_perc_post_process = Float.T(default=0.)
5128 t_perc_optimize = Float.T(default=0.)
5129 t_perc_stack = Float.T(default=0.)
5130 t_perc_static_get_store = Float.T(default=0.)
5131 t_perc_static_discretize_basesource = Float.T(default=0.)
5132 t_perc_static_sum_statics = Float.T(default=0.)
5133 t_perc_static_post_process = Float.T(default=0.)
5134 t_wallclock = Float.T(default=0.)
5135 t_cpu = Float.T(default=0.)
5136 n_read_blocks = Int.T(default=0)
5137 n_results = Int.T(default=0)
5138 n_subrequests = Int.T(default=0)
5139 n_stores = Int.T(default=0)
5140 n_records_stacked = Int.T(default=0)
5143class Response(Object):
5144 '''
5145 Resonse object to a synthetic seismogram computation request.
5146 '''
5148 request = Request.T()
5149 results_list = List.T(List.T(meta.SeismosizerResult.T()))
5150 stats = ProcessingStats.T()
5152 def pyrocko_traces(self):
5153 '''
5154 Return a list of requested
5155 :class:`~pyrocko.trace.Trace` instances.
5156 '''
5158 traces = []
5159 for results in self.results_list:
5160 for result in results:
5161 if not isinstance(result, meta.Result):
5162 continue
5163 traces.append(result.trace.pyrocko_trace())
5165 return traces
5167 def kite_scenes(self):
5168 '''
5169 Return a list of requested
5170 :class:`kite.Scene` instances.
5171 '''
5172 kite_scenes = []
5173 for results in self.results_list:
5174 for result in results:
5175 if isinstance(result, meta.KiteSceneResult):
5176 sc = result.get_scene()
5177 kite_scenes.append(sc)
5179 return kite_scenes
5181 def static_results(self):
5182 '''
5183 Return a list of requested
5184 :class:`~pyrocko.gf.meta.StaticResult` instances.
5185 '''
5186 statics = []
5187 for results in self.results_list:
5188 for result in results:
5189 if not isinstance(result, meta.StaticResult):
5190 continue
5191 statics.append(result)
5193 return statics
5195 def iter_results(self, get='pyrocko_traces'):
5196 '''
5197 Generator function to iterate over results of request.
5199 Yields associated :py:class:`Source`,
5200 :class:`~pyrocko.gf.targets.Target`,
5201 :class:`~pyrocko.trace.Trace` instances in each iteration.
5202 '''
5204 for isource, source in enumerate(self.request.sources):
5205 for itarget, target in enumerate(self.request.targets):
5206 result = self.results_list[isource][itarget]
5207 if get == 'pyrocko_traces':
5208 yield source, target, result.trace.pyrocko_trace()
5209 elif get == 'results':
5210 yield source, target, result
5212 def snuffle(self, **kwargs):
5213 '''
5214 Open *snuffler* with requested traces.
5215 '''
5217 trace.snuffle(self.pyrocko_traces(), **kwargs)
5220class Engine(Object):
5221 '''
5222 Base class for synthetic seismogram calculators.
5223 '''
5225 def get_store_ids(self):
5226 '''
5227 Get list of available GF store IDs
5228 '''
5230 return []
5233class Rule(object):
5234 pass
5237class VectorRule(Rule):
5239 def __init__(self, quantity, differentiate=0, integrate=0):
5240 self.components = [quantity + '.' + c for c in 'ned']
5241 self.differentiate = differentiate
5242 self.integrate = integrate
5244 def required_components(self, target):
5245 n, e, d = self.components
5246 sa, ca, sd, cd = target.get_sin_cos_factors()
5248 comps = []
5249 if nonzero(ca * cd):
5250 comps.append(n)
5252 if nonzero(sa * cd):
5253 comps.append(e)
5255 if nonzero(sd):
5256 comps.append(d)
5258 return tuple(comps)
5260 def apply_(self, target, base_seismogram):
5261 n, e, d = self.components
5262 sa, ca, sd, cd = target.get_sin_cos_factors()
5264 if nonzero(ca * cd):
5265 data = base_seismogram[n].data * (ca * cd)
5266 deltat = base_seismogram[n].deltat
5267 else:
5268 data = 0.0
5270 if nonzero(sa * cd):
5271 data = data + base_seismogram[e].data * (sa * cd)
5272 deltat = base_seismogram[e].deltat
5274 if nonzero(sd):
5275 data = data + base_seismogram[d].data * sd
5276 deltat = base_seismogram[d].deltat
5278 if self.differentiate:
5279 data = util.diff_fd(self.differentiate, 4, deltat, data)
5281 if self.integrate:
5282 raise NotImplementedError('Integration is not implemented yet.')
5284 return data
5287class HorizontalVectorRule(Rule):
5289 def __init__(self, quantity, differentiate=0, integrate=0):
5290 self.components = [quantity + '.' + c for c in 'ne']
5291 self.differentiate = differentiate
5292 self.integrate = integrate
5294 def required_components(self, target):
5295 n, e = self.components
5296 sa, ca, _, _ = target.get_sin_cos_factors()
5298 comps = []
5299 if nonzero(ca):
5300 comps.append(n)
5302 if nonzero(sa):
5303 comps.append(e)
5305 return tuple(comps)
5307 def apply_(self, target, base_seismogram):
5308 n, e = self.components
5309 sa, ca, _, _ = target.get_sin_cos_factors()
5311 if nonzero(ca):
5312 data = base_seismogram[n].data * ca
5313 else:
5314 data = 0.0
5316 if nonzero(sa):
5317 data = data + base_seismogram[e].data * sa
5319 if self.differentiate:
5320 deltat = base_seismogram[e].deltat
5321 data = util.diff_fd(self.differentiate, 4, deltat, data)
5323 if self.integrate:
5324 raise NotImplementedError('Integration is not implemented yet.')
5326 return data
5329class ScalarRule(Rule):
5331 def __init__(self, quantity, differentiate=0):
5332 self.c = quantity
5334 def required_components(self, target):
5335 return (self.c, )
5337 def apply_(self, target, base_seismogram):
5338 data = base_seismogram[self.c].data.copy()
5339 deltat = base_seismogram[self.c].deltat
5340 if self.differentiate:
5341 data = util.diff_fd(self.differentiate, 4, deltat, data)
5343 return data
5346class StaticDisplacement(Rule):
5348 def required_components(self, target):
5349 return tuple(['displacement.%s' % c for c in list('ned')])
5351 def apply_(self, target, base_statics):
5352 if isinstance(target, SatelliteTarget):
5353 los_fac = target.get_los_factors()
5354 base_statics['displacement.los'] =\
5355 (los_fac[:, 0] * -base_statics['displacement.d'] +
5356 los_fac[:, 1] * base_statics['displacement.e'] +
5357 los_fac[:, 2] * base_statics['displacement.n'])
5358 return base_statics
5361channel_rules = {
5362 'displacement': [VectorRule('displacement')],
5363 'rotation_displacement': [VectorRule('rotation_displacement')],
5364 'velocity': [
5365 VectorRule('velocity'),
5366 VectorRule('displacement', differentiate=1)],
5367 'acceleration': [
5368 VectorRule('acceleration'),
5369 VectorRule('velocity', differentiate=1),
5370 VectorRule('displacement', differentiate=2)],
5371 'pore_pressure': [ScalarRule('pore_pressure')],
5372 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')],
5373 'darcy_velocity': [VectorRule('darcy_velocity')],
5374}
5376static_rules = {
5377 'displacement': [StaticDisplacement()]
5378}
5381class OutOfBoundsContext(Object):
5382 source = Source.T()
5383 target = Target.T()
5384 distance = Float.T()
5385 components = List.T(String.T())
5388def process_dynamic_timeseries(work, psources, ptargets, engine, nthreads=0):
5389 dsource_cache = {}
5390 tcounters = list(range(6))
5392 store_ids = set()
5393 sources = set()
5394 targets = set()
5396 for itarget, target in enumerate(ptargets):
5397 target._id = itarget
5399 for w in work:
5400 _, _, isources, itargets = w
5402 sources.update([psources[isource] for isource in isources])
5403 targets.update([ptargets[itarget] for itarget in itargets])
5405 store_ids = set([t.store_id for t in targets])
5407 for isource, source in enumerate(psources):
5409 components = set()
5410 for itarget, target in enumerate(targets):
5411 rule = engine.get_rule(source, target)
5412 components.update(rule.required_components(target))
5414 for store_id in store_ids:
5415 store_targets = [t for t in targets if t.store_id == store_id]
5417 sample_rates = set([t.sample_rate for t in store_targets])
5418 interpolations = set([t.interpolation for t in store_targets])
5420 base_seismograms = []
5421 store_targets_out = []
5423 for samp_rate in sample_rates:
5424 for interp in interpolations:
5425 engine_targets = [
5426 t for t in store_targets if t.sample_rate == samp_rate
5427 and t.interpolation == interp]
5429 if not engine_targets:
5430 continue
5432 store_targets_out += engine_targets
5434 base_seismograms += engine.base_seismograms(
5435 source,
5436 engine_targets,
5437 components,
5438 dsource_cache,
5439 nthreads)
5441 for iseis, seismogram in enumerate(base_seismograms):
5442 for tr in seismogram.values():
5443 if tr.err != store.SeismosizerErrorEnum.SUCCESS:
5444 e = SeismosizerError(
5445 'Seismosizer failed with return code %i\n%s' % (
5446 tr.err, str(
5447 OutOfBoundsContext(
5448 source=source,
5449 target=store_targets[iseis],
5450 distance=source.distance_to(
5451 store_targets[iseis]),
5452 components=components))))
5453 raise e
5455 for seismogram, target in zip(base_seismograms, store_targets_out):
5457 try:
5458 result = engine._post_process_dynamic(
5459 seismogram, source, target)
5460 except SeismosizerError as e:
5461 result = e
5463 yield (isource, target._id, result), tcounters
5466def process_dynamic(work, psources, ptargets, engine, nthreads=0):
5467 dsource_cache = {}
5469 for w in work:
5470 _, _, isources, itargets = w
5472 sources = [psources[isource] for isource in isources]
5473 targets = [ptargets[itarget] for itarget in itargets]
5475 components = set()
5476 for target in targets:
5477 rule = engine.get_rule(sources[0], target)
5478 components.update(rule.required_components(target))
5480 for isource, source in zip(isources, sources):
5481 for itarget, target in zip(itargets, targets):
5483 try:
5484 base_seismogram, tcounters = engine.base_seismogram(
5485 source, target, components, dsource_cache, nthreads)
5486 except meta.OutOfBounds as e:
5487 e.context = OutOfBoundsContext(
5488 source=sources[0],
5489 target=targets[0],
5490 distance=sources[0].distance_to(targets[0]),
5491 components=components)
5492 raise
5494 n_records_stacked = 0
5495 t_optimize = 0.0
5496 t_stack = 0.0
5498 for _, tr in base_seismogram.items():
5499 n_records_stacked += tr.n_records_stacked
5500 t_optimize += tr.t_optimize
5501 t_stack += tr.t_stack
5503 try:
5504 result = engine._post_process_dynamic(
5505 base_seismogram, source, target)
5506 result.n_records_stacked = n_records_stacked
5507 result.n_shared_stacking = len(sources) *\
5508 len(targets)
5509 result.t_optimize = t_optimize
5510 result.t_stack = t_stack
5511 except SeismosizerError as e:
5512 result = e
5514 tcounters.append(xtime())
5515 yield (isource, itarget, result), tcounters
5518def process_static(work, psources, ptargets, engine, nthreads=0):
5519 for w in work:
5520 _, _, isources, itargets = w
5522 sources = [psources[isource] for isource in isources]
5523 targets = [ptargets[itarget] for itarget in itargets]
5525 for isource, source in zip(isources, sources):
5526 for itarget, target in zip(itargets, targets):
5527 components = engine.get_rule(source, target)\
5528 .required_components(target)
5530 try:
5531 base_statics, tcounters = engine.base_statics(
5532 source, target, components, nthreads)
5533 except meta.OutOfBounds as e:
5534 e.context = OutOfBoundsContext(
5535 source=sources[0],
5536 target=targets[0],
5537 distance=float('nan'),
5538 components=components)
5539 raise
5540 result = engine._post_process_statics(
5541 base_statics, source, target)
5542 tcounters.append(xtime())
5544 yield (isource, itarget, result), tcounters
5547class LocalEngine(Engine):
5548 '''
5549 Offline synthetic seismogram calculator.
5551 :param use_env: if ``True``, fill :py:attr:`store_superdirs` and
5552 :py:attr:`store_dirs` with paths set in environment variables
5553 GF_STORE_SUPERDIRS and GF_STORE_DIRS.
5554 :param use_config: if ``True``, fill :py:attr:`store_superdirs` and
5555 :py:attr:`store_dirs` with paths set in the user's config file.
5557 The config file can be found at :file:`~/.pyrocko/config.pf`
5559 .. code-block :: python
5561 gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/']
5562 gf_store_superdirs: ['/home/pyrocko/gf_stores/']
5563 '''
5565 store_superdirs = List.T(
5566 String.T(),
5567 help="directories which are searched for Green's function stores")
5569 store_dirs = List.T(
5570 String.T(),
5571 help="additional individual Green's function store directories")
5573 default_store_id = String.T(
5574 optional=True,
5575 help='default store ID to be used when a request does not provide '
5576 'one')
5578 nthreads = Int.T(
5579 default=1,
5580 help='default number of threads to utilize')
5582 def __init__(self, **kwargs):
5583 use_env = kwargs.pop('use_env', False)
5584 use_config = kwargs.pop('use_config', False)
5585 Engine.__init__(self, **kwargs)
5586 if use_env:
5587 env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '')
5588 env_store_dirs = os.environ.get('GF_STORE_DIRS', '')
5589 if env_store_superdirs:
5590 self.store_superdirs.extend(env_store_superdirs.split(':'))
5592 if env_store_dirs:
5593 self.store_dirs.extend(env_store_dirs.split(':'))
5595 if use_config:
5596 c = config.config()
5597 self.store_superdirs.extend(c.gf_store_superdirs)
5598 self.store_dirs.extend(c.gf_store_dirs)
5600 self._check_store_dirs_type()
5601 self._id_to_store_dir = {}
5602 self._open_stores = {}
5603 self._effective_default_store_id = None
5605 def _check_store_dirs_type(self):
5606 for sdir in ['store_dirs', 'store_superdirs']:
5607 if not isinstance(self.__getattribute__(sdir), list):
5608 raise TypeError('{} of {} is not of type list'.format(
5609 sdir, self.__class__.__name__))
5611 def _get_store_id(self, store_dir):
5612 store_ = store.Store(store_dir)
5613 store_id = store_.config.id
5614 store_.close()
5615 return store_id
5617 def _looks_like_store_dir(self, store_dir):
5618 return os.path.isdir(store_dir) and \
5619 all(os.path.isfile(pjoin(store_dir, x)) for x in
5620 ('index', 'traces', 'config'))
5622 def iter_store_dirs(self):
5623 store_dirs = set()
5624 for d in self.store_superdirs:
5625 if not os.path.exists(d):
5626 logger.warning('store_superdir not available: %s' % d)
5627 continue
5629 for entry in os.listdir(d):
5630 store_dir = os.path.realpath(pjoin(d, entry))
5631 if self._looks_like_store_dir(store_dir):
5632 store_dirs.add(store_dir)
5634 for store_dir in self.store_dirs:
5635 store_dirs.add(os.path.realpath(store_dir))
5637 return store_dirs
5639 def _scan_stores(self):
5640 for store_dir in self.iter_store_dirs():
5641 store_id = self._get_store_id(store_dir)
5642 if store_id not in self._id_to_store_dir:
5643 self._id_to_store_dir[store_id] = store_dir
5644 else:
5645 if store_dir != self._id_to_store_dir[store_id]:
5646 raise DuplicateStoreId(
5647 'GF store ID %s is used in (at least) two '
5648 'different stores. Locations are: %s and %s' %
5649 (store_id, self._id_to_store_dir[store_id], store_dir))
5651 def get_store_dir(self, store_id):
5652 '''
5653 Lookup directory given a GF store ID.
5654 '''
5656 if store_id not in self._id_to_store_dir:
5657 self._scan_stores()
5659 if store_id not in self._id_to_store_dir:
5660 raise NoSuchStore(
5661 store_id,
5662 self.store_dirs,
5663 self.store_superdirs,
5664 sorted(self._id_to_store_dir.keys()))
5666 return self._id_to_store_dir[store_id]
5668 def get_store_ids(self):
5669 '''
5670 Get list of available store IDs.
5671 '''
5673 self._scan_stores()
5674 return sorted(self._id_to_store_dir.keys())
5676 def effective_default_store_id(self):
5677 if self._effective_default_store_id is None:
5678 if self.default_store_id is None:
5679 store_ids = self.get_store_ids()
5680 if len(store_ids) == 1:
5681 self._effective_default_store_id = self.get_store_ids()[0]
5682 else:
5683 raise NoDefaultStoreSet()
5684 else:
5685 self._effective_default_store_id = self.default_store_id
5687 return self._effective_default_store_id
5689 def get_store(self, store_id=None):
5690 '''
5691 Get a store from the engine.
5693 :param store_id: identifier of the store (optional)
5694 :returns: :py:class:`~pyrocko.gf.store.Store` object
5696 If no ``store_id`` is provided the store
5697 associated with the :py:gattr:`default_store_id` is returned.
5698 Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is
5699 undefined.
5700 '''
5702 if store_id is None:
5703 store_id = self.effective_default_store_id()
5705 if store_id not in self._open_stores:
5706 store_dir = self.get_store_dir(store_id)
5707 self._open_stores[store_id] = store.Store(store_dir)
5709 return self._open_stores[store_id]
5711 def get_store_config(self, store_id):
5712 store = self.get_store(store_id)
5713 return store.config
5715 def get_store_extra(self, store_id, key):
5716 store = self.get_store(store_id)
5717 return store.get_extra(key)
5719 def close_cashed_stores(self):
5720 '''
5721 Close and remove ids from cashed stores.
5722 '''
5723 store_ids = []
5724 for store_id, store_ in self._open_stores.items():
5725 store_.close()
5726 store_ids.append(store_id)
5728 for store_id in store_ids:
5729 self._open_stores.pop(store_id)
5731 def get_rule(self, source, target):
5732 cprovided = self.get_store(target.store_id).get_provided_components()
5734 if isinstance(target, StaticTarget):
5735 quantity = target.quantity
5736 available_rules = static_rules
5737 elif isinstance(target, Target):
5738 quantity = target.effective_quantity()
5739 available_rules = channel_rules
5741 try:
5742 for rule in available_rules[quantity]:
5743 cneeded = rule.required_components(target)
5744 if all(c in cprovided for c in cneeded):
5745 return rule
5747 except KeyError:
5748 pass
5750 raise BadRequest(
5751 'No rule to calculate "%s" with GFs from store "%s" '
5752 'for source model "%s".' % (
5753 target.effective_quantity(),
5754 target.store_id,
5755 source.__class__.__name__))
5757 def _cached_discretize_basesource(self, source, store, cache, target):
5758 if (source, store) not in cache:
5759 cache[source, store] = source.discretize_basesource(store, target)
5761 return cache[source, store]
5763 def base_seismograms(self, source, targets, components, dsource_cache,
5764 nthreads=None):
5766 target = targets[0]
5768 interp = set([t.interpolation for t in targets])
5769 if len(interp) > 1:
5770 raise BadRequest('Targets have different interpolation schemes.')
5772 rates = set([t.sample_rate for t in targets])
5773 if len(rates) > 1:
5774 raise BadRequest('Targets have different sample rates.')
5776 store_ = self.get_store(target.store_id)
5777 receivers = [t.receiver(store_) for t in targets]
5779 if target.sample_rate is not None:
5780 deltat = 1. / target.sample_rate
5781 rate = target.sample_rate
5782 else:
5783 deltat = None
5784 rate = store_.config.sample_rate
5786 tmin = num.fromiter(
5787 (t.tmin for t in targets), dtype=float, count=len(targets))
5788 tmax = num.fromiter(
5789 (t.tmax for t in targets), dtype=float, count=len(targets))
5791 mask = num.logical_and(num.isfinite(tmin), num.isfinite(tmax))
5793 itmin = num.zeros_like(tmin, dtype=num.int64)
5794 itmax = num.zeros_like(tmin, dtype=num.int64)
5795 nsamples = num.full_like(tmin, -1, dtype=num.int64)
5797 itmin[mask] = num.floor(tmin[mask] * rate).astype(num.int64)
5798 itmax[mask] = num.ceil(tmax[mask] * rate).astype(num.int64)
5799 nsamples = itmax - itmin + 1
5800 nsamples[num.logical_not(mask)] = -1
5802 base_source = self._cached_discretize_basesource(
5803 source, store_, dsource_cache, target)
5805 base_seismograms = store_.calc_seismograms(
5806 base_source, receivers, components,
5807 deltat=deltat,
5808 itmin=itmin, nsamples=nsamples,
5809 interpolation=target.interpolation,
5810 optimization=target.optimization,
5811 nthreads=nthreads if nthreads is not None else 1)
5813 for i, base_seismogram in enumerate(base_seismograms):
5814 base_seismograms[i] = store.make_same_span(base_seismogram)
5816 return base_seismograms
5818 def base_seismogram(self, source, target, components, dsource_cache,
5819 nthreads):
5821 tcounters = [xtime()]
5823 store_ = self.get_store(target.store_id)
5824 receiver = target.receiver(store_)
5826 if target.tmin and target.tmax is not None:
5827 rate = store_.config.sample_rate
5828 itmin = int(num.floor(target.tmin * rate))
5829 itmax = int(num.ceil(target.tmax * rate))
5830 nsamples = itmax - itmin + 1
5831 else:
5832 itmin = None
5833 nsamples = None
5835 tcounters.append(xtime())
5836 base_source = self._cached_discretize_basesource(
5837 source, store_, dsource_cache, target)
5839 tcounters.append(xtime())
5841 if target.sample_rate is not None:
5842 deltat = 1. / target.sample_rate
5843 else:
5844 deltat = None
5846 base_seismogram = store_.seismogram(
5847 base_source, receiver, components,
5848 deltat=deltat,
5849 itmin=itmin, nsamples=nsamples,
5850 interpolation=target.interpolation,
5851 optimization=target.optimization,
5852 nthreads=nthreads)
5854 tcounters.append(xtime())
5856 base_seismogram = store.make_same_span(base_seismogram)
5858 tcounters.append(xtime())
5860 return base_seismogram, tcounters
5862 def base_statics(self, source, target, components, nthreads):
5863 tcounters = [xtime()]
5864 store_ = self.get_store(target.store_id)
5866 if target.tsnapshot is not None:
5867 rate = store_.config.sample_rate
5868 itsnapshot = int(num.floor(target.tsnapshot * rate))
5869 else:
5870 itsnapshot = None
5871 tcounters.append(xtime())
5873 base_source = source.discretize_basesource(store_, target=target)
5875 tcounters.append(xtime())
5877 base_statics = store_.statics(
5878 base_source,
5879 target,
5880 itsnapshot,
5881 components,
5882 target.interpolation,
5883 nthreads)
5885 tcounters.append(xtime())
5887 return base_statics, tcounters
5889 def _post_process_dynamic(self, base_seismogram, source, target):
5890 base_any = next(iter(base_seismogram.values()))
5891 deltat = base_any.deltat
5892 itmin = base_any.itmin
5894 rule = self.get_rule(source, target)
5895 data = rule.apply_(target, base_seismogram)
5897 factor = source.get_factor() * target.get_factor()
5898 if factor != 1.0:
5899 data = data * factor
5901 stf = source.effective_stf_post()
5903 times, amplitudes = stf.discretize_t(
5904 deltat, 0.0)
5906 # repeat end point to prevent boundary effects
5907 padded_data = num.empty(data.size + amplitudes.size, dtype=float)
5908 padded_data[:data.size] = data
5909 padded_data[data.size:] = data[-1]
5910 data = num.convolve(amplitudes, padded_data)
5912 tmin = itmin * deltat + times[0]
5914 tr = meta.SeismosizerTrace(
5915 codes=target.codes,
5916 data=data[:-amplitudes.size],
5917 deltat=deltat,
5918 tmin=tmin)
5920 return target.post_process(self, source, tr)
5922 def _post_process_statics(self, base_statics, source, starget):
5923 rule = self.get_rule(source, starget)
5924 data = rule.apply_(starget, base_statics)
5926 factor = source.get_factor()
5927 if factor != 1.0:
5928 for v in data.values():
5929 v *= factor
5931 return starget.post_process(self, source, base_statics)
5933 def process(self, *args, **kwargs):
5934 '''
5935 Process a request.
5937 ::
5939 process(**kwargs)
5940 process(request, **kwargs)
5941 process(sources, targets, **kwargs)
5943 The request can be given a a :py:class:`Request` object, or such an
5944 object is created using ``Request(**kwargs)`` for convenience.
5946 :returns: :py:class:`Response` object
5947 '''
5949 if len(args) not in (0, 1, 2):
5950 raise BadRequest('Invalid arguments.')
5952 if len(args) == 1:
5953 kwargs['request'] = args[0]
5955 elif len(args) == 2:
5956 kwargs.update(Request.args2kwargs(args))
5958 request = kwargs.pop('request', None)
5959 status_callback = kwargs.pop('status_callback', None)
5960 calc_timeseries = kwargs.pop('calc_timeseries', True)
5962 nprocs = kwargs.pop('nprocs', None)
5963 if nprocs is not None:
5964 raise BadRequest(
5965 'The `nprocs` keyword argument to process is no longer '
5966 'available. Please use `nthreads`.')
5968 nthreads = kwargs.pop('nthreads', self.nthreads)
5970 if request is None:
5971 request = Request(**kwargs)
5973 if resource:
5974 rs0 = resource.getrusage(resource.RUSAGE_SELF)
5975 rc0 = resource.getrusage(resource.RUSAGE_CHILDREN)
5976 tt0 = xtime()
5978 # make sure stores are open before fork()
5979 store_ids = set(target.store_id for target in request.targets)
5980 for store_id in store_ids:
5981 self.get_store(store_id)
5983 source_index = dict((x, i) for (i, x) in
5984 enumerate(request.sources))
5985 target_index = dict((x, i) for (i, x) in
5986 enumerate(request.targets))
5988 m = request.subrequest_map()
5990 skeys = sorted(m.keys(), key=cmp_to_key(cmp_none_aware))
5991 results_list = []
5993 for i in range(len(request.sources)):
5994 results_list.append([None] * len(request.targets))
5996 tcounters_dyn_list = []
5997 tcounters_static_list = []
5998 nsub = len(skeys)
5999 isub = 0
6001 # Processing dynamic targets through
6002 # parimap(process_subrequest_dynamic)
6004 if calc_timeseries:
6005 _process_dynamic = process_dynamic_timeseries
6006 else:
6007 _process_dynamic = process_dynamic
6009 if request.has_dynamic:
6010 work_dynamic = [
6011 (i, nsub,
6012 [source_index[source] for source in m[k][0]],
6013 [target_index[target] for target in m[k][1]
6014 if not isinstance(target, StaticTarget)])
6015 for (i, k) in enumerate(skeys)]
6017 for ii_results, tcounters_dyn in _process_dynamic(
6018 work_dynamic, request.sources, request.targets, self,
6019 nthreads):
6021 tcounters_dyn_list.append(num.diff(tcounters_dyn))
6022 isource, itarget, result = ii_results
6023 results_list[isource][itarget] = result
6025 if status_callback:
6026 status_callback(isub, nsub)
6028 isub += 1
6030 # Processing static targets through process_static
6031 if request.has_statics:
6032 work_static = [
6033 (i, nsub,
6034 [source_index[source] for source in m[k][0]],
6035 [target_index[target] for target in m[k][1]
6036 if isinstance(target, StaticTarget)])
6037 for (i, k) in enumerate(skeys)]
6039 for ii_results, tcounters_static in process_static(
6040 work_static, request.sources, request.targets, self,
6041 nthreads=nthreads):
6043 tcounters_static_list.append(num.diff(tcounters_static))
6044 isource, itarget, result = ii_results
6045 results_list[isource][itarget] = result
6047 if status_callback:
6048 status_callback(isub, nsub)
6050 isub += 1
6052 if status_callback:
6053 status_callback(nsub, nsub)
6055 tt1 = time.time()
6056 if resource:
6057 rs1 = resource.getrusage(resource.RUSAGE_SELF)
6058 rc1 = resource.getrusage(resource.RUSAGE_CHILDREN)
6060 s = ProcessingStats()
6062 if request.has_dynamic:
6063 tcumu_dyn = num.sum(num.vstack(tcounters_dyn_list), axis=0)
6064 t_dyn = float(num.sum(tcumu_dyn))
6065 perc_dyn = map(float, tcumu_dyn / t_dyn * 100.)
6066 (s.t_perc_get_store_and_receiver,
6067 s.t_perc_discretize_source,
6068 s.t_perc_make_base_seismogram,
6069 s.t_perc_make_same_span,
6070 s.t_perc_post_process) = perc_dyn
6071 else:
6072 t_dyn = 0.
6074 if request.has_statics:
6075 tcumu_static = num.sum(num.vstack(tcounters_static_list), axis=0)
6076 t_static = num.sum(tcumu_static)
6077 perc_static = map(float, tcumu_static / t_static * 100.)
6078 (s.t_perc_static_get_store,
6079 s.t_perc_static_discretize_basesource,
6080 s.t_perc_static_sum_statics,
6081 s.t_perc_static_post_process) = perc_static
6083 s.t_wallclock = tt1 - tt0
6084 if resource:
6085 s.t_cpu = (
6086 (rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) -
6087 (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime))
6088 s.n_read_blocks = (
6089 (rs1.ru_inblock + rc1.ru_inblock) -
6090 (rs0.ru_inblock + rc0.ru_inblock))
6092 n_records_stacked = 0.
6093 for results in results_list:
6094 for result in results:
6095 if not isinstance(result, meta.Result):
6096 continue
6097 shr = float(result.n_shared_stacking)
6098 n_records_stacked += result.n_records_stacked / shr
6099 s.t_perc_optimize += result.t_optimize / shr
6100 s.t_perc_stack += result.t_stack / shr
6101 s.n_records_stacked = int(n_records_stacked)
6102 if t_dyn != 0.:
6103 s.t_perc_optimize /= t_dyn * 100
6104 s.t_perc_stack /= t_dyn * 100
6106 return Response(
6107 request=request,
6108 results_list=results_list,
6109 stats=s)
6112class RemoteEngine(Engine):
6113 '''
6114 Client for remote synthetic seismogram calculator.
6115 '''
6117 site = String.T(default=ws.g_default_site, optional=True)
6118 url = String.T(default=ws.g_url, optional=True)
6120 def process(self, request=None, status_callback=None, **kwargs):
6122 if request is None:
6123 request = Request(**kwargs)
6125 return ws.seismosizer(url=self.url, site=self.site, request=request)
6128g_engine = None
6131def get_engine(store_superdirs=[]):
6132 global g_engine
6133 if g_engine is None:
6134 g_engine = LocalEngine(use_env=True, use_config=True)
6136 for d in store_superdirs:
6137 if d not in g_engine.store_superdirs:
6138 g_engine.store_superdirs.append(d)
6140 return g_engine
6143class SourceGroup(Object):
6145 def __getattr__(self, k):
6146 return num.fromiter((getattr(s, k) for s in self),
6147 dtype=float)
6149 def __iter__(self):
6150 raise NotImplementedError(
6151 'This method should be implemented in subclass.')
6153 def __len__(self):
6154 raise NotImplementedError(
6155 'This method should be implemented in subclass.')
6158class SourceList(SourceGroup):
6159 sources = List.T(Source.T())
6161 def append(self, s):
6162 self.sources.append(s)
6164 def __iter__(self):
6165 return iter(self.sources)
6167 def __len__(self):
6168 return len(self.sources)
6171class SourceGrid(SourceGroup):
6173 base = Source.T()
6174 variables = Dict.T(String.T(), Range.T())
6175 order = List.T(String.T())
6177 def __len__(self):
6178 n = 1
6179 for (k, v) in self.make_coords(self.base):
6180 n *= len(list(v))
6182 return n
6184 def __iter__(self):
6185 for items in permudef(self.make_coords(self.base)):
6186 s = self.base.clone(**{k: v for (k, v) in items})
6187 s.regularize()
6188 yield s
6190 def ordered_params(self):
6191 ks = list(self.variables.keys())
6192 for k in self.order + list(self.base.keys()):
6193 if k in ks:
6194 yield k
6195 ks.remove(k)
6196 if ks:
6197 raise Exception('Invalid parameter "%s" for source type "%s".' %
6198 (ks[0], self.base.__class__.__name__))
6200 def make_coords(self, base):
6201 return [(param, self.variables[param].make(base=base[param]))
6202 for param in self.ordered_params()]
6205source_classes = [
6206 Source,
6207 SourceWithMagnitude,
6208 SourceWithDerivedMagnitude,
6209 ExplosionSource,
6210 RectangularExplosionSource,
6211 DCSource,
6212 CLVDSource,
6213 VLVDSource,
6214 MTSource,
6215 RectangularSource,
6216 PseudoDynamicRupture,
6217 DoubleDCSource,
6218 RingfaultSource,
6219 CombiSource,
6220 CombiSFSource,
6221 SFSource,
6222 SimpleLandslideSource,
6223 PorePressurePointSource,
6224 PorePressureLineSource,
6225]
6227stf_classes = [
6228 STF,
6229 BoxcarSTF,
6230 TriangularSTF,
6231 HalfSinusoidSTF,
6232 ResonatorSTF,
6233 TremorSTF,
6234 SimpleLandslideSTF,
6235]
6237__all__ = '''
6238Cloneable
6239NoDefaultStoreSet
6240SeismosizerError
6241STFError
6242BadRequest
6243NoSuchStore
6244DerivedMagnitudeError
6245STFMode
6246'''.split() + [S.__name__ for S in source_classes + stf_classes] + '''
6247Request
6248ProcessingStats
6249Response
6250Engine
6251LocalEngine
6252RemoteEngine
6253source_classes
6254get_engine
6255Range
6256SourceGroup
6257SourceList
6258SourceGrid
6259map_anchor
6260'''.split()