1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
10.. _coordinate-system-names:
12Coordinate systems
13..................
15Coordinate system names commonly used in source models.
17================= ============================================
18Name Description
19================= ============================================
20``'xyz'`` northing, easting, depth in [m]
21``'xy'`` northing, easting in [m]
22``'latlon'`` latitude, longitude in [deg]
23``'lonlat'`` longitude, latitude in [deg]
24``'latlondepth'`` latitude, longitude in [deg], depth in [m]
25================= ============================================
26'''
29from __future__ import absolute_import, division, print_function
31from collections import defaultdict
32from functools import cmp_to_key
33import time
34import math
35import os
36import re
37import logging
38try:
39 import resource
40except ImportError:
41 resource = None
42from hashlib import sha1
44import numpy as num
45from scipy.interpolate import RegularGridInterpolator
47from pyrocko.guts import (Object, Float, String, StringChoice, List,
48 Timestamp, Int, SObject, ArgumentError, Dict,
49 ValidationError, Bool)
50from pyrocko.guts_array import Array
52from pyrocko import moment_tensor as pmt
53from pyrocko import trace, util, config, model, eikonal_ext
54from pyrocko.orthodrome import ne_to_latlon
55from pyrocko.model import Location
56from pyrocko.modelling import OkadaSource, make_okada_coefficient_matrix, \
57 okada_ext, invert_fault_dislocations_bem
59from . import meta, store, ws
60from .tractions import TractionField, DirectedTractions
61from .targets import Target, StaticTarget, SatelliteTarget
63pjoin = os.path.join
65guts_prefix = 'pf'
67d2r = math.pi / 180.
68r2d = 180. / math.pi
69km = 1e3
71logger = logging.getLogger('pyrocko.gf.seismosizer')
74def cmp_none_aware(a, b):
75 if isinstance(a, tuple) and isinstance(b, tuple):
76 for xa, xb in zip(a, b):
77 rv = cmp_none_aware(xa, xb)
78 if rv != 0:
79 return rv
81 return 0
83 anone = a is None
84 bnone = b is None
86 if anone and bnone:
87 return 0
89 if anone:
90 return -1
92 if bnone:
93 return 1
95 return bool(a > b) - bool(a < b)
98def xtime():
99 return time.time()
102class SeismosizerError(Exception):
103 pass
106class BadRequest(SeismosizerError):
107 pass
110class DuplicateStoreId(Exception):
111 pass
114class NoDefaultStoreSet(Exception):
115 pass
118class ConversionError(Exception):
119 pass
122class NoSuchStore(BadRequest):
124 def __init__(self, store_id=None, dirs=None):
125 BadRequest.__init__(self)
126 self.store_id = store_id
127 self.dirs = dirs
129 def __str__(self):
130 if self.store_id is not None:
131 rstr = 'no GF store with id "%s" found.' % self.store_id
132 else:
133 rstr = 'GF store not found.'
135 if self.dirs is not None:
136 rstr += ' Searched folders:\n %s' % '\n '.join(sorted(self.dirs))
137 return rstr
140def ufloat(s):
141 units = {
142 'k': 1e3,
143 'M': 1e6,
144 }
146 factor = 1.0
147 if s and s[-1] in units:
148 factor = units[s[-1]]
149 s = s[:-1]
150 if not s:
151 raise ValueError('unit without a number: \'%s\'' % s)
153 return float(s) * factor
156def ufloat_or_none(s):
157 if s:
158 return ufloat(s)
159 else:
160 return None
163def int_or_none(s):
164 if s:
165 return int(s)
166 else:
167 return None
170def nonzero(x, eps=1e-15):
171 return abs(x) > eps
174def permudef(ln, j=0):
175 if j < len(ln):
176 k, v = ln[j]
177 for y in v:
178 ln[j] = k, y
179 for s in permudef(ln, j + 1):
180 yield s
182 ln[j] = k, v
183 return
184 else:
185 yield ln
188def arr(x):
189 return num.atleast_1d(num.asarray(x))
192def discretize_rect_source(deltas, deltat, time, north, east, depth,
193 strike, dip, length, width,
194 anchor, velocity=None, stf=None,
195 nucleation_x=None, nucleation_y=None,
196 decimation_factor=1, pointsonly=False,
197 plane_coords=False,
198 aggressive_oversampling=False):
200 if stf is None:
201 stf = STF()
203 if not velocity and not pointsonly:
204 raise AttributeError('velocity is required in time mode')
206 mindeltagf = float(num.min(deltas))
207 if velocity:
208 mindeltagf = min(mindeltagf, deltat * velocity)
210 ln = length
211 wd = width
213 if aggressive_oversampling:
214 nl = int((2. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
215 nw = int((2. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
216 else:
217 nl = int((1. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
218 nw = int((1. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
220 n = int(nl * nw)
222 dl = ln / nl
223 dw = wd / nw
225 xl = num.linspace(-0.5 * (ln - dl), 0.5 * (ln - dl), nl)
226 xw = num.linspace(-0.5 * (wd - dw), 0.5 * (wd - dw), nw)
228 points = num.zeros((n, 3))
229 points[:, 0] = num.tile(xl, nw)
230 points[:, 1] = num.repeat(xw, nl)
232 if nucleation_x is not None:
233 dist_x = num.abs(nucleation_x - points[:, 0])
234 else:
235 dist_x = num.zeros(n)
237 if nucleation_y is not None:
238 dist_y = num.abs(nucleation_y - points[:, 1])
239 else:
240 dist_y = num.zeros(n)
242 dist = num.sqrt(dist_x**2 + dist_y**2)
243 times = dist / velocity
245 anch_x, anch_y = map_anchor[anchor]
247 points[:, 0] -= anch_x * 0.5 * length
248 points[:, 1] -= anch_y * 0.5 * width
250 if plane_coords:
251 return points, dl, dw, nl, nw
253 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
254 points = num.dot(rotmat.T, points.T).T
256 points[:, 0] += north
257 points[:, 1] += east
258 points[:, 2] += depth
260 if pointsonly:
261 return points, dl, dw, nl, nw
263 xtau, amplitudes = stf.discretize_t(deltat, time)
264 nt = xtau.size
266 points2 = num.repeat(points, nt, axis=0)
267 times2 = (times[:, num.newaxis] + xtau[num.newaxis, :]).ravel()
268 amplitudes2 = num.tile(amplitudes, n)
270 return points2, times2, amplitudes2, dl, dw, nl, nw
273def check_rect_source_discretisation(points2, nl, nw, store):
274 # We assume a non-rotated fault plane
275 N_CRITICAL = 8
276 points = points2.T.reshape((3, nl, nw))
277 if points.size <= N_CRITICAL:
278 logger.warning('RectangularSource is defined by only %d sub-sources!'
279 % points.size)
280 return True
282 distances = num.sqrt(
283 (points[0, 0, :] - points[0, 1, :])**2 +
284 (points[1, 0, :] - points[1, 1, :])**2 +
285 (points[2, 0, :] - points[2, 1, :])**2)
287 depths = points[2, 0, :]
288 vs_profile = store.config.get_vs(
289 lat=0., lon=0.,
290 points=num.repeat(depths[:, num.newaxis], 3, axis=1),
291 interpolation='multilinear')
293 min_wavelength = vs_profile * (store.config.deltat * 2)
294 if not num.all(min_wavelength > distances / 2):
295 return False
296 return True
299def outline_rect_source(strike, dip, length, width, anchor):
300 ln = length
301 wd = width
302 points = num.array(
303 [[-0.5 * ln, -0.5 * wd, 0.],
304 [0.5 * ln, -0.5 * wd, 0.],
305 [0.5 * ln, 0.5 * wd, 0.],
306 [-0.5 * ln, 0.5 * wd, 0.],
307 [-0.5 * ln, -0.5 * wd, 0.]])
309 anch_x, anch_y = map_anchor[anchor]
310 points[:, 0] -= anch_x * 0.5 * length
311 points[:, 1] -= anch_y * 0.5 * width
313 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
315 return num.dot(rotmat.T, points.T).T
318def from_plane_coords(
319 strike, dip, length, width, depth, x_plane_coords, y_plane_coords,
320 lat=0., lon=0.,
321 north_shift=0, east_shift=0,
322 anchor='top', cs='xy'):
324 ln = length
325 wd = width
326 x_abs = []
327 y_abs = []
328 if not isinstance(x_plane_coords, list):
329 x_plane_coords = [x_plane_coords]
330 y_plane_coords = [y_plane_coords]
332 for x_plane, y_plane in zip(x_plane_coords, y_plane_coords):
333 points = num.array(
334 [[-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
335 [0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
336 [0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
337 [-0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
338 [-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.]])
340 anch_x, anch_y = map_anchor[anchor]
341 points[:, 0] -= anch_x * 0.5 * length
342 points[:, 1] -= anch_y * 0.5 * width
344 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
346 points = num.dot(rotmat.T, points.T).T
347 points[:, 0] += north_shift
348 points[:, 1] += east_shift
349 points[:, 2] += depth
350 if cs in ('latlon', 'lonlat'):
351 latlon = ne_to_latlon(lat, lon,
352 points[:, 0], points[:, 1])
353 latlon = num.array(latlon).T
354 x_abs.append(latlon[1:2, 1])
355 y_abs.append(latlon[2:3, 0])
356 if cs == 'xy':
357 x_abs.append(points[1:2, 1])
358 y_abs.append(points[2:3, 0])
360 if cs == 'lonlat':
361 return y_abs, x_abs
362 else:
363 return x_abs, y_abs
366def points_on_rect_source(
367 strike, dip, length, width, anchor,
368 discretized_basesource=None, points_x=None, points_y=None):
370 ln = length
371 wd = width
373 if isinstance(points_x, list) or isinstance(points_x, float):
374 points_x = num.array([points_x])
375 if isinstance(points_y, list) or isinstance(points_y, float):
376 points_y = num.array([points_y])
378 if discretized_basesource:
379 ds = discretized_basesource
381 nl_patches = ds.nl + 1
382 nw_patches = ds.nw + 1
384 npoints = nl_patches * nw_patches
385 points = num.zeros((npoints, 3))
386 ln_patches = num.array([il for il in range(nl_patches)])
387 wd_patches = num.array([iw for iw in range(nw_patches)])
389 points_ln =\
390 2 * ((ln_patches - num.min(ln_patches)) / num.ptp(ln_patches)) - 1
391 points_wd =\
392 2 * ((wd_patches - num.min(wd_patches)) / num.ptp(wd_patches)) - 1
394 for il in range(nl_patches):
395 for iw in range(nw_patches):
396 points[il * nw_patches + iw, :] = num.array([
397 points_ln[il] * ln * 0.5,
398 points_wd[iw] * wd * 0.5, 0.0])
400 elif points_x.shape[0] > 0 and points_y.shape[0] > 0:
401 points = num.zeros(shape=((len(points_x), 3)))
402 for i, (x, y) in enumerate(zip(points_x, points_y)):
403 points[i, :] = num.array(
404 [x * 0.5 * ln, y * 0.5 * wd, 0.0])
406 anch_x, anch_y = map_anchor[anchor]
408 points[:, 0] -= anch_x * 0.5 * ln
409 points[:, 1] -= anch_y * 0.5 * wd
411 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
413 return num.dot(rotmat.T, points.T).T
416class InvalidGridDef(Exception):
417 pass
420class Range(SObject):
421 '''
422 Convenient range specification.
424 Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]::
426 Range('0 .. 10k : 1k')
427 Range(start=0., stop=10e3, step=1e3)
428 Range(0, 10e3, 1e3)
429 Range('0 .. 10k @ 11')
430 Range(start=0., stop=10*km, n=11)
432 Range(0, 10e3, n=11)
433 Range(values=[x*1e3 for x in range(11)])
435 Depending on the use context, it can be possible to omit any part of the
436 specification. E.g. in the context of extracting a subset of an already
437 existing range, the existing range's specification values would be filled
438 in where missing.
440 The values are distributed with equal spacing, unless the ``spacing``
441 argument is modified. The values can be created offset or relative to an
442 external base value with the ``relative`` argument if the use context
443 supports this.
445 The range specification can be expressed with a short string
446 representation::
448 'start .. stop @ num | spacing, relative'
449 'start .. stop : step | spacing, relative'
451 most parts of the expression can be omitted if not needed. Whitespace is
452 allowed for readability but can also be omitted.
453 '''
455 start = Float.T(optional=True)
456 stop = Float.T(optional=True)
457 step = Float.T(optional=True)
458 n = Int.T(optional=True)
459 values = Array.T(optional=True, dtype=float, shape=(None,))
461 spacing = StringChoice.T(
462 choices=['lin', 'log', 'symlog'],
463 default='lin',
464 optional=True)
466 relative = StringChoice.T(
467 choices=['', 'add', 'mult'],
468 default='',
469 optional=True)
471 pattern = re.compile(r'^((?P<start>.*)\.\.(?P<stop>[^@|:]*))?'
472 r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?'
473 r'(\|(?P<stuff>.+))?$')
475 def __init__(self, *args, **kwargs):
476 d = {}
477 if len(args) == 1:
478 d = self.parse(args[0])
479 elif len(args) in (2, 3):
480 d['start'], d['stop'] = [float(x) for x in args[:2]]
481 if len(args) == 3:
482 d['step'] = float(args[2])
484 for k, v in kwargs.items():
485 if k in d:
486 raise ArgumentError('%s specified more than once' % k)
488 d[k] = v
490 SObject.__init__(self, **d)
492 def __str__(self):
493 def sfloat(x):
494 if x is not None:
495 return '%g' % x
496 else:
497 return ''
499 if self.values:
500 return ','.join('%g' % x for x in self.values)
502 if self.start is None and self.stop is None:
503 s0 = ''
504 else:
505 s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop))
507 s1 = ''
508 if self.step is not None:
509 s1 = [' : %g', ':%g'][s0 == ''] % self.step
510 elif self.n is not None:
511 s1 = [' @ %i', '@%i'][s0 == ''] % self.n
513 if self.spacing == 'lin' and self.relative == '':
514 s2 = ''
515 else:
516 x = []
517 if self.spacing != 'lin':
518 x.append(self.spacing)
520 if self.relative != '':
521 x.append(self.relative)
523 s2 = ' | %s' % ','.join(x)
525 return s0 + s1 + s2
527 @classmethod
528 def parse(cls, s):
529 s = re.sub(r'\s+', '', s)
530 m = cls.pattern.match(s)
531 if not m:
532 try:
533 vals = [ufloat(x) for x in s.split(',')]
534 except Exception:
535 raise InvalidGridDef(
536 '"%s" is not a valid range specification' % s)
538 return dict(values=num.array(vals, dtype=float))
540 d = m.groupdict()
541 try:
542 start = ufloat_or_none(d['start'])
543 stop = ufloat_or_none(d['stop'])
544 step = ufloat_or_none(d['step'])
545 n = int_or_none(d['n'])
546 except Exception:
547 raise InvalidGridDef(
548 '"%s" is not a valid range specification' % s)
550 spacing = 'lin'
551 relative = ''
553 if d['stuff'] is not None:
554 t = d['stuff'].split(',')
555 for x in t:
556 if x in cls.spacing.choices:
557 spacing = x
558 elif x and x in cls.relative.choices:
559 relative = x
560 else:
561 raise InvalidGridDef(
562 '"%s" is not a valid range specification' % s)
564 return dict(start=start, stop=stop, step=step, n=n, spacing=spacing,
565 relative=relative)
567 def make(self, mi=None, ma=None, inc=None, base=None, eps=1e-5):
568 if self.values:
569 return self.values
571 start = self.start
572 stop = self.stop
573 step = self.step
574 n = self.n
576 swap = step is not None and step < 0.
577 if start is None:
578 start = [mi, ma][swap]
579 if stop is None:
580 stop = [ma, mi][swap]
581 if step is None and inc is not None:
582 step = [inc, -inc][ma < mi]
584 if start is None or stop is None:
585 raise InvalidGridDef(
586 'Cannot use range specification "%s" without start '
587 'and stop in this context' % self)
589 if step is None and n is None:
590 step = stop - start
592 if n is None:
593 if (step < 0) != (stop - start < 0):
594 raise InvalidGridDef(
595 'Range specification "%s" has inconsistent ordering '
596 '(step < 0 => stop > start)' % self)
598 n = int(round((stop - start) / step)) + 1
599 stop2 = start + (n - 1) * step
600 if abs(stop - stop2) > eps:
601 n = int(math.floor((stop - start) / step)) + 1
602 stop = start + (n - 1) * step
603 else:
604 stop = stop2
606 if start == stop:
607 n = 1
609 if self.spacing == 'lin':
610 vals = num.linspace(start, stop, n)
612 elif self.spacing in ('log', 'symlog'):
613 if start > 0. and stop > 0.:
614 vals = num.exp(num.linspace(num.log(start),
615 num.log(stop), n))
616 elif start < 0. and stop < 0.:
617 vals = -num.exp(num.linspace(num.log(-start),
618 num.log(-stop), n))
619 else:
620 raise InvalidGridDef(
621 'Log ranges should not include or cross zero '
622 '(in range specification "%s").' % self)
624 if self.spacing == 'symlog':
625 nvals = - vals
626 vals = num.concatenate((nvals[::-1], vals))
628 if self.relative in ('add', 'mult') and base is None:
629 raise InvalidGridDef(
630 'Cannot use relative range specification in this context.')
632 vals = self.make_relative(base, vals)
634 return list(map(float, vals))
636 def make_relative(self, base, vals):
637 if self.relative == 'add':
638 vals += base
640 if self.relative == 'mult':
641 vals *= base
643 return vals
646class GridDefElement(Object):
648 param = meta.StringID.T()
649 rs = Range.T()
651 def __init__(self, shorthand=None, **kwargs):
652 if shorthand is not None:
653 t = shorthand.split('=')
654 if len(t) != 2:
655 raise InvalidGridDef(
656 'Invalid grid specification element: %s' % shorthand)
658 sp, sr = t[0].strip(), t[1].strip()
660 kwargs['param'] = sp
661 kwargs['rs'] = Range(sr)
663 Object.__init__(self, **kwargs)
665 def shorthand(self):
666 return self.param + ' = ' + str(self.rs)
669class GridDef(Object):
671 elements = List.T(GridDefElement.T())
673 def __init__(self, shorthand=None, **kwargs):
674 if shorthand is not None:
675 t = shorthand.splitlines()
676 tt = []
677 for x in t:
678 x = x.strip()
679 if x:
680 tt.extend(x.split(';'))
682 elements = []
683 for se in tt:
684 elements.append(GridDef(se))
686 kwargs['elements'] = elements
688 Object.__init__(self, **kwargs)
690 def shorthand(self):
691 return '; '.join(str(x) for x in self.elements)
694class Cloneable(object):
696 def __iter__(self):
697 return iter(self.T.propnames)
699 def __getitem__(self, k):
700 if k not in self.keys():
701 raise KeyError(k)
703 return getattr(self, k)
705 def __setitem__(self, k, v):
706 if k not in self.keys():
707 raise KeyError(k)
709 return setattr(self, k, v)
711 def clone(self, **kwargs):
712 '''
713 Make a copy of the object.
715 A new object of the same class is created and initialized with the
716 parameters of the object on which this method is called on. If
717 ``kwargs`` are given, these are used to override any of the
718 initialization parameters.
719 '''
721 d = dict(self)
722 for k in d:
723 v = d[k]
724 if isinstance(v, Cloneable):
725 d[k] = v.clone()
727 d.update(kwargs)
728 return self.__class__(**d)
730 @classmethod
731 def keys(cls):
732 '''
733 Get list of the source model's parameter names.
734 '''
736 return cls.T.propnames
739class STF(Object, Cloneable):
741 '''
742 Base class for source time functions.
743 '''
745 def __init__(self, effective_duration=None, **kwargs):
746 if effective_duration is not None:
747 kwargs['duration'] = effective_duration / \
748 self.factor_duration_to_effective()
750 Object.__init__(self, **kwargs)
752 @classmethod
753 def factor_duration_to_effective(cls):
754 return 1.0
756 def centroid_time(self, tref):
757 return tref
759 @property
760 def effective_duration(self):
761 return self.duration * self.factor_duration_to_effective()
763 def discretize_t(self, deltat, tref):
764 tl = math.floor(tref / deltat) * deltat
765 th = math.ceil(tref / deltat) * deltat
766 if tl == th:
767 return num.array([tl], dtype=float), num.ones(1)
768 else:
769 return (
770 num.array([tl, th], dtype=float),
771 num.array([th - tref, tref - tl], dtype=float) / deltat)
773 def base_key(self):
774 return (type(self).__name__,)
777g_unit_pulse = STF()
780def sshift(times, amplitudes, tshift, deltat):
782 t0 = math.floor(tshift / deltat) * deltat
783 t1 = math.ceil(tshift / deltat) * deltat
784 if t0 == t1:
785 return times, amplitudes
787 amplitudes2 = num.zeros(amplitudes.size + 1, dtype=float)
789 amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes
790 amplitudes2[1:] += (tshift - t0) / deltat * amplitudes
792 times2 = num.arange(times.size + 1, dtype=float) * \
793 deltat + times[0] + t0
795 return times2, amplitudes2
798class BoxcarSTF(STF):
800 '''
801 Boxcar type source time function.
803 .. figure :: /static/stf-BoxcarSTF.svg
804 :width: 40%
805 :align: center
806 :alt: boxcar source time function
807 '''
809 duration = Float.T(
810 default=0.0,
811 help='duration of the boxcar')
813 anchor = Float.T(
814 default=0.0,
815 help='anchor point with respect to source.time: ('
816 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
817 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
818 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
820 @classmethod
821 def factor_duration_to_effective(cls):
822 return 1.0
824 def centroid_time(self, tref):
825 return tref - 0.5 * self.duration * self.anchor
827 def discretize_t(self, deltat, tref):
828 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
829 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
830 tmin = round(tmin_stf / deltat) * deltat
831 tmax = round(tmax_stf / deltat) * deltat
832 nt = int(round((tmax - tmin) / deltat)) + 1
833 times = num.linspace(tmin, tmax, nt)
834 amplitudes = num.ones_like(times)
835 if times.size > 1:
836 t_edges = num.linspace(
837 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
838 t = tmin_stf + self.duration * num.array(
839 [0.0, 0.0, 1.0, 1.0], dtype=float)
840 f = num.array([0., 1., 1., 0.], dtype=float)
841 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
842 amplitudes /= num.sum(amplitudes)
844 tshift = (num.sum(amplitudes * times) - self.centroid_time(tref))
846 return sshift(times, amplitudes, -tshift, deltat)
848 def base_key(self):
849 return (type(self).__name__, self.duration, self.anchor)
852class TriangularSTF(STF):
854 '''
855 Triangular type source time function.
857 .. figure :: /static/stf-TriangularSTF.svg
858 :width: 40%
859 :align: center
860 :alt: triangular source time function
861 '''
863 duration = Float.T(
864 default=0.0,
865 help='baseline of the triangle')
867 peak_ratio = Float.T(
868 default=0.5,
869 help='fraction of time compared to duration, '
870 'when the maximum amplitude is reached')
872 anchor = Float.T(
873 default=0.0,
874 help='anchor point with respect to source.time: ('
875 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
876 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
877 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
879 @classmethod
880 def factor_duration_to_effective(cls, peak_ratio=None):
881 if peak_ratio is None:
882 peak_ratio = cls.peak_ratio.default()
884 return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0)
886 def __init__(self, effective_duration=None, **kwargs):
887 if effective_duration is not None:
888 kwargs['duration'] = effective_duration / \
889 self.factor_duration_to_effective(
890 kwargs.get('peak_ratio', None))
892 STF.__init__(self, **kwargs)
894 @property
895 def centroid_ratio(self):
896 ra = self.peak_ratio
897 rb = 1.0 - ra
898 return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb)
900 def centroid_time(self, tref):
901 ca = self.centroid_ratio
902 cb = 1.0 - ca
903 if self.anchor <= 0.:
904 return tref - ca * self.duration * self.anchor
905 else:
906 return tref - cb * self.duration * self.anchor
908 @property
909 def effective_duration(self):
910 return self.duration * self.factor_duration_to_effective(
911 self.peak_ratio)
913 def tminmax_stf(self, tref):
914 ca = self.centroid_ratio
915 cb = 1.0 - ca
916 if self.anchor <= 0.:
917 tmin_stf = tref - ca * self.duration * (self.anchor + 1.)
918 tmax_stf = tmin_stf + self.duration
919 else:
920 tmax_stf = tref + cb * self.duration * (1. - self.anchor)
921 tmin_stf = tmax_stf - self.duration
923 return tmin_stf, tmax_stf
925 def discretize_t(self, deltat, tref):
926 tmin_stf, tmax_stf = self.tminmax_stf(tref)
928 tmin = round(tmin_stf / deltat) * deltat
929 tmax = round(tmax_stf / deltat) * deltat
930 nt = int(round((tmax - tmin) / deltat)) + 1
931 if nt > 1:
932 t_edges = num.linspace(
933 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
934 t = tmin_stf + self.duration * num.array(
935 [0.0, self.peak_ratio, 1.0], dtype=float)
936 f = num.array([0., 1., 0.], dtype=float)
937 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
938 amplitudes /= num.sum(amplitudes)
939 else:
940 amplitudes = num.ones(1)
942 times = num.linspace(tmin, tmax, nt)
943 return times, amplitudes
945 def base_key(self):
946 return (
947 type(self).__name__, self.duration, self.peak_ratio, self.anchor)
950class HalfSinusoidSTF(STF):
952 '''
953 Half sinusoid type source time function.
955 .. figure :: /static/stf-HalfSinusoidSTF.svg
956 :width: 40%
957 :align: center
958 :alt: half-sinusouid source time function
959 '''
961 duration = Float.T(
962 default=0.0,
963 help='duration of the half-sinusoid (baseline)')
965 anchor = Float.T(
966 default=0.0,
967 help='anchor point with respect to source.time: ('
968 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
969 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
970 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
972 exponent = Int.T(
973 default=1,
974 help='set to 2 to use square of the half-period sinusoidal function.')
976 def __init__(self, effective_duration=None, **kwargs):
977 if effective_duration is not None:
978 kwargs['duration'] = effective_duration / \
979 self.factor_duration_to_effective(
980 kwargs.get('exponent', 1))
982 STF.__init__(self, **kwargs)
984 @classmethod
985 def factor_duration_to_effective(cls, exponent):
986 if exponent == 1:
987 return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi
988 elif exponent == 2:
989 return math.sqrt(math.pi**2 - 6) / math.pi
990 else:
991 raise ValueError('Exponent for HalfSinusoidSTF must be 1 or 2.')
993 @property
994 def effective_duration(self):
995 return self.duration * self.factor_duration_to_effective(self.exponent)
997 def centroid_time(self, tref):
998 return tref - 0.5 * self.duration * self.anchor
1000 def discretize_t(self, deltat, tref):
1001 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1002 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1003 tmin = round(tmin_stf / deltat) * deltat
1004 tmax = round(tmax_stf / deltat) * deltat
1005 nt = int(round((tmax - tmin) / deltat)) + 1
1006 if nt > 1:
1007 t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace(
1008 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)))
1010 if self.exponent == 1:
1011 fint = -num.cos(
1012 (t_edges - tmin_stf) * (math.pi / self.duration))
1014 elif self.exponent == 2:
1015 fint = (t_edges - tmin_stf) / self.duration \
1016 - 1.0 / (2.0 * math.pi) * num.sin(
1017 (t_edges - tmin_stf) * (2.0 * math.pi / self.duration))
1018 else:
1019 raise ValueError(
1020 'Exponent for HalfSinusoidSTF must be 1 or 2.')
1022 amplitudes = fint[1:] - fint[:-1]
1023 amplitudes /= num.sum(amplitudes)
1024 else:
1025 amplitudes = num.ones(1)
1027 times = num.linspace(tmin, tmax, nt)
1028 return times, amplitudes
1030 def base_key(self):
1031 return (type(self).__name__, self.duration, self.anchor)
1034class SmoothRampSTF(STF):
1035 '''
1036 Smooth-ramp type source time function for near-field displacement.
1037 Based on moment function of double-couple point source proposed by Bruestle
1038 and Mueller (PEPI, 1983).
1040 .. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow
1041 earthquakes from Love-wave modelling for regional distances, PEPI 32,
1042 312-324.
1044 .. figure :: /static/stf-SmoothRampSTF.svg
1045 :width: 40%
1046 :alt: smooth ramp source time function
1047 '''
1048 duration = Float.T(
1049 default=0.0,
1050 help='duration of the ramp (baseline)')
1052 rise_ratio = Float.T(
1053 default=0.5,
1054 help='fraction of time compared to duration, '
1055 'when the maximum amplitude is reached')
1057 anchor = Float.T(
1058 default=0.0,
1059 help='anchor point with respect to source.time: ('
1060 '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, '
1061 '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, '
1062 '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)')
1064 def discretize_t(self, deltat, tref):
1065 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1066 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1067 tmin = round(tmin_stf / deltat) * deltat
1068 tmax = round(tmax_stf / deltat) * deltat
1069 D = round((tmax - tmin) / deltat) * deltat
1070 nt = int(round(D / deltat)) + 1
1071 times = num.linspace(tmin, tmax, nt)
1072 if nt > 1:
1073 rise_time = self.rise_ratio * self.duration
1074 amplitudes = num.ones_like(times)
1075 tp = tmin + rise_time
1076 ii = num.where(times <= tp)
1077 t_inc = times[ii]
1078 a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time)
1079 b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0
1080 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b)
1082 amplitudes /= num.sum(amplitudes)
1083 else:
1084 amplitudes = num.ones(1)
1086 return times, amplitudes
1088 def base_key(self):
1089 return (type(self).__name__,
1090 self.duration, self.rise_ratio, self.anchor)
1093class ResonatorSTF(STF):
1094 '''
1095 Simple resonator like source time function.
1097 .. math ::
1099 f(t) = 0 for t < 0
1100 f(t) = e^{-t/tau} * sin(2 * pi * f * t)
1103 .. figure :: /static/stf-SmoothRampSTF.svg
1104 :width: 40%
1105 :alt: smooth ramp source time function
1107 '''
1109 duration = Float.T(
1110 default=0.0,
1111 help='decay time')
1113 frequency = Float.T(
1114 default=1.0,
1115 help='resonance frequency')
1117 def discretize_t(self, deltat, tref):
1118 tmin_stf = tref
1119 tmax_stf = tref + self.duration * 3
1120 tmin = math.floor(tmin_stf / deltat) * deltat
1121 tmax = math.ceil(tmax_stf / deltat) * deltat
1122 times = util.arange2(tmin, tmax, deltat)
1123 amplitudes = num.exp(-(times - tref) / self.duration) \
1124 * num.sin(2.0 * num.pi * self.frequency * (times - tref))
1126 return times, amplitudes
1128 def base_key(self):
1129 return (type(self).__name__,
1130 self.duration, self.frequency)
1133class STFMode(StringChoice):
1134 choices = ['pre', 'post']
1137class Source(Location, Cloneable):
1138 '''
1139 Base class for all source models.
1140 '''
1142 name = String.T(optional=True, default='')
1144 time = Timestamp.T(
1145 default=Timestamp.D('1970-01-01 00:00:00'),
1146 help='source origin time.')
1148 stf = STF.T(
1149 optional=True,
1150 help='source time function.')
1152 stf_mode = STFMode.T(
1153 default='post',
1154 help='whether to apply source time function in pre or '
1155 'post-processing.')
1157 def __init__(self, **kwargs):
1158 Location.__init__(self, **kwargs)
1160 def update(self, **kwargs):
1161 '''
1162 Change some of the source models parameters.
1164 Example::
1166 >>> from pyrocko import gf
1167 >>> s = gf.DCSource()
1168 >>> s.update(strike=66., dip=33.)
1169 >>> print(s)
1170 --- !pf.DCSource
1171 depth: 0.0
1172 time: 1970-01-01 00:00:00
1173 magnitude: 6.0
1174 strike: 66.0
1175 dip: 33.0
1176 rake: 0.0
1178 '''
1180 for (k, v) in kwargs.items():
1181 self[k] = v
1183 def grid(self, **variables):
1184 '''
1185 Create grid of source model variations.
1187 :returns: :py:class:`SourceGrid` instance.
1189 Example::
1191 >>> from pyrocko import gf
1192 >>> base = DCSource()
1193 >>> R = gf.Range
1194 >>> for s in base.grid(R('
1196 '''
1197 return SourceGrid(base=self, variables=variables)
1199 def base_key(self):
1200 '''
1201 Get key to decide about source discretization / GF stack sharing.
1203 When two source models differ only in amplitude and origin time, the
1204 discretization and the GF stacking can be done only once for a unit
1205 amplitude and a zero origin time and the amplitude and origin times of
1206 the seismograms can be applied during post-processing of the synthetic
1207 seismogram.
1209 For any derived parameterized source model, this method is called to
1210 decide if discretization and stacking of the source should be shared.
1211 When two source models return an equal vector of values discretization
1212 is shared.
1213 '''
1214 return (self.depth, self.lat, self.north_shift,
1215 self.lon, self.east_shift, self.time, type(self).__name__) + \
1216 self.effective_stf_pre().base_key()
1218 def get_factor(self):
1219 '''
1220 Get the scaling factor to be applied during post-processing.
1222 Discretization of the base seismogram is usually done for a unit
1223 amplitude, because a common factor can be efficiently multiplied to
1224 final seismograms. This eliminates to do repeat the stacking when
1225 creating seismograms for a series of source models only differing in
1226 amplitude.
1228 This method should return the scaling factor to apply in the
1229 post-processing (often this is simply the scalar moment of the source).
1230 '''
1232 return 1.0
1234 def effective_stf_pre(self):
1235 '''
1236 Return the STF applied before stacking of the Green's functions.
1238 This STF is used during discretization of the parameterized source
1239 models, i.e. to produce a temporal distribution of point sources.
1241 Handling of the STF before stacking of the GFs is less efficient but
1242 allows to use different source time functions for different parts of
1243 the source.
1244 '''
1246 if self.stf is not None and self.stf_mode == 'pre':
1247 return self.stf
1248 else:
1249 return g_unit_pulse
1251 def effective_stf_post(self):
1252 '''
1253 Return the STF applied after stacking of the Green's fuctions.
1255 This STF is used in the post-processing of the synthetic seismograms.
1257 Handling of the STF after stacking of the GFs is usually more efficient
1258 but is only possible when a common STF is used for all subsources.
1259 '''
1261 if self.stf is not None and self.stf_mode == 'post':
1262 return self.stf
1263 else:
1264 return g_unit_pulse
1266 def _dparams_base(self):
1267 return dict(times=arr(self.time),
1268 lat=self.lat, lon=self.lon,
1269 north_shifts=arr(self.north_shift),
1270 east_shifts=arr(self.east_shift),
1271 depths=arr(self.depth))
1273 def _hash(self):
1274 sha = sha1()
1275 for k in self.base_key():
1276 sha.update(str(k).encode())
1277 return sha.hexdigest()
1279 def _dparams_base_repeated(self, times):
1280 if times is None:
1281 return self._dparams_base()
1283 nt = times.size
1284 north_shifts = num.repeat(self.north_shift, nt)
1285 east_shifts = num.repeat(self.east_shift, nt)
1286 depths = num.repeat(self.depth, nt)
1287 return dict(times=times,
1288 lat=self.lat, lon=self.lon,
1289 north_shifts=north_shifts,
1290 east_shifts=east_shifts,
1291 depths=depths)
1293 def pyrocko_event(self, store=None, target=None, **kwargs):
1294 duration = None
1295 if self.stf:
1296 duration = self.stf.effective_duration
1298 return model.Event(
1299 lat=self.lat,
1300 lon=self.lon,
1301 north_shift=self.north_shift,
1302 east_shift=self.east_shift,
1303 time=self.time,
1304 name=self.name,
1305 depth=self.depth,
1306 duration=duration,
1307 **kwargs)
1309 def outline(self, cs='xyz'):
1310 points = num.atleast_2d(num.zeros([1, 3]))
1312 points[:, 0] += self.north_shift
1313 points[:, 1] += self.east_shift
1314 points[:, 2] += self.depth
1315 if cs == 'xyz':
1316 return points
1317 elif cs == 'xy':
1318 return points[:, :2]
1319 elif cs in ('latlon', 'lonlat'):
1320 latlon = ne_to_latlon(
1321 self.lat, self.lon, points[:, 0], points[:, 1])
1323 latlon = num.array(latlon).T
1324 if cs == 'latlon':
1325 return latlon
1326 else:
1327 return latlon[:, ::-1]
1329 @classmethod
1330 def from_pyrocko_event(cls, ev, **kwargs):
1331 if ev.depth is None:
1332 raise ConversionError(
1333 'Cannot convert event object to source object: '
1334 'no depth information available')
1336 stf = None
1337 if ev.duration is not None:
1338 stf = HalfSinusoidSTF(effective_duration=ev.duration)
1340 d = dict(
1341 name=ev.name,
1342 time=ev.time,
1343 lat=ev.lat,
1344 lon=ev.lon,
1345 north_shift=ev.north_shift,
1346 east_shift=ev.east_shift,
1347 depth=ev.depth,
1348 stf=stf)
1349 d.update(kwargs)
1350 return cls(**d)
1352 def get_magnitude(self):
1353 raise NotImplementedError(
1354 '%s does not implement get_magnitude()'
1355 % self.__class__.__name__)
1358class SourceWithMagnitude(Source):
1359 '''
1360 Base class for sources containing a moment magnitude.
1361 '''
1363 magnitude = Float.T(
1364 default=6.0,
1365 help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1367 def __init__(self, **kwargs):
1368 if 'moment' in kwargs:
1369 mom = kwargs.pop('moment')
1370 if 'magnitude' not in kwargs:
1371 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1373 Source.__init__(self, **kwargs)
1375 @property
1376 def moment(self):
1377 return float(pmt.magnitude_to_moment(self.magnitude))
1379 @moment.setter
1380 def moment(self, value):
1381 self.magnitude = float(pmt.moment_to_magnitude(value))
1383 def pyrocko_event(self, store=None, target=None, **kwargs):
1384 return Source.pyrocko_event(
1385 self, store, target,
1386 magnitude=self.magnitude,
1387 **kwargs)
1389 @classmethod
1390 def from_pyrocko_event(cls, ev, **kwargs):
1391 d = {}
1392 if ev.magnitude:
1393 d.update(magnitude=ev.magnitude)
1395 d.update(kwargs)
1396 return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
1398 def get_magnitude(self):
1399 return self.magnitude
1402class DerivedMagnitudeError(ValidationError):
1403 pass
1406class SourceWithDerivedMagnitude(Source):
1408 class __T(Source.T):
1410 def validate_extra(self, val):
1411 Source.T.validate_extra(self, val)
1412 val.check_conflicts()
1414 def check_conflicts(self):
1415 '''
1416 Check for parameter conflicts.
1418 To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError`
1419 on conflicts.
1420 '''
1421 pass
1423 def get_magnitude(self, store=None, target=None):
1424 raise DerivedMagnitudeError('No magnitude set.')
1426 def get_moment(self, store=None, target=None):
1427 return float(pmt.magnitude_to_moment(
1428 self.get_magnitude(store, target)))
1430 def pyrocko_moment_tensor(self, store=None, target=None):
1431 raise NotImplementedError(
1432 '%s does not implement pyrocko_moment_tensor()'
1433 % self.__class__.__name__)
1435 def pyrocko_event(self, store=None, target=None, **kwargs):
1436 try:
1437 mt = self.pyrocko_moment_tensor(store, target)
1438 magnitude = self.get_magnitude()
1439 except (DerivedMagnitudeError, NotImplementedError):
1440 mt = None
1441 magnitude = None
1443 return Source.pyrocko_event(
1444 self, store, target,
1445 moment_tensor=mt,
1446 magnitude=magnitude,
1447 **kwargs)
1450class ExplosionSource(SourceWithDerivedMagnitude):
1451 '''
1452 An isotropic explosion point source.
1453 '''
1455 magnitude = Float.T(
1456 optional=True,
1457 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1459 volume_change = Float.T(
1460 optional=True,
1461 help='volume change of the explosion/implosion or '
1462 'the contracting/extending magmatic source. [m^3]')
1464 discretized_source_class = meta.DiscretizedExplosionSource
1466 def __init__(self, **kwargs):
1467 if 'moment' in kwargs:
1468 mom = kwargs.pop('moment')
1469 if 'magnitude' not in kwargs:
1470 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1472 SourceWithDerivedMagnitude.__init__(self, **kwargs)
1474 def base_key(self):
1475 return SourceWithDerivedMagnitude.base_key(self) + \
1476 (self.volume_change,)
1478 def check_conflicts(self):
1479 if self.magnitude is not None and self.volume_change is not None:
1480 raise DerivedMagnitudeError(
1481 'Magnitude and volume_change are both defined.')
1483 def get_magnitude(self, store=None, target=None):
1484 self.check_conflicts()
1486 if self.magnitude is not None:
1487 return self.magnitude
1489 elif self.volume_change is not None:
1490 moment = self.volume_change * \
1491 self.get_moment_to_volume_change_ratio(store, target)
1493 return float(pmt.moment_to_magnitude(abs(moment)))
1494 else:
1495 return float(pmt.moment_to_magnitude(1.0))
1497 def get_volume_change(self, store=None, target=None):
1498 self.check_conflicts()
1500 if self.volume_change is not None:
1501 return self.volume_change
1503 elif self.magnitude is not None:
1504 moment = float(pmt.magnitude_to_moment(self.magnitude))
1505 return moment / self.get_moment_to_volume_change_ratio(
1506 store, target)
1508 else:
1509 return 1.0 / self.get_moment_to_volume_change_ratio(store)
1511 def get_moment_to_volume_change_ratio(self, store, target=None):
1512 if store is None:
1513 raise DerivedMagnitudeError(
1514 'Need earth model to convert between volume change and '
1515 'magnitude.')
1517 points = num.array(
1518 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1520 interpolation = target.interpolation if target else 'multilinear'
1521 try:
1522 shear_moduli = store.config.get_shear_moduli(
1523 self.lat, self.lon,
1524 points=points,
1525 interpolation=interpolation)[0]
1527 bulk_moduli = store.config.get_bulk_moduli(
1528 self.lat, self.lon,
1529 points=points,
1530 interpolation=interpolation)[0]
1531 except meta.OutOfBounds:
1532 raise DerivedMagnitudeError(
1533 'Could not get shear modulus at source position.')
1535 return float(2. * shear_moduli + bulk_moduli)
1537 def get_factor(self):
1538 return 1.0
1540 def discretize_basesource(self, store, target=None):
1541 times, amplitudes = self.effective_stf_pre().discretize_t(
1542 store.config.deltat, self.time)
1544 amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.)
1546 if self.volume_change is not None:
1547 if self.volume_change < 0.:
1548 amplitudes *= -1
1550 return meta.DiscretizedExplosionSource(
1551 m0s=amplitudes,
1552 **self._dparams_base_repeated(times))
1554 def pyrocko_moment_tensor(self, store=None, target=None):
1555 a = self.get_moment(store, target) * math.sqrt(2. / 3.)
1556 return pmt.MomentTensor(m=pmt.symmat6(a, a, a, 0., 0., 0.))
1559class RectangularExplosionSource(ExplosionSource):
1560 '''
1561 Rectangular or line explosion source.
1562 '''
1564 discretized_source_class = meta.DiscretizedExplosionSource
1566 strike = Float.T(
1567 default=0.0,
1568 help='strike direction in [deg], measured clockwise from north')
1570 dip = Float.T(
1571 default=90.0,
1572 help='dip angle in [deg], measured downward from horizontal')
1574 length = Float.T(
1575 default=0.,
1576 help='length of rectangular source area [m]')
1578 width = Float.T(
1579 default=0.,
1580 help='width of rectangular source area [m]')
1582 anchor = StringChoice.T(
1583 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
1584 'bottom_left', 'bottom_right'],
1585 default='center',
1586 optional=True,
1587 help='Anchor point for positioning the plane, can be: top, center or'
1588 'bottom and also top_left, top_right,bottom_left,'
1589 'bottom_right, center_left and center right')
1591 nucleation_x = Float.T(
1592 optional=True,
1593 help='horizontal position of rupture nucleation in normalized fault '
1594 'plane coordinates (-1 = left edge, +1 = right edge)')
1596 nucleation_y = Float.T(
1597 optional=True,
1598 help='down-dip position of rupture nucleation in normalized fault '
1599 'plane coordinates (-1 = upper edge, +1 = lower edge)')
1601 velocity = Float.T(
1602 default=3500.,
1603 help='speed of explosion front [m/s]')
1605 aggressive_oversampling = Bool.T(
1606 default=False,
1607 help='Aggressive oversampling for basesource discretization. '
1608 'When using \'multilinear\' interpolation oversampling has'
1609 ' practically no effect.')
1611 def base_key(self):
1612 return Source.base_key(self) + (self.strike, self.dip, self.length,
1613 self.width, self.nucleation_x,
1614 self.nucleation_y, self.velocity,
1615 self.anchor)
1617 def discretize_basesource(self, store, target=None):
1619 if self.nucleation_x is not None:
1620 nucx = self.nucleation_x * 0.5 * self.length
1621 else:
1622 nucx = None
1624 if self.nucleation_y is not None:
1625 nucy = self.nucleation_y * 0.5 * self.width
1626 else:
1627 nucy = None
1629 stf = self.effective_stf_pre()
1631 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
1632 store.config.deltas, store.config.deltat,
1633 self.time, self.north_shift, self.east_shift, self.depth,
1634 self.strike, self.dip, self.length, self.width, self.anchor,
1635 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy)
1637 amplitudes /= num.sum(amplitudes)
1638 amplitudes *= self.get_moment(store, target)
1640 return meta.DiscretizedExplosionSource(
1641 lat=self.lat,
1642 lon=self.lon,
1643 times=times,
1644 north_shifts=points[:, 0],
1645 east_shifts=points[:, 1],
1646 depths=points[:, 2],
1647 m0s=amplitudes)
1649 def outline(self, cs='xyz'):
1650 points = outline_rect_source(self.strike, self.dip, self.length,
1651 self.width, self.anchor)
1653 points[:, 0] += self.north_shift
1654 points[:, 1] += self.east_shift
1655 points[:, 2] += self.depth
1656 if cs == 'xyz':
1657 return points
1658 elif cs == 'xy':
1659 return points[:, :2]
1660 elif cs in ('latlon', 'lonlat'):
1661 latlon = ne_to_latlon(
1662 self.lat, self.lon, points[:, 0], points[:, 1])
1664 latlon = num.array(latlon).T
1665 if cs == 'latlon':
1666 return latlon
1667 else:
1668 return latlon[:, ::-1]
1670 def get_nucleation_abs_coord(self, cs='xy'):
1672 if self.nucleation_x is None:
1673 return None, None
1675 coords = from_plane_coords(self.strike, self.dip, self.length,
1676 self.width, self.depth, self.nucleation_x,
1677 self.nucleation_y, lat=self.lat,
1678 lon=self.lon, north_shift=self.north_shift,
1679 east_shift=self.east_shift, cs=cs)
1680 return coords
1683class DCSource(SourceWithMagnitude):
1684 '''
1685 A double-couple point source.
1686 '''
1688 strike = Float.T(
1689 default=0.0,
1690 help='strike direction in [deg], measured clockwise from north')
1692 dip = Float.T(
1693 default=90.0,
1694 help='dip angle in [deg], measured downward from horizontal')
1696 rake = Float.T(
1697 default=0.0,
1698 help='rake angle in [deg], '
1699 'measured counter-clockwise from right-horizontal '
1700 'in on-plane view')
1702 discretized_source_class = meta.DiscretizedMTSource
1704 def base_key(self):
1705 return Source.base_key(self) + (self.strike, self.dip, self.rake)
1707 def get_factor(self):
1708 return float(pmt.magnitude_to_moment(self.magnitude))
1710 def discretize_basesource(self, store, target=None):
1711 mot = pmt.MomentTensor(
1712 strike=self.strike, dip=self.dip, rake=self.rake)
1714 times, amplitudes = self.effective_stf_pre().discretize_t(
1715 store.config.deltat, self.time)
1716 return meta.DiscretizedMTSource(
1717 m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis],
1718 **self._dparams_base_repeated(times))
1720 def pyrocko_moment_tensor(self, store=None, target=None):
1721 return pmt.MomentTensor(
1722 strike=self.strike,
1723 dip=self.dip,
1724 rake=self.rake,
1725 scalar_moment=self.moment)
1727 def pyrocko_event(self, store=None, target=None, **kwargs):
1728 return SourceWithMagnitude.pyrocko_event(
1729 self, store, target,
1730 moment_tensor=self.pyrocko_moment_tensor(store, target),
1731 **kwargs)
1733 @classmethod
1734 def from_pyrocko_event(cls, ev, **kwargs):
1735 d = {}
1736 mt = ev.moment_tensor
1737 if mt:
1738 (strike, dip, rake), _ = mt.both_strike_dip_rake()
1739 d.update(
1740 strike=float(strike),
1741 dip=float(dip),
1742 rake=float(rake),
1743 magnitude=float(mt.moment_magnitude()))
1745 d.update(kwargs)
1746 return super(DCSource, cls).from_pyrocko_event(ev, **d)
1749class CLVDSource(SourceWithMagnitude):
1750 '''
1751 A pure CLVD point source.
1752 '''
1754 discretized_source_class = meta.DiscretizedMTSource
1756 azimuth = Float.T(
1757 default=0.0,
1758 help='azimuth direction of largest dipole, clockwise from north [deg]')
1760 dip = Float.T(
1761 default=90.,
1762 help='dip direction of largest dipole, downward from horizontal [deg]')
1764 def base_key(self):
1765 return Source.base_key(self) + (self.azimuth, self.dip)
1767 def get_factor(self):
1768 return float(pmt.magnitude_to_moment(self.magnitude))
1770 @property
1771 def m6(self):
1772 a = math.sqrt(4. / 3.) * self.get_factor()
1773 m = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1774 rotmat1 = pmt.euler_to_matrix(
1775 d2r * (self.dip - 90.),
1776 d2r * (self.azimuth - 90.),
1777 0.)
1778 m = num.dot(rotmat1.T, num.dot(m, rotmat1))
1779 return pmt.to6(m)
1781 @property
1782 def m6_astuple(self):
1783 return tuple(self.m6.tolist())
1785 def discretize_basesource(self, store, target=None):
1786 factor = self.get_factor()
1787 times, amplitudes = self.effective_stf_pre().discretize_t(
1788 store.config.deltat, self.time)
1789 return meta.DiscretizedMTSource(
1790 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor,
1791 **self._dparams_base_repeated(times))
1793 def pyrocko_moment_tensor(self, store=None, target=None):
1794 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
1796 def pyrocko_event(self, store=None, target=None, **kwargs):
1797 mt = self.pyrocko_moment_tensor(store, target)
1798 return Source.pyrocko_event(
1799 self, store, target,
1800 moment_tensor=self.pyrocko_moment_tensor(store, target),
1801 magnitude=float(mt.moment_magnitude()),
1802 **kwargs)
1805class VLVDSource(SourceWithDerivedMagnitude):
1806 '''
1807 Volumetric linear vector dipole source.
1809 This source is a parameterization for a restricted moment tensor point
1810 source, useful to represent dyke or sill like inflation or deflation
1811 sources. The restriction is such that the moment tensor is rotational
1812 symmetric. It can be represented by a superposition of a linear vector
1813 dipole (here we use a CLVD for convenience) and an isotropic component. The
1814 restricted moment tensor has 4 degrees of freedom: 2 independent
1815 eigenvalues and 2 rotation angles orienting the the symmetry axis.
1817 In this parameterization, the isotropic component is controlled by
1818 ``volume_change``. To define the moment tensor, it must be converted to the
1819 scalar moment of the the MT's isotropic component. For the conversion, the
1820 shear modulus at the source's position must be known. This value is
1821 extracted from the earth model defined in the GF store in use.
1823 The CLVD part by controlled by its scalar moment :math:`M_0`:
1824 ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a
1825 positiv or negativ CLVD (the sign of the largest eigenvalue).
1826 '''
1828 discretized_source_class = meta.DiscretizedMTSource
1830 azimuth = Float.T(
1831 default=0.0,
1832 help='azimuth direction of symmetry axis, clockwise from north [deg].')
1834 dip = Float.T(
1835 default=90.,
1836 help='dip direction of symmetry axis, downward from horizontal [deg].')
1838 volume_change = Float.T(
1839 default=0.,
1840 help='volume change of the inflation/deflation [m^3].')
1842 clvd_moment = Float.T(
1843 default=0.,
1844 help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign '
1845 'controls the sign of the CLVD (the sign of its largest '
1846 'eigenvalue).')
1848 def get_moment_to_volume_change_ratio(self, store, target):
1849 if store is None or target is None:
1850 raise DerivedMagnitudeError(
1851 'Need earth model to convert between volume change and '
1852 'magnitude.')
1854 points = num.array(
1855 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1857 try:
1858 shear_moduli = store.config.get_shear_moduli(
1859 self.lat, self.lon,
1860 points=points,
1861 interpolation=target.interpolation)[0]
1863 bulk_moduli = store.config.get_bulk_moduli(
1864 self.lat, self.lon,
1865 points=points,
1866 interpolation=target.interpolation)[0]
1867 except meta.OutOfBounds:
1868 raise DerivedMagnitudeError(
1869 'Could not get shear modulus at source position.')
1871 return float(2. * shear_moduli + bulk_moduli)
1873 def base_key(self):
1874 return Source.base_key(self) + \
1875 (self.azimuth, self.dip, self.volume_change, self.clvd_moment)
1877 def get_magnitude(self, store=None, target=None):
1878 mt = self.pyrocko_moment_tensor(store, target)
1879 return float(pmt.moment_to_magnitude(mt.moment))
1881 def get_m6(self, store, target):
1882 a = math.sqrt(4. / 3.) * self.clvd_moment
1883 m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1885 rotmat1 = pmt.euler_to_matrix(
1886 d2r * (self.dip - 90.),
1887 d2r * (self.azimuth - 90.),
1888 0.)
1889 m_clvd = num.dot(rotmat1.T, num.dot(m_clvd, rotmat1))
1891 m_iso = self.volume_change * \
1892 self.get_moment_to_volume_change_ratio(store, target)
1894 m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0.,
1895 0., 0.,) * math.sqrt(2. / 3)
1897 m = pmt.to6(m_clvd) + pmt.to6(m_iso)
1898 return m
1900 def get_moment(self, store=None, target=None):
1901 return float(pmt.magnitude_to_moment(
1902 self.get_magnitude(store, target)))
1904 def get_m6_astuple(self, store, target):
1905 m6 = self.get_m6(store, target)
1906 return tuple(m6.tolist())
1908 def discretize_basesource(self, store, target=None):
1909 times, amplitudes = self.effective_stf_pre().discretize_t(
1910 store.config.deltat, self.time)
1912 m6 = self.get_m6(store, target)
1913 m6 *= amplitudes / self.get_factor()
1915 return meta.DiscretizedMTSource(
1916 m6s=m6[num.newaxis, :],
1917 **self._dparams_base_repeated(times))
1919 def pyrocko_moment_tensor(self, store=None, target=None):
1920 m6_astuple = self.get_m6_astuple(store, target)
1921 return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple))
1924class MTSource(Source):
1925 '''
1926 A moment tensor point source.
1927 '''
1929 discretized_source_class = meta.DiscretizedMTSource
1931 mnn = Float.T(
1932 default=1.,
1933 help='north-north component of moment tensor in [Nm]')
1935 mee = Float.T(
1936 default=1.,
1937 help='east-east component of moment tensor in [Nm]')
1939 mdd = Float.T(
1940 default=1.,
1941 help='down-down component of moment tensor in [Nm]')
1943 mne = Float.T(
1944 default=0.,
1945 help='north-east component of moment tensor in [Nm]')
1947 mnd = Float.T(
1948 default=0.,
1949 help='north-down component of moment tensor in [Nm]')
1951 med = Float.T(
1952 default=0.,
1953 help='east-down component of moment tensor in [Nm]')
1955 def __init__(self, **kwargs):
1956 if 'm6' in kwargs:
1957 for (k, v) in zip('mnn mee mdd mne mnd med'.split(),
1958 kwargs.pop('m6')):
1959 kwargs[k] = float(v)
1961 Source.__init__(self, **kwargs)
1963 @property
1964 def m6(self):
1965 return num.array(self.m6_astuple)
1967 @property
1968 def m6_astuple(self):
1969 return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)
1971 @m6.setter
1972 def m6(self, value):
1973 self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value
1975 def base_key(self):
1976 return Source.base_key(self) + self.m6_astuple
1978 def discretize_basesource(self, store, target=None):
1979 times, amplitudes = self.effective_stf_pre().discretize_t(
1980 store.config.deltat, self.time)
1981 return meta.DiscretizedMTSource(
1982 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis],
1983 **self._dparams_base_repeated(times))
1985 def get_magnitude(self, store=None, target=None):
1986 m6 = self.m6
1987 return pmt.moment_to_magnitude(
1988 math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) /
1989 math.sqrt(2.))
1991 def pyrocko_moment_tensor(self, store=None, target=None):
1992 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
1994 def pyrocko_event(self, store=None, target=None, **kwargs):
1995 mt = self.pyrocko_moment_tensor(store, target)
1996 return Source.pyrocko_event(
1997 self, store, target,
1998 moment_tensor=self.pyrocko_moment_tensor(store, target),
1999 magnitude=float(mt.moment_magnitude()),
2000 **kwargs)
2002 @classmethod
2003 def from_pyrocko_event(cls, ev, **kwargs):
2004 d = {}
2005 mt = ev.moment_tensor
2006 if mt:
2007 d.update(m6=tuple(map(float, mt.m6())))
2008 else:
2009 if ev.magnitude is not None:
2010 mom = pmt.magnitude_to_moment(ev.magnitude)
2011 v = math.sqrt(2. / 3.) * mom
2012 d.update(m6=(v, v, v, 0., 0., 0.))
2014 d.update(kwargs)
2015 return super(MTSource, cls).from_pyrocko_event(ev, **d)
2018map_anchor = {
2019 'center': (0.0, 0.0),
2020 'center_left': (-1.0, 0.0),
2021 'center_right': (1.0, 0.0),
2022 'top': (0.0, -1.0),
2023 'top_left': (-1.0, -1.0),
2024 'top_right': (1.0, -1.0),
2025 'bottom': (0.0, 1.0),
2026 'bottom_left': (-1.0, 1.0),
2027 'bottom_right': (1.0, 1.0)}
2030class RectangularSource(SourceWithDerivedMagnitude):
2031 '''
2032 Classical Haskell source model modified for bilateral rupture.
2033 '''
2035 discretized_source_class = meta.DiscretizedMTSource
2037 magnitude = Float.T(
2038 optional=True,
2039 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
2041 strike = Float.T(
2042 default=0.0,
2043 help='strike direction in [deg], measured clockwise from north')
2045 dip = Float.T(
2046 default=90.0,
2047 help='dip angle in [deg], measured downward from horizontal')
2049 rake = Float.T(
2050 default=0.0,
2051 help='rake angle in [deg], '
2052 'measured counter-clockwise from right-horizontal '
2053 'in on-plane view')
2055 length = Float.T(
2056 default=0.,
2057 help='length of rectangular source area [m]')
2059 width = Float.T(
2060 default=0.,
2061 help='width of rectangular source area [m]')
2063 anchor = StringChoice.T(
2064 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2065 'bottom_left', 'bottom_right'],
2066 default='center',
2067 optional=True,
2068 help='Anchor point for positioning the plane, can be: ``top, center '
2069 'bottom, top_left, top_right,bottom_left,'
2070 'bottom_right, center_left, center right``.')
2072 nucleation_x = Float.T(
2073 optional=True,
2074 help='horizontal position of rupture nucleation in normalized fault '
2075 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge)')
2077 nucleation_y = Float.T(
2078 optional=True,
2079 help='down-dip position of rupture nucleation in normalized fault '
2080 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge)')
2082 velocity = Float.T(
2083 default=3500.,
2084 help='speed of rupture front [m/s]')
2086 slip = Float.T(
2087 optional=True,
2088 help='Slip on the rectangular source area [m]')
2090 opening_fraction = Float.T(
2091 default=0.,
2092 help='Determines fraction of slip related to opening. '
2093 '(``-1``: pure tensile closing, '
2094 '``0``: pure shear, '
2095 '``1``: pure tensile opening)')
2097 decimation_factor = Int.T(
2098 optional=True,
2099 default=1,
2100 help='Sub-source decimation factor, a larger decimation will'
2101 ' make the result inaccurate but shorten the necessary'
2102 ' computation time (use for testing puposes only).')
2104 aggressive_oversampling = Bool.T(
2105 default=False,
2106 help='Aggressive oversampling for basesource discretization. '
2107 'When using \'multilinear\' interpolation oversampling has'
2108 ' practically no effect.')
2110 def __init__(self, **kwargs):
2111 if 'moment' in kwargs:
2112 mom = kwargs.pop('moment')
2113 if 'magnitude' not in kwargs:
2114 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
2116 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2118 def base_key(self):
2119 return SourceWithDerivedMagnitude.base_key(self) + (
2120 self.magnitude,
2121 self.slip,
2122 self.strike,
2123 self.dip,
2124 self.rake,
2125 self.length,
2126 self.width,
2127 self.nucleation_x,
2128 self.nucleation_y,
2129 self.velocity,
2130 self.decimation_factor,
2131 self.anchor)
2133 def check_conflicts(self):
2134 if self.magnitude is not None and self.slip is not None:
2135 raise DerivedMagnitudeError(
2136 'Magnitude and slip are both defined.')
2138 def get_magnitude(self, store=None, target=None):
2139 self.check_conflicts()
2140 if self.magnitude is not None:
2141 return self.magnitude
2143 elif self.slip is not None:
2144 if None in (store, target):
2145 raise DerivedMagnitudeError(
2146 'Magnitude for a rectangular source with slip defined '
2147 'can only be derived when earth model and target '
2148 'interpolation method are available.')
2150 amplitudes = self._discretize(store, target)[2]
2151 if amplitudes.ndim == 2:
2152 # CLVD component has no net moment, leave out
2153 return float(pmt.moment_to_magnitude(
2154 num.sum(num.abs(amplitudes[0:2, :]).sum())))
2155 else:
2156 return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
2158 else:
2159 return float(pmt.moment_to_magnitude(1.0))
2161 def get_factor(self):
2162 return 1.0
2164 def get_slip_tensile(self):
2165 return self.slip * self.opening_fraction
2167 def get_slip_shear(self):
2168 return self.slip - abs(self.get_slip_tensile)
2170 def _discretize(self, store, target):
2171 if self.nucleation_x is not None:
2172 nucx = self.nucleation_x * 0.5 * self.length
2173 else:
2174 nucx = None
2176 if self.nucleation_y is not None:
2177 nucy = self.nucleation_y * 0.5 * self.width
2178 else:
2179 nucy = None
2181 stf = self.effective_stf_pre()
2183 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
2184 store.config.deltas, store.config.deltat,
2185 self.time, self.north_shift, self.east_shift, self.depth,
2186 self.strike, self.dip, self.length, self.width, self.anchor,
2187 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy,
2188 decimation_factor=self.decimation_factor,
2189 aggressive_oversampling=self.aggressive_oversampling)
2191 if self.slip is not None:
2192 if target is not None:
2193 interpolation = target.interpolation
2194 else:
2195 interpolation = 'nearest_neighbor'
2196 logger.warning(
2197 'no target information available, will use '
2198 '"nearest_neighbor" interpolation when extracting shear '
2199 'modulus from earth model')
2201 shear_moduli = store.config.get_shear_moduli(
2202 self.lat, self.lon,
2203 points=points,
2204 interpolation=interpolation)
2206 tensile_slip = self.get_slip_tensile()
2207 shear_slip = self.slip - abs(tensile_slip)
2209 amplitudes_total = [shear_moduli * shear_slip]
2210 if tensile_slip != 0:
2211 bulk_moduli = store.config.get_bulk_moduli(
2212 self.lat, self.lon,
2213 points=points,
2214 interpolation=interpolation)
2216 tensile_iso = bulk_moduli * tensile_slip
2217 tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip
2219 amplitudes_total.extend([tensile_iso, tensile_clvd])
2221 amplitudes_total = num.vstack(amplitudes_total).squeeze() * \
2222 amplitudes * dl * dw
2224 else:
2225 # normalization to retain total moment
2226 amplitudes_norm = amplitudes / num.sum(amplitudes)
2227 moment = self.get_moment(store, target)
2229 amplitudes_total = [
2230 amplitudes_norm * moment * (1 - abs(self.opening_fraction))]
2231 if self.opening_fraction != 0.:
2232 amplitudes_total.append(
2233 amplitudes_norm * self.opening_fraction * moment)
2235 amplitudes_total = num.vstack(amplitudes_total).squeeze()
2237 return points, times, num.atleast_1d(amplitudes_total), dl, dw, nl, nw
2239 def discretize_basesource(self, store, target=None):
2241 points, times, amplitudes, dl, dw, nl, nw = self._discretize(
2242 store, target)
2244 mot = pmt.MomentTensor(
2245 strike=self.strike, dip=self.dip, rake=self.rake)
2246 m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0)
2248 if amplitudes.ndim == 1:
2249 m6s[:, :] *= amplitudes[:, num.newaxis]
2250 elif amplitudes.ndim == 2:
2251 # shear MT components
2252 rotmat1 = pmt.euler_to_matrix(
2253 d2r * self.dip, d2r * self.strike, d2r * -self.rake)
2254 m6s[:, :] *= amplitudes[0, :][:, num.newaxis]
2256 if amplitudes.shape[0] == 2:
2257 # tensile MT components - moment/magnitude input
2258 tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.)
2259 rot_tensile = pmt.to6(
2260 num.dot(rotmat1.T, num.dot(tensile, rotmat1)))
2262 m6s_tensile = rot_tensile[
2263 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2264 m6s += m6s_tensile
2266 elif amplitudes.shape[0] == 3:
2267 # tensile MT components - slip input
2268 iso = pmt.symmat6(1., 1., 1., 0., 0., 0.)
2269 clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.)
2271 rot_iso = pmt.to6(
2272 num.dot(rotmat1.T, num.dot(iso, rotmat1)))
2273 rot_clvd = pmt.to6(
2274 num.dot(rotmat1.T, num.dot(clvd, rotmat1)))
2276 m6s_iso = rot_iso[
2277 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2278 m6s_clvd = rot_clvd[
2279 num.newaxis, :] * amplitudes[2, :][:, num.newaxis]
2280 m6s += m6s_iso + m6s_clvd
2281 else:
2282 raise ValueError('Unknwown amplitudes shape!')
2283 else:
2284 raise ValueError(
2285 'Unexpected dimension of {}'.format(amplitudes.ndim))
2287 ds = meta.DiscretizedMTSource(
2288 lat=self.lat,
2289 lon=self.lon,
2290 times=times,
2291 north_shifts=points[:, 0],
2292 east_shifts=points[:, 1],
2293 depths=points[:, 2],
2294 m6s=m6s,
2295 dl=dl,
2296 dw=dw,
2297 nl=nl,
2298 nw=nw)
2300 return ds
2302 def outline(self, cs='xyz'):
2303 points = outline_rect_source(self.strike, self.dip, self.length,
2304 self.width, self.anchor)
2306 points[:, 0] += self.north_shift
2307 points[:, 1] += self.east_shift
2308 points[:, 2] += self.depth
2309 if cs == 'xyz':
2310 return points
2311 elif cs == 'xy':
2312 return points[:, :2]
2313 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2314 latlon = ne_to_latlon(
2315 self.lat, self.lon, points[:, 0], points[:, 1])
2317 latlon = num.array(latlon).T
2318 if cs == 'latlon':
2319 return latlon
2320 elif cs == 'lonlat':
2321 return latlon[:, ::-1]
2322 else:
2323 return num.concatenate(
2324 (latlon, points[:, 2].reshape((len(points), 1))),
2325 axis=1)
2327 def points_on_source(self, cs='xyz', **kwargs):
2329 points = points_on_rect_source(
2330 self.strike, self.dip, self.length, self.width,
2331 self.anchor, **kwargs)
2333 points[:, 0] += self.north_shift
2334 points[:, 1] += self.east_shift
2335 points[:, 2] += self.depth
2336 if cs == 'xyz':
2337 return points
2338 elif cs == 'xy':
2339 return points[:, :2]
2340 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2341 latlon = ne_to_latlon(
2342 self.lat, self.lon, points[:, 0], points[:, 1])
2344 latlon = num.array(latlon).T
2345 if cs == 'latlon':
2346 return latlon
2347 elif cs == 'lonlat':
2348 return latlon[:, ::-1]
2349 else:
2350 return num.concatenate(
2351 (latlon, points[:, 2].reshape((len(points), 1))),
2352 axis=1)
2354 def get_nucleation_abs_coord(self, cs='xy'):
2356 if self.nucleation_x is None:
2357 return None, None
2359 coords = from_plane_coords(self.strike, self.dip, self.length,
2360 self.width, self.depth, self.nucleation_x,
2361 self.nucleation_y, lat=self.lat,
2362 lon=self.lon, north_shift=self.north_shift,
2363 east_shift=self.east_shift, cs=cs)
2364 return coords
2366 def pyrocko_moment_tensor(self, store=None, target=None):
2367 return pmt.MomentTensor(
2368 strike=self.strike,
2369 dip=self.dip,
2370 rake=self.rake,
2371 scalar_moment=self.get_moment(store, target))
2373 def pyrocko_event(self, store=None, target=None, **kwargs):
2374 return SourceWithDerivedMagnitude.pyrocko_event(
2375 self, store, target,
2376 **kwargs)
2378 @classmethod
2379 def from_pyrocko_event(cls, ev, **kwargs):
2380 d = {}
2381 mt = ev.moment_tensor
2382 if mt:
2383 (strike, dip, rake), _ = mt.both_strike_dip_rake()
2384 d.update(
2385 strike=float(strike),
2386 dip=float(dip),
2387 rake=float(rake),
2388 magnitude=float(mt.moment_magnitude()))
2390 d.update(kwargs)
2391 return super(RectangularSource, cls).from_pyrocko_event(ev, **d)
2394class PseudoDynamicRupture(SourceWithDerivedMagnitude):
2395 '''
2396 Combined Eikonal and Okada quasi-dynamic rupture model.
2398 Details are described in :doc:`/topics/pseudo-dynamic-rupture`.
2399 Note: attribute `stf` is not used so far, but kept for future applications.
2400 '''
2402 discretized_source_class = meta.DiscretizedMTSource
2404 strike = Float.T(
2405 default=0.0,
2406 help='Strike direction in [deg], measured clockwise from north.')
2408 dip = Float.T(
2409 default=0.0,
2410 help='Dip angle in [deg], measured downward from horizontal.')
2412 length = Float.T(
2413 default=10. * km,
2414 help='Length of rectangular source area in [m].')
2416 width = Float.T(
2417 default=5. * km,
2418 help='Width of rectangular source area in [m].')
2420 anchor = StringChoice.T(
2421 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2422 'bottom_left', 'bottom_right'],
2423 default='center',
2424 optional=True,
2425 help='Anchor point for positioning the plane, can be: ``top, center, '
2426 'bottom, top_left, top_right, bottom_left, '
2427 'bottom_right, center_left, center_right``.')
2429 nucleation_x__ = Array.T(
2430 default=num.array([0.]),
2431 dtype=num.float64,
2432 serialize_as='list',
2433 help='Horizontal position of rupture nucleation in normalized fault '
2434 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge).')
2436 nucleation_y__ = Array.T(
2437 default=num.array([0.]),
2438 dtype=num.float64,
2439 serialize_as='list',
2440 help='Down-dip position of rupture nucleation in normalized fault '
2441 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge).')
2443 nucleation_time__ = Array.T(
2444 optional=True,
2445 help='Time in [s] after origin, when nucleation points defined by '
2446 '``nucleation_x`` and ``nucleation_y`` rupture.',
2447 dtype=num.float64,
2448 serialize_as='list')
2450 gamma = Float.T(
2451 default=0.8,
2452 help='Scaling factor between rupture velocity and S-wave velocity: '
2453 r':math:`v_r = \gamma * v_s`.')
2455 nx = Int.T(
2456 default=2,
2457 help='Number of discrete source patches in x direction (along '
2458 'strike).')
2460 ny = Int.T(
2461 default=2,
2462 help='Number of discrete source patches in y direction (down dip).')
2464 slip = Float.T(
2465 optional=True,
2466 help='Maximum slip of the rectangular source [m]. '
2467 'Setting the slip the tractions/stress field '
2468 'will be normalized to accomodate the desired maximum slip.')
2470 rake = Float.T(
2471 optional=True,
2472 help='Rake angle in [deg], '
2473 'measured counter-clockwise from right-horizontal '
2474 'in on-plane view. Rake is translated into homogenous tractions '
2475 'in strike and up-dip direction. ``rake`` is mutually exclusive '
2476 'with tractions parameter.')
2478 patches = List.T(
2479 OkadaSource.T(),
2480 optional=True,
2481 help='List of all boundary elements/sub faults/fault patches.')
2483 patch_mask__ = Array.T(
2484 dtype=bool,
2485 serialize_as='list',
2486 shape=(None,),
2487 optional=True,
2488 help='Mask for all boundary elements/sub faults/fault patches. True '
2489 'leaves the patch in the calculation, False excludes the patch.')
2491 tractions = TractionField.T(
2492 optional=True,
2493 help='Traction field the rupture plane is exposed to. See the'
2494 ':py:mod:`pyrocko.gf.tractions` module for more details. '
2495 'If ``tractions=None`` and ``rake`` is given'
2496 ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will'
2497 ' be used.')
2499 coef_mat = Array.T(
2500 optional=True,
2501 help='Coefficient matrix linking traction and dislocation field.',
2502 dtype=num.float64,
2503 shape=(None, None))
2505 eikonal_decimation = Int.T(
2506 optional=True,
2507 default=1,
2508 help='Sub-source eikonal factor, a smaller eikonal factor will'
2509 ' increase the accuracy of rupture front calculation but'
2510 ' increases also the computation time.')
2512 decimation_factor = Int.T(
2513 optional=True,
2514 default=1,
2515 help='Sub-source decimation factor, a larger decimation will'
2516 ' make the result inaccurate but shorten the necessary'
2517 ' computation time (use for testing puposes only).')
2519 nthreads = Int.T(
2520 optional=True,
2521 default=1,
2522 help='Number of threads for Okada forward modelling, '
2523 'matrix inversion and calculation of point subsources. '
2524 'Note: for small/medium matrices 1 thread is most efficient.')
2526 pure_shear = Bool.T(
2527 optional=True,
2528 default=False,
2529 help='Calculate only shear tractions and omit tensile tractions.')
2531 smooth_rupture = Bool.T(
2532 default=True,
2533 help='Smooth the tractions by weighting partially ruptured'
2534 ' fault patches.')
2536 aggressive_oversampling = Bool.T(
2537 default=False,
2538 help='Aggressive oversampling for basesource discretization. '
2539 'When using \'multilinear\' interpolation oversampling has'
2540 ' practically no effect.')
2542 def __init__(self, **kwargs):
2543 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2544 self._interpolators = {}
2545 self.check_conflicts()
2547 @property
2548 def nucleation_x(self):
2549 return self.nucleation_x__
2551 @nucleation_x.setter
2552 def nucleation_x(self, nucleation_x):
2553 if isinstance(nucleation_x, list):
2554 nucleation_x = num.array(nucleation_x)
2556 elif not isinstance(
2557 nucleation_x, num.ndarray) and nucleation_x is not None:
2559 nucleation_x = num.array([nucleation_x])
2560 self.nucleation_x__ = nucleation_x
2562 @property
2563 def nucleation_y(self):
2564 return self.nucleation_y__
2566 @nucleation_y.setter
2567 def nucleation_y(self, nucleation_y):
2568 if isinstance(nucleation_y, list):
2569 nucleation_y = num.array(nucleation_y)
2571 elif not isinstance(nucleation_y, num.ndarray) \
2572 and nucleation_y is not None:
2573 nucleation_y = num.array([nucleation_y])
2575 self.nucleation_y__ = nucleation_y
2577 @property
2578 def nucleation(self):
2579 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
2581 if (nucl_x is None) or (nucl_y is None):
2582 return None
2584 assert nucl_x.shape[0] == nucl_y.shape[0]
2586 return num.concatenate(
2587 (nucl_x[:, num.newaxis], nucl_y[:, num.newaxis]), axis=1)
2589 @nucleation.setter
2590 def nucleation(self, nucleation):
2591 if isinstance(nucleation, list):
2592 nucleation = num.array(nucleation)
2594 assert nucleation.shape[1] == 2
2596 self.nucleation_x = nucleation[:, 0]
2597 self.nucleation_y = nucleation[:, 1]
2599 @property
2600 def nucleation_time(self):
2601 return self.nucleation_time__
2603 @nucleation_time.setter
2604 def nucleation_time(self, nucleation_time):
2605 if not isinstance(nucleation_time, num.ndarray) \
2606 and nucleation_time is not None:
2607 nucleation_time = num.array([nucleation_time])
2609 self.nucleation_time__ = nucleation_time
2611 @property
2612 def patch_mask(self):
2613 if (self.patch_mask__ is not None and
2614 self.patch_mask__.shape == (self.nx * self.ny,)):
2616 return self.patch_mask__
2617 else:
2618 return num.ones(self.nx * self.ny, dtype=bool)
2620 @patch_mask.setter
2621 def patch_mask(self, patch_mask):
2622 if isinstance(patch_mask, list):
2623 patch_mask = num.array(patch_mask)
2625 self.patch_mask__ = patch_mask
2627 def get_tractions(self):
2628 '''
2629 Get source traction vectors.
2631 If :py:attr:`rake` is given, unit length directed traction vectors
2632 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are returned,
2633 else the given :py:attr:`tractions` are used.
2635 :returns:
2636 Traction vectors per patch.
2637 :rtype:
2638 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2639 '''
2641 if self.rake is not None:
2642 if num.isnan(self.rake):
2643 raise ValueError('Rake must be a real number, not NaN.')
2645 logger.warning(
2646 'Tractions are derived based on the given source rake.')
2647 tractions = DirectedTractions(rake=self.rake)
2648 else:
2649 tractions = self.tractions
2650 return tractions.get_tractions(self.nx, self.ny, self.patches)
2652 def base_key(self):
2653 return SourceWithDerivedMagnitude.base_key(self) + (
2654 self.slip,
2655 self.strike,
2656 self.dip,
2657 self.rake,
2658 self.length,
2659 self.width,
2660 float(self.nucleation_x.mean()),
2661 float(self.nucleation_y.mean()),
2662 self.decimation_factor,
2663 self.anchor,
2664 self.pure_shear,
2665 self.gamma,
2666 tuple(self.patch_mask))
2668 def check_conflicts(self):
2669 if self.tractions and self.rake:
2670 raise AttributeError(
2671 'Tractions and rake are mutually exclusive.')
2672 if self.tractions is None and self.rake is None:
2673 self.rake = 0.
2675 def get_magnitude(self, store=None, target=None):
2676 self.check_conflicts()
2677 if self.slip is not None or self.tractions is not None:
2678 if store is None:
2679 raise DerivedMagnitudeError(
2680 'Magnitude for a rectangular source with slip or '
2681 'tractions defined can only be derived when earth model '
2682 'is set.')
2684 moment_rate, calc_times = self.discretize_basesource(
2685 store, target=target).get_moment_rate(store.config.deltat)
2687 deltat = num.concatenate((
2688 (num.diff(calc_times)[0],),
2689 num.diff(calc_times)))
2691 return float(pmt.moment_to_magnitude(
2692 num.sum(moment_rate * deltat)))
2694 else:
2695 return float(pmt.moment_to_magnitude(1.0))
2697 def get_factor(self):
2698 return 1.0
2700 def outline(self, cs='xyz'):
2701 '''
2702 Get source outline corner coordinates.
2704 :param cs:
2705 :ref:`Output coordinate system <coordinate-system-names>`.
2706 :type cs:
2707 optional, str
2709 :returns:
2710 Corner points in desired coordinate system.
2711 :rtype:
2712 :py:class:`~numpy.ndarray`: ``(5, [2, 3])``.
2713 '''
2714 points = outline_rect_source(self.strike, self.dip, self.length,
2715 self.width, self.anchor)
2717 points[:, 0] += self.north_shift
2718 points[:, 1] += self.east_shift
2719 points[:, 2] += self.depth
2720 if cs == 'xyz':
2721 return points
2722 elif cs == 'xy':
2723 return points[:, :2]
2724 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2725 latlon = ne_to_latlon(
2726 self.lat, self.lon, points[:, 0], points[:, 1])
2728 latlon = num.array(latlon).T
2729 if cs == 'latlon':
2730 return latlon
2731 elif cs == 'lonlat':
2732 return latlon[:, ::-1]
2733 else:
2734 return num.concatenate(
2735 (latlon, points[:, 2].reshape((len(points), 1))),
2736 axis=1)
2738 def points_on_source(self, cs='xyz', **kwargs):
2739 '''
2740 Convert relative plane coordinates to geographical coordinates.
2742 Given x and y coordinates (relative source coordinates between -1.
2743 and 1.) are converted to desired geographical coordinates. Coordinates
2744 need to be given as :py:class:`~numpy.ndarray` arguments ``points_x``
2745 and ``points_y``.
2747 :param cs:
2748 :ref:`Output coordinate system <coordinate-system-names>`.
2749 :type cs:
2750 optional, str
2752 :returns:
2753 Point coordinates in desired coordinate system.
2754 :rtype:
2755 :py:class:`~numpy.ndarray`: ``(n_points, [2, 3])``.
2756 '''
2757 points = points_on_rect_source(
2758 self.strike, self.dip, self.length, self.width,
2759 self.anchor, **kwargs)
2761 points[:, 0] += self.north_shift
2762 points[:, 1] += self.east_shift
2763 points[:, 2] += self.depth
2764 if cs == 'xyz':
2765 return points
2766 elif cs == 'xy':
2767 return points[:, :2]
2768 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2769 latlon = ne_to_latlon(
2770 self.lat, self.lon, points[:, 0], points[:, 1])
2772 latlon = num.array(latlon).T
2773 if cs == 'latlon':
2774 return latlon
2775 elif cs == 'lonlat':
2776 return latlon[:, ::-1]
2777 else:
2778 return num.concatenate(
2779 (latlon, points[:, 2].reshape((len(points), 1))),
2780 axis=1)
2782 def pyrocko_moment_tensor(self, store=None, target=None):
2783 if store is not None:
2784 if not self.patches:
2785 self.discretize_patches(store)
2787 data = self.get_slip()
2788 else:
2789 data = self.get_tractions()
2791 weights = num.linalg.norm(data, axis=1)
2792 weights /= weights.sum()
2794 rakes = num.arctan2(data[:, 1], data[:, 0]) * r2d
2795 rake = num.average(rakes, weights=weights)
2797 return pmt.MomentTensor(
2798 strike=self.strike,
2799 dip=self.dip,
2800 rake=rake,
2801 scalar_moment=self.get_moment(store, target))
2803 def pyrocko_event(self, store=None, target=None, **kwargs):
2804 return SourceWithDerivedMagnitude.pyrocko_event(
2805 self, store, target,
2806 **kwargs)
2808 @classmethod
2809 def from_pyrocko_event(cls, ev, **kwargs):
2810 d = {}
2811 mt = ev.moment_tensor
2812 if mt:
2813 (strike, dip, rake), _ = mt.both_strike_dip_rake()
2814 d.update(
2815 strike=float(strike),
2816 dip=float(dip),
2817 rake=float(rake))
2819 d.update(kwargs)
2820 return super(PseudoDynamicRupture, cls).from_pyrocko_event(ev, **d)
2822 def _discretize_points(self, store, *args, **kwargs):
2823 '''
2824 Discretize source plane with equal vertical and horizontal spacing.
2826 Additional ``*args`` and ``**kwargs`` are passed to
2827 :py:meth:`points_on_source`.
2829 :param store:
2830 Green's function database (needs to cover whole region of the
2831 source).
2832 :type store:
2833 :py:class:`~pyrocko.gf.store.Store`
2835 :returns:
2836 Number of points in strike and dip direction, distance
2837 between adjacent points, coordinates (latlondepth) and coordinates
2838 (xy on fault) for discrete points.
2839 :rtype:
2840 (int, int, float, :py:class:`~numpy.ndarray`,
2841 :py:class:`~numpy.ndarray`).
2842 '''
2843 anch_x, anch_y = map_anchor[self.anchor]
2845 npoints = int(self.width // km) + 1
2846 points = num.zeros((npoints, 3))
2847 points[:, 1] = num.linspace(-1., 1., npoints)
2848 points[:, 1] = (points[:, 1] - anch_y) * self.width/2
2850 rotmat = pmt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.0)
2851 points = num.dot(rotmat.T, points.T).T
2852 points[:, 2] += self.depth
2854 vs_min = store.config.get_vs(
2855 self.lat, self.lon, points,
2856 interpolation='nearest_neighbor')
2857 vr_min = max(vs_min.min(), .5*km) * self.gamma
2859 oversampling = 10.
2860 delta_l = self.length / (self.nx * oversampling)
2861 delta_w = self.width / (self.ny * oversampling)
2863 delta = self.eikonal_decimation * num.min([
2864 store.config.deltat * vr_min / oversampling,
2865 delta_l, delta_w] + [
2866 deltas for deltas in store.config.deltas])
2868 delta = delta_w / num.ceil(delta_w / delta)
2870 nx = int(num.ceil(self.length / delta)) + 1
2871 ny = int(num.ceil(self.width / delta)) + 1
2873 rem_l = (nx-1)*delta - self.length
2874 lim_x = rem_l / self.length
2876 points_xy = num.zeros((nx * ny, 2))
2877 points_xy[:, 0] = num.repeat(
2878 num.linspace(-1.-lim_x, 1.+lim_x, nx), ny)
2879 points_xy[:, 1] = num.tile(
2880 num.linspace(-1., 1., ny), nx)
2882 points = self.points_on_source(
2883 points_x=points_xy[:, 0],
2884 points_y=points_xy[:, 1],
2885 **kwargs)
2887 return nx, ny, delta, points, points_xy
2889 def _discretize_rupture_v(self, store, interpolation='nearest_neighbor',
2890 points=None):
2891 '''
2892 Get rupture velocity for discrete points on source plane.
2894 :param store:
2895 Green's function database (needs to cover the whole region of the
2896 source)
2897 :type store:
2898 optional, :py:class:`~pyrocko.gf.store.Store`
2900 :param interpolation:
2901 Interpolation method to use (choose between ``'nearest_neighbor'``
2902 and ``'multilinear'``).
2903 :type interpolation:
2904 optional, str
2906 :param points:
2907 Coordinates on fault (-1.:1.) of discrete points.
2908 :type points:
2909 optional, :py:class:`~numpy.ndarray`: ``(n_points, 2)``
2911 :returns:
2912 Rupture velocity assumed as :math:`v_s * \\gamma` for discrete
2913 points.
2914 :rtype:
2915 :py:class:`~numpy.ndarray`: ``(n_points, )``.
2916 '''
2918 if points is None:
2919 _, _, _, points, _ = self._discretize_points(store, cs='xyz')
2921 return store.config.get_vs(
2922 self.lat, self.lon,
2923 points=points,
2924 interpolation=interpolation) * self.gamma
2926 def discretize_time(
2927 self, store, interpolation='nearest_neighbor',
2928 vr=None, times=None, *args, **kwargs):
2929 '''
2930 Get rupture start time for discrete points on source plane.
2932 :param store:
2933 Green's function database (needs to cover whole region of the
2934 source)
2935 :type store:
2936 :py:class:`~pyrocko.gf.store.Store`
2938 :param interpolation:
2939 Interpolation method to use (choose between ``'nearest_neighbor'``
2940 and ``'multilinear'``).
2941 :type interpolation:
2942 optional, str
2944 :param vr:
2945 Array, containing rupture user defined rupture velocity values.
2946 :type vr:
2947 optional, :py:class:`~numpy.ndarray`
2949 :param times:
2950 Array, containing zeros, where rupture is starting, real positive
2951 numbers at later secondary nucleation points and -1, where time
2952 will be calculated. If not given, rupture starts at nucleation_x,
2953 nucleation_y. Times are given for discrete points with equal
2954 horizontal and vertical spacing.
2955 :type times:
2956 optional, :py:class:`~numpy.ndarray`
2958 :returns:
2959 Coordinates (latlondepth), coordinates (xy), rupture velocity,
2960 rupture propagation time of discrete points.
2961 :rtype:
2962 :py:class:`~numpy.ndarray`: ``(n_points, 3)``,
2963 :py:class:`~numpy.ndarray`: ``(n_points, 2)``,
2964 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``,
2965 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``.
2966 '''
2967 nx, ny, delta, points, points_xy = self._discretize_points(
2968 store, cs='xyz')
2970 if vr is None or vr.shape != tuple((nx, ny)):
2971 if vr:
2972 logger.warning(
2973 'Given rupture velocities are not in right shape: '
2974 '(%i, %i), but needed is (%i, %i).', *vr.shape + (nx, ny))
2975 vr = self._discretize_rupture_v(store, interpolation, points)\
2976 .reshape(nx, ny)
2978 if vr.shape != tuple((nx, ny)):
2979 logger.warning(
2980 'Given rupture velocities are not in right shape. Therefore'
2981 ' standard rupture velocity array is used.')
2983 def initialize_times():
2984 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
2986 if nucl_x.shape != nucl_y.shape:
2987 raise ValueError(
2988 'Nucleation coordinates have different shape.')
2990 dist_points = num.array([
2991 num.linalg.norm(points_xy - num.array([x, y]).ravel(), axis=1)
2992 for x, y in zip(nucl_x, nucl_y)])
2993 nucl_indices = num.argmin(dist_points, axis=1)
2995 if self.nucleation_time is None:
2996 nucl_times = num.zeros_like(nucl_indices)
2997 else:
2998 if self.nucleation_time.shape == nucl_x.shape:
2999 nucl_times = self.nucleation_time
3000 else:
3001 raise ValueError(
3002 'Nucleation coordinates and times have different '
3003 'shapes')
3005 t = num.full(nx * ny, -1.)
3006 t[nucl_indices] = nucl_times
3007 return t.reshape(nx, ny)
3009 if times is None:
3010 times = initialize_times()
3011 elif times.shape != tuple((nx, ny)):
3012 times = initialize_times()
3013 logger.warning(
3014 'Given times are not in right shape. Therefore standard time'
3015 ' array is used.')
3017 eikonal_ext.eikonal_solver_fmm_cartesian(
3018 speeds=vr, times=times, delta=delta)
3020 return points, points_xy, vr, times
3022 def get_vr_time_interpolators(
3023 self, store, interpolation='nearest_neighbor', force=False,
3024 **kwargs):
3025 '''
3026 Get interpolators for rupture velocity and rupture time.
3028 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3030 :param store:
3031 Green's function database (needs to cover whole region of the
3032 source).
3033 :type store:
3034 :py:class:`~pyrocko.gf.store.Store`
3036 :param interpolation:
3037 Interpolation method to use (choose between ``'nearest_neighbor'``
3038 and ``'multilinear'``).
3039 :type interpolation:
3040 optional, str
3042 :param force:
3043 Force recalculation of the interpolators (e.g. after change of
3044 nucleation point locations/times). Default is ``False``.
3045 :type force:
3046 optional, bool
3047 '''
3048 interp_map = {'multilinear': 'linear', 'nearest_neighbor': 'nearest'}
3049 if interpolation not in interp_map:
3050 raise TypeError(
3051 'Interpolation method %s not available' % interpolation)
3053 if not self._interpolators.get(interpolation, False) or force:
3054 _, points_xy, vr, times = self.discretize_time(
3055 store, **kwargs)
3057 if self.length <= 0.:
3058 raise ValueError(
3059 'length must be larger then 0. not %g' % self.length)
3061 if self.width <= 0.:
3062 raise ValueError(
3063 'width must be larger then 0. not %g' % self.width)
3065 nx, ny = times.shape
3066 anch_x, anch_y = map_anchor[self.anchor]
3068 points_xy[:, 0] = (points_xy[:, 0] - anch_x) * self.length / 2.
3069 points_xy[:, 1] = (points_xy[:, 1] - anch_y) * self.width / 2.
3071 self._interpolators[interpolation] = (
3072 nx, ny, times, vr,
3073 RegularGridInterpolator(
3074 (points_xy[::ny, 0], points_xy[:ny, 1]), times,
3075 method=interp_map[interpolation]),
3076 RegularGridInterpolator(
3077 (points_xy[::ny, 0], points_xy[:ny, 1]), vr,
3078 method=interp_map[interpolation]))
3079 return self._interpolators[interpolation]
3081 def discretize_patches(
3082 self, store, interpolation='nearest_neighbor', force=False,
3083 grid_shape=(),
3084 **kwargs):
3085 '''
3086 Get rupture start time and OkadaSource elements for points on rupture.
3088 All source elements and their corresponding center points are
3089 calculated and stored in the :py:attr:`patches` attribute.
3091 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3093 :param store:
3094 Green's function database (needs to cover whole region of the
3095 source).
3096 :type store:
3097 :py:class:`~pyrocko.gf.store.Store`
3099 :param interpolation:
3100 Interpolation method to use (choose between ``'nearest_neighbor'``
3101 and ``'multilinear'``).
3102 :type interpolation:
3103 optional, str
3105 :param force:
3106 Force recalculation of the vr and time interpolators ( e.g. after
3107 change of nucleation point locations/times). Default is ``False``.
3108 :type force:
3109 optional, bool
3111 :param grid_shape:
3112 Desired sub fault patch grid size (nlength, nwidth). Either factor
3113 or grid_shape should be set.
3114 :type grid_shape:
3115 optional, tuple of int
3116 '''
3117 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3118 self.get_vr_time_interpolators(
3119 store,
3120 interpolation=interpolation, force=force, **kwargs)
3121 anch_x, anch_y = map_anchor[self.anchor]
3123 al = self.length / 2.
3124 aw = self.width / 2.
3125 al1 = -(al + anch_x * al)
3126 al2 = al - anch_x * al
3127 aw1 = -aw + anch_y * aw
3128 aw2 = aw + anch_y * aw
3129 assert num.abs([al1, al2]).sum() == self.length
3130 assert num.abs([aw1, aw2]).sum() == self.width
3132 def get_lame(*a, **kw):
3133 shear_mod = store.config.get_shear_moduli(*a, **kw)
3134 lamb = store.config.get_vp(*a, **kw)**2 \
3135 * store.config.get_rho(*a, **kw) - 2. * shear_mod
3136 return shear_mod, lamb / (2. * (lamb + shear_mod))
3138 shear_mod, poisson = get_lame(
3139 self.lat, self.lon,
3140 num.array([[self.north_shift, self.east_shift, self.depth]]),
3141 interpolation=interpolation)
3143 okada_src = OkadaSource(
3144 lat=self.lat, lon=self.lon,
3145 strike=self.strike, dip=self.dip,
3146 north_shift=self.north_shift, east_shift=self.east_shift,
3147 depth=self.depth,
3148 al1=al1, al2=al2, aw1=aw1, aw2=aw2,
3149 poisson=poisson.mean(),
3150 shearmod=shear_mod.mean(),
3151 opening=kwargs.get('opening', 0.))
3153 if not (self.nx and self.ny):
3154 if grid_shape:
3155 self.nx, self.ny = grid_shape
3156 else:
3157 self.nx = nx
3158 self.ny = ny
3160 source_disc, source_points = okada_src.discretize(self.nx, self.ny)
3162 shear_mod, poisson = get_lame(
3163 self.lat, self.lon,
3164 num.array([src.source_patch()[:3] for src in source_disc]),
3165 interpolation=interpolation)
3167 if (self.nx, self.ny) != (nx, ny):
3168 times_interp = time_interpolator(source_points[:, :2])
3169 vr_interp = vr_interpolator(source_points[:, :2])
3170 else:
3171 times_interp = times.T.ravel()
3172 vr_interp = vr.T.ravel()
3174 for isrc, src in enumerate(source_disc):
3175 src.vr = vr_interp[isrc]
3176 src.time = times_interp[isrc] + self.time
3178 self.patches = source_disc
3180 def discretize_basesource(self, store, target=None):
3181 '''
3182 Prepare source for synthetic waveform calculation.
3184 :param store:
3185 Green's function database (needs to cover whole region of the
3186 source).
3187 :type store:
3188 :py:class:`~pyrocko.gf.store.Store`
3190 :param target:
3191 Target information.
3192 :type target:
3193 optional, :py:class:`~pyrocko.gf.targets.Target`
3195 :returns:
3196 Source discretized by a set of moment tensors and times.
3197 :rtype:
3198 :py:class:`~pyrocko.gf.meta.DiscretizedMTSource`
3199 '''
3200 if not target:
3201 interpolation = 'nearest_neighbor'
3202 else:
3203 interpolation = target.interpolation
3205 if not self.patches:
3206 self.discretize_patches(store, interpolation)
3208 if self.coef_mat is None:
3209 self.calc_coef_mat()
3211 delta_slip, slip_times = self.get_delta_slip(store)
3212 npatches = self.nx * self.ny
3213 ntimes = slip_times.size
3215 anch_x, anch_y = map_anchor[self.anchor]
3217 pln = self.length / self.nx
3218 pwd = self.width / self.ny
3220 patch_coords = num.array([
3221 (p.ix, p.iy)
3222 for p in self.patches]).reshape(self.nx, self.ny, 2)
3224 # boundary condition is zero-slip
3225 # is not valid to avoid unwished interpolation effects
3226 slip_grid = num.zeros((self.nx + 2, self.ny + 2, ntimes, 3))
3227 slip_grid[1:-1, 1:-1, :, :] = \
3228 delta_slip.reshape(self.nx, self.ny, ntimes, 3)
3230 slip_grid[0, 0, :, :] = slip_grid[1, 1, :, :]
3231 slip_grid[0, -1, :, :] = slip_grid[1, -2, :, :]
3232 slip_grid[-1, 0, :, :] = slip_grid[-2, 1, :, :]
3233 slip_grid[-1, -1, :, :] = slip_grid[-2, -2, :, :]
3235 slip_grid[1:-1, 0, :, :] = slip_grid[1:-1, 1, :, :]
3236 slip_grid[1:-1, -1, :, :] = slip_grid[1:-1, -2, :, :]
3237 slip_grid[0, 1:-1, :, :] = slip_grid[1, 1:-1, :, :]
3238 slip_grid[-1, 1:-1, :, :] = slip_grid[-2, 1:-1, :, :]
3240 def make_grid(patch_parameter):
3241 grid = num.zeros((self.nx + 2, self.ny + 2))
3242 grid[1:-1, 1:-1] = patch_parameter.reshape(self.nx, self.ny)
3244 grid[0, 0] = grid[1, 1]
3245 grid[0, -1] = grid[1, -2]
3246 grid[-1, 0] = grid[-2, 1]
3247 grid[-1, -1] = grid[-2, -2]
3249 grid[1:-1, 0] = grid[1:-1, 1]
3250 grid[1:-1, -1] = grid[1:-1, -2]
3251 grid[0, 1:-1] = grid[1, 1:-1]
3252 grid[-1, 1:-1] = grid[-2, 1:-1]
3254 return grid
3256 lamb = self.get_patch_attribute('lamb')
3257 mu = self.get_patch_attribute('shearmod')
3259 lamb_grid = make_grid(lamb)
3260 mu_grid = make_grid(mu)
3262 coords_x = num.zeros(self.nx + 2)
3263 coords_x[1:-1] = patch_coords[:, 0, 0]
3264 coords_x[0] = coords_x[1] - pln / 2
3265 coords_x[-1] = coords_x[-2] + pln / 2
3267 coords_y = num.zeros(self.ny + 2)
3268 coords_y[1:-1] = patch_coords[0, :, 1]
3269 coords_y[0] = coords_y[1] - pwd / 2
3270 coords_y[-1] = coords_y[-2] + pwd / 2
3272 slip_interp = RegularGridInterpolator(
3273 (coords_x, coords_y, slip_times),
3274 slip_grid, method='nearest')
3276 lamb_interp = RegularGridInterpolator(
3277 (coords_x, coords_y),
3278 lamb_grid, method='nearest')
3280 mu_interp = RegularGridInterpolator(
3281 (coords_x, coords_y),
3282 mu_grid, method='nearest')
3284 # discretize basesources
3285 mindeltagf = min(tuple(
3286 (self.length / self.nx, self.width / self.ny) +
3287 tuple(store.config.deltas)))
3289 nl = int((1. / self.decimation_factor) *
3290 num.ceil(pln / mindeltagf)) + 1
3291 nw = int((1. / self.decimation_factor) *
3292 num.ceil(pwd / mindeltagf)) + 1
3293 nsrc_patch = int(nl * nw)
3294 dl = pln / nl
3295 dw = pwd / nw
3297 patch_area = dl * dw
3299 xl = num.linspace(-0.5 * (pln - dl), 0.5 * (pln - dl), nl)
3300 xw = num.linspace(-0.5 * (pwd - dw), 0.5 * (pwd - dw), nw)
3302 base_coords = num.zeros((nsrc_patch, 3))
3303 base_coords[:, 0] = num.tile(xl, nw)
3304 base_coords[:, 1] = num.repeat(xw, nl)
3305 base_coords = num.tile(base_coords, (npatches, 1))
3307 center_coords = num.zeros((npatches, 3))
3308 center_coords[:, 0] = num.repeat(
3309 num.arange(self.nx) * pln + pln / 2, self.ny) - self.length / 2
3310 center_coords[:, 1] = num.tile(
3311 num.arange(self.ny) * pwd + pwd / 2, self.nx) - self.width / 2
3313 base_coords -= center_coords.repeat(nsrc_patch, axis=0)
3314 nbaselocs = base_coords.shape[0]
3316 base_interp = base_coords.repeat(ntimes, axis=0)
3318 base_times = num.tile(slip_times, nbaselocs)
3319 base_interp[:, 0] -= anch_x * self.length / 2
3320 base_interp[:, 1] -= anch_y * self.width / 2
3321 base_interp[:, 2] = base_times
3323 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3324 store, interpolation=interpolation)
3326 time_eikonal_max = time_interpolator.values.max()
3328 nbasesrcs = base_interp.shape[0]
3329 delta_slip = slip_interp(base_interp).reshape(nbaselocs, ntimes, 3)
3330 lamb = lamb_interp(base_interp[:, :2]).ravel()
3331 mu = mu_interp(base_interp[:, :2]).ravel()
3333 if False:
3334 try:
3335 import matplotlib.pyplot as plt
3336 coords = base_coords.copy()
3337 norm = num.sum(num.linalg.norm(delta_slip, axis=2), axis=1)
3338 plt.scatter(coords[:, 0], coords[:, 1], c=norm)
3339 plt.show()
3340 except AttributeError:
3341 pass
3343 base_interp[:, 2] = 0.
3344 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
3345 base_interp = num.dot(rotmat.T, base_interp.T).T
3346 base_interp[:, 0] += self.north_shift
3347 base_interp[:, 1] += self.east_shift
3348 base_interp[:, 2] += self.depth
3350 slip_strike = delta_slip[:, :, 0].ravel()
3351 slip_dip = delta_slip[:, :, 1].ravel()
3352 slip_norm = delta_slip[:, :, 2].ravel()
3354 slip_shear = num.linalg.norm([slip_strike, slip_dip], axis=0)
3355 slip_rake = r2d * num.arctan2(slip_dip, slip_strike)
3357 m6s = okada_ext.patch2m6(
3358 strikes=num.full(nbasesrcs, self.strike, dtype=float),
3359 dips=num.full(nbasesrcs, self.dip, dtype=float),
3360 rakes=slip_rake,
3361 disl_shear=slip_shear,
3362 disl_norm=slip_norm,
3363 lamb=lamb,
3364 mu=mu,
3365 nthreads=self.nthreads)
3367 m6s *= patch_area
3369 dl = -self.patches[0].al1 + self.patches[0].al2
3370 dw = -self.patches[0].aw1 + self.patches[0].aw2
3372 base_times[base_times > time_eikonal_max] = time_eikonal_max
3374 ds = meta.DiscretizedMTSource(
3375 lat=self.lat,
3376 lon=self.lon,
3377 times=base_times + self.time,
3378 north_shifts=base_interp[:, 0],
3379 east_shifts=base_interp[:, 1],
3380 depths=base_interp[:, 2],
3381 m6s=m6s,
3382 dl=dl,
3383 dw=dw,
3384 nl=self.nx,
3385 nw=self.ny)
3387 return ds
3389 def calc_coef_mat(self):
3390 '''
3391 Calculate coefficients connecting tractions and dislocations.
3392 '''
3393 if not self.patches:
3394 raise ValueError(
3395 'Patches are needed. Please calculate them first.')
3397 self.coef_mat = make_okada_coefficient_matrix(
3398 self.patches, nthreads=self.nthreads, pure_shear=self.pure_shear)
3400 def get_patch_attribute(self, attr):
3401 '''
3402 Get patch attributes.
3404 :param attr:
3405 Name of selected attribute (see
3406 :py:class`pyrocko.modelling.okada.OkadaSource`).
3407 :type attr:
3408 str
3410 :returns:
3411 Array with attribute value for each fault patch.
3412 :rtype:
3413 :py:class:`~numpy.ndarray`
3415 '''
3416 if not self.patches:
3417 raise ValueError(
3418 'Patches are needed. Please calculate them first.')
3419 return num.array([getattr(p, attr) for p in self.patches])
3421 def get_slip(
3422 self,
3423 time=None,
3424 scale_slip=True,
3425 interpolation='nearest_neighbor',
3426 **kwargs):
3427 '''
3428 Get slip per subfault patch for given time after rupture start.
3430 :param time:
3431 Time after origin [s], for which slip is computed. If not
3432 given, final static slip is returned.
3433 :type time:
3434 optional, float > 0.
3436 :param scale_slip:
3437 If ``True`` and :py:attr:`slip` given, all slip values are scaled
3438 to fit the given maximum slip.
3439 :type scale_slip:
3440 optional, bool
3442 :param interpolation:
3443 Interpolation method to use (choose between ``'nearest_neighbor'``
3444 and ``'multilinear'``).
3445 :type interpolation:
3446 optional, str
3448 :returns:
3449 Inverted dislocations (:math:`u_{strike}, u_{dip}, u_{tensile}`)
3450 for each source patch.
3451 :rtype:
3452 :py:class:`~numpy.ndarray`: ``(n_sources, 3)``
3453 '''
3455 if self.patches is None:
3456 raise ValueError(
3457 'Please discretize the source first (discretize_patches())')
3458 npatches = len(self.patches)
3459 tractions = self.get_tractions()
3460 time_patch_max = self.get_patch_attribute('time').max() - self.time
3462 time_patch = time
3463 if time is None:
3464 time_patch = time_patch_max
3466 if self.coef_mat is None:
3467 self.calc_coef_mat()
3469 if tractions.shape != (npatches, 3):
3470 raise AttributeError(
3471 'The traction vector is of invalid shape.'
3472 ' Required shape is (npatches, 3)')
3474 patch_mask = num.ones(npatches, dtype=bool)
3475 if self.patch_mask is not None:
3476 patch_mask = self.patch_mask
3478 times = self.get_patch_attribute('time') - self.time
3479 times[~patch_mask] = time_patch + 1. # exlcude unmasked patches
3480 relevant_sources = num.nonzero(times <= time_patch)[0]
3481 disloc_est = num.zeros_like(tractions)
3483 if self.smooth_rupture:
3484 patch_activation = num.zeros(npatches)
3486 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3487 self.get_vr_time_interpolators(
3488 store, interpolation=interpolation)
3490 # Getting the native Eikonal grid, bit hackish
3491 points_x = num.round(time_interpolator.grid[0], decimals=2)
3492 points_y = num.round(time_interpolator.grid[1], decimals=2)
3493 times_eikonal = time_interpolator.values
3495 time_max = time
3496 if time is None:
3497 time_max = times_eikonal.max()
3499 for ip, p in enumerate(self.patches):
3500 ul = num.round((p.ix + p.al1, p.iy + p.aw1), decimals=2)
3501 lr = num.round((p.ix + p.al2, p.iy + p.aw2), decimals=2)
3503 idx_length = num.logical_and(
3504 points_x >= ul[0], points_x <= lr[0])
3505 idx_width = num.logical_and(
3506 points_y >= ul[1], points_y <= lr[1])
3508 times_patch = times_eikonal[num.ix_(idx_length, idx_width)]
3509 if times_patch.size == 0:
3510 raise AttributeError('could not use smooth_rupture')
3512 patch_activation[ip] = \
3513 (times_patch <= time_max).sum() / times_patch.size
3515 if time_patch == 0 and time_patch != time_patch_max:
3516 patch_activation[ip] = 0.
3518 patch_activation[~patch_mask] = 0. # exlcude unmasked patches
3520 relevant_sources = num.nonzero(patch_activation > 0.)[0]
3522 if relevant_sources.size == 0:
3523 return disloc_est
3525 indices_disl = num.repeat(relevant_sources * 3, 3)
3526 indices_disl[1::3] += 1
3527 indices_disl[2::3] += 2
3529 disloc_est[relevant_sources] = invert_fault_dislocations_bem(
3530 stress_field=tractions[relevant_sources, :].ravel(),
3531 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3532 pure_shear=self.pure_shear, nthreads=self.nthreads,
3533 epsilon=None,
3534 **kwargs)
3536 if self.smooth_rupture:
3537 disloc_est *= patch_activation[:, num.newaxis]
3539 if scale_slip and self.slip is not None:
3540 disloc_tmax = num.zeros(npatches)
3542 indices_disl = num.repeat(num.nonzero(patch_mask)[0] * 3, 3)
3543 indices_disl[1::3] += 1
3544 indices_disl[2::3] += 2
3546 disloc_tmax[patch_mask] = num.linalg.norm(
3547 invert_fault_dislocations_bem(
3548 stress_field=tractions[patch_mask, :].ravel(),
3549 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3550 pure_shear=self.pure_shear, nthreads=self.nthreads,
3551 epsilon=None,
3552 **kwargs), axis=1)
3554 disloc_tmax_max = disloc_tmax.max()
3555 if disloc_tmax_max == 0.:
3556 logger.warning(
3557 'slip scaling not performed. Maximum slip is 0.')
3559 disloc_est *= self.slip / disloc_tmax_max
3561 return disloc_est
3563 def get_delta_slip(
3564 self,
3565 store=None,
3566 deltat=None,
3567 delta=True,
3568 interpolation='nearest_neighbor',
3569 **kwargs):
3570 '''
3571 Get slip change snapshots.
3573 The time interval, within which the slip changes are computed is
3574 determined by the sampling rate of the Green's function database or
3575 ``deltat``. Additional ``**kwargs`` are passed to :py:meth:`get_slip`.
3577 :param store:
3578 Green's function database (needs to cover whole region of of the
3579 source). Its sampling interval is used as time increment for slip
3580 difference calculation. Either ``deltat`` or ``store`` should be
3581 given.
3582 :type store:
3583 optional, :py:class:`~pyrocko.gf.store.Store`
3585 :param deltat:
3586 Time interval for slip difference calculation [s]. Either
3587 ``deltat`` or ``store`` should be given.
3588 :type deltat:
3589 optional, float
3591 :param delta:
3592 If ``True``, slip differences between two time steps are given. If
3593 ``False``, cumulative slip for all time steps.
3594 :type delta:
3595 optional, bool
3597 :param interpolation:
3598 Interpolation method to use (choose between ``'nearest_neighbor'``
3599 and ``'multilinear'``).
3600 :type interpolation:
3601 optional, str
3603 :returns:
3604 Displacement changes(:math:`\\Delta u_{strike},
3605 \\Delta u_{dip} , \\Delta u_{tensile}`) for each source patch and
3606 time; corner times, for which delta slip is computed. The order of
3607 displacement changes array is:
3609 .. math::
3611 &[[\\\\
3612 &[\\Delta u_{strike, patch1, t1},
3613 \\Delta u_{dip, patch1, t1},
3614 \\Delta u_{tensile, patch1, t1}],\\\\
3615 &[\\Delta u_{strike, patch1, t2},
3616 \\Delta u_{dip, patch1, t2},
3617 \\Delta u_{tensile, patch1, t2}]\\\\
3618 &], [\\\\
3619 &[\\Delta u_{strike, patch2, t1}, ...],\\\\
3620 &[\\Delta u_{strike, patch2, t2}, ...]]]\\\\
3622 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3623 :py:class:`~numpy.ndarray`: ``(n_times, )``
3624 '''
3625 if store and deltat:
3626 raise AttributeError(
3627 'Argument collision. '
3628 'Please define only the store or the deltat argument.')
3630 if store:
3631 deltat = store.config.deltat
3633 if not deltat:
3634 raise AttributeError('Please give a GF store or set deltat.')
3636 npatches = len(self.patches)
3638 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3639 store, interpolation=interpolation)
3640 tmax = time_interpolator.values.max()
3642 calc_times = num.arange(0., tmax + deltat, deltat)
3643 calc_times[calc_times > tmax] = tmax
3645 disloc_est = num.zeros((npatches, calc_times.size, 3))
3647 for itime, t in enumerate(calc_times):
3648 disloc_est[:, itime, :] = self.get_slip(
3649 time=t, scale_slip=False, **kwargs)
3651 if self.slip:
3652 disloc_tmax = num.linalg.norm(
3653 self.get_slip(scale_slip=False, time=tmax),
3654 axis=1)
3656 disloc_tmax_max = disloc_tmax.max()
3657 if disloc_tmax_max == 0.:
3658 logger.warning(
3659 'Slip scaling not performed. Maximum slip is 0.')
3660 else:
3661 disloc_est *= self.slip / disloc_tmax_max
3663 if not delta:
3664 return disloc_est, calc_times
3666 # if we have only one timestep there is no gradient
3667 if calc_times.size > 1:
3668 disloc_init = disloc_est[:, 0, :]
3669 disloc_est = num.diff(disloc_est, axis=1)
3670 disloc_est = num.concatenate((
3671 disloc_init[:, num.newaxis, :], disloc_est), axis=1)
3673 calc_times = calc_times
3675 return disloc_est, calc_times
3677 def get_slip_rate(self, *args, **kwargs):
3678 '''
3679 Get slip rate inverted from patches.
3681 The time interval, within which the slip rates are computed is
3682 determined by the sampling rate of the Green's function database or
3683 ``deltat``. Additional ``*args`` and ``**kwargs`` are passed to
3684 :py:meth:`get_delta_slip`.
3686 :returns:
3687 Slip rates (:math:`\\Delta u_{strike}/\\Delta t`,
3688 :math:`\\Delta u_{dip}/\\Delta t, \\Delta u_{tensile}/\\Delta t`)
3689 for each source patch and time; corner times, for which slip rate
3690 is computed. The order of sliprate array is:
3692 .. math::
3694 &[[\\\\
3695 &[\\Delta u_{strike, patch1, t1}/\\Delta t,
3696 \\Delta u_{dip, patch1, t1}/\\Delta t,
3697 \\Delta u_{tensile, patch1, t1}/\\Delta t],\\\\
3698 &[\\Delta u_{strike, patch1, t2}/\\Delta t,
3699 \\Delta u_{dip, patch1, t2}/\\Delta t,
3700 \\Delta u_{tensile, patch1, t2}/\\Delta t]], [\\\\
3701 &[\\Delta u_{strike, patch2, t1}/\\Delta t, ...],\\\\
3702 &[\\Delta u_{strike, patch2, t2}/\\Delta t, ...]]]\\\\
3704 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3705 :py:class:`~numpy.ndarray`: ``(n_times, )``
3706 '''
3707 ddisloc_est, calc_times = self.get_delta_slip(
3708 *args, delta=True, **kwargs)
3710 dt = num.concatenate(
3711 [(num.diff(calc_times)[0], ), num.diff(calc_times)])
3712 slip_rate = num.linalg.norm(ddisloc_est, axis=2) / dt
3714 return slip_rate, calc_times
3716 def get_moment_rate_patches(self, *args, **kwargs):
3717 '''
3718 Get scalar seismic moment rate for each patch individually.
3720 Additional ``*args`` and ``**kwargs`` are passed to
3721 :py:meth:`get_slip_rate`.
3723 :returns:
3724 Seismic moment rate for each source patch and time; corner times,
3725 for which patch moment rate is computed based on slip rate. The
3726 order of the moment rate array is:
3728 .. math::
3730 &[\\\\
3731 &[(\\Delta M / \\Delta t)_{patch1, t1},
3732 (\\Delta M / \\Delta t)_{patch1, t2}, ...],\\\\
3733 &[(\\Delta M / \\Delta t)_{patch2, t1},
3734 (\\Delta M / \\Delta t)_{patch, t2}, ...],\\\\
3735 &[...]]\\\\
3737 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times)``,
3738 :py:class:`~numpy.ndarray`: ``(n_times, )``
3739 '''
3740 slip_rate, calc_times = self.get_slip_rate(*args, **kwargs)
3742 shear_mod = self.get_patch_attribute('shearmod')
3743 p_length = self.get_patch_attribute('length')
3744 p_width = self.get_patch_attribute('width')
3746 dA = p_length * p_width
3748 mom_rate = shear_mod[:, num.newaxis] * slip_rate * dA[:, num.newaxis]
3750 return mom_rate, calc_times
3752 def get_moment_rate(self, store, target=None, deltat=None):
3753 '''
3754 Get seismic source moment rate for the total source (STF).
3756 :param store:
3757 Green's function database (needs to cover whole region of of the
3758 source). Its ``deltat`` [s] is used as time increment for slip
3759 difference calculation. Either ``deltat`` or ``store`` should be
3760 given.
3761 :type store:
3762 :py:class:`~pyrocko.gf.store.Store`
3764 :param target:
3765 Target information, needed for interpolation method.
3766 :type target:
3767 optional, :py:class:`~pyrocko.gf.targets.Target`
3769 :param deltat:
3770 Time increment for slip difference calculation [s]. If not given
3771 ``store.deltat`` is used.
3772 :type deltat:
3773 optional, float
3775 :return:
3776 Seismic moment rate [Nm/s] for each time; corner times, for which
3777 moment rate is computed. The order of the moment rate array is:
3779 .. math::
3781 &[\\\\
3782 &(\\Delta M / \\Delta t)_{t1},\\\\
3783 &(\\Delta M / \\Delta t)_{t2},\\\\
3784 &...]\\\\
3786 :rtype:
3787 :py:class:`~numpy.ndarray`: ``(n_times, )``,
3788 :py:class:`~numpy.ndarray`: ``(n_times, )``
3789 '''
3790 if not deltat:
3791 deltat = store.config.deltat
3792 return self.discretize_basesource(
3793 store, target=target).get_moment_rate(deltat)
3795 def get_moment(self, *args, **kwargs):
3796 '''
3797 Get seismic cumulative moment.
3799 Additional ``*args`` and ``**kwargs`` are passed to
3800 :py:meth:`get_magnitude`.
3802 :returns:
3803 Cumulative seismic moment in [Nm].
3804 :rtype:
3805 float
3806 '''
3807 return float(pmt.magnitude_to_moment(self.get_magnitude(
3808 *args, **kwargs)))
3810 def rescale_slip(self, magnitude=None, moment=None, **kwargs):
3811 '''
3812 Rescale source slip based on given target magnitude or seismic moment.
3814 Rescale the maximum source slip to fit the source moment magnitude or
3815 seismic moment to the given target values. Either ``magnitude`` or
3816 ``moment`` need to be given. Additional ``**kwargs`` are passed to
3817 :py:meth:`get_moment`.
3819 :param magnitude:
3820 Target moment magnitude :math:`M_\\mathrm{w}` as in
3821 [Hanks and Kanamori, 1979]
3822 :type magnitude:
3823 optional, float
3825 :param moment:
3826 Target seismic moment :math:`M_0` [Nm].
3827 :type moment:
3828 optional, float
3829 '''
3830 if self.slip is None:
3831 self.slip = 1.
3832 logger.warning('No slip found for rescaling. '
3833 'An initial slip of 1 m is assumed.')
3835 if magnitude is None and moment is None:
3836 raise ValueError(
3837 'Either target magnitude or moment need to be given.')
3839 moment_init = self.get_moment(**kwargs)
3841 if magnitude is not None:
3842 moment = pmt.magnitude_to_moment(magnitude)
3844 self.slip *= moment / moment_init
3846 def get_centroid(self, store, *args, **kwargs):
3847 '''
3848 Centroid of the pseudo dynamic rupture model.
3850 The centroid location and time are derived from the locations and times
3851 of the individual patches weighted with their moment contribution.
3852 Additional ``**kwargs`` are passed to :py:meth:`pyrocko_moment_tensor`.
3854 :param store:
3855 Green's function database (needs to cover whole region of of the
3856 source). Its ``deltat`` [s] is used as time increment for slip
3857 difference calculation. Either ``deltat`` or ``store`` should be
3858 given.
3859 :type store:
3860 :py:class:`~pyrocko.gf.store.Store`
3862 :returns:
3863 The centroid location and associated moment tensor.
3864 :rtype:
3865 :py:class:`pyrocko.model.Event`
3866 '''
3867 _, _, _, _, time, _ = self.get_vr_time_interpolators(store)
3868 t_max = time.values.max()
3870 moment_rate, times = self.get_moment_rate_patches(deltat=t_max)
3872 moment = num.sum(moment_rate * times, axis=1)
3873 weights = moment / moment.sum()
3875 norths = self.get_patch_attribute('north_shift')
3876 easts = self.get_patch_attribute('east_shift')
3877 depths = self.get_patch_attribute('depth')
3879 centroid_n = num.sum(weights * norths)
3880 centroid_e = num.sum(weights * easts)
3881 centroid_d = num.sum(weights * depths)
3883 centroid_lat, centroid_lon = ne_to_latlon(
3884 self.lat, self.lon, centroid_n, centroid_e)
3886 moment_rate_, times = self.get_moment_rate(store)
3887 delta_times = num.concatenate((
3888 [times[1] - times[0]],
3889 num.diff(times)))
3890 moment_src = delta_times * moment_rate
3892 centroid_t = num.sum(
3893 moment_src / num.sum(moment_src) * times) + self.time
3895 mt = self.pyrocko_moment_tensor(store, *args, **kwargs)
3897 return model.Event(
3898 lat=centroid_lat,
3899 lon=centroid_lon,
3900 depth=centroid_d,
3901 time=centroid_t,
3902 moment_tensor=mt,
3903 magnitude=mt.magnitude,
3904 duration=t_max)
3907class DoubleDCSource(SourceWithMagnitude):
3908 '''
3909 Two double-couple point sources separated in space and time.
3910 Moment share between the sub-sources is controlled by the
3911 parameter mix.
3912 The position of the subsources is dependent on the moment
3913 distribution between the two sources. Depth, east and north
3914 shift are given for the centroid between the two double-couples.
3915 The subsources will positioned according to their moment shares
3916 around this centroid position.
3917 This is done according to their delta parameters, which are
3918 therefore in relation to that centroid.
3919 Note that depth of the subsources therefore can be
3920 depth+/-delta_depth. For shallow earthquakes therefore
3921 the depth has to be chosen deeper to avoid sampling
3922 above surface.
3923 '''
3925 strike1 = Float.T(
3926 default=0.0,
3927 help='strike direction in [deg], measured clockwise from north')
3929 dip1 = Float.T(
3930 default=90.0,
3931 help='dip angle in [deg], measured downward from horizontal')
3933 azimuth = Float.T(
3934 default=0.0,
3935 help='azimuth to second double-couple [deg], '
3936 'measured at first, clockwise from north')
3938 rake1 = Float.T(
3939 default=0.0,
3940 help='rake angle in [deg], '
3941 'measured counter-clockwise from right-horizontal '
3942 'in on-plane view')
3944 strike2 = Float.T(
3945 default=0.0,
3946 help='strike direction in [deg], measured clockwise from north')
3948 dip2 = Float.T(
3949 default=90.0,
3950 help='dip angle in [deg], measured downward from horizontal')
3952 rake2 = Float.T(
3953 default=0.0,
3954 help='rake angle in [deg], '
3955 'measured counter-clockwise from right-horizontal '
3956 'in on-plane view')
3958 delta_time = Float.T(
3959 default=0.0,
3960 help='separation of double-couples in time (t2-t1) [s]')
3962 delta_depth = Float.T(
3963 default=0.0,
3964 help='difference in depth (z2-z1) [m]')
3966 distance = Float.T(
3967 default=0.0,
3968 help='distance between the two double-couples [m]')
3970 mix = Float.T(
3971 default=0.5,
3972 help='how to distribute the moment to the two doublecouples '
3973 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1')
3975 stf1 = STF.T(
3976 optional=True,
3977 help='Source time function of subsource 1 '
3978 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
3980 stf2 = STF.T(
3981 optional=True,
3982 help='Source time function of subsource 2 '
3983 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
3985 discretized_source_class = meta.DiscretizedMTSource
3987 def base_key(self):
3988 return (
3989 self.time, self.depth, self.lat, self.north_shift,
3990 self.lon, self.east_shift, type(self).__name__) + \
3991 self.effective_stf1_pre().base_key() + \
3992 self.effective_stf2_pre().base_key() + (
3993 self.strike1, self.dip1, self.rake1,
3994 self.strike2, self.dip2, self.rake2,
3995 self.delta_time, self.delta_depth,
3996 self.azimuth, self.distance, self.mix)
3998 def get_factor(self):
3999 return self.moment
4001 def effective_stf1_pre(self):
4002 return self.stf1 or self.stf or g_unit_pulse
4004 def effective_stf2_pre(self):
4005 return self.stf2 or self.stf or g_unit_pulse
4007 def effective_stf_post(self):
4008 return g_unit_pulse
4010 def split(self):
4011 a1 = 1.0 - self.mix
4012 a2 = self.mix
4013 delta_north = math.cos(self.azimuth * d2r) * self.distance
4014 delta_east = math.sin(self.azimuth * d2r) * self.distance
4016 dc1 = DCSource(
4017 lat=self.lat,
4018 lon=self.lon,
4019 time=self.time - self.delta_time * a2,
4020 north_shift=self.north_shift - delta_north * a2,
4021 east_shift=self.east_shift - delta_east * a2,
4022 depth=self.depth - self.delta_depth * a2,
4023 moment=self.moment * a1,
4024 strike=self.strike1,
4025 dip=self.dip1,
4026 rake=self.rake1,
4027 stf=self.stf1 or self.stf)
4029 dc2 = DCSource(
4030 lat=self.lat,
4031 lon=self.lon,
4032 time=self.time + self.delta_time * a1,
4033 north_shift=self.north_shift + delta_north * a1,
4034 east_shift=self.east_shift + delta_east * a1,
4035 depth=self.depth + self.delta_depth * a1,
4036 moment=self.moment * a2,
4037 strike=self.strike2,
4038 dip=self.dip2,
4039 rake=self.rake2,
4040 stf=self.stf2 or self.stf)
4042 return [dc1, dc2]
4044 def discretize_basesource(self, store, target=None):
4045 a1 = 1.0 - self.mix
4046 a2 = self.mix
4047 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4048 rake=self.rake1, scalar_moment=a1)
4049 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4050 rake=self.rake2, scalar_moment=a2)
4052 delta_north = math.cos(self.azimuth * d2r) * self.distance
4053 delta_east = math.sin(self.azimuth * d2r) * self.distance
4055 times1, amplitudes1 = self.effective_stf1_pre().discretize_t(
4056 store.config.deltat, self.time - self.delta_time * a2)
4058 times2, amplitudes2 = self.effective_stf2_pre().discretize_t(
4059 store.config.deltat, self.time + self.delta_time * a1)
4061 nt1 = times1.size
4062 nt2 = times2.size
4064 ds = meta.DiscretizedMTSource(
4065 lat=self.lat,
4066 lon=self.lon,
4067 times=num.concatenate((times1, times2)),
4068 north_shifts=num.concatenate((
4069 num.repeat(self.north_shift - delta_north * a2, nt1),
4070 num.repeat(self.north_shift + delta_north * a1, nt2))),
4071 east_shifts=num.concatenate((
4072 num.repeat(self.east_shift - delta_east * a2, nt1),
4073 num.repeat(self.east_shift + delta_east * a1, nt2))),
4074 depths=num.concatenate((
4075 num.repeat(self.depth - self.delta_depth * a2, nt1),
4076 num.repeat(self.depth + self.delta_depth * a1, nt2))),
4077 m6s=num.vstack((
4078 mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis],
4079 mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis])))
4081 return ds
4083 def pyrocko_moment_tensor(self, store=None, target=None):
4084 a1 = 1.0 - self.mix
4085 a2 = self.mix
4086 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4087 rake=self.rake1,
4088 scalar_moment=a1 * self.moment)
4089 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4090 rake=self.rake2,
4091 scalar_moment=a2 * self.moment)
4092 return pmt.MomentTensor(m=mot1.m() + mot2.m())
4094 def pyrocko_event(self, store=None, target=None, **kwargs):
4095 return SourceWithMagnitude.pyrocko_event(
4096 self, store, target,
4097 moment_tensor=self.pyrocko_moment_tensor(store, target),
4098 **kwargs)
4100 @classmethod
4101 def from_pyrocko_event(cls, ev, **kwargs):
4102 d = {}
4103 mt = ev.moment_tensor
4104 if mt:
4105 (strike, dip, rake), _ = mt.both_strike_dip_rake()
4106 d.update(
4107 strike1=float(strike),
4108 dip1=float(dip),
4109 rake1=float(rake),
4110 strike2=float(strike),
4111 dip2=float(dip),
4112 rake2=float(rake),
4113 mix=0.0,
4114 magnitude=float(mt.moment_magnitude()))
4116 d.update(kwargs)
4117 source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d)
4118 source.stf1 = source.stf
4119 source.stf2 = HalfSinusoidSTF(effective_duration=0.)
4120 source.stf = None
4121 return source
4124class RingfaultSource(SourceWithMagnitude):
4125 '''
4126 A ring fault with vertical doublecouples.
4127 '''
4129 diameter = Float.T(
4130 default=1.0,
4131 help='diameter of the ring in [m]')
4133 sign = Float.T(
4134 default=1.0,
4135 help='inside of the ring moves up (+1) or down (-1)')
4137 strike = Float.T(
4138 default=0.0,
4139 help='strike direction of the ring plane, clockwise from north,'
4140 ' in [deg]')
4142 dip = Float.T(
4143 default=0.0,
4144 help='dip angle of the ring plane from horizontal in [deg]')
4146 npointsources = Int.T(
4147 default=360,
4148 help='number of point sources to use')
4150 discretized_source_class = meta.DiscretizedMTSource
4152 def base_key(self):
4153 return Source.base_key(self) + (
4154 self.strike, self.dip, self.diameter, self.npointsources)
4156 def get_factor(self):
4157 return self.sign * self.moment
4159 def discretize_basesource(self, store=None, target=None):
4160 n = self.npointsources
4161 phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False)
4163 points = num.zeros((n, 3))
4164 points[:, 0] = num.cos(phi) * 0.5 * self.diameter
4165 points[:, 1] = num.sin(phi) * 0.5 * self.diameter
4167 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
4168 points = num.dot(rotmat.T, points.T).T # !!! ?
4170 points[:, 0] += self.north_shift
4171 points[:, 1] += self.east_shift
4172 points[:, 2] += self.depth
4174 m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90.,
4175 scalar_moment=1.0 / n).m())
4177 rotmats = num.transpose(
4178 [[num.cos(phi), num.sin(phi), num.zeros(n)],
4179 [-num.sin(phi), num.cos(phi), num.zeros(n)],
4180 [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1))
4182 ms = num.zeros((n, 3, 3))
4183 for i in range(n):
4184 mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i]))
4185 ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat))
4187 m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2],
4188 ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T
4190 times, amplitudes = self.effective_stf_pre().discretize_t(
4191 store.config.deltat, self.time)
4193 nt = times.size
4195 return meta.DiscretizedMTSource(
4196 times=num.tile(times, n),
4197 lat=self.lat,
4198 lon=self.lon,
4199 north_shifts=num.repeat(points[:, 0], nt),
4200 east_shifts=num.repeat(points[:, 1], nt),
4201 depths=num.repeat(points[:, 2], nt),
4202 m6s=num.repeat(m6s, nt, axis=0) * num.tile(
4203 amplitudes, n)[:, num.newaxis])
4206class CombiSource(Source):
4207 '''
4208 Composite source model.
4209 '''
4211 discretized_source_class = meta.DiscretizedMTSource
4213 subsources = List.T(Source.T())
4215 def __init__(self, subsources=[], **kwargs):
4216 if not subsources:
4217 raise BadRequest(
4218 'Need at least one sub-source to create a CombiSource object.')
4220 lats = num.array(
4221 [subsource.lat for subsource in subsources], dtype=float)
4222 lons = num.array(
4223 [subsource.lon for subsource in subsources], dtype=float)
4225 lat, lon = lats[0], lons[0]
4226 if not num.all(lats == lat) and num.all(lons == lon):
4227 subsources = [s.clone() for s in subsources]
4228 for subsource in subsources[1:]:
4229 subsource.set_origin(lat, lon)
4231 depth = float(num.mean([p.depth for p in subsources]))
4232 time = float(num.mean([p.time for p in subsources]))
4233 north_shift = float(num.mean([p.north_shift for p in subsources]))
4234 east_shift = float(num.mean([p.east_shift for p in subsources]))
4235 kwargs.update(
4236 time=time,
4237 lat=float(lat),
4238 lon=float(lon),
4239 north_shift=north_shift,
4240 east_shift=east_shift,
4241 depth=depth)
4243 Source.__init__(self, subsources=subsources, **kwargs)
4245 def get_factor(self):
4246 return 1.0
4248 def discretize_basesource(self, store, target=None):
4249 dsources = []
4250 for sf in self.subsources:
4251 ds = sf.discretize_basesource(store, target)
4252 ds.m6s *= sf.get_factor()
4253 dsources.append(ds)
4255 return meta.DiscretizedMTSource.combine(dsources)
4258class SFSource(Source):
4259 '''
4260 A single force point source.
4262 Supported GF schemes: `'elastic5'`.
4263 '''
4265 discretized_source_class = meta.DiscretizedSFSource
4267 fn = Float.T(
4268 default=0.,
4269 help='northward component of single force [N]')
4271 fe = Float.T(
4272 default=0.,
4273 help='eastward component of single force [N]')
4275 fd = Float.T(
4276 default=0.,
4277 help='downward component of single force [N]')
4279 def __init__(self, **kwargs):
4280 Source.__init__(self, **kwargs)
4282 def base_key(self):
4283 return Source.base_key(self) + (self.fn, self.fe, self.fd)
4285 def get_factor(self):
4286 return 1.0
4288 def discretize_basesource(self, store, target=None):
4289 times, amplitudes = self.effective_stf_pre().discretize_t(
4290 store.config.deltat, self.time)
4291 forces = amplitudes[:, num.newaxis] * num.array(
4292 [[self.fn, self.fe, self.fd]], dtype=float)
4294 return meta.DiscretizedSFSource(forces=forces,
4295 **self._dparams_base_repeated(times))
4297 def pyrocko_event(self, store=None, target=None, **kwargs):
4298 return Source.pyrocko_event(
4299 self, store, target,
4300 **kwargs)
4302 @classmethod
4303 def from_pyrocko_event(cls, ev, **kwargs):
4304 d = {}
4305 d.update(kwargs)
4306 return super(SFSource, cls).from_pyrocko_event(ev, **d)
4309class PorePressurePointSource(Source):
4310 '''
4311 Excess pore pressure point source.
4313 For poro-elastic initial value problem where an excess pore pressure is
4314 brought into a small source volume.
4315 '''
4317 discretized_source_class = meta.DiscretizedPorePressureSource
4319 pp = Float.T(
4320 default=1.0,
4321 help='initial excess pore pressure in [Pa]')
4323 def base_key(self):
4324 return Source.base_key(self)
4326 def get_factor(self):
4327 return self.pp
4329 def discretize_basesource(self, store, target=None):
4330 return meta.DiscretizedPorePressureSource(pp=arr(1.0),
4331 **self._dparams_base())
4334class PorePressureLineSource(Source):
4335 '''
4336 Excess pore pressure line source.
4338 The line source is centered at (north_shift, east_shift, depth).
4339 '''
4341 discretized_source_class = meta.DiscretizedPorePressureSource
4343 pp = Float.T(
4344 default=1.0,
4345 help='initial excess pore pressure in [Pa]')
4347 length = Float.T(
4348 default=0.0,
4349 help='length of the line source [m]')
4351 azimuth = Float.T(
4352 default=0.0,
4353 help='azimuth direction, clockwise from north [deg]')
4355 dip = Float.T(
4356 default=90.,
4357 help='dip direction, downward from horizontal [deg]')
4359 def base_key(self):
4360 return Source.base_key(self) + (self.azimuth, self.dip, self.length)
4362 def get_factor(self):
4363 return self.pp
4365 def discretize_basesource(self, store, target=None):
4367 n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1
4369 a = num.linspace(-0.5 * self.length, 0.5 * self.length, n)
4371 sa = math.sin(self.azimuth * d2r)
4372 ca = math.cos(self.azimuth * d2r)
4373 sd = math.sin(self.dip * d2r)
4374 cd = math.cos(self.dip * d2r)
4376 points = num.zeros((n, 3))
4377 points[:, 0] = self.north_shift + a * ca * cd
4378 points[:, 1] = self.east_shift + a * sa * cd
4379 points[:, 2] = self.depth + a * sd
4381 return meta.DiscretizedPorePressureSource(
4382 times=util.num_full(n, self.time),
4383 lat=self.lat,
4384 lon=self.lon,
4385 north_shifts=points[:, 0],
4386 east_shifts=points[:, 1],
4387 depths=points[:, 2],
4388 pp=num.ones(n) / n)
4391class Request(Object):
4392 '''
4393 Synthetic seismogram computation request.
4395 ::
4397 Request(**kwargs)
4398 Request(sources, targets, **kwargs)
4399 '''
4401 sources = List.T(
4402 Source.T(),
4403 help='list of sources for which to produce synthetics.')
4405 targets = List.T(
4406 Target.T(),
4407 help='list of targets for which to produce synthetics.')
4409 @classmethod
4410 def args2kwargs(cls, args):
4411 if len(args) not in (0, 2, 3):
4412 raise BadRequest('Invalid arguments.')
4414 if len(args) == 2:
4415 return dict(sources=args[0], targets=args[1])
4416 else:
4417 return {}
4419 def __init__(self, *args, **kwargs):
4420 kwargs.update(self.args2kwargs(args))
4421 sources = kwargs.pop('sources', [])
4422 targets = kwargs.pop('targets', [])
4424 if isinstance(sources, Source):
4425 sources = [sources]
4427 if isinstance(targets, Target) or isinstance(targets, StaticTarget):
4428 targets = [targets]
4430 Object.__init__(self, sources=sources, targets=targets, **kwargs)
4432 @property
4433 def targets_dynamic(self):
4434 return [t for t in self.targets if isinstance(t, Target)]
4436 @property
4437 def targets_static(self):
4438 return [t for t in self.targets if isinstance(t, StaticTarget)]
4440 @property
4441 def has_dynamic(self):
4442 return True if len(self.targets_dynamic) > 0 else False
4444 @property
4445 def has_statics(self):
4446 return True if len(self.targets_static) > 0 else False
4448 def subsources_map(self):
4449 m = defaultdict(list)
4450 for source in self.sources:
4451 m[source.base_key()].append(source)
4453 return m
4455 def subtargets_map(self):
4456 m = defaultdict(list)
4457 for target in self.targets:
4458 m[target.base_key()].append(target)
4460 return m
4462 def subrequest_map(self):
4463 ms = self.subsources_map()
4464 mt = self.subtargets_map()
4465 m = {}
4466 for (ks, ls) in ms.items():
4467 for (kt, lt) in mt.items():
4468 m[ks, kt] = (ls, lt)
4470 return m
4473class ProcessingStats(Object):
4474 t_perc_get_store_and_receiver = Float.T(default=0.)
4475 t_perc_discretize_source = Float.T(default=0.)
4476 t_perc_make_base_seismogram = Float.T(default=0.)
4477 t_perc_make_same_span = Float.T(default=0.)
4478 t_perc_post_process = Float.T(default=0.)
4479 t_perc_optimize = Float.T(default=0.)
4480 t_perc_stack = Float.T(default=0.)
4481 t_perc_static_get_store = Float.T(default=0.)
4482 t_perc_static_discretize_basesource = Float.T(default=0.)
4483 t_perc_static_sum_statics = Float.T(default=0.)
4484 t_perc_static_post_process = Float.T(default=0.)
4485 t_wallclock = Float.T(default=0.)
4486 t_cpu = Float.T(default=0.)
4487 n_read_blocks = Int.T(default=0)
4488 n_results = Int.T(default=0)
4489 n_subrequests = Int.T(default=0)
4490 n_stores = Int.T(default=0)
4491 n_records_stacked = Int.T(default=0)
4494class Response(Object):
4495 '''
4496 Resonse object to a synthetic seismogram computation request.
4497 '''
4499 request = Request.T()
4500 results_list = List.T(List.T(meta.SeismosizerResult.T()))
4501 stats = ProcessingStats.T()
4503 def pyrocko_traces(self):
4504 '''
4505 Return a list of requested
4506 :class:`~pyrocko.trace.Trace` instances.
4507 '''
4509 traces = []
4510 for results in self.results_list:
4511 for result in results:
4512 if not isinstance(result, meta.Result):
4513 continue
4514 traces.append(result.trace.pyrocko_trace())
4516 return traces
4518 def kite_scenes(self):
4519 '''
4520 Return a list of requested
4521 :class:`~kite.scenes` instances.
4522 '''
4523 kite_scenes = []
4524 for results in self.results_list:
4525 for result in results:
4526 if isinstance(result, meta.KiteSceneResult):
4527 sc = result.get_scene()
4528 kite_scenes.append(sc)
4530 return kite_scenes
4532 def static_results(self):
4533 '''
4534 Return a list of requested
4535 :class:`~pyrocko.gf.meta.StaticResult` instances.
4536 '''
4537 statics = []
4538 for results in self.results_list:
4539 for result in results:
4540 if not isinstance(result, meta.StaticResult):
4541 continue
4542 statics.append(result)
4544 return statics
4546 def iter_results(self, get='pyrocko_traces'):
4547 '''
4548 Generator function to iterate over results of request.
4550 Yields associated :py:class:`Source`,
4551 :class:`~pyrocko.gf.targets.Target`,
4552 :class:`~pyrocko.trace.Trace` instances in each iteration.
4553 '''
4555 for isource, source in enumerate(self.request.sources):
4556 for itarget, target in enumerate(self.request.targets):
4557 result = self.results_list[isource][itarget]
4558 if get == 'pyrocko_traces':
4559 yield source, target, result.trace.pyrocko_trace()
4560 elif get == 'results':
4561 yield source, target, result
4563 def snuffle(self, **kwargs):
4564 '''
4565 Open *snuffler* with requested traces.
4566 '''
4568 trace.snuffle(self.pyrocko_traces(), **kwargs)
4571class Engine(Object):
4572 '''
4573 Base class for synthetic seismogram calculators.
4574 '''
4576 def get_store_ids(self):
4577 '''
4578 Get list of available GF store IDs
4579 '''
4581 return []
4584class Rule(object):
4585 pass
4588class VectorRule(Rule):
4590 def __init__(self, quantity, differentiate=0, integrate=0):
4591 self.components = [quantity + '.' + c for c in 'ned']
4592 self.differentiate = differentiate
4593 self.integrate = integrate
4595 def required_components(self, target):
4596 n, e, d = self.components
4597 sa, ca, sd, cd = target.get_sin_cos_factors()
4599 comps = []
4600 if nonzero(ca * cd):
4601 comps.append(n)
4603 if nonzero(sa * cd):
4604 comps.append(e)
4606 if nonzero(sd):
4607 comps.append(d)
4609 return tuple(comps)
4611 def apply_(self, target, base_seismogram):
4612 n, e, d = self.components
4613 sa, ca, sd, cd = target.get_sin_cos_factors()
4615 if nonzero(ca * cd):
4616 data = base_seismogram[n].data * (ca * cd)
4617 deltat = base_seismogram[n].deltat
4618 else:
4619 data = 0.0
4621 if nonzero(sa * cd):
4622 data = data + base_seismogram[e].data * (sa * cd)
4623 deltat = base_seismogram[e].deltat
4625 if nonzero(sd):
4626 data = data + base_seismogram[d].data * sd
4627 deltat = base_seismogram[d].deltat
4629 if self.differentiate:
4630 data = util.diff_fd(self.differentiate, 4, deltat, data)
4632 if self.integrate:
4633 raise NotImplementedError('Integration is not implemented yet.')
4635 return data
4638class HorizontalVectorRule(Rule):
4640 def __init__(self, quantity, differentiate=0, integrate=0):
4641 self.components = [quantity + '.' + c for c in 'ne']
4642 self.differentiate = differentiate
4643 self.integrate = integrate
4645 def required_components(self, target):
4646 n, e = self.components
4647 sa, ca, _, _ = target.get_sin_cos_factors()
4649 comps = []
4650 if nonzero(ca):
4651 comps.append(n)
4653 if nonzero(sa):
4654 comps.append(e)
4656 return tuple(comps)
4658 def apply_(self, target, base_seismogram):
4659 n, e = self.components
4660 sa, ca, _, _ = target.get_sin_cos_factors()
4662 if nonzero(ca):
4663 data = base_seismogram[n].data * ca
4664 else:
4665 data = 0.0
4667 if nonzero(sa):
4668 data = data + base_seismogram[e].data * sa
4670 if self.differentiate:
4671 deltat = base_seismogram[e].deltat
4672 data = util.diff_fd(self.differentiate, 4, deltat, data)
4674 if self.integrate:
4675 raise NotImplementedError('Integration is not implemented yet.')
4677 return data
4680class ScalarRule(Rule):
4682 def __init__(self, quantity, differentiate=0):
4683 self.c = quantity
4685 def required_components(self, target):
4686 return (self.c, )
4688 def apply_(self, target, base_seismogram):
4689 data = base_seismogram[self.c].data.copy()
4690 deltat = base_seismogram[self.c].deltat
4691 if self.differentiate:
4692 data = util.diff_fd(self.differentiate, 4, deltat, data)
4694 return data
4697class StaticDisplacement(Rule):
4699 def required_components(self, target):
4700 return tuple(['displacement.%s' % c for c in list('ned')])
4702 def apply_(self, target, base_statics):
4703 if isinstance(target, SatelliteTarget):
4704 los_fac = target.get_los_factors()
4705 base_statics['displacement.los'] =\
4706 (los_fac[:, 0] * -base_statics['displacement.d'] +
4707 los_fac[:, 1] * base_statics['displacement.e'] +
4708 los_fac[:, 2] * base_statics['displacement.n'])
4709 return base_statics
4712channel_rules = {
4713 'displacement': [VectorRule('displacement')],
4714 'rotation': [VectorRule('rotation')],
4715 'velocity': [
4716 VectorRule('velocity'),
4717 VectorRule('displacement', differentiate=1)],
4718 'acceleration': [
4719 VectorRule('acceleration'),
4720 VectorRule('velocity', differentiate=1),
4721 VectorRule('displacement', differentiate=2)],
4722 'pore_pressure': [ScalarRule('pore_pressure')],
4723 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')],
4724 'darcy_velocity': [VectorRule('darcy_velocity')],
4725}
4727static_rules = {
4728 'displacement': [StaticDisplacement()]
4729}
4732class OutOfBoundsContext(Object):
4733 source = Source.T()
4734 target = Target.T()
4735 distance = Float.T()
4736 components = List.T(String.T())
4739def process_dynamic_timeseries(work, psources, ptargets, engine, nthreads=0):
4740 dsource_cache = {}
4741 tcounters = list(range(6))
4743 store_ids = set()
4744 sources = set()
4745 targets = set()
4747 for itarget, target in enumerate(ptargets):
4748 target._id = itarget
4750 for w in work:
4751 _, _, isources, itargets = w
4753 sources.update([psources[isource] for isource in isources])
4754 targets.update([ptargets[itarget] for itarget in itargets])
4756 store_ids = set([t.store_id for t in targets])
4758 for isource, source in enumerate(psources):
4760 components = set()
4761 for itarget, target in enumerate(targets):
4762 rule = engine.get_rule(source, target)
4763 components.update(rule.required_components(target))
4765 for store_id in store_ids:
4766 store_targets = [t for t in targets if t.store_id == store_id]
4768 sample_rates = set([t.sample_rate for t in store_targets])
4769 interpolations = set([t.interpolation for t in store_targets])
4771 base_seismograms = []
4772 store_targets_out = []
4774 for samp_rate in sample_rates:
4775 for interp in interpolations:
4776 engine_targets = [
4777 t for t in store_targets if t.sample_rate == samp_rate
4778 and t.interpolation == interp]
4780 if not engine_targets:
4781 continue
4783 store_targets_out += engine_targets
4785 base_seismograms += engine.base_seismograms(
4786 source,
4787 engine_targets,
4788 components,
4789 dsource_cache,
4790 nthreads)
4792 for iseis, seismogram in enumerate(base_seismograms):
4793 for tr in seismogram.values():
4794 if tr.err != store.SeismosizerErrorEnum.SUCCESS:
4795 e = SeismosizerError(
4796 'Seismosizer failed with return code %i\n%s' % (
4797 tr.err, str(
4798 OutOfBoundsContext(
4799 source=source,
4800 target=store_targets[iseis],
4801 distance=source.distance_to(
4802 store_targets[iseis]),
4803 components=components))))
4804 raise e
4806 for seismogram, target in zip(base_seismograms, store_targets_out):
4808 try:
4809 result = engine._post_process_dynamic(
4810 seismogram, source, target)
4811 except SeismosizerError as e:
4812 result = e
4814 yield (isource, target._id, result), tcounters
4817def process_dynamic(work, psources, ptargets, engine, nthreads=0):
4818 dsource_cache = {}
4820 for w in work:
4821 _, _, isources, itargets = w
4823 sources = [psources[isource] for isource in isources]
4824 targets = [ptargets[itarget] for itarget in itargets]
4826 components = set()
4827 for target in targets:
4828 rule = engine.get_rule(sources[0], target)
4829 components.update(rule.required_components(target))
4831 for isource, source in zip(isources, sources):
4832 for itarget, target in zip(itargets, targets):
4834 try:
4835 base_seismogram, tcounters = engine.base_seismogram(
4836 source, target, components, dsource_cache, nthreads)
4837 except meta.OutOfBounds as e:
4838 e.context = OutOfBoundsContext(
4839 source=sources[0],
4840 target=targets[0],
4841 distance=sources[0].distance_to(targets[0]),
4842 components=components)
4843 raise
4845 n_records_stacked = 0
4846 t_optimize = 0.0
4847 t_stack = 0.0
4849 for _, tr in base_seismogram.items():
4850 n_records_stacked += tr.n_records_stacked
4851 t_optimize += tr.t_optimize
4852 t_stack += tr.t_stack
4854 try:
4855 result = engine._post_process_dynamic(
4856 base_seismogram, source, target)
4857 result.n_records_stacked = n_records_stacked
4858 result.n_shared_stacking = len(sources) *\
4859 len(targets)
4860 result.t_optimize = t_optimize
4861 result.t_stack = t_stack
4862 except SeismosizerError as e:
4863 result = e
4865 tcounters.append(xtime())
4866 yield (isource, itarget, result), tcounters
4869def process_static(work, psources, ptargets, engine, nthreads=0):
4870 for w in work:
4871 _, _, isources, itargets = w
4873 sources = [psources[isource] for isource in isources]
4874 targets = [ptargets[itarget] for itarget in itargets]
4876 for isource, source in zip(isources, sources):
4877 for itarget, target in zip(itargets, targets):
4878 components = engine.get_rule(source, target)\
4879 .required_components(target)
4881 try:
4882 base_statics, tcounters = engine.base_statics(
4883 source, target, components, nthreads)
4884 except meta.OutOfBounds as e:
4885 e.context = OutOfBoundsContext(
4886 source=sources[0],
4887 target=targets[0],
4888 distance=float('nan'),
4889 components=components)
4890 raise
4891 result = engine._post_process_statics(
4892 base_statics, source, target)
4893 tcounters.append(xtime())
4895 yield (isource, itarget, result), tcounters
4898class LocalEngine(Engine):
4899 '''
4900 Offline synthetic seismogram calculator.
4902 :param use_env: if ``True``, fill :py:attr:`store_superdirs` and
4903 :py:attr:`store_dirs` with paths set in environment variables
4904 GF_STORE_SUPERDIRS and GF_STORE_DIRS.
4905 :param use_config: if ``True``, fill :py:attr:`store_superdirs` and
4906 :py:attr:`store_dirs` with paths set in the user's config file.
4908 The config file can be found at :file:`~/.pyrocko/config.pf`
4910 .. code-block :: python
4912 gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/']
4913 gf_store_superdirs: ['/home/pyrocko/gf_stores/']
4914 '''
4916 store_superdirs = List.T(
4917 String.T(),
4918 help='directories which are searched for Green\'s function stores')
4920 store_dirs = List.T(
4921 String.T(),
4922 help='additional individual Green\'s function store directories')
4924 default_store_id = String.T(
4925 optional=True,
4926 help='default store ID to be used when a request does not provide '
4927 'one')
4929 def __init__(self, **kwargs):
4930 use_env = kwargs.pop('use_env', False)
4931 use_config = kwargs.pop('use_config', False)
4932 Engine.__init__(self, **kwargs)
4933 if use_env:
4934 env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '')
4935 env_store_dirs = os.environ.get('GF_STORE_DIRS', '')
4936 if env_store_superdirs:
4937 self.store_superdirs.extend(env_store_superdirs.split(':'))
4939 if env_store_dirs:
4940 self.store_dirs.extend(env_store_dirs.split(':'))
4942 if use_config:
4943 c = config.config()
4944 self.store_superdirs.extend(c.gf_store_superdirs)
4945 self.store_dirs.extend(c.gf_store_dirs)
4947 self._check_store_dirs_type()
4948 self._id_to_store_dir = {}
4949 self._open_stores = {}
4950 self._effective_default_store_id = None
4952 def _check_store_dirs_type(self):
4953 for sdir in ['store_dirs', 'store_superdirs']:
4954 if not isinstance(self.__getattribute__(sdir), list):
4955 raise TypeError("{} of {} is not of type list".format(
4956 sdir, self.__class__.__name__))
4958 def _get_store_id(self, store_dir):
4959 store_ = store.Store(store_dir)
4960 store_id = store_.config.id
4961 store_.close()
4962 return store_id
4964 def _looks_like_store_dir(self, store_dir):
4965 return os.path.isdir(store_dir) and \
4966 all(os.path.isfile(pjoin(store_dir, x)) for x in
4967 ('index', 'traces', 'config'))
4969 def iter_store_dirs(self):
4970 store_dirs = set()
4971 for d in self.store_superdirs:
4972 if not os.path.exists(d):
4973 logger.warning('store_superdir not available: %s' % d)
4974 continue
4976 for entry in os.listdir(d):
4977 store_dir = os.path.realpath(pjoin(d, entry))
4978 if self._looks_like_store_dir(store_dir):
4979 store_dirs.add(store_dir)
4981 for store_dir in self.store_dirs:
4982 store_dirs.add(os.path.realpath(store_dir))
4984 return store_dirs
4986 def _scan_stores(self):
4987 for store_dir in self.iter_store_dirs():
4988 store_id = self._get_store_id(store_dir)
4989 if store_id not in self._id_to_store_dir:
4990 self._id_to_store_dir[store_id] = store_dir
4991 else:
4992 if store_dir != self._id_to_store_dir[store_id]:
4993 raise DuplicateStoreId(
4994 'GF store ID %s is used in (at least) two '
4995 'different stores. Locations are: %s and %s' %
4996 (store_id, self._id_to_store_dir[store_id], store_dir))
4998 def get_store_dir(self, store_id):
4999 '''
5000 Lookup directory given a GF store ID.
5001 '''
5003 if store_id not in self._id_to_store_dir:
5004 self._scan_stores()
5006 if store_id not in self._id_to_store_dir:
5007 raise NoSuchStore(store_id, self.iter_store_dirs())
5009 return self._id_to_store_dir[store_id]
5011 def get_store_ids(self):
5012 '''
5013 Get list of available store IDs.
5014 '''
5016 self._scan_stores()
5017 return sorted(self._id_to_store_dir.keys())
5019 def effective_default_store_id(self):
5020 if self._effective_default_store_id is None:
5021 if self.default_store_id is None:
5022 store_ids = self.get_store_ids()
5023 if len(store_ids) == 1:
5024 self._effective_default_store_id = self.get_store_ids()[0]
5025 else:
5026 raise NoDefaultStoreSet()
5027 else:
5028 self._effective_default_store_id = self.default_store_id
5030 return self._effective_default_store_id
5032 def get_store(self, store_id=None):
5033 '''
5034 Get a store from the engine.
5036 :param store_id: identifier of the store (optional)
5037 :returns: :py:class:`~pyrocko.gf.store.Store` object
5039 If no ``store_id`` is provided the store
5040 associated with the :py:gattr:`default_store_id` is returned.
5041 Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is
5042 undefined.
5043 '''
5045 if store_id is None:
5046 store_id = self.effective_default_store_id()
5048 if store_id not in self._open_stores:
5049 store_dir = self.get_store_dir(store_id)
5050 self._open_stores[store_id] = store.Store(store_dir)
5052 return self._open_stores[store_id]
5054 def get_store_config(self, store_id):
5055 store = self.get_store(store_id)
5056 return store.config
5058 def get_store_extra(self, store_id, key):
5059 store = self.get_store(store_id)
5060 return store.get_extra(key)
5062 def close_cashed_stores(self):
5063 '''
5064 Close and remove ids from cashed stores.
5065 '''
5066 store_ids = []
5067 for store_id, store_ in self._open_stores.items():
5068 store_.close()
5069 store_ids.append(store_id)
5071 for store_id in store_ids:
5072 self._open_stores.pop(store_id)
5074 def get_rule(self, source, target):
5075 cprovided = self.get_store(target.store_id).get_provided_components()
5077 if isinstance(target, StaticTarget):
5078 quantity = target.quantity
5079 available_rules = static_rules
5080 elif isinstance(target, Target):
5081 quantity = target.effective_quantity()
5082 available_rules = channel_rules
5084 try:
5085 for rule in available_rules[quantity]:
5086 cneeded = rule.required_components(target)
5087 if all(c in cprovided for c in cneeded):
5088 return rule
5090 except KeyError:
5091 pass
5093 raise BadRequest(
5094 'No rule to calculate "%s" with GFs from store "%s" '
5095 'for source model "%s".' % (
5096 target.effective_quantity(),
5097 target.store_id,
5098 source.__class__.__name__))
5100 def _cached_discretize_basesource(self, source, store, cache, target):
5101 if (source, store) not in cache:
5102 cache[source, store] = source.discretize_basesource(store, target)
5104 return cache[source, store]
5106 def base_seismograms(self, source, targets, components, dsource_cache,
5107 nthreads=0):
5109 target = targets[0]
5111 interp = set([t.interpolation for t in targets])
5112 if len(interp) > 1:
5113 raise BadRequest('Targets have different interpolation schemes.')
5115 rates = set([t.sample_rate for t in targets])
5116 if len(rates) > 1:
5117 raise BadRequest('Targets have different sample rates.')
5119 store_ = self.get_store(target.store_id)
5120 receivers = [t.receiver(store_) for t in targets]
5122 if target.sample_rate is not None:
5123 deltat = 1. / target.sample_rate
5124 rate = target.sample_rate
5125 else:
5126 deltat = None
5127 rate = store_.config.sample_rate
5129 tmin = num.fromiter(
5130 (t.tmin for t in targets), dtype=float, count=len(targets))
5131 tmax = num.fromiter(
5132 (t.tmax for t in targets), dtype=float, count=len(targets))
5134 itmin = num.floor(tmin * rate).astype(num.int64)
5135 itmax = num.ceil(tmax * rate).astype(num.int64)
5136 nsamples = itmax - itmin + 1
5138 mask = num.isnan(tmin)
5139 itmin[mask] = 0
5140 nsamples[mask] = -1
5142 base_source = self._cached_discretize_basesource(
5143 source, store_, dsource_cache, target)
5145 base_seismograms = store_.calc_seismograms(
5146 base_source, receivers, components,
5147 deltat=deltat,
5148 itmin=itmin, nsamples=nsamples,
5149 interpolation=target.interpolation,
5150 optimization=target.optimization,
5151 nthreads=nthreads)
5153 for i, base_seismogram in enumerate(base_seismograms):
5154 base_seismograms[i] = store.make_same_span(base_seismogram)
5156 return base_seismograms
5158 def base_seismogram(self, source, target, components, dsource_cache,
5159 nthreads):
5161 tcounters = [xtime()]
5163 store_ = self.get_store(target.store_id)
5164 receiver = target.receiver(store_)
5166 if target.tmin and target.tmax is not None:
5167 rate = store_.config.sample_rate
5168 itmin = int(num.floor(target.tmin * rate))
5169 itmax = int(num.ceil(target.tmax * rate))
5170 nsamples = itmax - itmin + 1
5171 else:
5172 itmin = None
5173 nsamples = None
5175 tcounters.append(xtime())
5176 base_source = self._cached_discretize_basesource(
5177 source, store_, dsource_cache, target)
5179 tcounters.append(xtime())
5181 if target.sample_rate is not None:
5182 deltat = 1. / target.sample_rate
5183 else:
5184 deltat = None
5186 base_seismogram = store_.seismogram(
5187 base_source, receiver, components,
5188 deltat=deltat,
5189 itmin=itmin, nsamples=nsamples,
5190 interpolation=target.interpolation,
5191 optimization=target.optimization,
5192 nthreads=nthreads)
5194 tcounters.append(xtime())
5196 base_seismogram = store.make_same_span(base_seismogram)
5198 tcounters.append(xtime())
5200 return base_seismogram, tcounters
5202 def base_statics(self, source, target, components, nthreads):
5203 tcounters = [xtime()]
5204 store_ = self.get_store(target.store_id)
5206 if target.tsnapshot is not None:
5207 rate = store_.config.sample_rate
5208 itsnapshot = int(num.floor(target.tsnapshot * rate))
5209 else:
5210 itsnapshot = None
5211 tcounters.append(xtime())
5213 base_source = source.discretize_basesource(store_, target=target)
5215 tcounters.append(xtime())
5217 base_statics = store_.statics(
5218 base_source,
5219 target,
5220 itsnapshot,
5221 components,
5222 target.interpolation,
5223 nthreads)
5225 tcounters.append(xtime())
5227 return base_statics, tcounters
5229 def _post_process_dynamic(self, base_seismogram, source, target):
5230 base_any = next(iter(base_seismogram.values()))
5231 deltat = base_any.deltat
5232 itmin = base_any.itmin
5234 rule = self.get_rule(source, target)
5235 data = rule.apply_(target, base_seismogram)
5237 factor = source.get_factor() * target.get_factor()
5238 if factor != 1.0:
5239 data = data * factor
5241 stf = source.effective_stf_post()
5243 times, amplitudes = stf.discretize_t(
5244 deltat, 0.0)
5246 # repeat end point to prevent boundary effects
5247 padded_data = num.empty(data.size + amplitudes.size, dtype=float)
5248 padded_data[:data.size] = data
5249 padded_data[data.size:] = data[-1]
5250 data = num.convolve(amplitudes, padded_data)
5252 tmin = itmin * deltat + times[0]
5254 tr = meta.SeismosizerTrace(
5255 codes=target.codes,
5256 data=data[:-amplitudes.size],
5257 deltat=deltat,
5258 tmin=tmin)
5260 return target.post_process(self, source, tr)
5262 def _post_process_statics(self, base_statics, source, starget):
5263 rule = self.get_rule(source, starget)
5264 data = rule.apply_(starget, base_statics)
5266 factor = source.get_factor()
5267 if factor != 1.0:
5268 for v in data.values():
5269 v *= factor
5271 return starget.post_process(self, source, base_statics)
5273 def process(self, *args, **kwargs):
5274 '''
5275 Process a request.
5277 ::
5279 process(**kwargs)
5280 process(request, **kwargs)
5281 process(sources, targets, **kwargs)
5283 The request can be given a a :py:class:`Request` object, or such an
5284 object is created using ``Request(**kwargs)`` for convenience.
5286 :returns: :py:class:`Response` object
5287 '''
5289 if len(args) not in (0, 1, 2):
5290 raise BadRequest('Invalid arguments.')
5292 if len(args) == 1:
5293 kwargs['request'] = args[0]
5295 elif len(args) == 2:
5296 kwargs.update(Request.args2kwargs(args))
5298 request = kwargs.pop('request', None)
5299 status_callback = kwargs.pop('status_callback', None)
5300 calc_timeseries = kwargs.pop('calc_timeseries', True)
5302 nprocs = kwargs.pop('nprocs', None)
5303 nthreads = kwargs.pop('nthreads', 1)
5304 if nprocs is not None:
5305 nthreads = nprocs
5307 if request is None:
5308 request = Request(**kwargs)
5310 if resource:
5311 rs0 = resource.getrusage(resource.RUSAGE_SELF)
5312 rc0 = resource.getrusage(resource.RUSAGE_CHILDREN)
5313 tt0 = xtime()
5315 # make sure stores are open before fork()
5316 store_ids = set(target.store_id for target in request.targets)
5317 for store_id in store_ids:
5318 self.get_store(store_id)
5320 source_index = dict((x, i) for (i, x) in
5321 enumerate(request.sources))
5322 target_index = dict((x, i) for (i, x) in
5323 enumerate(request.targets))
5325 m = request.subrequest_map()
5327 skeys = sorted(m.keys(), key=cmp_to_key(cmp_none_aware))
5328 results_list = []
5330 for i in range(len(request.sources)):
5331 results_list.append([None] * len(request.targets))
5333 tcounters_dyn_list = []
5334 tcounters_static_list = []
5335 nsub = len(skeys)
5336 isub = 0
5338 # Processing dynamic targets through
5339 # parimap(process_subrequest_dynamic)
5341 if calc_timeseries:
5342 _process_dynamic = process_dynamic_timeseries
5343 else:
5344 _process_dynamic = process_dynamic
5346 if request.has_dynamic:
5347 work_dynamic = [
5348 (i, nsub,
5349 [source_index[source] for source in m[k][0]],
5350 [target_index[target] for target in m[k][1]
5351 if not isinstance(target, StaticTarget)])
5352 for (i, k) in enumerate(skeys)]
5354 for ii_results, tcounters_dyn in _process_dynamic(
5355 work_dynamic, request.sources, request.targets, self,
5356 nthreads):
5358 tcounters_dyn_list.append(num.diff(tcounters_dyn))
5359 isource, itarget, result = ii_results
5360 results_list[isource][itarget] = result
5362 if status_callback:
5363 status_callback(isub, nsub)
5365 isub += 1
5367 # Processing static targets through process_static
5368 if request.has_statics:
5369 work_static = [
5370 (i, nsub,
5371 [source_index[source] for source in m[k][0]],
5372 [target_index[target] for target in m[k][1]
5373 if isinstance(target, StaticTarget)])
5374 for (i, k) in enumerate(skeys)]
5376 for ii_results, tcounters_static in process_static(
5377 work_static, request.sources, request.targets, self,
5378 nthreads=nthreads):
5380 tcounters_static_list.append(num.diff(tcounters_static))
5381 isource, itarget, result = ii_results
5382 results_list[isource][itarget] = result
5384 if status_callback:
5385 status_callback(isub, nsub)
5387 isub += 1
5389 if status_callback:
5390 status_callback(nsub, nsub)
5392 tt1 = time.time()
5393 if resource:
5394 rs1 = resource.getrusage(resource.RUSAGE_SELF)
5395 rc1 = resource.getrusage(resource.RUSAGE_CHILDREN)
5397 s = ProcessingStats()
5399 if request.has_dynamic:
5400 tcumu_dyn = num.sum(num.vstack(tcounters_dyn_list), axis=0)
5401 t_dyn = float(num.sum(tcumu_dyn))
5402 perc_dyn = map(float, tcumu_dyn / t_dyn * 100.)
5403 (s.t_perc_get_store_and_receiver,
5404 s.t_perc_discretize_source,
5405 s.t_perc_make_base_seismogram,
5406 s.t_perc_make_same_span,
5407 s.t_perc_post_process) = perc_dyn
5408 else:
5409 t_dyn = 0.
5411 if request.has_statics:
5412 tcumu_static = num.sum(num.vstack(tcounters_static_list), axis=0)
5413 t_static = num.sum(tcumu_static)
5414 perc_static = map(float, tcumu_static / t_static * 100.)
5415 (s.t_perc_static_get_store,
5416 s.t_perc_static_discretize_basesource,
5417 s.t_perc_static_sum_statics,
5418 s.t_perc_static_post_process) = perc_static
5420 s.t_wallclock = tt1 - tt0
5421 if resource:
5422 s.t_cpu = (
5423 (rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) -
5424 (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime))
5425 s.n_read_blocks = (
5426 (rs1.ru_inblock + rc1.ru_inblock) -
5427 (rs0.ru_inblock + rc0.ru_inblock))
5429 n_records_stacked = 0.
5430 for results in results_list:
5431 for result in results:
5432 if not isinstance(result, meta.Result):
5433 continue
5434 shr = float(result.n_shared_stacking)
5435 n_records_stacked += result.n_records_stacked / shr
5436 s.t_perc_optimize += result.t_optimize / shr
5437 s.t_perc_stack += result.t_stack / shr
5438 s.n_records_stacked = int(n_records_stacked)
5439 if t_dyn != 0.:
5440 s.t_perc_optimize /= t_dyn * 100
5441 s.t_perc_stack /= t_dyn * 100
5443 return Response(
5444 request=request,
5445 results_list=results_list,
5446 stats=s)
5449class RemoteEngine(Engine):
5450 '''
5451 Client for remote synthetic seismogram calculator.
5452 '''
5454 site = String.T(default=ws.g_default_site, optional=True)
5455 url = String.T(default=ws.g_url, optional=True)
5457 def process(self, request=None, status_callback=None, **kwargs):
5459 if request is None:
5460 request = Request(**kwargs)
5462 return ws.seismosizer(url=self.url, site=self.site, request=request)
5465g_engine = None
5468def get_engine(store_superdirs=[]):
5469 global g_engine
5470 if g_engine is None:
5471 g_engine = LocalEngine(use_env=True, use_config=True)
5473 for d in store_superdirs:
5474 if d not in g_engine.store_superdirs:
5475 g_engine.store_superdirs.append(d)
5477 return g_engine
5480class SourceGroup(Object):
5482 def __getattr__(self, k):
5483 return num.fromiter((getattr(s, k) for s in self),
5484 dtype=float)
5486 def __iter__(self):
5487 raise NotImplementedError(
5488 'This method should be implemented in subclass.')
5490 def __len__(self):
5491 raise NotImplementedError(
5492 'This method should be implemented in subclass.')
5495class SourceList(SourceGroup):
5496 sources = List.T(Source.T())
5498 def append(self, s):
5499 self.sources.append(s)
5501 def __iter__(self):
5502 return iter(self.sources)
5504 def __len__(self):
5505 return len(self.sources)
5508class SourceGrid(SourceGroup):
5510 base = Source.T()
5511 variables = Dict.T(String.T(), Range.T())
5512 order = List.T(String.T())
5514 def __len__(self):
5515 n = 1
5516 for (k, v) in self.make_coords(self.base):
5517 n *= len(list(v))
5519 return n
5521 def __iter__(self):
5522 for items in permudef(self.make_coords(self.base)):
5523 s = self.base.clone(**{k: v for (k, v) in items})
5524 s.regularize()
5525 yield s
5527 def ordered_params(self):
5528 ks = list(self.variables.keys())
5529 for k in self.order + list(self.base.keys()):
5530 if k in ks:
5531 yield k
5532 ks.remove(k)
5533 if ks:
5534 raise Exception('Invalid parameter "%s" for source type "%s".' %
5535 (ks[0], self.base.__class__.__name__))
5537 def make_coords(self, base):
5538 return [(param, self.variables[param].make(base=base[param]))
5539 for param in self.ordered_params()]
5542source_classes = [
5543 Source,
5544 SourceWithMagnitude,
5545 SourceWithDerivedMagnitude,
5546 ExplosionSource,
5547 RectangularExplosionSource,
5548 DCSource,
5549 CLVDSource,
5550 VLVDSource,
5551 MTSource,
5552 RectangularSource,
5553 PseudoDynamicRupture,
5554 DoubleDCSource,
5555 RingfaultSource,
5556 CombiSource,
5557 SFSource,
5558 PorePressurePointSource,
5559 PorePressureLineSource,
5560]
5562stf_classes = [
5563 STF,
5564 BoxcarSTF,
5565 TriangularSTF,
5566 HalfSinusoidSTF,
5567 ResonatorSTF,
5568]
5570__all__ = '''
5571SeismosizerError
5572BadRequest
5573NoSuchStore
5574DerivedMagnitudeError
5575STFMode
5576'''.split() + [S.__name__ for S in source_classes + stf_classes] + '''
5577Request
5578ProcessingStats
5579Response
5580Engine
5581LocalEngine
5582RemoteEngine
5583source_classes
5584get_engine
5585Range
5586SourceGroup
5587SourceList
5588SourceGrid
5589map_anchor
5590'''.split()