# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
Timestamp, Int, SObject, ArgumentError, Dict, ValidationError)
if isinstance(a, tuple) and isinstance(b, tuple): for xa, xb in zip(a, b): rv = cmp_none_aware(xa, xb) if rv != 0: return rv
return 0
anone = a is None bnone = b is None
if anone and bnone: return 0
if anone: return -1
if bnone: return 1
return bool(a > b) - bool(a < b)
BadRequest.__init__(self) self.store_id = store_id self.dirs = dirs
def __str__(self): if self.store_id is not None: rstr = 'no GF store with id "%s" found.' % self.store_id else: rstr = 'GF store not found.'
if self.dirs is not None: rstr += ' Searched folders:\n %s' % '\n '.join(sorted(self.dirs)) return rstr
units = { 'k': 1e3, 'M': 1e6, }
factor = 1.0 if s and s[-1] in units: factor = units[s[-1]] s = s[:-1] if not s: raise ValueError('unit without a number: \'%s\'' % s)
return float(s) * factor
if s: return ufloat(s) else: return None
if s: return int(s) else: return None
if j < len(ln): k, v = ln[j] for y in v: ln[j] = k, y for s in permudef(ln, j + 1): yield s
ln[j] = k, v return else: yield ln
return num.atleast_1d(num.asarray(x))
strike, dip, length, width, anchor, velocity, stf=None, nucleation_x=None, nucleation_y=None, decimation_factor=1):
stf = STF()
dist_x = num.abs(nucleation_x - points[:, 0]) else:
dist_y = num.abs(nucleation_y - points[:, 1]) else:
pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
# We assume a non-rotated fault plane N_CRITICAL = 8 points = points2.T.reshape((3, nl, nw)) if points.size <= N_CRITICAL: logger.warning('RectangularSource is defined by only %d sub-sources!' % points.size) return True
distances = num.sqrt( (points[0, 0, :] - points[0, 1, :])**2 + (points[1, 0, :] - points[1, 1, :])**2 + (points[2, 0, :] - points[2, 1, :])**2)
depths = points[2, 0, :] vs_profile = store.config.get_vs( lat=0., lon=0., points=num.repeat(depths[:, num.newaxis], 3, axis=1), interpolation='multilinear')
min_wavelength = vs_profile * (store.config.deltat * 2) if not num.all(min_wavelength > distances/2): return False return True
[[-0.5 * ln, -0.5 * wd, 0.], [0.5 * ln, -0.5 * wd, 0.], [0.5 * ln, 0.5 * wd, 0.], [-0.5 * ln, 0.5 * wd, 0.], [-0.5 * ln, -0.5 * wd, 0.]])
pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
strike, dip, length, width, depth, x_plane_coords, y_plane_coords, lat=0., lon=0., north_shift=0, east_shift=0, anchor='top', cs='xy'):
ln = length wd = width x_abs = [] y_abs = [] if not isinstance(x_plane_coords, list): x_plane_coords = [x_plane_coords] y_plane_coords = [y_plane_coords]
for x_plane, y_plane in zip(x_plane_coords, y_plane_coords): points = num.array( [[-0.5 * ln*x_plane, -0.5 * wd*y_plane, 0.], [0.5 * ln*x_plane, -0.5 * wd*y_plane, 0.], [0.5 * ln*x_plane, 0.5 * wd*y_plane, 0.], [-0.5 * ln*x_plane, 0.5 * wd*y_plane, 0.], [-0.5 * ln*x_plane, -0.5 * wd*y_plane, 0.]])
anch_x, anch_y = map_anchor[anchor] points[:, 0] -= anch_x * 0.5 * length points[:, 1] -= anch_y * 0.5 * width
rotmat = num.asarray( pmt.euler_to_matrix(dip * d2r, strike * d2r, 0.0))
points = num.dot(rotmat.T, points.T).T points[:, 0] += north_shift points[:, 1] += east_shift points[:, 2] += depth if cs in ('latlon', 'lonlat'): latlon = ne_to_latlon(lat, lon, points[:, 0], points[:, 1]) latlon = num.array(latlon).T x_abs.append(latlon[1:2, 1]) y_abs.append(latlon[2:3, 0]) if cs == 'xy': x_abs.append(points[1:2, 1]) y_abs.append(points[2:3, 0])
if cs == 'lonlat': return y_abs, x_abs else: return x_abs, y_abs
''' Convenient range specification.
Equivalent ways to sepecify the range [ 0., 1000., ... 10000. ]::
Range('0 .. 10k : 1k') Range(start=0., stop=10e3, step=1e3) Range(0, 10e3, 1e3) Range('0 .. 10k @ 11') Range(start=0., stop=10*km, n=11)
Range(0, 10e3, n=11) Range(values=[x*1e3 for x in range(11)])
Depending on the use context, it can be possible to omit any part of the specification. E.g. in the context of extracting a subset of an already existing range, the existing range's specification values would be filled in where missing.
The values are distributed with equal spacing, unless the ``spacing`` argument is modified. The values can be created offset or relative to an external base value with the ``relative`` argument if the use context supports this.
The range specification can be expressed with a short string representation::
'start .. stop @ num | spacing, relative' 'start .. stop : step | spacing, relative'
most parts of the expression can be omitted if not needed. Whitespace is allowed for readability but can also be omitted. '''
choices=['lin', 'log', 'symlog'], default='lin', optional=True)
choices=['', 'add', 'mult'], default='', optional=True)
r'(@(?P<n>[^|]+)|:(?P<step>[^|]+))?' r'(\|(?P<stuff>.+))?$')
d = {} if len(args) == 1: d = self.parse(args[0]) elif len(args) in (2, 3): d['start'], d['stop'] = [float(x) for x in args[:2]] if len(args) == 3: d['step'] = float(args[2])
for k, v in kwargs.items(): if k in d: raise ArgumentError('%s specified more than once' % k)
d[k] = v
SObject.__init__(self, **d)
def __str__(self): def sfloat(x): if x is not None: return '%g' % x else: return ''
if self.values: return ','.join('%g' % x for x in self.values)
if self.start is None and self.stop is None: s0 = '' else: s0 = '%s .. %s' % (sfloat(self.start), sfloat(self.stop))
s1 = '' if self.step is not None: s1 = [' : %g', ':%g'][s0 == ''] % self.step elif self.n is not None: s1 = [' @ %i', '@%i'][s0 == ''] % self.n
if self.spacing == 'lin' and self.relative == '': s2 = '' else: x = [] if self.spacing != 'lin': x.append(self.spacing)
if self.relative != '': x.append(self.relative)
s2 = ' | %s' % ','.join(x)
return s0 + s1 + s2
def parse(cls, s): s = re.sub(r'\s+', '', s) m = cls.pattern.match(s) if not m: try: vals = [ufloat(x) for x in s.split(',')] except Exception: raise InvalidGridDef( '"%s" is not a valid range specification' % s)
return dict(values=num.array(vals, dtype=num.float))
d = m.groupdict() try: start = ufloat_or_none(d['start']) stop = ufloat_or_none(d['stop']) step = ufloat_or_none(d['step']) n = int_or_none(d['n']) except Exception: raise InvalidGridDef( '"%s" is not a valid range specification' % s)
spacing = 'lin' relative = ''
if d['stuff'] is not None: t = d['stuff'].split(',') for x in t: if x in cls.spacing.choices: spacing = x elif x and x in cls.relative.choices: relative = x else: raise InvalidGridDef( '"%s" is not a valid range specification' % s)
return dict(start=start, stop=stop, step=step, n=n, spacing=spacing, relative=relative)
if self.values: return self.values
start = self.start stop = self.stop step = self.step n = self.n
swap = step is not None and step < 0. if start is None: start = [mi, ma][swap] if stop is None: stop = [ma, mi][swap] if step is None and inc is not None: step = [inc, -inc][ma < mi]
if start is None or stop is None: raise InvalidGridDef( 'Cannot use range specification "%s" without start ' 'and stop in this context' % self)
if step is None and n is None: step = stop - start
if n is None: if (step < 0) != (stop - start < 0): raise InvalidGridDef( 'Range specification "%s" has inconsistent ordering ' '(step < 0 => stop > start)' % self)
n = int(round((stop - start) / step)) + 1 stop2 = start + (n - 1) * step if abs(stop - stop2) > eps: n = int(math.floor((stop - start) / step)) + 1 stop = start + (n - 1) * step else: stop = stop2
if start == stop: n = 1
if self.spacing == 'lin': vals = num.linspace(start, stop, n)
elif self.spacing in ('log', 'symlog'): if start > 0. and stop > 0.: vals = num.exp(num.linspace(num.log(start), num.log(stop), n)) elif start < 0. and stop < 0.: vals = -num.exp(num.linspace(num.log(-start), num.log(-stop), n)) else: raise InvalidGridDef( 'log ranges should not include or cross zero ' '(in range specification "%s")' % self)
if self.spacing == 'symlog': nvals = - vals vals = num.concatenate((nvals[::-1], vals))
if self.relative in ('add', 'mult') and base is None: raise InvalidGridDef( 'cannot use relative range specification in this context')
vals = self.make_relative(base, vals)
return list(map(float, vals))
if self.relative == 'add': vals += base
if self.relative == 'mult': vals *= base
return vals
if shorthand is not None: t = shorthand.split('=') if len(t) != 2: raise InvalidGridDef( 'invalid grid specification element: %s' % shorthand)
sp, sr = t[0].strip(), t[1].strip()
kwargs['param'] = sp kwargs['rs'] = Range(sr)
Object.__init__(self, **kwargs)
return self.param + ' = ' + str(self.rs)
if shorthand is not None: t = shorthand.splitlines() tt = [] for x in t: x = x.strip() if x: tt.extend(x.split(';'))
elements = [] for se in tt: elements.append(GridDef(se))
kwargs['elements'] = elements
Object.__init__(self, **kwargs)
return '; '.join(str(x) for x in self.elements)
return iter(self.T.propnames)
if k not in self.keys(): raise KeyError(k)
return getattr(self, k)
if k not in self.keys(): raise KeyError(k)
return setattr(self, k, v)
''' Make a copy of the object.
A new object of the same class is created and initialized with the parameters of the object on which this method is called on. If ``kwargs`` are given, these are used to override any of the initialization parameters. '''
d = dict(self) for k in d: v = d[k] if isinstance(v, Cloneable): d[k] = v.clone()
d.update(kwargs) return self.__class__(**d)
def keys(cls): ''' Get list of the source model's parameter names. '''
return cls.T.propnames
''' Base class for source time functions. '''
kwargs['duration'] = effective_duration / \ self.factor_duration_to_effective()
def factor_duration_to_effective(cls): return 1.0
return tref
def effective_duration(self): return self.duration * self.factor_duration_to_effective()
else: return ( num.array([tl, th], dtype=num.float), num.array([th - tref, tref - tl], dtype=num.float) / deltat)
t0 = math.floor(tshift / deltat) * deltat t1 = math.ceil(tshift / deltat) * deltat if t0 == t1: return times, amplitudes
amplitudes2 = num.zeros(amplitudes.size + 1, dtype=num.float)
amplitudes2[:-1] += (t1 - tshift) / deltat * amplitudes amplitudes2[1:] += (tshift - t0) / deltat * amplitudes
times2 = num.arange(times.size + 1, dtype=num.float) * \ deltat + times[0] + t0
return times2, amplitudes2
''' Boxcar type source time function.
.. figure :: /static/stf-BoxcarSTF.svg :width: 40% :align: center :alt: boxcar source time function '''
default=0.0, help='duration of the boxcar')
default=0.0, help='anchor point with respect to source.time: (' '-1.0: left -> source duration [0, T] ~ hypocenter time, ' ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
def factor_duration_to_effective(cls): return 1.0
return tref - 0.5 * self.duration * self.anchor
tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 tmin = round(tmin_stf / deltat) * deltat tmax = round(tmax_stf / deltat) * deltat nt = int(round((tmax - tmin) / deltat)) + 1 times = num.linspace(tmin, tmax, nt) amplitudes = num.ones_like(times) if times.size > 1: t_edges = num.linspace( tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1) t = tmin_stf + self.duration * num.array( [0.0, 0.0, 1.0, 1.0], dtype=num.float) f = num.array([0., 1., 1., 0.], dtype=num.float) amplitudes = util.plf_integrate_piecewise(t_edges, t, f) amplitudes /= num.sum(amplitudes)
tshift = (num.sum(amplitudes * times) - self.centroid_time(tref))
return sshift(times, amplitudes, -tshift, deltat)
return (type(self).__name__, self.duration, self.anchor)
''' Triangular type source time function.
.. figure :: /static/stf-TriangularSTF.svg :width: 40% :align: center :alt: triangular source time function '''
default=0.0, help='baseline of the triangle')
default=0.5, help='fraction of time compared to duration, ' 'when the maximum amplitude is reached')
default=0.0, help='anchor point with respect to source.time: (' '-1.0: left -> source duration [0, T] ~ hypocenter time, ' ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
if peak_ratio is None: peak_ratio = cls.peak_ratio.default()
return math.sqrt((peak_ratio**2 - peak_ratio + 1.0) * 2.0 / 3.0)
kwargs['duration'] = effective_duration / \ self.factor_duration_to_effective( kwargs.get('peak_ratio', None))
def centroid_ratio(self): ra = self.peak_ratio rb = 1.0 - ra return self.peak_ratio + (rb**2 / 3. - ra**2 / 3.) / (ra + rb)
ca = self.centroid_ratio cb = 1.0 - ca if self.anchor <= 0.: return tref - ca * self.duration * self.anchor else: return tref - cb * self.duration * self.anchor
def effective_duration(self): return self.duration * self.factor_duration_to_effective( self.peak_ratio)
ca = self.centroid_ratio cb = 1.0 - ca if self.anchor <= 0.: tmin_stf = tref - ca * self.duration * (self.anchor + 1.) tmax_stf = tmin_stf + self.duration else: tmax_stf = tref + cb * self.duration * (1. - self.anchor) tmin_stf = tmax_stf - self.duration
return tmin_stf, tmax_stf
tmin_stf, tmax_stf = self.tminmax_stf(tref)
tmin = round(tmin_stf / deltat) * deltat tmax = round(tmax_stf / deltat) * deltat nt = int(round((tmax - tmin) / deltat)) + 1 if nt > 1: t_edges = num.linspace( tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1) t = tmin_stf + self.duration * num.array( [0.0, self.peak_ratio, 1.0], dtype=num.float) f = num.array([0., 1., 0.], dtype=num.float) amplitudes = util.plf_integrate_piecewise(t_edges, t, f) amplitudes /= num.sum(amplitudes) else: amplitudes = num.ones(1)
times = num.linspace(tmin, tmax, nt) return times, amplitudes
return ( type(self).__name__, self.duration, self.peak_ratio, self.anchor)
''' Half sinusoid type source time function.
.. figure :: /static/stf-HalfSinusoidSTF.svg :width: 40% :align: center :alt: half-sinusouid source time function '''
default=0.0, help='duration of the half-sinusoid (baseline)')
default=0.0, help='anchor point with respect to source.time: (' '-1.0: left -> source duration [0, T] ~ hypocenter time, ' ' 0.0: center -> source duration [-T/2, T/2] ~ centroid time, ' '+1.0: right -> source duration [-T, 0] ~ rupture end time)')
default=1, help='set to 2 to use square of the half-period sinusoidal function.')
kwargs['duration'] = effective_duration / \ self.factor_duration_to_effective( kwargs.get('exponent', 1))
def factor_duration_to_effective(cls, exponent): if exponent == 1: return math.sqrt(3.0 * math.pi**2 - 24.0) / math.pi elif exponent == 2: return math.sqrt(math.pi**2 - 6) / math.pi else: raise ValueError('exponent for HalfSinusoidSTF must be 1 or 2')
def effective_duration(self): return self.duration * self.factor_duration_to_effective(self.exponent)
return tref - 0.5 * self.duration * self.anchor
tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 tmin = round(tmin_stf / deltat) * deltat tmax = round(tmax_stf / deltat) * deltat nt = int(round((tmax - tmin) / deltat)) + 1 if nt > 1: t_edges = num.maximum(tmin_stf, num.minimum(tmax_stf, num.linspace( tmin - 0.5 * deltat, tmax + 0.5 * deltat, nt + 1)))
if self.exponent == 1: fint = -num.cos( (t_edges - tmin_stf) * (math.pi / self.duration))
elif self.exponent == 2: fint = (t_edges - tmin_stf) / self.duration \ - 1.0 / (2.0 * math.pi) * num.sin( (t_edges - tmin_stf) * (2.0 * math.pi / self.duration)) else: raise ValueError('exponent for HalfSinusoidSTF must be 1 or 2')
amplitudes = fint[1:] - fint[:-1] amplitudes /= num.sum(amplitudes) else: amplitudes = num.ones(1)
times = num.linspace(tmin, tmax, nt) return times, amplitudes
return (type(self).__name__, self.duration, self.anchor)
'''Smooth-ramp type source time function for near-field displacement. Based on moment function of double-couple point source proposed by Bruestle and Mueller (PEPI, 1983).
.. [1] W. Bruestle, G. Mueller (1983), Moment and duration of shallow earthquakes from Love-wave modelling for regional distances, PEPI 32, 312-324.
.. figure :: /static/stf-SmoothRampSTF.svg :width: 40% :alt: smooth ramp source time function ''' default=0.0, help='duration of the ramp (baseline)')
default=0.5, help='fraction of time compared to duration, ' 'when the maximum amplitude is reached')
default=0.0, help='anchor point with respect to source.time: (' '-1.0: left -> source duration ``[0, T]`` ~ hypocenter time, ' '0.0: center -> source duration ``[-T/2, T/2]`` ~ centroid time, ' '+1.0: right -> source duration ``[-T, 0]`` ~ rupture end time)')
tmin_stf = tref - self.duration * (self.anchor + 1.) * 0.5 tmax_stf = tref + self.duration * (1. - self.anchor) * 0.5 tmin = round(tmin_stf / deltat) * deltat tmax = round(tmax_stf / deltat) * deltat D = round((tmax - tmin) / deltat) * deltat nt = int(round(D / deltat)) + 1 times = num.linspace(tmin, tmax, nt) if nt > 1: rise_time = self.rise_ratio * self.duration amplitudes = num.ones_like(times) tp = tmin + rise_time ii = num.where(times <= tp) t_inc = times[ii] a = num.cos(num.pi * (t_inc - tmin_stf) / rise_time) b = num.cos(3 * num.pi * (t_inc - tmin_stf) / rise_time) - 1.0 amplitudes[ii] = (9. / 16.) * (1 - a + (1. / 9.) * b)
amplitudes /= num.sum(amplitudes) else: amplitudes = num.ones(1)
return times, amplitudes
return (type(self).__name__, self.duration, self.rise_ratio, self.anchor)
''' Simple resonator like source time function.
.. math ::
f(t) = 0 for t < 0 f(t) = e^{-t/tau} * sin(2 * pi * f * t)
.. figure :: /static/stf-SmoothRampSTF.svg :width: 40% :alt: smooth ramp source time function
'''
default=0.0, help='decay time')
default=1.0, help='resonance frequency')
tmin_stf = tref tmax_stf = tref + self.duration * 3 tmin = math.floor(tmin_stf / deltat) * deltat tmax = math.ceil(tmax_stf / deltat) * deltat times = util.arange2(tmin, tmax, deltat) amplitudes = num.exp(-(times-tref)/self.duration) \ * num.sin(2.0 * num.pi * self.frequency * (times-tref))
return times, amplitudes
return (type(self).__name__, self.duration, self.frequency)
''' Base class for all source models. '''
default=0., help='source origin time.')
optional=True, help='source time function.')
default='post', help='whether to apply source time function in pre or ' 'post-processing.')
''' Change some of the source models parameters.
Example::
>>> from pyrocko import gf >>> s = gf.DCSource() >>> s.update(strike=66., dip=33.) >>> print s --- !pf.DCSource depth: 0.0 time: 1970-01-01 00:00:00 magnitude: 6.0 strike: 66.0 dip: 33.0 rake: 0.0
'''
for (k, v) in kwargs.items(): self[k] = v
''' Create grid of source model variations.
:returns: :py:class:`SourceGrid` instance.
Example::
>>> from pyrocko import gf >>> base = DCSource() >>> R = gf.Range >>> for s in base.grid(R('
''' return SourceGrid(base=self, variables=variables)
''' Get key to decide about source discretization / GF stack sharing.
When two source models differ only in amplitude and origin time, the discretization and the GF stacking can be done only once for a unit amplitude and a zero origin time and the amplitude and origin times of the seismograms can be applied during post-processing of the synthetic seismogram.
For any derived parameterized source model, this method is called to decide if discretization and stacking of the source should be shared. When two source models return an equal vector of values discretization is shared. ''' self.lon, self.east_shift, self.time, type(self).__name__) + \ self.effective_stf_pre().base_key()
''' Get the scaling factor to be applied during post-processing.
Discretization of the base seismogram is usually done for a unit amplitude, because a common factor can be efficiently multiplied to final seismograms. This eliminates to do repeat the stacking when creating seismograms for a series of source models only differing in amplitude.
This method should return the scaling factor to apply in the post-processing (often this is simply the scalar moment of the source). '''
return 1.0
''' Return the STF applied before stacking of the Green's functions.
This STF is used during discretization of the parameterized source models, i.e. to produce a temporal distribution of point sources.
Handling of the STF before stacking of the GFs is less efficient but allows to use different source time functions for different parts of the source. '''
return self.stf else:
''' Return the STF applied after stacking of the Green's fuctions.
This STF is used in the post-processing of the synthetic seismograms.
Handling of the STF after stacking of the GFs is usually more efficient but is only possible when a common STF is used for all subsources. '''
return self.stf else:
return dict(times=arr(self.time), lat=self.lat, lon=self.lon, north_shifts=arr(self.north_shift), east_shifts=arr(self.east_shift), depths=arr(self.depth))
return self._dparams_base()
lat=self.lat, lon=self.lon, north_shifts=north_shifts, east_shifts=east_shifts, depths=depths)
duration = None if self.stf: duration = self.stf.effective_duration
return model.Event( lat=self.lat, lon=self.lon, north_shift=self.north_shift, east_shift=self.east_shift, time=self.time, name=self.name, depth=self.depth, duration=duration, **kwargs)
points = num.atleast_2d(num.zeros([1, 3]))
points[:, 0] += self.north_shift points[:, 1] += self.east_shift points[:, 2] += self.depth if cs == 'xyz': return points elif cs == 'xy': return points[:, :2] elif cs in ('latlon', 'lonlat'): latlon = ne_to_latlon( self.lat, self.lon, points[:, 0], points[:, 1])
latlon = num.array(latlon).T if cs == 'latlon': return latlon else: return latlon[:, ::-1]
def from_pyrocko_event(cls, ev, **kwargs): if ev.depth is None: raise ConversionError( 'cannot convert event object to source object: ' 'no depth information available')
stf = None if ev.duration is not None: stf = HalfSinusoidSTF(effective_duration=ev.duration)
d = dict( name=ev.name, time=ev.time, lat=ev.lat, lon=ev.lon, north_shift=ev.north_shift, east_shift=ev.east_shift, depth=ev.depth, stf=stf) d.update(kwargs) return cls(**d)
raise NotImplementedError( '%s does not implement get_magnitude()' % self.__class__.__name__)
''' Base class for sources containing a moment magnitude. '''
default=6.0, help='Moment magnitude Mw as in [Hanks and Kanamori, 1979]')
def moment(self):
def moment(self, value): self.magnitude = float(pmt.moment_to_magnitude(value))
return Source.pyrocko_event( self, store, target, magnitude=self.magnitude, **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} if ev.magnitude: d.update(magnitude=ev.magnitude)
d.update(kwargs) return super(SourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
return self.magnitude
Source.T.validate_extra(self, val) val.check_conflicts()
''' Check for parameter conflicts.
To be overloaded in subclasses. Raises :py:exc:`DerivedMagnitudeError` on conflicts. ''' pass
raise DerivedMagnitudeError('no magnitude set')
self.get_magnitude(store, target)))
raise NotImplementedError( '%s does not implement pyrocko_moment_tensor()' % self.__class__.__name__)
try: mt = self.pyrocko_moment_tensor(store, target) magnitude = self.get_magnitude() except (DerivedMagnitudeError, NotImplementedError): mt = None magnitude = None
return Source.pyrocko_event( self, store, target, moment_tensor=mt, magnitude=magnitude, **kwargs)
''' An isotropic explosion point source. '''
optional=True, help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
optional=True, help='volume change of the explosion/implosion or ' 'the contracting/extending magmatic source. [m^3]')
mom = kwargs.pop('moment') if 'magnitude' not in kwargs: kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
return SourceWithDerivedMagnitude.base_key(self) + \ (self.volume_change,)
raise DerivedMagnitudeError( 'magnitude and volume_change are both defined')
return self.magnitude
moment = self.volume_change * \ self.get_moment_to_volume_change_ratio(store, target)
return float(pmt.moment_to_magnitude(abs(moment))) else:
self.check_conflicts()
if self.volume_change is not None: return self.volume_change
elif self.magnitude is not None: moment = float(pmt.magnitude_to_moment(self.magnitude)) return moment / self.get_moment_to_volume_change_ratio( store, target)
else: return 1.0 / self.get_moment_to_volume_change_ratio(store)
if store is None: raise DerivedMagnitudeError( 'need earth model to convert between volume change and ' 'magnitude')
points = num.array( [[self.north_shift, self.east_shift, self.depth]], dtype=num.float)
interpolation = target.interpolation if target else 'multilinear' try: shear_moduli = store.config.get_shear_moduli( self.lat, self.lon, points=points, interpolation=interpolation)[0] except meta.OutOfBounds: raise DerivedMagnitudeError( 'could not get shear modulus at source position')
return float(3. * shear_moduli)
return 1.0
times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time)
amplitudes *= self.get_moment(store, target) * math.sqrt(2. / 3.)
if self.volume_change is not None: if self.volume_change < 0.: amplitudes *= -1
return meta.DiscretizedExplosionSource( m0s=amplitudes, **self._dparams_base_repeated(times))
''' Rectangular or line explosion source. '''
default=0.0, help='strike direction in [deg], measured clockwise from north')
default=90.0, help='dip angle in [deg], measured downward from horizontal')
default=0., help='length of rectangular source area [m]')
default=0., help='width of rectangular source area [m]')
choices=['top', 'top_left', 'top_right', 'center', 'bottom', 'bottom_left', 'bottom_right'], default='center', optional=True, help='Anchor point for positioning the plane, can be: top, center or' 'bottom and also top_left, top_right,bottom_left,' 'bottom_right, center_left and center right')
optional=True, help='horizontal position of rupture nucleation in normalized fault ' 'plane coordinates (-1 = left edge, +1 = right edge)')
optional=True, help='down-dip position of rupture nucleation in normalized fault ' 'plane coordinates (-1 = upper edge, +1 = lower edge)')
default=3500., help='speed of explosion front [m/s]')
return Source.base_key(self) + (self.strike, self.dip, self.length, self.width, self.nucleation_x, self.nucleation_y, self.velocity, self.anchor)
if self.nucleation_x is not None: nucx = self.nucleation_x * 0.5 * self.length else: nucx = None
if self.nucleation_y is not None: nucy = self.nucleation_y * 0.5 * self.width else: nucy = None
stf = self.effective_stf_pre()
points, times, amplitudes, dl, dw, nl, nw = discretize_rect_source( store.config.deltas, store.config.deltat, self.time, self.north_shift, self.east_shift, self.depth, self.strike, self.dip, self.length, self.width, self.anchor, self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy)
amplitudes /= num.sum(amplitudes) amplitudes *= self.get_moment(store, target)
return meta.DiscretizedExplosionSource( lat=self.lat, lon=self.lon, times=times, north_shifts=points[:, 0], east_shifts=points[:, 1], depths=points[:, 2], m0s=amplitudes)
points = outline_rect_source(self.strike, self.dip, self.length, self.width, self.anchor)
points[:, 0] += self.north_shift points[:, 1] += self.east_shift points[:, 2] += self.depth if cs == 'xyz': return points elif cs == 'xy': return points[:, :2] elif cs in ('latlon', 'lonlat'): latlon = ne_to_latlon( self.lat, self.lon, points[:, 0], points[:, 1])
latlon = num.array(latlon).T if cs == 'latlon': return latlon else: return latlon[:, ::-1]
if self.nucleation_x is None: return None, None
coords = from_plane_coords(self.strike, self.dip, self.length, self.width, self.depth, self.nucleation_x, self.nucleation_y, lat=self.lat, lon=self.lon, north_shift=self.north_shift, east_shift=self.east_shift, cs=cs) return coords
''' A double-couple point source. '''
default=0.0, help='strike direction in [deg], measured clockwise from north')
default=90.0, help='dip angle in [deg], measured downward from horizontal')
default=0.0, help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view')
strike=self.strike, dip=self.dip, rake=self.rake)
store.config.deltat, self.time) m6s=mot.m6()[num.newaxis, :] * amplitudes[:, num.newaxis], **self._dparams_base_repeated(times))
strike=self.strike, dip=self.dip, rake=self.rake, scalar_moment=self.moment)
return SourceWithMagnitude.pyrocko_event( self, store, target, moment_tensor=self.pyrocko_moment_tensor(store, target), **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} mt = ev.moment_tensor if mt: (strike, dip, rake), _ = mt.both_strike_dip_rake() d.update( strike=float(strike), dip=float(dip), rake=float(rake), magnitude=float(mt.moment_magnitude()))
d.update(kwargs) return super(DCSource, cls).from_pyrocko_event(ev, **d)
''' A pure CLVD point source. '''
default=0.0, help='azimuth direction of largest dipole, clockwise from north [deg]')
default=90., help='dip direction of largest dipole, downward from horizontal [deg]')
return Source.base_key(self) + (self.azimuth, self.dip)
def m6(self): d2r * (self.dip - 90.), d2r * (self.azimuth - 90.), 0.)
def m6_astuple(self):
factor = self.get_factor() times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time) return meta.DiscretizedMTSource( m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis] / factor, **self._dparams_base_repeated(times))
mt = self.pyrocko_moment_tensor(store, target) return Source.pyrocko_event( self, store, target, moment_tensor=self.pyrocko_moment_tensor(store, target), magnitude=float(mt.moment_magnitude()), **kwargs)
''' Volumetric linear vector dipole source.
This source is a parameterization for a restricted moment tensor point source, useful to represent dyke or sill like inflation or deflation sources. The restriction is such that the moment tensor is rotational symmetric. It can be represented by a superposition of a linear vector dipole (here we use a CLVD for convenience) and an isotropic component. The restricted moment tensor has 4 degrees of freedom: 2 independent eigenvalues and 2 rotation angles orienting the the symmetry axis.
In this parameterization, the isotropic component is controlled by ``volume_change``. To define the moment tensor, it must be converted to the scalar moment of the the MT's isotropic component. For the conversion, the shear modulus at the source's position must be known. This value is extracted from the earth model defined in the GF store in use.
The CLVD part by controlled by its scalar moment :math:`M_0`: ``clvd_moment``. The sign of ``clvd_moment`` is used to switch between a positiv or negativ CLVD (the sign of the largest eigenvalue). '''
default=0.0, help='azimuth direction of symmetry axis, clockwise from north [deg].')
default=90., help='dip direction of symmetry axis, downward from horizontal [deg].')
default=0., help='volume change of the inflation/deflation [m^3].')
default=0., help='scalar moment :math:`M_0` of the CLVD component [Nm]. The sign ' 'controls the sign of the CLVD (the sign of its largest ' 'eigenvalue).')
if store is None or target is None: raise DerivedMagnitudeError( 'need earth model to convert between volume change and ' 'magnitude')
points = num.array( [[self.north_shift, self.east_shift, self.depth]], dtype=num.float)
try: shear_moduli = store.config.get_shear_moduli( self.lat, self.lon, points=points, interpolation=target.interpolation)[0] except meta.OutOfBounds: raise DerivedMagnitudeError( 'could not get shear modulus at source position')
return float(3. * shear_moduli)
return Source.base_key(self) + \ (self.azimuth, self.dip, self.volume_change, self.clvd_moment)
mt = self.pyrocko_moment_tensor(store, target) return float(pmt.moment_to_magnitude(mt.moment))
a = math.sqrt(4. / 3.) * self.clvd_moment m_clvd = pmt.symmat6(-0.5 * a, -0.5 * a, a, 0., 0., 0.)
rotmat1 = pmt.euler_to_matrix( d2r * (self.dip - 90.), d2r * (self.azimuth - 90.), 0.) m_clvd = rotmat1.T * m_clvd * rotmat1
m_iso = self.volume_change * \ self.get_moment_to_volume_change_ratio(store, target)
m_iso = pmt.symmat6(m_iso, m_iso, m_iso, 0., 0., 0.,) * math.sqrt(2./3)
m = pmt.to6(m_clvd) + pmt.to6(m_iso) return m
return float(pmt.magnitude_to_moment( self.get_magnitude(store, target)))
m6 = self.get_m6(store, target) return tuple(m6.tolist())
times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time)
m6 = self.get_m6(store, target) m6 *= amplitudes / self.get_factor()
return meta.DiscretizedMTSource( m6s=m6[num.newaxis, :], **self._dparams_base_repeated(times))
m6_astuple = self.get_m6_astuple(store, target) return pmt.MomentTensor(m=pmt.symmat6(*m6_astuple))
''' A moment tensor point source. '''
default=1., help='north-north component of moment tensor in [Nm]')
default=1., help='east-east component of moment tensor in [Nm]')
default=1., help='down-down component of moment tensor in [Nm]')
default=0., help='north-east component of moment tensor in [Nm]')
default=0., help='north-down component of moment tensor in [Nm]')
default=0., help='east-down component of moment tensor in [Nm]')
if 'm6' in kwargs: for (k, v) in zip('mnn mee mdd mne mnd med'.split(), kwargs.pop('m6')): kwargs[k] = float(v)
Source.__init__(self, **kwargs)
def m6(self): return num.array(self.m6_astuple)
def m6_astuple(self): return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)
def m6(self, value): self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value
return Source.base_key(self) + self.m6_astuple
times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time) return meta.DiscretizedMTSource( m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis], **self._dparams_base_repeated(times))
m6 = self.m6 return pmt.moment_to_magnitude( math.sqrt(num.sum(m6[0:3]**2) + 2.0 * num.sum(m6[3:6]**2)) / math.sqrt(2.))
return pmt.MomentTensor(m=pmt.symmat6(*self.m6_astuple))
mt = self.pyrocko_moment_tensor(store, target) return Source.pyrocko_event( self, store, target, moment_tensor=self.pyrocko_moment_tensor(store, target), magnitude=float(mt.moment_magnitude()), **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} mt = ev.moment_tensor if mt: d.update(m6=tuple(map(float, mt.m6()))) else: if ev.magnitude is not None: mom = pmt.magnitude_to_moment(ev.magnitude) v = math.sqrt(2./3.) * mom d.update(m6=(v, v, v, 0., 0., 0.))
d.update(kwargs) return super(MTSource, cls).from_pyrocko_event(ev, **d)
'center': (0.0, 0.0), 'center_left': (-1.0, 0.0), 'center_right': (1.0, 0.0), 'top': (0.0, -1.0), 'top_left': (-1.0, -1.0), 'top_right': (1.0, -1.0), 'bottom': (0.0, 1.0), 'bottom_left': (-1.0, 1.0), 'bottom_right': (1.0, 1.0)}
''' Classical Haskell source model modified for bilateral rupture. '''
optional=True, help='moment magnitude Mw as in [Hanks and Kanamori, 1979]')
default=0.0, help='strike direction in [deg], measured clockwise from north')
default=90.0, help='dip angle in [deg], measured downward from horizontal')
default=0.0, help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view')
default=0., help='length of rectangular source area [m]')
default=0., help='width of rectangular source area [m]')
choices=['top', 'top_left', 'top_right', 'center', 'bottom', 'bottom_left', 'bottom_right'], default='center', optional=True, help='Anchor point for positioning the plane, can be: top, center or' 'bottom and also top_left, top_right,bottom_left,' 'bottom_right, center_left and center right')
optional=True, help='horizontal position of rupture nucleation in normalized fault ' 'plane coordinates (-1 = left edge, +1 = right edge)')
optional=True, help='down-dip position of rupture nucleation in normalized fault ' 'plane coordinates (-1 = upper edge, +1 = lower edge)')
default=3500., help='speed of rupture front [m/s]')
optional=True, help='Slip on the rectangular source area [m]')
optional=True, default=1, help='Sub-source decimation factor, a larger decimation will' ' make the result inaccurate but shorten the necessary' ' computation time (use for testing puposes only).')
mom = kwargs.pop('moment') if 'magnitude' not in kwargs: kwargs['magnitude'] = float(pmt.moment_to_magnitude(mom))
self.magnitude, self.slip, self.strike, self.dip, self.rake, self.length, self.width, self.nucleation_x, self.nucleation_y, self.velocity, self.decimation_factor, self.anchor)
if self.magnitude is not None and self.slip is not None: raise DerivedMagnitudeError( 'magnitude and slip are both defined')
self.check_conflicts() if self.magnitude is not None: return self.magnitude
elif self.slip is not None: if None in (store, target): raise DerivedMagnitudeError( 'magnitude for a rectangular source with slip defined ' 'can only be derived when earth model and target ' 'interpolation method are available')
amplitudes = self._discretize(store, target)[2] return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
else: return float(pmt.moment_to_magnitude(1.0))
nucx = self.nucleation_x * 0.5 * self.length else:
nucy = self.nucleation_y * 0.5 * self.width else:
store.config.deltas, store.config.deltat, self.time, self.north_shift, self.east_shift, self.depth, self.strike, self.dip, self.length, self.width, self.anchor, self.velocity, stf=stf, nucleation_x=nucx, nucleation_y=nucy, decimation_factor=self.decimation_factor)
else: interpolation = 'nearest_neighbor' logger.warn( 'no target information available, will use ' '"nearest_neighbor" interpolation when extracting shear ' 'modulus from earth model')
self.lat, self.lon, points=points, interpolation=interpolation)
else: # normalization to retain total moment amplitudes /= num.sum(amplitudes) amplitudes *= self.get_moment(store, target)
strike=self.strike, dip=self.dip, rake=self.rake)
lat=self.lat, lon=self.lon, times=times, north_shifts=points[:, 0], east_shifts=points[:, 1], depths=points[:, 2], m6s=m6s)
self.width, self.anchor)
return points elif cs in ('latlon', 'lonlat'): latlon = ne_to_latlon( self.lat, self.lon, points[:, 0], points[:, 1])
latlon = num.array(latlon).T if cs == 'latlon': return latlon else: return latlon[:, ::-1]
if self.nucleation_x is None: return None, None
coords = from_plane_coords(self.strike, self.dip, self.length, self.width, self.depth, self.nucleation_x, self.nucleation_y, lat=self.lat, lon=self.lon, north_shift=self.north_shift, east_shift=self.east_shift, cs=cs) return coords
return pmt.MomentTensor( strike=self.strike, dip=self.dip, rake=self.rake, scalar_moment=self.get_moment(store, target))
return SourceWithDerivedMagnitude.pyrocko_event( self, store, target, **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} mt = ev.moment_tensor if mt: (strike, dip, rake), _ = mt.both_strike_dip_rake() d.update( strike=float(strike), dip=float(dip), rake=float(rake), magnitude=float(mt.moment_magnitude()))
d.update(kwargs) return super(RectangularSource, cls).from_pyrocko_event(ev, **d)
''' Two double-couple point sources separated in space and time. Moment share between the sub-sources is controlled by the parameter mix. The position of the subsources is dependent on the moment distribution between the two sources. Depth, east and north shift are given for the centroid between the two double-couples. The subsources will positioned according to their moment shares around this centroid position. This is done according to their delta parameters, which are therefore in relation to that centroid. Note that depth of the subsources therefore can be depth+/-delta_depth. For shallow earthquakes therefore the depth has to be chosen deeper to avoid sampling above surface. '''
default=0.0, help='strike direction in [deg], measured clockwise from north')
default=90.0, help='dip angle in [deg], measured downward from horizontal')
default=0.0, help='azimuth to second double-couple [deg], ' 'measured at first, clockwise from north')
default=0.0, help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view')
default=0.0, help='strike direction in [deg], measured clockwise from north')
default=90.0, help='dip angle in [deg], measured downward from horizontal')
default=0.0, help='rake angle in [deg], ' 'measured counter-clockwise from right-horizontal ' 'in on-plane view')
default=0.0, help='separation of double-couples in time (t2-t1) [s]')
default=0.0, help='difference in depth (z2-z1) [m]')
default=0.0, help='distance between the two double-couples [m]')
default=0.5, help='how to distribute the moment to the two doublecouples ' 'mix=0 -> m1=1 and m2=0; mix=1 -> m1=0, m2=1')
optional=True, help='Source time function of subsource 1 ' '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
optional=True, help='Source time function of subsource 2 ' '(if given, overrides STF from attribute :py:gattr:`Source.stf`)')
return ( self.time, self.depth, self.lat, self.north_shift, self.lon, self.east_shift, type(self).__name__) + \ self.effective_stf1_pre().base_key() + \ self.effective_stf2_pre().base_key() + ( self.strike1, self.dip1, self.rake1, self.strike2, self.dip2, self.rake2, self.delta_time, self.delta_depth, self.azimuth, self.distance, self.mix)
return self.moment
return self.stf1 or self.stf or g_unit_pulse
return self.stf2 or self.stf or g_unit_pulse
return g_unit_pulse
a1 = 1.0 - self.mix a2 = self.mix delta_north = math.cos(self.azimuth * d2r) * self.distance delta_east = math.sin(self.azimuth * d2r) * self.distance
dc1 = DCSource( lat=self.lat, lon=self.lon, time=self.time - self.delta_time * a1, north_shift=self.north_shift - delta_north * a1, east_shift=self.east_shift - delta_east * a1, depth=self.depth - self.delta_depth * a1, moment=self.moment * a1, strike=self.strike1, dip=self.dip1, rake=self.rake1, stf=self.stf1 or self.stf)
dc2 = DCSource( lat=self.lat, lon=self.lon, time=self.time + self.delta_time * a2, north_shift=self.north_shift + delta_north * a2, east_shift=self.east_shift + delta_east * a2, depth=self.depth + self.delta_depth * a2, moment=self.moment * a2, strike=self.strike2, dip=self.dip2, rake=self.rake2, stf=self.stf2 or self.stf)
return [dc1, dc2]
a1 = 1.0 - self.mix a2 = self.mix mot1 = pmt.MomentTensor(strike=self.strike1, dip=self.dip1, rake=self.rake1, scalar_moment=a1) mot2 = pmt.MomentTensor(strike=self.strike2, dip=self.dip2, rake=self.rake2, scalar_moment=a2)
delta_north = math.cos(self.azimuth * d2r) * self.distance delta_east = math.sin(self.azimuth * d2r) * self.distance
times1, amplitudes1 = self.effective_stf1_pre().discretize_t( store.config.deltat, self.time - self.delta_time * a1)
times2, amplitudes2 = self.effective_stf2_pre().discretize_t( store.config.deltat, self.time + self.delta_time * a2)
nt1 = times1.size nt2 = times2.size
ds = meta.DiscretizedMTSource( lat=self.lat, lon=self.lon, times=num.concatenate((times1, times2)), north_shifts=num.concatenate(( num.repeat(self.north_shift - delta_north * a1, nt1), num.repeat(self.north_shift + delta_north * a2, nt2))), east_shifts=num.concatenate(( num.repeat(self.east_shift - delta_east * a1, nt1), num.repeat(self.east_shift + delta_east * a2, nt2))), depths=num.concatenate(( num.repeat(self.depth - self.delta_depth * a1, nt1), num.repeat(self.depth + self.delta_depth * a2, nt2))), m6s=num.vstack(( mot1.m6()[num.newaxis, :] * amplitudes1[:, num.newaxis], mot2.m6()[num.newaxis, :] * amplitudes2[:, num.newaxis])))
return ds
rake=self.rake1, scalar_moment=a1 * self.moment) rake=self.rake2, scalar_moment=a2 * self.moment)
return SourceWithMagnitude.pyrocko_event( self, store, target, moment_tensor=self.pyrocko_moment_tensor(store, target), **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} mt = ev.moment_tensor if mt: (strike, dip, rake), _ = mt.both_strike_dip_rake() d.update( strike1=float(strike), dip1=float(dip), rake1=float(rake), strike2=float(strike), dip2=float(dip), rake2=float(rake), mix=0.0, magnitude=float(mt.moment_magnitude()))
d.update(kwargs) source = super(DoubleDCSource, cls).from_pyrocko_event(ev, **d) source.stf1 = source.stf source.stf2 = HalfSinusoidSTF(effective_duration=0.) source.stf = None return source
''' A ring fault with vertical doublecouples. '''
default=1.0, help='diameter of the ring in [m]')
default=1.0, help='inside of the ring moves up (+1) or down (-1)')
default=0.0, help='strike direction of the ring plane, clockwise from north,' ' in [deg]')
default=0.0, help='dip angle of the ring plane from horizontal in [deg]')
default=360, help='number of point sources to use')
return Source.base_key(self) + ( self.strike, self.dip, self.diameter, self.npointsources)
return self.sign * self.moment
n = self.npointsources phi = num.linspace(0, 2.0 * num.pi, n, endpoint=False)
points = num.zeros((n, 3)) points[:, 0] = num.cos(phi) * 0.5 * self.diameter points[:, 1] = num.sin(phi) * 0.5 * self.diameter
rotmat = num.array(pmt.euler_to_matrix( self.dip * d2r, self.strike * d2r, 0.0)) points = num.dot(rotmat.T, points.T).T # !!! ?
points[:, 0] += self.north_shift points[:, 1] += self.east_shift points[:, 2] += self.depth
m = num.array(pmt.MomentTensor(strike=90., dip=90., rake=-90., scalar_moment=1.0 / n).m())
rotmats = num.transpose( [[num.cos(phi), num.sin(phi), num.zeros(n)], [-num.sin(phi), num.cos(phi), num.zeros(n)], [num.zeros(n), num.zeros(n), num.ones(n)]], (2, 0, 1))
ms = num.zeros((n, 3, 3)) for i in range(n): mtemp = num.dot(rotmats[i].T, num.dot(m, rotmats[i])) ms[i, :, :] = num.dot(rotmat.T, num.dot(mtemp, rotmat))
m6s = num.vstack((ms[:, 0, 0], ms[:, 1, 1], ms[:, 2, 2], ms[:, 0, 1], ms[:, 0, 2], ms[:, 1, 2])).T
times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time)
nt = times.size
return meta.DiscretizedMTSource( times=num.tile(times, n), lat=self.lat, lon=self.lon, north_shifts=num.repeat(points[:, 0], nt), east_shifts=num.repeat(points[:, 1], nt), depths=num.repeat(points[:, 2], nt), m6s=num.repeat(m6s, nt, axis=0) * num.tile( amplitudes, n)[:, num.newaxis])
'''Composite source model.'''
if not subsources: raise BadRequest( 'need at least one sub-source to create a CombiSource object.')
lats = num.array( [subsource.lat for subsource in subsources], dtype=num.float) lons = num.array( [subsource.lon for subsource in subsources], dtype=num.float)
lat, lon = lats[0], lons[0] if not num.all(lats == lat) and num.all(lons == lon): subsources = [s.clone() for s in subsources] for subsource in subsources[1:]: subsource.set_origin(lat, lon)
depth = float(num.mean([p.depth for p in subsources])) time = float(num.mean([p.time for p in subsources])) north_shift = float(num.mean([p.north_shift for p in subsources])) east_shift = float(num.mean([p.east_shift for p in subsources])) kwargs.update( time=time, lat=float(lat), lon=float(lon), north_shift=north_shift, east_shift=east_shift, depth=depth)
Source.__init__(self, subsources=subsources, **kwargs)
return 1.0
dsources = [] for sf in self.subsources: ds = sf.discretize_basesource(store, target) ds.m6s *= sf.get_factor() dsources.append(ds)
return meta.DiscretizedMTSource.combine(dsources)
''' A single force point source. '''
default=0., help='northward component of single force [N]')
default=0., help='eastward component of single force [N]')
default=0., help='downward component of single force [N]')
Source.__init__(self, **kwargs)
return Source.base_key(self) + (self.fn, self.fe, self.fd)
return 1.0
times, amplitudes = self.effective_stf_pre().discretize_t( store.config.deltat, self.time) forces = num.array([[self.fn, self.fe, self.fd]], dtype=num.float) forces *= amplitudes[:, num.newaxis] return meta.DiscretizedSFSource(forces=forces, **self._dparams_base_repeated(times))
return Source.pyrocko_event( self, store, target, **kwargs)
def from_pyrocko_event(cls, ev, **kwargs): d = {} d.update(kwargs) return super(SFSource, cls).from_pyrocko_event(ev, **d)
''' Excess pore pressure point source.
For poro-elastic initial value problem where an excess pore pressure is brought into a small source volume. '''
default=1.0, help='initial excess pore pressure in [Pa]')
return Source.base_key(self)
return self.pp
return meta.DiscretizedPorePressureSource(pp=arr(1.0), **self._dparams_base())
''' Excess pore pressure line source.
The line source is centered at (north_shift, east_shift, depth). '''
default=1.0, help='initial excess pore pressure in [Pa]')
default=0.0, help='length of the line source [m]')
default=0.0, help='azimuth direction, clockwise from north [deg]')
default=90., help='dip direction, downward from horizontal [deg]')
return Source.base_key(self) + (self.azimuth, self.dip, self.length)
return self.pp
n = 2 * int(math.ceil(self.length / num.min(store.config.deltas))) + 1
a = num.linspace(-0.5 * self.length, 0.5 * self.length, n)
sa = math.sin(self.azimuth * d2r) ca = math.cos(self.azimuth * d2r) sd = math.sin(self.dip * d2r) cd = math.cos(self.dip * d2r)
points = num.zeros((n, 3)) points[:, 0] = self.north_shift + a * ca * cd points[:, 1] = self.east_shift + a * sa * cd points[:, 2] = self.depth + a * sd
return meta.DiscretizedPorePressureSource( times=util.num_full(n, self.time), lat=self.lat, lon=self.lon, north_shifts=points[:, 0], east_shifts=points[:, 1], depths=points[:, 2], pp=num.ones(n) / n)
''' Synthetic seismogram computation request.
::
Request(**kwargs) Request(sources, targets, **kwargs) '''
Source.T(), help='list of sources for which to produce synthetics.')
Target.T(), help='list of targets for which to produce synthetics.')
def args2kwargs(cls, args): raise BadRequest('invalid arguments')
else:
targets = [targets]
def targets_dynamic(self):
def targets_static(self):
def has_dynamic(self):
def has_statics(self):
''' Resonse object to a synthetic seismogram computation request. '''
''' Return a list of requested :class:`~pyrocko.trace.Trace` instances. '''
continue
''' Return a list of requested :class:`~kite.scenes` instances. ''' kite_scenes = [] for results in self.results_list: for result in results: if isinstance(result, meta.KiteSceneResult): sc = result.get_scene() kite_scenes.append(sc)
return kite_scenes
''' Return a list of requested :class:`~pyrocko.gf.meta.StaticResult` instances. ''' statics = [] for results in self.results_list: for result in results: if not isinstance(result, meta.StaticResult): continue statics.append(result)
return statics
''' Generator function to iterate over results of request.
Yields associated :py:class:`Source`, :class:`~pyrocko.gf.targets.Target`, :class:`~pyrocko.trace.Trace` instances in each iteration. '''
for isource, source in enumerate(self.request.sources): for itarget, target in enumerate(self.request.targets): result = self.results_list[isource][itarget] if get == 'pyrocko_traces': yield source, target, result.trace.pyrocko_trace() elif get == 'results': yield source, target, result
''' Open *snuffler* with requested traces. '''
trace.snuffle(self.pyrocko_traces(), **kwargs)
''' Base class for synthetic seismogram calculators. '''
''' Get list of available GF store IDs '''
return []
else:
data = num.diff(data)
n, e = self.components sa, ca, _, _ = target.get_sin_cos_factors()
comps = [] if nonzero(ca): comps.append(n)
if nonzero(sa): comps.append(e)
return tuple(comps)
n, e = self.components sa, ca, _, _ = target.get_sin_cos_factors()
if nonzero(ca): data = base_seismogram[n].data * ca else: data = 0.0
if nonzero(sa): data = data + base_seismogram[e].data * sa
return data
return (self.c, )
return base_seismogram[self.c].data.copy()
(los_fac[:, 0] * -base_statics['displacement.d'] + los_fac[:, 1] * base_statics['displacement.e'] + los_fac[:, 2] * base_statics['displacement.n'])
'displacement': [VectorRule('displacement')], 'velocity': [ VectorRule('velocity'), VectorRule('displacement', differentiate=1)], 'acceleration': [ VectorRule('acceleration'), VectorRule('velocity', differentiate=1), VectorRule('displacement', differentiate=2)], 'pore_pressure': [ScalarRule('pore_pressure')], 'vertical_tilt': [HorizontalVectorRule('vertical_tilt')], 'darcy_velocity': [VectorRule('darcy_velocity')], }
'displacement': [StaticDisplacement()] }
source, store_targets, components, dsource_cache, nthreads)
e = SeismosizerError( 'Seismosizer failed with return code %i\n%s' % ( tr.err, str( OutOfBoundsContext( source=source, target=store_targets[iseis], distance=source.distance_to( store_targets[iseis]), components=components)))) raise e
seismogram, source, target) except SeismosizerError as e: result = e
dsource_cache = {}
for w in work: _, _, isources, itargets = w
sources = [psources[isource] for isource in isources] targets = [ptargets[itarget] for itarget in itargets]
components = set() for target in targets: rule = engine.get_rule(sources[0], target) components.update(rule.required_components(target))
for isource, source in zip(isources, sources): for itarget, target in zip(itargets, targets):
try: base_seismogram, tcounters = engine.base_seismogram( source, target, components, dsource_cache, nthreads) except meta.OutOfBounds as e: e.context = OutOfBoundsContext( source=sources[0], target=targets[0], distance=sources[0].distance_to(targets[0]), components=components) raise
n_records_stacked = 0 t_optimize = 0.0 t_stack = 0.0
for _, tr in base_seismogram.items(): n_records_stacked += tr.n_records_stacked t_optimize += tr.t_optimize t_stack += tr.t_stack
try: result = engine._post_process_dynamic( base_seismogram, source, target) result.n_records_stacked = n_records_stacked result.n_shared_stacking = len(sources) *\ len(targets) result.t_optimize = t_optimize result.t_stack = t_stack except SeismosizerError as e: result = e
tcounters.append(xtime()) yield (isource, itarget, result), tcounters
.required_components(target)
source, target, components, nthreads) except meta.OutOfBounds as e: e.context = OutOfBoundsContext( source=sources[0], target=targets[0], distance=float('nan'), components=components) raise base_statics, source, target)
''' Offline synthetic seismogram calculator.
:param use_env: if ``True``, fill :py:attr:`store_superdirs` and :py:attr:`store_dirs` with paths set in environment variables GF_STORE_SUPERDIRS and GF_STORE_DIRS. :param use_config: if ``True``, fill :py:attr:`store_superdirs` and :py:attr:`store_dirs` with paths set in the user's config file.
The config file can be found at :file:`~/.pyrocko/config.pf`
.. code-block :: python
gf_store_dirs: ['/home/pyrocko/gf_stores/ak135/'] gf_store_superdirs: ['/home/pyrocko/gf_stores/'] '''
String.T(), help='directories which are searched for Green\'s function stores')
String.T(), help='additional individual Green\'s function store directories')
optional=True, help='default store ID to be used when a request does not provide ' 'one')
env_store_superdirs = os.environ.get('GF_STORE_SUPERDIRS', '') env_store_dirs = os.environ.get('GF_STORE_DIRS', '') if env_store_superdirs: self.store_superdirs.extend(env_store_superdirs.split(':'))
if env_store_dirs: self.store_dirs.extend(env_store_dirs.split(':'))
raise TypeError("{} of {} is not of type list".format( sdir, self.__class__.__name__))
all(os.path.isfile(pjoin(store_dir, x)) for x in ('index', 'traces', 'config'))
logger.warning('store_superdir not available: %s' % d) continue
store_dirs.add(os.path.realpath(store_dir))
else: if store_dir != self._id_to_store_dir[store_id]: raise DuplicateStoreId( 'GF store ID %s is used in (at least) two ' 'different stores. Locations are: %s and %s' % (store_id, self._id_to_store_dir[store_id], store_dir))
''' Lookup directory given a GF store ID. '''
raise NoSuchStore(store_id, self.iter_store_dirs())
''' Get list of available store IDs. '''
if self._effective_default_store_id is None: if self.default_store_id is None: store_ids = self.get_store_ids() if len(store_ids) == 1: self._effective_default_store_id = self.get_store_ids()[0] else: raise NoDefaultStoreSet() else: self._effective_default_store_id = self.default_store_id
return self._effective_default_store_id
''' Get a store from the engine.
:param store_id: identifier of the store (optional) :returns: :py:class:`~pyrocko.gf.store.Store` object
If no ``store_id`` is provided the store associated with the :py:gattr:`default_store_id` is returned. Raises :py:exc:`NoDefaultStoreSet` if :py:gattr:`default_store_id` is undefined. '''
store_id = self.effective_default_store_id()
store = self.get_store(store_id) return store.config
store = self.get_store(store_id) return store.get_extra(key)
''' Close and remove ids from cashed stores. ''' store_ids = [] for store_id, store_ in self._open_stores.items(): store_.close() store_ids.append(store_id)
for store_id in store_ids: self._open_stores.pop(store_id)
except KeyError: pass
raise BadRequest( 'no rule to calculate "%s" with GFs from store "%s" ' 'for source model "%s"' % ( target.effective_quantity(), target.store_id, source.__class__.__name__))
nthreads=0):
logging.warning('Targets have different interpolation schemes!' ' Choosing %s for all targets.' % target.interpolation)
(t.tmin for t in targets), dtype=num.float, count=len(targets)) (t.tmax for t in targets), dtype=num.float, count=len(targets))
source, store_, dsource_cache, target)
deltat = 1. / target.sample_rate else:
base_source, receivers, components, deltat=deltat, itmin=itmin, nsamples=nsamples, interpolation=target.interpolation, optimization=target.optimization, nthreads=nthreads)
nthreads):
tcounters = [xtime()]
store_ = self.get_store(target.store_id) receiver = target.receiver(store_)
if target.tmin and target.tmax is not None: rate = store_.config.sample_rate itmin = int(num.floor(target.tmin * rate)) itmax = int(num.ceil(target.tmax * rate)) nsamples = itmax - itmin + 1 else: itmin = None nsamples = None
tcounters.append(xtime()) base_source = self._cached_discretize_basesource( source, store_, dsource_cache, target)
tcounters.append(xtime())
if target.sample_rate is not None: deltat = 1. / target.sample_rate else: deltat = None
base_seismogram = store_.seismogram( base_source, receiver, components, deltat=deltat, itmin=itmin, nsamples=nsamples, interpolation=target.interpolation, optimization=target.optimization, nthreads=nthreads)
tcounters.append(xtime())
base_seismogram = store.make_same_span(base_seismogram)
tcounters.append(xtime())
return base_seismogram, tcounters
else: itsnapshot = None
base_source, target, itsnapshot, components, target.interpolation, nthreads)
deltat, 0.0)
# repeat end point to prevent boundary effects
codes=target.codes, data=data[:-amplitudes.size], deltat=deltat, tmin=tmin)
for v in data.values(): v *= factor
''' Process a request.
::
process(**kwargs) process(request, **kwargs) process(sources, targets, **kwargs)
The request can be given a a :py:class:`Request` object, or such an object is created using ``Request(**kwargs)`` for convenience.
:returns: :py:class:`Response` object '''
raise BadRequest('invalid arguments')
kwargs['request'] = args[0]
nthreads = nprocs
# make sure stores are open before fork()
enumerate(request.sources)) enumerate(request.targets))
# Processing dynamic targets through # parimap(process_subrequest_dynamic)
else: _process_dynamic = process_dynamic
(i, nsub, [source_index[source] for source in m[k][0]], [target_index[target] for target in m[k][1] if not isinstance(target, StaticTarget)]) for (i, k) in enumerate(skeys)]
work_dynamic, request.sources, request.targets, self, nthreads):
status_callback(isub, nsub)
# Processing static targets through process_static (i, nsub, [source_index[source] for source in m[k][0]], [target_index[target] for target in m[k][1] if isinstance(target, StaticTarget)]) for (i, k) in enumerate(skeys)]
work_static, request.sources, request.targets, self, nthreads=nthreads):
status_callback(isub, nsub)
status_callback(nsub, nsub)
s.t_perc_discretize_source, s.t_perc_make_base_seismogram, s.t_perc_make_same_span, s.t_perc_post_process) = perc_dyn else:
s.t_perc_static_discretize_basesource, s.t_perc_static_sum_statics, s.t_perc_static_post_process) = perc_static
(rs1.ru_utime + rs1.ru_stime + rc1.ru_utime + rc1.ru_stime) - (rs0.ru_utime + rs0.ru_stime + rc0.ru_utime + rc0.ru_stime)) (rs1.ru_inblock + rc1.ru_inblock) - (rs0.ru_inblock + rc0.ru_inblock))
request=request, results_list=results_list, stats=s)
''' Client for remote synthetic seismogram calculator. '''
if request is None: request = Request(**kwargs)
return ws.seismosizer(url=self.url, site=self.site, request=request)
global g_engine if g_engine is None: g_engine = LocalEngine(use_env=True, use_config=True)
for d in store_superdirs: if d not in g_engine.store_superdirs: g_engine.store_superdirs.append(d)
return g_engine
return num.fromiter((getattr(s, k) for s in self), dtype=num.float)
raise NotImplementedError( 'this method should be implemented in subclass')
raise NotImplementedError( 'this method should be implemented in subclass')
self.sources.append(s)
return iter(self.sources)
return len(self.sources)
n = 1 for (k, v) in self.make_coords(self.base): n *= len(list(v))
return n
for items in permudef(self.make_coords(self.base)): s = self.base.clone(**{k: v for (k, v) in items}) s.regularize() yield s
ks = list(self.variables.keys()) for k in self.order + list(self.base.keys()): if k in ks: yield k ks.remove(k) if ks: raise Exception('Invalid parameter "%s" for source type "%s"' % (ks[0], self.base.__class__.__name__))
return [(param, self.variables[param].make(base=base[param])) for param in self.ordered_params()]
Source, SourceWithMagnitude, SourceWithDerivedMagnitude, ExplosionSource, RectangularExplosionSource, DCSource, CLVDSource, VLVDSource, MTSource, RectangularSource, DoubleDCSource, RingfaultSource, CombiSource, SFSource, PorePressurePointSource, PorePressureLineSource, ]
STF, BoxcarSTF, TriangularSTF, HalfSinusoidSTF, ResonatorSTF, ]
SeismosizerError BadRequest NoSuchStore DerivedMagnitudeError STFMode '''.split() + [S.__name__ for S in source_classes + stf_classes] + ''' Request ProcessingStats Response Engine LocalEngine RemoteEngine source_classes get_engine Range SourceGroup SourceList SourceGrid '''.split() |