Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/seismosizer.py: 83%
2580 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-10-02 07:18 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7High level synthetic seismogram synthesis.
9.. _coordinate-system-names:
11Coordinate systems
12..................
14Coordinate system names commonly used in source models.
16================= ============================================
17Name Description
18================= ============================================
19``'xyz'`` northing, easting, depth in [m]
20``'xy'`` northing, easting in [m]
21``'latlon'`` latitude, longitude in [deg]
22``'lonlat'`` longitude, latitude in [deg]
23``'latlondepth'`` latitude, longitude in [deg], depth in [m]
24================= ============================================
25'''
27from collections import defaultdict
28from functools import cmp_to_key
29import time
30import math
31import os
32import re
33import logging
34try:
35 import resource
36except ImportError:
37 resource = None
38from hashlib import sha1
40import numpy as num
41from scipy.interpolate import RegularGridInterpolator
43from pyrocko.guts import (Object, Float, String, StringChoice, List,
44 Timestamp, Int, SObject, ArgumentError, Dict,
45 ValidationError, Bool)
46from pyrocko.guts_array import Array
48from pyrocko import moment_tensor as pmt
49from pyrocko import trace, util, config, model, eikonal_ext
50from pyrocko.orthodrome import ne_to_latlon
51from pyrocko.model import Location
52from pyrocko.modelling import OkadaSource, make_okada_coefficient_matrix, \
53 okada_ext, invert_fault_dislocations_bem
55from . import meta, store, ws
56from .tractions import TractionField, DirectedTractions
57from .targets import Target, StaticTarget, SatelliteTarget
59pjoin = os.path.join
61guts_prefix = 'pf'
63d2r = math.pi / 180.
64r2d = 180. / math.pi
65km = 1e3
67logger = logging.getLogger('pyrocko.gf.seismosizer')
70def cmp_none_aware(a, b):
71 if isinstance(a, tuple) and isinstance(b, tuple):
72 for xa, xb in zip(a, b):
73 rv = cmp_none_aware(xa, xb)
74 if rv != 0:
75 return rv
77 return 0
79 anone = a is None
80 bnone = b is None
82 if anone and bnone:
83 return 0
85 if anone:
86 return -1
88 if bnone:
89 return 1
91 return bool(a > b) - bool(a < b)
94def xtime():
95 return time.time()
98class SeismosizerError(Exception):
99 pass
102class BadRequest(SeismosizerError):
103 pass
106class DuplicateStoreId(Exception):
107 pass
110class NoDefaultStoreSet(Exception):
111 '''
112 Raised, when a default store would be used but none is set.
113 '''
114 pass
117class ConversionError(Exception):
118 pass
121class STFError(SeismosizerError):
122 pass
125class NoSuchStore(BadRequest):
127 def __init__(self, store_id=None, dirs=None):
128 BadRequest.__init__(self)
129 self.store_id = store_id
130 self.dirs = dirs
132 def __str__(self):
133 if self.store_id is not None:
134 rstr = 'no GF store with id "%s" found.' % self.store_id
135 else:
136 rstr = 'GF store not found.'
138 if self.dirs is not None:
139 rstr += ' Searched folders:\n %s' % '\n '.join(sorted(self.dirs))
140 return rstr
143def ufloat(s):
144 units = {
145 'k': 1e3,
146 'M': 1e6,
147 }
149 factor = 1.0
150 if s and s[-1] in units:
151 factor = units[s[-1]]
152 s = s[:-1]
153 if not s:
154 raise ValueError("unit without a number: '%s'" % s)
156 return float(s) * factor
159def ufloat_or_none(s):
160 if s:
161 return ufloat(s)
162 else:
163 return None
166def int_or_none(s):
167 if s:
168 return int(s)
169 else:
170 return None
173def nonzero(x, eps=1e-15):
174 return abs(x) > eps
177def permudef(ln, j=0):
178 if j < len(ln):
179 k, v = ln[j]
180 for y in v:
181 ln[j] = k, y
182 for s in permudef(ln, j + 1):
183 yield s
185 ln[j] = k, v
186 return
187 else:
188 yield ln
191def arr(x):
192 return num.atleast_1d(num.asarray(x))
195def discretize_rect_source(deltas, deltat, time, north, east, depth,
196 strike, dip, length, width,
197 anchor, velocity=None, stf=None,
198 nucleation_x=None, nucleation_y=None,
199 decimation_factor=1, pointsonly=False,
200 plane_coords=False,
201 aggressive_oversampling=False):
203 if stf is None:
204 stf = STF()
206 if not velocity and not pointsonly:
207 raise AttributeError('velocity is required in time mode')
209 mindeltagf = float(num.min(deltas))
210 if velocity:
211 mindeltagf = min(mindeltagf, deltat * velocity)
213 ln = length
214 wd = width
216 if aggressive_oversampling:
217 nl = int((2. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
218 nw = int((2. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
219 else:
220 nl = int((1. / decimation_factor) * num.ceil(ln / mindeltagf)) + 1
221 nw = int((1. / decimation_factor) * num.ceil(wd / mindeltagf)) + 1
223 n = int(nl * nw)
225 dl = ln / nl
226 dw = wd / nw
228 xl = num.linspace(-0.5 * (ln - dl), 0.5 * (ln - dl), nl)
229 xw = num.linspace(-0.5 * (wd - dw), 0.5 * (wd - dw), nw)
231 points = num.zeros((n, 3))
232 points[:, 0] = num.tile(xl, nw)
233 points[:, 1] = num.repeat(xw, nl)
235 if nucleation_x is not None:
236 dist_x = num.abs(nucleation_x - points[:, 0])
237 else:
238 dist_x = num.zeros(n)
240 if nucleation_y is not None:
241 dist_y = num.abs(nucleation_y - points[:, 1])
242 else:
243 dist_y = num.zeros(n)
245 dist = num.sqrt(dist_x**2 + dist_y**2)
246 times = dist / velocity
248 anch_x, anch_y = map_anchor[anchor]
250 points[:, 0] -= anch_x * 0.5 * length
251 points[:, 1] -= anch_y * 0.5 * width
253 if plane_coords:
254 return points, dl, dw, nl, nw
256 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
257 points = num.dot(rotmat.T, points.T).T
259 points[:, 0] += north
260 points[:, 1] += east
261 points[:, 2] += depth
263 if pointsonly:
264 return points, dl, dw, nl, nw
266 xtau, amplitudes = stf.discretize_t(deltat, time)
267 nt = xtau.size
269 points2 = num.repeat(points, nt, axis=0)
270 times2 = (times[:, num.newaxis] + xtau[num.newaxis, :]).ravel()
271 amplitudes2 = num.tile(amplitudes, n)
273 return points2, times2, amplitudes2, dl, dw, nl, nw
276def check_rect_source_discretisation(points2, nl, nw, store):
277 # We assume a non-rotated fault plane
278 N_CRITICAL = 8
279 points = points2.T.reshape((3, nl, nw))
280 if points.size <= N_CRITICAL:
281 logger.warning('RectangularSource is defined by only %d sub-sources!'
282 % points.size)
283 return True
285 distances = num.sqrt(
286 (points[0, 0, :] - points[0, 1, :])**2 +
287 (points[1, 0, :] - points[1, 1, :])**2 +
288 (points[2, 0, :] - points[2, 1, :])**2)
290 depths = points[2, 0, :]
291 vs_profile = store.config.get_vs(
292 lat=0., lon=0.,
293 points=num.repeat(depths[:, num.newaxis], 3, axis=1),
294 interpolation='multilinear')
296 min_wavelength = vs_profile * (store.config.deltat * 2)
297 if not num.all(min_wavelength > distances / 2):
298 return False
299 return True
302def outline_rect_source(strike, dip, length, width, anchor):
303 ln = length
304 wd = width
305 points = num.array(
306 [[-0.5 * ln, -0.5 * wd, 0.],
307 [0.5 * ln, -0.5 * wd, 0.],
308 [0.5 * ln, 0.5 * wd, 0.],
309 [-0.5 * ln, 0.5 * wd, 0.],
310 [-0.5 * ln, -0.5 * wd, 0.]])
312 anch_x, anch_y = map_anchor[anchor]
313 points[:, 0] -= anch_x * 0.5 * length
314 points[:, 1] -= anch_y * 0.5 * width
316 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
318 return num.dot(rotmat.T, points.T).T
321def from_plane_coords(
322 strike, dip, length, width, depth, x_plane_coords, y_plane_coords,
323 lat=0., lon=0.,
324 north_shift=0, east_shift=0,
325 anchor='top', cs='xy'):
327 ln = length
328 wd = width
329 x_abs = []
330 y_abs = []
331 if not isinstance(x_plane_coords, list):
332 x_plane_coords = [x_plane_coords]
333 y_plane_coords = [y_plane_coords]
335 for x_plane, y_plane in zip(x_plane_coords, y_plane_coords):
336 points = num.array(
337 [[-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
338 [0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.],
339 [0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
340 [-0.5 * ln * x_plane, 0.5 * wd * y_plane, 0.],
341 [-0.5 * ln * x_plane, -0.5 * wd * y_plane, 0.]])
343 anch_x, anch_y = map_anchor[anchor]
344 points[:, 0] -= anch_x * 0.5 * length
345 points[:, 1] -= anch_y * 0.5 * width
347 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
349 points = num.dot(rotmat.T, points.T).T
350 points[:, 0] += north_shift
351 points[:, 1] += east_shift
352 points[:, 2] += depth
353 if cs in ('latlon', 'lonlat'):
354 latlon = ne_to_latlon(lat, lon,
355 points[:, 0], points[:, 1])
356 latlon = num.array(latlon).T
357 x_abs.append(latlon[1:2, 1])
358 y_abs.append(latlon[2:3, 0])
359 if cs == 'xy':
360 x_abs.append(points[1:2, 1])
361 y_abs.append(points[2:3, 0])
363 if cs == 'lonlat':
364 return y_abs, x_abs
365 else:
366 return x_abs, y_abs
369def points_on_rect_source(
370 strike, dip, length, width, anchor,
371 discretized_basesource=None, points_x=None, points_y=None):
373 ln = length
374 wd = width
376 if isinstance(points_x, list) or isinstance(points_x, float):
377 points_x = num.array([points_x])
378 if isinstance(points_y, list) or isinstance(points_y, float):
379 points_y = num.array([points_y])
381 if discretized_basesource:
382 ds = discretized_basesource
384 nl_patches = ds.nl + 1
385 nw_patches = ds.nw + 1
387 npoints = nl_patches * nw_patches
388 points = num.zeros((npoints, 3))
389 ln_patches = num.array([il for il in range(nl_patches)])
390 wd_patches = num.array([iw for iw in range(nw_patches)])
392 points_ln =\
393 2 * ((ln_patches - num.min(ln_patches)) / num.ptp(ln_patches)) - 1
394 points_wd =\
395 2 * ((wd_patches - num.min(wd_patches)) / num.ptp(wd_patches)) - 1
397 for il in range(nl_patches):
398 for iw in range(nw_patches):
399 points[il * nw_patches + iw, :] = num.array([
400 points_ln[il] * ln * 0.5,
401 points_wd[iw] * wd * 0.5, 0.0])
403 elif points_x.shape[0] > 0 and points_y.shape[0] > 0:
404 points = num.zeros(shape=((len(points_x), 3)))
405 for i, (x, y) in enumerate(zip(points_x, points_y)):
406 points[i, :] = num.array(
407 [x * 0.5 * ln, y * 0.5 * wd, 0.0])
409 anch_x, anch_y = map_anchor[anchor]
411 points[:, 0] -= anch_x * 0.5 * ln
412 points[:, 1] -= anch_y * 0.5 * wd
414 rotmat = pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0)
416 return num.dot(rotmat.T, points.T).T
419class InvalidGridDef(Exception):
420 pass
423class Range(SObject):
424 '''
425 Convenient range specification.
427 Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]::
429 Range('0 .. 10k : 1k')
430 Range(start=0., stop=10e3, step=1e3)
431 Range(0, 10e3, 1e3)
432 Range('0 .. 10k @ 11')
433 Range(start=0., stop=10*km, n=11)
435 Range(0, 10e3, n=11)
436 Range(values=[x*1e3 for x in range(11)])
438 Depending on the use context, it can be possible to omit any part of the
439 specification. E.g. in the context of extracting a subset of an already
440 existing range, the existing range's specification values would be filled
441 in where missing.
443 The values are distributed with equal spacing, unless the ``spacing``
444 argument is modified. The values can be created offset or relative to an
445 external base value with the ``relative`` argument if the use context
446 supports this.
448 The range specification can be expressed with a short string
449 representation::
451 'start .. stop @ num | spacing, relative'
452 'start .. stop : step | spacing, relative'
454 most parts of the expression can be omitted if not needed. Whitespace is
455 allowed for readability but can also be omitted.
456 '''
458 start = Float.T(optional=True)
459 stop = Float.T(optional=True)
460 step = Float.T(optional=True)
461 n = Int.T(optional=True)
462 values = Array.T(optional=True, dtype=float, shape=(None,))
464 spacing = StringChoice.T(
465 choices=['lin', 'log', 'symlog'],
466 default='lin',
467 optional=True)
469 relative = StringChoice.T(
470 choices=['', 'add', 'mult'],
471 default='',
472 optional=True)
474 pattern = re.compile(r'^((?P<start>.*)\.\.(?P<stop>[^@|:]*))?'
475 r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?'
476 r'(\|(?P<stuff>.+))?$')
478 def __init__(self, *args, **kwargs):
479 d = {}
480 if len(args) == 1:
481 d = self.parse(args[0])
482 elif len(args) in (2, 3):
483 d['start'], d['stop'] = [float(x) for x in args[:2]]
484 if len(args) == 3:
485 d['step'] = float(args[2])
487 for k, v in kwargs.items():
488 if k in d:
489 raise ArgumentError('%s specified more than once' % k)
491 d[k] = v
493 SObject.__init__(self, **d)
495 def __str__(self):
496 def sfloat(x):
497 if x is not None:
498 return '%g' % x
499 else:
500 return ''
502 if self.values:
503 return ','.join('%g' % x for x in self.values)
505 if self.start is None and self.stop is None:
506 s0 = ''
507 else:
508 s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop))
510 s1 = ''
511 if self.step is not None:
512 s1 = [' : %g', ':%g'][s0 == ''] % self.step
513 elif self.n is not None:
514 s1 = [' @ %i', '@%i'][s0 == ''] % self.n
516 if self.spacing == 'lin' and self.relative == '':
517 s2 = ''
518 else:
519 x = []
520 if self.spacing != 'lin':
521 x.append(self.spacing)
523 if self.relative != '':
524 x.append(self.relative)
526 s2 = ' | %s' % ','.join(x)
528 return s0 + s1 + s2
530 @classmethod
531 def parse(cls, s):
532 s = re.sub(r'\s+', '', s)
533 m = cls.pattern.match(s)
534 if not m:
535 try:
536 vals = [ufloat(x) for x in s.split(',')]
537 except Exception:
538 raise InvalidGridDef(
539 '"%s" is not a valid range specification' % s)
541 return dict(values=num.array(vals, dtype=float))
543 d = m.groupdict()
544 try:
545 start = ufloat_or_none(d['start'])
546 stop = ufloat_or_none(d['stop'])
547 step = ufloat_or_none(d['step'])
548 n = int_or_none(d['n'])
549 except Exception:
550 raise InvalidGridDef(
551 '"%s" is not a valid range specification' % s)
553 spacing = 'lin'
554 relative = ''
556 if d['stuff'] is not None:
557 t = d['stuff'].split(',')
558 for x in t:
559 if x in cls.spacing.choices:
560 spacing = x
561 elif x and x in cls.relative.choices:
562 relative = x
563 else:
564 raise InvalidGridDef(
565 '"%s" is not a valid range specification' % s)
567 return dict(start=start, stop=stop, step=step, n=n, spacing=spacing,
568 relative=relative)
570 def make(self, mi=None, ma=None, inc=None, base=None, eps=1e-5):
571 if self.values:
572 return self.values
574 start = self.start
575 stop = self.stop
576 step = self.step
577 n = self.n
579 swap = step is not None and step < 0.
580 if start is None:
581 start = [mi, ma][swap]
582 if stop is None:
583 stop = [ma, mi][swap]
584 if step is None and inc is not None:
585 step = [inc, -inc][ma < mi]
587 if start is None or stop is None:
588 raise InvalidGridDef(
589 'Cannot use range specification "%s" without start '
590 'and stop in this context' % self)
592 if step is None and n is None:
593 step = stop - start
595 if n is None:
596 if (step < 0) != (stop - start < 0):
597 raise InvalidGridDef(
598 'Range specification "%s" has inconsistent ordering '
599 '(step < 0 => stop > start)' % self)
601 n = int(round((stop - start) / step)) + 1
602 stop2 = start + (n - 1) * step
603 if abs(stop - stop2) > eps:
604 n = int(math.floor((stop - start) / step)) + 1
605 stop = start + (n - 1) * step
606 else:
607 stop = stop2
609 if start == stop:
610 n = 1
612 if self.spacing == 'lin':
613 vals = num.linspace(start, stop, n)
615 elif self.spacing in ('log', 'symlog'):
616 if start > 0. and stop > 0.:
617 vals = num.exp(num.linspace(num.log(start),
618 num.log(stop), n))
619 elif start < 0. and stop < 0.:
620 vals = -num.exp(num.linspace(num.log(-start),
621 num.log(-stop), n))
622 else:
623 raise InvalidGridDef(
624 'Log ranges should not include or cross zero '
625 '(in range specification "%s").' % self)
627 if self.spacing == 'symlog':
628 nvals = - vals
629 vals = num.concatenate((nvals[::-1], vals))
631 if self.relative in ('add', 'mult') and base is None:
632 raise InvalidGridDef(
633 'Cannot use relative range specification in this context.')
635 vals = self.make_relative(base, vals)
637 return list(map(float, vals))
639 def make_relative(self, base, vals):
640 if self.relative == 'add':
641 vals += base
643 if self.relative == 'mult':
644 vals *= base
646 return vals
649class GridDefElement(Object):
651 param = meta.StringID.T()
652 rs = Range.T()
654 def __init__(self, shorthand=None, **kwargs):
655 if shorthand is not None:
656 t = shorthand.split('=')
657 if len(t) != 2:
658 raise InvalidGridDef(
659 'Invalid grid specification element: %s' % shorthand)
661 sp, sr = t[0].strip(), t[1].strip()
663 kwargs['param'] = sp
664 kwargs['rs'] = Range(sr)
666 Object.__init__(self, **kwargs)
668 def shorthand(self):
669 return self.param + ' = ' + str(self.rs)
672class GridDef(Object):
674 elements = List.T(GridDefElement.T())
676 def __init__(self, shorthand=None, **kwargs):
677 if shorthand is not None:
678 t = shorthand.splitlines()
679 tt = []
680 for x in t:
681 x = x.strip()
682 if x:
683 tt.extend(x.split(';'))
685 elements = []
686 for se in tt:
687 elements.append(GridDef(se))
689 kwargs['elements'] = elements
691 Object.__init__(self, **kwargs)
693 def shorthand(self):
694 return '; '.join(str(x) for x in self.elements)
697class Cloneable(object):
698 '''
699 Mix-in class for Guts objects, providing dict-like access and cloning.
700 '''
702 def __iter__(self):
703 return iter(self.T.propnames)
705 def __getitem__(self, k):
706 if k not in self.keys():
707 raise KeyError(k)
709 return getattr(self, k)
711 def __setitem__(self, k, v):
712 if k not in self.keys():
713 raise KeyError(k)
715 return setattr(self, k, v)
717 def clone(self, **kwargs):
718 '''
719 Make a copy of the object.
721 A new object of the same class is created and initialized with the
722 parameters of the object on which this method is called on. If
723 ``kwargs`` are given, these are used to override any of the
724 initialization parameters.
725 '''
727 d = dict(self)
728 for k in d:
729 v = d[k]
730 if isinstance(v, Cloneable):
731 d[k] = v.clone()
733 d.update(kwargs)
734 return self.__class__(**d)
736 @classmethod
737 def keys(cls):
738 '''
739 Get list of the source model's parameter names.
740 '''
742 return cls.T.propnames
745class STF(Object, Cloneable):
747 '''
748 Base class for source time functions.
749 '''
751 def __init__(self, effective_duration=None, **kwargs):
752 if effective_duration is not None:
753 kwargs['duration'] = effective_duration / \
754 self.factor_duration_to_effective()
756 Object.__init__(self, **kwargs)
758 @classmethod
759 def factor_duration_to_effective(cls):
760 return 1.0
762 def centroid_time(self, tref):
763 return tref
765 @property
766 def effective_duration(self):
767 return self.duration * self.factor_duration_to_effective()
769 def discretize_t(self, deltat, tref):
770 tl = math.floor(tref / deltat) * deltat
771 th = math.ceil(tref / deltat) * deltat
772 if tl == th:
773 return num.array([tl], dtype=float), num.ones(1)
774 else:
775 return (
776 num.array([tl, th], dtype=float),
777 num.array([th - tref, tref - tl], dtype=float) / deltat)
779 def base_key(self):
780 return (type(self).__name__,)
783g_unit_pulse = STF()
786def sshift(times, amplitudes, tshift, deltat):
788 t0 = math.floor(tshift / deltat) * deltat
789 t1 = math.ceil(tshift / deltat) * deltat
790 if t0 == t1:
791 return times, amplitudes
793 amplitudes2 = num.zeros(amplitudes.size + 1, dtype=float)
795 amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes
796 amplitudes2[1:] += (tshift - t0) / deltat * amplitudes
798 times2 = num.arange(times.size + 1, dtype=float) * \
799 deltat + times[0] + t0
801 return times2, amplitudes2
804class BoxcarSTF(STF):
806 '''
807 Boxcar type source time function.
809 .. figure :: /static/stf-BoxcarSTF.svg
810 :width: 40%
811 :align: center
812 :alt: boxcar source time function
813 '''
815 duration = Float.T(
816 default=0.0,
817 help='duration of the boxcar')
819 anchor = Float.T(
820 default=0.0,
821 help='anchor point with respect to source.time: ('
822 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
823 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
824 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
826 @classmethod
827 def factor_duration_to_effective(cls):
828 return 1.0
830 def centroid_time(self, tref):
831 return tref - 0.5 * self.duration * self.anchor
833 def discretize_t(self, deltat, tref):
834 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
835 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
836 tmin = round(tmin_stf / deltat) * deltat
837 tmax = round(tmax_stf / deltat) * deltat
838 nt = int(round((tmax - tmin) / deltat)) + 1
839 times = num.linspace(tmin, tmax, nt)
840 amplitudes = num.ones_like(times)
841 if times.size > 1:
842 t_edges = num.linspace(
843 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
844 t = tmin_stf + self.duration * num.array(
845 [0.0, 0.0, 1.0, 1.0], dtype=float)
846 f = num.array([0., 1., 1., 0.], dtype=float)
847 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
848 amplitudes /= num.sum(amplitudes)
850 tshift = (num.sum(amplitudes * times) - self.centroid_time(tref))
852 return sshift(times, amplitudes, -tshift, deltat)
854 def base_key(self):
855 return (type(self).__name__, self.duration, self.anchor)
858class TriangularSTF(STF):
860 '''
861 Triangular type source time function.
863 .. figure :: /static/stf-TriangularSTF.svg
864 :width: 40%
865 :align: center
866 :alt: triangular source time function
867 '''
869 duration = Float.T(
870 default=0.0,
871 help='baseline of the triangle')
873 peak_ratio = Float.T(
874 default=0.5,
875 help='fraction of time compared to duration, '
876 'when the maximum amplitude is reached')
878 anchor = Float.T(
879 default=0.0,
880 help='anchor point with respect to source.time: ('
881 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
882 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
883 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
885 @classmethod
886 def factor_duration_to_effective(cls, peak_ratio=None):
887 if peak_ratio is None:
888 peak_ratio = cls.peak_ratio.default()
890 return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0)
892 def __init__(self, effective_duration=None, **kwargs):
893 if effective_duration is not None:
894 kwargs['duration'] = effective_duration / \
895 self.factor_duration_to_effective(
896 kwargs.get('peak_ratio', None))
898 STF.__init__(self, **kwargs)
900 @property
901 def centroid_ratio(self):
902 ra = self.peak_ratio
903 rb = 1.0 - ra
904 return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb)
906 def centroid_time(self, tref):
907 ca = self.centroid_ratio
908 cb = 1.0 - ca
909 if self.anchor <= 0.:
910 return tref - ca * self.duration * self.anchor
911 else:
912 return tref - cb * self.duration * self.anchor
914 @property
915 def effective_duration(self):
916 return self.duration * self.factor_duration_to_effective(
917 self.peak_ratio)
919 def tminmax_stf(self, tref):
920 ca = self.centroid_ratio
921 cb = 1.0 - ca
922 if self.anchor <= 0.:
923 tmin_stf = tref - ca * self.duration * (self.anchor + 1.)
924 tmax_stf = tmin_stf + self.duration
925 else:
926 tmax_stf = tref + cb * self.duration * (1. - self.anchor)
927 tmin_stf = tmax_stf - self.duration
929 return tmin_stf, tmax_stf
931 def discretize_t(self, deltat, tref):
932 tmin_stf, tmax_stf = self.tminmax_stf(tref)
934 tmin = round(tmin_stf / deltat) * deltat
935 tmax = round(tmax_stf / deltat) * deltat
936 nt = int(round((tmax - tmin) / deltat)) + 1
937 if nt > 1:
938 t_edges = num.linspace(
939 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)
940 t = tmin_stf + self.duration * num.array(
941 [0.0, self.peak_ratio, 1.0], dtype=float)
942 f = num.array([0., 1., 0.], dtype=float)
943 amplitudes = util.plf_integrate_piecewise(t_edges, t, f)
944 amplitudes /= num.sum(amplitudes)
945 else:
946 amplitudes = num.ones(1)
948 times = num.linspace(tmin, tmax, nt)
949 return times, amplitudes
951 def base_key(self):
952 return (
953 type(self).__name__, self.duration, self.peak_ratio, self.anchor)
956class HalfSinusoidSTF(STF):
958 '''
959 Half sinusoid type source time function.
961 .. figure :: /static/stf-HalfSinusoidSTF.svg
962 :width: 40%
963 :align: center
964 :alt: half-sinusouid source time function
965 '''
967 duration = Float.T(
968 default=0.0,
969 help='duration of the half-sinusoid (baseline)')
971 anchor = Float.T(
972 default=0.0,
973 help='anchor point with respect to source.time: ('
974 '-1.0: left -> source duration [0, T] ~ hypocenter time, '
975 ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, '
976 '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
978 exponent = Int.T(
979 default=1,
980 help='set to 2 to use square of the half-period sinusoidal function.')
982 def __init__(self, effective_duration=None, **kwargs):
983 if effective_duration is not None:
984 kwargs['duration'] = effective_duration / \
985 self.factor_duration_to_effective(
986 kwargs.get('exponent', 1))
988 STF.__init__(self, **kwargs)
990 @classmethod
991 def factor_duration_to_effective(cls, exponent):
992 if exponent == 1:
993 return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi
994 elif exponent == 2:
995 return math.sqrt(math.pi**2 - 6) / math.pi
996 else:
997 raise ValueError('Exponent for HalfSinusoidSTF must be 1 or 2.')
999 @property
1000 def effective_duration(self):
1001 return self.duration * self.factor_duration_to_effective(self.exponent)
1003 def centroid_time(self, tref):
1004 return tref - 0.5 * self.duration * self.anchor
1006 def discretize_t(self, deltat, tref):
1007 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1008 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1009 tmin = round(tmin_stf / deltat) * deltat
1010 tmax = round(tmax_stf / deltat) * deltat
1011 nt = int(round((tmax - tmin) / deltat)) + 1
1012 if nt > 1:
1013 t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace(
1014 tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)))
1016 if self.exponent == 1:
1017 fint = -num.cos(
1018 (t_edges - tmin_stf) * (math.pi / self.duration))
1020 elif self.exponent == 2:
1021 fint = (t_edges - tmin_stf) / self.duration \
1022 - 1.0 / (2.0 * math.pi) * num.sin(
1023 (t_edges - tmin_stf) * (2.0 * math.pi / self.duration))
1024 else:
1025 raise ValueError(
1026 'Exponent for HalfSinusoidSTF must be 1 or 2.')
1028 amplitudes = fint[1:] - fint[:-1]
1029 amplitudes /= num.sum(amplitudes)
1030 else:
1031 amplitudes = num.ones(1)
1033 times = num.linspace(tmin, tmax, nt)
1034 return times, amplitudes
1036 def base_key(self):
1037 return (type(self).__name__, self.duration, self.anchor)
1040class SmoothRampSTF(STF):
1041 '''
1042 Smooth-ramp type source time function for near-field displacement.
1043 Based on moment function of double-couple point source proposed by [1]_.
1045 .. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow
1046 earthquakes from Love-wave modelling for regional distances, PEPI 32,
1047 312-324.
1049 .. figure :: /static/stf-SmoothRampSTF.svg
1050 :width: 40%
1051 :alt: smooth ramp source time function
1052 '''
1053 duration = Float.T(
1054 default=0.0,
1055 help='duration of the ramp (baseline)')
1057 rise_ratio = Float.T(
1058 default=0.5,
1059 help='fraction of time compared to duration, '
1060 'when the maximum amplitude is reached')
1062 anchor = Float.T(
1063 default=0.0,
1064 help='anchor point with respect to source.time: ('
1065 '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, '
1066 '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, '
1067 '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)')
1069 def discretize_t(self, deltat, tref):
1070 tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5
1071 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5
1072 tmin = round(tmin_stf / deltat) * deltat
1073 tmax = round(tmax_stf / deltat) * deltat
1074 D = round((tmax - tmin) / deltat) * deltat
1075 nt = int(round(D / deltat)) + 1
1076 times = num.linspace(tmin, tmax, nt)
1077 if nt > 1:
1078 rise_time = self.rise_ratio * self.duration
1079 amplitudes = num.ones_like(times)
1080 tp = tmin + rise_time
1081 ii = num.where(times <= tp)
1082 t_inc = times[ii]
1083 a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time)
1084 b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0
1085 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b)
1087 amplitudes /= num.sum(amplitudes)
1088 else:
1089 amplitudes = num.ones(1)
1091 return times, amplitudes
1093 def base_key(self):
1094 return (type(self).__name__,
1095 self.duration, self.rise_ratio, self.anchor)
1098class ResonatorSTF(STF):
1099 '''
1100 Simple resonator like source time function.
1102 .. math ::
1104 f(t) = 0 for t < 0
1105 f(t) = e^{-t/tau} * sin(2 * pi * f * t)
1108 .. figure :: /static/stf-SmoothRampSTF.svg
1109 :width: 40%
1110 :alt: smooth ramp source time function
1112 '''
1114 duration = Float.T(
1115 default=0.0,
1116 help='decay time')
1118 frequency = Float.T(
1119 default=1.0,
1120 help='resonance frequency')
1122 def discretize_t(self, deltat, tref):
1123 tmin_stf = tref
1124 tmax_stf = tref + self.duration * 3
1125 tmin = math.floor(tmin_stf / deltat) * deltat
1126 tmax = math.ceil(tmax_stf / deltat) * deltat
1127 times = util.arange2(tmin, tmax, deltat)
1128 amplitudes = num.exp(-(times - tref) / self.duration) \
1129 * num.sin(2.0 * num.pi * self.frequency * (times - tref))
1131 return times, amplitudes
1133 def base_key(self):
1134 return (type(self).__name__,
1135 self.duration, self.frequency)
1138class TremorSTF(STF):
1139 '''
1140 Oszillating source time function.
1142 .. math ::
1144 f(t) = 0 for t < -tau/2 or t > tau/2
1145 f(t) = cos(pi/tau*t) * sin(2 * pi * f * t)
1147 '''
1149 duration = Float.T(
1150 default=0.0,
1151 help='Tremor duration [s]')
1153 frequency = Float.T(
1154 default=1.0,
1155 help='Frequency [Hz]')
1157 def discretize_t(self, deltat, tref):
1158 tmin_stf = tref - 0.5 * self.duration
1159 tmax_stf = tref + 0.5 * self.duration
1160 tmin = math.floor(tmin_stf / deltat) * deltat
1161 tmax = math.ceil(tmax_stf / deltat) * deltat
1162 times = util.arange2(tmin, tmax, deltat)
1163 mask = num.logical_and(
1164 tref - 0.5 * self.duration < times,
1165 times < tref + 0.5 * self.duration)
1167 amplitudes = num.zeros_like(times)
1168 amplitudes[mask] = num.cos(num.pi/self.duration*(times[mask] - tref)) \
1169 * num.sin(2.0 * num.pi * self.frequency * (times[mask] - tref))
1170 amplitudes[mask] *= deltat
1172 return times, amplitudes
1174 def base_key(self):
1175 return (type(self).__name__,
1176 self.duration, self.frequency)
1179class SimpleLandslideSTF(STF):
1181 '''
1182 Doublepulse land-slide STF which respects conservation of momentum.
1183 '''
1185 duration_acceleration = Float.T(
1186 default=1.0,
1187 help='Duratian of the acceleration phase [s].')
1189 duration_deceleration = Float.T(
1190 default=1.0,
1191 help='Duration of the deceleration phase [s].')
1193 mute_acceleration = Bool.T(
1194 default=False,
1195 help='set acceleration to zero (for testing)')
1197 mute_deceleration = Bool.T(
1198 default=False,
1199 help='set acceleration to zero (for testing)')
1201 def discretize_t(self, deltat, tref):
1203 d_acc = self.duration_acceleration
1204 d_dec = self.duration_deceleration
1206 tmin_stf = tref
1207 tmax_stf = tref + d_acc + d_dec
1208 tmin = math.floor(tmin_stf / deltat) * deltat
1209 tmax = math.ceil(tmax_stf / deltat) * deltat
1210 times = util.arange2(tmin, tmax, deltat, epsilon=1e-3)
1212 mask_acc = num.logical_and(
1213 tref <= times,
1214 times < tref + d_acc)
1216 mask_dec = num.logical_and(
1217 tref + d_acc <= times,
1218 times < tref + d_acc + d_dec)
1220 n_acc = num.sum(mask_acc)
1221 if n_acc < 1:
1222 raise STFError(
1223 'SimpleLandslideSTF: `duration_acceleration` must be longer '
1224 'than sampling interval.')
1226 n_dec = num.sum(mask_dec)
1227 if n_dec < 1:
1228 raise STFError(
1229 'SimpleLandslideSTF: `duration_deceleration` must be longer '
1230 'than sampling interval.')
1232 amplitudes = num.zeros_like(times)
1234 amplitudes[mask_acc] = - num.sin(
1235 (times[mask_acc] - tref) / d_acc * num.pi)
1237 amplitudes[mask_dec] = num.sin(
1238 (times[mask_dec] - (tref + d_acc)) / d_dec * num.pi)
1240 sum_acc = num.abs(num.sum(amplitudes[mask_acc]))
1241 if sum_acc != 0.0:
1242 amplitudes[mask_acc] /= sum_acc
1243 else:
1244 amplitudes[mask_acc] = 1.0 / n_acc
1246 sum_dec = num.abs(num.sum(amplitudes[mask_dec]))
1247 if sum_dec != 0.0:
1248 amplitudes[mask_dec] /= sum_dec
1249 else:
1250 amplitudes[mask_dec] = 1.0 / n_dec
1252 if self.mute_acceleration:
1253 amplitudes[mask_acc] = 0.0
1255 if self.mute_deceleration:
1256 amplitudes[mask_dec] = 0.0
1258 return times, amplitudes
1261class STFMode(StringChoice):
1262 choices = ['pre', 'post']
1265class Source(Location, Cloneable):
1266 '''
1267 Base class for all source models.
1268 '''
1270 name = String.T(optional=True, default='')
1272 time = Timestamp.T(
1273 default=Timestamp.D('1970-01-01 00:00:00'),
1274 help='source origin time.')
1276 stf = STF.T(
1277 optional=True,
1278 help='source time function.')
1280 stf_mode = STFMode.T(
1281 default='post',
1282 help='whether to apply source time function in pre or '
1283 'post-processing.')
1285 def __init__(self, **kwargs):
1286 Location.__init__(self, **kwargs)
1288 def update(self, **kwargs):
1289 '''
1290 Change some of the source models parameters.
1292 Example::
1294 >>> from pyrocko import gf
1295 >>> s = gf.DCSource()
1296 >>> s.update(strike=66., dip=33.)
1297 >>> print(s)
1298 --- !pf.DCSource
1299 depth: 0.0
1300 time: 1970-01-01 00:00:00
1301 magnitude: 6.0
1302 strike: 66.0
1303 dip: 33.0
1304 rake: 0.0
1306 '''
1308 for (k, v) in kwargs.items():
1309 self[k] = v
1311 def grid(self, **variables):
1312 '''
1313 Create grid of source model variations.
1315 :returns: :py:class:`SourceGrid` instance.
1317 Example::
1319 >>> from pyrocko import gf
1320 >>> base = DCSource()
1321 >>> R = gf.Range
1322 >>> for s in base.grid(R('
1324 '''
1325 return SourceGrid(base=self, variables=variables)
1327 def base_key(self):
1328 '''
1329 Get key to decide about source discretization / GF stack sharing.
1331 When two source models differ only in amplitude and origin time, the
1332 discretization and the GF stacking can be done only once for a unit
1333 amplitude and a zero origin time and the amplitude and origin times of
1334 the seismograms can be applied during post-processing of the synthetic
1335 seismogram.
1337 For any derived parameterized source model, this method is called to
1338 decide if discretization and stacking of the source should be shared.
1339 When two source models return an equal vector of values discretization
1340 is shared.
1341 '''
1342 return (self.depth, self.lat, self.north_shift,
1343 self.lon, self.east_shift, self.time, type(self).__name__) + \
1344 self.effective_stf_pre().base_key()
1346 def get_factor(self):
1347 '''
1348 Get the scaling factor to be applied during post-processing.
1350 Discretization of the base seismogram is usually done for a unit
1351 amplitude, because a common factor can be efficiently multiplied to
1352 final seismograms. This eliminates to do repeat the stacking when
1353 creating seismograms for a series of source models only differing in
1354 amplitude.
1356 This method should return the scaling factor to apply in the
1357 post-processing (often this is simply the scalar moment of the source).
1358 '''
1360 return 1.0
1362 def effective_stf_pre(self):
1363 '''
1364 Return the STF applied before stacking of the Green's functions.
1366 This STF is used during discretization of the parameterized source
1367 models, i.e. to produce a temporal distribution of point sources.
1369 Handling of the STF before stacking of the GFs is less efficient but
1370 allows to use different source time functions for different parts of
1371 the source.
1372 '''
1374 if self.stf is not None and self.stf_mode == 'pre':
1375 return self.stf
1376 else:
1377 return g_unit_pulse
1379 def effective_stf_post(self):
1380 '''
1381 Return the STF applied after stacking of the Green's fuctions.
1383 This STF is used in the post-processing of the synthetic seismograms.
1385 Handling of the STF after stacking of the GFs is usually more efficient
1386 but is only possible when a common STF is used for all subsources.
1387 '''
1389 if self.stf is not None and self.stf_mode == 'post':
1390 return self.stf
1391 else:
1392 return g_unit_pulse
1394 def _dparams_base(self):
1395 return dict(times=arr(self.time),
1396 lat=self.lat, lon=self.lon,
1397 north_shifts=arr(self.north_shift),
1398 east_shifts=arr(self.east_shift),
1399 depths=arr(self.depth))
1401 def _hash(self):
1402 sha = sha1()
1403 for k in self.base_key():
1404 sha.update(str(k).encode())
1405 return sha.hexdigest()
1407 def _dparams_base_repeated(self, times):
1408 if times is None:
1409 return self._dparams_base()
1411 nt = times.size
1412 north_shifts = num.repeat(self.north_shift, nt)
1413 east_shifts = num.repeat(self.east_shift, nt)
1414 depths = num.repeat(self.depth, nt)
1415 return dict(times=times,
1416 lat=self.lat, lon=self.lon,
1417 north_shifts=north_shifts,
1418 east_shifts=east_shifts,
1419 depths=depths)
1421 def pyrocko_event(self, store=None, target=None, **kwargs):
1422 duration = None
1423 if self.stf:
1424 duration = self.stf.effective_duration
1426 return model.Event(
1427 lat=self.lat,
1428 lon=self.lon,
1429 north_shift=self.north_shift,
1430 east_shift=self.east_shift,
1431 time=self.time,
1432 name=self.name,
1433 depth=self.depth,
1434 duration=duration,
1435 **kwargs)
1437 def geometry(self, **kwargs):
1438 raise NotImplementedError
1440 def outline(self, cs='xyz'):
1441 points = num.atleast_2d(num.zeros([1, 3]))
1443 points[:, 0] += self.north_shift
1444 points[:, 1] += self.east_shift
1445 points[:, 2] += self.depth
1446 if cs == 'xyz':
1447 return points
1448 elif cs == 'xy':
1449 return points[:, :2]
1450 elif cs in ('latlon', 'lonlat'):
1451 latlon = ne_to_latlon(
1452 self.lat, self.lon, points[:, 0], points[:, 1])
1454 latlon = num.array(latlon).T
1455 if cs == 'latlon':
1456 return latlon
1457 else:
1458 return latlon[:, ::-1]
1460 @classmethod
1461 def from_pyrocko_event(cls, ev, **kwargs):
1462 if ev.depth is None:
1463 raise ConversionError(
1464 'Cannot convert event object to source object: '
1465 'no depth information available')
1467 stf = None
1468 if ev.duration is not None:
1469 stf = HalfSinusoidSTF(effective_duration=ev.duration)
1471 d = dict(
1472 name=ev.name,
1473 time=ev.time,
1474 lat=ev.lat,
1475 lon=ev.lon,
1476 north_shift=ev.north_shift,
1477 east_shift=ev.east_shift,
1478 depth=ev.depth,
1479 stf=stf)
1480 d.update(kwargs)
1481 return cls(**d)
1483 def get_magnitude(self):
1484 raise NotImplementedError(
1485 '%s does not implement get_magnitude()'
1486 % self.__class__.__name__)
1489class SourceWithMagnitude(Source):
1490 '''
1491 Base class for sources containing a moment magnitude.
1492 '''
1494 magnitude = Float.T(
1495 default=6.0,
1496 help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1498 def __init__(self, **kwargs):
1499 if 'moment' in kwargs:
1500 mom = kwargs.pop('moment')
1501 if 'magnitude' not in kwargs:
1502 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1504 Source.__init__(self, **kwargs)
1506 @property
1507 def moment(self):
1508 return float(pmt.magnitude_to_moment(self.magnitude))
1510 @moment.setter
1511 def moment(self, value):
1512 self.magnitude = float(pmt.moment_to_magnitude(value))
1514 def pyrocko_event(self, store=None, target=None, **kwargs):
1515 return Source.pyrocko_event(
1516 self, store, target,
1517 magnitude=self.magnitude,
1518 **kwargs)
1520 @classmethod
1521 def from_pyrocko_event(cls, ev, **kwargs):
1522 d = {}
1523 if ev.magnitude:
1524 d.update(magnitude=ev.magnitude)
1526 d.update(kwargs)
1527 return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
1529 def get_magnitude(self):
1530 return self.magnitude
1533class DerivedMagnitudeError(ValidationError):
1534 '''
1535 Raised when conversion between magnitude, moment, volume change or
1536 displacement failed.
1537 '''
1538 pass
1541class SourceWithDerivedMagnitude(Source):
1543 class __T(Source.T):
1545 def validate_extra(self, val):
1546 Source.T.validate_extra(self, val)
1547 val.check_conflicts()
1549 def check_conflicts(self):
1550 '''
1551 Check for parameter conflicts.
1553 To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError`
1554 on conflicts.
1555 '''
1556 pass
1558 def get_magnitude(self, store=None, target=None):
1559 raise DerivedMagnitudeError('No magnitude set.')
1561 def get_moment(self, store=None, target=None):
1562 return float(pmt.magnitude_to_moment(
1563 self.get_magnitude(store, target)))
1565 def pyrocko_moment_tensor(self, store=None, target=None):
1566 raise NotImplementedError(
1567 '%s does not implement pyrocko_moment_tensor()'
1568 % self.__class__.__name__)
1570 def pyrocko_event(self, store=None, target=None, **kwargs):
1571 try:
1572 mt = self.pyrocko_moment_tensor(store, target)
1573 magnitude = self.get_magnitude()
1574 except (DerivedMagnitudeError, NotImplementedError):
1575 mt = None
1576 magnitude = None
1578 return Source.pyrocko_event(
1579 self, store, target,
1580 moment_tensor=mt,
1581 magnitude=magnitude,
1582 **kwargs)
1585class ExplosionSource(SourceWithDerivedMagnitude):
1586 '''
1587 An isotropic explosion point source.
1588 '''
1590 magnitude = Float.T(
1591 optional=True,
1592 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
1594 volume_change = Float.T(
1595 optional=True,
1596 help='volume change of the explosion/implosion or '
1597 'the contracting/extending magmatic source. [m^3]')
1599 discretized_source_class = meta.DiscretizedExplosionSource
1601 def __init__(self, **kwargs):
1602 if 'moment' in kwargs:
1603 mom = kwargs.pop('moment')
1604 if 'magnitude' not in kwargs:
1605 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
1607 SourceWithDerivedMagnitude.__init__(self, **kwargs)
1609 def base_key(self):
1610 return SourceWithDerivedMagnitude.base_key(self) + \
1611 (self.volume_change,)
1613 def check_conflicts(self):
1614 if self.magnitude is not None and self.volume_change is not None:
1615 raise DerivedMagnitudeError(
1616 'Magnitude and volume_change are both defined.')
1618 def get_magnitude(self, store=None, target=None):
1619 self.check_conflicts()
1621 if self.magnitude is not None:
1622 return self.magnitude
1624 elif self.volume_change is not None:
1625 moment = self.volume_change * \
1626 self.get_moment_to_volume_change_ratio(store, target)
1628 return float(pmt.moment_to_magnitude(abs(moment)))
1629 else:
1630 return float(pmt.moment_to_magnitude(1.0))
1632 def get_volume_change(self, store=None, target=None):
1633 self.check_conflicts()
1635 if self.volume_change is not None:
1636 return self.volume_change
1638 elif self.magnitude is not None:
1639 moment = float(pmt.magnitude_to_moment(self.magnitude))
1640 return moment / self.get_moment_to_volume_change_ratio(
1641 store, target)
1643 else:
1644 return 1.0 / self.get_moment_to_volume_change_ratio(store)
1646 def get_moment_to_volume_change_ratio(self, store, target=None):
1647 if store is None:
1648 raise DerivedMagnitudeError(
1649 'Need earth model to convert between volume change and '
1650 'magnitude.')
1652 points = num.array(
1653 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1655 interpolation = target.interpolation if target else 'multilinear'
1656 try:
1657 shear_moduli = store.config.get_shear_moduli(
1658 self.lat, self.lon,
1659 points=points,
1660 interpolation=interpolation)[0]
1662 bulk_moduli = store.config.get_bulk_moduli(
1663 self.lat, self.lon,
1664 points=points,
1665 interpolation=interpolation)[0]
1666 except meta.OutOfBounds:
1667 raise DerivedMagnitudeError(
1668 'Could not get shear modulus at source position.')
1670 return float(2. * shear_moduli + bulk_moduli)
1672 def get_factor(self):
1673 return 1.0
1675 def discretize_basesource(self, store, target=None):
1676 times, amplitudes = self.effective_stf_pre().discretize_t(
1677 store.config.deltat, self.time)
1679 amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.)
1681 if self.volume_change is not None:
1682 if self.volume_change < 0.:
1683 amplitudes *= -1
1685 return meta.DiscretizedExplosionSource(
1686 m0s=amplitudes,
1687 **self._dparams_base_repeated(times))
1689 def pyrocko_moment_tensor(self, store=None, target=None):
1690 a = self.get_moment(store, target) * math.sqrt(2. / 3.)
1691 return pmt.MomentTensor(m=pmt.symmat6(a, a, a, 0., 0., 0.))
1694class RectangularExplosionSource(ExplosionSource):
1695 '''
1696 Rectangular or line explosion source.
1697 '''
1699 discretized_source_class = meta.DiscretizedExplosionSource
1701 strike = Float.T(
1702 default=0.0,
1703 help='strike direction in [deg], measured clockwise from north')
1705 dip = Float.T(
1706 default=90.0,
1707 help='dip angle in [deg], measured downward from horizontal')
1709 length = Float.T(
1710 default=0.,
1711 help='length of rectangular source area [m]')
1713 width = Float.T(
1714 default=0.,
1715 help='width of rectangular source area [m]')
1717 anchor = StringChoice.T(
1718 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
1719 'bottom_left', 'bottom_right'],
1720 default='center',
1721 optional=True,
1722 help='Anchor point for positioning the plane, can be: top, center or'
1723 'bottom and also top_left, top_right,bottom_left,'
1724 'bottom_right, center_left and center right')
1726 nucleation_x = Float.T(
1727 optional=True,
1728 help='horizontal position of rupture nucleation in normalized fault '
1729 'plane coordinates (-1 = left edge, +1 = right edge)')
1731 nucleation_y = Float.T(
1732 optional=True,
1733 help='down-dip position of rupture nucleation in normalized fault '
1734 'plane coordinates (-1 = upper edge, +1 = lower edge)')
1736 velocity = Float.T(
1737 default=3500.,
1738 help='speed of explosion front [m/s]')
1740 aggressive_oversampling = Bool.T(
1741 default=False,
1742 help='Aggressive oversampling for basesource discretization. '
1743 "When using 'multilinear' interpolation oversampling has"
1744 ' practically no effect.')
1746 def base_key(self):
1747 return Source.base_key(self) + (self.strike, self.dip, self.length,
1748 self.width, self.nucleation_x,
1749 self.nucleation_y, self.velocity,
1750 self.anchor)
1752 def discretize_basesource(self, store, target=None):
1754 if self.nucleation_x is not None:
1755 nucx = self.nucleation_x * 0.5 * self.length
1756 else:
1757 nucx = None
1759 if self.nucleation_y is not None:
1760 nucy = self.nucleation_y * 0.5 * self.width
1761 else:
1762 nucy = None
1764 stf = self.effective_stf_pre()
1766 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
1767 store.config.deltas, store.config.deltat,
1768 self.time, self.north_shift, self.east_shift, self.depth,
1769 self.strike, self.dip, self.length, self.width, self.anchor,
1770 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy)
1772 amplitudes /= num.sum(amplitudes)
1773 amplitudes *= self.get_moment(store, target)
1775 return meta.DiscretizedExplosionSource(
1776 lat=self.lat,
1777 lon=self.lon,
1778 times=times,
1779 north_shifts=points[:, 0],
1780 east_shifts=points[:, 1],
1781 depths=points[:, 2],
1782 m0s=amplitudes)
1784 def outline(self, cs='xyz'):
1785 points = outline_rect_source(self.strike, self.dip, self.length,
1786 self.width, self.anchor)
1788 points[:, 0] += self.north_shift
1789 points[:, 1] += self.east_shift
1790 points[:, 2] += self.depth
1791 if cs == 'xyz':
1792 return points
1793 elif cs == 'xy':
1794 return points[:, :2]
1795 elif cs in ('latlon', 'lonlat'):
1796 latlon = ne_to_latlon(
1797 self.lat, self.lon, points[:, 0], points[:, 1])
1799 latlon = num.array(latlon).T
1800 if cs == 'latlon':
1801 return latlon
1802 else:
1803 return latlon[:, ::-1]
1805 def get_nucleation_abs_coord(self, cs='xy'):
1807 if self.nucleation_x is None:
1808 return None, None
1810 coords = from_plane_coords(self.strike, self.dip, self.length,
1811 self.width, self.depth, self.nucleation_x,
1812 self.nucleation_y, lat=self.lat,
1813 lon=self.lon, north_shift=self.north_shift,
1814 east_shift=self.east_shift, cs=cs)
1815 return coords
1818class DCSource(SourceWithMagnitude):
1819 '''
1820 A double-couple point source.
1821 '''
1823 strike = Float.T(
1824 default=0.0,
1825 help='strike direction in [deg], measured clockwise from north')
1827 dip = Float.T(
1828 default=90.0,
1829 help='dip angle in [deg], measured downward from horizontal')
1831 rake = Float.T(
1832 default=0.0,
1833 help='rake angle in [deg], '
1834 'measured counter-clockwise from right-horizontal '
1835 'in on-plane view')
1837 discretized_source_class = meta.DiscretizedMTSource
1839 def base_key(self):
1840 return Source.base_key(self) + (self.strike, self.dip, self.rake)
1842 def get_factor(self):
1843 return float(pmt.magnitude_to_moment(self.magnitude))
1845 def discretize_basesource(self, store, target=None):
1846 mot = pmt.MomentTensor(
1847 strike=self.strike, dip=self.dip, rake=self.rake)
1849 times, amplitudes = self.effective_stf_pre().discretize_t(
1850 store.config.deltat, self.time)
1851 return meta.DiscretizedMTSource(
1852 m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis],
1853 **self._dparams_base_repeated(times))
1855 def pyrocko_moment_tensor(self, store=None, target=None):
1856 return pmt.MomentTensor(
1857 strike=self.strike,
1858 dip=self.dip,
1859 rake=self.rake,
1860 scalar_moment=self.moment)
1862 def pyrocko_event(self, store=None, target=None, **kwargs):
1863 return SourceWithMagnitude.pyrocko_event(
1864 self, store, target,
1865 moment_tensor=self.pyrocko_moment_tensor(store, target),
1866 **kwargs)
1868 @classmethod
1869 def from_pyrocko_event(cls, ev, **kwargs):
1870 d = {}
1871 mt = ev.moment_tensor
1872 if mt:
1873 (strike, dip, rake), _ = mt.both_strike_dip_rake()
1874 d.update(
1875 strike=float(strike),
1876 dip=float(dip),
1877 rake=float(rake),
1878 magnitude=float(mt.moment_magnitude()))
1880 d.update(kwargs)
1881 return super(DCSource, cls).from_pyrocko_event(ev, **d)
1884class CLVDSource(SourceWithMagnitude):
1885 '''
1886 A pure CLVD point source.
1887 '''
1889 discretized_source_class = meta.DiscretizedMTSource
1891 azimuth = Float.T(
1892 default=0.0,
1893 help='azimuth direction of largest dipole, clockwise from north [deg]')
1895 dip = Float.T(
1896 default=90.,
1897 help='dip direction of largest dipole, downward from horizontal [deg]')
1899 def base_key(self):
1900 return Source.base_key(self) + (self.azimuth, self.dip)
1902 def get_factor(self):
1903 return float(pmt.magnitude_to_moment(self.magnitude))
1905 @property
1906 def m6(self):
1907 a = math.sqrt(4. / 3.) * self.get_factor()
1908 m = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
1909 rotmat1 = pmt.euler_to_matrix(
1910 d2r * (self.dip - 90.),
1911 d2r * (self.azimuth - 90.),
1912 0.)
1913 m = num.dot(rotmat1.T, num.dot(m, rotmat1))
1914 return pmt.to6(m)
1916 @property
1917 def m6_astuple(self):
1918 return tuple(self.m6.tolist())
1920 def discretize_basesource(self, store, target=None):
1921 factor = self.get_factor()
1922 times, amplitudes = self.effective_stf_pre().discretize_t(
1923 store.config.deltat, self.time)
1924 return meta.DiscretizedMTSource(
1925 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor,
1926 **self._dparams_base_repeated(times))
1928 def pyrocko_moment_tensor(self, store=None, target=None):
1929 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
1931 def pyrocko_event(self, store=None, target=None, **kwargs):
1932 mt = self.pyrocko_moment_tensor(store, target)
1933 return Source.pyrocko_event(
1934 self, store, target,
1935 moment_tensor=self.pyrocko_moment_tensor(store, target),
1936 magnitude=float(mt.moment_magnitude()),
1937 **kwargs)
1940class VLVDSource(SourceWithDerivedMagnitude):
1941 '''
1942 Volumetric linear vector dipole source.
1944 This source is a parameterization for a restricted moment tensor point
1945 source, useful to represent dyke or sill like inflation or deflation
1946 sources. The restriction is such that the moment tensor is rotational
1947 symmetric. It can be represented by a superposition of a linear vector
1948 dipole (here we use a CLVD for convenience) and an isotropic component. The
1949 restricted moment tensor has 4 degrees of freedom: 2 independent
1950 eigenvalues and 2 rotation angles orienting the the symmetry axis.
1952 In this parameterization, the isotropic component is controlled by
1953 ``volume_change``. To define the moment tensor, it must be converted to the
1954 scalar moment of the the MT's isotropic component. For the conversion, the
1955 shear modulus at the source's position must be known. This value is
1956 extracted from the earth model defined in the GF store in use.
1958 The CLVD part by controlled by its scalar moment :math:`M_0`:
1959 ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a
1960 positiv or negativ CLVD (the sign of the largest eigenvalue).
1961 '''
1963 discretized_source_class = meta.DiscretizedMTSource
1965 azimuth = Float.T(
1966 default=0.0,
1967 help='azimuth direction of symmetry axis, clockwise from north [deg].')
1969 dip = Float.T(
1970 default=90.,
1971 help='dip direction of symmetry axis, downward from horizontal [deg].')
1973 volume_change = Float.T(
1974 default=0.,
1975 help='volume change of the inflation/deflation [m^3].')
1977 clvd_moment = Float.T(
1978 default=0.,
1979 help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign '
1980 'controls the sign of the CLVD (the sign of its largest '
1981 'eigenvalue).')
1983 def get_moment_to_volume_change_ratio(self, store, target):
1984 if store is None or target is None:
1985 raise DerivedMagnitudeError(
1986 'Need earth model to convert between volume change and '
1987 'magnitude.')
1989 points = num.array(
1990 [[self.north_shift, self.east_shift, self.depth]], dtype=float)
1992 try:
1993 shear_moduli = store.config.get_shear_moduli(
1994 self.lat, self.lon,
1995 points=points,
1996 interpolation=target.interpolation)[0]
1998 bulk_moduli = store.config.get_bulk_moduli(
1999 self.lat, self.lon,
2000 points=points,
2001 interpolation=target.interpolation)[0]
2002 except meta.OutOfBounds:
2003 raise DerivedMagnitudeError(
2004 'Could not get shear modulus at source position.')
2006 return float(2. * shear_moduli + bulk_moduli)
2008 def base_key(self):
2009 return Source.base_key(self) + \
2010 (self.azimuth, self.dip, self.volume_change, self.clvd_moment)
2012 def get_magnitude(self, store=None, target=None):
2013 mt = self.pyrocko_moment_tensor(store, target)
2014 return float(pmt.moment_to_magnitude(mt.moment))
2016 def get_m6(self, store, target):
2017 a = math.sqrt(4. / 3.) * self.clvd_moment
2018 m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
2020 rotmat1 = pmt.euler_to_matrix(
2021 d2r * (self.dip - 90.),
2022 d2r * (self.azimuth - 90.),
2023 0.)
2024 m_clvd = num.dot(rotmat1.T, num.dot(m_clvd, rotmat1))
2026 m_iso = self.volume_change * \
2027 self.get_moment_to_volume_change_ratio(store, target)
2029 m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0.,
2030 0., 0.,) * math.sqrt(2. / 3)
2032 m = pmt.to6(m_clvd) + pmt.to6(m_iso)
2033 return m
2035 def get_moment(self, store=None, target=None):
2036 return float(pmt.magnitude_to_moment(
2037 self.get_magnitude(store, target)))
2039 def get_m6_astuple(self, store, target):
2040 m6 = self.get_m6(store, target)
2041 return tuple(m6.tolist())
2043 def discretize_basesource(self, store, target=None):
2044 times, amplitudes = self.effective_stf_pre().discretize_t(
2045 store.config.deltat, self.time)
2047 m6 = self.get_m6(store, target)
2048 m6 *= amplitudes / self.get_factor()
2050 return meta.DiscretizedMTSource(
2051 m6s=m6[num.newaxis, :],
2052 **self._dparams_base_repeated(times))
2054 def pyrocko_moment_tensor(self, store=None, target=None):
2055 m6_astuple = self.get_m6_astuple(store, target)
2056 return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple))
2059class MTSource(Source):
2060 '''
2061 A moment tensor point source.
2062 '''
2064 discretized_source_class = meta.DiscretizedMTSource
2066 mnn = Float.T(
2067 default=1.,
2068 help='north-north component of moment tensor in [Nm]')
2070 mee = Float.T(
2071 default=1.,
2072 help='east-east component of moment tensor in [Nm]')
2074 mdd = Float.T(
2075 default=1.,
2076 help='down-down component of moment tensor in [Nm]')
2078 mne = Float.T(
2079 default=0.,
2080 help='north-east component of moment tensor in [Nm]')
2082 mnd = Float.T(
2083 default=0.,
2084 help='north-down component of moment tensor in [Nm]')
2086 med = Float.T(
2087 default=0.,
2088 help='east-down component of moment tensor in [Nm]')
2090 def __init__(self, **kwargs):
2091 if 'm6' in kwargs:
2092 for (k, v) in zip('mnn mee mdd mne mnd med'.split(),
2093 kwargs.pop('m6')):
2094 kwargs[k] = float(v)
2096 Source.__init__(self, **kwargs)
2098 @property
2099 def m6(self):
2100 return num.array(self.m6_astuple)
2102 @property
2103 def m6_astuple(self):
2104 return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)
2106 @m6.setter
2107 def m6(self, value):
2108 self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value
2110 def base_key(self):
2111 return Source.base_key(self) + self.m6_astuple
2113 def discretize_basesource(self, store, target=None):
2114 times, amplitudes = self.effective_stf_pre().discretize_t(
2115 store.config.deltat, self.time)
2116 return meta.DiscretizedMTSource(
2117 m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis],
2118 **self._dparams_base_repeated(times))
2120 def get_magnitude(self, store=None, target=None):
2121 m6 = self.m6
2122 return pmt.moment_to_magnitude(
2123 math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) /
2124 math.sqrt(2.))
2126 def pyrocko_moment_tensor(self, store=None, target=None):
2127 return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
2129 def pyrocko_event(self, store=None, target=None, **kwargs):
2130 mt = self.pyrocko_moment_tensor(store, target)
2131 return Source.pyrocko_event(
2132 self, store, target,
2133 moment_tensor=self.pyrocko_moment_tensor(store, target),
2134 magnitude=float(mt.moment_magnitude()),
2135 **kwargs)
2137 @classmethod
2138 def from_pyrocko_event(cls, ev, **kwargs):
2139 d = {}
2140 mt = ev.moment_tensor
2141 if mt:
2142 d.update(m6=tuple(map(float, mt.m6())))
2143 else:
2144 if ev.magnitude is not None:
2145 mom = pmt.magnitude_to_moment(ev.magnitude)
2146 v = math.sqrt(2. / 3.) * mom
2147 d.update(m6=(v, v, v, 0., 0., 0.))
2149 d.update(kwargs)
2150 return super(MTSource, cls).from_pyrocko_event(ev, **d)
2153map_anchor = {
2154 'center': (0.0, 0.0),
2155 'center_left': (-1.0, 0.0),
2156 'center_right': (1.0, 0.0),
2157 'top': (0.0, -1.0),
2158 'top_left': (-1.0, -1.0),
2159 'top_right': (1.0, -1.0),
2160 'bottom': (0.0, 1.0),
2161 'bottom_left': (-1.0, 1.0),
2162 'bottom_right': (1.0, 1.0)}
2165class RectangularSource(SourceWithDerivedMagnitude):
2166 '''
2167 Classical Haskell source model modified for bilateral rupture.
2168 '''
2170 discretized_source_class = meta.DiscretizedMTSource
2172 magnitude = Float.T(
2173 optional=True,
2174 help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
2176 strike = Float.T(
2177 default=0.0,
2178 help='strike direction in [deg], measured clockwise from north')
2180 dip = Float.T(
2181 default=90.0,
2182 help='dip angle in [deg], measured downward from horizontal')
2184 rake = Float.T(
2185 default=0.0,
2186 help='rake angle in [deg], '
2187 'measured counter-clockwise from right-horizontal '
2188 'in on-plane view')
2190 length = Float.T(
2191 default=0.,
2192 help='length of rectangular source area [m]')
2194 width = Float.T(
2195 default=0.,
2196 help='width of rectangular source area [m]')
2198 anchor = StringChoice.T(
2199 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2200 'bottom_left', 'bottom_right'],
2201 default='center',
2202 optional=True,
2203 help='Anchor point for positioning the plane, can be: ``top, center '
2204 'bottom, top_left, top_right,bottom_left,'
2205 'bottom_right, center_left, center right``.')
2207 nucleation_x = Float.T(
2208 optional=True,
2209 help='horizontal position of rupture nucleation in normalized fault '
2210 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge)')
2212 nucleation_y = Float.T(
2213 optional=True,
2214 help='down-dip position of rupture nucleation in normalized fault '
2215 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge)')
2217 velocity = Float.T(
2218 default=3500.,
2219 help='speed of rupture front [m/s]')
2221 slip = Float.T(
2222 optional=True,
2223 help='Slip on the rectangular source area [m]')
2225 opening_fraction = Float.T(
2226 default=0.,
2227 help='Determines fraction of slip related to opening. '
2228 '(``-1``: pure tensile closing, '
2229 '``0``: pure shear, '
2230 '``1``: pure tensile opening)')
2232 decimation_factor = Int.T(
2233 optional=True,
2234 default=1,
2235 help='Sub-source decimation factor, a larger decimation will'
2236 ' make the result inaccurate but shorten the necessary'
2237 ' computation time (use for testing puposes only).')
2239 aggressive_oversampling = Bool.T(
2240 default=False,
2241 help='Aggressive oversampling for basesource discretization. '
2242 "When using 'multilinear' interpolation oversampling has"
2243 ' practically no effect.')
2245 def __init__(self, **kwargs):
2246 if 'moment' in kwargs:
2247 mom = kwargs.pop('moment')
2248 if 'magnitude' not in kwargs:
2249 kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
2251 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2253 def base_key(self):
2254 return SourceWithDerivedMagnitude.base_key(self) + (
2255 self.magnitude,
2256 self.slip,
2257 self.strike,
2258 self.dip,
2259 self.rake,
2260 self.length,
2261 self.width,
2262 self.nucleation_x,
2263 self.nucleation_y,
2264 self.velocity,
2265 self.decimation_factor,
2266 self.anchor)
2268 def check_conflicts(self):
2269 if self.magnitude is not None and self.slip is not None:
2270 raise DerivedMagnitudeError(
2271 'Magnitude and slip are both defined.')
2273 def get_magnitude(self, store=None, target=None):
2274 self.check_conflicts()
2275 if self.magnitude is not None:
2276 return self.magnitude
2278 elif self.slip is not None:
2279 if None in (store, target):
2280 raise DerivedMagnitudeError(
2281 'Magnitude for a rectangular source with slip defined '
2282 'can only be derived when earth model and target '
2283 'interpolation method are available.')
2285 amplitudes = self._discretize(store, target)[2]
2286 if amplitudes.ndim == 2:
2287 # CLVD component has no net moment, leave out
2288 return float(pmt.moment_to_magnitude(
2289 num.sum(num.abs(amplitudes[0:2, :]).sum())))
2290 else:
2291 return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
2293 else:
2294 return float(pmt.moment_to_magnitude(1.0))
2296 def get_factor(self):
2297 return 1.0
2299 def get_slip_tensile(self):
2300 return self.slip * self.opening_fraction
2302 def get_slip_shear(self):
2303 return self.slip - abs(self.get_slip_tensile)
2305 def _discretize(self, store, target):
2306 if self.nucleation_x is not None:
2307 nucx = self.nucleation_x * 0.5 * self.length
2308 else:
2309 nucx = None
2311 if self.nucleation_y is not None:
2312 nucy = self.nucleation_y * 0.5 * self.width
2313 else:
2314 nucy = None
2316 stf = self.effective_stf_pre()
2318 points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source(
2319 store.config.deltas, store.config.deltat,
2320 self.time, self.north_shift, self.east_shift, self.depth,
2321 self.strike, self.dip, self.length, self.width, self.anchor,
2322 self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy,
2323 decimation_factor=self.decimation_factor,
2324 aggressive_oversampling=self.aggressive_oversampling)
2326 if self.slip is not None:
2327 if target is not None:
2328 interpolation = target.interpolation
2329 else:
2330 interpolation = 'nearest_neighbor'
2331 logger.warning(
2332 'no target information available, will use '
2333 '"nearest_neighbor" interpolation when extracting shear '
2334 'modulus from earth model')
2336 shear_moduli = store.config.get_shear_moduli(
2337 self.lat, self.lon,
2338 points=points,
2339 interpolation=interpolation)
2341 tensile_slip = self.get_slip_tensile()
2342 shear_slip = self.slip - abs(tensile_slip)
2344 amplitudes_total = [shear_moduli * shear_slip]
2345 if tensile_slip != 0:
2346 bulk_moduli = store.config.get_bulk_moduli(
2347 self.lat, self.lon,
2348 points=points,
2349 interpolation=interpolation)
2351 tensile_iso = bulk_moduli * tensile_slip
2352 tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip
2354 amplitudes_total.extend([tensile_iso, tensile_clvd])
2356 amplitudes_total = num.vstack(amplitudes_total).squeeze() * \
2357 amplitudes * dl * dw
2359 else:
2360 # normalization to retain total moment
2361 amplitudes_norm = amplitudes / num.sum(amplitudes)
2362 moment = self.get_moment(store, target)
2364 amplitudes_total = [
2365 amplitudes_norm * moment * (1 - abs(self.opening_fraction))]
2366 if self.opening_fraction != 0.:
2367 amplitudes_total.append(
2368 amplitudes_norm * self.opening_fraction * moment)
2370 amplitudes_total = num.vstack(amplitudes_total).squeeze()
2372 return points, times, num.atleast_1d(amplitudes_total), dl, dw, nl, nw
2374 def discretize_basesource(self, store, target=None):
2376 points, times, amplitudes, dl, dw, nl, nw = self._discretize(
2377 store, target)
2379 mot = pmt.MomentTensor(
2380 strike=self.strike, dip=self.dip, rake=self.rake)
2381 m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0)
2383 if amplitudes.ndim == 1:
2384 m6s[:, :] *= amplitudes[:, num.newaxis]
2385 elif amplitudes.ndim == 2:
2386 # shear MT components
2387 rotmat1 = pmt.euler_to_matrix(
2388 d2r * self.dip, d2r * self.strike, d2r * -self.rake)
2389 m6s[:, :] *= amplitudes[0, :][:, num.newaxis]
2391 if amplitudes.shape[0] == 2:
2392 # tensile MT components - moment/magnitude input
2393 tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.)
2394 rot_tensile = pmt.to6(
2395 num.dot(rotmat1.T, num.dot(tensile, rotmat1)))
2397 m6s_tensile = rot_tensile[
2398 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2399 m6s += m6s_tensile
2401 elif amplitudes.shape[0] == 3:
2402 # tensile MT components - slip input
2403 iso = pmt.symmat6(1., 1., 1., 0., 0., 0.)
2404 clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.)
2406 rot_iso = pmt.to6(
2407 num.dot(rotmat1.T, num.dot(iso, rotmat1)))
2408 rot_clvd = pmt.to6(
2409 num.dot(rotmat1.T, num.dot(clvd, rotmat1)))
2411 m6s_iso = rot_iso[
2412 num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
2413 m6s_clvd = rot_clvd[
2414 num.newaxis, :] * amplitudes[2, :][:, num.newaxis]
2415 m6s += m6s_iso + m6s_clvd
2416 else:
2417 raise ValueError('Unknwown amplitudes shape!')
2418 else:
2419 raise ValueError(
2420 'Unexpected dimension of {}'.format(amplitudes.ndim))
2422 ds = meta.DiscretizedMTSource(
2423 lat=self.lat,
2424 lon=self.lon,
2425 times=times,
2426 north_shifts=points[:, 0],
2427 east_shifts=points[:, 1],
2428 depths=points[:, 2],
2429 m6s=m6s,
2430 dl=dl,
2431 dw=dw,
2432 nl=nl,
2433 nw=nw)
2435 return ds
2437 def xy_to_coord(self, x, y, cs='xyz'):
2438 ln, wd = self.length, self.width
2439 strike, dip = self.strike, self.dip
2441 def array_check(variable):
2442 if not isinstance(variable, num.ndarray):
2443 return num.array(variable)
2444 else:
2445 return variable
2447 x, y = array_check(x), array_check(y)
2449 if x.shape[0] != y.shape[0]:
2450 raise ValueError('Shapes of x and y mismatch')
2452 x, y = x * 0.5 * ln, y * 0.5 * wd
2454 points = num.hstack((
2455 x.reshape(-1, 1), y.reshape(-1, 1), num.zeros((x.shape[0], 1))))
2457 anch_x, anch_y = map_anchor[self.anchor]
2458 points[:, 0] -= anch_x * 0.5 * ln
2459 points[:, 1] -= anch_y * 0.5 * wd
2461 rotmat = num.asarray(
2462 pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
2464 points_rot = num.dot(rotmat.T, points.T).T
2466 points_rot[:, 0] += self.north_shift
2467 points_rot[:, 1] += self.east_shift
2468 points_rot[:, 2] += self.depth
2470 if cs == 'xyz':
2471 return points_rot
2472 elif cs == 'xy':
2473 return points_rot[:, :2]
2474 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2475 latlon = ne_to_latlon(
2476 self.lat, self.lon, points_rot[:, 0], points_rot[:, 1])
2477 latlon = num.array(latlon).T
2478 if cs == 'latlon':
2479 return latlon
2480 elif cs == 'lonlat':
2481 return latlon[:, ::-1]
2482 else:
2483 return num.concatenate(
2484 (latlon, points_rot[:, 2].reshape((len(points_rot), 1))),
2485 axis=1)
2487 def outline(self, cs='xyz'):
2488 x = num.array([-1., 1., 1., -1., -1.])
2489 y = num.array([-1., -1., 1., 1., -1.])
2491 return self.xy_to_coord(x, y, cs=cs)
2493 def points_on_source(self, cs='xyz', **kwargs):
2495 points = points_on_rect_source(
2496 self.strike, self.dip, self.length, self.width,
2497 self.anchor, **kwargs)
2499 points[:, 0] += self.north_shift
2500 points[:, 1] += self.east_shift
2501 points[:, 2] += self.depth
2502 if cs == 'xyz':
2503 return points
2504 elif cs == 'xy':
2505 return points[:, :2]
2506 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2507 latlon = ne_to_latlon(
2508 self.lat, self.lon, points[:, 0], points[:, 1])
2510 latlon = num.array(latlon).T
2511 if cs == 'latlon':
2512 return latlon
2513 elif cs == 'lonlat':
2514 return latlon[:, ::-1]
2515 else:
2516 return num.concatenate(
2517 (latlon, points[:, 2].reshape((len(points), 1))),
2518 axis=1)
2520 def geometry(self, *args, **kwargs):
2521 from pyrocko.model import Geometry
2523 ds = self.discretize_basesource(*args, **kwargs)
2524 nx, ny = ds.nl, ds.nw
2526 def patch_outlines_xy(nx, ny):
2527 points = num.zeros((nx * ny, 2))
2528 points[:, 0] = num.tile(num.linspace(-1., 1., nx), ny)
2529 points[:, 1] = num.repeat(num.linspace(-1., 1., ny), nx)
2531 return points
2533 points_ds = patch_outlines_xy(nx + 1, ny + 1)
2534 npoints = (nx + 1) * (ny + 1)
2536 vertices = num.hstack((
2537 num.ones((npoints, 1)) * self.lat,
2538 num.ones((npoints, 1)) * self.lon,
2539 self.xy_to_coord(points_ds[:, 0], points_ds[:, 1], cs='xyz')))
2541 faces = num.array([[
2542 iy * (nx + 1) + ix,
2543 iy * (nx + 1) + ix + 1,
2544 (iy + 1) * (nx + 1) + ix + 1,
2545 (iy + 1) * (nx + 1) + ix,
2546 iy * (nx + 1) + ix]
2547 for iy in range(ny) for ix in range(nx)])
2549 xyz = self.outline('xyz')
2550 latlon = num.ones((5, 2)) * num.array([self.lat, self.lon])
2551 patchverts = num.hstack((latlon, xyz))
2553 geom = Geometry()
2554 geom.setup(vertices, faces)
2555 geom.set_outlines([patchverts])
2557 if self.stf:
2558 geom.times = num.unique(ds.times)
2560 if self.nucleation_x is not None and self.nucleation_y is not None:
2561 geom.add_property('t_arrival', ds.times)
2563 geom.add_property(
2564 'moment', ds.moments().reshape(ds.nl*ds.nw, -1))
2566 geom.add_property(
2567 'slip', num.ones_like(ds.times) * self.slip)
2569 return geom
2571 def get_nucleation_abs_coord(self, cs='xy'):
2573 if self.nucleation_x is None:
2574 return None, None
2576 coords = from_plane_coords(self.strike, self.dip, self.length,
2577 self.width, self.depth, self.nucleation_x,
2578 self.nucleation_y, lat=self.lat,
2579 lon=self.lon, north_shift=self.north_shift,
2580 east_shift=self.east_shift, cs=cs)
2581 return coords
2583 def pyrocko_moment_tensor(self, store=None, target=None):
2584 return pmt.MomentTensor(
2585 strike=self.strike,
2586 dip=self.dip,
2587 rake=self.rake,
2588 scalar_moment=self.get_moment(store, target))
2590 def pyrocko_event(self, store=None, target=None, **kwargs):
2591 return SourceWithDerivedMagnitude.pyrocko_event(
2592 self, store, target,
2593 **kwargs)
2595 @classmethod
2596 def from_pyrocko_event(cls, ev, **kwargs):
2597 d = {}
2598 mt = ev.moment_tensor
2599 if mt:
2600 (strike, dip, rake), _ = mt.both_strike_dip_rake()
2601 d.update(
2602 strike=float(strike),
2603 dip=float(dip),
2604 rake=float(rake),
2605 magnitude=float(mt.moment_magnitude()))
2607 d.update(kwargs)
2608 return super(RectangularSource, cls).from_pyrocko_event(ev, **d)
2611class PseudoDynamicRupture(SourceWithDerivedMagnitude):
2612 '''
2613 Combined Eikonal and Okada quasi-dynamic rupture model.
2615 Details are described in :doc:`/topics/pseudo-dynamic-rupture`.
2616 Note: attribute `stf` is not used so far, but kept for future applications.
2617 '''
2619 discretized_source_class = meta.DiscretizedMTSource
2621 strike = Float.T(
2622 default=0.0,
2623 help='Strike direction in [deg], measured clockwise from north.')
2625 dip = Float.T(
2626 default=0.0,
2627 help='Dip angle in [deg], measured downward from horizontal.')
2629 length = Float.T(
2630 default=10. * km,
2631 help='Length of rectangular source area in [m].')
2633 width = Float.T(
2634 default=5. * km,
2635 help='Width of rectangular source area in [m].')
2637 anchor = StringChoice.T(
2638 choices=['top', 'top_left', 'top_right', 'center', 'bottom',
2639 'bottom_left', 'bottom_right'],
2640 default='center',
2641 optional=True,
2642 help='Anchor point for positioning the plane, can be: ``top, center, '
2643 'bottom, top_left, top_right, bottom_left, '
2644 'bottom_right, center_left, center_right``.')
2646 nucleation_x__ = Array.T(
2647 default=num.array([0.]),
2648 dtype=num.float64,
2649 serialize_as='list',
2650 help='Horizontal position of rupture nucleation in normalized fault '
2651 'plane coordinates (``-1.`` = left edge, ``+1.`` = right edge).')
2653 nucleation_y__ = Array.T(
2654 default=num.array([0.]),
2655 dtype=num.float64,
2656 serialize_as='list',
2657 help='Down-dip position of rupture nucleation in normalized fault '
2658 'plane coordinates (``-1.`` = upper edge, ``+1.`` = lower edge).')
2660 nucleation_time__ = Array.T(
2661 optional=True,
2662 help='Time in [s] after origin, when nucleation points defined by '
2663 '``nucleation_x`` and ``nucleation_y`` rupture.',
2664 dtype=num.float64,
2665 serialize_as='list')
2667 gamma = Float.T(
2668 default=0.8,
2669 help='Scaling factor between rupture velocity and S-wave velocity: '
2670 r':math:`v_r = \gamma * v_s`.')
2672 nx = Int.T(
2673 default=2,
2674 help='Number of discrete source patches in x direction (along '
2675 'strike).')
2677 ny = Int.T(
2678 default=2,
2679 help='Number of discrete source patches in y direction (down dip).')
2681 slip = Float.T(
2682 optional=True,
2683 help='Maximum slip of the rectangular source [m]. '
2684 'Setting the slip the tractions/stress field '
2685 'will be normalized to accomodate the desired maximum slip.')
2687 rake = Float.T(
2688 optional=True,
2689 help='Rake angle in [deg], '
2690 'measured counter-clockwise from right-horizontal '
2691 'in on-plane view. Rake is translated into homogenous tractions '
2692 'in strike and up-dip direction. ``rake`` is mutually exclusive '
2693 'with tractions parameter.')
2695 patches = List.T(
2696 OkadaSource.T(),
2697 optional=True,
2698 help='List of all boundary elements/sub faults/fault patches.')
2700 patch_mask__ = Array.T(
2701 dtype=bool,
2702 serialize_as='list',
2703 shape=(None,),
2704 optional=True,
2705 help='Mask for all boundary elements/sub faults/fault patches. True '
2706 'leaves the patch in the calculation, False excludes the patch.')
2708 tractions = TractionField.T(
2709 optional=True,
2710 help='Traction field the rupture plane is exposed to. See the '
2711 ':py:mod:`pyrocko.gf.tractions` module for more details. '
2712 'If ``tractions=None`` and ``rake`` is given'
2713 ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will'
2714 ' be used.')
2716 coef_mat = Array.T(
2717 optional=True,
2718 help='Coefficient matrix linking traction and dislocation field.',
2719 dtype=num.float64,
2720 shape=(None, None))
2722 eikonal_decimation = Int.T(
2723 optional=True,
2724 default=1,
2725 help='Sub-source eikonal factor, a smaller eikonal factor will'
2726 ' increase the accuracy of rupture front calculation but'
2727 ' increases also the computation time.')
2729 decimation_factor = Int.T(
2730 optional=True,
2731 default=1,
2732 help='Sub-source decimation factor, a larger decimation will'
2733 ' make the result inaccurate but shorten the necessary'
2734 ' computation time (use for testing puposes only).')
2736 nthreads = Int.T(
2737 optional=True,
2738 default=1,
2739 help='Number of threads for Okada forward modelling, '
2740 'matrix inversion and calculation of point subsources. '
2741 'Note: for small/medium matrices 1 thread is most efficient.')
2743 pure_shear = Bool.T(
2744 optional=True,
2745 default=False,
2746 help='Calculate only shear tractions and omit tensile tractions.')
2748 smooth_rupture = Bool.T(
2749 default=True,
2750 help='Smooth the tractions by weighting partially ruptured'
2751 ' fault patches.')
2753 aggressive_oversampling = Bool.T(
2754 default=False,
2755 help='Aggressive oversampling for basesource discretization. '
2756 "When using 'multilinear' interpolation oversampling has"
2757 ' practically no effect.')
2759 def __init__(self, **kwargs):
2760 SourceWithDerivedMagnitude.__init__(self, **kwargs)
2761 self._interpolators = {}
2762 self.check_conflicts()
2764 @property
2765 def nucleation_x(self):
2766 return self.nucleation_x__
2768 @nucleation_x.setter
2769 def nucleation_x(self, nucleation_x):
2770 if isinstance(nucleation_x, list):
2771 nucleation_x = num.array(nucleation_x)
2773 elif not isinstance(
2774 nucleation_x, num.ndarray) and nucleation_x is not None:
2776 nucleation_x = num.array([nucleation_x])
2777 self.nucleation_x__ = nucleation_x
2779 @property
2780 def nucleation_y(self):
2781 return self.nucleation_y__
2783 @nucleation_y.setter
2784 def nucleation_y(self, nucleation_y):
2785 if isinstance(nucleation_y, list):
2786 nucleation_y = num.array(nucleation_y)
2788 elif not isinstance(nucleation_y, num.ndarray) \
2789 and nucleation_y is not None:
2790 nucleation_y = num.array([nucleation_y])
2792 self.nucleation_y__ = nucleation_y
2794 @property
2795 def nucleation(self):
2796 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
2798 if (nucl_x is None) or (nucl_y is None):
2799 return None
2801 assert nucl_x.shape[0] == nucl_y.shape[0]
2803 return num.concatenate(
2804 (nucl_x[:, num.newaxis], nucl_y[:, num.newaxis]), axis=1)
2806 @nucleation.setter
2807 def nucleation(self, nucleation):
2808 if isinstance(nucleation, list):
2809 nucleation = num.array(nucleation)
2811 assert nucleation.shape[1] == 2
2813 self.nucleation_x = nucleation[:, 0]
2814 self.nucleation_y = nucleation[:, 1]
2816 @property
2817 def nucleation_time(self):
2818 return self.nucleation_time__
2820 @nucleation_time.setter
2821 def nucleation_time(self, nucleation_time):
2822 if not isinstance(nucleation_time, num.ndarray) \
2823 and nucleation_time is not None:
2824 nucleation_time = num.array([nucleation_time])
2826 self.nucleation_time__ = nucleation_time
2828 @property
2829 def patch_mask(self):
2830 if (self.patch_mask__ is not None and
2831 self.patch_mask__.shape == (self.nx * self.ny,)):
2833 return self.patch_mask__
2834 else:
2835 return num.ones(self.nx * self.ny, dtype=bool)
2837 @patch_mask.setter
2838 def patch_mask(self, patch_mask):
2839 if isinstance(patch_mask, list):
2840 patch_mask = num.array(patch_mask)
2842 self.patch_mask__ = patch_mask
2844 def get_tractions(self):
2845 '''
2846 Get source traction vectors.
2848 If :py:attr:`rake` is given, unit length directed traction vectors
2849 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are returned,
2850 else the given :py:attr:`tractions` are used.
2852 :returns:
2853 Traction vectors per patch.
2854 :rtype:
2855 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2856 '''
2858 if self.rake is not None:
2859 if num.isnan(self.rake):
2860 raise ValueError('Rake must be a real number, not NaN.')
2862 logger.warning(
2863 'Tractions are derived based on the given source rake.')
2864 tractions = DirectedTractions(rake=self.rake)
2865 else:
2866 tractions = self.tractions
2867 return tractions.get_tractions(self.nx, self.ny, self.patches)
2869 def get_scaled_tractions(self, store):
2870 '''
2871 Get traction vectors rescaled to given slip.
2873 Opposing to :py:meth:`get_tractions` traction vectors
2874 (:py:class:`~pyrocko.gf.tractions.DirectedTractions`) are rescaled to
2875 the given :py:attr:`slip` before returning. If no :py:attr:`slip` and
2876 :py:attr:`rake` are provided, the given :py:attr:`tractions` are
2877 returned without scaling.
2879 :param store:
2880 Green's function database (needs to cover whole region of the
2881 source).
2882 :type store:
2883 :py:class:`~pyrocko.gf.store.Store`
2885 :returns:
2886 Traction vectors per patch.
2887 :rtype:
2888 :py:class:`~numpy.ndarray`: ``(n_patches, 3)``.
2889 '''
2890 tractions = self.tractions
2891 factor = 1.
2893 if self.rake is not None and self.slip is not None:
2894 if num.isnan(self.rake):
2895 raise ValueError('Rake must be a real number, not NaN.')
2897 self.discretize_patches(store)
2898 slip_0t = max(num.linalg.norm(
2899 self.get_slip(scale_slip=False),
2900 axis=1))
2902 factor = self.slip / slip_0t
2903 tractions = DirectedTractions(rake=self.rake)
2905 return tractions.get_tractions(self.nx, self.ny, self.patches) * factor
2907 def base_key(self):
2908 return SourceWithDerivedMagnitude.base_key(self) + (
2909 self.slip,
2910 self.strike,
2911 self.dip,
2912 self.rake,
2913 self.length,
2914 self.width,
2915 float(self.nucleation_x.mean()),
2916 float(self.nucleation_y.mean()),
2917 self.decimation_factor,
2918 self.anchor,
2919 self.pure_shear,
2920 self.gamma,
2921 tuple(self.patch_mask))
2923 def check_conflicts(self):
2924 if self.tractions and self.rake:
2925 raise AttributeError(
2926 'Tractions and rake are mutually exclusive.')
2927 if self.tractions is None and self.rake is None:
2928 self.rake = 0.
2930 def get_magnitude(self, store=None, target=None):
2931 '''
2932 Get total seismic moment magnitude Mw.
2934 :param store:
2935 GF store to guide the discretization and providing the earthmodel
2936 which is needed to calculate moment from slip.
2937 :type store:
2938 :py:class:`~pyrocko.gf.store.Store`
2940 :param target:
2941 Target, used to get GF interpolation settings.
2942 :type target:
2943 :py:class:`pyrocko.gf.targets.Target`
2945 :returns:
2946 Moment magnitude
2947 :rtype:
2948 float
2949 '''
2950 self.check_conflicts()
2951 if self.slip is not None or self.tractions is not None:
2952 if store is None:
2953 raise DerivedMagnitudeError(
2954 'Magnitude for a rectangular source with slip or '
2955 'tractions defined can only be derived when earth model '
2956 'is set.')
2958 moment_rate, calc_times = self.discretize_basesource(
2959 store, target=target).get_moment_rate(store.config.deltat)
2961 deltat = num.concatenate((
2962 (num.diff(calc_times)[0],),
2963 num.diff(calc_times)))
2965 return float(pmt.moment_to_magnitude(
2966 num.sum(moment_rate * deltat)))
2968 else:
2969 return float(pmt.moment_to_magnitude(1.0))
2971 def get_factor(self):
2972 return 1.0
2974 def outline(self, cs='xyz'):
2975 '''
2976 Get source outline corner coordinates.
2978 :param cs:
2979 :ref:`Output coordinate system <coordinate-system-names>`.
2980 :type cs:
2981 str
2983 :returns:
2984 Corner points in desired coordinate system.
2985 :rtype:
2986 :py:class:`~numpy.ndarray`: ``(5, [2, 3])``.
2987 '''
2988 points = outline_rect_source(self.strike, self.dip, self.length,
2989 self.width, self.anchor)
2991 points[:, 0] += self.north_shift
2992 points[:, 1] += self.east_shift
2993 points[:, 2] += self.depth
2994 if cs == 'xyz':
2995 return points
2996 elif cs == 'xy':
2997 return points[:, :2]
2998 elif cs in ('latlon', 'lonlat', 'latlondepth'):
2999 latlon = ne_to_latlon(
3000 self.lat, self.lon, points[:, 0], points[:, 1])
3002 latlon = num.array(latlon).T
3003 if cs == 'latlon':
3004 return latlon
3005 elif cs == 'lonlat':
3006 return latlon[:, ::-1]
3007 else:
3008 return num.concatenate(
3009 (latlon, points[:, 2].reshape((len(points), 1))),
3010 axis=1)
3012 def points_on_source(self, cs='xyz', **kwargs):
3013 '''
3014 Convert relative plane coordinates to geographical coordinates.
3016 Given x and y coordinates (relative source coordinates between -1.
3017 and 1.) are converted to desired geographical coordinates. Coordinates
3018 need to be given as :py:class:`~numpy.ndarray` arguments ``points_x``
3019 and ``points_y``.
3021 :param cs:
3022 :ref:`Output coordinate system <coordinate-system-names>`.
3023 :type cs:
3024 str
3026 :returns:
3027 Point coordinates in desired coordinate system.
3028 :rtype:
3029 :py:class:`~numpy.ndarray`: ``(n_points, [2, 3])``.
3030 '''
3031 points = points_on_rect_source(
3032 self.strike, self.dip, self.length, self.width,
3033 self.anchor, **kwargs)
3035 points[:, 0] += self.north_shift
3036 points[:, 1] += self.east_shift
3037 points[:, 2] += self.depth
3038 if cs == 'xyz':
3039 return points
3040 elif cs == 'xy':
3041 return points[:, :2]
3042 elif cs in ('latlon', 'lonlat', 'latlondepth'):
3043 latlon = ne_to_latlon(
3044 self.lat, self.lon, points[:, 0], points[:, 1])
3046 latlon = num.array(latlon).T
3047 if cs == 'latlon':
3048 return latlon
3049 elif cs == 'lonlat':
3050 return latlon[:, ::-1]
3051 else:
3052 return num.concatenate(
3053 (latlon, points[:, 2].reshape((len(points), 1))),
3054 axis=1)
3056 def pyrocko_moment_tensor(self, store=None, target=None):
3057 '''
3058 Get overall moment tensor of the rupture.
3060 :param store:
3061 GF store to guide the discretization and providing the earthmodel
3062 which is needed to calculate moment from slip.
3063 :type store:
3064 :py:class:`~pyrocko.gf.store.Store`
3066 :param target:
3067 Target, used to get GF interpolation settings.
3068 :type target:
3069 :py:class:`pyrocko.gf.targets.Target`
3071 :returns:
3072 Moment tensor.
3073 :rtype:
3074 :py:class:`~pyrocko.moment_tensor.MomentTensor`
3075 '''
3076 if store is not None:
3077 if not self.patches:
3078 self.discretize_patches(store)
3080 data = self.get_slip()
3081 else:
3082 data = self.get_tractions()
3084 weights = num.linalg.norm(data, axis=1)
3085 weights /= weights.sum()
3087 rakes = num.arctan2(data[:, 1], data[:, 0]) * r2d
3088 rake = num.average(rakes, weights=weights)
3090 return pmt.MomentTensor(
3091 strike=self.strike,
3092 dip=self.dip,
3093 rake=rake,
3094 scalar_moment=self.get_moment(store, target))
3096 def pyrocko_event(self, store=None, target=None, **kwargs):
3097 return SourceWithDerivedMagnitude.pyrocko_event(
3098 self, store, target,
3099 **kwargs)
3101 @classmethod
3102 def from_pyrocko_event(cls, ev, **kwargs):
3103 d = {}
3104 mt = ev.moment_tensor
3105 if mt:
3106 (strike, dip, rake), _ = mt.both_strike_dip_rake()
3107 d.update(
3108 strike=float(strike),
3109 dip=float(dip),
3110 rake=float(rake))
3112 d.update(kwargs)
3113 return super(PseudoDynamicRupture, cls).from_pyrocko_event(ev, **d)
3115 def _discretize_points(self, store, *args, **kwargs):
3116 '''
3117 Discretize source plane with equal vertical and horizontal spacing.
3119 Additional ``*args`` and ``**kwargs`` are passed to
3120 :py:meth:`points_on_source`.
3122 :param store:
3123 Green's function database (needs to cover whole region of the
3124 source).
3125 :type store:
3126 :py:class:`~pyrocko.gf.store.Store`
3128 :returns:
3129 Number of points in strike and dip direction, distance
3130 between adjacent points, coordinates (latlondepth) and coordinates
3131 (xy on fault) for discrete points.
3132 :rtype:
3133 (int, int, float, :py:class:`~numpy.ndarray`,
3134 :py:class:`~numpy.ndarray`).
3135 '''
3136 anch_x, anch_y = map_anchor[self.anchor]
3138 npoints = int(self.width // km) + 1
3139 points = num.zeros((npoints, 3))
3140 points[:, 1] = num.linspace(-1., 1., npoints)
3141 points[:, 1] = (points[:, 1] - anch_y) * self.width/2
3143 rotmat = pmt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.0)
3144 points = num.dot(rotmat.T, points.T).T
3145 points[:, 2] += self.depth
3147 vs_min = store.config.get_vs(
3148 self.lat, self.lon, points,
3149 interpolation='nearest_neighbor')
3150 vr_min = max(vs_min.min(), .5*km) * self.gamma
3152 oversampling = 10.
3153 delta_l = self.length / (self.nx * oversampling)
3154 delta_w = self.width / (self.ny * oversampling)
3156 delta = self.eikonal_decimation * num.min([
3157 store.config.deltat * vr_min / oversampling,
3158 delta_l, delta_w] + [
3159 deltas for deltas in store.config.deltas])
3161 delta = delta_w / num.ceil(delta_w / delta)
3163 nx = int(num.ceil(self.length / delta)) + 1
3164 ny = int(num.ceil(self.width / delta)) + 1
3166 rem_l = (nx-1)*delta - self.length
3167 lim_x = rem_l / self.length
3169 points_xy = num.zeros((nx * ny, 2))
3170 points_xy[:, 0] = num.repeat(
3171 num.linspace(-1.-lim_x, 1.+lim_x, nx), ny)
3172 points_xy[:, 1] = num.tile(
3173 num.linspace(-1., 1., ny), nx)
3175 points = self.points_on_source(
3176 points_x=points_xy[:, 0],
3177 points_y=points_xy[:, 1],
3178 **kwargs)
3180 return nx, ny, delta, points, points_xy
3182 def _discretize_rupture_v(self, store, interpolation='nearest_neighbor',
3183 points=None):
3184 '''
3185 Get rupture velocity for discrete points on source plane.
3187 :param store:
3188 Green's function database (needs to cover the whole region of the
3189 source)
3190 :type store:
3191 :py:class:`~pyrocko.gf.store.Store`
3193 :param interpolation:
3194 Interpolation method to use (choose between ``'nearest_neighbor'``
3195 and ``'multilinear'``).
3196 :type interpolation:
3197 str
3199 :param points:
3200 Coordinates on fault (-1.:1.) of discrete points.
3201 :type points:
3202 :py:class:`~numpy.ndarray`: ``(n_points, 2)``
3204 :returns:
3205 Rupture velocity assumed as :math:`v_s * \\gamma` for discrete
3206 points.
3207 :rtype:
3208 :py:class:`~numpy.ndarray`: ``(n_points, )``.
3209 '''
3211 if points is None:
3212 _, _, _, points, _ = self._discretize_points(store, cs='xyz')
3214 return store.config.get_vs(
3215 self.lat, self.lon,
3216 points=points,
3217 interpolation=interpolation) * self.gamma
3219 def discretize_time(
3220 self, store, interpolation='nearest_neighbor',
3221 vr=None, times=None, *args, **kwargs):
3222 '''
3223 Get rupture start time for discrete points on source plane.
3225 :param store:
3226 Green's function database (needs to cover whole region of the
3227 source)
3228 :type store:
3229 :py:class:`~pyrocko.gf.store.Store`
3231 :param interpolation:
3232 Interpolation method to use (choose between ``'nearest_neighbor'``
3233 and ``'multilinear'``).
3234 :type interpolation:
3235 str
3237 :param vr:
3238 Array, containing rupture user defined rupture velocity values.
3239 :type vr:
3240 :py:class:`~numpy.ndarray`
3242 :param times:
3243 Array, containing zeros, where rupture is starting, real positive
3244 numbers at later secondary nucleation points and -1, where time
3245 will be calculated. If not given, rupture starts at nucleation_x,
3246 nucleation_y. Times are given for discrete points with equal
3247 horizontal and vertical spacing.
3248 :type times:
3249 :py:class:`~numpy.ndarray`
3251 :returns:
3252 Coordinates (latlondepth), coordinates (xy), rupture velocity,
3253 rupture propagation time of discrete points.
3254 :rtype:
3255 :py:class:`~numpy.ndarray`: ``(n_points, 3)``,
3256 :py:class:`~numpy.ndarray`: ``(n_points, 2)``,
3257 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``,
3258 :py:class:`~numpy.ndarray`: ``(n_points_dip, n_points_strike)``.
3259 '''
3260 nx, ny, delta, points, points_xy = self._discretize_points(
3261 store, cs='xyz')
3263 if vr is None or vr.shape != tuple((nx, ny)):
3264 if vr:
3265 logger.warning(
3266 'Given rupture velocities are not in right shape: '
3267 '(%i, %i), but needed is (%i, %i).', *vr.shape + (nx, ny))
3268 vr = self._discretize_rupture_v(store, interpolation, points)\
3269 .reshape(nx, ny)
3271 if vr.shape != tuple((nx, ny)):
3272 logger.warning(
3273 'Given rupture velocities are not in right shape. Therefore'
3274 ' standard rupture velocity array is used.')
3276 def initialize_times():
3277 nucl_x, nucl_y = self.nucleation_x, self.nucleation_y
3279 if nucl_x.shape != nucl_y.shape:
3280 raise ValueError(
3281 'Nucleation coordinates have different shape.')
3283 dist_points = num.array([
3284 num.linalg.norm(points_xy - num.array([x, y]).ravel(), axis=1)
3285 for x, y in zip(nucl_x, nucl_y)])
3286 nucl_indices = num.argmin(dist_points, axis=1)
3288 if self.nucleation_time is None:
3289 nucl_times = num.zeros_like(nucl_indices)
3290 else:
3291 if self.nucleation_time.shape == nucl_x.shape:
3292 nucl_times = self.nucleation_time
3293 else:
3294 raise ValueError(
3295 'Nucleation coordinates and times have different '
3296 'shapes')
3298 t = num.full(nx * ny, -1.)
3299 t[nucl_indices] = nucl_times
3300 return t.reshape(nx, ny)
3302 if times is None:
3303 times = initialize_times()
3304 elif times.shape != tuple((nx, ny)):
3305 times = initialize_times()
3306 logger.warning(
3307 'Given times are not in right shape. Therefore standard time'
3308 ' array is used.')
3310 eikonal_ext.eikonal_solver_fmm_cartesian(
3311 speeds=vr, times=times, delta=delta)
3313 return points, points_xy, vr, times
3315 def get_vr_time_interpolators(
3316 self, store, interpolation='nearest_neighbor', force=False,
3317 **kwargs):
3318 '''
3319 Get interpolators for rupture velocity and rupture time.
3321 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3323 :param store:
3324 Green's function database (needs to cover whole region of the
3325 source).
3326 :type store:
3327 :py:class:`~pyrocko.gf.store.Store`
3329 :param interpolation:
3330 Interpolation method to use (choose between ``'nearest_neighbor'``
3331 and ``'multilinear'``).
3332 :type interpolation:
3333 str
3335 :param force:
3336 Force recalculation of the interpolators (e.g. after change of
3337 nucleation point locations/times). Default is ``False``.
3338 :type force:
3339 bool
3340 '''
3341 interp_map = {'multilinear': 'linear', 'nearest_neighbor': 'nearest'}
3342 if interpolation not in interp_map:
3343 raise TypeError(
3344 'Interpolation method %s not available' % interpolation)
3346 if not self._interpolators.get(interpolation, False) or force:
3347 _, points_xy, vr, times = self.discretize_time(
3348 store, **kwargs)
3350 if self.length <= 0.:
3351 raise ValueError(
3352 'length must be larger then 0. not %g' % self.length)
3354 if self.width <= 0.:
3355 raise ValueError(
3356 'width must be larger then 0. not %g' % self.width)
3358 nx, ny = times.shape
3359 anch_x, anch_y = map_anchor[self.anchor]
3361 points_xy[:, 0] = (points_xy[:, 0] - anch_x) * self.length / 2.
3362 points_xy[:, 1] = (points_xy[:, 1] - anch_y) * self.width / 2.
3364 ascont = num.ascontiguousarray
3366 self._interpolators[interpolation] = (
3367 nx, ny, times, vr,
3368 RegularGridInterpolator(
3369 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3370 times,
3371 method=interp_map[interpolation]),
3372 RegularGridInterpolator(
3373 (ascont(points_xy[::ny, 0]), ascont(points_xy[:ny, 1])),
3374 vr,
3375 method=interp_map[interpolation]))
3377 return self._interpolators[interpolation]
3379 def discretize_patches(
3380 self, store, interpolation='nearest_neighbor', force=False,
3381 grid_shape=(),
3382 **kwargs):
3383 '''
3384 Get rupture start time and OkadaSource elements for points on rupture.
3386 All source elements and their corresponding center points are
3387 calculated and stored in the :py:attr:`patches` attribute.
3389 Additional ``**kwargs`` are passed to :py:meth:`discretize_time`.
3391 :param store:
3392 Green's function database (needs to cover whole region of the
3393 source).
3394 :type store:
3395 :py:class:`~pyrocko.gf.store.Store`
3397 :param interpolation:
3398 Interpolation method to use (choose between ``'nearest_neighbor'``
3399 and ``'multilinear'``).
3400 :type interpolation:
3401 str
3403 :param force:
3404 Force recalculation of the vr and time interpolators ( e.g. after
3405 change of nucleation point locations/times). Default is ``False``.
3406 :type force:
3407 bool
3409 :param grid_shape:
3410 Desired sub fault patch grid size (nlength, nwidth). Either factor
3411 or grid_shape should be set.
3412 :type grid_shape:
3413 :py:class:`tuple` of :py:class:`int`
3414 '''
3415 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3416 self.get_vr_time_interpolators(
3417 store,
3418 interpolation=interpolation, force=force, **kwargs)
3419 anch_x, anch_y = map_anchor[self.anchor]
3421 al = self.length / 2.
3422 aw = self.width / 2.
3423 al1 = -(al + anch_x * al)
3424 al2 = al - anch_x * al
3425 aw1 = -aw + anch_y * aw
3426 aw2 = aw + anch_y * aw
3427 assert num.abs([al1, al2]).sum() == self.length
3428 assert num.abs([aw1, aw2]).sum() == self.width
3430 def get_lame(*a, **kw):
3431 shear_mod = store.config.get_shear_moduli(*a, **kw)
3432 lamb = store.config.get_vp(*a, **kw)**2 \
3433 * store.config.get_rho(*a, **kw) - 2. * shear_mod
3434 return shear_mod, lamb / (2. * (lamb + shear_mod))
3436 shear_mod, poisson = get_lame(
3437 self.lat, self.lon,
3438 num.array([[self.north_shift, self.east_shift, self.depth]]),
3439 interpolation=interpolation)
3441 okada_src = OkadaSource(
3442 lat=self.lat, lon=self.lon,
3443 strike=self.strike, dip=self.dip,
3444 north_shift=self.north_shift, east_shift=self.east_shift,
3445 depth=self.depth,
3446 al1=al1, al2=al2, aw1=aw1, aw2=aw2,
3447 poisson=poisson.mean(),
3448 shearmod=shear_mod.mean(),
3449 opening=kwargs.get('opening', 0.))
3451 if not (self.nx and self.ny):
3452 if grid_shape:
3453 self.nx, self.ny = grid_shape
3454 else:
3455 self.nx = nx
3456 self.ny = ny
3458 source_disc, source_points = okada_src.discretize(self.nx, self.ny)
3460 shear_mod, poisson = get_lame(
3461 self.lat, self.lon,
3462 num.array([src.source_patch()[:3] for src in source_disc]),
3463 interpolation=interpolation)
3465 if (self.nx, self.ny) != (nx, ny):
3466 times_interp = time_interpolator(
3467 num.ascontiguousarray(source_points[:, :2]))
3468 vr_interp = vr_interpolator(
3469 num.ascontiguousarray(source_points[:, :2]))
3470 else:
3471 times_interp = times.T.ravel()
3472 vr_interp = vr.T.ravel()
3474 for isrc, src in enumerate(source_disc):
3475 src.vr = vr_interp[isrc]
3476 src.time = times_interp[isrc] + self.time
3478 self.patches = source_disc
3480 def discretize_basesource(self, store, target=None):
3481 '''
3482 Prepare source for synthetic waveform calculation.
3484 :param store:
3485 Green's function database (needs to cover whole region of the
3486 source).
3487 :type store:
3488 :py:class:`~pyrocko.gf.store.Store`
3490 :param target:
3491 Target information.
3492 :type target:
3493 :py:class:`~pyrocko.gf.targets.Target`
3495 :returns:
3496 Source discretized by a set of moment tensors and times.
3497 :rtype:
3498 :py:class:`~pyrocko.gf.meta.DiscretizedMTSource`
3499 '''
3500 if not target:
3501 interpolation = 'nearest_neighbor'
3502 else:
3503 interpolation = target.interpolation
3505 if not self.patches:
3506 self.discretize_patches(store, interpolation)
3508 if self.coef_mat is None:
3509 self.calc_coef_mat()
3511 delta_slip, slip_times = self.get_delta_slip(store)
3512 npatches = self.nx * self.ny
3513 ntimes = slip_times.size
3515 anch_x, anch_y = map_anchor[self.anchor]
3517 pln = self.length / self.nx
3518 pwd = self.width / self.ny
3520 patch_coords = num.array([
3521 (p.ix, p.iy)
3522 for p in self.patches]).reshape(self.nx, self.ny, 2)
3524 # boundary condition is zero-slip
3525 # is not valid to avoid unwished interpolation effects
3526 slip_grid = num.zeros((self.nx + 2, self.ny + 2, ntimes, 3))
3527 slip_grid[1:-1, 1:-1, :, :] = \
3528 delta_slip.reshape(self.nx, self.ny, ntimes, 3)
3530 slip_grid[0, 0, :, :] = slip_grid[1, 1, :, :]
3531 slip_grid[0, -1, :, :] = slip_grid[1, -2, :, :]
3532 slip_grid[-1, 0, :, :] = slip_grid[-2, 1, :, :]
3533 slip_grid[-1, -1, :, :] = slip_grid[-2, -2, :, :]
3535 slip_grid[1:-1, 0, :, :] = slip_grid[1:-1, 1, :, :]
3536 slip_grid[1:-1, -1, :, :] = slip_grid[1:-1, -2, :, :]
3537 slip_grid[0, 1:-1, :, :] = slip_grid[1, 1:-1, :, :]
3538 slip_grid[-1, 1:-1, :, :] = slip_grid[-2, 1:-1, :, :]
3540 def make_grid(patch_parameter):
3541 grid = num.zeros((self.nx + 2, self.ny + 2))
3542 grid[1:-1, 1:-1] = patch_parameter.reshape(self.nx, self.ny)
3544 grid[0, 0] = grid[1, 1]
3545 grid[0, -1] = grid[1, -2]
3546 grid[-1, 0] = grid[-2, 1]
3547 grid[-1, -1] = grid[-2, -2]
3549 grid[1:-1, 0] = grid[1:-1, 1]
3550 grid[1:-1, -1] = grid[1:-1, -2]
3551 grid[0, 1:-1] = grid[1, 1:-1]
3552 grid[-1, 1:-1] = grid[-2, 1:-1]
3554 return grid
3556 lamb = self.get_patch_attribute('lamb')
3557 mu = self.get_patch_attribute('shearmod')
3559 lamb_grid = make_grid(lamb)
3560 mu_grid = make_grid(mu)
3562 coords_x = num.zeros(self.nx + 2)
3563 coords_x[1:-1] = patch_coords[:, 0, 0]
3564 coords_x[0] = coords_x[1] - pln / 2
3565 coords_x[-1] = coords_x[-2] + pln / 2
3567 coords_y = num.zeros(self.ny + 2)
3568 coords_y[1:-1] = patch_coords[0, :, 1]
3569 coords_y[0] = coords_y[1] - pwd / 2
3570 coords_y[-1] = coords_y[-2] + pwd / 2
3572 slip_interp = RegularGridInterpolator(
3573 (coords_x, coords_y, slip_times),
3574 slip_grid, method='nearest')
3576 lamb_interp = RegularGridInterpolator(
3577 (coords_x, coords_y),
3578 lamb_grid, method='nearest')
3580 mu_interp = RegularGridInterpolator(
3581 (coords_x, coords_y),
3582 mu_grid, method='nearest')
3584 # discretize basesources
3585 mindeltagf = min(tuple(
3586 (self.length / self.nx, self.width / self.ny) +
3587 tuple(store.config.deltas)))
3589 nl = int((1. / self.decimation_factor) *
3590 num.ceil(pln / mindeltagf)) + 1
3591 nw = int((1. / self.decimation_factor) *
3592 num.ceil(pwd / mindeltagf)) + 1
3593 nsrc_patch = int(nl * nw)
3594 dl = pln / nl
3595 dw = pwd / nw
3597 patch_area = dl * dw
3599 xl = num.linspace(-0.5 * (pln - dl), 0.5 * (pln - dl), nl)
3600 xw = num.linspace(-0.5 * (pwd - dw), 0.5 * (pwd - dw), nw)
3602 base_coords = num.zeros((nsrc_patch, 3))
3603 base_coords[:, 0] = num.tile(xl, nw)
3604 base_coords[:, 1] = num.repeat(xw, nl)
3605 base_coords = num.tile(base_coords, (npatches, 1))
3607 center_coords = num.zeros((npatches, 3))
3608 center_coords[:, 0] = num.repeat(
3609 num.arange(self.nx) * pln + pln / 2, self.ny) - self.length / 2
3610 center_coords[:, 1] = num.tile(
3611 num.arange(self.ny) * pwd + pwd / 2, self.nx) - self.width / 2
3613 base_coords -= center_coords.repeat(nsrc_patch, axis=0)
3614 nbaselocs = base_coords.shape[0]
3616 base_interp = base_coords.repeat(ntimes, axis=0)
3618 base_times = num.tile(slip_times, nbaselocs)
3619 base_interp[:, 0] -= anch_x * self.length / 2
3620 base_interp[:, 1] -= anch_y * self.width / 2
3621 base_interp[:, 2] = base_times
3623 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3624 store, interpolation=interpolation)
3626 time_eikonal_max = time_interpolator.values.max()
3628 nbasesrcs = base_interp.shape[0]
3629 delta_slip = slip_interp(base_interp).reshape(nbaselocs, ntimes, 3)
3630 lamb = lamb_interp(base_interp[:, :2]).ravel()
3631 mu = mu_interp(base_interp[:, :2]).ravel()
3633 if False:
3634 try:
3635 import matplotlib.pyplot as plt
3636 coords = base_coords.copy()
3637 norm = num.sum(num.linalg.norm(delta_slip, axis=2), axis=1)
3638 plt.scatter(coords[:, 0], coords[:, 1], c=norm)
3639 plt.show()
3640 except AttributeError:
3641 pass
3643 base_interp[:, 2] = 0.
3644 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
3645 base_interp = num.dot(rotmat.T, base_interp.T).T
3646 base_interp[:, 0] += self.north_shift
3647 base_interp[:, 1] += self.east_shift
3648 base_interp[:, 2] += self.depth
3650 slip_strike = delta_slip[:, :, 0].ravel()
3651 slip_dip = delta_slip[:, :, 1].ravel()
3652 slip_norm = delta_slip[:, :, 2].ravel()
3654 slip_shear = num.linalg.norm([slip_strike, slip_dip], axis=0)
3655 slip_rake = r2d * num.arctan2(slip_dip, slip_strike)
3657 m6s = okada_ext.patch2m6(
3658 strikes=num.full(nbasesrcs, self.strike, dtype=float),
3659 dips=num.full(nbasesrcs, self.dip, dtype=float),
3660 rakes=slip_rake,
3661 disl_shear=slip_shear,
3662 disl_norm=slip_norm,
3663 lamb=lamb,
3664 mu=mu,
3665 nthreads=self.nthreads)
3667 m6s *= patch_area
3669 dl = -self.patches[0].al1 + self.patches[0].al2
3670 dw = -self.patches[0].aw1 + self.patches[0].aw2
3672 base_times[base_times > time_eikonal_max] = time_eikonal_max
3674 ds = meta.DiscretizedMTSource(
3675 lat=self.lat,
3676 lon=self.lon,
3677 times=base_times + self.time,
3678 north_shifts=base_interp[:, 0],
3679 east_shifts=base_interp[:, 1],
3680 depths=base_interp[:, 2],
3681 m6s=m6s,
3682 dl=dl,
3683 dw=dw,
3684 nl=self.nx,
3685 nw=self.ny)
3687 return ds
3689 def calc_coef_mat(self):
3690 '''
3691 Calculate coefficients connecting tractions and dislocations.
3692 '''
3693 if not self.patches:
3694 raise ValueError(
3695 'Patches are needed. Please calculate them first.')
3697 self.coef_mat = make_okada_coefficient_matrix(
3698 self.patches, nthreads=self.nthreads, pure_shear=self.pure_shear)
3700 def get_patch_attribute(self, attr):
3701 '''
3702 Get patch attributes.
3704 :param attr:
3705 Name of selected attribute (see
3706 :py:class`pyrocko.modelling.okada.OkadaSource`).
3707 :type attr:
3708 str
3710 :returns:
3711 Array with attribute value for each fault patch.
3712 :rtype:
3713 :py:class:`~numpy.ndarray`
3715 '''
3716 if not self.patches:
3717 raise ValueError(
3718 'Patches are needed. Please calculate them first.')
3719 return num.array([getattr(p, attr) for p in self.patches])
3721 def get_slip(
3722 self,
3723 time=None,
3724 scale_slip=True,
3725 interpolation='nearest_neighbor',
3726 **kwargs):
3727 '''
3728 Get slip per subfault patch for given time after rupture start.
3730 :param time:
3731 Time after origin [s], for which slip is computed. If not
3732 given, final static slip is returned.
3733 :type time:
3734 float
3736 :param scale_slip:
3737 If ``True`` and :py:attr:`slip` given, all slip values are scaled
3738 to fit the given maximum slip.
3739 :type scale_slip:
3740 bool
3742 :param interpolation:
3743 Interpolation method to use (choose between ``'nearest_neighbor'``
3744 and ``'multilinear'``).
3745 :type interpolation:
3746 str
3748 :returns:
3749 Inverted dislocations (:math:`u_{strike}, u_{dip}, u_{tensile}`)
3750 for each source patch.
3751 :rtype:
3752 :py:class:`~numpy.ndarray`: ``(n_sources, 3)``
3753 '''
3755 if self.patches is None:
3756 raise ValueError(
3757 'Please discretize the source first (discretize_patches())')
3758 npatches = len(self.patches)
3759 tractions = self.get_tractions()
3760 time_patch_max = self.get_patch_attribute('time').max() - self.time
3762 time_patch = time
3763 if time is None:
3764 time_patch = time_patch_max
3766 if self.coef_mat is None:
3767 self.calc_coef_mat()
3769 if tractions.shape != (npatches, 3):
3770 raise AttributeError(
3771 'The traction vector is of invalid shape.'
3772 ' Required shape is (npatches, 3)')
3774 patch_mask = num.ones(npatches, dtype=bool)
3775 if self.patch_mask is not None:
3776 patch_mask = self.patch_mask
3778 times = self.get_patch_attribute('time') - self.time
3779 times[~patch_mask] = time_patch + 1. # exlcude unmasked patches
3780 relevant_sources = num.nonzero(times <= time_patch)[0]
3781 disloc_est = num.zeros_like(tractions)
3783 if self.smooth_rupture:
3784 patch_activation = num.zeros(npatches)
3786 nx, ny, times, vr, time_interpolator, vr_interpolator = \
3787 self.get_vr_time_interpolators(
3788 store, interpolation=interpolation)
3790 # Getting the native Eikonal grid, bit hackish
3791 points_x = num.round(time_interpolator.grid[0], decimals=2)
3792 points_y = num.round(time_interpolator.grid[1], decimals=2)
3793 times_eikonal = time_interpolator.values
3795 time_max = time
3796 if time is None:
3797 time_max = times_eikonal.max()
3799 for ip, p in enumerate(self.patches):
3800 ul = num.round((p.ix + p.al1, p.iy + p.aw1), decimals=2)
3801 lr = num.round((p.ix + p.al2, p.iy + p.aw2), decimals=2)
3803 idx_length = num.logical_and(
3804 points_x >= ul[0], points_x <= lr[0])
3805 idx_width = num.logical_and(
3806 points_y >= ul[1], points_y <= lr[1])
3808 times_patch = times_eikonal[num.ix_(idx_length, idx_width)]
3809 if times_patch.size == 0:
3810 raise AttributeError('could not use smooth_rupture')
3812 patch_activation[ip] = \
3813 (times_patch <= time_max).sum() / times_patch.size
3815 if time_patch == 0 and time_patch != time_patch_max:
3816 patch_activation[ip] = 0.
3818 patch_activation[~patch_mask] = 0. # exlcude unmasked patches
3820 relevant_sources = num.nonzero(patch_activation > 0.)[0]
3822 if relevant_sources.size == 0:
3823 return disloc_est
3825 indices_disl = num.repeat(relevant_sources * 3, 3)
3826 indices_disl[1::3] += 1
3827 indices_disl[2::3] += 2
3829 disloc_est[relevant_sources] = invert_fault_dislocations_bem(
3830 stress_field=tractions[relevant_sources, :].ravel(),
3831 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3832 pure_shear=self.pure_shear, nthreads=self.nthreads,
3833 epsilon=None,
3834 **kwargs)
3836 if self.smooth_rupture:
3837 disloc_est *= patch_activation[:, num.newaxis]
3839 if scale_slip and self.slip is not None:
3840 disloc_tmax = num.zeros(npatches)
3842 indices_disl = num.repeat(num.nonzero(patch_mask)[0] * 3, 3)
3843 indices_disl[1::3] += 1
3844 indices_disl[2::3] += 2
3846 disloc_tmax[patch_mask] = num.linalg.norm(
3847 invert_fault_dislocations_bem(
3848 stress_field=tractions[patch_mask, :].ravel(),
3849 coef_mat=self.coef_mat[indices_disl, :][:, indices_disl],
3850 pure_shear=self.pure_shear, nthreads=self.nthreads,
3851 epsilon=None,
3852 **kwargs), axis=1)
3854 disloc_tmax_max = disloc_tmax.max()
3855 if disloc_tmax_max == 0.:
3856 logger.warning(
3857 'slip scaling not performed. Maximum slip is 0.')
3859 disloc_est *= self.slip / disloc_tmax_max
3861 return disloc_est
3863 def get_delta_slip(
3864 self,
3865 store=None,
3866 deltat=None,
3867 delta=True,
3868 interpolation='nearest_neighbor',
3869 **kwargs):
3870 '''
3871 Get slip change snapshots.
3873 The time interval, within which the slip changes are computed is
3874 determined by the sampling rate of the Green's function database or
3875 ``deltat``. Additional ``**kwargs`` are passed to :py:meth:`get_slip`.
3877 :param store:
3878 Green's function database (needs to cover whole region of of the
3879 source). Its sampling interval is used as time increment for slip
3880 difference calculation. Either ``deltat`` or ``store`` should be
3881 given.
3882 :type store:
3883 :py:class:`~pyrocko.gf.store.Store`
3885 :param deltat:
3886 Time interval for slip difference calculation [s]. Either
3887 ``deltat`` or ``store`` should be given.
3888 :type deltat:
3889 float
3891 :param delta:
3892 If ``True``, slip differences between two time steps are given. If
3893 ``False``, cumulative slip for all time steps.
3894 :type delta:
3895 bool
3897 :param interpolation:
3898 Interpolation method to use (choose between ``'nearest_neighbor'``
3899 and ``'multilinear'``).
3900 :type interpolation:
3901 str
3903 :returns:
3904 Displacement changes(:math:`\\Delta u_{strike},
3905 \\Delta u_{dip} , \\Delta u_{tensile}`) for each source patch and
3906 time; corner times, for which delta slip is computed. The order of
3907 displacement changes array is:
3909 .. math::
3911 &[[\\\\
3912 &[\\Delta u_{strike, patch1, t1},
3913 \\Delta u_{dip, patch1, t1},
3914 \\Delta u_{tensile, patch1, t1}],\\\\
3915 &[\\Delta u_{strike, patch1, t2},
3916 \\Delta u_{dip, patch1, t2},
3917 \\Delta u_{tensile, patch1, t2}]\\\\
3918 &], [\\\\
3919 &[\\Delta u_{strike, patch2, t1}, ...],\\\\
3920 &[\\Delta u_{strike, patch2, t2}, ...]]]\\\\
3922 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
3923 :py:class:`~numpy.ndarray`: ``(n_times, )``
3924 '''
3925 if store and deltat:
3926 raise AttributeError(
3927 'Argument collision. '
3928 'Please define only the store or the deltat argument.')
3930 if store:
3931 deltat = store.config.deltat
3933 if not deltat:
3934 raise AttributeError('Please give a GF store or set deltat.')
3936 npatches = len(self.patches)
3938 _, _, _, _, time_interpolator, _ = self.get_vr_time_interpolators(
3939 store, interpolation=interpolation)
3940 tmax = time_interpolator.values.max()
3942 calc_times = num.arange(0., tmax + deltat, deltat)
3943 calc_times[calc_times > tmax] = tmax
3945 disloc_est = num.zeros((npatches, calc_times.size, 3))
3947 for itime, t in enumerate(calc_times):
3948 disloc_est[:, itime, :] = self.get_slip(
3949 time=t, scale_slip=False, **kwargs)
3951 if self.slip:
3952 disloc_tmax = num.linalg.norm(
3953 self.get_slip(scale_slip=False, time=tmax),
3954 axis=1)
3956 disloc_tmax_max = disloc_tmax.max()
3957 if disloc_tmax_max == 0.:
3958 logger.warning(
3959 'Slip scaling not performed. Maximum slip is 0.')
3960 else:
3961 disloc_est *= self.slip / disloc_tmax_max
3963 if not delta:
3964 return disloc_est, calc_times
3966 # if we have only one timestep there is no gradient
3967 if calc_times.size > 1:
3968 disloc_init = disloc_est[:, 0, :]
3969 disloc_est = num.diff(disloc_est, axis=1)
3970 disloc_est = num.concatenate((
3971 disloc_init[:, num.newaxis, :], disloc_est), axis=1)
3973 calc_times = calc_times
3975 return disloc_est, calc_times
3977 def get_slip_rate(self, *args, **kwargs):
3978 '''
3979 Get slip rate inverted from patches.
3981 The time interval, within which the slip rates are computed is
3982 determined by the sampling rate of the Green's function database or
3983 ``deltat``. Additional ``*args`` and ``**kwargs`` are passed to
3984 :py:meth:`get_delta_slip`.
3986 :returns:
3987 Slip rates (:math:`\\Delta u_{strike}/\\Delta t`,
3988 :math:`\\Delta u_{dip}/\\Delta t, \\Delta u_{tensile}/\\Delta t`)
3989 for each source patch and time; corner times, for which slip rate
3990 is computed. The order of sliprate array is:
3992 .. math::
3994 &[[\\\\
3995 &[\\Delta u_{strike, patch1, t1}/\\Delta t,
3996 \\Delta u_{dip, patch1, t1}/\\Delta t,
3997 \\Delta u_{tensile, patch1, t1}/\\Delta t],\\\\
3998 &[\\Delta u_{strike, patch1, t2}/\\Delta t,
3999 \\Delta u_{dip, patch1, t2}/\\Delta t,
4000 \\Delta u_{tensile, patch1, t2}/\\Delta t]], [\\\\
4001 &[\\Delta u_{strike, patch2, t1}/\\Delta t, ...],\\\\
4002 &[\\Delta u_{strike, patch2, t2}/\\Delta t, ...]]]\\\\
4004 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times, 3)``,
4005 :py:class:`~numpy.ndarray`: ``(n_times, )``
4006 '''
4007 ddisloc_est, calc_times = self.get_delta_slip(
4008 *args, delta=True, **kwargs)
4010 dt = num.concatenate(
4011 [(num.diff(calc_times)[0], ), num.diff(calc_times)])
4012 slip_rate = num.linalg.norm(ddisloc_est, axis=2) / dt
4014 return slip_rate, calc_times
4016 def get_moment_rate_patches(self, *args, **kwargs):
4017 '''
4018 Get scalar seismic moment rate for each patch individually.
4020 Additional ``*args`` and ``**kwargs`` are passed to
4021 :py:meth:`get_slip_rate`.
4023 :returns:
4024 Seismic moment rate for each source patch and time; corner times,
4025 for which patch moment rate is computed based on slip rate. The
4026 order of the moment rate array is:
4028 .. math::
4030 &[\\\\
4031 &[(\\Delta M / \\Delta t)_{patch1, t1},
4032 (\\Delta M / \\Delta t)_{patch1, t2}, ...],\\\\
4033 &[(\\Delta M / \\Delta t)_{patch2, t1},
4034 (\\Delta M / \\Delta t)_{patch, t2}, ...],\\\\
4035 &[...]]\\\\
4037 :rtype: :py:class:`~numpy.ndarray`: ``(n_sources, n_times)``,
4038 :py:class:`~numpy.ndarray`: ``(n_times, )``
4039 '''
4040 slip_rate, calc_times = self.get_slip_rate(*args, **kwargs)
4042 shear_mod = self.get_patch_attribute('shearmod')
4043 p_length = self.get_patch_attribute('length')
4044 p_width = self.get_patch_attribute('width')
4046 dA = p_length * p_width
4048 mom_rate = shear_mod[:, num.newaxis] * slip_rate * dA[:, num.newaxis]
4050 return mom_rate, calc_times
4052 def get_moment_rate(self, store, target=None, deltat=None):
4053 '''
4054 Get seismic source moment rate for the total source (STF).
4056 :param store:
4057 Green's function database (needs to cover whole region of of the
4058 source). Its ``deltat`` [s] is used as time increment for slip
4059 difference calculation. Either ``deltat`` or ``store`` should be
4060 given.
4061 :type store:
4062 :py:class:`~pyrocko.gf.store.Store`
4064 :param target:
4065 Target information, needed for interpolation method.
4066 :type target:
4067 :py:class:`~pyrocko.gf.targets.Target`
4069 :param deltat:
4070 Time increment for slip difference calculation [s]. If not given
4071 ``store.deltat`` is used.
4072 :type deltat:
4073 float
4075 :return:
4076 Seismic moment rate [Nm/s] for each time; corner times, for which
4077 moment rate is computed. The order of the moment rate array is:
4079 .. math::
4081 &[\\\\
4082 &(\\Delta M / \\Delta t)_{t1},\\\\
4083 &(\\Delta M / \\Delta t)_{t2},\\\\
4084 &...]\\\\
4086 :rtype:
4087 :py:class:`~numpy.ndarray`: ``(n_times, )``,
4088 :py:class:`~numpy.ndarray`: ``(n_times, )``
4089 '''
4090 if not deltat:
4091 deltat = store.config.deltat
4092 return self.discretize_basesource(
4093 store, target=target).get_moment_rate(deltat)
4095 def get_moment(self, *args, **kwargs):
4096 '''
4097 Get cumulative seismic moment.
4099 Additional ``*args`` and ``**kwargs`` are passed to
4100 :py:meth:`get_magnitude`.
4102 :returns:
4103 Cumulative seismic moment in [Nm].
4104 :rtype:
4105 float
4106 '''
4107 return float(pmt.magnitude_to_moment(self.get_magnitude(
4108 *args, **kwargs)))
4110 def rescale_slip(self, magnitude=None, moment=None, **kwargs):
4111 '''
4112 Rescale source slip based on given target magnitude or seismic moment.
4114 Rescale the maximum source slip to fit the source moment magnitude or
4115 seismic moment to the given target values. Either ``magnitude`` or
4116 ``moment`` need to be given. Additional ``**kwargs`` are passed to
4117 :py:meth:`get_moment`.
4119 :param magnitude:
4120 Target moment magnitude :math:`M_\\mathrm{w}` as in
4121 [Hanks and Kanamori, 1979]
4122 :type magnitude:
4123 float
4125 :param moment:
4126 Target seismic moment :math:`M_0` [Nm].
4127 :type moment:
4128 float
4129 '''
4130 if self.slip is None:
4131 self.slip = 1.
4132 logger.warning('No slip found for rescaling. '
4133 'An initial slip of 1 m is assumed.')
4135 if magnitude is None and moment is None:
4136 raise ValueError(
4137 'Either target magnitude or moment need to be given.')
4139 moment_init = self.get_moment(**kwargs)
4141 if magnitude is not None:
4142 moment = pmt.magnitude_to_moment(magnitude)
4144 self.slip *= moment / moment_init
4146 def get_centroid(self, store, *args, **kwargs):
4147 '''
4148 Centroid of the pseudo dynamic rupture model.
4150 The centroid location and time are derived from the locations and times
4151 of the individual patches weighted with their moment contribution.
4152 Additional ``**kwargs`` are passed to :py:meth:`pyrocko_moment_tensor`.
4154 :param store:
4155 Green's function database (needs to cover whole region of of the
4156 source). Its ``deltat`` [s] is used as time increment for slip
4157 difference calculation. Either ``deltat`` or ``store`` should be
4158 given.
4159 :type store:
4160 :py:class:`~pyrocko.gf.store.Store`
4162 :returns:
4163 The centroid location and associated moment tensor.
4164 :rtype:
4165 :py:class:`pyrocko.model.event.Event`
4166 '''
4167 _, _, _, _, time, _ = self.get_vr_time_interpolators(store)
4168 t_max = time.values.max()
4170 moment_rate, times = self.get_moment_rate_patches(deltat=t_max)
4172 moment = num.sum(moment_rate * times, axis=1)
4173 weights = moment / moment.sum()
4175 norths = self.get_patch_attribute('north_shift')
4176 easts = self.get_patch_attribute('east_shift')
4177 depths = self.get_patch_attribute('depth')
4179 centroid_n = num.sum(weights * norths)
4180 centroid_e = num.sum(weights * easts)
4181 centroid_d = num.sum(weights * depths)
4183 centroid_lat, centroid_lon = ne_to_latlon(
4184 self.lat, self.lon, centroid_n, centroid_e)
4186 moment_rate_, times = self.get_moment_rate(store)
4187 delta_times = num.concatenate((
4188 [times[1] - times[0]],
4189 num.diff(times)))
4190 moment_src = delta_times * moment_rate
4192 centroid_t = num.sum(
4193 moment_src / num.sum(moment_src) * times) + self.time
4195 mt = self.pyrocko_moment_tensor(store, *args, **kwargs)
4197 return model.Event(
4198 lat=centroid_lat,
4199 lon=centroid_lon,
4200 depth=centroid_d,
4201 time=centroid_t,
4202 moment_tensor=mt,
4203 magnitude=mt.magnitude,
4204 duration=t_max)
4206 def get_coulomb_failure_stress(
4207 self,
4208 receiver_points,
4209 friction,
4210 pressure,
4211 strike,
4212 dip,
4213 rake,
4214 time=None,
4215 *args,
4216 **kwargs):
4217 '''
4218 Calculate Coulomb failure stress change CFS.
4220 The function obtains the Coulomb failure stress change :math:`\\Delta
4221 \\sigma_C` at arbitrary receiver points with a commonly oriented
4222 receiver plane assuming:
4224 .. math::
4226 \\Delta \\sigma_C = \\sigma_S - \\mu (\\sigma_N - \\Delta p)
4228 with the shear stress :math:`\\sigma_S`, the coefficient of friction
4229 :math:`\\mu`, the normal stress :math:`\\sigma_N`, and the pore fluid
4230 pressure change :math:`\\Delta p`. Each receiver point is characterized
4231 by its geographical coordinates, and depth. The required receiver plane
4232 orientation is defined by ``strike``, ``dip``, and ``rake``. The
4233 Coulomb failure stress change is calculated for a given time after
4234 rupture origin time.
4236 :param receiver_points:
4237 Location of the receiver points in Northing, Easting, and depth in
4238 [m].
4239 :type receiver_points:
4240 :py:class:`~numpy.ndarray`: ``(n_receiver, 3)``
4242 :param friction:
4243 Coefficient of friction.
4244 :type friction:
4245 float
4247 :param pressure:
4248 Pore pressure change in [Pa].
4249 :type pressure:
4250 float
4252 :param strike:
4253 Strike of the receiver plane in [deg].
4254 :type strike:
4255 float
4257 :param dip:
4258 Dip of the receiver plane in [deg].
4259 :type dip:
4260 float
4262 :param rake:
4263 Rake of the receiver plane in [deg].
4264 :type rake:
4265 float
4267 :param time:
4268 Time after origin [s], for which the resulting :math:`\\Delta
4269 \\Sigma_c` is computed. If not given, :math:`\\Delta \\Sigma_c` is
4270 derived based on the final static slip.
4271 :type time:
4272 float
4274 :returns:
4275 The Coulomb failure stress change :math:`\\Delta \\Sigma_c` at each
4276 receiver point in [Pa].
4277 :rtype:
4278 :py:class:`~numpy.ndarray`: ``(n_receiver,)``
4279 '''
4280 # dislocation at given time
4281 source_slip = self.get_slip(time=time, scale_slip=True)
4283 # source planes
4284 source_patches = num.array([
4285 src.source_patch() for src in self.patches])
4287 # earth model
4288 lambda_mean = num.mean([src.lamb for src in self.patches])
4289 mu_mean = num.mean([src.shearmod for src in self.patches])
4291 # Dislocation and spatial derivatives from okada in NED
4292 results = okada_ext.okada(
4293 source_patches,
4294 source_slip,
4295 receiver_points,
4296 lambda_mean,
4297 mu_mean,
4298 rotate_sdn=False, # TODO Check
4299 stack_sources=0, # TODO Check
4300 *args, **kwargs)
4302 # resolve stress tensor (sum!)
4303 diag_ind = [0, 4, 8]
4304 kron = num.zeros(9)
4305 kron[diag_ind] = 1.
4306 kron = kron[num.newaxis, num.newaxis, :]
4308 eps = 0.5 * (
4309 results[:, :, 3:] +
4310 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
4312 dilatation \
4313 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
4315 stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps
4317 # superposed stress of all sources at receiver locations
4318 stress_sum = num.sum(stress, axis=0)
4320 # get shear and normal stress from stress tensor
4321 strike_rad = d2r * strike
4322 dip_rad = d2r * dip
4323 rake_rad = d2r * rake
4325 n_rec = receiver_points.shape[0]
4326 stress_normal = num.zeros(n_rec)
4327 tau = num.zeros(n_rec)
4329 # Get vectors in receiver fault normal (ns), strike (rst) and
4330 # dip (rdi) direction
4331 ns = num.zeros(3)
4332 rst = num.zeros(3)
4333 rdi = num.zeros(3)
4335 ns[0] = num.sin(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4336 ns[1] = num.sin(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4337 ns[2] = -num.cos(dip_rad)
4339 rst[0] = num.cos(strike_rad)
4340 rst[1] = num.sin(strike_rad)
4341 rst[2] = 0.0
4343 rdi[0] = num.cos(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
4344 rdi[1] = num.cos(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
4345 rdi[2] = num.sin(dip_rad)
4347 ts = rst * num.cos(rake_rad) - rdi * num.sin(rake_rad)
4349 stress_normal = num.sum(
4350 num.tile(ns, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4352 tau = num.sum(
4353 num.tile(ts, 3) * stress_sum * num.repeat(ns, 3), axis=1)
4355 # calculate cfs using formula above and return
4356 return tau + friction * (stress_normal + pressure)
4359class DoubleDCSource(SourceWithMagnitude):
4360 '''
4361 Two double-couple point sources separated in space and time.
4362 Moment share between the sub-sources is controlled by the
4363 parameter mix.
4364 The position of the subsources is dependent on the moment
4365 distribution between the two sources. Depth, east and north
4366 shift are given for the centroid between the two double-couples.
4367 The subsources will positioned according to their moment shares
4368 around this centroid position.
4369 This is done according to their delta parameters, which are
4370 therefore in relation to that centroid.
4371 Note that depth of the subsources therefore can be
4372 depth+/-delta_depth. For shallow earthquakes therefore
4373 the depth has to be chosen deeper to avoid sampling
4374 above surface.
4375 '''
4377 strike1 = Float.T(
4378 default=0.0,
4379 help='strike direction in [deg], measured clockwise from north')
4381 dip1 = Float.T(
4382 default=90.0,
4383 help='dip angle in [deg], measured downward from horizontal')
4385 azimuth = Float.T(
4386 default=0.0,
4387 help='azimuth to second double-couple [deg], '
4388 'measured at first, clockwise from north')
4390 rake1 = Float.T(
4391 default=0.0,
4392 help='rake angle in [deg], '
4393 'measured counter-clockwise from right-horizontal '
4394 'in on-plane view')
4396 strike2 = Float.T(
4397 default=0.0,
4398 help='strike direction in [deg], measured clockwise from north')
4400 dip2 = Float.T(
4401 default=90.0,
4402 help='dip angle in [deg], measured downward from horizontal')
4404 rake2 = Float.T(
4405 default=0.0,
4406 help='rake angle in [deg], '
4407 'measured counter-clockwise from right-horizontal '
4408 'in on-plane view')
4410 delta_time = Float.T(
4411 default=0.0,
4412 help='separation of double-couples in time (t2-t1) [s]')
4414 delta_depth = Float.T(
4415 default=0.0,
4416 help='difference in depth (z2-z1) [m]')
4418 distance = Float.T(
4419 default=0.0,
4420 help='distance between the two double-couples [m]')
4422 mix = Float.T(
4423 default=0.5,
4424 help='how to distribute the moment to the two doublecouples '
4425 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1')
4427 stf1 = STF.T(
4428 optional=True,
4429 help='Source time function of subsource 1 '
4430 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4432 stf2 = STF.T(
4433 optional=True,
4434 help='Source time function of subsource 2 '
4435 '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
4437 discretized_source_class = meta.DiscretizedMTSource
4439 def base_key(self):
4440 return (
4441 self.time, self.depth, self.lat, self.north_shift,
4442 self.lon, self.east_shift, type(self).__name__) + \
4443 self.effective_stf1_pre().base_key() + \
4444 self.effective_stf2_pre().base_key() + (
4445 self.strike1, self.dip1, self.rake1,
4446 self.strike2, self.dip2, self.rake2,
4447 self.delta_time, self.delta_depth,
4448 self.azimuth, self.distance, self.mix)
4450 def get_factor(self):
4451 return self.moment
4453 def effective_stf1_pre(self):
4454 return self.stf1 or self.stf or g_unit_pulse
4456 def effective_stf2_pre(self):
4457 return self.stf2 or self.stf or g_unit_pulse
4459 def effective_stf_post(self):
4460 return g_unit_pulse
4462 def split(self):
4463 a1 = 1.0 - self.mix
4464 a2 = self.mix
4465 delta_north = math.cos(self.azimuth * d2r) * self.distance
4466 delta_east = math.sin(self.azimuth * d2r) * self.distance
4468 dc1 = DCSource(
4469 lat=self.lat,
4470 lon=self.lon,
4471 time=self.time - self.delta_time * a2,
4472 north_shift=self.north_shift - delta_north * a2,
4473 east_shift=self.east_shift - delta_east * a2,
4474 depth=self.depth - self.delta_depth * a2,
4475 moment=self.moment * a1,
4476 strike=self.strike1,
4477 dip=self.dip1,
4478 rake=self.rake1,
4479 stf=self.stf1 or self.stf)
4481 dc2 = DCSource(
4482 lat=self.lat,
4483 lon=self.lon,
4484 time=self.time + self.delta_time * a1,
4485 north_shift=self.north_shift + delta_north * a1,
4486 east_shift=self.east_shift + delta_east * a1,
4487 depth=self.depth + self.delta_depth * a1,
4488 moment=self.moment * a2,
4489 strike=self.strike2,
4490 dip=self.dip2,
4491 rake=self.rake2,
4492 stf=self.stf2 or self.stf)
4494 return [dc1, dc2]
4496 def discretize_basesource(self, store, target=None):
4497 a1 = 1.0 - self.mix
4498 a2 = self.mix
4499 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4500 rake=self.rake1, scalar_moment=a1)
4501 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4502 rake=self.rake2, scalar_moment=a2)
4504 delta_north = math.cos(self.azimuth * d2r) * self.distance
4505 delta_east = math.sin(self.azimuth * d2r) * self.distance
4507 times1, amplitudes1 = self.effective_stf1_pre().discretize_t(
4508 store.config.deltat, self.time - self.delta_time * a2)
4510 times2, amplitudes2 = self.effective_stf2_pre().discretize_t(
4511 store.config.deltat, self.time + self.delta_time * a1)
4513 nt1 = times1.size
4514 nt2 = times2.size
4516 ds = meta.DiscretizedMTSource(
4517 lat=self.lat,
4518 lon=self.lon,
4519 times=num.concatenate((times1, times2)),
4520 north_shifts=num.concatenate((
4521 num.repeat(self.north_shift - delta_north * a2, nt1),
4522 num.repeat(self.north_shift + delta_north * a1, nt2))),
4523 east_shifts=num.concatenate((
4524 num.repeat(self.east_shift - delta_east * a2, nt1),
4525 num.repeat(self.east_shift + delta_east * a1, nt2))),
4526 depths=num.concatenate((
4527 num.repeat(self.depth - self.delta_depth * a2, nt1),
4528 num.repeat(self.depth + self.delta_depth * a1, nt2))),
4529 m6s=num.vstack((
4530 mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis],
4531 mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis])))
4533 return ds
4535 def pyrocko_moment_tensor(self, store=None, target=None):
4536 a1 = 1.0 - self.mix
4537 a2 = self.mix
4538 mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1,
4539 rake=self.rake1,
4540 scalar_moment=a1 * self.moment)
4541 mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2,
4542 rake=self.rake2,
4543 scalar_moment=a2 * self.moment)
4544 return pmt.MomentTensor(m=mot1.m() + mot2.m())
4546 def pyrocko_event(self, store=None, target=None, **kwargs):
4547 return SourceWithMagnitude.pyrocko_event(
4548 self, store, target,
4549 moment_tensor=self.pyrocko_moment_tensor(store, target),
4550 **kwargs)
4552 @classmethod
4553 def from_pyrocko_event(cls, ev, **kwargs):
4554 d = {}
4555 mt = ev.moment_tensor
4556 if mt:
4557 (strike, dip, rake), _ = mt.both_strike_dip_rake()
4558 d.update(
4559 strike1=float(strike),
4560 dip1=float(dip),
4561 rake1=float(rake),
4562 strike2=float(strike),
4563 dip2=float(dip),
4564 rake2=float(rake),
4565 mix=0.0,
4566 magnitude=float(mt.moment_magnitude()))
4568 d.update(kwargs)
4569 source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d)
4570 source.stf1 = source.stf
4571 source.stf2 = HalfSinusoidSTF(effective_duration=0.)
4572 source.stf = None
4573 return source
4576class RingfaultSource(SourceWithMagnitude):
4577 '''
4578 A ring fault with vertical doublecouples.
4579 '''
4581 diameter = Float.T(
4582 default=1.0,
4583 help='diameter of the ring in [m]')
4585 sign = Float.T(
4586 default=1.0,
4587 help='inside of the ring moves up (+1) or down (-1)')
4589 strike = Float.T(
4590 default=0.0,
4591 help='strike direction of the ring plane, clockwise from north,'
4592 ' in [deg]')
4594 dip = Float.T(
4595 default=0.0,
4596 help='dip angle of the ring plane from horizontal in [deg]')
4598 npointsources = Int.T(
4599 default=360,
4600 help='number of point sources to use')
4602 discretized_source_class = meta.DiscretizedMTSource
4604 def base_key(self):
4605 return Source.base_key(self) + (
4606 self.strike, self.dip, self.diameter, self.npointsources)
4608 def get_factor(self):
4609 return self.sign * self.moment
4611 def discretize_basesource(self, store=None, target=None):
4612 n = self.npointsources
4613 phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False)
4615 points = num.zeros((n, 3))
4616 points[:, 0] = num.cos(phi) * 0.5 * self.diameter
4617 points[:, 1] = num.sin(phi) * 0.5 * self.diameter
4619 rotmat = pmt.euler_to_matrix(self.dip * d2r, self.strike * d2r, 0.0)
4620 points = num.dot(rotmat.T, points.T).T # !!! ?
4622 points[:, 0] += self.north_shift
4623 points[:, 1] += self.east_shift
4624 points[:, 2] += self.depth
4626 m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90.,
4627 scalar_moment=1.0 / n).m())
4629 rotmats = num.transpose(
4630 [[num.cos(phi), num.sin(phi), num.zeros(n)],
4631 [-num.sin(phi), num.cos(phi), num.zeros(n)],
4632 [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1))
4634 ms = num.zeros((n, 3, 3))
4635 for i in range(n):
4636 mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i]))
4637 ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat))
4639 m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2],
4640 ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T
4642 times, amplitudes = self.effective_stf_pre().discretize_t(
4643 store.config.deltat, self.time)
4645 nt = times.size
4647 return meta.DiscretizedMTSource(
4648 times=num.tile(times, n),
4649 lat=self.lat,
4650 lon=self.lon,
4651 north_shifts=num.repeat(points[:, 0], nt),
4652 east_shifts=num.repeat(points[:, 1], nt),
4653 depths=num.repeat(points[:, 2], nt),
4654 m6s=num.repeat(m6s, nt, axis=0) * num.tile(
4655 amplitudes, n)[:, num.newaxis])
4658class CombiSource(Source):
4659 '''
4660 Composite source model.
4661 '''
4663 discretized_source_class = meta.DiscretizedMTSource
4665 subsources = List.T(Source.T())
4667 def __init__(self, subsources=[], **kwargs):
4668 if not subsources:
4669 raise BadRequest(
4670 'Need at least one sub-source to create a CombiSource object.')
4672 lats = num.array(
4673 [subsource.lat for subsource in subsources], dtype=float)
4674 lons = num.array(
4675 [subsource.lon for subsource in subsources], dtype=float)
4677 lat, lon = lats[0], lons[0]
4678 if not num.all(lats == lat) and num.all(lons == lon):
4679 subsources = [s.clone() for s in subsources]
4680 for subsource in subsources[1:]:
4681 subsource.set_origin(lat, lon)
4683 depth = float(num.mean([p.depth for p in subsources]))
4684 time = float(num.mean([p.time for p in subsources]))
4685 north_shift = float(num.mean([p.north_shift for p in subsources]))
4686 east_shift = float(num.mean([p.east_shift for p in subsources]))
4687 kwargs.update(
4688 time=time,
4689 lat=float(lat),
4690 lon=float(lon),
4691 north_shift=north_shift,
4692 east_shift=east_shift,
4693 depth=depth)
4695 Source.__init__(self, subsources=subsources, **kwargs)
4697 def get_factor(self):
4698 return 1.0
4700 def discretize_basesource(self, store, target=None):
4701 dsources = []
4702 for sf in self.subsources:
4703 ds = sf.discretize_basesource(store, target)
4704 ds.m6s *= sf.get_factor()
4705 dsources.append(ds)
4707 return meta.DiscretizedMTSource.combine(dsources)
4710class CombiSFSource(Source):
4711 '''
4712 Composite source model.
4713 '''
4715 discretized_source_class = meta.DiscretizedSFSource
4717 subsources = List.T(Source.T())
4719 def __init__(self, subsources=[], **kwargs):
4720 if not subsources:
4721 raise BadRequest(
4722 'Need at least one sub-source to create a CombiSFSource '
4723 'object.')
4725 lats = num.array(
4726 [subsource.lat for subsource in subsources], dtype=float)
4727 lons = num.array(
4728 [subsource.lon for subsource in subsources], dtype=float)
4730 lat, lon = lats[0], lons[0]
4731 if not num.all(lats == lat) and num.all(lons == lon):
4732 subsources = [s.clone() for s in subsources]
4733 for subsource in subsources[1:]:
4734 subsource.set_origin(lat, lon)
4736 depth = float(num.mean([p.depth for p in subsources]))
4737 time = float(num.mean([p.time for p in subsources]))
4738 north_shift = float(num.mean([p.north_shift for p in subsources]))
4739 east_shift = float(num.mean([p.east_shift for p in subsources]))
4740 kwargs.update(
4741 time=time,
4742 lat=float(lat),
4743 lon=float(lon),
4744 north_shift=north_shift,
4745 east_shift=east_shift,
4746 depth=depth)
4748 Source.__init__(self, subsources=subsources, **kwargs)
4750 def get_factor(self):
4751 return 1.0
4753 def discretize_basesource(self, store, target=None):
4754 dsources = []
4755 for sf in self.subsources:
4756 ds = sf.discretize_basesource(store, target)
4757 ds.forces *= sf.get_factor()
4758 dsources.append(ds)
4760 return meta.DiscretizedSFSource.combine(dsources)
4763class SFSource(Source):
4764 '''
4765 A single force point source.
4767 Supported GF schemes: `'elastic5'`.
4768 '''
4770 discretized_source_class = meta.DiscretizedSFSource
4772 fn = Float.T(
4773 default=0.,
4774 help='northward component of single force [N]')
4776 fe = Float.T(
4777 default=0.,
4778 help='eastward component of single force [N]')
4780 fd = Float.T(
4781 default=0.,
4782 help='downward component of single force [N]')
4784 def __init__(self, **kwargs):
4785 Source.__init__(self, **kwargs)
4787 def base_key(self):
4788 return Source.base_key(self) + (self.fn, self.fe, self.fd)
4790 def get_factor(self):
4791 return 1.0
4793 @property
4794 def force(self):
4795 return math.sqrt(self.fn**2 + self.fe**2 + self.fd**2)
4797 def discretize_basesource(self, store, target=None):
4798 times, amplitudes = self.effective_stf_pre().discretize_t(
4799 store.config.deltat, self.time)
4800 forces = amplitudes[:, num.newaxis] * num.array(
4801 [[self.fn, self.fe, self.fd]], dtype=float)
4803 return meta.DiscretizedSFSource(forces=forces,
4804 **self._dparams_base_repeated(times))
4806 def pyrocko_event(self, store=None, target=None, **kwargs):
4807 return Source.pyrocko_event(
4808 self, store, target,
4809 **kwargs)
4811 @classmethod
4812 def from_pyrocko_event(cls, ev, **kwargs):
4813 d = {}
4814 d.update(kwargs)
4815 return super(SFSource, cls).from_pyrocko_event(ev, **d)
4818class SimpleLandslideSource(Source):
4819 '''
4820 A single force landslide source respecting conservation of momentum.
4822 The landslide is modelled point-like in space but with individual source
4823 time functions for each force component. The source time functions
4824 :py:class:`SimpleLandslideSTF` impose the constraint that everything is at
4825 rest before and after the event but are allowed to have different
4826 acceleration and deceleration durations. It should thus be suitable as a
4827 first order approximation of a rock fall or landslide.
4828 For realistic landslides, the horizontal accelerations and decelerations
4829 can be delayed with respect to the vertical ones but cannot precede them.
4831 Supported GF schemes: `'elastic5'`.
4832 '''
4834 discretized_source_class = meta.DiscretizedSFSource
4836 stf_mode = STFMode.T(
4837 default='pre',
4838 help='SimpleLandslideSource only works with `stf_mode == "pre"`.')
4840 impulse_n = Float.T(
4841 default=0.,
4842 help='northward component of impulse [Ns]')
4844 impulse_e = Float.T(
4845 default=0.,
4846 help='eastward component of impulse [Ns]')
4848 impulse_d = Float.T(
4849 default=0.,
4850 help='downward component of impulse [Ns]')
4852 azimuth = Float.T(
4853 default=0.,
4854 help='azimuth direction of the mass movement [deg]')
4856 stf_v = SimpleLandslideSTF.T(
4857 default=SimpleLandslideSTF.D(),
4858 help='source time function for vertical force component')
4860 stf_h = SimpleLandslideSTF.T(
4861 default=SimpleLandslideSTF.D(),
4862 help='source time function for horizontal force component')
4864 anchor_stf = StringChoice.T(
4865 choices=['onset', 'centroid'],
4866 default='onset',
4867 help='``"onset"``: STFs start at origin time ``"centroid"``: STFs all '
4868 'start at the same time but so that the centroid is at the given '
4869 'origin time.')
4871 def __init__(self, **kwargs):
4872 Source.__init__(self, **kwargs)
4874 def base_key(self):
4875 return Source.base_key(self) + (
4876 self.impulse_n, self.impulse_e, self.impulse_d) \
4877 + self.stf_v.base_key() + self.stf_h.base_key()
4879 def get_factor(self):
4880 return 1.0
4882 def discretize_basesource(self, store, target=None):
4883 if self.stf_mode != 'pre':
4884 raise Exception(
4885 'SimpleLandslideSource: '
4886 'Only works with `stf_mode == "pre"`.')
4888 if self.stf is not None:
4889 raise Exception(
4890 'SimpleLandslideSource: '
4891 'Setting `stf` is not supported: use `stf_v` and `stf_h`.')
4893 if self.anchor_stf == 'centroid':
4894 duration_acc = num.array([
4895 self.stf_h.duration_acceleration,
4896 self.stf_h.duration_acceleration,
4897 self.stf_v.duration_acceleration], dtype=float)
4899 impulse = num.array([
4900 self.impulse_n,
4901 self.impulse_e,
4902 self.impulse_d], dtype=float)
4904 tshift_centroid = \
4905 - num.sum(duration_acc * impulse**2) \
4906 / num.sum(impulse**2)
4908 elif self.anchor_stf == 'onset':
4909 tshift_centroid = 0.0
4911 times, amplitudes = self.stf_v.discretize_t(
4912 store.config.deltat,
4913 self.time + tshift_centroid)
4915 forces_v = num.zeros((times.size, 3))
4916 forces_v[:, 2] = amplitudes * self.impulse_d
4918 dsource_v = meta.DiscretizedSFSource(
4919 forces=forces_v,
4920 **self._dparams_base_repeated(times))
4922 times, amplitudes = self.stf_h.discretize_t(
4923 store.config.deltat,
4924 self.time + tshift_centroid)
4926 forces_h = num.zeros((times.size, 3))
4927 forces_h[:, 0] = \
4928 amplitudes * self.impulse_n
4929 forces_h[:, 1] = \
4930 amplitudes * self.impulse_e
4932 dsource_h = meta.DiscretizedSFSource(
4933 forces=forces_h,
4934 **self._dparams_base_repeated(times))
4936 return meta.DiscretizedSFSource.combine([dsource_v, dsource_h])
4938 def pyrocko_event(self, store=None, target=None, **kwargs):
4939 return Source.pyrocko_event(
4940 self, store, target,
4941 **kwargs)
4943 @classmethod
4944 def from_pyrocko_event(cls, ev, **kwargs):
4945 d = {}
4946 d.update(kwargs)
4947 return super(SimpleLandslideSource, cls).from_pyrocko_event(ev, **d)
4950class PorePressurePointSource(Source):
4951 '''
4952 Excess pore pressure point source.
4954 For poro-elastic initial value problem where an excess pore pressure is
4955 brought into a small source volume.
4956 '''
4958 discretized_source_class = meta.DiscretizedPorePressureSource
4960 pp = Float.T(
4961 default=1.0,
4962 help='initial excess pore pressure in [Pa]')
4964 def base_key(self):
4965 return Source.base_key(self)
4967 def get_factor(self):
4968 return self.pp
4970 def discretize_basesource(self, store, target=None):
4971 return meta.DiscretizedPorePressureSource(pp=arr(1.0),
4972 **self._dparams_base())
4975class PorePressureLineSource(Source):
4976 '''
4977 Excess pore pressure line source.
4979 The line source is centered at (north_shift, east_shift, depth).
4980 '''
4982 discretized_source_class = meta.DiscretizedPorePressureSource
4984 pp = Float.T(
4985 default=1.0,
4986 help='initial excess pore pressure in [Pa]')
4988 length = Float.T(
4989 default=0.0,
4990 help='length of the line source [m]')
4992 azimuth = Float.T(
4993 default=0.0,
4994 help='azimuth direction, clockwise from north [deg]')
4996 dip = Float.T(
4997 default=90.,
4998 help='dip direction, downward from horizontal [deg]')
5000 def base_key(self):
5001 return Source.base_key(self) + (self.azimuth, self.dip, self.length)
5003 def get_factor(self):
5004 return self.pp
5006 def discretize_basesource(self, store, target=None):
5008 n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1
5010 a = num.linspace(-0.5 * self.length, 0.5 * self.length, n)
5012 sa = math.sin(self.azimuth * d2r)
5013 ca = math.cos(self.azimuth * d2r)
5014 sd = math.sin(self.dip * d2r)
5015 cd = math.cos(self.dip * d2r)
5017 points = num.zeros((n, 3))
5018 points[:, 0] = self.north_shift + a * ca * cd
5019 points[:, 1] = self.east_shift + a * sa * cd
5020 points[:, 2] = self.depth + a * sd
5022 return meta.DiscretizedPorePressureSource(
5023 times=num.full(n, self.time),
5024 lat=self.lat,
5025 lon=self.lon,
5026 north_shifts=points[:, 0],
5027 east_shifts=points[:, 1],
5028 depths=points[:, 2],
5029 pp=num.ones(n) / n)
5032class Request(Object):
5033 '''
5034 Synthetic seismogram computation request.
5036 ::
5038 Request(**kwargs)
5039 Request(sources, targets, **kwargs)
5040 '''
5042 sources = List.T(
5043 Source.T(),
5044 help='list of sources for which to produce synthetics.')
5046 targets = List.T(
5047 Target.T(),
5048 help='list of targets for which to produce synthetics.')
5050 @classmethod
5051 def args2kwargs(cls, args):
5052 if len(args) not in (0, 2, 3):
5053 raise BadRequest('Invalid arguments.')
5055 if len(args) == 2:
5056 return dict(sources=args[0], targets=args[1])
5057 else:
5058 return {}
5060 def __init__(self, *args, **kwargs):
5061 kwargs.update(self.args2kwargs(args))
5062 sources = kwargs.pop('sources', [])
5063 targets = kwargs.pop('targets', [])
5065 if isinstance(sources, Source):
5066 sources = [sources]
5068 if isinstance(targets, Target) or isinstance(targets, StaticTarget):
5069 targets = [targets]
5071 Object.__init__(self, sources=sources, targets=targets, **kwargs)
5073 @property
5074 def targets_dynamic(self):
5075 return [t for t in self.targets if isinstance(t, Target)]
5077 @property
5078 def targets_static(self):
5079 return [t for t in self.targets if isinstance(t, StaticTarget)]
5081 @property
5082 def has_dynamic(self):
5083 return True if len(self.targets_dynamic) > 0 else False
5085 @property
5086 def has_statics(self):
5087 return True if len(self.targets_static) > 0 else False
5089 def subsources_map(self):
5090 m = defaultdict(list)
5091 for source in self.sources:
5092 m[source.base_key()].append(source)
5094 return m
5096 def subtargets_map(self):
5097 m = defaultdict(list)
5098 for target in self.targets:
5099 m[target.base_key()].append(target)
5101 return m
5103 def subrequest_map(self):
5104 ms = self.subsources_map()
5105 mt = self.subtargets_map()
5106 m = {}
5107 for (ks, ls) in ms.items():
5108 for (kt, lt) in mt.items():
5109 m[ks, kt] = (ls, lt)
5111 return m
5114class ProcessingStats(Object):
5115 t_perc_get_store_and_receiver = Float.T(default=0.)
5116 t_perc_discretize_source = Float.T(default=0.)
5117 t_perc_make_base_seismogram = Float.T(default=0.)
5118 t_perc_make_same_span = Float.T(default=0.)
5119 t_perc_post_process = Float.T(default=0.)
5120 t_perc_optimize = Float.T(default=0.)
5121 t_perc_stack = Float.T(default=0.)
5122 t_perc_static_get_store = Float.T(default=0.)
5123 t_perc_static_discretize_basesource = Float.T(default=0.)
5124 t_perc_static_sum_statics = Float.T(default=0.)
5125 t_perc_static_post_process = Float.T(default=0.)
5126 t_wallclock = Float.T(default=0.)
5127 t_cpu = Float.T(default=0.)
5128 n_read_blocks = Int.T(default=0)
5129 n_results = Int.T(default=0)
5130 n_subrequests = Int.T(default=0)
5131 n_stores = Int.T(default=0)
5132 n_records_stacked = Int.T(default=0)
5135class Response(Object):
5136 '''
5137 Resonse object to a synthetic seismogram computation request.
5138 '''
5140 request = Request.T()
5141 results_list = List.T(List.T(meta.SeismosizerResult.T()))
5142 stats = ProcessingStats.T()
5144 def pyrocko_traces(self):
5145 '''
5146 Return a list of requested
5147 :class:`~pyrocko.trace.Trace` instances.
5148 '''
5150 traces = []
5151 for results in self.results_list:
5152 for result in results:
5153 if not isinstance(result, meta.Result):
5154 continue
5155 traces.append(result.trace.pyrocko_trace())
5157 return traces
5159 def kite_scenes(self):
5160 '''
5161 Return a list of requested
5162 :class:`kite.Scene` instances.
5163 '''
5164 kite_scenes = []
5165 for results in self.results_list:
5166 for result in results:
5167 if isinstance(result, meta.KiteSceneResult):
5168 sc = result.get_scene()
5169 kite_scenes.append(sc)
5171 return kite_scenes
5173 def static_results(self):
5174 '''
5175 Return a list of requested
5176 :class:`~pyrocko.gf.meta.StaticResult` instances.
5177 '''
5178 statics = []
5179 for results in self.results_list:
5180 for result in results:
5181 if not isinstance(result, meta.StaticResult):
5182 continue
5183 statics.append(result)
5185 return statics
5187 def iter_results(self, get='pyrocko_traces'):
5188 '''
5189 Generator function to iterate over results of request.
5191 Yields associated :py:class:`Source`,
5192 :class:`~pyrocko.gf.targets.Target`,
5193 :class:`~pyrocko.trace.Trace` instances in each iteration.
5194 '''
5196 for isource, source in enumerate(self.request.sources):
5197 for itarget, target in enumerate(self.request.targets):
5198 result = self.results_list[isource][itarget]
5199 if get == 'pyrocko_traces':
5200 yield source, target, result.trace.pyrocko_trace()
5201 elif get == 'results':
5202 yield source, target, result
5204 def snuffle(self, **kwargs):
5205 '''
5206 Open *snuffler* with requested traces.
5207 '''
5209 trace.snuffle(self.pyrocko_traces(), **kwargs)
5212class Engine(Object):
5213 '''
5214 Base class for synthetic seismogram calculators.
5215 '''
5217 def get_store_ids(self):
5218 '''
5219 Get list of available GF store IDs
5220 '''
5222 return []
5225class Rule(object):
5226 pass
5229class VectorRule(Rule):
5231 def __init__(self, quantity, differentiate=0, integrate=0):
5232 self.components = [quantity + '.' + c for c in 'ned']
5233 self.differentiate = differentiate
5234 self.integrate = integrate
5236 def required_components(self, target):
5237 n, e, d = self.components
5238 sa, ca, sd, cd = target.get_sin_cos_factors()
5240 comps = []
5241 if nonzero(ca * cd):
5242 comps.append(n)
5244 if nonzero(sa * cd):
5245 comps.append(e)
5247 if nonzero(sd):
5248 comps.append(d)
5250 return tuple(comps)
5252 def apply_(self, target, base_seismogram):
5253 n, e, d = self.components
5254 sa, ca, sd, cd = target.get_sin_cos_factors()
5256 if nonzero(ca * cd):
5257 data = base_seismogram[n].data * (ca * cd)
5258 deltat = base_seismogram[n].deltat
5259 else:
5260 data = 0.0
5262 if nonzero(sa * cd):
5263 data = data + base_seismogram[e].data * (sa * cd)
5264 deltat = base_seismogram[e].deltat
5266 if nonzero(sd):
5267 data = data + base_seismogram[d].data * sd
5268 deltat = base_seismogram[d].deltat
5270 if self.differentiate:
5271 data = util.diff_fd(self.differentiate, 4, deltat, data)
5273 if self.integrate:
5274 raise NotImplementedError('Integration is not implemented yet.')
5276 return data
5279class HorizontalVectorRule(Rule):
5281 def __init__(self, quantity, differentiate=0, integrate=0):
5282 self.components = [quantity + '.' + c for c in 'ne']
5283 self.differentiate = differentiate
5284 self.integrate = integrate
5286 def required_components(self, target):
5287 n, e = self.components
5288 sa, ca, _, _ = target.get_sin_cos_factors()
5290 comps = []
5291 if nonzero(ca):
5292 comps.append(n)
5294 if nonzero(sa):
5295 comps.append(e)
5297 return tuple(comps)
5299 def apply_(self, target, base_seismogram):
5300 n, e = self.components
5301 sa, ca, _, _ = target.get_sin_cos_factors()
5303 if nonzero(ca):
5304 data = base_seismogram[n].data * ca
5305 else:
5306 data = 0.0
5308 if nonzero(sa):
5309 data = data + base_seismogram[e].data * sa
5311 if self.differentiate:
5312 deltat = base_seismogram[e].deltat
5313 data = util.diff_fd(self.differentiate, 4, deltat, data)
5315 if self.integrate:
5316 raise NotImplementedError('Integration is not implemented yet.')
5318 return data
5321class ScalarRule(Rule):
5323 def __init__(self, quantity, differentiate=0):
5324 self.c = quantity
5326 def required_components(self, target):
5327 return (self.c, )
5329 def apply_(self, target, base_seismogram):
5330 data = base_seismogram[self.c].data.copy()
5331 deltat = base_seismogram[self.c].deltat
5332 if self.differentiate:
5333 data = util.diff_fd(self.differentiate, 4, deltat, data)
5335 return data
5338class StaticDisplacement(Rule):
5340 def required_components(self, target):
5341 return tuple(['displacement.%s' % c for c in list('ned')])
5343 def apply_(self, target, base_statics):
5344 if isinstance(target, SatelliteTarget):
5345 los_fac = target.get_los_factors()
5346 base_statics['displacement.los'] =\
5347 (los_fac[:, 0] * -base_statics['displacement.d'] +
5348 los_fac[:, 1] * base_statics['displacement.e'] +
5349 los_fac[:, 2] * base_statics['displacement.n'])
5350 return base_statics
5353channel_rules = {
5354 'displacement': [VectorRule('displacement')],
5355 'rotation_displacement': [VectorRule('rotation_displacement')],
5356 'velocity': [
5357 VectorRule('velocity'),
5358 VectorRule('displacement', differentiate=1)],
5359 'acceleration': [
5360 VectorRule('acceleration'),
5361 VectorRule('velocity', differentiate=1),
5362 VectorRule('displacement', differentiate=2)],
5363 'pore_pressure': [ScalarRule('pore_pressure')],
5364 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')],
5365 'darcy_velocity': [VectorRule('darcy_velocity')],
5366}
5368static_rules = {
5369 'displacement': [StaticDisplacement()]
5370}
5373class OutOfBoundsContext(Object):
5374 source = Source.T()
5375 target = Target.T()
5376 distance = Float.T()
5377 components = List.T(String.T())
5380def process_dynamic_timeseries(work, psources, ptargets, engine, nthreads=0):
5381 dsource_cache = {}
5382 tcounters = list(range(6))
5384 store_ids = set()
5385 sources = set()
5386 targets = set()
5388 for itarget, target in enumerate(ptargets):
5389 target._id = itarget
5391 for w in work:
5392 _, _, isources, itargets = w
5394 sources.update([psources[isource] for isource in isources])
5395 targets.update([ptargets[itarget] for itarget in itargets])
5397 store_ids = set([t.store_id for t in targets])
5399 for isource, source in enumerate(psources):
5401 components = set()
5402 for itarget, target in enumerate(targets):
5403 rule = engine.get_rule(source, target)
5404 components.update(rule.required_components(target))
5406 for store_id in store_ids:
5407 store_targets = [t for t in targets if t.store_id == store_id]
5409 sample_rates = set([t.sample_rate for t in store_targets])
5410 interpolations = set([t.interpolation for t in store_targets])
5412 base_seismograms = []
5413 store_targets_out = []
5415 for samp_rate in sample_rates:
5416 for interp in interpolations:
5417 engine_targets = [
5418 t for t in store_targets if t.sample_rate == samp_rate
5419 and t.interpolation == interp]
5421 if not engine_targets:
5422 continue
5424 store_targets_out += engine_targets
5426 base_seismograms += engine.base_seismograms(
5427 source,
5428 engine_targets,
5429 components,
5430 dsource_cache,
5431 nthreads)
5433 for iseis, seismogram in enumerate(base_seismograms):
5434 for tr in seismogram.values():
5435 if tr.err != store.SeismosizerErrorEnum.SUCCESS:
5436 e = SeismosizerError(
5437 'Seismosizer failed with return code %i\n%s' % (
5438 tr.err, str(
5439 OutOfBoundsContext(
5440 source=source,
5441 target=store_targets[iseis],
5442 distance=source.distance_to(
5443 store_targets[iseis]),
5444 components=components))))
5445 raise e
5447 for seismogram, target in zip(base_seismograms, store_targets_out):
5449 try:
5450 result = engine._post_process_dynamic(
5451 seismogram, source, target)
5452 except SeismosizerError as e:
5453 result = e
5455 yield (isource, target._id, result), tcounters
5458def process_dynamic(work, psources, ptargets, engine, nthreads=0):
5459 dsource_cache = {}
5461 for w in work:
5462 _, _, isources, itargets = w
5464 sources = [psources[isource] for isource in isources]
5465 targets = [ptargets[itarget] for itarget in itargets]
5467 components = set()
5468 for target in targets:
5469 rule = engine.get_rule(sources[0], target)
5470 components.update(rule.required_components(target))
5472 for isource, source in zip(isources, sources):
5473 for itarget, target in zip(itargets, targets):
5475 try:
5476 base_seismogram, tcounters = engine.base_seismogram(
5477 source, target, components, dsource_cache, nthreads)
5478 except meta.OutOfBounds as e:
5479 e.context = OutOfBoundsContext(
5480 source=sources[0],
5481 target=targets[0],
5482 distance=sources[0].distance_to(targets[0]),
5483 components=components)
5484 raise
5486 n_records_stacked = 0
5487 t_optimize = 0.0
5488 t_stack = 0.0
5490 for _, tr in base_seismogram.items():
5491 n_records_stacked += tr.n_records_stacked
5492 t_optimize += tr.t_optimize
5493 t_stack += tr.t_stack
5495 try:
5496 result = engine._post_process_dynamic(
5497 base_seismogram, source, target)
5498 result.n_records_stacked = n_records_stacked
5499 result.n_shared_stacking = len(sources) *\
5500 len(targets)
5501 result.t_optimize = t_optimize
5502 result.t_stack = t_stack
5503 except SeismosizerError as e:
5504 result = e
5506 tcounters.append(xtime())
5507 yield (isource, itarget, result), tcounters
5510def process_static(work, psources, ptargets, engine, nthreads=0):
5511 for w in work:
5512 _, _, isources, itargets = w
5514 sources = [psources[isource] for isource in isources]
5515 targets = [ptargets[itarget] for itarget in itargets]
5517 for isource, source in zip(isources, sources):
5518 for itarget, target in zip(itargets, targets):
5519 components = engine.get_rule(source, target)\
5520 .required_components(target)
5522 try:
5523 base_statics, tcounters = engine.base_statics(
5524 source, target, components, nthreads)
5525 except meta.OutOfBounds as e:
5526 e.context = OutOfBoundsContext(
5527 source=sources[0],
5528 target=targets[0],
5529 distance=float('nan'),
5530 components=components)
5531 raise
5532 result = engine._post_process_statics(
5533 base_statics, source, target)
5534 tcounters.append(xtime())
5536 yield (isource, itarget, result), tcounters
5539class LocalEngine(Engine):
5540 '''
5541 Offline synthetic seismogram calculator.
5543 :param use_env: if ``True``, fill :py:attr:`store_superdirs` and
5544 :py:attr:`store_dirs` with paths set in environment variables
5545 GF_STORE_SUPERDIRS and GF_STORE_DIRS.
5546 :param use_config: if ``True``, fill :py:attr:`store_superdirs` and
5547 :py:attr:`store_dirs` with paths set in the user's config file.
5549 The config file can be found at :file:`~/.pyrocko/config.pf`
5551 .. code-block :: python
5553 gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/']
5554 gf_store_superdirs: ['/home/pyrocko/gf_stores/']
5555 '''
5557 store_superdirs = List.T(
5558 String.T(),
5559 help="directories which are searched for Green's function stores")
5561 store_dirs = List.T(
5562 String.T(),
5563 help="additional individual Green's function store directories")
5565 default_store_id = String.T(
5566 optional=True,
5567 help='default store ID to be used when a request does not provide '
5568 'one')
5570 def __init__(self, **kwargs):
5571 use_env = kwargs.pop('use_env', False)
5572 use_config = kwargs.pop('use_config', False)
5573 Engine.__init__(self, **kwargs)
5574 if use_env:
5575 env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '')
5576 env_store_dirs = os.environ.get('GF_STORE_DIRS', '')
5577 if env_store_superdirs:
5578 self.store_superdirs.extend(env_store_superdirs.split(':'))
5580 if env_store_dirs:
5581 self.store_dirs.extend(env_store_dirs.split(':'))
5583 if use_config:
5584 c = config.config()
5585 self.store_superdirs.extend(c.gf_store_superdirs)
5586 self.store_dirs.extend(c.gf_store_dirs)
5588 self._check_store_dirs_type()
5589 self._id_to_store_dir = {}
5590 self._open_stores = {}
5591 self._effective_default_store_id = None
5593 def _check_store_dirs_type(self):
5594 for sdir in ['store_dirs', 'store_superdirs']:
5595 if not isinstance(self.__getattribute__(sdir), list):
5596 raise TypeError('{} of {} is not of type list'.format(
5597 sdir, self.__class__.__name__))
5599 def _get_store_id(self, store_dir):
5600 store_ = store.Store(store_dir)
5601 store_id = store_.config.id
5602 store_.close()
5603 return store_id
5605 def _looks_like_store_dir(self, store_dir):
5606 return os.path.isdir(store_dir) and \
5607 all(os.path.isfile(pjoin(store_dir, x)) for x in
5608 ('index', 'traces', 'config'))
5610 def iter_store_dirs(self):
5611 store_dirs = set()
5612 for d in self.store_superdirs:
5613 if not os.path.exists(d):
5614 logger.warning('store_superdir not available: %s' % d)
5615 continue
5617 for entry in os.listdir(d):
5618 store_dir = os.path.realpath(pjoin(d, entry))
5619 if self._looks_like_store_dir(store_dir):
5620 store_dirs.add(store_dir)
5622 for store_dir in self.store_dirs:
5623 store_dirs.add(os.path.realpath(store_dir))
5625 return store_dirs
5627 def _scan_stores(self):
5628 for store_dir in self.iter_store_dirs():
5629 store_id = self._get_store_id(store_dir)
5630 if store_id not in self._id_to_store_dir:
5631 self._id_to_store_dir[store_id] = store_dir
5632 else:
5633 if store_dir != self._id_to_store_dir[store_id]:
5634 raise DuplicateStoreId(
5635 'GF store ID %s is used in (at least) two '
5636 'different stores. Locations are: %s and %s' %
5637 (store_id, self._id_to_store_dir[store_id], store_dir))
5639 def get_store_dir(self, store_id):
5640 '''
5641 Lookup directory given a GF store ID.
5642 '''
5644 if store_id not in self._id_to_store_dir:
5645 self._scan_stores()
5647 if store_id not in self._id_to_store_dir:
5648 raise NoSuchStore(store_id, self.iter_store_dirs())
5650 return self._id_to_store_dir[store_id]
5652 def get_store_ids(self):
5653 '''
5654 Get list of available store IDs.
5655 '''
5657 self._scan_stores()
5658 return sorted(self._id_to_store_dir.keys())
5660 def effective_default_store_id(self):
5661 if self._effective_default_store_id is None:
5662 if self.default_store_id is None:
5663 store_ids = self.get_store_ids()
5664 if len(store_ids) == 1:
5665 self._effective_default_store_id = self.get_store_ids()[0]
5666 else:
5667 raise NoDefaultStoreSet()
5668 else:
5669 self._effective_default_store_id = self.default_store_id
5671 return self._effective_default_store_id
5673 def get_store(self, store_id=None):
5674 '''
5675 Get a store from the engine.
5677 :param store_id: identifier of the store (optional)
5678 :returns: :py:class:`~pyrocko.gf.store.Store` object
5680 If no ``store_id`` is provided the store
5681 associated with the :py:gattr:`default_store_id` is returned.
5682 Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is
5683 undefined.
5684 '''
5686 if store_id is None:
5687 store_id = self.effective_default_store_id()
5689 if store_id not in self._open_stores:
5690 store_dir = self.get_store_dir(store_id)
5691 self._open_stores[store_id] = store.Store(store_dir)
5693 return self._open_stores[store_id]
5695 def get_store_config(self, store_id):
5696 store = self.get_store(store_id)
5697 return store.config
5699 def get_store_extra(self, store_id, key):
5700 store = self.get_store(store_id)
5701 return store.get_extra(key)
5703 def close_cashed_stores(self):
5704 '''
5705 Close and remove ids from cashed stores.
5706 '''
5707 store_ids = []
5708 for store_id, store_ in self._open_stores.items():
5709 store_.close()
5710 store_ids.append(store_id)
5712 for store_id in store_ids:
5713 self._open_stores.pop(store_id)
5715 def get_rule(self, source, target):
5716 cprovided = self.get_store(target.store_id).get_provided_components()
5718 if isinstance(target, StaticTarget):
5719 quantity = target.quantity
5720 available_rules = static_rules
5721 elif isinstance(target, Target):
5722 quantity = target.effective_quantity()
5723 available_rules = channel_rules
5725 try:
5726 for rule in available_rules[quantity]:
5727 cneeded = rule.required_components(target)
5728 if all(c in cprovided for c in cneeded):
5729 return rule
5731 except KeyError:
5732 pass
5734 raise BadRequest(
5735 'No rule to calculate "%s" with GFs from store "%s" '
5736 'for source model "%s".' % (
5737 target.effective_quantity(),
5738 target.store_id,
5739 source.__class__.__name__))
5741 def _cached_discretize_basesource(self, source, store, cache, target):
5742 if (source, store) not in cache:
5743 cache[source, store] = source.discretize_basesource(store, target)
5745 return cache[source, store]
5747 def base_seismograms(self, source, targets, components, dsource_cache,
5748 nthreads=0):
5750 target = targets[0]
5752 interp = set([t.interpolation for t in targets])
5753 if len(interp) > 1:
5754 raise BadRequest('Targets have different interpolation schemes.')
5756 rates = set([t.sample_rate for t in targets])
5757 if len(rates) > 1:
5758 raise BadRequest('Targets have different sample rates.')
5760 store_ = self.get_store(target.store_id)
5761 receivers = [t.receiver(store_) for t in targets]
5763 if target.sample_rate is not None:
5764 deltat = 1. / target.sample_rate
5765 rate = target.sample_rate
5766 else:
5767 deltat = None
5768 rate = store_.config.sample_rate
5770 tmin = num.fromiter(
5771 (t.tmin for t in targets), dtype=float, count=len(targets))
5772 tmax = num.fromiter(
5773 (t.tmax for t in targets), dtype=float, count=len(targets))
5775 mask = num.logical_and(num.isfinite(tmin), num.isfinite(tmax))
5777 itmin = num.zeros_like(tmin, dtype=num.int64)
5778 itmax = num.zeros_like(tmin, dtype=num.int64)
5779 nsamples = num.full_like(tmin, -1, dtype=num.int64)
5781 itmin[mask] = num.floor(tmin[mask] * rate).astype(num.int64)
5782 itmax[mask] = num.ceil(tmax[mask] * rate).astype(num.int64)
5783 nsamples = itmax - itmin + 1
5784 nsamples[num.logical_not(mask)] = -1
5786 base_source = self._cached_discretize_basesource(
5787 source, store_, dsource_cache, target)
5789 base_seismograms = store_.calc_seismograms(
5790 base_source, receivers, components,
5791 deltat=deltat,
5792 itmin=itmin, nsamples=nsamples,
5793 interpolation=target.interpolation,
5794 optimization=target.optimization,
5795 nthreads=nthreads)
5797 for i, base_seismogram in enumerate(base_seismograms):
5798 base_seismograms[i] = store.make_same_span(base_seismogram)
5800 return base_seismograms
5802 def base_seismogram(self, source, target, components, dsource_cache,
5803 nthreads):
5805 tcounters = [xtime()]
5807 store_ = self.get_store(target.store_id)
5808 receiver = target.receiver(store_)
5810 if target.tmin and target.tmax is not None:
5811 rate = store_.config.sample_rate
5812 itmin = int(num.floor(target.tmin * rate))
5813 itmax = int(num.ceil(target.tmax * rate))
5814 nsamples = itmax - itmin + 1
5815 else:
5816 itmin = None
5817 nsamples = None
5819 tcounters.append(xtime())
5820 base_source = self._cached_discretize_basesource(
5821 source, store_, dsource_cache, target)
5823 tcounters.append(xtime())
5825 if target.sample_rate is not None:
5826 deltat = 1. / target.sample_rate
5827 else:
5828 deltat = None
5830 base_seismogram = store_.seismogram(
5831 base_source, receiver, components,
5832 deltat=deltat,
5833 itmin=itmin, nsamples=nsamples,
5834 interpolation=target.interpolation,
5835 optimization=target.optimization,
5836 nthreads=nthreads)
5838 tcounters.append(xtime())
5840 base_seismogram = store.make_same_span(base_seismogram)
5842 tcounters.append(xtime())
5844 return base_seismogram, tcounters
5846 def base_statics(self, source, target, components, nthreads):
5847 tcounters = [xtime()]
5848 store_ = self.get_store(target.store_id)
5850 if target.tsnapshot is not None:
5851 rate = store_.config.sample_rate
5852 itsnapshot = int(num.floor(target.tsnapshot * rate))
5853 else:
5854 itsnapshot = None
5855 tcounters.append(xtime())
5857 base_source = source.discretize_basesource(store_, target=target)
5859 tcounters.append(xtime())
5861 base_statics = store_.statics(
5862 base_source,
5863 target,
5864 itsnapshot,
5865 components,
5866 target.interpolation,
5867 nthreads)
5869 tcounters.append(xtime())
5871 return base_statics, tcounters
5873 def _post_process_dynamic(self, base_seismogram, source, target):
5874 base_any = next(iter(base_seismogram.values()))
5875 deltat = base_any.deltat
5876 itmin = base_any.itmin
5878 rule = self.get_rule(source, target)
5879 data = rule.apply_(target, base_seismogram)
5881 factor = source.get_factor() * target.get_factor()
5882 if factor != 1.0:
5883 data = data * factor
5885 stf = source.effective_stf_post()
5887 times, amplitudes = stf.discretize_t(
5888 deltat, 0.0)
5890 # repeat end point to prevent boundary effects
5891 padded_data = num.empty(data.size + amplitudes.size, dtype=float)
5892 padded_data[:data.size] = data
5893 padded_data[data.size:] = data[-1]
5894 data = num.convolve(amplitudes, padded_data)
5896 tmin = itmin * deltat + times[0]
5898 tr = meta.SeismosizerTrace(
5899 codes=target.codes,
5900 data=data[:-amplitudes.size],
5901 deltat=deltat,
5902 tmin=tmin)
5904 return target.post_process(self, source, tr)
5906 def _post_process_statics(self, base_statics, source, starget):
5907 rule = self.get_rule(source, starget)
5908 data = rule.apply_(starget, base_statics)
5910 factor = source.get_factor()
5911 if factor != 1.0:
5912 for v in data.values():
5913 v *= factor
5915 return starget.post_process(self, source, base_statics)
5917 def process(self, *args, **kwargs):
5918 '''
5919 Process a request.
5921 ::
5923 process(**kwargs)
5924 process(request, **kwargs)
5925 process(sources, targets, **kwargs)
5927 The request can be given a a :py:class:`Request` object, or such an
5928 object is created using ``Request(**kwargs)`` for convenience.
5930 :returns: :py:class:`Response` object
5931 '''
5933 if len(args) not in (0, 1, 2):
5934 raise BadRequest('Invalid arguments.')
5936 if len(args) == 1:
5937 kwargs['request'] = args[0]
5939 elif len(args) == 2:
5940 kwargs.update(Request.args2kwargs(args))
5942 request = kwargs.pop('request', None)
5943 status_callback = kwargs.pop('status_callback', None)
5944 calc_timeseries = kwargs.pop('calc_timeseries', True)
5946 nprocs = kwargs.pop('nprocs', None)
5947 nthreads = kwargs.pop('nthreads', 1)
5948 if nprocs is not None:
5949 nthreads = nprocs
5951 if request is None:
5952 request = Request(**kwargs)
5954 if resource:
5955 rs0 = resource.getrusage(resource.RUSAGE_SELF)
5956 rc0 = resource.getrusage(resource.RUSAGE_CHILDREN)
5957 tt0 = xtime()
5959 # make sure stores are open before fork()
5960 store_ids = set(target.store_id for target in request.targets)
5961 for store_id in store_ids:
5962 self.get_store(store_id)
5964 source_index = dict((x, i) for (i, x) in
5965 enumerate(request.sources))
5966 target_index = dict((x, i) for (i, x) in
5967 enumerate(request.targets))
5969 m = request.subrequest_map()
5971 skeys = sorted(m.keys(), key=cmp_to_key(cmp_none_aware))
5972 results_list = []
5974 for i in range(len(request.sources)):
5975 results_list.append([None] * len(request.targets))
5977 tcounters_dyn_list = []
5978 tcounters_static_list = []
5979 nsub = len(skeys)
5980 isub = 0
5982 # Processing dynamic targets through
5983 # parimap(process_subrequest_dynamic)
5985 if calc_timeseries:
5986 _process_dynamic = process_dynamic_timeseries
5987 else:
5988 _process_dynamic = process_dynamic
5990 if request.has_dynamic:
5991 work_dynamic = [
5992 (i, nsub,
5993 [source_index[source] for source in m[k][0]],
5994 [target_index[target] for target in m[k][1]
5995 if not isinstance(target, StaticTarget)])
5996 for (i, k) in enumerate(skeys)]
5998 for ii_results, tcounters_dyn in _process_dynamic(
5999 work_dynamic, request.sources, request.targets, self,
6000 nthreads):
6002 tcounters_dyn_list.append(num.diff(tcounters_dyn))
6003 isource, itarget, result = ii_results
6004 results_list[isource][itarget] = result
6006 if status_callback:
6007 status_callback(isub, nsub)
6009 isub += 1
6011 # Processing static targets through process_static
6012 if request.has_statics:
6013 work_static = [
6014 (i, nsub,
6015 [source_index[source] for source in m[k][0]],
6016 [target_index[target] for target in m[k][1]
6017 if isinstance(target, StaticTarget)])
6018 for (i, k) in enumerate(skeys)]
6020 for ii_results, tcounters_static in process_static(
6021 work_static, request.sources, request.targets, self,
6022 nthreads=nthreads):
6024 tcounters_static_list.append(num.diff(tcounters_static))
6025 isource, itarget, result = ii_results
6026 results_list[isource][itarget] = result
6028 if status_callback:
6029 status_callback(isub, nsub)
6031 isub += 1
6033 if status_callback:
6034 status_callback(nsub, nsub)
6036 tt1 = time.time()
6037 if resource:
6038 rs1 = resource.getrusage(resource.RUSAGE_SELF)
6039 rc1 = resource.getrusage(resource.RUSAGE_CHILDREN)
6041 s = ProcessingStats()
6043 if request.has_dynamic:
6044 tcumu_dyn = num.sum(num.vstack(tcounters_dyn_list), axis=0)
6045 t_dyn = float(num.sum(tcumu_dyn))
6046 perc_dyn = map(float, tcumu_dyn / t_dyn * 100.)
6047 (s.t_perc_get_store_and_receiver,
6048 s.t_perc_discretize_source,
6049 s.t_perc_make_base_seismogram,
6050 s.t_perc_make_same_span,
6051 s.t_perc_post_process) = perc_dyn
6052 else:
6053 t_dyn = 0.
6055 if request.has_statics:
6056 tcumu_static = num.sum(num.vstack(tcounters_static_list), axis=0)
6057 t_static = num.sum(tcumu_static)
6058 perc_static = map(float, tcumu_static / t_static * 100.)
6059 (s.t_perc_static_get_store,
6060 s.t_perc_static_discretize_basesource,
6061 s.t_perc_static_sum_statics,
6062 s.t_perc_static_post_process) = perc_static
6064 s.t_wallclock = tt1 - tt0
6065 if resource:
6066 s.t_cpu = (
6067 (rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) -
6068 (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime))
6069 s.n_read_blocks = (
6070 (rs1.ru_inblock + rc1.ru_inblock) -
6071 (rs0.ru_inblock + rc0.ru_inblock))
6073 n_records_stacked = 0.
6074 for results in results_list:
6075 for result in results:
6076 if not isinstance(result, meta.Result):
6077 continue
6078 shr = float(result.n_shared_stacking)
6079 n_records_stacked += result.n_records_stacked / shr
6080 s.t_perc_optimize += result.t_optimize / shr
6081 s.t_perc_stack += result.t_stack / shr
6082 s.n_records_stacked = int(n_records_stacked)
6083 if t_dyn != 0.:
6084 s.t_perc_optimize /= t_dyn * 100
6085 s.t_perc_stack /= t_dyn * 100
6087 return Response(
6088 request=request,
6089 results_list=results_list,
6090 stats=s)
6093class RemoteEngine(Engine):
6094 '''
6095 Client for remote synthetic seismogram calculator.
6096 '''
6098 site = String.T(default=ws.g_default_site, optional=True)
6099 url = String.T(default=ws.g_url, optional=True)
6101 def process(self, request=None, status_callback=None, **kwargs):
6103 if request is None:
6104 request = Request(**kwargs)
6106 return ws.seismosizer(url=self.url, site=self.site, request=request)
6109g_engine = None
6112def get_engine(store_superdirs=[]):
6113 global g_engine
6114 if g_engine is None:
6115 g_engine = LocalEngine(use_env=True, use_config=True)
6117 for d in store_superdirs:
6118 if d not in g_engine.store_superdirs:
6119 g_engine.store_superdirs.append(d)
6121 return g_engine
6124class SourceGroup(Object):
6126 def __getattr__(self, k):
6127 return num.fromiter((getattr(s, k) for s in self),
6128 dtype=float)
6130 def __iter__(self):
6131 raise NotImplementedError(
6132 'This method should be implemented in subclass.')
6134 def __len__(self):
6135 raise NotImplementedError(
6136 'This method should be implemented in subclass.')
6139class SourceList(SourceGroup):
6140 sources = List.T(Source.T())
6142 def append(self, s):
6143 self.sources.append(s)
6145 def __iter__(self):
6146 return iter(self.sources)
6148 def __len__(self):
6149 return len(self.sources)
6152class SourceGrid(SourceGroup):
6154 base = Source.T()
6155 variables = Dict.T(String.T(), Range.T())
6156 order = List.T(String.T())
6158 def __len__(self):
6159 n = 1
6160 for (k, v) in self.make_coords(self.base):
6161 n *= len(list(v))
6163 return n
6165 def __iter__(self):
6166 for items in permudef(self.make_coords(self.base)):
6167 s = self.base.clone(**{k: v for (k, v) in items})
6168 s.regularize()
6169 yield s
6171 def ordered_params(self):
6172 ks = list(self.variables.keys())
6173 for k in self.order + list(self.base.keys()):
6174 if k in ks:
6175 yield k
6176 ks.remove(k)
6177 if ks:
6178 raise Exception('Invalid parameter "%s" for source type "%s".' %
6179 (ks[0], self.base.__class__.__name__))
6181 def make_coords(self, base):
6182 return [(param, self.variables[param].make(base=base[param]))
6183 for param in self.ordered_params()]
6186source_classes = [
6187 Source,
6188 SourceWithMagnitude,
6189 SourceWithDerivedMagnitude,
6190 ExplosionSource,
6191 RectangularExplosionSource,
6192 DCSource,
6193 CLVDSource,
6194 VLVDSource,
6195 MTSource,
6196 RectangularSource,
6197 PseudoDynamicRupture,
6198 DoubleDCSource,
6199 RingfaultSource,
6200 CombiSource,
6201 CombiSFSource,
6202 SFSource,
6203 SimpleLandslideSource,
6204 PorePressurePointSource,
6205 PorePressureLineSource,
6206]
6208stf_classes = [
6209 STF,
6210 BoxcarSTF,
6211 TriangularSTF,
6212 HalfSinusoidSTF,
6213 ResonatorSTF,
6214 TremorSTF,
6215 SimpleLandslideSTF,
6216]
6218__all__ = '''
6219Cloneable
6220NoDefaultStoreSet
6221SeismosizerError
6222STFError
6223BadRequest
6224NoSuchStore
6225DerivedMagnitudeError
6226STFMode
6227'''.split() + [S.__name__ for S in source_classes + stf_classes] + '''
6228Request
6229ProcessingStats
6230Response
6231Engine
6232LocalEngine
6233RemoteEngine
6234source_classes
6235get_engine
6236Range
6237SourceGroup
6238SourceList
6239SourceGrid
6240map_anchor
6241'''.split()