Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/seismosizer.py: 83%
2470 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-12 12:00 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7High level synthetic seismogram synthesis.
9.. _coordinate-system-names:
11Coordinate systems
12..................
14Coordinate system names commonly used in source models.
16================= ============================================
17Name Description
18================= ============================================
19``'xyz'`` northing, easting, depth in [m]
20``'xy'`` northing, easting in [m]
21``'latlon'`` latitude, longitude in [deg]
22``'lonlat'`` longitude, latitude in [deg]
23``'latlondepth'`` latitude, longitude in [deg], depth in [m]
24================= ============================================
25'''
27from collections import defaultdict
28from functools import cmp_to_key
29import time
30import math
31import os
32import re
33import logging
34try:
35 import resource
36except ImportError:
37 resource = None
38from hashlib import sha1
40import numpy as num
41from scipy.interpolate import RegularGridInterpolator
43from pyrocko.guts import (Object, Float, String, StringChoice, List,
44 Timestamp, Int, SObject, ArgumentError, Dict,
45 ValidationError, Bool)
46from pyrocko.guts_array import Array
48from pyrocko import moment_tensor as pmt
49from pyrocko import trace, util, config, model, eikonal_ext
50from pyrocko.orthodrome import ne_to_latlon
51from pyrocko.model import Location
52from pyrocko.modelling import OkadaSource, make_okada_coefficient_matrix, \
53 okada_ext, invert_fault_dislocations_bem
55from . import meta, store, ws
56from .tractions import TractionField, DirectedTractions
57from .targets import Target, StaticTarget, SatelliteTarget
59pjoin = os.path.join
61guts_prefix = 'pf'
63d2r = math.pi / 180.
64r2d = 180. / math.pi
65km = 1e3
67logger = logging.getLogger('pyrocko.gf.seismosizer')
70def cmp_none_aware(a, b):
71 if isinstance(a, tuple) and isinstance(b, tuple):
72 for xa, xb in zip(a, b):
73 rv = cmp_none_aware(xa, xb)
74 if rv != 0:
75 return rv
77 return 0
79 anone = a is None
80 bnone = b is None
82 if anone and bnone:
83 return 0
85 if anone:
86 return -1
88 if bnone:
89 return 1
91 return bool(a > b) - bool(a < b)
94def xtime():
95 return time.time()
98class SeismosizerError(Exception):
99 pass
102class BadRequest(SeismosizerError):
103 pass
106class DuplicateStoreId(Exception):
107 pass
110class NoDefaultStoreSet(Exception):
111 '''
112 Raised, when a default store would be used but none is set.
113 '''
114 pass
117class ConversionError(Exception):
118 pass
121class NoSuchStore(BadRequest):
123 def __init__(self, store_id=None, dirs=None):
124 BadRequest.__init__(self)
125 self.store_id = store_id
126 self.dirs = dirs
128 def __str__(self):
129 if self.store_id is not None:
130 rstr = 'no GF store with id "%s" found.' % self.store_id
131 else:
132 rstr = 'GF store not found.'
134 if self.dirs is not None:
135 rstr += ' Searched folders:\n %s' % '\n '.join(sorted(self.dirs))
136 return rstr
139def ufloat(s):
140 units = {
141 'k': 1e3,
142 'M': 1e6,
143 }
145 factor = 1.0
146 if s and s[-1] in units:
147 factor = units[s[-1]]
148 s = s[:-1]
149 if not s:
150 raise ValueError("unit without a number: '%s'" % s)
152 return float(s) * factor
155def ufloat_or_none(s):
156 if s:
157 return ufloat(s)
158 else:
159 return None
162def int_or_none(s):
163 if s:
164 return int(s)
165 else:
166 return None
169def nonzero(x, eps=1e-15):
170 return abs(x) > eps
173def permudef(ln, j=0):
174 if j < len(ln):
175 k, v = ln[j]
176 for y in v:
177 ln[j] = k, y
178 for s in permudef(ln, j + 1):
179 yield s
181 ln[j] = k, v
182 return
183 else:
184 yield ln
187def arr(x):
188 return num.atleast_1d(num.asarray(x))
191def discretize_rect_source(deltas, deltat, time, north, east, depth,
192 strike, dip, length, width,
193 anchor, velocity=None, stf=None,
194 nucleation_x=None, nucleation_y=None,
195 decimation_factor=1, pointsonly=False,
196 plane_coords=False,
197 aggressive_oversampling=False):
199 if stf is None:
200 stf = STF()
202 if not velocity and not pointsonly:
203 raise AttributeError('velocity is required in time mode')
205 mindeltagf = float(num.min(deltas))
206 if velocity:
207 mindeltagf = min(mindeltagf, deltat * velocity)
209 ln = length
210 wd = width
212 if aggressive_oversampling:
213 nl = int((2. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
214 nw = int((2. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
215 else:
216 nl = int((1. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
217 nw = int((1. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
219 n = int(nl * nw)
221 dl = ln / nl
222 dw = wd / nw
224 xl = num.linspace(-0.5 * (ln - dl), 0.5 * (ln - dl), nl)
225 xw = num.linspace(-0.5 * (wd - dw), 0.5 * (wd - dw), nw)
227 points = num.zeros((n, 3))
228 points[:, 0] = num.tile(xl, nw)
229 points[:, 1] = num.repeat(xw, nl)
231 if nucleation_x is not None:
232 dist_x = num.abs(nucleation_x - points[:, 0])
233 else:
234 dist_x = num.zeros(n)
236 if nucleation_y is not None:
237 dist_y = num.abs(nucleation_y - points[:, 1])
238 else:
239 dist_y = num.zeros(n)
241 dist = num.sqrt(dist_x**2 + dist_y**2)
242 times = dist / velocity
244 anch_x, anch_y = map_anchor[anchor]
246 points[:, 0] -= anch_x * 0.5 * length
247 points[:, 1] -= anch_y * 0.5 * width
249 if plane_coords:
250 return points, dl, dw, nl, nw
252 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
253 points = num.dot(rotmat.T, points.T).T
255 points[:, 0] += north
256 points[:, 1] += east
257 points[:, 2] += depth
259 if pointsonly:
260 return points, dl, dw, nl, nw
262 xtau, amplitudes = stf.discretize_t(deltat, time)
263 nt = xtau.size
265 points2 = num.repeat(points, nt, axis=0)
266 times2 = (times[:, num.newaxis] + xtau[num.newaxis, :]).ravel()
267 amplitudes2 = num.tile(amplitudes, n)
269 return points2, times2, amplitudes2, dl, dw, nl, nw
272def check_rect_source_discretisation(points2, nl, nw, store):
273 # We assume a non-rotated fault plane
274 N_CRITICAL = 8
275 points = points2.T.reshape((3, nl, nw))
276 if points.size <= N_CRITICAL:
277 logger.warning('RectangularSource is defined by only %d sub-sources!'
278 % points.size)
279 return True
281 distances = num.sqrt(
282 (points[0, 0, :] - points[0, 1, :])**2 +
283 (points[1, 0, :] - points[1, 1, :])**2 +
284 (points[2, 0, :] - points[2, 1, :])**2)
286 depths = points[2, 0, :]
287 vs_profile = store.config.get_vs(
288 lat=0., lon=0.,
289 points=num.repeat(depths[:, num.newaxis], 3, axis=1),
290 interpolation='multilinear')
292 min_wavelength = vs_profile * (store.config.deltat * 2)
293 if not num.all(min_wavelength > distances / 2):
294 return False
295 return True
298def outline_rect_source(strike, dip, length, width, anchor):
299 ln = length
300 wd = width
301 points = num.array(
302 [[-0.5 * ln, -0.5 * wd, 0.],
303 [0.5 * ln, -0.5 * wd, 0.],
304 [0.5 * ln, 0.5 * wd, 0.],
305 [-0.5 * ln, 0.5 * wd, 0.],
306 [-0.5 * ln, -0.5 * wd, 0.]])
308 anch_x, anch_y = map_anchor[anchor]
309 points[:, 0] -= anch_x * 0.5 * length
310 points[:, 1] -= anch_y * 0.5 * width
312 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
314 return num.dot(rotmat.T, points.T).T
317def from_plane_coords(
318 strike, dip, length, width, depth, x_plane_coords, y_plane_coords,
319 lat=0., lon=0.,
320 north_shift=0, east_shift=0,
321 anchor='top', cs='xy'):
323 ln = length
324 wd = width
325 x_abs = []
326 y_abs = []
327 if not isinstance(x_plane_coords, list):
328 x_plane_coords = [x_plane_coords]
329 y_plane_coords = [y_plane_coords]
331 for x_plane, y_plane in zip(x_plane_coords, y_plane_coords):
332 points = num.array(
333 [[-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
334 [0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
335 [0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
336 [-0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
337 [-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.]])
339 anch_x, anch_y = map_anchor[anchor]
340 points[:, 0] -= anch_x * 0.5 * length
341 points[:, 1] -= anch_y * 0.5 * width
343 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
345 points = num.dot(rotmat.T, points.T).T
346 points[:, 0] += north_shift
347 points[:, 1] += east_shift
348 points[:, 2] += depth
349 if cs in ('latlon', 'lonlat'):
350 latlon = ne_to_latlon(lat, lon,
351 points[:, 0], points[:, 1])
352 latlon = num.array(latlon).T
353 x_abs.append(latlon[1:2, 1])
354 y_abs.append(latlon[2:3, 0])
355 if cs == 'xy':
356 x_abs.append(points[1:2, 1])
357 y_abs.append(points[2:3, 0])
359 if cs == 'lonlat':
360 return y_abs, x_abs
361 else:
362 return x_abs, y_abs
365def points_on_rect_source(
366 strike, dip, length, width, anchor,
367 discretized_basesource=None, points_x=None, points_y=None):
369 ln = length
370 wd = width
372 if isinstance(points_x, list) or isinstance(points_x, float):
373 points_x = num.array([points_x])
374 if isinstance(points_y, list) or isinstance(points_y, float):
375 points_y = num.array([points_y])
377 if discretized_basesource:
378 ds = discretized_basesource
380 nl_patches = ds.nl + 1
381 nw_patches = ds.nw + 1
383 npoints = nl_patches * nw_patches
384 points = num.zeros((npoints, 3))
385 ln_patches = num.array([il for il in range(nl_patches)])
386 wd_patches = num.array([iw for iw in range(nw_patches)])
388 points_ln =\
389 2 * ((ln_patches - num.min(ln_patches)) / num.ptp(ln_patches)) - 1
390 points_wd =\
391 2 * ((wd_patches - num.min(wd_patches)) / num.ptp(wd_patches)) - 1
393 for il in range(nl_patches):
394 for iw in range(nw_patches):
395 points[il * nw_patches + iw, :] = num.array([
396 points_ln[il] * ln * 0.5,
397 points_wd[iw] * wd * 0.5, 0.0])
399 elif points_x.shape[0] > 0 and points_y.shape[0] > 0:
400 points = num.zeros(shape=((len(points_x), 3)))
401 for i, (x, y) in enumerate(zip(points_x, points_y)):
402 points[i, :] = num.array(
403 [x * 0.5 * ln, y * 0.5 * wd, 0.0])
405 anch_x, anch_y = map_anchor[anchor]
407 points[:, 0] -= anch_x * 0.5 * ln
408 points[:, 1] -= anch_y * 0.5 * wd
410 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
412 return num.dot(rotmat.T, points.T).T
415class InvalidGridDef(Exception):
416 pass
419class Range(SObject):
420 '''
421 Convenient range specification.
423 Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]::
425 Range('0 .. 10k : 1k')
426 Range(start=0., stop=10e3, step=1e3)
427 Range(0, 10e3, 1e3)
428 Range('0 .. 10k @ 11')
429 Range(start=0., stop=10*km, n=11)
431 Range(0, 10e3, n=11)
432 Range(values=[x*1e3 for x in range(11)])
434 Depending on the use context, it can be possible to omit any part of the
435 specification. E.g. in the context of extracting a subset of an already
436 existing range, the existing range's specification values would be filled
437 in where missing.
439 The values are distributed with equal spacing, unless the ``spacing``
440 argument is modified. The values can be created offset or relative to an
441 external base value with the ``relative`` argument if the use context
442 supports this.
444 The range specification can be expressed with a short string
445 representation::
447 'start .. stop @ num | spacing, relative'
448 'start .. stop : step | spacing, relative'
450 most parts of the expression can be omitted if not needed. Whitespace is
451 allowed for readability but can also be omitted.
452 '''
454 start = Float.T(optional=True)
455 stop = Float.T(optional=True)
456 step = Float.T(optional=True)
457 n = Int.T(optional=True)
458 values = Array.T(optional=True, dtype=float, shape=(None,))
460 spacing = StringChoice.T(
461 choices=['lin', 'log', 'symlog'],
462 default='lin',
463 optional=True)
465 relative = StringChoice.T(
466 choices=['', 'add', 'mult'],
467 default='',
468 optional=True)
470 pattern = re.compile(r'^((?P<start>.*)\.\.(?P<stop>[^@|:]*))?'
471 r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?'
472 r'(\|(?P<stuff>.+))?$')
474 def __init__(self, *args, **kwargs):
475 d = {}
476 if len(args) == 1:
477 d = self.parse(args[0])
478 elif len(args) in (2, 3):
479 d['start'], d['stop'] = [float(x) for x in args[:2]]
480 if len(args) == 3:
481 d['step'] = float(args[2])
483 for k, v in kwargs.items():
484 if k in d:
485 raise ArgumentError('%s specified more than once' % k)
487 d[k] = v
489 SObject.__init__(self, **d)
491 def __str__(self):
492 def sfloat(x):
493 if x is not None:
494 return '%g' % x
495 else:
496 return ''
498 if self.values:
499 return ','.join('%g' % x for x in self.values)
501 if self.start is None and self.stop is None:
502 s0 = ''
503 else:
504 s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop))
506 s1 = ''
507 if self.step is not None:
508 s1 = [' : %g', ':%g'][s0 == ''] % self.step
509 elif self.n is not None:
510 s1 = [' @ %i', '@%i'][s0 == ''] % self.n
512 if self.spacing == 'lin' and self.relative == '':
513 s2 = ''
514 else:
515 x = []
516 if self.spacing != 'lin':
517 x.append(self.spacing)
519 if self.relative != '':
520 x.append(self.relative)
522 s2 = ' | %s' % ','.join(x)
524 return s0 + s1 + s2
526 @classmethod
527 def parse(cls, s):
528 s = re.sub(r'\s+', '', s)
529 m = cls.pattern.match(s)
530 if not m:
531 try:
532 vals = [ufloat(x) for x in s.split(',')]
533 except Exception:
534 raise InvalidGridDef(
535 '"%s" is not a valid range specification' % s)
537 return dict(values=num.array(vals, dtype=float))
539 d = m.groupdict()
540 try:
541 start = ufloat_or_none(d['start'])
542 stop = ufloat_or_none(d['stop'])
543 step = ufloat_or_none(d['step'])
544 n = int_or_none(d['n'])
545 except Exception:
546 raise InvalidGridDef(
547 '"%s" is not a valid range specification' % s)
549 spacing = 'lin'
550 relative = ''
552 if d['stuff'] is not None:
553 t = d['stuff'].split(',')
554 for x in t:
555 if x in cls.spacing.choices:
556 spacing = x
557 elif x and x in cls.relative.choices:
558 relative = x
559 else:
560 raise InvalidGridDef(
561 '"%s" is not a valid range specification' % s)
563 return dict(start=start, stop=stop, step=step, n=n, spacing=spacing,
564 relative=relative)
566 def make(self, mi=None, ma=None, inc=None, base=None, eps=1e-5):
567 if self.values:
568 return self.values
570 start = self.start
571 stop = self.stop
572 step = self.step
573 n = self.n
575 swap = step is not None and step < 0.
576 if start is None:
577 start = [mi, ma][swap]
578 if stop is None:
579 stop = [ma, mi][swap]
580 if step is None and inc is not None:
581 step = [inc, -inc][ma < mi]
583 if start is None or stop is None:
584 raise InvalidGridDef(
585 'Cannot use range specification "%s" without start '
586 'and stop in this context' % self)
588 if step is None and n is None:
589 step = stop - start
591 if n is None:
592 if (step < 0) != (stop - start < 0):
593 raise InvalidGridDef(
594 'Range specification "%s" has inconsistent ordering '
595 '(step < 0 => stop > start)' % self)
597 n = int(round((stop - start) / step)) + 1
598 stop2 = start + (n - 1) * step
599 if abs(stop - stop2) > eps:
600 n = int(math.floor((stop - start) / step)) + 1
601 stop = start + (n - 1) * step
602 else:
603 stop = stop2
605 if start == stop:
606 n = 1
608 if self.spacing == 'lin':
609 vals = num.linspace(start, stop, n)
611 elif self.spacing in ('log', 'symlog'):
612 if start > 0. and stop > 0.:
613 vals = num.exp(num.linspace(num.log(start),
614 num.log(stop), n))
615 elif start < 0. and stop < 0.:
616 vals = -num.exp(num.linspace(num.log(-start),
617 num.log(-stop), n))
618 else:
619 raise InvalidGridDef(
620 'Log ranges should not include or cross zero '
621 '(in range specification "%s").' % self)
623 if self.spacing == 'symlog':
624 nvals = - vals
625 vals = num.concatenate((nvals[::-1], vals))
627 if self.relative in ('add', 'mult') and base is None:
628 raise InvalidGridDef(
629 'Cannot use relative range specification in this context.')
631 vals = self.make_relative(base, vals)
633 return list(map(float, vals))
635 def make_relative(self, base, vals):
636 if self.relative == 'add':
637 vals += base
639 if self.relative == 'mult':
640 vals *= base
642 return vals
645class GridDefElement(Object):
647 param = meta.StringID.T()
648 rs = Range.T()
650 def __init__(self, shorthand=None, **kwargs):
651 if shorthand is not None:
652 t = shorthand.split('=')
653 if len(t) != 2:
654 raise InvalidGridDef(
655 'Invalid grid specification element: %s' % shorthand)
657 sp, sr = t[0].strip(), t[1].strip()
659 kwargs['param'] = sp
660 kwargs['rs'] = Range(sr)
662 Object.__init__(self, **kwargs)
664 def shorthand(self):
665 return self.param + ' = ' + str(self.rs)
668class GridDef(Object):
670 elements = List.T(GridDefElement.T())
672 def __init__(self, shorthand=None, **kwargs):
673 if shorthand is not None:
674 t = shorthand.splitlines()
675 tt = []
676 for x in t:
677 x = x.strip()
678 if x:
679 tt.extend(x.split(';'))
681 elements = []
682 for se in tt:
683 elements.append(GridDef(se))
685 kwargs['elements'] = elements
687 Object.__init__(self, **kwargs)
689 def shorthand(self):
690 return '; '.join(str(x) for x in self.elements)
693class Cloneable(object):
694 '''
695 Mix-in class for Guts objects, providing dict-like access and cloning.
696 '''
698 def __iter__(self):
699 return iter(self.T.propnames)
701 def __getitem__(self, k):
702 if k not in self.keys():
703 raise KeyError(k)
705 return getattr(self, k)
707 def __setitem__(self, k, v):
708 if k not in self.keys():
709 raise KeyError(k)
711 return setattr(self, k, v)
713 def clone(self, **kwargs):
714 '''
715 Make a copy of the object.
717 A new object of the same class is created and initialized with the
718 parameters of the object on which this method is called on. If
719 ``kwargs`` are given, these are used to override any of the
720 initialization parameters.
721 '''
723 d = dict(self)
724 for k in d:
725 v = d[k]
726 if isinstance(v, Cloneable):
727 d[k] = v.clone()
729 d.update(kwargs)
730 return self.__class__(**d)
732 @classmethod
733 def keys(cls):
734 '''
735 Get list of the source model's parameter names.
736 '''
738 return cls.T.propnames
741class STF(Object, Cloneable):
743 '''
744 Base class for source time functions.
745 '''
747 def __init__(self, effective_duration=None, **kwargs):
748 if effective_duration is not None:
749 kwargs['duration'] = effective_duration / \
750 self.factor_duration_to_effective()
752 Object.__init__(self, **kwargs)
754 @classmethod
755 def factor_duration_to_effective(cls):
756 return 1.0
758 def centroid_time(self, tref):
759 return tref
761 @property
762 def effective_duration(self):
763 return self.duration * self.factor_duration_to_effective()
765 def discretize_t(self, deltat, tref):
766 tl = math.floor(tref / deltat) * deltat
767 th = math.ceil(tref / deltat) * deltat
768 if tl == th:
769 return num.array([tl], dtype=float), num.ones(1)
770 else:
771 return (
772 num.array([tl, th], dtype=float),
773 num.array([th - tref, tref - tl], dtype=float) / deltat)
775 def base_key(self):
776 return (type(self).__name__,)
779g_unit_pulse = STF()
782def sshift(times, amplitudes, tshift, deltat):
784 t0 = math.floor(tshift / deltat) * deltat
785 t1 = math.ceil(tshift / deltat) * deltat
786 if t0 == t1:
787 return times, amplitudes
789 amplitudes2 = num.zeros(amplitudes.size + 1, dtype=float)
791 amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes
792 amplitudes2[1:] += (tshift - t0) / deltat * amplitudes
794 times2 = num.arange(times.size + 1, dtype=float) * \
795 deltat + times[0] + t0
797 return times2, amplitudes2
800class BoxcarSTF(STF):
802 '''
803 Boxcar type source time function.
805 .. figure :: /static/stf-BoxcarSTF.svg
806 :width: 40%
807 :align: center
808 :alt: boxcar source time function
809 '''
811 duration = Float.T(
812 default=0.0,
813 help='duration of the boxcar')
815 anchor = Float.T(
816 default=0.0,
817 help='anchor point with respect to source.time: ('
818 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
819 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
820 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
822 @classmethod
823 def factor_duration_to_effective(cls):
824 return 1.0
826 def centroid_time(self, tref):
827 return tref - 0.5 * self.duration * self.anchor
829 def discretize_t(self, deltat, tref):
830 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
831 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
832 tmin = round(tmin_stf / deltat) * deltat
833 tmax = round(tmax_stf / deltat) * deltat
834 nt = int(round((tmax - tmin) / deltat)) + 1
835 times = num.linspace(tmin, tmax, nt)
836 amplitudes = num.ones_like(times)
837 if times.size > 1:
838 t_edges = num.linspace(
839 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
840 t = tmin_stf + self.duration * num.array(
841 [0.0, 0.0, 1.0, 1.0], dtype=float)
842 f = num.array([0., 1., 1., 0.], dtype=float)
843 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
844 amplitudes /= num.sum(amplitudes)
846 tshift = (num.sum(amplitudes * times) - self.centroid_time(tref))
848 return sshift(times, amplitudes, -tshift, deltat)
850 def base_key(self):
851 return (type(self).__name__, self.duration, self.anchor)
854class TriangularSTF(STF):
856 '''
857 Triangular type source time function.
859 .. figure :: /static/stf-TriangularSTF.svg
860 :width: 40%
861 :align: center
862 :alt: triangular source time function
863 '''
865 duration = Float.T(
866 default=0.0,
867 help='baseline of the triangle')
869 peak_ratio = Float.T(
870 default=0.5,
871 help='fraction of time compared to duration, '
872 'when the maximum amplitude is reached')
874 anchor = Float.T(
875 default=0.0,
876 help='anchor point with respect to source.time: ('
877 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
878 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
879 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
881 @classmethod
882 def factor_duration_to_effective(cls, peak_ratio=None):
883 if peak_ratio is None:
884 peak_ratio = cls.peak_ratio.default()
886 return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0)
888 def __init__(self, effective_duration=None, **kwargs):
889 if effective_duration is not None:
890 kwargs['duration'] = effective_duration / \
891 self.factor_duration_to_effective(
892 kwargs.get('peak_ratio', None))
894 STF.__init__(self, **kwargs)
896 @property
897 def centroid_ratio(self):
898 ra = self.peak_ratio
899 rb = 1.0 - ra
900 return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb)
902 def centroid_time(self, tref):
903 ca = self.centroid_ratio
904 cb = 1.0 - ca
905 if self.anchor <= 0.:
906 return tref - ca * self.duration * self.anchor
907 else:
908 return tref - cb * self.duration * self.anchor
910 @property
911 def effective_duration(self):
912 return self.duration * self.factor_duration_to_effective(
913 self.peak_ratio)
915 def tminmax_stf(self, tref):
916 ca = self.centroid_ratio
917 cb = 1.0 - ca
918 if self.anchor <= 0.:
919 tmin_stf = tref - ca * self.duration * (self.anchor + 1.)
920 tmax_stf = tmin_stf + self.duration
921 else:
922 tmax_stf = tref + cb * self.duration * (1. - self.anchor)
923 tmin_stf = tmax_stf - self.duration
925 return tmin_stf, tmax_stf
927 def discretize_t(self, deltat, tref):
928 tmin_stf, tmax_stf = self.tminmax_stf(tref)
930 tmin = round(tmin_stf / deltat) * deltat
931 tmax = round(tmax_stf / deltat) * deltat
932 nt = int(round((tmax - tmin) / deltat)) + 1
933 if nt > 1:
934 t_edges = num.linspace(
935 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
936 t = tmin_stf + self.duration * num.array(
937 [0.0, self.peak_ratio, 1.0], dtype=float)
938 f = num.array([0., 1., 0.], dtype=float)
939 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
940 amplitudes /= num.sum(amplitudes)
941 else:
942 amplitudes = num.ones(1)
944 times = num.linspace(tmin, tmax, nt)
945 return times, amplitudes
947 def base_key(self):
948 return (
949 type(self).__name__, self.duration, self.peak_ratio, self.anchor)
952class HalfSinusoidSTF(STF):
954 '''
955 Half sinusoid type source time function.
957 .. figure :: /static/stf-HalfSinusoidSTF.svg
958 :width: 40%
959 :align: center
960 :alt: half-sinusouid source time function
961 '''
963 duration = Float.T(
964 default=0.0,
965 help='duration of the half-sinusoid (baseline)')
967 anchor = Float.T(
968 default=0.0,
969 help='anchor point with respect to source.time: ('
970 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
971 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
972 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
974 exponent = Int.T(
975 default=1,
976 help='set to 2 to use square of the half-period sinusoidal function.')
978 def __init__(self, effective_duration=None, **kwargs):
979 if effective_duration is not None:
980 kwargs['duration'] = effective_duration / \
981 self.factor_duration_to_effective(
982 kwargs.get('exponent', 1))
984 STF.__init__(self, **kwargs)
986 @classmethod
987 def factor_duration_to_effective(cls, exponent):
988 if exponent == 1:
989 return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi
990 elif exponent == 2:
991 return math.sqrt(math.pi**2 - 6) / math.pi
992 else:
993 raise ValueError('Exponent for HalfSinusoidSTF must be 1 or 2.')
995 @property
996 def effective_duration(self):
997 return self.duration * self.factor_duration_to_effective(self.exponent)
999 def centroid_time(self, tref):
1000 return tref - 0.5 * self.duration * self.anchor
1002 def discretize_t(self, deltat, tref):
1003 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1004 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1005 tmin = round(tmin_stf / deltat) * deltat
1006 tmax = round(tmax_stf / deltat) * deltat
1007 nt = int(round((tmax - tmin) / deltat)) + 1
1008 if nt > 1:
1009 t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace(
1010 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)))
1012 if self.exponent == 1:
1013 fint = -num.cos(
1014 (t_edges - tmin_stf) * (math.pi / self.duration))
1016 elif self.exponent == 2:
1017 fint = (t_edges - tmin_stf) / self.duration \
1018 - 1.0 / (2.0 * math.pi) * num.sin(
1019 (t_edges - tmin_stf) * (2.0 * math.pi / self.duration))
1020 else:
1021 raise ValueError(
1022 'Exponent for HalfSinusoidSTF must be 1 or 2.')
1024 amplitudes = fint[1:] - fint[:-1]
1025 amplitudes /= num.sum(amplitudes)
1026 else:
1027 amplitudes = num.ones(1)
1029 times = num.linspace(tmin, tmax, nt)
1030 return times, amplitudes
1032 def base_key(self):
1033 return (type(self).__name__, self.duration, self.anchor)
1036class SmoothRampSTF(STF):
1037 '''
1038 Smooth-ramp type source time function for near-field displacement.
1039 Based on moment function of double-couple point source proposed by [1]_.
1041 .. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow
1042 earthquakes from Love-wave modelling for regional distances, PEPI 32,
1043 312-324.
1045 .. figure :: /static/stf-SmoothRampSTF.svg
1046 :width: 40%
1047 :alt: smooth ramp source time function
1048 '''
1049 duration = Float.T(
1050 default=0.0,
1051 help='duration of the ramp (baseline)')
1053 rise_ratio = Float.T(
1054 default=0.5,
1055 help='fraction of time compared to duration, '
1056 'when the maximum amplitude is reached')
1058 anchor = Float.T(
1059 default=0.0,
1060 help='anchor point with respect to source.time: ('
1061 '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, '
1062 '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, '
1063 '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)')
1065 def discretize_t(self, deltat, tref):
1066 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1067 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1068 tmin = round(tmin_stf / deltat) * deltat
1069 tmax = round(tmax_stf / deltat) * deltat
1070 D = round((tmax - tmin) / deltat) * deltat
1071 nt = int(round(D / deltat)) + 1
1072 times = num.linspace(tmin, tmax, nt)
1073 if nt > 1:
1074 rise_time = self.rise_ratio * self.duration
1075 amplitudes = num.ones_like(times)
1076 tp = tmin + rise_time
1077 ii = num.where(times <= tp)
1078 t_inc = times[ii]
1079 a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time)
1080 b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0
1081 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b)
1083 amplitudes /= num.sum(amplitudes)
1084 else:
1085 amplitudes = num.ones(1)
1087 return times, amplitudes
1089 def base_key(self):
1090 return (type(self).__name__,
1091 self.duration, self.rise_ratio, self.anchor)
1094class ResonatorSTF(STF):
1095 '''
1096 Simple resonator like source time function.
1098 .. math ::
1100 f(t) = 0 for t < 0
1101 f(t) = e^{-t/tau} * sin(2 * pi * f * t)
1104 .. figure :: /static/stf-SmoothRampSTF.svg
1105 :width: 40%
1106 :alt: smooth ramp source time function
1108 '''
1110 duration = Float.T(
1111 default=0.0,
1112 help='decay time')
1114 frequency = Float.T(
1115 default=1.0,
1116 help='resonance frequency')
1118 def discretize_t(self, deltat, tref):
1119 tmin_stf = tref
1120 tmax_stf = tref + self.duration * 3
1121 tmin = math.floor(tmin_stf / deltat) * deltat
1122 tmax = math.ceil(tmax_stf / deltat) * deltat
1123 times = util.arange2(tmin, tmax, deltat)
1124 amplitudes = num.exp(-(times - tref) / self.duration) \
1125 * num.sin(2.0 * num.pi * self.frequency * (times - tref))
1127 return times, amplitudes
1129 def base_key(self):
1130 return (type(self).__name__,
1131 self.duration, self.frequency)
1134class TremorSTF(STF):
1135 '''
1136 Oszillating source time function.
1138 .. math ::
1140 f(t) = 0 for t < -tau/2 or t > tau/2
1141 f(t) = cos(pi/tau*t) * sin(2 * pi * f * t)
1143 '''
1145 duration = Float.T(
1146 default=0.0,
1147 help='Tremor duration [s]')
1149 frequency = Float.T(
1150 default=1.0,
1151 help='Frequency [Hz]')
1153 def discretize_t(self, deltat, tref):
1154 tmin_stf = tref - 0.5 * self.duration
1155 tmax_stf = tref + 0.5 * self.duration
1156 tmin = math.floor(tmin_stf / deltat) * deltat
1157 tmax = math.ceil(tmax_stf / deltat) * deltat
1158 times = util.arange2(tmin, tmax, deltat)
1159 mask = num.logical_and(
1160 tref - 0.5 * self.duration < times,
1161 times < tref + 0.5 * self.duration)
1163 amplitudes = num.zeros_like(times)
1164 amplitudes[mask] = num.cos(num.pi/self.duration*(times[mask] - tref)) \
1165 * num.sin(2.0 * num.pi * self.frequency * (times[mask] - tref))
1167 return times, amplitudes
1169 def base_key(self):
1170 return (type(self).__name__,
1171 self.duration, self.frequency)
1174class STFMode(StringChoice):
1175 choices = ['pre', 'post']
1178class Source(Location, Cloneable):
1179 '''
1180 Base class for all source models.
1181 '''
1183 name = String.T(optional=True, default='')
1185 time = Timestamp.T(
1186 default=Timestamp.D('1970-01-01 00:00:00'),
1187 help='source origin time.')
1189 stf = STF.T(
1190 optional=True,
1191 help='source time function.')
1193 stf_mode = STFMode.T(
1194 default='post',
1195 help='whether to apply source time function in pre or '
1196 'post-processing.')
1198 def __init__(self, **kwargs):
1199 Location.__init__(self, **kwargs)
1201 def update(self, **kwargs):
1202 '''
1203 Change some of the source models parameters.
1205 Example::
1207 >>> from pyrocko import gf
1208 >>> s = gf.DCSource()
1209 >>> s.update(strike=66., dip=33.)
1210 >>> print(s)
1211 --- !pf.DCSource
1212 depth: 0.0
1213 time: 1970-01-01 00:00:00
1214 magnitude: 6.0
1215 strike: 66.0
1216 dip: 33.0
1217 rake: 0.0
1219 '''
1221 for (k, v) in kwargs.items():
1222 self[k] = v
1224 def grid(self, **variables):
1225 '''
1226 Create grid of source model variations.
1228 :returns: :py:class:`SourceGrid` instance.
1230 Example::
1232 >>> from pyrocko import gf
1233 >>> base = DCSource()
1234 >>> R = gf.Range
1235 >>> for s in base.grid(R('
1237 '''
1238 return SourceGrid(base=self, variables=variables)
1240 def base_key(self):
1241 '''
1242 Get key to decide about source discretization / GF stack sharing.
1244 When two source models differ only in amplitude and origin time, the
1245 discretization and the GF stacking can be done only once for a unit
1246 amplitude and a zero origin time and the amplitude and origin times of
1247 the seismograms can be applied during post-processing of the synthetic
1248 seismogram.
1250 For any derived parameterized source model, this method is called to
1251 decide if discretization and stacking of the source should be shared.
1252 When two source models return an equal vector of values discretization
1253 is shared.
1254 '''
1255 return (self.depth, self.lat, self.north_shift,
1256 self.lon, self.east_shift, self.time, type(self).__name__) + \
1257 self.effective_stf_pre().base_key()
1259 def get_factor(self):
1260 '''
1261 Get the scaling factor to be applied during post-processing.
1263 Discretization of the base seismogram is usually done for a unit
1264 amplitude, because a common factor can be efficiently multiplied to
1265 final seismograms. This eliminates to do repeat the stacking when
1266 creating seismograms for a series of source models only differing in
1267 amplitude.
1269 This method should return the scaling factor to apply in the
1270 post-processing (often this is simply the scalar moment of the source).
1271 '''
1273 return 1.0
1275 def effective_stf_pre(self):
1276 '''
1277 Return the STF applied before stacking of the Green's functions.
1279 This STF is used during discretization of the parameterized source
1280 models, i.e. to produce a temporal distribution of point sources.
1282 Handling of the STF before stacking of the GFs is less efficient but
1283 allows to use different source time functions for different parts of
1284 the source.
1285 '''
1287 if self.stf is not None and self.stf_mode == 'pre':
1288 return self.stf
1289 else:
1290 return g_unit_pulse
1292 def effective_stf_post(self):
1293 '''
1294 Return the STF applied after stacking of the Green's fuctions.
1296 This STF is used in the post-processing of the synthetic seismograms.
1298 Handling of the STF after stacking of the GFs is usually more efficient
1299 but is only possible when a common STF is used for all subsources.
1300 '''
1302 if self.stf is not None and self.stf_mode == 'post':
1303 return self.stf
1304 else:
1305 return g_unit_pulse
1307 def _dparams_base(self):
1308 return dict(times=arr(self.time),
1309 lat=self.lat, lon=self.lon,
1310 north_shifts=arr(self.north_shift),
1311 east_shifts=arr(self.east_shift),
1312 depths=arr(self.depth))
1314 def _hash(self):
1315 sha = sha1()
1316 for k in self.base_key():
1317 sha.update(str(k).encode())
1318 return sha.hexdigest()
1320 def _dparams_base_repeated(self, times):
1321 if times is None:
1322 return self._dparams_base()
1324 nt = times.size
1325 north_shifts = num.repeat(self.north_shift, nt)
1326 east_shifts = num.repeat(self.east_shift, nt)
1327 depths = num.repeat(self.depth, nt)
1328 return dict(times=times,
1329 lat=self.lat, lon=self.lon,
1330 north_shifts=north_shifts,
1331 east_shifts=east_shifts,
1332 depths=depths)
1334 def pyrocko_event(self, store=None, target=None, **kwargs):
1335 duration = None
1336 if self.stf:
1337 duration = self.stf.effective_duration
1339 return model.Event(
1340 lat=self.lat,
1341 lon=self.lon,
1342 north_shift=self.north_shift,
1343 east_shift=self.east_shift,
1344 time=self.time,
1345 name=self.name,
1346 depth=self.depth,
1347 duration=duration,
1348 **kwargs)
1350 def geometry(self, **kwargs):
1351 raise NotImplementedError
1353 def outline(self, cs='xyz'):
1354 points = num.atleast_2d(num.zeros([1, 3]))
1356 points[:, 0] += self.north_shift
1357 points[:, 1] += self.east_shift
1358 points[:, 2] += self.depth
1359 if cs == 'xyz':
1360 return points
1361 elif cs == 'xy':
1362 return points[:, :2]
1363 elif cs in ('latlon', 'lonlat'):
1364 latlon = ne_to_latlon(
1365 self.lat, self.lon, points[:, 0], points[:, 1])
1367 latlon = num.array(latlon).T
1368 if cs == 'latlon':
1369 return latlon
1370 else:
1371 return latlon[:, ::-1]
1373 @classmethod
1374 def from_pyrocko_event(cls, ev, **kwargs):
1375 if ev.depth is None:
1376 raise ConversionError(
1377 'Cannot convert event object to source object: '
1378 'no depth information available')
1380 stf = None
1381 if ev.duration is not None:
1382 stf = HalfSinusoidSTF(effective_duration=ev.duration)
1384 d = dict(
1385 name=ev.name,
1386 time=ev.time,
1387 lat=ev.lat,
1388 lon=ev.lon,
1389 north_shift=ev.north_shift,
1390 east_shift=ev.east_shift,
1391 depth=ev.depth,
1392 stf=stf)
1393 d.update(kwargs)
1394 return cls(**d)
1396 def get_magnitude(self):
1397 raise NotImplementedError(
1398 '%s does not implement get_magnitude()'
1399 % self.__class__.__name__)
1402class SourceWithMagnitude(Source):
1403 '''
1404 Base class for sources containing a moment magnitude.
1405 '''
1407 magnitude = Float.T(
1408 default=6.0,
1409 help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1411 def __init__(self, **kwargs):
1412 if 'moment' in kwargs:
1413 mom = kwargs.pop('moment')
1414 if 'magnitude' not in kwargs:
1415 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1417 Source.__init__(self, **kwargs)
1419 @property
1420 def moment(self):
1421 return float(pmt.magnitude_to_moment(self.magnitude))
1423 @moment.setter
1424 def moment(self, value):
1425 self.magnitude = float(pmt.moment_to_magnitude(value))
1427 def pyrocko_event(self, store=None, target=None, **kwargs):
1428 return Source.pyrocko_event(
1429 self, store, target,
1430 magnitude=self.magnitude,
1431 **kwargs)
1433 @classmethod
1434 def from_pyrocko_event(cls, ev, **kwargs):
1435 d = {}
1436 if ev.magnitude:
1437 d.update(magnitude=ev.magnitude)
1439 d.update(kwargs)
1440 return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
1442 def get_magnitude(self):
1443 return self.magnitude
1446class DerivedMagnitudeError(ValidationError):
1447 '''
1448 Raised when conversion between magnitude, moment, volume change or
1449 displacement failed.
1450 '''
1451 pass
1454class SourceWithDerivedMagnitude(Source):
1456 class __T(Source.T):
1458 def validate_extra(self, val):
1459 Source.T.validate_extra(self, val)
1460 val.check_conflicts()
1462 def check_conflicts(self):
1463 '''
1464 Check for parameter conflicts.
1466 To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError`
1467 on conflicts.
1468 '''
1469 pass
1471 def get_magnitude(self, store=None, target=None):
1472 raise DerivedMagnitudeError('No magnitude set.')
1474 def get_moment(self, store=None, target=None):
1475 return float(pmt.magnitude_to_moment(
1476 self.get_magnitude(store, target)))
1478 def pyrocko_moment_tensor(self, store=None, target=None):
1479 raise NotImplementedError(
1480 '%s does not implement pyrocko_moment_tensor()'
1481 % self.__class__.__name__)
1483 def pyrocko_event(self, store=None, target=None, **kwargs):
1484 try:
1485 mt = self.pyrocko_moment_tensor(store, target)
1486 magnitude = self.get_magnitude()
1487 except (DerivedMagnitudeError, NotImplementedError):
1488 mt = None
1489 magnitude = None
1491 return Source.pyrocko_event(
1492 self, store, target,
1493 moment_tensor=mt,
1494 magnitude=magnitude,
1495 **kwargs)
1498class ExplosionSource(SourceWithDerivedMagnitude):
1499 '''
1500 An isotropic explosion point source.
1501 '''
1503 magnitude = Float.T(
1504 optional=True,
1505 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1507 volume_change = Float.T(
1508 optional=True,
1509 help='volume change of the explosion/implosion or '
1510 'the contracting/extending magmatic source. [m^3]')
1512 discretized_source_class = meta.DiscretizedExplosionSource
1514 def __init__(self, **kwargs):
1515 if 'moment' in kwargs:
1516 mom = kwargs.pop('moment')
1517 if 'magnitude' not in kwargs:
1518 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1520 SourceWithDerivedMagnitude.__init__(self, **kwargs)
1522 def base_key(self):
1523 return SourceWithDerivedMagnitude.base_key(self) + \
1524 (self.volume_change,)
1526 def check_conflicts(self):
1527 if self.magnitude is not None and self.volume_change is not None:
1528 raise DerivedMagnitudeError(
1529 'Magnitude and volume_change are both defined.')
1531 def get_magnitude(self, store=None, target=None):
1532 self.check_conflicts()
1534 if self.magnitude is not None:
1535 return self.magnitude
1537 elif self.volume_change is not None:
1538 moment = self.volume_change * \
1539 self.get_moment_to_volume_change_ratio(store, target)
1541 return float(pmt.moment_to_magnitude(abs(moment)))
1542 else:
1543 return float(pmt.moment_to_magnitude(1.0))
1545 def get_volume_change(self, store=None, target=None):
1546 self.check_conflicts()
1548 if self.volume_change is not None:
1549 return self.volume_change
1551 elif self.magnitude is not None:
1552 moment = float(pmt.magnitude_to_moment(self.magnitude))
1553 return moment / self.get_moment_to_volume_change_ratio(
1554 store, target)
1556 else:
1557 return 1.0 / self.get_moment_to_volume_change_ratio(store)
1559 def get_moment_to_volume_change_ratio(self, store, target=None):
1560 if store is None:
1561 raise DerivedMagnitudeError(
1562 'Need earth model to convert between volume change and '
1563 'magnitude.')
1565 points = num.array(
1566 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1568 interpolation = target.interpolation if target else 'multilinear'
1569 try:
1570 shear_moduli = store.config.get_shear_moduli(
1571 self.lat, self.lon,
1572 points=points,
1573 interpolation=interpolation)[0]
1575 bulk_moduli = store.config.get_bulk_moduli(
1576 self.lat, self.lon,
1577 points=points,
1578 interpolation=interpolation)[0]
1579 except meta.OutOfBounds:
1580 raise DerivedMagnitudeError(
1581 'Could not get shear modulus at source position.')
1583 return float(2. * shear_moduli + bulk_moduli)
1585 def get_factor(self):
1586 return 1.0
1588 def discretize_basesource(self, store, target=None):
1589 times, amplitudes = self.effective_stf_pre().discretize_t(
1590 store.config.deltat, self.time)
1592 amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.)
1594 if self.volume_change is not None:
1595 if self.volume_change < 0.:
1596 amplitudes *= -1
1598 return meta.DiscretizedExplosionSource(
1599 m0s=amplitudes,
1600 **self._dparams_base_repeated(times))
1602 def pyrocko_moment_tensor(self, store=None, target=None):
1603 a = self.get_moment(store, target) * math.sqrt(2. / 3.)
1604 return pmt.MomentTensor(m=pmt.symmat6(a, a, a, 0., 0., 0.))
1607class RectangularExplosionSource(ExplosionSource):
1608 '''
1609 Rectangular or line explosion source.
1610 '''
1612 discretized_source_class = meta.DiscretizedExplosionSource
1614 strike = Float.T(
1615 default=0.0,
1616 help='strike direction in [deg], measured clockwise from north')
1618 dip = Float.T(
1619 default=90.0,
1620 help='dip angle in [deg], measured downward from horizontal')
1622 length = Float.T(
1623 default=0.,
1624 help='length of rectangular source area [m]')
1626 width = Float.T(
1627 default=0.,
1628 help='width of rectangular source area [m]')
1630 anchor = StringChoice.T(
1631 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
1632 'bottom_left', 'bottom_right'],
1633 default='center',
1634 optional=True,
1635 help='Anchor point for positioning the plane, can be: top, center or'
1636 'bottom and also top_left, top_right,bottom_left,'
1637 'bottom_right, center_left and center right')
1639 nucleation_x = Float.T(
1640 optional=True,
1641 help='horizontal position of rupture nucleation in normalized fault '
1642 'plane coordinates (-1 = left edge, +1 = right edge)')
1644 nucleation_y = Float.T(
1645 optional=True,
1646 help='down-dip position of rupture nucleation in normalized fault '
1647 'plane coordinates (-1 = upper edge, +1 = lower edge)')
1649 velocity = Float.T(
1650 default=3500.,
1651 help='speed of explosion front [m/s]')
1653 aggressive_oversampling = Bool.T(
1654 default=False,
1655 help='Aggressive oversampling for basesource discretization. '
1656 "When using 'multilinear' interpolation oversampling has"
1657 ' practically no effect.')
1659 def base_key(self):
1660 return Source.base_key(self) + (self.strike, self.dip, self.length,
1661 self.width, self.nucleation_x,
1662 self.nucleation_y, self.velocity,
1663 self.anchor)
1665 def discretize_basesource(self, store, target=None):
1667 if self.nucleation_x is not None:
1668 nucx = self.nucleation_x * 0.5 * self.length
1669 else:
1670 nucx = None
1672 if self.nucleation_y is not None:
1673 nucy = self.nucleation_y * 0.5 * self.width
1674 else:
1675 nucy = None
1677 stf = self.effective_stf_pre()
1679 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
1680 store.config.deltas, store.config.deltat,
1681 self.time, self.north_shift, self.east_shift, self.depth,
1682 self.strike, self.dip, self.length, self.width, self.anchor,
1683 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy)
1685 amplitudes /= num.sum(amplitudes)
1686 amplitudes *= self.get_moment(store, target)
1688 return meta.DiscretizedExplosionSource(
1689 lat=self.lat,
1690 lon=self.lon,
1691 times=times,
1692 north_shifts=points[:, 0],
1693 east_shifts=points[:, 1],
1694 depths=points[:, 2],
1695 m0s=amplitudes)
1697 def outline(self, cs='xyz'):
1698 points = outline_rect_source(self.strike, self.dip, self.length,
1699 self.width, self.anchor)
1701 points[:, 0] += self.north_shift
1702 points[:, 1] += self.east_shift
1703 points[:, 2] += self.depth
1704 if cs == 'xyz':
1705 return points
1706 elif cs == 'xy':
1707 return points[:, :2]
1708 elif cs in ('latlon', 'lonlat'):
1709 latlon = ne_to_latlon(
1710 self.lat, self.lon, points[:, 0], points[:, 1])
1712 latlon = num.array(latlon).T
1713 if cs == 'latlon':
1714 return latlon
1715 else:
1716 return latlon[:, ::-1]
1718 def get_nucleation_abs_coord(self, cs='xy'):
1720 if self.nucleation_x is None:
1721 return None, None
1723 coords = from_plane_coords(self.strike, self.dip, self.length,
1724 self.width, self.depth, self.nucleation_x,
1725 self.nucleation_y, lat=self.lat,
1726 lon=self.lon, north_shift=self.north_shift,
1727 east_shift=self.east_shift, cs=cs)
1728 return coords
1731class DCSource(SourceWithMagnitude):
1732 '''
1733 A double-couple point source.
1734 '''
1736 strike = Float.T(
1737 default=0.0,
1738 help='strike direction in [deg], measured clockwise from north')
1740 dip = Float.T(
1741 default=90.0,
1742 help='dip angle in [deg], measured downward from horizontal')
1744 rake = Float.T(
1745 default=0.0,
1746 help='rake angle in [deg], '
1747 'measured counter-clockwise from right-horizontal '
1748 'in on-plane view')
1750 discretized_source_class = meta.DiscretizedMTSource
1752 def base_key(self):
1753 return Source.base_key(self) + (self.strike, self.dip, self.rake)
1755 def get_factor(self):
1756 return float(pmt.magnitude_to_moment(self.magnitude))
1758 def discretize_basesource(self, store, target=None):
1759 mot = pmt.MomentTensor(
1760 strike=self.strike, dip=self.dip, rake=self.rake)
1762 times, amplitudes = self.effective_stf_pre().discretize_t(
1763 store.config.deltat, self.time)
1764 return meta.DiscretizedMTSource(
1765 m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis],
1766 **self._dparams_base_repeated(times))
1768 def pyrocko_moment_tensor(self, store=None, target=None):
1769 return pmt.MomentTensor(
1770 strike=self.strike,
1771 dip=self.dip,
1772 rake=self.rake,
1773 scalar_moment=self.moment)
1775 def pyrocko_event(self, store=None, target=None, **kwargs):
1776 return SourceWithMagnitude.pyrocko_event(
1777 self, store, target,
1778 moment_tensor=self.pyrocko_moment_tensor(store, target),
1779 **kwargs)
1781 @classmethod
1782 def from_pyrocko_event(cls, ev, **kwargs):
1783 d = {}
1784 mt = ev.moment_tensor
1785 if mt:
1786 (strike, dip, rake), _ = mt.both_strike_dip_rake()
1787 d.update(
1788 strike=float(strike),
1789 dip=float(dip),
1790 rake=float(rake),
1791 magnitude=float(mt.moment_magnitude()))
1793 d.update(kwargs)
1794 return super(DCSource, cls).from_pyrocko_event(ev, **d)
1797class CLVDSource(SourceWithMagnitude):
1798 '''
1799 A pure CLVD point source.
1800 '''
1802 discretized_source_class = meta.DiscretizedMTSource
1804 azimuth = Float.T(
1805 default=0.0,
1806 help='azimuth direction of largest dipole, clockwise from north [deg]')
1808 dip = Float.T(
1809 default=90.,
1810 help='dip direction of largest dipole, downward from horizontal [deg]')
1812 def base_key(self):
1813 return Source.base_key(self) + (self.azimuth, self.dip)
1815 def get_factor(self):
1816 return float(pmt.magnitude_to_moment(self.magnitude))
1818 @property
1819 def m6(self):
1820 a = math.sqrt(4. / 3.) * self.get_factor()
1821 m = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1822 rotmat1 = pmt.euler_to_matrix(
1823 d2r * (self.dip - 90.),
1824 d2r * (self.azimuth - 90.),
1825 0.)
1826 m = num.dot(rotmat1.T, num.dot(m, rotmat1))
1827 return pmt.to6(m)
1829 @property
1830 def m6_astuple(self):
1831 return tuple(self.m6.tolist())
1833 def discretize_basesource(self, store, target=None):
1834 factor = self.get_factor()
1835 times, amplitudes = self.effective_stf_pre().discretize_t(
1836 store.config.deltat, self.time)
1837 return meta.DiscretizedMTSource(
1838 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor,
1839 **self._dparams_base_repeated(times))
1841 def pyrocko_moment_tensor(self, store=None, target=None):
1842 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
1844 def pyrocko_event(self, store=None, target=None, **kwargs):
1845 mt = self.pyrocko_moment_tensor(store, target)
1846 return Source.pyrocko_event(
1847 self, store, target,
1848 moment_tensor=self.pyrocko_moment_tensor(store, target),
1849 magnitude=float(mt.moment_magnitude()),
1850 **kwargs)
1853class VLVDSource(SourceWithDerivedMagnitude):
1854 '''
1855 Volumetric linear vector dipole source.
1857 This source is a parameterization for a restricted moment tensor point
1858 source, useful to represent dyke or sill like inflation or deflation
1859 sources. The restriction is such that the moment tensor is rotational
1860 symmetric. It can be represented by a superposition of a linear vector
1861 dipole (here we use a CLVD for convenience) and an isotropic component. The
1862 restricted moment tensor has 4 degrees of freedom: 2 independent
1863 eigenvalues and 2 rotation angles orienting the the symmetry axis.
1865 In this parameterization, the isotropic component is controlled by
1866 ``volume_change``. To define the moment tensor, it must be converted to the
1867 scalar moment of the the MT's isotropic component. For the conversion, the
1868 shear modulus at the source's position must be known. This value is
1869 extracted from the earth model defined in the GF store in use.
1871 The CLVD part by controlled by its scalar moment :math:`M_0`:
1872 ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a
1873 positiv or negativ CLVD (the sign of the largest eigenvalue).
1874 '''
1876 discretized_source_class = meta.DiscretizedMTSource
1878 azimuth = Float.T(
1879 default=0.0,
1880 help='azimuth direction of symmetry axis, clockwise from north [deg].')
1882 dip = Float.T(
1883 default=90.,
1884 help='dip direction of symmetry axis, downward from horizontal [deg].')
1886 volume_change = Float.T(
1887 default=0.,
1888 help='volume change of the inflation/deflation [m^3].')
1890 clvd_moment = Float.T(
1891 default=0.,
1892 help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign '
1893 'controls the sign of the CLVD (the sign of its largest '
1894 'eigenvalue).')
1896 def get_moment_to_volume_change_ratio(self, store, target):
1897 if store is None or target is None:
1898 raise DerivedMagnitudeError(
1899 'Need earth model to convert between volume change and '
1900 'magnitude.')
1902 points = num.array(
1903 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1905 try:
1906 shear_moduli = store.config.get_shear_moduli(
1907 self.lat, self.lon,
1908 points=points,
1909 interpolation=target.interpolation)[0]
1911 bulk_moduli = store.config.get_bulk_moduli(
1912 self.lat, self.lon,
1913 points=points,
1914 interpolation=target.interpolation)[0]
1915 except meta.OutOfBounds:
1916 raise DerivedMagnitudeError(
1917 'Could not get shear modulus at source position.')
1919 return float(2. * shear_moduli + bulk_moduli)
1921 def base_key(self):
1922 return Source.base_key(self) + \
1923 (self.azimuth, self.dip, self.volume_change, self.clvd_moment)
1925 def get_magnitude(self, store=None, target=None):
1926 mt = self.pyrocko_moment_tensor(store, target)
1927 return float(pmt.moment_to_magnitude(mt.moment))
1929 def get_m6(self, store, target):
1930 a = math.sqrt(4. / 3.) * self.clvd_moment
1931 m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1933 rotmat1 = pmt.euler_to_matrix(
1934 d2r * (self.dip - 90.),
1935 d2r * (self.azimuth - 90.),
1936 0.)
1937 m_clvd = num.dot(rotmat1.T, num.dot(m_clvd, rotmat1))
1939 m_iso = self.volume_change * \
1940 self.get_moment_to_volume_change_ratio(store, target)
1942 m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0.,
1943 0., 0.,) * math.sqrt(2. / 3)
1945 m = pmt.to6(m_clvd) + pmt.to6(m_iso)
1946 return m
1948 def get_moment(self, store=None, target=None):
1949 return float(pmt.magnitude_to_moment(
1950 self.get_magnitude(store, target)))
1952 def get_m6_astuple(self, store, target):
1953 m6 = self.get_m6(store, target)
1954 return tuple(m6.tolist())
1956 def discretize_basesource(self, store, target=None):
1957 times, amplitudes = self.effective_stf_pre().discretize_t(
1958 store.config.deltat, self.time)
1960 m6 = self.get_m6(store, target)
1961 m6 *= amplitudes / self.get_factor()
1963 return meta.DiscretizedMTSource(
1964 m6s=m6[num.newaxis, :],
1965 **self._dparams_base_repeated(times))
1967 def pyrocko_moment_tensor(self, store=None, target=None):
1968 m6_astuple = self.get_m6_astuple(store, target)
1969 return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple))
1972class MTSource(Source):
1973 '''
1974 A moment tensor point source.
1975 '''
1977 discretized_source_class = meta.DiscretizedMTSource
1979 mnn = Float.T(
1980 default=1.,
1981 help='north-north component of moment tensor in [Nm]')
1983 mee = Float.T(
1984 default=1.,
1985 help='east-east component of moment tensor in [Nm]')
1987 mdd = Float.T(
1988 default=1.,
1989 help='down-down component of moment tensor in [Nm]')
1991 mne = Float.T(
1992 default=0.,
1993 help='north-east component of moment tensor in [Nm]')
1995 mnd = Float.T(
1996 default=0.,
1997 help='north-down component of moment tensor in [Nm]')
1999 med = Float.T(
2000 default=0.,
2001 help='east-down component of moment tensor in [Nm]')
2003 def __init__(self, **kwargs):
2004 if 'm6' in kwargs:
2005 for (k, v) in zip('mnn mee mdd mne mnd med'.split(),
2006 kwargs.pop('m6')):
2007 kwargs[k] = float(v)
2009 Source.__init__(self, **kwargs)
2011 @property
2012 def m6(self):
2013 return num.array(self.m6_astuple)
2015 @property
2016 def m6_astuple(self):
2017 return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)
2019 @m6.setter
2020 def m6(self, value):
2021 self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value
2023 def base_key(self):
2024 return Source.base_key(self) + self.m6_astuple
2026 def discretize_basesource(self, store, target=None):
2027 times, amplitudes = self.effective_stf_pre().discretize_t(
2028 store.config.deltat, self.time)
2029 return meta.DiscretizedMTSource(
2030 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis],
2031 **self._dparams_base_repeated(times))
2033 def get_magnitude(self, store=None, target=None):
2034 m6 = self.m6
2035 return pmt.moment_to_magnitude(
2036 math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) /
2037 math.sqrt(2.))
2039 def pyrocko_moment_tensor(self, store=None, target=None):
2040 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
2042 def pyrocko_event(self, store=None, target=None, **kwargs):
2043 mt = self.pyrocko_moment_tensor(store, target)
2044 return Source.pyrocko_event(
2045 self, store, target,
2046 moment_tensor=self.pyrocko_moment_tensor(store, target),
2047 magnitude=float(mt.moment_magnitude()),
2048 **kwargs)
2050 @classmethod
2051 def from_pyrocko_event(cls, ev, **kwargs):
2052 d = {}
2053 mt = ev.moment_tensor
2054 if mt:
2055 d.update(m6=tuple(map(float, mt.m6())))
2056 else:
2057 if ev.magnitude is not None:
2058 mom = pmt.magnitude_to_moment(ev.magnitude)
2059 v = math.sqrt(2. / 3.) * mom
2060 d.update(m6=(v, v, v, 0., 0., 0.))
2062 d.update(kwargs)
2063 return super(MTSource, cls).from_pyrocko_event(ev, **d)
2066map_anchor = {
2067 'center': (0.0, 0.0),
2068 'center_left': (-1.0, 0.0),
2069 'center_right': (1.0, 0.0),
2070 'top': (0.0, -1.0),
2071 'top_left': (-1.0, -1.0),
2072 'top_right': (1.0, -1.0),
2073 'bottom': (0.0, 1.0),
2074 'bottom_left': (-1.0, 1.0),
2075 'bottom_right': (1.0, 1.0)}
2078class RectangularSource(SourceWithDerivedMagnitude):
2079 '''
2080 Classical Haskell source model modified for bilateral rupture.
2081 '''
2083 discretized_source_class = meta.DiscretizedMTSource
2085 magnitude = Float.T(
2086 optional=True,
2087 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
2089 strike = Float.T(
2090 default=0.0,
2091 help='strike direction in [deg], measured clockwise from north')
2093 dip = Float.T(
2094 default=90.0,
2095 help='dip angle in [deg], measured downward from horizontal')
2097 rake = Float.T(
2098 default=0.0,
2099 help='rake angle in [deg], '
2100 'measured counter-clockwise from right-horizontal '
2101 'in on-plane view')
2103 length = Float.T(
2104 default=0.,
2105 help='length of rectangular source area [m]')
2107 width = Float.T(
2108 default=0.,
2109 help='width of rectangular source area [m]')
2111 anchor = StringChoice.T(
2112 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2113 'bottom_left', 'bottom_right'],
2114 default='center',
2115 optional=True,
2116 help='Anchor point for positioning the plane, can be: ``top, center '
2117 'bottom, top_left, top_right,bottom_left,'
2118 'bottom_right, center_left, center right``.')
2120 nucleation_x = Float.T(
2121 optional=True,
2122 help='horizontal position of rupture nucleation in normalized fault '
2123 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge)')
2125 nucleation_y = Float.T(
2126 optional=True,
2127 help='down-dip position of rupture nucleation in normalized fault '
2128 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge)')
2130 velocity = Float.T(
2131 default=3500.,
2132 help='speed of rupture front [m/s]')
2134 slip = Float.T(
2135 optional=True,
2136 help='Slip on the rectangular source area [m]')
2138 opening_fraction = Float.T(
2139 default=0.,
2140 help='Determines fraction of slip related to opening. '
2141 '(``-1``: pure tensile closing, '
2142 '``0``: pure shear, '
2143 '``1``: pure tensile opening)')
2145 decimation_factor = Int.T(
2146 optional=True,
2147 default=1,
2148 help='Sub-source decimation factor, a larger decimation will'
2149 ' make the result inaccurate but shorten the necessary'
2150 ' computation time (use for testing puposes only).')
2152 aggressive_oversampling = Bool.T(
2153 default=False,
2154 help='Aggressive oversampling for basesource discretization. '
2155 "When using 'multilinear' interpolation oversampling has"
2156 ' practically no effect.')
2158 def __init__(self, **kwargs):
2159 if 'moment' in kwargs:
2160 mom = kwargs.pop('moment')
2161 if 'magnitude' not in kwargs:
2162 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
2164 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2166 def base_key(self):
2167 return SourceWithDerivedMagnitude.base_key(self) + (
2168 self.magnitude,
2169 self.slip,
2170 self.strike,
2171 self.dip,
2172 self.rake,
2173 self.length,
2174 self.width,
2175 self.nucleation_x,
2176 self.nucleation_y,
2177 self.velocity,
2178 self.decimation_factor,
2179 self.anchor)
2181 def check_conflicts(self):
2182 if self.magnitude is not None and self.slip is not None:
2183 raise DerivedMagnitudeError(
2184 'Magnitude and slip are both defined.')
2186 def get_magnitude(self, store=None, target=None):
2187 self.check_conflicts()
2188 if self.magnitude is not None:
2189 return self.magnitude
2191 elif self.slip is not None:
2192 if None in (store, target):
2193 raise DerivedMagnitudeError(
2194 'Magnitude for a rectangular source with slip defined '
2195 'can only be derived when earth model and target '
2196 'interpolation method are available.')
2198 amplitudes = self._discretize(store, target)[2]
2199 if amplitudes.ndim == 2:
2200 # CLVD component has no net moment, leave out
2201 return float(pmt.moment_to_magnitude(
2202 num.sum(num.abs(amplitudes[0:2, :]).sum())))
2203 else:
2204 return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
2206 else:
2207 return float(pmt.moment_to_magnitude(1.0))
2209 def get_factor(self):
2210 return 1.0
2212 def get_slip_tensile(self):
2213 return self.slip * self.opening_fraction
2215 def get_slip_shear(self):
2216 return self.slip - abs(self.get_slip_tensile)
2218 def _discretize(self, store, target):
2219 if self.nucleation_x is not None:
2220 nucx = self.nucleation_x * 0.5 * self.length
2221 else:
2222 nucx = None
2224 if self.nucleation_y is not None:
2225 nucy = self.nucleation_y * 0.5 * self.width
2226 else:
2227 nucy = None
2229 stf = self.effective_stf_pre()
2231 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
2232 store.config.deltas, store.config.deltat,
2233 self.time, self.north_shift, self.east_shift, self.depth,
2234 self.strike, self.dip, self.length, self.width, self.anchor,
2235 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy,
2236 decimation_factor=self.decimation_factor,
2237 aggressive_oversampling=self.aggressive_oversampling)
2239 if self.slip is not None:
2240 if target is not None:
2241 interpolation = target.interpolation
2242 else:
2243 interpolation = 'nearest_neighbor'
2244 logger.warning(
2245 'no target information available, will use '
2246 '"nearest_neighbor" interpolation when extracting shear '
2247 'modulus from earth model')
2249 shear_moduli = store.config.get_shear_moduli(
2250 self.lat, self.lon,
2251 points=points,
2252 interpolation=interpolation)
2254 tensile_slip = self.get_slip_tensile()
2255 shear_slip = self.slip - abs(tensile_slip)
2257 amplitudes_total = [shear_moduli * shear_slip]
2258 if tensile_slip != 0:
2259 bulk_moduli = store.config.get_bulk_moduli(
2260 self.lat, self.lon,
2261 points=points,
2262 interpolation=interpolation)
2264 tensile_iso = bulk_moduli * tensile_slip
2265 tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip
2267 amplitudes_total.extend([tensile_iso, tensile_clvd])
2269 amplitudes_total = num.vstack(amplitudes_total).squeeze() * \
2270 amplitudes * dl * dw
2272 else:
2273 # normalization to retain total moment
2274 amplitudes_norm = amplitudes / num.sum(amplitudes)
2275 moment = self.get_moment(store, target)
2277 amplitudes_total = [
2278 amplitudes_norm * moment * (1 - abs(self.opening_fraction))]
2279 if self.opening_fraction != 0.:
2280 amplitudes_total.append(
2281 amplitudes_norm * self.opening_fraction * moment)
2283 amplitudes_total = num.vstack(amplitudes_total).squeeze()
2285 return points, times, num.atleast_1d(amplitudes_total), dl, dw, nl, nw
2287 def discretize_basesource(self, store, target=None):
2289 points, times, amplitudes, dl, dw, nl, nw = self._discretize(
2290 store, target)
2292 mot = pmt.MomentTensor(
2293 strike=self.strike, dip=self.dip, rake=self.rake)
2294 m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0)
2296 if amplitudes.ndim == 1:
2297 m6s[:, :] *= amplitudes[:, num.newaxis]
2298 elif amplitudes.ndim == 2:
2299 # shear MT components
2300 rotmat1 = pmt.euler_to_matrix(
2301 d2r * self.dip, d2r * self.strike, d2r * -self.rake)
2302 m6s[:, :] *= amplitudes[0, :][:, num.newaxis]
2304 if amplitudes.shape[0] == 2:
2305 # tensile MT components - moment/magnitude input
2306 tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.)
2307 rot_tensile = pmt.to6(
2308 num.dot(rotmat1.T, num.dot(tensile, rotmat1)))
2310 m6s_tensile = rot_tensile[
2311 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2312 m6s += m6s_tensile
2314 elif amplitudes.shape[0] == 3:
2315 # tensile MT components - slip input
2316 iso = pmt.symmat6(1., 1., 1., 0., 0., 0.)
2317 clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.)
2319 rot_iso = pmt.to6(
2320 num.dot(rotmat1.T, num.dot(iso, rotmat1)))
2321 rot_clvd = pmt.to6(
2322 num.dot(rotmat1.T, num.dot(clvd, rotmat1)))
2324 m6s_iso = rot_iso[
2325 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2326 m6s_clvd = rot_clvd[
2327 num.newaxis, :] * amplitudes[2, :][:, num.newaxis]
2328 m6s += m6s_iso + m6s_clvd
2329 else:
2330 raise ValueError('Unknwown amplitudes shape!')
2331 else:
2332 raise ValueError(
2333 'Unexpected dimension of {}'.format(amplitudes.ndim))
2335 ds = meta.DiscretizedMTSource(
2336 lat=self.lat,
2337 lon=self.lon,
2338 times=times,
2339 north_shifts=points[:, 0],
2340 east_shifts=points[:, 1],
2341 depths=points[:, 2],
2342 m6s=m6s,
2343 dl=dl,
2344 dw=dw,
2345 nl=nl,
2346 nw=nw)
2348 return ds
2350 def xy_to_coord(self, x, y, cs='xyz'):
2351 ln, wd = self.length, self.width
2352 strike, dip = self.strike, self.dip
2354 def array_check(variable):
2355 if not isinstance(variable, num.ndarray):
2356 return num.array(variable)
2357 else:
2358 return variable
2360 x, y = array_check(x), array_check(y)
2362 if x.shape[0] != y.shape[0]:
2363 raise ValueError('Shapes of x and y mismatch')
2365 x, y = x * 0.5 * ln, y * 0.5 * wd
2367 points = num.hstack((
2368 x.reshape(-1, 1), y.reshape(-1, 1), num.zeros((x.shape[0], 1))))
2370 anch_x, anch_y = map_anchor[self.anchor]
2371 points[:, 0] -= anch_x * 0.5 * ln
2372 points[:, 1] -= anch_y * 0.5 * wd
2374 rotmat = num.asarray(
2375 pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
2377 points_rot = num.dot(rotmat.T, points.T).T
2379 points_rot[:, 0] += self.north_shift
2380 points_rot[:, 1] += self.east_shift
2381 points_rot[:, 2] += self.depth
2383 if cs == 'xyz':
2384 return points_rot
2385 elif cs == 'xy':
2386 return points_rot[:, :2]
2387 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2388 latlon = ne_to_latlon(
2389 self.lat, self.lon, points_rot[:, 0], points_rot[:, 1])
2390 latlon = num.array(latlon).T
2391 if cs == 'latlon':
2392 return latlon
2393 elif cs == 'lonlat':
2394 return latlon[:, ::-1]
2395 else:
2396 return num.concatenate(
2397 (latlon, points_rot[:, 2].reshape((len(points_rot), 1))),
2398 axis=1)
2400 def outline(self, cs='xyz'):
2401 x = num.array([-1., 1., 1., -1., -1.])
2402 y = num.array([-1., -1., 1., 1., -1.])
2404 return self.xy_to_coord(x, y, cs=cs)
2406 def points_on_source(self, cs='xyz', **kwargs):
2408 points = points_on_rect_source(
2409 self.strike, self.dip, self.length, self.width,
2410 self.anchor, **kwargs)
2412 points[:, 0] += self.north_shift
2413 points[:, 1] += self.east_shift
2414 points[:, 2] += self.depth
2415 if cs == 'xyz':
2416 return points
2417 elif cs == 'xy':
2418 return points[:, :2]
2419 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2420 latlon = ne_to_latlon(
2421 self.lat, self.lon, points[:, 0], points[:, 1])
2423 latlon = num.array(latlon).T
2424 if cs == 'latlon':
2425 return latlon
2426 elif cs == 'lonlat':
2427 return latlon[:, ::-1]
2428 else:
2429 return num.concatenate(
2430 (latlon, points[:, 2].reshape((len(points), 1))),
2431 axis=1)
2433 def geometry(self, *args, **kwargs):
2434 from pyrocko.model import Geometry
2436 ds = self.discretize_basesource(*args, **kwargs)
2437 nx, ny = ds.nl, ds.nw
2439 def patch_outlines_xy(nx, ny):
2440 points = num.zeros((nx * ny, 2))
2441 points[:, 0] = num.tile(num.linspace(-1., 1., nx), ny)
2442 points[:, 1] = num.repeat(num.linspace(-1., 1., ny), nx)
2444 return points
2446 points_ds = patch_outlines_xy(nx + 1, ny + 1)
2447 npoints = (nx + 1) * (ny + 1)
2449 vertices = num.hstack((
2450 num.ones((npoints, 1)) * self.lat,
2451 num.ones((npoints, 1)) * self.lon,
2452 self.xy_to_coord(points_ds[:, 0], points_ds[:, 1], cs='xyz')))
2454 faces = num.array([[
2455 iy * (nx + 1) + ix,
2456 iy * (nx + 1) + ix + 1,
2457 (iy + 1) * (nx + 1) + ix + 1,
2458 (iy + 1) * (nx + 1) + ix,
2459 iy * (nx + 1) + ix]
2460 for iy in range(ny) for ix in range(nx)])
2462 xyz = self.outline('xyz')
2463 latlon = num.ones((5, 2)) * num.array([self.lat, self.lon])
2464 patchverts = num.hstack((latlon, xyz))
2466 geom = Geometry()
2467 geom.setup(vertices, faces)
2468 geom.set_outlines([patchverts])
2470 if self.stf:
2471 geom.times = num.unique(ds.times)
2473 if self.nucleation_x is not None and self.nucleation_y is not None:
2474 geom.add_property('t_arrival', ds.times)
2476 geom.add_property(
2477 'moment', ds.moments().reshape(ds.nl*ds.nw, -1))
2479 geom.add_property(
2480 'slip', num.ones_like(ds.times) * self.slip)
2482 return geom
2484 def get_nucleation_abs_coord(self, cs='xy'):
2486 if self.nucleation_x is None:
2487 return None, None
2489 coords = from_plane_coords(self.strike, self.dip, self.length,
2490 self.width, self.depth, self.nucleation_x,
2491 self.nucleation_y, lat=self.lat,
2492 lon=self.lon, north_shift=self.north_shift,
2493 east_shift=self.east_shift, cs=cs)
2494 return coords
2496 def pyrocko_moment_tensor(self, store=None, target=None):
2497 return pmt.MomentTensor(
2498 strike=self.strike,
2499 dip=self.dip,
2500 rake=self.rake,
2501 scalar_moment=self.get_moment(store, target))
2503 def pyrocko_event(self, store=None, target=None, **kwargs):
2504 return SourceWithDerivedMagnitude.pyrocko_event(
2505 self, store, target,
2506 **kwargs)
2508 @classmethod
2509 def from_pyrocko_event(cls, ev, **kwargs):
2510 d = {}
2511 mt = ev.moment_tensor
2512 if mt:
2513 (strike, dip, rake), _ = mt.both_strike_dip_rake()
2514 d.update(
2515 strike=float(strike),
2516 dip=float(dip),
2517 rake=float(rake),
2518 magnitude=float(mt.moment_magnitude()))
2520 d.update(kwargs)
2521 return super(RectangularSource, cls).from_pyrocko_event(ev, **d)
2524class PseudoDynamicRupture(SourceWithDerivedMagnitude):
2525 '''
2526 Combined Eikonal and Okada quasi-dynamic rupture model.
2528 Details are described in :doc:`/topics/pseudo-dynamic-rupture`.
2529 Note: attribute `stf` is not used so far, but kept for future applications.
2530 '''
2532 discretized_source_class = meta.DiscretizedMTSource
2534 strike = Float.T(
2535 default=0.0,
2536 help='Strike direction in [deg], measured clockwise from north.')
2538 dip = Float.T(
2539 default=0.0,
2540 help='Dip angle in [deg], measured downward from horizontal.')
2542 length = Float.T(
2543 default=10. * km,
2544 help='Length of rectangular source area in [m].')
2546 width = Float.T(
2547 default=5. * km,
2548 help='Width of rectangular source area in [m].')
2550 anchor = StringChoice.T(
2551 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2552 'bottom_left', 'bottom_right'],
2553 default='center',
2554 optional=True,
2555 help='Anchor point for positioning the plane, can be: ``top, center, '
2556 'bottom, top_left, top_right, bottom_left, '
2557 'bottom_right, center_left, center_right``.')
2559 nucleation_x__ = Array.T(
2560 default=num.array([0.]),
2561 dtype=num.float64,
2562 serialize_as='list',
2563 help='Horizontal position of rupture nucleation in normalized fault '
2564 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge).')
2566 nucleation_y__ = Array.T(
2567 default=num.array([0.]),
2568 dtype=num.float64,
2569 serialize_as='list',
2570 help='Down-dip position of rupture nucleation in normalized fault '
2571 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge).')
2573 nucleation_time__ = Array.T(
2574 optional=True,
2575 help='Time in [s] after origin, when nucleation points defined by '
2576 '``nucleation_x`` and ``nucleation_y`` rupture.',
2577 dtype=num.float64,
2578 serialize_as='list')
2580 gamma = Float.T(
2581 default=0.8,
2582 help='Scaling factor between rupture velocity and S-wave velocity: '
2583 r':math:`v_r = \gamma * v_s`.')
2585 nx = Int.T(
2586 default=2,
2587 help='Number of discrete source patches in x direction (along '
2588 'strike).')
2590 ny = Int.T(
2591 default=2,
2592 help='Number of discrete source patches in y direction (down dip).')
2594 slip = Float.T(
2595 optional=True,
2596 help='Maximum slip of the rectangular source [m]. '
2597 'Setting the slip the tractions/stress field '
2598 'will be normalized to accomodate the desired maximum slip.')
2600 rake = Float.T(
2601 optional=True,
2602 help='Rake angle in [deg], '
2603 'measured counter-clockwise from right-horizontal '
2604 'in on-plane view. Rake is translated into homogenous tractions '
2605 'in strike and up-dip direction. ``rake`` is mutually exclusive '
2606 'with tractions parameter.')
2608 patches = List.T(
2609 OkadaSource.T(),
2610 optional=True,
2611 help='List of all boundary elements/sub faults/fault patches.')
2613 patch_mask__ = Array.T(
2614 dtype=bool,
2615 serialize_as='list',
2616 shape=(None,),
2617 optional=True,
2618 help='Mask for all boundary elements/sub faults/fault patches. True '
2619 'leaves the patch in the calculation, False excludes the patch.')
2621 tractions = TractionField.T(
2622 optional=True,
2623 help='Traction field the rupture plane is exposed to. See the '
2624 ':py:mod:`pyrocko.gf.tractions` module for more details. '
2625 'If ``tractions=None`` and ``rake`` is given'
2626 ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will'
2627 ' be used.')
2629 coef_mat = Array.T(
2630 optional=True,
2631 help='Coefficient matrix linking traction and dislocation field.',
2632 dtype=num.float64,
2633 shape=(None, None))
2635 eikonal_decimation = Int.T(
2636 optional=True,
2637 default=1,
2638 help='Sub-source eikonal factor, a smaller eikonal factor will'
2639 ' increase the accuracy of rupture front calculation but'
2640 ' increases also the computation time.')
2642 decimation_factor = Int.T(
2643 optional=True,
2644 default=1,
2645 help='Sub-source decimation factor, a larger decimation will'
2646 ' make the result inaccurate but shorten the necessary'
2647 ' computation time (use for testing puposes only).')
2649 nthreads = Int.T(
2650 optional=True,
2651 default=1,
2652 help='Number of threads for Okada forward modelling, '
2653 'matrix inversion and calculation of point subsources. '
2654 'Note: for small/medium matrices 1 thread is most efficient.')
2656 pure_shear = Bool.T(
2657 optional=True,
2658 default=False,
2659 help='Calculate only shear tractions and omit tensile tractions.')
2661 smooth_rupture = Bool.T(
2662 default=True,
2663 help='Smooth the tractions by weighting partially ruptured'
2664 ' fault patches.')
2666 aggressive_oversampling = Bool.T(
2667 default=False,
2668 help='Aggressive oversampling for basesource discretization. '
2669 "When using 'multilinear' interpolation oversampling has"
2670 ' practically no effect.')
2672 def __init__(self, **kwargs):
2673 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2674 self._interpolators = {}
2675 self.check_conflicts()
2677 @property
2678 def nucleation_x(self):
2679 return self.nucleation_x__
2681 @nucleation_x.setter
2682 def nucleation_x(self, nucleation_x):
2683 if isinstance(nucleation_x, list):
2684 nucleation_x = num.array(nucleation_x)
2686 elif not isinstance(
2687 nucleation_x, num.ndarray) and nucleation_x is not None:
2689 nucleation_x = num.array([nucleation_x])
2690 self.nucleation_x__ = nucleation_x
2692 @property
2693 def nucleation_y(self):
2694 return self.nucleation_y__
2696 @nucleation_y.setter
2697 def nucleation_y(self, nucleation_y):
2698 if isinstance(nucleation_y, list):
2699 nucleation_y = num.array(nucleation_y)
2701 elif not isinstance(nucleation_y, num.ndarray) \
2702 and nucleation_y is not None:
2703 nucleation_y = num.array([nucleation_y])
2705 self.nucleation_y__ = nucleation_y
2707 @property
2708 def nucleation(self):
2709 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
2711 if (nucl_x is None) or (nucl_y is None):
2712 return None
2714 assert nucl_x.shape[0] == nucl_y.shape[0]
2716 return num.concatenate(
2717 (nucl_x[:, num.newaxis], nucl_y[:, num.newaxis]), axis=1)
2719 @nucleation.setter
2720 def nucleation(self, nucleation):
2721 if isinstance(nucleation, list):
2722 nucleation = num.array(nucleation)
2724 assert nucleation.shape[1] == 2
2726 self.nucleation_x = nucleation[:, 0]
2727 self.nucleation_y = nucleation[:, 1]
2729 @property
2730 def nucleation_time(self):
2731 return self.nucleation_time__
2733 @nucleation_time.setter
2734 def nucleation_time(self, nucleation_time):
2735 if not isinstance(nucleation_time, num.ndarray) \
2736 and nucleation_time is not None:
2737 nucleation_time = num.array([nucleation_time])
2739 self.nucleation_time__ = nucleation_time
2741 @property
2742 def patch_mask(self):
2743 if (self.patch_mask__ is not None and
2744 self.patch_mask__.shape == (self.nx * self.ny,)):
2746 return self.patch_mask__
2747 else:
2748 return num.ones(self.nx * self.ny, dtype=bool)
2750 @patch_mask.setter
2751 def patch_mask(self, patch_mask):
2752 if isinstance(patch_mask, list):
2753 patch_mask = num.array(patch_mask)
2755 self.patch_mask__ = patch_mask
2757 def get_tractions(self):
2758 '''
2759 Get source traction vectors.
2761 If :py:attr:`rake` is given, unit length directed traction vectors
2762 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are returned,
2763 else the given :py:attr:`tractions` are used.
2765 :returns:
2766 Traction vectors per patch.
2767 :rtype:
2768 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2769 '''
2771 if self.rake is not None:
2772 if num.isnan(self.rake):
2773 raise ValueError('Rake must be a real number, not NaN.')
2775 logger.warning(
2776 'Tractions are derived based on the given source rake.')
2777 tractions = DirectedTractions(rake=self.rake)
2778 else:
2779 tractions = self.tractions
2780 return tractions.get_tractions(self.nx, self.ny, self.patches)
2782 def get_scaled_tractions(self, store):
2783 '''
2784 Get traction vectors rescaled to given slip.
2786 Opposing to :py:meth:`get_tractions` traction vectors
2787 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are rescaled to
2788 the given :py:attr:`slip` before returning. If no :py:attr:`slip` and
2789 :py:attr:`rake` are provided, the given :py:attr:`tractions` are
2790 returned without scaling.
2792 :param store:
2793 Green's function database (needs to cover whole region of the
2794 source).
2795 :type store:
2796 :py:class:`~pyrocko.gf.store.Store`
2798 :returns:
2799 Traction vectors per patch.
2800 :rtype:
2801 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2802 '''
2803 tractions = self.tractions
2804 factor = 1.
2806 if self.rake is not None and self.slip is not None:
2807 if num.isnan(self.rake):
2808 raise ValueError('Rake must be a real number, not NaN.')
2810 self.discretize_patches(store)
2811 slip_0t = max(num.linalg.norm(
2812 self.get_slip(scale_slip=False),
2813 axis=1))
2815 factor = self.slip / slip_0t
2816 tractions = DirectedTractions(rake=self.rake)
2818 return tractions.get_tractions(self.nx, self.ny, self.patches) * factor
2820 def base_key(self):
2821 return SourceWithDerivedMagnitude.base_key(self) + (
2822 self.slip,
2823 self.strike,
2824 self.dip,
2825 self.rake,
2826 self.length,
2827 self.width,
2828 float(self.nucleation_x.mean()),
2829 float(self.nucleation_y.mean()),
2830 self.decimation_factor,
2831 self.anchor,
2832 self.pure_shear,
2833 self.gamma,
2834 tuple(self.patch_mask))
2836 def check_conflicts(self):
2837 if self.tractions and self.rake:
2838 raise AttributeError(
2839 'Tractions and rake are mutually exclusive.')
2840 if self.tractions is None and self.rake is None:
2841 self.rake = 0.
2843 def get_magnitude(self, store=None, target=None):
2844 '''
2845 Get total seismic moment magnitude Mw.
2847 :param store:
2848 GF store to guide the discretization and providing the earthmodel
2849 which is needed to calculate moment from slip.
2850 :type store:
2851 :py:class:`~pyrocko.gf.store.Store`
2853 :param target:
2854 Target, used to get GF interpolation settings.
2855 :type target:
2856 :py:class:`pyrocko.gf.targets.Target`
2858 :returns:
2859 Moment magnitude
2860 :rtype:
2861 float
2862 '''
2863 self.check_conflicts()
2864 if self.slip is not None or self.tractions is not None:
2865 if store is None:
2866 raise DerivedMagnitudeError(
2867 'Magnitude for a rectangular source with slip or '
2868 'tractions defined can only be derived when earth model '
2869 'is set.')
2871 moment_rate, calc_times = self.discretize_basesource(
2872 store, target=target).get_moment_rate(store.config.deltat)
2874 deltat = num.concatenate((
2875 (num.diff(calc_times)[0],),
2876 num.diff(calc_times)))
2878 return float(pmt.moment_to_magnitude(
2879 num.sum(moment_rate * deltat)))
2881 else:
2882 return float(pmt.moment_to_magnitude(1.0))
2884 def get_factor(self):
2885 return 1.0
2887 def outline(self, cs='xyz'):
2888 '''
2889 Get source outline corner coordinates.
2891 :param cs:
2892 :ref:`Output coordinate system <coordinate-system-names>`.
2893 :type cs:
2894 str
2896 :returns:
2897 Corner points in desired coordinate system.
2898 :rtype:
2899 :py:class:`~numpy.ndarray`: ``(5, [2, 3])``.
2900 '''
2901 points = outline_rect_source(self.strike, self.dip, self.length,
2902 self.width, self.anchor)
2904 points[:, 0] += self.north_shift
2905 points[:, 1] += self.east_shift
2906 points[:, 2] += self.depth
2907 if cs == 'xyz':
2908 return points
2909 elif cs == 'xy':
2910 return points[:, :2]
2911 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2912 latlon = ne_to_latlon(
2913 self.lat, self.lon, points[:, 0], points[:, 1])
2915 latlon = num.array(latlon).T
2916 if cs == 'latlon':
2917 return latlon
2918 elif cs == 'lonlat':
2919 return latlon[:, ::-1]
2920 else:
2921 return num.concatenate(
2922 (latlon, points[:, 2].reshape((len(points), 1))),
2923 axis=1)
2925 def points_on_source(self, cs='xyz', **kwargs):
2926 '''
2927 Convert relative plane coordinates to geographical coordinates.
2929 Given x and y coordinates (relative source coordinates between -1.
2930 and 1.) are converted to desired geographical coordinates. Coordinates
2931 need to be given as :py:class:`~numpy.ndarray` arguments ``points_x``
2932 and ``points_y``.
2934 :param cs:
2935 :ref:`Output coordinate system <coordinate-system-names>`.
2936 :type cs:
2937 str
2939 :returns:
2940 Point coordinates in desired coordinate system.
2941 :rtype:
2942 :py:class:`~numpy.ndarray`: ``(n_points, [2, 3])``.
2943 '''
2944 points = points_on_rect_source(
2945 self.strike, self.dip, self.length, self.width,
2946 self.anchor, **kwargs)
2948 points[:, 0] += self.north_shift
2949 points[:, 1] += self.east_shift
2950 points[:, 2] += self.depth
2951 if cs == 'xyz':
2952 return points
2953 elif cs == 'xy':
2954 return points[:, :2]
2955 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2956 latlon = ne_to_latlon(
2957 self.lat, self.lon, points[:, 0], points[:, 1])
2959 latlon = num.array(latlon).T
2960 if cs == 'latlon':
2961 return latlon
2962 elif cs == 'lonlat':
2963 return latlon[:, ::-1]
2964 else:
2965 return num.concatenate(
2966 (latlon, points[:, 2].reshape((len(points), 1))),
2967 axis=1)
2969 def pyrocko_moment_tensor(self, store=None, target=None):
2970 '''
2971 Get overall moment tensor of the rupture.
2973 :param store:
2974 GF store to guide the discretization and providing the earthmodel
2975 which is needed to calculate moment from slip.
2976 :type store:
2977 :py:class:`~pyrocko.gf.store.Store`
2979 :param target:
2980 Target, used to get GF interpolation settings.
2981 :type target:
2982 :py:class:`pyrocko.gf.targets.Target`
2984 :returns:
2985 Moment tensor.
2986 :rtype:
2987 :py:class:`~pyrocko.moment_tensor.MomentTensor`
2988 '''
2989 if store is not None:
2990 if not self.patches:
2991 self.discretize_patches(store)
2993 data = self.get_slip()
2994 else:
2995 data = self.get_tractions()
2997 weights = num.linalg.norm(data, axis=1)
2998 weights /= weights.sum()
3000 rakes = num.arctan2(data[:, 1], data[:, 0]) * r2d
3001 rake = num.average(rakes, weights=weights)
3003 return pmt.MomentTensor(
3004 strike=self.strike,
3005 dip=self.dip,
3006 rake=rake,
3007 scalar_moment=self.get_moment(store, target))
3009 def pyrocko_event(self, store=None, target=None, **kwargs):
3010 return SourceWithDerivedMagnitude.pyrocko_event(
3011 self, store, target,
3012 **kwargs)
3014 @classmethod
3015 def from_pyrocko_event(cls, ev, **kwargs):
3016 d = {}
3017 mt = ev.moment_tensor
3018 if mt:
3019 (strike, dip, rake), _ = mt.both_strike_dip_rake()
3020 d.update(
3021 strike=float(strike),
3022 dip=float(dip),
3023 rake=float(rake))
3025 d.update(kwargs)
3026 return super(PseudoDynamicRupture, cls).from_pyrocko_event(ev, **d)
3028 def _discretize_points(self, store, *args, **kwargs):
3029 '''
3030 Discretize source plane with equal vertical and horizontal spacing.
3032 Additional ``*args`` and ``**kwargs`` are passed to
3033 :py:meth:`points_on_source`.
3035 :param store:
3036 Green's function database (needs to cover whole region of the
3037 source).
3038 :type store:
3039 :py:class:`~pyrocko.gf.store.Store`
3041 :returns:
3042 Number of points in strike and dip direction, distance
3043 between adjacent points, coordinates (latlondepth) and coordinates
3044 (xy on fault) for discrete points.
3045 :rtype:
3046 (int, int, float, :py:class:`~numpy.ndarray`,
3047 :py:class:`~numpy.ndarray`).
3048 '''
3049 anch_x, anch_y = map_anchor[self.anchor]
3051 npoints = int(self.width // km) + 1
3052 points = num.zeros((npoints, 3))
3053 points[:, 1] = num.linspace(-1., 1., npoints)
3054 points[:, 1] = (points[:, 1] - anch_y) * self.width/2
3056 rotmat = pmt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.0)
3057 points = num.dot(rotmat.T, points.T).T
3058 points[:, 2] += self.depth
3060 vs_min = store.config.get_vs(
3061 self.lat, self.lon, points,
3062 interpolation='nearest_neighbor')
3063 vr_min = max(vs_min.min(), .5*km) * self.gamma
3065 oversampling = 10.
3066 delta_l = self.length / (self.nx * oversampling)
3067 delta_w = self.width / (self.ny * oversampling)
3069 delta = self.eikonal_decimation * num.min([
3070 store.config.deltat * vr_min / oversampling,
3071 delta_l, delta_w] + [
3072 deltas for deltas in store.config.deltas])
3074 delta = delta_w / num.ceil(delta_w / delta)
3076 nx = int(num.ceil(self.length / delta)) + 1
3077 ny = int(num.ceil(self.width / delta)) + 1
3079 rem_l = (nx-1)*delta - self.length
3080 lim_x = rem_l / self.length
3082 points_xy = num.zeros((nx * ny, 2))
3083 points_xy[:, 0] = num.repeat(
3084 num.linspace(-1.-lim_x, 1.+lim_x, nx), ny)
3085 points_xy[:, 1] = num.tile(
3086 num.linspace(-1., 1., ny), nx)
3088 points = self.points_on_source(
3089 points_x=points_xy[:, 0],
3090 points_y=points_xy[:, 1],
3091 **kwargs)
3093 return nx, ny, delta, points, points_xy
3095 def _discretize_rupture_v(self, store, interpolation='nearest_neighbor',
3096 points=None):
3097 '''
3098 Get rupture velocity for discrete points on source plane.
3100 :param store:
3101 Green's function database (needs to cover the whole region of the
3102 source)
3103 :type store:
3104 :py:class:`~pyrocko.gf.store.Store`
3106 :param interpolation:
3107 Interpolation method to use (choose between ``'nearest_neighbor'``
3108 and ``'multilinear'``).
3109 :type interpolation:
3110 str
3112 :param points:
3113 Coordinates on fault (-1.:1.) of discrete points.
3114 :type points:
3115 :py:class:`~numpy.ndarray`: ``(n_points, 2)``
3117 :returns:
3118 Rupture velocity assumed as :math:`v_s * \\gamma` for discrete
3119 points.
3120 :rtype:
3121 :py:class:`~numpy.ndarray`: ``(n_points, )``.
3122 '''
3124 if points is None:
3125 _, _, _, points, _ = self._discretize_points(store, cs='xyz')
3127 return store.config.get_vs(
3128 self.lat, self.lon,
3129 points=points,
3130 interpolation=interpolation) * self.gamma
3132 def discretize_time(
3133 self, store, interpolation='nearest_neighbor',
3134 vr=None, times=None, *args, **kwargs):
3135 '''
3136 Get rupture start time for discrete points on source plane.
3138 :param store:
3139 Green's function database (needs to cover whole region of the
3140 source)
3141 :type store:
3142 :py:class:`~pyrocko.gf.store.Store`
3144 :param interpolation:
3145 Interpolation method to use (choose between ``'nearest_neighbor'``
3146 and ``'multilinear'``).
3147 :type interpolation:
3148 str
3150 :param vr:
3151 Array, containing rupture user defined rupture velocity values.
3152 :type vr:
3153 :py:class:`~numpy.ndarray`
3155 :param times:
3156 Array, containing zeros, where rupture is starting, real positive
3157 numbers at later secondary nucleation points and -1, where time
3158 will be calculated. If not given, rupture starts at nucleation_x,
3159 nucleation_y. Times are given for discrete points with equal
3160 horizontal and vertical spacing.
3161 :type times:
3162 :py:class:`~numpy.ndarray`
3164 :returns:
3165 Coordinates (latlondepth), coordinates (xy), rupture velocity,
3166 rupture propagation time of discrete points.
3167 :rtype:
3168 :py:class:`~numpy.ndarray`: ``(n_points, 3)``,
3169 :py:class:`~numpy.ndarray`: ``(n_points, 2)``,
3170 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``,
3171 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``.
3172 '''
3173 nx, ny, delta, points, points_xy = self._discretize_points(
3174 store, cs='xyz')
3176 if vr is None or vr.shape != tuple((nx, ny)):
3177 if vr:
3178 logger.warning(
3179 'Given rupture velocities are not in right shape: '
3180 '(%i, %i), but needed is (%i, %i).', *vr.shape + (nx, ny))
3181 vr = self._discretize_rupture_v(store, interpolation, points)\
3182 .reshape(nx, ny)
3184 if vr.shape != tuple((nx, ny)):
3185 logger.warning(
3186 'Given rupture velocities are not in right shape. Therefore'
3187 ' standard rupture velocity array is used.')
3189 def initialize_times():
3190 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
3192 if nucl_x.shape != nucl_y.shape:
3193 raise ValueError(
3194 'Nucleation coordinates have different shape.')
3196 dist_points = num.array([
3197 num.linalg.norm(points_xy - num.array([x, y]).ravel(), axis=1)
3198 for x, y in zip(nucl_x, nucl_y)])
3199 nucl_indices = num.argmin(dist_points, axis=1)
3201 if self.nucleation_time is None:
3202 nucl_times = num.zeros_like(nucl_indices)
3203 else:
3204 if self.nucleation_time.shape == nucl_x.shape:
3205 nucl_times = self.nucleation_time
3206 else:
3207 raise ValueError(
3208 'Nucleation coordinates and times have different '
3209 'shapes')
3211 t = num.full(nx * ny, -1.)
3212 t[nucl_indices] = nucl_times
3213 return t.reshape(nx, ny)
3215 if times is None:
3216 times = initialize_times()
3217 elif times.shape != tuple((nx, ny)):
3218 times = initialize_times()
3219 logger.warning(
3220 'Given times are not in right shape. Therefore standard time'
3221 ' array is used.')
3223 eikonal_ext.eikonal_solver_fmm_cartesian(
3224 speeds=vr, times=times, delta=delta)
3226 return points, points_xy, vr, times
3228 def get_vr_time_interpolators(
3229 self, store, interpolation='nearest_neighbor', force=False,
3230 **kwargs):
3231 '''
3232 Get interpolators for rupture velocity and rupture time.
3234 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3236 :param store:
3237 Green's function database (needs to cover whole region of the
3238 source).
3239 :type store:
3240 :py:class:`~pyrocko.gf.store.Store`
3242 :param interpolation:
3243 Interpolation method to use (choose between ``'nearest_neighbor'``
3244 and ``'multilinear'``).
3245 :type interpolation:
3246 str
3248 :param force:
3249 Force recalculation of the interpolators (e.g. after change of
3250 nucleation point locations/times). Default is ``False``.
3251 :type force:
3252 bool
3253 '''
3254 interp_map = {'multilinear': 'linear', 'nearest_neighbor': 'nearest'}
3255 if interpolation not in interp_map:
3256 raise TypeError(
3257 'Interpolation method %s not available' % interpolation)
3259 if not self._interpolators.get(interpolation, False) or force:
3260 _, points_xy, vr, times = self.discretize_time(
3261 store, **kwargs)
3263 if self.length <= 0.:
3264 raise ValueError(
3265 'length must be larger then 0. not %g' % self.length)
3267 if self.width <= 0.:
3268 raise ValueError(
3269 'width must be larger then 0. not %g' % self.width)
3271 nx, ny = times.shape
3272 anch_x, anch_y = map_anchor[self.anchor]
3274 points_xy[:, 0] = (points_xy[:, 0] - anch_x) * self.length / 2.
3275 points_xy[:, 1] = (points_xy[:, 1] - anch_y) * self.width / 2.
3277 ascont = num.ascontiguousarray
3279 self._interpolators[interpolation] = (
3280 nx, ny, times, vr,
3281 RegularGridInterpolator(
3282 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3283 times,
3284 method=interp_map[interpolation]),
3285 RegularGridInterpolator(
3286 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3287 vr,
3288 method=interp_map[interpolation]))
3290 return self._interpolators[interpolation]
3292 def discretize_patches(
3293 self, store, interpolation='nearest_neighbor', force=False,
3294 grid_shape=(),
3295 **kwargs):
3296 '''
3297 Get rupture start time and OkadaSource elements for points on rupture.
3299 All source elements and their corresponding center points are
3300 calculated and stored in the :py:attr:`patches` attribute.
3302 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3304 :param store:
3305 Green's function database (needs to cover whole region of the
3306 source).
3307 :type store:
3308 :py:class:`~pyrocko.gf.store.Store`
3310 :param interpolation:
3311 Interpolation method to use (choose between ``'nearest_neighbor'``
3312 and ``'multilinear'``).
3313 :type interpolation:
3314 str
3316 :param force:
3317 Force recalculation of the vr and time interpolators ( e.g. after
3318 change of nucleation point locations/times). Default is ``False``.
3319 :type force:
3320 bool
3322 :param grid_shape:
3323 Desired sub fault patch grid size (nlength, nwidth). Either factor
3324 or grid_shape should be set.
3325 :type grid_shape:
3326 :py:class:`tuple` of :py:class:`int`
3327 '''
3328 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3329 self.get_vr_time_interpolators(
3330 store,
3331 interpolation=interpolation, force=force, **kwargs)
3332 anch_x, anch_y = map_anchor[self.anchor]
3334 al = self.length / 2.
3335 aw = self.width / 2.
3336 al1 = -(al + anch_x * al)
3337 al2 = al - anch_x * al
3338 aw1 = -aw + anch_y * aw
3339 aw2 = aw + anch_y * aw
3340 assert num.abs([al1, al2]).sum() == self.length
3341 assert num.abs([aw1, aw2]).sum() == self.width
3343 def get_lame(*a, **kw):
3344 shear_mod = store.config.get_shear_moduli(*a, **kw)
3345 lamb = store.config.get_vp(*a, **kw)**2 \
3346 * store.config.get_rho(*a, **kw) - 2. * shear_mod
3347 return shear_mod, lamb / (2. * (lamb + shear_mod))
3349 shear_mod, poisson = get_lame(
3350 self.lat, self.lon,
3351 num.array([[self.north_shift, self.east_shift, self.depth]]),
3352 interpolation=interpolation)
3354 okada_src = OkadaSource(
3355 lat=self.lat, lon=self.lon,
3356 strike=self.strike, dip=self.dip,
3357 north_shift=self.north_shift, east_shift=self.east_shift,
3358 depth=self.depth,
3359 al1=al1, al2=al2, aw1=aw1, aw2=aw2,
3360 poisson=poisson.mean(),
3361 shearmod=shear_mod.mean(),
3362 opening=kwargs.get('opening', 0.))
3364 if not (self.nx and self.ny):
3365 if grid_shape:
3366 self.nx, self.ny = grid_shape
3367 else:
3368 self.nx = nx
3369 self.ny = ny
3371 source_disc, source_points = okada_src.discretize(self.nx, self.ny)
3373 shear_mod, poisson = get_lame(
3374 self.lat, self.lon,
3375 num.array([src.source_patch()[:3] for src in source_disc]),
3376 interpolation=interpolation)
3378 if (self.nx, self.ny) != (nx, ny):
3379 times_interp = time_interpolator(
3380 num.ascontiguousarray(source_points[:, :2]))
3381 vr_interp = vr_interpolator(
3382 num.ascontiguousarray(source_points[:, :2]))
3383 else:
3384 times_interp = times.T.ravel()
3385 vr_interp = vr.T.ravel()
3387 for isrc, src in enumerate(source_disc):
3388 src.vr = vr_interp[isrc]
3389 src.time = times_interp[isrc] + self.time
3391 self.patches = source_disc
3393 def discretize_basesource(self, store, target=None):
3394 '''
3395 Prepare source for synthetic waveform calculation.
3397 :param store:
3398 Green's function database (needs to cover whole region of the
3399 source).
3400 :type store:
3401 :py:class:`~pyrocko.gf.store.Store`
3403 :param target:
3404 Target information.
3405 :type target:
3406 :py:class:`~pyrocko.gf.targets.Target`
3408 :returns:
3409 Source discretized by a set of moment tensors and times.
3410 :rtype:
3411 :py:class:`~pyrocko.gf.meta.DiscretizedMTSource`
3412 '''
3413 if not target:
3414 interpolation = 'nearest_neighbor'
3415 else:
3416 interpolation = target.interpolation
3418 if not self.patches:
3419 self.discretize_patches(store, interpolation)
3421 if self.coef_mat is None:
3422 self.calc_coef_mat()
3424 delta_slip, slip_times = self.get_delta_slip(store)
3425 npatches = self.nx * self.ny
3426 ntimes = slip_times.size
3428 anch_x, anch_y = map_anchor[self.anchor]
3430 pln = self.length / self.nx
3431 pwd = self.width / self.ny
3433 patch_coords = num.array([
3434 (p.ix, p.iy)
3435 for p in self.patches]).reshape(self.nx, self.ny, 2)
3437 # boundary condition is zero-slip
3438 # is not valid to avoid unwished interpolation effects
3439 slip_grid = num.zeros((self.nx + 2, self.ny + 2, ntimes, 3))
3440 slip_grid[1:-1, 1:-1, :, :] = \
3441 delta_slip.reshape(self.nx, self.ny, ntimes, 3)
3443 slip_grid[0, 0, :, :] = slip_grid[1, 1, :, :]
3444 slip_grid[0, -1, :, :] = slip_grid[1, -2, :, :]
3445 slip_grid[-1, 0, :, :] = slip_grid[-2, 1, :, :]
3446 slip_grid[-1, -1, :, :] = slip_grid[-2, -2, :, :]
3448 slip_grid[1:-1, 0, :, :] = slip_grid[1:-1, 1, :, :]
3449 slip_grid[1:-1, -1, :, :] = slip_grid[1:-1, -2, :, :]
3450 slip_grid[0, 1:-1, :, :] = slip_grid[1, 1:-1, :, :]
3451 slip_grid[-1, 1:-1, :, :] = slip_grid[-2, 1:-1, :, :]
3453 def make_grid(patch_parameter):
3454 grid = num.zeros((self.nx + 2, self.ny + 2))
3455 grid[1:-1, 1:-1] = patch_parameter.reshape(self.nx, self.ny)
3457 grid[0, 0] = grid[1, 1]
3458 grid[0, -1] = grid[1, -2]
3459 grid[-1, 0] = grid[-2, 1]
3460 grid[-1, -1] = grid[-2, -2]
3462 grid[1:-1, 0] = grid[1:-1, 1]
3463 grid[1:-1, -1] = grid[1:-1, -2]
3464 grid[0, 1:-1] = grid[1, 1:-1]
3465 grid[-1, 1:-1] = grid[-2, 1:-1]
3467 return grid
3469 lamb = self.get_patch_attribute('lamb')
3470 mu = self.get_patch_attribute('shearmod')
3472 lamb_grid = make_grid(lamb)
3473 mu_grid = make_grid(mu)
3475 coords_x = num.zeros(self.nx + 2)
3476 coords_x[1:-1] = patch_coords[:, 0, 0]
3477 coords_x[0] = coords_x[1] - pln / 2
3478 coords_x[-1] = coords_x[-2] + pln / 2
3480 coords_y = num.zeros(self.ny + 2)
3481 coords_y[1:-1] = patch_coords[0, :, 1]
3482 coords_y[0] = coords_y[1] - pwd / 2
3483 coords_y[-1] = coords_y[-2] + pwd / 2
3485 slip_interp = RegularGridInterpolator(
3486 (coords_x, coords_y, slip_times),
3487 slip_grid, method='nearest')
3489 lamb_interp = RegularGridInterpolator(
3490 (coords_x, coords_y),
3491 lamb_grid, method='nearest')
3493 mu_interp = RegularGridInterpolator(
3494 (coords_x, coords_y),
3495 mu_grid, method='nearest')
3497 # discretize basesources
3498 mindeltagf = min(tuple(
3499 (self.length / self.nx, self.width / self.ny) +
3500 tuple(store.config.deltas)))
3502 nl = int((1. / self.decimation_factor) *
3503 num.ceil(pln / mindeltagf)) + 1
3504 nw = int((1. / self.decimation_factor) *
3505 num.ceil(pwd / mindeltagf)) + 1
3506 nsrc_patch = int(nl * nw)
3507 dl = pln / nl
3508 dw = pwd / nw
3510 patch_area = dl * dw
3512 xl = num.linspace(-0.5 * (pln - dl), 0.5 * (pln - dl), nl)
3513 xw = num.linspace(-0.5 * (pwd - dw), 0.5 * (pwd - dw), nw)
3515 base_coords = num.zeros((nsrc_patch, 3))
3516 base_coords[:, 0] = num.tile(xl, nw)
3517 base_coords[:, 1] = num.repeat(xw, nl)
3518 base_coords = num.tile(base_coords, (npatches, 1))
3520 center_coords = num.zeros((npatches, 3))
3521 center_coords[:, 0] = num.repeat(
3522 num.arange(self.nx) * pln + pln / 2, self.ny) - self.length / 2
3523 center_coords[:, 1] = num.tile(
3524 num.arange(self.ny) * pwd + pwd / 2, self.nx) - self.width / 2
3526 base_coords -= center_coords.repeat(nsrc_patch, axis=0)
3527 nbaselocs = base_coords.shape[0]
3529 base_interp = base_coords.repeat(ntimes, axis=0)
3531 base_times = num.tile(slip_times, nbaselocs)
3532 base_interp[:, 0] -= anch_x * self.length / 2
3533 base_interp[:, 1] -= anch_y * self.width / 2
3534 base_interp[:, 2] = base_times
3536 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3537 store, interpolation=interpolation)
3539 time_eikonal_max = time_interpolator.values.max()
3541 nbasesrcs = base_interp.shape[0]
3542 delta_slip = slip_interp(base_interp).reshape(nbaselocs, ntimes, 3)
3543 lamb = lamb_interp(base_interp[:, :2]).ravel()
3544 mu = mu_interp(base_interp[:, :2]).ravel()
3546 if False:
3547 try:
3548 import matplotlib.pyplot as plt
3549 coords = base_coords.copy()
3550 norm = num.sum(num.linalg.norm(delta_slip, axis=2), axis=1)
3551 plt.scatter(coords[:, 0], coords[:, 1], c=norm)
3552 plt.show()
3553 except AttributeError:
3554 pass
3556 base_interp[:, 2] = 0.
3557 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
3558 base_interp = num.dot(rotmat.T, base_interp.T).T
3559 base_interp[:, 0] += self.north_shift
3560 base_interp[:, 1] += self.east_shift
3561 base_interp[:, 2] += self.depth
3563 slip_strike = delta_slip[:, :, 0].ravel()
3564 slip_dip = delta_slip[:, :, 1].ravel()
3565 slip_norm = delta_slip[:, :, 2].ravel()
3567 slip_shear = num.linalg.norm([slip_strike, slip_dip], axis=0)
3568 slip_rake = r2d * num.arctan2(slip_dip, slip_strike)
3570 m6s = okada_ext.patch2m6(
3571 strikes=num.full(nbasesrcs, self.strike, dtype=float),
3572 dips=num.full(nbasesrcs, self.dip, dtype=float),
3573 rakes=slip_rake,
3574 disl_shear=slip_shear,
3575 disl_norm=slip_norm,
3576 lamb=lamb,
3577 mu=mu,
3578 nthreads=self.nthreads)
3580 m6s *= patch_area
3582 dl = -self.patches[0].al1 + self.patches[0].al2
3583 dw = -self.patches[0].aw1 + self.patches[0].aw2
3585 base_times[base_times > time_eikonal_max] = time_eikonal_max
3587 ds = meta.DiscretizedMTSource(
3588 lat=self.lat,
3589 lon=self.lon,
3590 times=base_times + self.time,
3591 north_shifts=base_interp[:, 0],
3592 east_shifts=base_interp[:, 1],
3593 depths=base_interp[:, 2],
3594 m6s=m6s,
3595 dl=dl,
3596 dw=dw,
3597 nl=self.nx,
3598 nw=self.ny)
3600 return ds
3602 def calc_coef_mat(self):
3603 '''
3604 Calculate coefficients connecting tractions and dislocations.
3605 '''
3606 if not self.patches:
3607 raise ValueError(
3608 'Patches are needed. Please calculate them first.')
3610 self.coef_mat = make_okada_coefficient_matrix(
3611 self.patches, nthreads=self.nthreads, pure_shear=self.pure_shear)
3613 def get_patch_attribute(self, attr):
3614 '''
3615 Get patch attributes.
3617 :param attr:
3618 Name of selected attribute (see
3619 :py:class`pyrocko.modelling.okada.OkadaSource`).
3620 :type attr:
3621 str
3623 :returns:
3624 Array with attribute value for each fault patch.
3625 :rtype:
3626 :py:class:`~numpy.ndarray`
3628 '''
3629 if not self.patches:
3630 raise ValueError(
3631 'Patches are needed. Please calculate them first.')
3632 return num.array([getattr(p, attr) for p in self.patches])
3634 def get_slip(
3635 self,
3636 time=None,
3637 scale_slip=True,
3638 interpolation='nearest_neighbor',
3639 **kwargs):
3640 '''
3641 Get slip per subfault patch for given time after rupture start.
3643 :param time:
3644 Time after origin [s], for which slip is computed. If not
3645 given, final static slip is returned.
3646 :type time:
3647 float
3649 :param scale_slip:
3650 If ``True`` and :py:attr:`slip` given, all slip values are scaled
3651 to fit the given maximum slip.
3652 :type scale_slip:
3653 bool
3655 :param interpolation:
3656 Interpolation method to use (choose between ``'nearest_neighbor'``
3657 and ``'multilinear'``).
3658 :type interpolation:
3659 str
3661 :returns:
3662 Inverted dislocations (:math:`u_{strike}, u_{dip}, u_{tensile}`)
3663 for each source patch.
3664 :rtype:
3665 :py:class:`~numpy.ndarray`: ``(n_sources, 3)``
3666 '''
3668 if self.patches is None:
3669 raise ValueError(
3670 'Please discretize the source first (discretize_patches())')
3671 npatches = len(self.patches)
3672 tractions = self.get_tractions()
3673 time_patch_max = self.get_patch_attribute('time').max() - self.time
3675 time_patch = time
3676 if time is None:
3677 time_patch = time_patch_max
3679 if self.coef_mat is None:
3680 self.calc_coef_mat()
3682 if tractions.shape != (npatches, 3):
3683 raise AttributeError(
3684 'The traction vector is of invalid shape.'
3685 ' Required shape is (npatches, 3)')
3687 patch_mask = num.ones(npatches, dtype=bool)
3688 if self.patch_mask is not None:
3689 patch_mask = self.patch_mask
3691 times = self.get_patch_attribute('time') - self.time
3692 times[~patch_mask] = time_patch + 1. # exlcude unmasked patches
3693 relevant_sources = num.nonzero(times <= time_patch)[0]
3694 disloc_est = num.zeros_like(tractions)
3696 if self.smooth_rupture:
3697 patch_activation = num.zeros(npatches)
3699 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3700 self.get_vr_time_interpolators(
3701 store, interpolation=interpolation)
3703 # Getting the native Eikonal grid, bit hackish
3704 points_x = num.round(time_interpolator.grid[0], decimals=2)
3705 points_y = num.round(time_interpolator.grid[1], decimals=2)
3706 times_eikonal = time_interpolator.values
3708 time_max = time
3709 if time is None:
3710 time_max = times_eikonal.max()
3712 for ip, p in enumerate(self.patches):
3713 ul = num.round((p.ix + p.al1, p.iy + p.aw1), decimals=2)
3714 lr = num.round((p.ix + p.al2, p.iy + p.aw2), decimals=2)
3716 idx_length = num.logical_and(
3717 points_x >= ul[0], points_x <= lr[0])
3718 idx_width = num.logical_and(
3719 points_y >= ul[1], points_y <= lr[1])
3721 times_patch = times_eikonal[num.ix_(idx_length, idx_width)]
3722 if times_patch.size == 0:
3723 raise AttributeError('could not use smooth_rupture')
3725 patch_activation[ip] = \
3726 (times_patch <= time_max).sum() / times_patch.size
3728 if time_patch == 0 and time_patch != time_patch_max:
3729 patch_activation[ip] = 0.
3731 patch_activation[~patch_mask] = 0. # exlcude unmasked patches
3733 relevant_sources = num.nonzero(patch_activation > 0.)[0]
3735 if relevant_sources.size == 0:
3736 return disloc_est
3738 indices_disl = num.repeat(relevant_sources * 3, 3)
3739 indices_disl[1::3] += 1
3740 indices_disl[2::3] += 2
3742 disloc_est[relevant_sources] = invert_fault_dislocations_bem(
3743 stress_field=tractions[relevant_sources, :].ravel(),
3744 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3745 pure_shear=self.pure_shear, nthreads=self.nthreads,
3746 epsilon=None,
3747 **kwargs)
3749 if self.smooth_rupture:
3750 disloc_est *= patch_activation[:, num.newaxis]
3752 if scale_slip and self.slip is not None:
3753 disloc_tmax = num.zeros(npatches)
3755 indices_disl = num.repeat(num.nonzero(patch_mask)[0] * 3, 3)
3756 indices_disl[1::3] += 1
3757 indices_disl[2::3] += 2
3759 disloc_tmax[patch_mask] = num.linalg.norm(
3760 invert_fault_dislocations_bem(
3761 stress_field=tractions[patch_mask, :].ravel(),
3762 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3763 pure_shear=self.pure_shear, nthreads=self.nthreads,
3764 epsilon=None,
3765 **kwargs), axis=1)
3767 disloc_tmax_max = disloc_tmax.max()
3768 if disloc_tmax_max == 0.:
3769 logger.warning(
3770 'slip scaling not performed. Maximum slip is 0.')
3772 disloc_est *= self.slip / disloc_tmax_max
3774 return disloc_est
3776 def get_delta_slip(
3777 self,
3778 store=None,
3779 deltat=None,
3780 delta=True,
3781 interpolation='nearest_neighbor',
3782 **kwargs):
3783 '''
3784 Get slip change snapshots.
3786 The time interval, within which the slip changes are computed is
3787 determined by the sampling rate of the Green's function database or
3788 ``deltat``. Additional ``**kwargs`` are passed to :py:meth:`get_slip`.
3790 :param store:
3791 Green's function database (needs to cover whole region of of the
3792 source). Its sampling interval is used as time increment for slip
3793 difference calculation. Either ``deltat`` or ``store`` should be
3794 given.
3795 :type store:
3796 :py:class:`~pyrocko.gf.store.Store`
3798 :param deltat:
3799 Time interval for slip difference calculation [s]. Either
3800 ``deltat`` or ``store`` should be given.
3801 :type deltat:
3802 float
3804 :param delta:
3805 If ``True``, slip differences between two time steps are given. If
3806 ``False``, cumulative slip for all time steps.
3807 :type delta:
3808 bool
3810 :param interpolation:
3811 Interpolation method to use (choose between ``'nearest_neighbor'``
3812 and ``'multilinear'``).
3813 :type interpolation:
3814 str
3816 :returns:
3817 Displacement changes(:math:`\\Delta u_{strike},
3818 \\Delta u_{dip} , \\Delta u_{tensile}`) for each source patch and
3819 time; corner times, for which delta slip is computed. The order of
3820 displacement changes array is:
3822 .. math::
3824 &[[\\\\
3825 &[\\Delta u_{strike, patch1, t1},
3826 \\Delta u_{dip, patch1, t1},
3827 \\Delta u_{tensile, patch1, t1}],\\\\
3828 &[\\Delta u_{strike, patch1, t2},
3829 \\Delta u_{dip, patch1, t2},
3830 \\Delta u_{tensile, patch1, t2}]\\\\
3831 &], [\\\\
3832 &[\\Delta u_{strike, patch2, t1}, ...],\\\\
3833 &[\\Delta u_{strike, patch2, t2}, ...]]]\\\\
3835 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3836 :py:class:`~numpy.ndarray`: ``(n_times, )``
3837 '''
3838 if store and deltat:
3839 raise AttributeError(
3840 'Argument collision. '
3841 'Please define only the store or the deltat argument.')
3843 if store:
3844 deltat = store.config.deltat
3846 if not deltat:
3847 raise AttributeError('Please give a GF store or set deltat.')
3849 npatches = len(self.patches)
3851 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3852 store, interpolation=interpolation)
3853 tmax = time_interpolator.values.max()
3855 calc_times = num.arange(0., tmax + deltat, deltat)
3856 calc_times[calc_times > tmax] = tmax
3858 disloc_est = num.zeros((npatches, calc_times.size, 3))
3860 for itime, t in enumerate(calc_times):
3861 disloc_est[:, itime, :] = self.get_slip(
3862 time=t, scale_slip=False, **kwargs)
3864 if self.slip:
3865 disloc_tmax = num.linalg.norm(
3866 self.get_slip(scale_slip=False, time=tmax),
3867 axis=1)
3869 disloc_tmax_max = disloc_tmax.max()
3870 if disloc_tmax_max == 0.:
3871 logger.warning(
3872 'Slip scaling not performed. Maximum slip is 0.')
3873 else:
3874 disloc_est *= self.slip / disloc_tmax_max
3876 if not delta:
3877 return disloc_est, calc_times
3879 # if we have only one timestep there is no gradient
3880 if calc_times.size > 1:
3881 disloc_init = disloc_est[:, 0, :]
3882 disloc_est = num.diff(disloc_est, axis=1)
3883 disloc_est = num.concatenate((
3884 disloc_init[:, num.newaxis, :], disloc_est), axis=1)
3886 calc_times = calc_times
3888 return disloc_est, calc_times
3890 def get_slip_rate(self, *args, **kwargs):
3891 '''
3892 Get slip rate inverted from patches.
3894 The time interval, within which the slip rates are computed is
3895 determined by the sampling rate of the Green's function database or
3896 ``deltat``. Additional ``*args`` and ``**kwargs`` are passed to
3897 :py:meth:`get_delta_slip`.
3899 :returns:
3900 Slip rates (:math:`\\Delta u_{strike}/\\Delta t`,
3901 :math:`\\Delta u_{dip}/\\Delta t, \\Delta u_{tensile}/\\Delta t`)
3902 for each source patch and time; corner times, for which slip rate
3903 is computed. The order of sliprate array is:
3905 .. math::
3907 &[[\\\\
3908 &[\\Delta u_{strike, patch1, t1}/\\Delta t,
3909 \\Delta u_{dip, patch1, t1}/\\Delta t,
3910 \\Delta u_{tensile, patch1, t1}/\\Delta t],\\\\
3911 &[\\Delta u_{strike, patch1, t2}/\\Delta t,
3912 \\Delta u_{dip, patch1, t2}/\\Delta t,
3913 \\Delta u_{tensile, patch1, t2}/\\Delta t]], [\\\\
3914 &[\\Delta u_{strike, patch2, t1}/\\Delta t, ...],\\\\
3915 &[\\Delta u_{strike, patch2, t2}/\\Delta t, ...]]]\\\\
3917 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3918 :py:class:`~numpy.ndarray`: ``(n_times, )``
3919 '''
3920 ddisloc_est, calc_times = self.get_delta_slip(
3921 *args, delta=True, **kwargs)
3923 dt = num.concatenate(
3924 [(num.diff(calc_times)[0], ), num.diff(calc_times)])
3925 slip_rate = num.linalg.norm(ddisloc_est, axis=2) / dt
3927 return slip_rate, calc_times
3929 def get_moment_rate_patches(self, *args, **kwargs):
3930 '''
3931 Get scalar seismic moment rate for each patch individually.
3933 Additional ``*args`` and ``**kwargs`` are passed to
3934 :py:meth:`get_slip_rate`.
3936 :returns:
3937 Seismic moment rate for each source patch and time; corner times,
3938 for which patch moment rate is computed based on slip rate. The
3939 order of the moment rate array is:
3941 .. math::
3943 &[\\\\
3944 &[(\\Delta M / \\Delta t)_{patch1, t1},
3945 (\\Delta M / \\Delta t)_{patch1, t2}, ...],\\\\
3946 &[(\\Delta M / \\Delta t)_{patch2, t1},
3947 (\\Delta M / \\Delta t)_{patch, t2}, ...],\\\\
3948 &[...]]\\\\
3950 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times)``,
3951 :py:class:`~numpy.ndarray`: ``(n_times, )``
3952 '''
3953 slip_rate, calc_times = self.get_slip_rate(*args, **kwargs)
3955 shear_mod = self.get_patch_attribute('shearmod')
3956 p_length = self.get_patch_attribute('length')
3957 p_width = self.get_patch_attribute('width')
3959 dA = p_length * p_width
3961 mom_rate = shear_mod[:, num.newaxis] * slip_rate * dA[:, num.newaxis]
3963 return mom_rate, calc_times
3965 def get_moment_rate(self, store, target=None, deltat=None):
3966 '''
3967 Get seismic source moment rate for the total source (STF).
3969 :param store:
3970 Green's function database (needs to cover whole region of of the
3971 source). Its ``deltat`` [s] is used as time increment for slip
3972 difference calculation. Either ``deltat`` or ``store`` should be
3973 given.
3974 :type store:
3975 :py:class:`~pyrocko.gf.store.Store`
3977 :param target:
3978 Target information, needed for interpolation method.
3979 :type target:
3980 :py:class:`~pyrocko.gf.targets.Target`
3982 :param deltat:
3983 Time increment for slip difference calculation [s]. If not given
3984 ``store.deltat`` is used.
3985 :type deltat:
3986 float
3988 :return:
3989 Seismic moment rate [Nm/s] for each time; corner times, for which
3990 moment rate is computed. The order of the moment rate array is:
3992 .. math::
3994 &[\\\\
3995 &(\\Delta M / \\Delta t)_{t1},\\\\
3996 &(\\Delta M / \\Delta t)_{t2},\\\\
3997 &...]\\\\
3999 :rtype:
4000 :py:class:`~numpy.ndarray`: ``(n_times, )``,
4001 :py:class:`~numpy.ndarray`: ``(n_times, )``
4002 '''
4003 if not deltat:
4004 deltat = store.config.deltat
4005 return self.discretize_basesource(
4006 store, target=target).get_moment_rate(deltat)
4008 def get_moment(self, *args, **kwargs):
4009 '''
4010 Get cumulative seismic moment.
4012 Additional ``*args`` and ``**kwargs`` are passed to
4013 :py:meth:`get_magnitude`.
4015 :returns:
4016 Cumulative seismic moment in [Nm].
4017 :rtype:
4018 float
4019 '''
4020 return float(pmt.magnitude_to_moment(self.get_magnitude(
4021 *args, **kwargs)))
4023 def rescale_slip(self, magnitude=None, moment=None, **kwargs):
4024 '''
4025 Rescale source slip based on given target magnitude or seismic moment.
4027 Rescale the maximum source slip to fit the source moment magnitude or
4028 seismic moment to the given target values. Either ``magnitude`` or
4029 ``moment`` need to be given. Additional ``**kwargs`` are passed to
4030 :py:meth:`get_moment`.
4032 :param magnitude:
4033 Target moment magnitude :math:`M_\\mathrm{w}` as in
4034 [Hanks and Kanamori, 1979]
4035 :type magnitude:
4036 float
4038 :param moment:
4039 Target seismic moment :math:`M_0` [Nm].
4040 :type moment:
4041 float
4042 '''
4043 if self.slip is None:
4044 self.slip = 1.
4045 logger.warning('No slip found for rescaling. '
4046 'An initial slip of 1 m is assumed.')
4048 if magnitude is None and moment is None:
4049 raise ValueError(
4050 'Either target magnitude or moment need to be given.')
4052 moment_init = self.get_moment(**kwargs)
4054 if magnitude is not None:
4055 moment = pmt.magnitude_to_moment(magnitude)
4057 self.slip *= moment / moment_init
4059 def get_centroid(self, store, *args, **kwargs):
4060 '''
4061 Centroid of the pseudo dynamic rupture model.
4063 The centroid location and time are derived from the locations and times
4064 of the individual patches weighted with their moment contribution.
4065 Additional ``**kwargs`` are passed to :py:meth:`pyrocko_moment_tensor`.
4067 :param store:
4068 Green's function database (needs to cover whole region of of the
4069 source). Its ``deltat`` [s] is used as time increment for slip
4070 difference calculation. Either ``deltat`` or ``store`` should be
4071 given.
4072 :type store:
4073 :py:class:`~pyrocko.gf.store.Store`
4075 :returns:
4076 The centroid location and associated moment tensor.
4077 :rtype:
4078 :py:class:`pyrocko.model.event.Event`
4079 '''
4080 _, _, _, _, time, _ = self.get_vr_time_interpolators(store)
4081 t_max = time.values.max()
4083 moment_rate, times = self.get_moment_rate_patches(deltat=t_max)
4085 moment = num.sum(moment_rate * times, axis=1)
4086 weights = moment / moment.sum()
4088 norths = self.get_patch_attribute('north_shift')
4089 easts = self.get_patch_attribute('east_shift')
4090 depths = self.get_patch_attribute('depth')
4092 centroid_n = num.sum(weights * norths)
4093 centroid_e = num.sum(weights * easts)
4094 centroid_d = num.sum(weights * depths)
4096 centroid_lat, centroid_lon = ne_to_latlon(
4097 self.lat, self.lon, centroid_n, centroid_e)
4099 moment_rate_, times = self.get_moment_rate(store)
4100 delta_times = num.concatenate((
4101 [times[1] - times[0]],
4102 num.diff(times)))
4103 moment_src = delta_times * moment_rate
4105 centroid_t = num.sum(
4106 moment_src / num.sum(moment_src) * times) + self.time
4108 mt = self.pyrocko_moment_tensor(store, *args, **kwargs)
4110 return model.Event(
4111 lat=centroid_lat,
4112 lon=centroid_lon,
4113 depth=centroid_d,
4114 time=centroid_t,
4115 moment_tensor=mt,
4116 magnitude=mt.magnitude,
4117 duration=t_max)
4119 def get_coulomb_failure_stress(
4120 self,
4121 receiver_points,
4122 friction,
4123 pressure,
4124 strike,
4125 dip,
4126 rake,
4127 time=None,
4128 *args,
4129 **kwargs):
4130 '''
4131 Calculate Coulomb failure stress change CFS.
4133 The function obtains the Coulomb failure stress change :math:`\\Delta
4134 \\sigma_C` at arbitrary receiver points with a commonly oriented
4135 receiver plane assuming:
4137 .. math::
4139 \\Delta \\sigma_C = \\sigma_S - \\mu (\\sigma_N - \\Delta p)
4141 with the shear stress :math:`\\sigma_S`, the coefficient of friction
4142 :math:`\\mu`, the normal stress :math:`\\sigma_N`, and the pore fluid
4143 pressure change :math:`\\Delta p`. Each receiver point is characterized
4144 by its geographical coordinates, and depth. The required receiver plane
4145 orientation is defined by ``strike``, ``dip``, and ``rake``. The
4146 Coulomb failure stress change is calculated for a given time after
4147 rupture origin time.
4149 :param receiver_points:
4150 Location of the receiver points in Northing, Easting, and depth in
4151 [m].
4152 :type receiver_points:
4153 :py:class:`~numpy.ndarray`: ``(n_receiver, 3)``
4155 :param friction:
4156 Coefficient of friction.
4157 :type friction:
4158 float
4160 :param pressure:
4161 Pore pressure change in [Pa].
4162 :type pressure:
4163 float
4165 :param strike:
4166 Strike of the receiver plane in [deg].
4167 :type strike:
4168 float
4170 :param dip:
4171 Dip of the receiver plane in [deg].
4172 :type dip:
4173 float
4175 :param rake:
4176 Rake of the receiver plane in [deg].
4177 :type rake:
4178 float
4180 :param time:
4181 Time after origin [s], for which the resulting :math:`\\Delta
4182 \\Sigma_c` is computed. If not given, :math:`\\Delta \\Sigma_c` is
4183 derived based on the final static slip.
4184 :type time:
4185 float
4187 :returns:
4188 The Coulomb failure stress change :math:`\\Delta \\Sigma_c` at each
4189 receiver point in [Pa].
4190 :rtype:
4191 :py:class:`~numpy.ndarray`: ``(n_receiver,)``
4192 '''
4193 # dislocation at given time
4194 source_slip = self.get_slip(time=time, scale_slip=True)
4196 # source planes
4197 source_patches = num.array([
4198 src.source_patch() for src in self.patches])
4200 # earth model
4201 lambda_mean = num.mean([src.lamb for src in self.patches])
4202 mu_mean = num.mean([src.shearmod for src in self.patches])
4204 # Dislocation and spatial derivatives from okada in NED
4205 results = okada_ext.okada(
4206 source_patches,
4207 source_slip,
4208 receiver_points,
4209 lambda_mean,
4210 mu_mean,
4211 rotate_sdn=False, # TODO Check
4212 stack_sources=0, # TODO Check
4213 *args, **kwargs)
4215 # resolve stress tensor (sum!)
4216 diag_ind = [0, 4, 8]
4217 kron = num.zeros(9)
4218 kron[diag_ind] = 1.
4219 kron = kron[num.newaxis, num.newaxis, :]
4221 eps = 0.5 * (
4222 results[:, :, 3:] +
4223 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
4225 dilatation \
4226 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
4228 stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps
4230 # superposed stress of all sources at receiver locations
4231 stress_sum = num.sum(stress, axis=0)
4233 # get shear and normal stress from stress tensor
4234 strike_rad = d2r * strike
4235 dip_rad = d2r * dip
4236 rake_rad = d2r * rake
4238 n_rec = receiver_points.shape[0]
4239 stress_normal = num.zeros(n_rec)
4240 tau = num.zeros(n_rec)
4242 # Get vectors in receiver fault normal (ns), strike (rst) and
4243 # dip (rdi) direction
4244 ns = num.zeros(3)
4245 rst = num.zeros(3)
4246 rdi = num.zeros(3)
4248 ns[0] = num.sin(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4249 ns[1] = num.sin(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4250 ns[2] = -num.cos(dip_rad)
4252 rst[0] = num.cos(strike_rad)
4253 rst[1] = num.sin(strike_rad)
4254 rst[2] = 0.0
4256 rdi[0] = num.cos(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4257 rdi[1] = num.cos(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4258 rdi[2] = num.sin(dip_rad)
4260 ts = rst * num.cos(rake_rad) - rdi * num.sin(rake_rad)
4262 stress_normal = num.sum(
4263 num.tile(ns, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4265 tau = num.sum(
4266 num.tile(ts, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4268 # calculate cfs using formula above and return
4269 return tau + friction * (stress_normal + pressure)
4272class DoubleDCSource(SourceWithMagnitude):
4273 '''
4274 Two double-couple point sources separated in space and time.
4275 Moment share between the sub-sources is controlled by the
4276 parameter mix.
4277 The position of the subsources is dependent on the moment
4278 distribution between the two sources. Depth, east and north
4279 shift are given for the centroid between the two double-couples.
4280 The subsources will positioned according to their moment shares
4281 around this centroid position.
4282 This is done according to their delta parameters, which are
4283 therefore in relation to that centroid.
4284 Note that depth of the subsources therefore can be
4285 depth+/-delta_depth. For shallow earthquakes therefore
4286 the depth has to be chosen deeper to avoid sampling
4287 above surface.
4288 '''
4290 strike1 = Float.T(
4291 default=0.0,
4292 help='strike direction in [deg], measured clockwise from north')
4294 dip1 = Float.T(
4295 default=90.0,
4296 help='dip angle in [deg], measured downward from horizontal')
4298 azimuth = Float.T(
4299 default=0.0,
4300 help='azimuth to second double-couple [deg], '
4301 'measured at first, clockwise from north')
4303 rake1 = Float.T(
4304 default=0.0,
4305 help='rake angle in [deg], '
4306 'measured counter-clockwise from right-horizontal '
4307 'in on-plane view')
4309 strike2 = Float.T(
4310 default=0.0,
4311 help='strike direction in [deg], measured clockwise from north')
4313 dip2 = Float.T(
4314 default=90.0,
4315 help='dip angle in [deg], measured downward from horizontal')
4317 rake2 = Float.T(
4318 default=0.0,
4319 help='rake angle in [deg], '
4320 'measured counter-clockwise from right-horizontal '
4321 'in on-plane view')
4323 delta_time = Float.T(
4324 default=0.0,
4325 help='separation of double-couples in time (t2-t1) [s]')
4327 delta_depth = Float.T(
4328 default=0.0,
4329 help='difference in depth (z2-z1) [m]')
4331 distance = Float.T(
4332 default=0.0,
4333 help='distance between the two double-couples [m]')
4335 mix = Float.T(
4336 default=0.5,
4337 help='how to distribute the moment to the two doublecouples '
4338 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1')
4340 stf1 = STF.T(
4341 optional=True,
4342 help='Source time function of subsource 1 '
4343 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4345 stf2 = STF.T(
4346 optional=True,
4347 help='Source time function of subsource 2 '
4348 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4350 discretized_source_class = meta.DiscretizedMTSource
4352 def base_key(self):
4353 return (
4354 self.time, self.depth, self.lat, self.north_shift,
4355 self.lon, self.east_shift, type(self).__name__) + \
4356 self.effective_stf1_pre().base_key() + \
4357 self.effective_stf2_pre().base_key() + (
4358 self.strike1, self.dip1, self.rake1,
4359 self.strike2, self.dip2, self.rake2,
4360 self.delta_time, self.delta_depth,
4361 self.azimuth, self.distance, self.mix)
4363 def get_factor(self):
4364 return self.moment
4366 def effective_stf1_pre(self):
4367 return self.stf1 or self.stf or g_unit_pulse
4369 def effective_stf2_pre(self):
4370 return self.stf2 or self.stf or g_unit_pulse
4372 def effective_stf_post(self):
4373 return g_unit_pulse
4375 def split(self):
4376 a1 = 1.0 - self.mix
4377 a2 = self.mix
4378 delta_north = math.cos(self.azimuth * d2r) * self.distance
4379 delta_east = math.sin(self.azimuth * d2r) * self.distance
4381 dc1 = DCSource(
4382 lat=self.lat,
4383 lon=self.lon,
4384 time=self.time - self.delta_time * a2,
4385 north_shift=self.north_shift - delta_north * a2,
4386 east_shift=self.east_shift - delta_east * a2,
4387 depth=self.depth - self.delta_depth * a2,
4388 moment=self.moment * a1,
4389 strike=self.strike1,
4390 dip=self.dip1,
4391 rake=self.rake1,
4392 stf=self.stf1 or self.stf)
4394 dc2 = DCSource(
4395 lat=self.lat,
4396 lon=self.lon,
4397 time=self.time + self.delta_time * a1,
4398 north_shift=self.north_shift + delta_north * a1,
4399 east_shift=self.east_shift + delta_east * a1,
4400 depth=self.depth + self.delta_depth * a1,
4401 moment=self.moment * a2,
4402 strike=self.strike2,
4403 dip=self.dip2,
4404 rake=self.rake2,
4405 stf=self.stf2 or self.stf)
4407 return [dc1, dc2]
4409 def discretize_basesource(self, store, target=None):
4410 a1 = 1.0 - self.mix
4411 a2 = self.mix
4412 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4413 rake=self.rake1, scalar_moment=a1)
4414 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4415 rake=self.rake2, scalar_moment=a2)
4417 delta_north = math.cos(self.azimuth * d2r) * self.distance
4418 delta_east = math.sin(self.azimuth * d2r) * self.distance
4420 times1, amplitudes1 = self.effective_stf1_pre().discretize_t(
4421 store.config.deltat, self.time - self.delta_time * a2)
4423 times2, amplitudes2 = self.effective_stf2_pre().discretize_t(
4424 store.config.deltat, self.time + self.delta_time * a1)
4426 nt1 = times1.size
4427 nt2 = times2.size
4429 ds = meta.DiscretizedMTSource(
4430 lat=self.lat,
4431 lon=self.lon,
4432 times=num.concatenate((times1, times2)),
4433 north_shifts=num.concatenate((
4434 num.repeat(self.north_shift - delta_north * a2, nt1),
4435 num.repeat(self.north_shift + delta_north * a1, nt2))),
4436 east_shifts=num.concatenate((
4437 num.repeat(self.east_shift - delta_east * a2, nt1),
4438 num.repeat(self.east_shift + delta_east * a1, nt2))),
4439 depths=num.concatenate((
4440 num.repeat(self.depth - self.delta_depth * a2, nt1),
4441 num.repeat(self.depth + self.delta_depth * a1, nt2))),
4442 m6s=num.vstack((
4443 mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis],
4444 mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis])))
4446 return ds
4448 def pyrocko_moment_tensor(self, store=None, target=None):
4449 a1 = 1.0 - self.mix
4450 a2 = self.mix
4451 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4452 rake=self.rake1,
4453 scalar_moment=a1 * self.moment)
4454 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4455 rake=self.rake2,
4456 scalar_moment=a2 * self.moment)
4457 return pmt.MomentTensor(m=mot1.m() + mot2.m())
4459 def pyrocko_event(self, store=None, target=None, **kwargs):
4460 return SourceWithMagnitude.pyrocko_event(
4461 self, store, target,
4462 moment_tensor=self.pyrocko_moment_tensor(store, target),
4463 **kwargs)
4465 @classmethod
4466 def from_pyrocko_event(cls, ev, **kwargs):
4467 d = {}
4468 mt = ev.moment_tensor
4469 if mt:
4470 (strike, dip, rake), _ = mt.both_strike_dip_rake()
4471 d.update(
4472 strike1=float(strike),
4473 dip1=float(dip),
4474 rake1=float(rake),
4475 strike2=float(strike),
4476 dip2=float(dip),
4477 rake2=float(rake),
4478 mix=0.0,
4479 magnitude=float(mt.moment_magnitude()))
4481 d.update(kwargs)
4482 source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d)
4483 source.stf1 = source.stf
4484 source.stf2 = HalfSinusoidSTF(effective_duration=0.)
4485 source.stf = None
4486 return source
4489class RingfaultSource(SourceWithMagnitude):
4490 '''
4491 A ring fault with vertical doublecouples.
4492 '''
4494 diameter = Float.T(
4495 default=1.0,
4496 help='diameter of the ring in [m]')
4498 sign = Float.T(
4499 default=1.0,
4500 help='inside of the ring moves up (+1) or down (-1)')
4502 strike = Float.T(
4503 default=0.0,
4504 help='strike direction of the ring plane, clockwise from north,'
4505 ' in [deg]')
4507 dip = Float.T(
4508 default=0.0,
4509 help='dip angle of the ring plane from horizontal in [deg]')
4511 npointsources = Int.T(
4512 default=360,
4513 help='number of point sources to use')
4515 discretized_source_class = meta.DiscretizedMTSource
4517 def base_key(self):
4518 return Source.base_key(self) + (
4519 self.strike, self.dip, self.diameter, self.npointsources)
4521 def get_factor(self):
4522 return self.sign * self.moment
4524 def discretize_basesource(self, store=None, target=None):
4525 n = self.npointsources
4526 phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False)
4528 points = num.zeros((n, 3))
4529 points[:, 0] = num.cos(phi) * 0.5 * self.diameter
4530 points[:, 1] = num.sin(phi) * 0.5 * self.diameter
4532 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
4533 points = num.dot(rotmat.T, points.T).T # !!! ?
4535 points[:, 0] += self.north_shift
4536 points[:, 1] += self.east_shift
4537 points[:, 2] += self.depth
4539 m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90.,
4540 scalar_moment=1.0 / n).m())
4542 rotmats = num.transpose(
4543 [[num.cos(phi), num.sin(phi), num.zeros(n)],
4544 [-num.sin(phi), num.cos(phi), num.zeros(n)],
4545 [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1))
4547 ms = num.zeros((n, 3, 3))
4548 for i in range(n):
4549 mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i]))
4550 ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat))
4552 m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2],
4553 ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T
4555 times, amplitudes = self.effective_stf_pre().discretize_t(
4556 store.config.deltat, self.time)
4558 nt = times.size
4560 return meta.DiscretizedMTSource(
4561 times=num.tile(times, n),
4562 lat=self.lat,
4563 lon=self.lon,
4564 north_shifts=num.repeat(points[:, 0], nt),
4565 east_shifts=num.repeat(points[:, 1], nt),
4566 depths=num.repeat(points[:, 2], nt),
4567 m6s=num.repeat(m6s, nt, axis=0) * num.tile(
4568 amplitudes, n)[:, num.newaxis])
4571class CombiSource(Source):
4572 '''
4573 Composite source model.
4574 '''
4576 discretized_source_class = meta.DiscretizedMTSource
4578 subsources = List.T(Source.T())
4580 def __init__(self, subsources=[], **kwargs):
4581 if not subsources:
4582 raise BadRequest(
4583 'Need at least one sub-source to create a CombiSource object.')
4585 lats = num.array(
4586 [subsource.lat for subsource in subsources], dtype=float)
4587 lons = num.array(
4588 [subsource.lon for subsource in subsources], dtype=float)
4590 lat, lon = lats[0], lons[0]
4591 if not num.all(lats == lat) and num.all(lons == lon):
4592 subsources = [s.clone() for s in subsources]
4593 for subsource in subsources[1:]:
4594 subsource.set_origin(lat, lon)
4596 depth = float(num.mean([p.depth for p in subsources]))
4597 time = float(num.mean([p.time for p in subsources]))
4598 north_shift = float(num.mean([p.north_shift for p in subsources]))
4599 east_shift = float(num.mean([p.east_shift for p in subsources]))
4600 kwargs.update(
4601 time=time,
4602 lat=float(lat),
4603 lon=float(lon),
4604 north_shift=north_shift,
4605 east_shift=east_shift,
4606 depth=depth)
4608 Source.__init__(self, subsources=subsources, **kwargs)
4610 def get_factor(self):
4611 return 1.0
4613 def discretize_basesource(self, store, target=None):
4614 dsources = []
4615 for sf in self.subsources:
4616 ds = sf.discretize_basesource(store, target)
4617 ds.m6s *= sf.get_factor()
4618 dsources.append(ds)
4620 return meta.DiscretizedMTSource.combine(dsources)
4623class SFSource(Source):
4624 '''
4625 A single force point source.
4627 Supported GF schemes: `'elastic5'`.
4628 '''
4630 discretized_source_class = meta.DiscretizedSFSource
4632 fn = Float.T(
4633 default=0.,
4634 help='northward component of single force [N]')
4636 fe = Float.T(
4637 default=0.,
4638 help='eastward component of single force [N]')
4640 fd = Float.T(
4641 default=0.,
4642 help='downward component of single force [N]')
4644 def __init__(self, **kwargs):
4645 Source.__init__(self, **kwargs)
4647 def base_key(self):
4648 return Source.base_key(self) + (self.fn, self.fe, self.fd)
4650 def get_factor(self):
4651 return 1.0
4653 def discretize_basesource(self, store, target=None):
4654 times, amplitudes = self.effective_stf_pre().discretize_t(
4655 store.config.deltat, self.time)
4656 forces = amplitudes[:, num.newaxis] * num.array(
4657 [[self.fn, self.fe, self.fd]], dtype=float)
4659 return meta.DiscretizedSFSource(forces=forces,
4660 **self._dparams_base_repeated(times))
4662 def pyrocko_event(self, store=None, target=None, **kwargs):
4663 return Source.pyrocko_event(
4664 self, store, target,
4665 **kwargs)
4667 @classmethod
4668 def from_pyrocko_event(cls, ev, **kwargs):
4669 d = {}
4670 d.update(kwargs)
4671 return super(SFSource, cls).from_pyrocko_event(ev, **d)
4674class PorePressurePointSource(Source):
4675 '''
4676 Excess pore pressure point source.
4678 For poro-elastic initial value problem where an excess pore pressure is
4679 brought into a small source volume.
4680 '''
4682 discretized_source_class = meta.DiscretizedPorePressureSource
4684 pp = Float.T(
4685 default=1.0,
4686 help='initial excess pore pressure in [Pa]')
4688 def base_key(self):
4689 return Source.base_key(self)
4691 def get_factor(self):
4692 return self.pp
4694 def discretize_basesource(self, store, target=None):
4695 return meta.DiscretizedPorePressureSource(pp=arr(1.0),
4696 **self._dparams_base())
4699class PorePressureLineSource(Source):
4700 '''
4701 Excess pore pressure line source.
4703 The line source is centered at (north_shift, east_shift, depth).
4704 '''
4706 discretized_source_class = meta.DiscretizedPorePressureSource
4708 pp = Float.T(
4709 default=1.0,
4710 help='initial excess pore pressure in [Pa]')
4712 length = Float.T(
4713 default=0.0,
4714 help='length of the line source [m]')
4716 azimuth = Float.T(
4717 default=0.0,
4718 help='azimuth direction, clockwise from north [deg]')
4720 dip = Float.T(
4721 default=90.,
4722 help='dip direction, downward from horizontal [deg]')
4724 def base_key(self):
4725 return Source.base_key(self) + (self.azimuth, self.dip, self.length)
4727 def get_factor(self):
4728 return self.pp
4730 def discretize_basesource(self, store, target=None):
4732 n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1
4734 a = num.linspace(-0.5 * self.length, 0.5 * self.length, n)
4736 sa = math.sin(self.azimuth * d2r)
4737 ca = math.cos(self.azimuth * d2r)
4738 sd = math.sin(self.dip * d2r)
4739 cd = math.cos(self.dip * d2r)
4741 points = num.zeros((n, 3))
4742 points[:, 0] = self.north_shift + a * ca * cd
4743 points[:, 1] = self.east_shift + a * sa * cd
4744 points[:, 2] = self.depth + a * sd
4746 return meta.DiscretizedPorePressureSource(
4747 times=util.num_full(n, self.time),
4748 lat=self.lat,
4749 lon=self.lon,
4750 north_shifts=points[:, 0],
4751 east_shifts=points[:, 1],
4752 depths=points[:, 2],
4753 pp=num.ones(n) / n)
4756class Request(Object):
4757 '''
4758 Synthetic seismogram computation request.
4760 ::
4762 Request(**kwargs)
4763 Request(sources, targets, **kwargs)
4764 '''
4766 sources = List.T(
4767 Source.T(),
4768 help='list of sources for which to produce synthetics.')
4770 targets = List.T(
4771 Target.T(),
4772 help='list of targets for which to produce synthetics.')
4774 @classmethod
4775 def args2kwargs(cls, args):
4776 if len(args) not in (0, 2, 3):
4777 raise BadRequest('Invalid arguments.')
4779 if len(args) == 2:
4780 return dict(sources=args[0], targets=args[1])
4781 else:
4782 return {}
4784 def __init__(self, *args, **kwargs):
4785 kwargs.update(self.args2kwargs(args))
4786 sources = kwargs.pop('sources', [])
4787 targets = kwargs.pop('targets', [])
4789 if isinstance(sources, Source):
4790 sources = [sources]
4792 if isinstance(targets, Target) or isinstance(targets, StaticTarget):
4793 targets = [targets]
4795 Object.__init__(self, sources=sources, targets=targets, **kwargs)
4797 @property
4798 def targets_dynamic(self):
4799 return [t for t in self.targets if isinstance(t, Target)]
4801 @property
4802 def targets_static(self):
4803 return [t for t in self.targets if isinstance(t, StaticTarget)]
4805 @property
4806 def has_dynamic(self):
4807 return True if len(self.targets_dynamic) > 0 else False
4809 @property
4810 def has_statics(self):
4811 return True if len(self.targets_static) > 0 else False
4813 def subsources_map(self):
4814 m = defaultdict(list)
4815 for source in self.sources:
4816 m[source.base_key()].append(source)
4818 return m
4820 def subtargets_map(self):
4821 m = defaultdict(list)
4822 for target in self.targets:
4823 m[target.base_key()].append(target)
4825 return m
4827 def subrequest_map(self):
4828 ms = self.subsources_map()
4829 mt = self.subtargets_map()
4830 m = {}
4831 for (ks, ls) in ms.items():
4832 for (kt, lt) in mt.items():
4833 m[ks, kt] = (ls, lt)
4835 return m
4838class ProcessingStats(Object):
4839 t_perc_get_store_and_receiver = Float.T(default=0.)
4840 t_perc_discretize_source = Float.T(default=0.)
4841 t_perc_make_base_seismogram = Float.T(default=0.)
4842 t_perc_make_same_span = Float.T(default=0.)
4843 t_perc_post_process = Float.T(default=0.)
4844 t_perc_optimize = Float.T(default=0.)
4845 t_perc_stack = Float.T(default=0.)
4846 t_perc_static_get_store = Float.T(default=0.)
4847 t_perc_static_discretize_basesource = Float.T(default=0.)
4848 t_perc_static_sum_statics = Float.T(default=0.)
4849 t_perc_static_post_process = Float.T(default=0.)
4850 t_wallclock = Float.T(default=0.)
4851 t_cpu = Float.T(default=0.)
4852 n_read_blocks = Int.T(default=0)
4853 n_results = Int.T(default=0)
4854 n_subrequests = Int.T(default=0)
4855 n_stores = Int.T(default=0)
4856 n_records_stacked = Int.T(default=0)
4859class Response(Object):
4860 '''
4861 Resonse object to a synthetic seismogram computation request.
4862 '''
4864 request = Request.T()
4865 results_list = List.T(List.T(meta.SeismosizerResult.T()))
4866 stats = ProcessingStats.T()
4868 def pyrocko_traces(self):
4869 '''
4870 Return a list of requested
4871 :class:`~pyrocko.trace.Trace` instances.
4872 '''
4874 traces = []
4875 for results in self.results_list:
4876 for result in results:
4877 if not isinstance(result, meta.Result):
4878 continue
4879 traces.append(result.trace.pyrocko_trace())
4881 return traces
4883 def kite_scenes(self):
4884 '''
4885 Return a list of requested
4886 :class:`kite.Scene` instances.
4887 '''
4888 kite_scenes = []
4889 for results in self.results_list:
4890 for result in results:
4891 if isinstance(result, meta.KiteSceneResult):
4892 sc = result.get_scene()
4893 kite_scenes.append(sc)
4895 return kite_scenes
4897 def static_results(self):
4898 '''
4899 Return a list of requested
4900 :class:`~pyrocko.gf.meta.StaticResult` instances.
4901 '''
4902 statics = []
4903 for results in self.results_list:
4904 for result in results:
4905 if not isinstance(result, meta.StaticResult):
4906 continue
4907 statics.append(result)
4909 return statics
4911 def iter_results(self, get='pyrocko_traces'):
4912 '''
4913 Generator function to iterate over results of request.
4915 Yields associated :py:class:`Source`,
4916 :class:`~pyrocko.gf.targets.Target`,
4917 :class:`~pyrocko.trace.Trace` instances in each iteration.
4918 '''
4920 for isource, source in enumerate(self.request.sources):
4921 for itarget, target in enumerate(self.request.targets):
4922 result = self.results_list[isource][itarget]
4923 if get == 'pyrocko_traces':
4924 yield source, target, result.trace.pyrocko_trace()
4925 elif get == 'results':
4926 yield source, target, result
4928 def snuffle(self, **kwargs):
4929 '''
4930 Open *snuffler* with requested traces.
4931 '''
4933 trace.snuffle(self.pyrocko_traces(), **kwargs)
4936class Engine(Object):
4937 '''
4938 Base class for synthetic seismogram calculators.
4939 '''
4941 def get_store_ids(self):
4942 '''
4943 Get list of available GF store IDs
4944 '''
4946 return []
4949class Rule(object):
4950 pass
4953class VectorRule(Rule):
4955 def __init__(self, quantity, differentiate=0, integrate=0):
4956 self.components = [quantity + '.' + c for c in 'ned']
4957 self.differentiate = differentiate
4958 self.integrate = integrate
4960 def required_components(self, target):
4961 n, e, d = self.components
4962 sa, ca, sd, cd = target.get_sin_cos_factors()
4964 comps = []
4965 if nonzero(ca * cd):
4966 comps.append(n)
4968 if nonzero(sa * cd):
4969 comps.append(e)
4971 if nonzero(sd):
4972 comps.append(d)
4974 return tuple(comps)
4976 def apply_(self, target, base_seismogram):
4977 n, e, d = self.components
4978 sa, ca, sd, cd = target.get_sin_cos_factors()
4980 if nonzero(ca * cd):
4981 data = base_seismogram[n].data * (ca * cd)
4982 deltat = base_seismogram[n].deltat
4983 else:
4984 data = 0.0
4986 if nonzero(sa * cd):
4987 data = data + base_seismogram[e].data * (sa * cd)
4988 deltat = base_seismogram[e].deltat
4990 if nonzero(sd):
4991 data = data + base_seismogram[d].data * sd
4992 deltat = base_seismogram[d].deltat
4994 if self.differentiate:
4995 data = util.diff_fd(self.differentiate, 4, deltat, data)
4997 if self.integrate:
4998 raise NotImplementedError('Integration is not implemented yet.')
5000 return data
5003class HorizontalVectorRule(Rule):
5005 def __init__(self, quantity, differentiate=0, integrate=0):
5006 self.components = [quantity + '.' + c for c in 'ne']
5007 self.differentiate = differentiate
5008 self.integrate = integrate
5010 def required_components(self, target):
5011 n, e = self.components
5012 sa, ca, _, _ = target.get_sin_cos_factors()
5014 comps = []
5015 if nonzero(ca):
5016 comps.append(n)
5018 if nonzero(sa):
5019 comps.append(e)
5021 return tuple(comps)
5023 def apply_(self, target, base_seismogram):
5024 n, e = self.components
5025 sa, ca, _, _ = target.get_sin_cos_factors()
5027 if nonzero(ca):
5028 data = base_seismogram[n].data * ca
5029 else:
5030 data = 0.0
5032 if nonzero(sa):
5033 data = data + base_seismogram[e].data * sa
5035 if self.differentiate:
5036 deltat = base_seismogram[e].deltat
5037 data = util.diff_fd(self.differentiate, 4, deltat, data)
5039 if self.integrate:
5040 raise NotImplementedError('Integration is not implemented yet.')
5042 return data
5045class ScalarRule(Rule):
5047 def __init__(self, quantity, differentiate=0):
5048 self.c = quantity
5050 def required_components(self, target):
5051 return (self.c, )
5053 def apply_(self, target, base_seismogram):
5054 data = base_seismogram[self.c].data.copy()
5055 deltat = base_seismogram[self.c].deltat
5056 if self.differentiate:
5057 data = util.diff_fd(self.differentiate, 4, deltat, data)
5059 return data
5062class StaticDisplacement(Rule):
5064 def required_components(self, target):
5065 return tuple(['displacement.%s' % c for c in list('ned')])
5067 def apply_(self, target, base_statics):
5068 if isinstance(target, SatelliteTarget):
5069 los_fac = target.get_los_factors()
5070 base_statics['displacement.los'] =\
5071 (los_fac[:, 0] * -base_statics['displacement.d'] +
5072 los_fac[:, 1] * base_statics['displacement.e'] +
5073 los_fac[:, 2] * base_statics['displacement.n'])
5074 return base_statics
5077channel_rules = {
5078 'displacement': [VectorRule('displacement')],
5079 'rotation': [VectorRule('rotation')],
5080 'velocity': [
5081 VectorRule('velocity'),
5082 VectorRule('displacement', differentiate=1)],
5083 'acceleration': [
5084 VectorRule('acceleration'),
5085 VectorRule('velocity', differentiate=1),
5086 VectorRule('displacement', differentiate=2)],
5087 'pore_pressure': [ScalarRule('pore_pressure')],
5088 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')],
5089 'darcy_velocity': [VectorRule('darcy_velocity')],
5090}
5092static_rules = {
5093 'displacement': [StaticDisplacement()]
5094}
5097class OutOfBoundsContext(Object):
5098 source = Source.T()
5099 target = Target.T()
5100 distance = Float.T()
5101 components = List.T(String.T())
5104def process_dynamic_timeseries(work, psources, ptargets, engine, nthreads=0):
5105 dsource_cache = {}
5106 tcounters = list(range(6))
5108 store_ids = set()
5109 sources = set()
5110 targets = set()
5112 for itarget, target in enumerate(ptargets):
5113 target._id = itarget
5115 for w in work:
5116 _, _, isources, itargets = w
5118 sources.update([psources[isource] for isource in isources])
5119 targets.update([ptargets[itarget] for itarget in itargets])
5121 store_ids = set([t.store_id for t in targets])
5123 for isource, source in enumerate(psources):
5125 components = set()
5126 for itarget, target in enumerate(targets):
5127 rule = engine.get_rule(source, target)
5128 components.update(rule.required_components(target))
5130 for store_id in store_ids:
5131 store_targets = [t for t in targets if t.store_id == store_id]
5133 sample_rates = set([t.sample_rate for t in store_targets])
5134 interpolations = set([t.interpolation for t in store_targets])
5136 base_seismograms = []
5137 store_targets_out = []
5139 for samp_rate in sample_rates:
5140 for interp in interpolations:
5141 engine_targets = [
5142 t for t in store_targets if t.sample_rate == samp_rate
5143 and t.interpolation == interp]
5145 if not engine_targets:
5146 continue
5148 store_targets_out += engine_targets
5150 base_seismograms += engine.base_seismograms(
5151 source,
5152 engine_targets,
5153 components,
5154 dsource_cache,
5155 nthreads)
5157 for iseis, seismogram in enumerate(base_seismograms):
5158 for tr in seismogram.values():
5159 if tr.err != store.SeismosizerErrorEnum.SUCCESS:
5160 e = SeismosizerError(
5161 'Seismosizer failed with return code %i\n%s' % (
5162 tr.err, str(
5163 OutOfBoundsContext(
5164 source=source,
5165 target=store_targets[iseis],
5166 distance=source.distance_to(
5167 store_targets[iseis]),
5168 components=components))))
5169 raise e
5171 for seismogram, target in zip(base_seismograms, store_targets_out):
5173 try:
5174 result = engine._post_process_dynamic(
5175 seismogram, source, target)
5176 except SeismosizerError as e:
5177 result = e
5179 yield (isource, target._id, result), tcounters
5182def process_dynamic(work, psources, ptargets, engine, nthreads=0):
5183 dsource_cache = {}
5185 for w in work:
5186 _, _, isources, itargets = w
5188 sources = [psources[isource] for isource in isources]
5189 targets = [ptargets[itarget] for itarget in itargets]
5191 components = set()
5192 for target in targets:
5193 rule = engine.get_rule(sources[0], target)
5194 components.update(rule.required_components(target))
5196 for isource, source in zip(isources, sources):
5197 for itarget, target in zip(itargets, targets):
5199 try:
5200 base_seismogram, tcounters = engine.base_seismogram(
5201 source, target, components, dsource_cache, nthreads)
5202 except meta.OutOfBounds as e:
5203 e.context = OutOfBoundsContext(
5204 source=sources[0],
5205 target=targets[0],
5206 distance=sources[0].distance_to(targets[0]),
5207 components=components)
5208 raise
5210 n_records_stacked = 0
5211 t_optimize = 0.0
5212 t_stack = 0.0
5214 for _, tr in base_seismogram.items():
5215 n_records_stacked += tr.n_records_stacked
5216 t_optimize += tr.t_optimize
5217 t_stack += tr.t_stack
5219 try:
5220 result = engine._post_process_dynamic(
5221 base_seismogram, source, target)
5222 result.n_records_stacked = n_records_stacked
5223 result.n_shared_stacking = len(sources) *\
5224 len(targets)
5225 result.t_optimize = t_optimize
5226 result.t_stack = t_stack
5227 except SeismosizerError as e:
5228 result = e
5230 tcounters.append(xtime())
5231 yield (isource, itarget, result), tcounters
5234def process_static(work, psources, ptargets, engine, nthreads=0):
5235 for w in work:
5236 _, _, isources, itargets = w
5238 sources = [psources[isource] for isource in isources]
5239 targets = [ptargets[itarget] for itarget in itargets]
5241 for isource, source in zip(isources, sources):
5242 for itarget, target in zip(itargets, targets):
5243 components = engine.get_rule(source, target)\
5244 .required_components(target)
5246 try:
5247 base_statics, tcounters = engine.base_statics(
5248 source, target, components, nthreads)
5249 except meta.OutOfBounds as e:
5250 e.context = OutOfBoundsContext(
5251 source=sources[0],
5252 target=targets[0],
5253 distance=float('nan'),
5254 components=components)
5255 raise
5256 result = engine._post_process_statics(
5257 base_statics, source, target)
5258 tcounters.append(xtime())
5260 yield (isource, itarget, result), tcounters
5263class LocalEngine(Engine):
5264 '''
5265 Offline synthetic seismogram calculator.
5267 :param use_env: if ``True``, fill :py:attr:`store_superdirs` and
5268 :py:attr:`store_dirs` with paths set in environment variables
5269 GF_STORE_SUPERDIRS and GF_STORE_DIRS.
5270 :param use_config: if ``True``, fill :py:attr:`store_superdirs` and
5271 :py:attr:`store_dirs` with paths set in the user's config file.
5273 The config file can be found at :file:`~/.pyrocko/config.pf`
5275 .. code-block :: python
5277 gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/']
5278 gf_store_superdirs: ['/home/pyrocko/gf_stores/']
5279 '''
5281 store_superdirs = List.T(
5282 String.T(),
5283 help="directories which are searched for Green's function stores")
5285 store_dirs = List.T(
5286 String.T(),
5287 help="additional individual Green's function store directories")
5289 default_store_id = String.T(
5290 optional=True,
5291 help='default store ID to be used when a request does not provide '
5292 'one')
5294 def __init__(self, **kwargs):
5295 use_env = kwargs.pop('use_env', False)
5296 use_config = kwargs.pop('use_config', False)
5297 Engine.__init__(self, **kwargs)
5298 if use_env:
5299 env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '')
5300 env_store_dirs = os.environ.get('GF_STORE_DIRS', '')
5301 if env_store_superdirs:
5302 self.store_superdirs.extend(env_store_superdirs.split(':'))
5304 if env_store_dirs:
5305 self.store_dirs.extend(env_store_dirs.split(':'))
5307 if use_config:
5308 c = config.config()
5309 self.store_superdirs.extend(c.gf_store_superdirs)
5310 self.store_dirs.extend(c.gf_store_dirs)
5312 self._check_store_dirs_type()
5313 self._id_to_store_dir = {}
5314 self._open_stores = {}
5315 self._effective_default_store_id = None
5317 def _check_store_dirs_type(self):
5318 for sdir in ['store_dirs', 'store_superdirs']:
5319 if not isinstance(self.__getattribute__(sdir), list):
5320 raise TypeError('{} of {} is not of type list'.format(
5321 sdir, self.__class__.__name__))
5323 def _get_store_id(self, store_dir):
5324 store_ = store.Store(store_dir)
5325 store_id = store_.config.id
5326 store_.close()
5327 return store_id
5329 def _looks_like_store_dir(self, store_dir):
5330 return os.path.isdir(store_dir) and \
5331 all(os.path.isfile(pjoin(store_dir, x)) for x in
5332 ('index', 'traces', 'config'))
5334 def iter_store_dirs(self):
5335 store_dirs = set()
5336 for d in self.store_superdirs:
5337 if not os.path.exists(d):
5338 logger.warning('store_superdir not available: %s' % d)
5339 continue
5341 for entry in os.listdir(d):
5342 store_dir = os.path.realpath(pjoin(d, entry))
5343 if self._looks_like_store_dir(store_dir):
5344 store_dirs.add(store_dir)
5346 for store_dir in self.store_dirs:
5347 store_dirs.add(os.path.realpath(store_dir))
5349 return store_dirs
5351 def _scan_stores(self):
5352 for store_dir in self.iter_store_dirs():
5353 store_id = self._get_store_id(store_dir)
5354 if store_id not in self._id_to_store_dir:
5355 self._id_to_store_dir[store_id] = store_dir
5356 else:
5357 if store_dir != self._id_to_store_dir[store_id]:
5358 raise DuplicateStoreId(
5359 'GF store ID %s is used in (at least) two '
5360 'different stores. Locations are: %s and %s' %
5361 (store_id, self._id_to_store_dir[store_id], store_dir))
5363 def get_store_dir(self, store_id):
5364 '''
5365 Lookup directory given a GF store ID.
5366 '''
5368 if store_id not in self._id_to_store_dir:
5369 self._scan_stores()
5371 if store_id not in self._id_to_store_dir:
5372 raise NoSuchStore(store_id, self.iter_store_dirs())
5374 return self._id_to_store_dir[store_id]
5376 def get_store_ids(self):
5377 '''
5378 Get list of available store IDs.
5379 '''
5381 self._scan_stores()
5382 return sorted(self._id_to_store_dir.keys())
5384 def effective_default_store_id(self):
5385 if self._effective_default_store_id is None:
5386 if self.default_store_id is None:
5387 store_ids = self.get_store_ids()
5388 if len(store_ids) == 1:
5389 self._effective_default_store_id = self.get_store_ids()[0]
5390 else:
5391 raise NoDefaultStoreSet()
5392 else:
5393 self._effective_default_store_id = self.default_store_id
5395 return self._effective_default_store_id
5397 def get_store(self, store_id=None):
5398 '''
5399 Get a store from the engine.
5401 :param store_id: identifier of the store (optional)
5402 :returns: :py:class:`~pyrocko.gf.store.Store` object
5404 If no ``store_id`` is provided the store
5405 associated with the :py:gattr:`default_store_id` is returned.
5406 Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is
5407 undefined.
5408 '''
5410 if store_id is None:
5411 store_id = self.effective_default_store_id()
5413 if store_id not in self._open_stores:
5414 store_dir = self.get_store_dir(store_id)
5415 self._open_stores[store_id] = store.Store(store_dir)
5417 return self._open_stores[store_id]
5419 def get_store_config(self, store_id):
5420 store = self.get_store(store_id)
5421 return store.config
5423 def get_store_extra(self, store_id, key):
5424 store = self.get_store(store_id)
5425 return store.get_extra(key)
5427 def close_cashed_stores(self):
5428 '''
5429 Close and remove ids from cashed stores.
5430 '''
5431 store_ids = []
5432 for store_id, store_ in self._open_stores.items():
5433 store_.close()
5434 store_ids.append(store_id)
5436 for store_id in store_ids:
5437 self._open_stores.pop(store_id)
5439 def get_rule(self, source, target):
5440 cprovided = self.get_store(target.store_id).get_provided_components()
5442 if isinstance(target, StaticTarget):
5443 quantity = target.quantity
5444 available_rules = static_rules
5445 elif isinstance(target, Target):
5446 quantity = target.effective_quantity()
5447 available_rules = channel_rules
5449 try:
5450 for rule in available_rules[quantity]:
5451 cneeded = rule.required_components(target)
5452 if all(c in cprovided for c in cneeded):
5453 return rule
5455 except KeyError:
5456 pass
5458 raise BadRequest(
5459 'No rule to calculate "%s" with GFs from store "%s" '
5460 'for source model "%s".' % (
5461 target.effective_quantity(),
5462 target.store_id,
5463 source.__class__.__name__))
5465 def _cached_discretize_basesource(self, source, store, cache, target):
5466 if (source, store) not in cache:
5467 cache[source, store] = source.discretize_basesource(store, target)
5469 return cache[source, store]
5471 def base_seismograms(self, source, targets, components, dsource_cache,
5472 nthreads=0):
5474 target = targets[0]
5476 interp = set([t.interpolation for t in targets])
5477 if len(interp) > 1:
5478 raise BadRequest('Targets have different interpolation schemes.')
5480 rates = set([t.sample_rate for t in targets])
5481 if len(rates) > 1:
5482 raise BadRequest('Targets have different sample rates.')
5484 store_ = self.get_store(target.store_id)
5485 receivers = [t.receiver(store_) for t in targets]
5487 if target.sample_rate is not None:
5488 deltat = 1. / target.sample_rate
5489 rate = target.sample_rate
5490 else:
5491 deltat = None
5492 rate = store_.config.sample_rate
5494 tmin = num.fromiter(
5495 (t.tmin for t in targets), dtype=float, count=len(targets))
5496 tmax = num.fromiter(
5497 (t.tmax for t in targets), dtype=float, count=len(targets))
5499 mask = num.logical_and(num.isfinite(tmin), num.isfinite(tmax))
5501 itmin = num.zeros_like(tmin, dtype=num.int64)
5502 itmax = num.zeros_like(tmin, dtype=num.int64)
5503 nsamples = num.full_like(tmin, -1, dtype=num.int64)
5505 itmin[mask] = num.floor(tmin[mask] * rate).astype(num.int64)
5506 itmax[mask] = num.ceil(tmax[mask] * rate).astype(num.int64)
5507 nsamples = itmax - itmin + 1
5508 nsamples[num.logical_not(mask)] = -1
5510 base_source = self._cached_discretize_basesource(
5511 source, store_, dsource_cache, target)
5513 base_seismograms = store_.calc_seismograms(
5514 base_source, receivers, components,
5515 deltat=deltat,
5516 itmin=itmin, nsamples=nsamples,
5517 interpolation=target.interpolation,
5518 optimization=target.optimization,
5519 nthreads=nthreads)
5521 for i, base_seismogram in enumerate(base_seismograms):
5522 base_seismograms[i] = store.make_same_span(base_seismogram)
5524 return base_seismograms
5526 def base_seismogram(self, source, target, components, dsource_cache,
5527 nthreads):
5529 tcounters = [xtime()]
5531 store_ = self.get_store(target.store_id)
5532 receiver = target.receiver(store_)
5534 if target.tmin and target.tmax is not None:
5535 rate = store_.config.sample_rate
5536 itmin = int(num.floor(target.tmin * rate))
5537 itmax = int(num.ceil(target.tmax * rate))
5538 nsamples = itmax - itmin + 1
5539 else:
5540 itmin = None
5541 nsamples = None
5543 tcounters.append(xtime())
5544 base_source = self._cached_discretize_basesource(
5545 source, store_, dsource_cache, target)
5547 tcounters.append(xtime())
5549 if target.sample_rate is not None:
5550 deltat = 1. / target.sample_rate
5551 else:
5552 deltat = None
5554 base_seismogram = store_.seismogram(
5555 base_source, receiver, components,
5556 deltat=deltat,
5557 itmin=itmin, nsamples=nsamples,
5558 interpolation=target.interpolation,
5559 optimization=target.optimization,
5560 nthreads=nthreads)
5562 tcounters.append(xtime())
5564 base_seismogram = store.make_same_span(base_seismogram)
5566 tcounters.append(xtime())
5568 return base_seismogram, tcounters
5570 def base_statics(self, source, target, components, nthreads):
5571 tcounters = [xtime()]
5572 store_ = self.get_store(target.store_id)
5574 if target.tsnapshot is not None:
5575 rate = store_.config.sample_rate
5576 itsnapshot = int(num.floor(target.tsnapshot * rate))
5577 else:
5578 itsnapshot = None
5579 tcounters.append(xtime())
5581 base_source = source.discretize_basesource(store_, target=target)
5583 tcounters.append(xtime())
5585 base_statics = store_.statics(
5586 base_source,
5587 target,
5588 itsnapshot,
5589 components,
5590 target.interpolation,
5591 nthreads)
5593 tcounters.append(xtime())
5595 return base_statics, tcounters
5597 def _post_process_dynamic(self, base_seismogram, source, target):
5598 base_any = next(iter(base_seismogram.values()))
5599 deltat = base_any.deltat
5600 itmin = base_any.itmin
5602 rule = self.get_rule(source, target)
5603 data = rule.apply_(target, base_seismogram)
5605 factor = source.get_factor() * target.get_factor()
5606 if factor != 1.0:
5607 data = data * factor
5609 stf = source.effective_stf_post()
5611 times, amplitudes = stf.discretize_t(
5612 deltat, 0.0)
5614 # repeat end point to prevent boundary effects
5615 padded_data = num.empty(data.size + amplitudes.size, dtype=float)
5616 padded_data[:data.size] = data
5617 padded_data[data.size:] = data[-1]
5618 data = num.convolve(amplitudes, padded_data)
5620 tmin = itmin * deltat + times[0]
5622 tr = meta.SeismosizerTrace(
5623 codes=target.codes,
5624 data=data[:-amplitudes.size],
5625 deltat=deltat,
5626 tmin=tmin)
5628 return target.post_process(self, source, tr)
5630 def _post_process_statics(self, base_statics, source, starget):
5631 rule = self.get_rule(source, starget)
5632 data = rule.apply_(starget, base_statics)
5634 factor = source.get_factor()
5635 if factor != 1.0:
5636 for v in data.values():
5637 v *= factor
5639 return starget.post_process(self, source, base_statics)
5641 def process(self, *args, **kwargs):
5642 '''
5643 Process a request.
5645 ::
5647 process(**kwargs)
5648 process(request, **kwargs)
5649 process(sources, targets, **kwargs)
5651 The request can be given a a :py:class:`Request` object, or such an
5652 object is created using ``Request(**kwargs)`` for convenience.
5654 :returns: :py:class:`Response` object
5655 '''
5657 if len(args) not in (0, 1, 2):
5658 raise BadRequest('Invalid arguments.')
5660 if len(args) == 1:
5661 kwargs['request'] = args[0]
5663 elif len(args) == 2:
5664 kwargs.update(Request.args2kwargs(args))
5666 request = kwargs.pop('request', None)
5667 status_callback = kwargs.pop('status_callback', None)
5668 calc_timeseries = kwargs.pop('calc_timeseries', True)
5670 nprocs = kwargs.pop('nprocs', None)
5671 nthreads = kwargs.pop('nthreads', 1)
5672 if nprocs is not None:
5673 nthreads = nprocs
5675 if request is None:
5676 request = Request(**kwargs)
5678 if resource:
5679 rs0 = resource.getrusage(resource.RUSAGE_SELF)
5680 rc0 = resource.getrusage(resource.RUSAGE_CHILDREN)
5681 tt0 = xtime()
5683 # make sure stores are open before fork()
5684 store_ids = set(target.store_id for target in request.targets)
5685 for store_id in store_ids:
5686 self.get_store(store_id)
5688 source_index = dict((x, i) for (i, x) in
5689 enumerate(request.sources))
5690 target_index = dict((x, i) for (i, x) in
5691 enumerate(request.targets))
5693 m = request.subrequest_map()
5695 skeys = sorted(m.keys(), key=cmp_to_key(cmp_none_aware))
5696 results_list = []
5698 for i in range(len(request.sources)):
5699 results_list.append([None] * len(request.targets))
5701 tcounters_dyn_list = []
5702 tcounters_static_list = []
5703 nsub = len(skeys)
5704 isub = 0
5706 # Processing dynamic targets through
5707 # parimap(process_subrequest_dynamic)
5709 if calc_timeseries:
5710 _process_dynamic = process_dynamic_timeseries
5711 else:
5712 _process_dynamic = process_dynamic
5714 if request.has_dynamic:
5715 work_dynamic = [
5716 (i, nsub,
5717 [source_index[source] for source in m[k][0]],
5718 [target_index[target] for target in m[k][1]
5719 if not isinstance(target, StaticTarget)])
5720 for (i, k) in enumerate(skeys)]
5722 for ii_results, tcounters_dyn in _process_dynamic(
5723 work_dynamic, request.sources, request.targets, self,
5724 nthreads):
5726 tcounters_dyn_list.append(num.diff(tcounters_dyn))
5727 isource, itarget, result = ii_results
5728 results_list[isource][itarget] = result
5730 if status_callback:
5731 status_callback(isub, nsub)
5733 isub += 1
5735 # Processing static targets through process_static
5736 if request.has_statics:
5737 work_static = [
5738 (i, nsub,
5739 [source_index[source] for source in m[k][0]],
5740 [target_index[target] for target in m[k][1]
5741 if isinstance(target, StaticTarget)])
5742 for (i, k) in enumerate(skeys)]
5744 for ii_results, tcounters_static in process_static(
5745 work_static, request.sources, request.targets, self,
5746 nthreads=nthreads):
5748 tcounters_static_list.append(num.diff(tcounters_static))
5749 isource, itarget, result = ii_results
5750 results_list[isource][itarget] = result
5752 if status_callback:
5753 status_callback(isub, nsub)
5755 isub += 1
5757 if status_callback:
5758 status_callback(nsub, nsub)
5760 tt1 = time.time()
5761 if resource:
5762 rs1 = resource.getrusage(resource.RUSAGE_SELF)
5763 rc1 = resource.getrusage(resource.RUSAGE_CHILDREN)
5765 s = ProcessingStats()
5767 if request.has_dynamic:
5768 tcumu_dyn = num.sum(num.vstack(tcounters_dyn_list), axis=0)
5769 t_dyn = float(num.sum(tcumu_dyn))
5770 perc_dyn = map(float, tcumu_dyn / t_dyn * 100.)
5771 (s.t_perc_get_store_and_receiver,
5772 s.t_perc_discretize_source,
5773 s.t_perc_make_base_seismogram,
5774 s.t_perc_make_same_span,
5775 s.t_perc_post_process) = perc_dyn
5776 else:
5777 t_dyn = 0.
5779 if request.has_statics:
5780 tcumu_static = num.sum(num.vstack(tcounters_static_list), axis=0)
5781 t_static = num.sum(tcumu_static)
5782 perc_static = map(float, tcumu_static / t_static * 100.)
5783 (s.t_perc_static_get_store,
5784 s.t_perc_static_discretize_basesource,
5785 s.t_perc_static_sum_statics,
5786 s.t_perc_static_post_process) = perc_static
5788 s.t_wallclock = tt1 - tt0
5789 if resource:
5790 s.t_cpu = (
5791 (rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) -
5792 (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime))
5793 s.n_read_blocks = (
5794 (rs1.ru_inblock + rc1.ru_inblock) -
5795 (rs0.ru_inblock + rc0.ru_inblock))
5797 n_records_stacked = 0.
5798 for results in results_list:
5799 for result in results:
5800 if not isinstance(result, meta.Result):
5801 continue
5802 shr = float(result.n_shared_stacking)
5803 n_records_stacked += result.n_records_stacked / shr
5804 s.t_perc_optimize += result.t_optimize / shr
5805 s.t_perc_stack += result.t_stack / shr
5806 s.n_records_stacked = int(n_records_stacked)
5807 if t_dyn != 0.:
5808 s.t_perc_optimize /= t_dyn * 100
5809 s.t_perc_stack /= t_dyn * 100
5811 return Response(
5812 request=request,
5813 results_list=results_list,
5814 stats=s)
5817class RemoteEngine(Engine):
5818 '''
5819 Client for remote synthetic seismogram calculator.
5820 '''
5822 site = String.T(default=ws.g_default_site, optional=True)
5823 url = String.T(default=ws.g_url, optional=True)
5825 def process(self, request=None, status_callback=None, **kwargs):
5827 if request is None:
5828 request = Request(**kwargs)
5830 return ws.seismosizer(url=self.url, site=self.site, request=request)
5833g_engine = None
5836def get_engine(store_superdirs=[]):
5837 global g_engine
5838 if g_engine is None:
5839 g_engine = LocalEngine(use_env=True, use_config=True)
5841 for d in store_superdirs:
5842 if d not in g_engine.store_superdirs:
5843 g_engine.store_superdirs.append(d)
5845 return g_engine
5848class SourceGroup(Object):
5850 def __getattr__(self, k):
5851 return num.fromiter((getattr(s, k) for s in self),
5852 dtype=float)
5854 def __iter__(self):
5855 raise NotImplementedError(
5856 'This method should be implemented in subclass.')
5858 def __len__(self):
5859 raise NotImplementedError(
5860 'This method should be implemented in subclass.')
5863class SourceList(SourceGroup):
5864 sources = List.T(Source.T())
5866 def append(self, s):
5867 self.sources.append(s)
5869 def __iter__(self):
5870 return iter(self.sources)
5872 def __len__(self):
5873 return len(self.sources)
5876class SourceGrid(SourceGroup):
5878 base = Source.T()
5879 variables = Dict.T(String.T(), Range.T())
5880 order = List.T(String.T())
5882 def __len__(self):
5883 n = 1
5884 for (k, v) in self.make_coords(self.base):
5885 n *= len(list(v))
5887 return n
5889 def __iter__(self):
5890 for items in permudef(self.make_coords(self.base)):
5891 s = self.base.clone(**{k: v for (k, v) in items})
5892 s.regularize()
5893 yield s
5895 def ordered_params(self):
5896 ks = list(self.variables.keys())
5897 for k in self.order + list(self.base.keys()):
5898 if k in ks:
5899 yield k
5900 ks.remove(k)
5901 if ks:
5902 raise Exception('Invalid parameter "%s" for source type "%s".' %
5903 (ks[0], self.base.__class__.__name__))
5905 def make_coords(self, base):
5906 return [(param, self.variables[param].make(base=base[param]))
5907 for param in self.ordered_params()]
5910source_classes = [
5911 Source,
5912 SourceWithMagnitude,
5913 SourceWithDerivedMagnitude,
5914 ExplosionSource,
5915 RectangularExplosionSource,
5916 DCSource,
5917 CLVDSource,
5918 VLVDSource,
5919 MTSource,
5920 RectangularSource,
5921 PseudoDynamicRupture,
5922 DoubleDCSource,
5923 RingfaultSource,
5924 CombiSource,
5925 SFSource,
5926 PorePressurePointSource,
5927 PorePressureLineSource,
5928]
5930stf_classes = [
5931 STF,
5932 BoxcarSTF,
5933 TriangularSTF,
5934 HalfSinusoidSTF,
5935 ResonatorSTF,
5936 TremorSTF,
5937]
5939__all__ = '''
5940Cloneable
5941NoDefaultStoreSet
5942SeismosizerError
5943BadRequest
5944NoSuchStore
5945DerivedMagnitudeError
5946STFMode
5947'''.split() + [S.__name__ for S in source_classes + stf_classes] + '''
5948Request
5949ProcessingStats
5950Response
5951Engine
5952LocalEngine
5953RemoteEngine
5954source_classes
5955get_engine
5956Range
5957SourceGroup
5958SourceList
5959SourceGrid
5960map_anchor
5961'''.split()