1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import re
8import fnmatch
9import logging
11import numpy as num
12from scipy.interpolate import interp1d
14from pyrocko.guts import (Object, SObject, String, StringChoice,
15 StringPattern, Unicode, Float, Bool, Int, TBase,
16 List, ValidationError, Timestamp, Tuple, Dict)
17from pyrocko.guts import dump, load # noqa
18from pyrocko.guts_array import literal, Array
19from pyrocko.model import Location, gnss
21from pyrocko import cake, orthodrome, spit, moment_tensor, trace
22from pyrocko.util import num_full
24from .error import StoreError
26try:
27 new_str = unicode
28except NameError:
29 new_str = str
31guts_prefix = 'pf'
33d2r = math.pi / 180.
34r2d = 1.0 / d2r
35km = 1000.
36vicinity_eps = 1e-5
38logger = logging.getLogger('pyrocko.gf.meta')
41def fmt_choices(cls):
42 return 'choices: %s' % ', '.join(
43 "``'%s'``" % choice for choice in cls.choices)
46class InterpolationMethod(StringChoice):
47 choices = ['nearest_neighbor', 'multilinear']
50class SeismosizerTrace(Object):
52 codes = Tuple.T(
53 4, String.T(),
54 default=('', 'STA', '', 'Z'),
55 help='network, station, location and channel codes')
57 data = Array.T(
58 shape=(None,),
59 dtype=num.float32,
60 serialize_as='base64',
61 serialize_dtype=num.dtype('<f4'),
62 help='numpy array with data samples')
64 deltat = Float.T(
65 default=1.0,
66 help='sampling interval [s]')
68 tmin = Timestamp.T(
69 default=Timestamp.D('1970-01-01 00:00:00'),
70 help='time of first sample as a system timestamp [s]')
72 def pyrocko_trace(self):
73 c = self.codes
74 return trace.Trace(c[0], c[1], c[2], c[3],
75 ydata=self.data,
76 deltat=self.deltat,
77 tmin=self.tmin)
79 @classmethod
80 def from_pyrocko_trace(cls, tr, **kwargs):
81 d = dict(
82 codes=tr.nslc_id,
83 tmin=tr.tmin,
84 deltat=tr.deltat,
85 data=num.asarray(tr.get_ydata(), dtype=num.float32))
86 d.update(kwargs)
87 return cls(**d)
90class SeismosizerResult(Object):
91 n_records_stacked = Int.T(optional=True, default=1)
92 t_stack = Float.T(optional=True, default=0.)
95class Result(SeismosizerResult):
96 trace = SeismosizerTrace.T(optional=True)
97 n_shared_stacking = Int.T(optional=True, default=1)
98 t_optimize = Float.T(optional=True, default=0.)
101class StaticResult(SeismosizerResult):
102 result = Dict.T(
103 String.T(),
104 Array.T(shape=(None,), dtype=float, serialize_as='base64'))
107class GNSSCampaignResult(StaticResult):
108 campaign = gnss.GNSSCampaign.T(
109 optional=True)
112class SatelliteResult(StaticResult):
114 theta = Array.T(
115 optional=True,
116 shape=(None,), dtype=float, serialize_as='base64')
118 phi = Array.T(
119 optional=True,
120 shape=(None,), dtype=float, serialize_as='base64')
123class KiteSceneResult(SatelliteResult):
125 shape = Tuple.T()
127 def get_scene(self, component='displacement.los'):
128 try:
129 from kite import Scene
130 except ImportError:
131 raise ImportError('Kite not installed')
133 def reshape(arr):
134 return arr.reshape(self.shape)
136 displacement = self.result[component]
138 displacement, theta, phi = map(
139 reshape, (displacement, self.theta, self.phi))
141 sc = Scene(
142 displacement=displacement,
143 theta=theta, phi=phi,
144 config=self.config)
146 return sc
149class ComponentSchemeDescription(Object):
150 name = String.T()
151 source_terms = List.T(String.T())
152 ncomponents = Int.T()
153 default_stored_quantity = String.T(optional=True)
154 provided_components = List.T(String.T())
157component_scheme_descriptions = [
158 ComponentSchemeDescription(
159 name='scalar1',
160 source_terms=['m0'],
161 ncomponents=1,
162 default_stored_quantity=None,
163 provided_components=[
164 'scalar']),
165 ComponentSchemeDescription(
166 name='elastic2',
167 source_terms=['m0'],
168 ncomponents=2,
169 default_stored_quantity='displacement',
170 provided_components=[
171 'n', 'e', 'd']),
172 ComponentSchemeDescription(
173 name='elastic8',
174 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
175 ncomponents=8,
176 default_stored_quantity='displacement',
177 provided_components=[
178 'n', 'e', 'd']),
179 ComponentSchemeDescription(
180 name='elastic10',
181 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
182 ncomponents=10,
183 default_stored_quantity='displacement',
184 provided_components=[
185 'n', 'e', 'd']),
186 ComponentSchemeDescription(
187 name='elastic18',
188 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
189 ncomponents=18,
190 default_stored_quantity='displacement',
191 provided_components=[
192 'n', 'e', 'd']),
193 ComponentSchemeDescription(
194 name='elastic5',
195 source_terms=['fn', 'fe', 'fd'],
196 ncomponents=5,
197 default_stored_quantity='displacement',
198 provided_components=[
199 'n', 'e', 'd']),
200 ComponentSchemeDescription(
201 name='poroelastic10',
202 source_terms=['pore_pressure'],
203 ncomponents=10,
204 default_stored_quantity=None,
205 provided_components=[
206 'displacement.n', 'displacement.e', 'displacement.d',
207 'vertical_tilt.n', 'vertical_tilt.e',
208 'pore_pressure',
209 'darcy_velocity.n', 'darcy_velocity.e', 'darcy_velocity.d'])]
212# new names?
213# 'mt_to_displacement_1d'
214# 'mt_to_displacement_1d_ff_only'
215# 'mt_to_gravity_1d'
216# 'mt_to_stress_1d'
217# 'explosion_to_displacement_1d'
218# 'sf_to_displacement_1d'
219# 'mt_to_displacement_3d'
220# 'mt_to_gravity_3d'
221# 'mt_to_stress_3d'
222# 'pore_pressure_to_displacement_1d'
223# 'pore_pressure_to_vertical_tilt_1d'
224# 'pore_pressure_to_pore_pressure_1d'
225# 'pore_pressure_to_darcy_velocity_1d'
228component_schemes = [c.name for c in component_scheme_descriptions]
229component_scheme_to_description = dict(
230 (c.name, c) for c in component_scheme_descriptions)
233class ComponentScheme(StringChoice):
234 '''
235 Different Green's Function component schemes are available:
237 ================= =========================================================
238 Name Description
239 ================= =========================================================
240 ``elastic10`` Elastodynamic for
241 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
242 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT
243 sources only
244 ``elastic8`` Elastodynamic for far-field only
245 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
246 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores,
247 MT sources only
248 ``elastic2`` Elastodynamic for
249 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
250 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, purely
251 isotropic sources only
252 ``elastic5`` Elastodynamic for
253 :py:class:`~pyrocko.gf.meta.ConfigTypeA`
254 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, SF
255 sources only
256 ``elastic18`` Elastodynamic for
257 :py:class:`~pyrocko.gf.meta.ConfigTypeC` stores, MT
258 sources only
259 ``poroelastic10`` Poroelastic for :py:class:`~pyrocko.gf.meta.ConfigTypeA`
260 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores
261 ``scalar1`` Scalar store e.g. acoustic pressure
262 for :py:class:`~pyrocko.gf.meta.ConfigTypeA`
263 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores
264 ================= =========================================================
265 '''
267 choices = component_schemes
270class Earthmodel1D(Object):
271 dummy_for = cake.LayeredModel
273 class __T(TBase):
274 def regularize_extra(self, val):
275 if isinstance(val, str):
276 val = cake.LayeredModel.from_scanlines(
277 cake.read_nd_model_str(val))
279 return val
281 def to_save(self, val):
282 return literal(cake.write_nd_model_str(val))
285class StringID(StringPattern):
286 pattern = r'^[A-Za-z][A-Za-z0-9._]{0,64}$'
289class ScopeType(StringChoice):
290 choices = [
291 'global',
292 'regional',
293 'local',
294 ]
297class WaveType(StringChoice):
298 choices = [
299 'full waveform',
300 'bodywave',
301 'P wave',
302 'S wave',
303 'surface wave',
304 ]
307class NearfieldTermsType(StringChoice):
308 choices = [
309 'complete',
310 'incomplete',
311 'missing',
312 ]
315class QuantityType(StringChoice):
316 choices = [
317 'displacement',
318 'rotation',
319 'velocity',
320 'acceleration',
321 'pressure',
322 'volume_change',
323 'tilt',
324 'pore_pressure',
325 'pressure',
326 'darcy_velocity',
327 'vertical_tilt']
330class Reference(Object):
331 id = StringID.T()
332 type = String.T()
333 title = Unicode.T()
334 journal = Unicode.T(optional=True)
335 volume = Unicode.T(optional=True)
336 number = Unicode.T(optional=True)
337 pages = Unicode.T(optional=True)
338 year = String.T()
339 note = Unicode.T(optional=True)
340 issn = String.T(optional=True)
341 doi = String.T(optional=True)
342 url = String.T(optional=True)
343 eprint = String.T(optional=True)
344 authors = List.T(Unicode.T())
345 publisher = Unicode.T(optional=True)
346 keywords = Unicode.T(optional=True)
347 note = Unicode.T(optional=True)
348 abstract = Unicode.T(optional=True)
350 @classmethod
351 def from_bibtex(cls, filename=None, stream=None):
353 from pybtex.database.input import bibtex
355 parser = bibtex.Parser()
357 if filename is not None:
358 bib_data = parser.parse_file(filename)
359 elif stream is not None:
360 bib_data = parser.parse_stream(stream)
362 references = []
364 for id_, entry in bib_data.entries.items():
365 d = {}
366 avail = entry.fields.keys()
367 for prop in cls.T.properties:
368 if prop.name == 'authors' or prop.name not in avail:
369 continue
371 d[prop.name] = entry.fields[prop.name]
373 if 'author' in entry.persons:
374 d['authors'] = []
375 for person in entry.persons['author']:
376 d['authors'].append(new_str(person))
378 c = Reference(id=id_, type=entry.type, **d)
379 references.append(c)
381 return references
384_fpat = r'[+-]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?'
385_spat = StringID.pattern[1:-1]
386_cake_pat = cake.PhaseDef.allowed_characters_pattern
387_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic
389_ppat = r'(' \
390 r'cake:' + _cake_pat + \
391 r'|iaspei:' + _iaspei_pat + \
392 r'|vel_surface:' + _fpat + \
393 r'|vel:' + _fpat + \
394 r'|stored:' + _spat + \
395 r')'
398timing_regex_old = re.compile(
399 r'^((first|last)?\((' + _spat + r'(\|' + _spat + r')*)\)|(' +
400 _spat + r'))?(' + _fpat + ')?$')
403timing_regex = re.compile(
404 r'^((first|last)?\{(' + _ppat + r'(\|' + _ppat + r')*)\}|(' +
405 _ppat + r'))?(' + _fpat + '(S|%)?)?$')
408class PhaseSelect(StringChoice):
409 choices = ['', 'first', 'last']
412class InvalidTimingSpecification(ValidationError):
413 pass
416class Timing(SObject):
417 '''
418 Definition of a time instant relative to one or more named phase arrivals.
420 Instances of this class can be used e.g. in cutting and tapering
421 operations. They can hold an absolute time or an offset to a named phase
422 arrival or group of such arrivals.
424 Timings can be instantiated from a simple string defintion i.e. with
425 ``Timing(str)`` where ``str`` is something like
426 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``,
427 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase
428 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'``
429 is appended, it is interpreted as percent. If the an ``'S'`` is appended
430 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km].
432 Phase definitions can be specified in either of the following ways:
434 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table
435 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake
436 (see :py:class:`pyrocko.cake.PhaseDef`)
437 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance /
438 velocity [km/s]
439 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s]
441 **Examples:**
443 * ``'100'`` : absolute time; 100 s
444 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to
445 stored travel time table named ``'P'``
446 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to
447 stored travel time table named ``'P'``
448 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to
449 stored travel time table named ``'P'``
450 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is
451 undefined for a given geometry
452 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A
453 and B is chosen, if both phases are defined for a given geometry
454 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen
455 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of
456 arrivals A, B, and C
457 '''
459 def __init__(self, s=None, **kwargs):
461 if s is not None:
462 offset_is = None
463 s = re.sub(r'\s+', '', s)
464 try:
465 offset = float(s.rstrip('S'))
467 if s.endswith('S'):
468 offset_is = 'slowness'
470 phase_defs = []
471 select = ''
473 except ValueError:
474 matched = False
475 m = timing_regex.match(s)
476 if m:
477 if m.group(3):
478 phase_defs = m.group(3).split('|')
479 elif m.group(19):
480 phase_defs = [m.group(19)]
481 else:
482 phase_defs = []
484 select = m.group(2) or ''
486 offset = 0.0
487 soff = m.group(27)
488 if soff:
489 offset = float(soff.rstrip('S%'))
490 if soff.endswith('S'):
491 offset_is = 'slowness'
492 elif soff.endswith('%'):
493 offset_is = 'percent'
495 matched = True
497 else:
498 m = timing_regex_old.match(s)
499 if m:
500 if m.group(3):
501 phase_defs = m.group(3).split('|')
502 elif m.group(5):
503 phase_defs = [m.group(5)]
504 else:
505 phase_defs = []
507 select = m.group(2) or ''
509 offset = 0.0
510 if m.group(6):
511 offset = float(m.group(6))
513 phase_defs = [
514 'stored:' + phase_def for phase_def in phase_defs]
516 matched = True
518 if not matched:
519 raise InvalidTimingSpecification(s)
521 kwargs = dict(
522 phase_defs=phase_defs,
523 select=select,
524 offset=offset,
525 offset_is=offset_is)
527 SObject.__init__(self, **kwargs)
529 def __str__(self):
530 s = []
531 if self.phase_defs:
532 sphases = '|'.join(self.phase_defs)
533 # if len(self.phase_defs) > 1 or self.select:
534 sphases = '{%s}' % sphases
536 if self.select:
537 sphases = self.select + sphases
539 s.append(sphases)
541 if self.offset != 0.0 or not self.phase_defs:
542 s.append('%+g' % self.offset)
543 if self.offset_is == 'slowness':
544 s.append('S')
545 elif self.offset_is == 'percent':
546 s.append('%')
548 return ''.join(s)
550 def evaluate(self, get_phase, args):
551 try:
552 if self.offset_is == 'slowness' and self.offset != 0.0:
553 phase_offset = get_phase(
554 'vel_surface:%g' % (1.0/self.offset))
555 offset = phase_offset(args)
556 else:
557 offset = self.offset
559 if self.phase_defs:
560 phases = [
561 get_phase(phase_def) for phase_def in self.phase_defs]
562 times = [phase(args) for phase in phases]
563 if self.offset_is == 'percent':
564 times = [t*(1.+offset/100.)
565 for t in times if t is not None]
566 else:
567 times = [t+offset for t in times if t is not None]
569 if not times:
570 return None
571 elif self.select == 'first':
572 return min(times)
573 elif self.select == 'last':
574 return max(times)
575 else:
576 return times[0]
577 else:
578 return offset
580 except spit.OutOfBounds:
581 raise OutOfBounds(args)
583 phase_defs = List.T(String.T())
584 offset = Float.T(default=0.0)
585 offset_is = String.T(optional=True)
586 select = PhaseSelect.T(
587 default='',
588 help=('Can be either ``\'%s\'``, ``\'%s\'``, or ``\'%s\'``. ' %
589 tuple(PhaseSelect.choices)))
592def mkdefs(s):
593 defs = []
594 for x in s.split(','):
595 try:
596 defs.append(float(x))
597 except ValueError:
598 if x.startswith('!'):
599 defs.extend(cake.PhaseDef.classic(x[1:]))
600 else:
601 defs.append(cake.PhaseDef(x))
603 return defs
606class TPDef(Object):
607 '''
608 Maps an arrival phase identifier to an arrival phase definition.
609 '''
611 id = StringID.T(
612 help='name used to identify the phase')
613 definition = String.T(
614 help='definition of the phase in either cake syntax as defined in '
615 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an '
616 '``!``, as a *classic phase name*, or, if it is a simple '
617 'number, as a constant horizontal velocity.')
619 @property
620 def phases(self):
621 return [x for x in mkdefs(self.definition)
622 if isinstance(x, cake.PhaseDef)]
624 @property
625 def horizontal_velocities(self):
626 return [x for x in mkdefs(self.definition) if isinstance(x, float)]
629class OutOfBounds(Exception):
630 def __init__(self, values=None, reason=None):
631 Exception.__init__(self)
632 self.values = values
633 self.reason = reason
634 self.context = None
636 def __str__(self):
637 scontext = ''
638 if self.context:
639 scontext = '\n%s' % str(self.context)
641 if self.reason:
642 scontext += '\n%s' % self.reason
644 if self.values:
645 return 'out of bounds: (%s)%s' % (
646 ', '.join('%g' % x for x in self.values), scontext)
647 else:
648 return 'out of bounds%s' % scontext
651class MultiLocation(Object):
653 lats = Array.T(
654 optional=True, shape=(None,), dtype=float,
655 help='Latitudes of targets.')
656 lons = Array.T(
657 optional=True, shape=(None,), dtype=float,
658 help='Longitude of targets.')
659 north_shifts = Array.T(
660 optional=True, shape=(None,), dtype=float,
661 help='North shifts of targets.')
662 east_shifts = Array.T(
663 optional=True, shape=(None,), dtype=float,
664 help='East shifts of targets.')
665 elevation = Array.T(
666 optional=True, shape=(None,), dtype=float,
667 help='Elevations of targets.')
669 def __init__(self, *args, **kwargs):
670 self._coords5 = None
671 Object.__init__(self, *args, **kwargs)
673 @property
674 def coords5(self):
675 if self._coords5 is not None:
676 return self._coords5
677 props = [self.lats, self.lons, self.north_shifts, self.east_shifts,
678 self.elevation]
679 sizes = [p.size for p in props if p is not None]
680 if not sizes:
681 raise AttributeError('Empty StaticTarget')
682 if num.unique(sizes).size != 1:
683 raise AttributeError('Inconsistent coordinate sizes.')
685 self._coords5 = num.zeros((sizes[0], 5))
686 for idx, p in enumerate(props):
687 if p is not None:
688 self._coords5[:, idx] = p.flatten()
690 return self._coords5
692 @property
693 def ncoords(self):
694 if self.coords5.shape[0] is None:
695 return 0
696 return int(self.coords5.shape[0])
698 def get_latlon(self):
699 '''
700 Get all coordinates as lat/lon.
702 :returns: Coordinates in Latitude, Longitude
703 :rtype: :class:`numpy.ndarray`, (N, 2)
704 '''
705 latlons = num.empty((self.ncoords, 2))
706 for i in range(self.ncoords):
707 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4])
708 return latlons
710 def get_corner_coords(self):
711 '''
712 Returns the corner coordinates of the multi-location object.
714 :returns: In lat/lon: lower left, upper left, upper right, lower right
715 :rtype: tuple
716 '''
717 latlon = self.get_latlon()
718 ll = (latlon[:, 0].min(), latlon[:, 1].min())
719 ur = (latlon[:, 0].max(), latlon[:, 1].max())
720 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1]))
723class Receiver(Location):
724 codes = Tuple.T(
725 3, String.T(),
726 optional=True,
727 help='network, station, and location codes')
729 def pyrocko_station(self):
730 from pyrocko import model
732 lat, lon = self.effective_latlon
734 return model.Station(
735 network=self.codes[0],
736 station=self.codes[1],
737 location=self.codes[2],
738 lat=lat,
739 lon=lon,
740 depth=self.depth)
743def g(x, d):
744 if x is None:
745 return d
746 else:
747 return x
750class UnavailableScheme(Exception):
751 pass
754class InvalidNComponents(Exception):
755 pass
758class DiscretizedSource(Object):
759 '''
760 Base class for discretized sources.
762 To compute synthetic seismograms, the parameterized source models (see of
763 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first
764 discretized into a number of point sources. These spacio-temporal point
765 source distributions are represented by subclasses of the
766 :py:class:`DiscretizedSource`. For elastodynamic problems there is the
767 :py:class:`DiscretizedMTSource` for moment tensor point source
768 distributions and the :py:class:`DiscretizedExplosionSource` for pure
769 explosion/implosion type source distributions. The geometry calculations
770 are implemented in the base class. How Green's function components have to
771 be weighted and summed is defined in the derived classes.
773 Like in the :py:class:`Location` class, the positions of the point sources
774 contained in the discretized source are defined by a common reference point
775 (:py:attr:`lat`, :py:attr:`lon`) and cartesian offsets to this
776 (:py:attr:`north_shifts`, :py:attr:`east_shifts`, :py:attr:`depths`).
777 Alternatively latitude and longitude of each contained point source can be
778 specified directly (:py:attr:`lats`, :py:attr:`lons`).
779 '''
781 times = Array.T(shape=(None,), dtype=float)
782 lats = Array.T(shape=(None,), dtype=float, optional=True)
783 lons = Array.T(shape=(None,), dtype=float, optional=True)
784 lat = Float.T(optional=True)
785 lon = Float.T(optional=True)
786 north_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
787 east_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
788 depths = Array.T(shape=(None,), dtype=num.float64)
789 dl = Float.T(optional=True)
790 dw = Float.T(optional=True)
791 nl = Float.T(optional=True)
792 nw = Float.T(optional=True)
794 @classmethod
795 def check_scheme(cls, scheme):
796 '''
797 Check if given GF component scheme is supported by the class.
799 Raises :py:class:`UnavailableScheme` if the given scheme is not
800 supported by this discretized source class.
801 '''
803 if scheme not in cls.provided_schemes:
804 raise UnavailableScheme(
805 'source type "%s" does not support GF component scheme "%s"' %
806 (cls.__name__, scheme))
808 def __init__(self, **kwargs):
809 Object.__init__(self, **kwargs)
810 self._latlons = None
812 def __setattr__(self, name, value):
813 if name in ('lat', 'lon', 'north_shifts', 'east_shifts',
814 'lats', 'lons'):
815 self.__dict__['_latlons'] = None
817 Object.__setattr__(self, name, value)
819 def get_source_terms(self, scheme):
820 raise NotImplementedError()
822 def make_weights(self, receiver, scheme):
823 raise NotImplementedError()
825 @property
826 def effective_latlons(self):
827 '''
828 Property holding the offest-corrected lats and lons of all points.
829 '''
831 if self._latlons is None:
832 if self.lats is not None and self.lons is not None:
833 if (self.north_shifts is not None and
834 self.east_shifts is not None):
835 self._latlons = orthodrome.ne_to_latlon(
836 self.lats, self.lons,
837 self.north_shifts, self.east_shifts)
838 else:
839 self._latlons = self.lats, self.lons
840 else:
841 lat = g(self.lat, 0.0)
842 lon = g(self.lon, 0.0)
843 self._latlons = orthodrome.ne_to_latlon(
844 lat, lon, self.north_shifts, self.east_shifts)
846 return self._latlons
848 @property
849 def effective_north_shifts(self):
850 if self.north_shifts is not None:
851 return self.north_shifts
852 else:
853 return num.zeros(self.times.size)
855 @property
856 def effective_east_shifts(self):
857 if self.east_shifts is not None:
858 return self.east_shifts
859 else:
860 return num.zeros(self.times.size)
862 def same_origin(self, receiver):
863 '''
864 Check if receiver has same reference point.
865 '''
867 return (g(self.lat, 0.0) == receiver.lat and
868 g(self.lon, 0.0) == receiver.lon and
869 self.lats is None and self.lons is None)
871 def azibazis_to(self, receiver):
872 '''
873 Compute azimuths and backaziumuths to/from receiver, for all contained
874 points.
875 '''
877 if self.same_origin(receiver):
878 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts,
879 receiver.north_shift - self.north_shifts)
880 bazis = azis + 180.
881 else:
882 slats, slons = self.effective_latlons
883 rlat, rlon = receiver.effective_latlon
884 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon)
885 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons)
887 return azis, bazis
889 def distances_to(self, receiver):
890 '''
891 Compute distances to receiver for all contained points.
892 '''
893 if self.same_origin(receiver):
894 return num.sqrt((self.north_shifts - receiver.north_shift)**2 +
895 (self.east_shifts - receiver.east_shift)**2)
897 else:
898 slats, slons = self.effective_latlons
899 rlat, rlon = receiver.effective_latlon
900 return orthodrome.distance_accurate50m_numpy(slats, slons,
901 rlat, rlon)
903 def element_coords(self, i):
904 if self.lats is not None and self.lons is not None:
905 lat = float(self.lats[i])
906 lon = float(self.lons[i])
907 else:
908 lat = self.lat
909 lon = self.lon
911 if self.north_shifts is not None and self.east_shifts is not None:
912 north_shift = float(self.north_shifts[i])
913 east_shift = float(self.east_shifts[i])
915 else:
916 north_shift = east_shift = 0.0
918 return lat, lon, north_shift, east_shift
920 def coords5(self):
921 xs = num.zeros((self.nelements, 5))
923 if self.lats is not None and self.lons is not None:
924 xs[:, 0] = self.lats
925 xs[:, 1] = self.lons
926 else:
927 xs[:, 0] = g(self.lat, 0.0)
928 xs[:, 1] = g(self.lon, 0.0)
930 if self.north_shifts is not None and self.east_shifts is not None:
931 xs[:, 2] = self.north_shifts
932 xs[:, 3] = self.east_shifts
934 xs[:, 4] = self.depths
936 return xs
938 @property
939 def nelements(self):
940 return self.times.size
942 @classmethod
943 def combine(cls, sources, **kwargs):
944 '''
945 Combine several discretized source models.
947 Concatenenates all point sources in the given discretized ``sources``.
948 Care must be taken when using this function that the external amplitude
949 factors and reference times of the parameterized (undiscretized)
950 sources match or are accounted for.
951 '''
953 first = sources[0]
955 if not all(type(s) == type(first) for s in sources):
956 raise Exception('DiscretizedSource.combine must be called with '
957 'sources of same type.')
959 latlons = []
960 for s in sources:
961 latlons.append(s.effective_latlons)
963 lats, lons = num.hstack(latlons)
965 if all((s.lats is None and s.lons is None) for s in sources):
966 rlats = num.array([s.lat for s in sources], dtype=float)
967 rlons = num.array([s.lon for s in sources], dtype=float)
968 same_ref = num.all(
969 rlats == rlats[0]) and num.all(rlons == rlons[0])
970 else:
971 same_ref = False
973 cat = num.concatenate
974 times = cat([s.times for s in sources])
975 depths = cat([s.depths for s in sources])
977 if same_ref:
978 lat = first.lat
979 lon = first.lon
980 north_shifts = cat([s.effective_north_shifts for s in sources])
981 east_shifts = cat([s.effective_east_shifts for s in sources])
982 lats = None
983 lons = None
984 else:
985 lat = None
986 lon = None
987 north_shifts = None
988 east_shifts = None
990 return cls(
991 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
992 north_shifts=north_shifts, east_shifts=east_shifts,
993 depths=depths, **kwargs)
995 def centroid_position(self):
996 moments = self.moments()
997 norm = num.sum(moments)
998 if norm != 0.0:
999 w = moments / num.sum(moments)
1000 else:
1001 w = num.ones(moments.size)
1003 if self.lats is not None and self.lons is not None:
1004 lats, lons = self.effective_latlons
1005 rlat, rlon = num.mean(lats), num.mean(lons)
1006 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
1007 else:
1008 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
1009 n, e = self.north_shifts, self.east_shifts
1011 cn = num.sum(n*w)
1012 ce = num.sum(e*w)
1013 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1015 if self.lats is not None and self.lons is not None:
1016 lat = clat
1017 lon = clon
1018 north_shift = 0.
1019 east_shift = 0.
1020 else:
1021 lat = g(self.lat, 0.0)
1022 lon = g(self.lon, 0.0)
1023 north_shift = cn
1024 east_shift = ce
1026 depth = num.sum(self.depths*w)
1027 time = num.sum(self.times*w)
1028 return tuple(float(x) for x in
1029 (time, lat, lon, north_shift, east_shift, depth))
1032class DiscretizedExplosionSource(DiscretizedSource):
1033 m0s = Array.T(shape=(None,), dtype=float)
1035 provided_schemes = (
1036 'scalar1',
1037 'elastic2',
1038 'elastic8',
1039 'elastic10',
1040 )
1042 def get_source_terms(self, scheme):
1043 self.check_scheme(scheme)
1045 if scheme in ('scalar1', 'elastic2'):
1046 return self.m0s[:, num.newaxis].copy()
1048 elif scheme in ('elastic8', 'elastic10'):
1049 m6s = num.zeros((self.m0s.size, 6))
1050 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1051 return m6s
1052 else:
1053 assert False
1055 def make_weights(self, receiver, scheme):
1056 self.check_scheme(scheme)
1058 azis, bazis = self.azibazis_to(receiver)
1060 sb = num.sin(bazis*d2r-num.pi)
1061 cb = num.cos(bazis*d2r-num.pi)
1063 m0s = self.m0s
1064 n = azis.size
1066 cat = num.concatenate
1067 rep = num.repeat
1069 if scheme == 'elastic2':
1070 w_n = cb*m0s
1071 g_n = filledi(0, n)
1072 w_e = sb*m0s
1073 g_e = filledi(0, n)
1074 w_d = m0s
1075 g_d = filledi(1, n)
1077 elif scheme == 'elastic8':
1078 w_n = cat((cb*m0s, cb*m0s))
1079 g_n = rep((0, 2), n)
1080 w_e = cat((sb*m0s, sb*m0s))
1081 g_e = rep((0, 2), n)
1082 w_d = cat((m0s, m0s))
1083 g_d = rep((5, 7), n)
1085 elif scheme == 'elastic10':
1086 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1087 g_n = rep((0, 2, 8), n)
1088 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1089 g_e = rep((0, 2, 8), n)
1090 w_d = cat((m0s, m0s, m0s))
1091 g_d = rep((5, 7, 9), n)
1093 else:
1094 assert False
1096 return (
1097 ('displacement.n', w_n, g_n),
1098 ('displacement.e', w_e, g_e),
1099 ('displacement.d', w_d, g_d),
1100 )
1102 def split(self):
1103 from pyrocko.gf.seismosizer import ExplosionSource
1104 sources = []
1105 for i in range(self.nelements):
1106 lat, lon, north_shift, east_shift = self.element_coords(i)
1107 sources.append(ExplosionSource(
1108 time=float(self.times[i]),
1109 lat=lat,
1110 lon=lon,
1111 north_shift=north_shift,
1112 east_shift=east_shift,
1113 depth=float(self.depths[i]),
1114 moment=float(self.m0s[i])))
1116 return sources
1118 def moments(self):
1119 return self.m0s
1121 def centroid(self):
1122 from pyrocko.gf.seismosizer import ExplosionSource
1123 time, lat, lon, north_shift, east_shift, depth = \
1124 self.centroid_position()
1126 return ExplosionSource(
1127 time=time,
1128 lat=lat,
1129 lon=lon,
1130 north_shift=north_shift,
1131 east_shift=east_shift,
1132 depth=depth,
1133 moment=float(num.sum(self.m0s)))
1135 @classmethod
1136 def combine(cls, sources, **kwargs):
1137 '''
1138 Combine several discretized source models.
1140 Concatenenates all point sources in the given discretized ``sources``.
1141 Care must be taken when using this function that the external amplitude
1142 factors and reference times of the parameterized (undiscretized)
1143 sources match or are accounted for.
1144 '''
1146 if 'm0s' not in kwargs:
1147 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1149 return super(DiscretizedExplosionSource, cls).combine(sources,
1150 **kwargs)
1153class DiscretizedSFSource(DiscretizedSource):
1154 forces = Array.T(shape=(None, 3), dtype=float)
1156 provided_schemes = (
1157 'elastic5',
1158 )
1160 def get_source_terms(self, scheme):
1161 self.check_scheme(scheme)
1163 return self.forces
1165 def make_weights(self, receiver, scheme):
1166 self.check_scheme(scheme)
1168 azis, bazis = self.azibazis_to(receiver)
1170 sa = num.sin(azis*d2r)
1171 ca = num.cos(azis*d2r)
1172 sb = num.sin(bazis*d2r-num.pi)
1173 cb = num.cos(bazis*d2r-num.pi)
1175 forces = self.forces
1176 fn = forces[:, 0]
1177 fe = forces[:, 1]
1178 fd = forces[:, 2]
1180 f0 = fd
1181 f1 = ca * fn + sa * fe
1182 f2 = ca * fe - sa * fn
1184 n = azis.size
1186 if scheme == 'elastic5':
1187 ioff = 0
1189 cat = num.concatenate
1190 rep = num.repeat
1192 w_n = cat((cb*f0, cb*f1, -sb*f2))
1193 g_n = ioff + rep((0, 1, 2), n)
1194 w_e = cat((sb*f0, sb*f1, cb*f2))
1195 g_e = ioff + rep((0, 1, 2), n)
1196 w_d = cat((f0, f1))
1197 g_d = ioff + rep((3, 4), n)
1199 return (
1200 ('displacement.n', w_n, g_n),
1201 ('displacement.e', w_e, g_e),
1202 ('displacement.d', w_d, g_d),
1203 )
1205 @classmethod
1206 def combine(cls, sources, **kwargs):
1207 '''
1208 Combine several discretized source models.
1210 Concatenenates all point sources in the given discretized ``sources``.
1211 Care must be taken when using this function that the external amplitude
1212 factors and reference times of the parameterized (undiscretized)
1213 sources match or are accounted for.
1214 '''
1216 if 'forces' not in kwargs:
1217 kwargs['forces'] = num.vstack([s.forces for s in sources])
1219 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1221 def moments(self):
1222 return num.sum(self.forces**2, axis=1)
1224 def centroid(self):
1225 from pyrocko.gf.seismosizer import SFSource
1226 time, lat, lon, north_shift, east_shift, depth = \
1227 self.centroid_position()
1229 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1230 return SFSource(
1231 time=time,
1232 lat=lat,
1233 lon=lon,
1234 north_shift=north_shift,
1235 east_shift=east_shift,
1236 depth=depth,
1237 fn=fn,
1238 fe=fe,
1239 fd=fd)
1242class DiscretizedMTSource(DiscretizedSource):
1243 m6s = Array.T(
1244 shape=(None, 6), dtype=float,
1245 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1247 provided_schemes = (
1248 'elastic8',
1249 'elastic10',
1250 'elastic18',
1251 )
1253 def get_source_terms(self, scheme):
1254 self.check_scheme(scheme)
1255 return self.m6s
1257 def make_weights(self, receiver, scheme):
1258 self.check_scheme(scheme)
1260 m6s = self.m6s
1261 n = m6s.shape[0]
1262 rep = num.repeat
1264 if scheme == 'elastic18':
1265 w_n = m6s.flatten()
1266 g_n = num.tile(num.arange(0, 6), n)
1267 w_e = m6s.flatten()
1268 g_e = num.tile(num.arange(6, 12), n)
1269 w_d = m6s.flatten()
1270 g_d = num.tile(num.arange(12, 18), n)
1272 else:
1273 azis, bazis = self.azibazis_to(receiver)
1275 sa = num.sin(azis*d2r)
1276 ca = num.cos(azis*d2r)
1277 s2a = num.sin(2.*azis*d2r)
1278 c2a = num.cos(2.*azis*d2r)
1279 sb = num.sin(bazis*d2r-num.pi)
1280 cb = num.cos(bazis*d2r-num.pi)
1282 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1283 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1284 f2 = m6s[:, 2]
1285 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1286 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1287 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1289 cat = num.concatenate
1291 if scheme == 'elastic8':
1292 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1293 g_n = rep((0, 1, 2, 3, 4), n)
1294 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1295 g_e = rep((0, 1, 2, 3, 4), n)
1296 w_d = cat((f0, f1, f2))
1297 g_d = rep((5, 6, 7), n)
1299 elif scheme == 'elastic10':
1300 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1301 g_n = rep((0, 1, 2, 8, 3, 4), n)
1302 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1303 g_e = rep((0, 1, 2, 8, 3, 4), n)
1304 w_d = cat((f0, f1, f2, f5))
1305 g_d = rep((5, 6, 7, 9), n)
1307 else:
1308 assert False
1310 return (
1311 ('displacement.n', w_n, g_n),
1312 ('displacement.e', w_e, g_e),
1313 ('displacement.d', w_d, g_d),
1314 )
1316 def split(self):
1317 from pyrocko.gf.seismosizer import MTSource
1318 sources = []
1319 for i in range(self.nelements):
1320 lat, lon, north_shift, east_shift = self.element_coords(i)
1321 sources.append(MTSource(
1322 time=float(self.times[i]),
1323 lat=lat,
1324 lon=lon,
1325 north_shift=north_shift,
1326 east_shift=east_shift,
1327 depth=float(self.depths[i]),
1328 m6=self.m6s[i]))
1330 return sources
1332 def moments(self):
1333 moments = num.array(
1334 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))
1335 for m6 in self.m6s])
1336 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)
1338 def get_moment_rate(self, deltat=None):
1339 moments = self.moments()
1340 times = self.times
1341 times -= times.min()
1343 t_max = times.max()
1344 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat
1345 mom_times[mom_times > t_max] = t_max
1347 # Right open histrogram bins
1348 mom, _ = num.histogram(
1349 times,
1350 bins=mom_times,
1351 weights=moments)
1353 deltat = num.diff(mom_times)
1354 mom_rate = mom / deltat
1356 return mom_rate, mom_times[1:]
1358 def centroid(self):
1359 from pyrocko.gf.seismosizer import MTSource
1360 time, lat, lon, north_shift, east_shift, depth = \
1361 self.centroid_position()
1363 return MTSource(
1364 time=time,
1365 lat=lat,
1366 lon=lon,
1367 north_shift=north_shift,
1368 east_shift=east_shift,
1369 depth=depth,
1370 m6=num.sum(self.m6s, axis=0))
1372 @classmethod
1373 def combine(cls, sources, **kwargs):
1374 '''
1375 Combine several discretized source models.
1377 Concatenenates all point sources in the given discretized ``sources``.
1378 Care must be taken when using this function that the external amplitude
1379 factors and reference times of the parameterized (undiscretized)
1380 sources match or are accounted for.
1381 '''
1383 if 'm6s' not in kwargs:
1384 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1386 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1389class DiscretizedPorePressureSource(DiscretizedSource):
1390 pp = Array.T(shape=(None,), dtype=float)
1392 provided_schemes = (
1393 'poroelastic10',
1394 )
1396 def get_source_terms(self, scheme):
1397 self.check_scheme(scheme)
1398 return self.pp[:, num.newaxis].copy()
1400 def make_weights(self, receiver, scheme):
1401 self.check_scheme(scheme)
1403 azis, bazis = self.azibazis_to(receiver)
1405 sb = num.sin(bazis*d2r-num.pi)
1406 cb = num.cos(bazis*d2r-num.pi)
1408 pp = self.pp
1409 n = bazis.size
1411 w_un = cb*pp
1412 g_un = filledi(1, n)
1413 w_ue = sb*pp
1414 g_ue = filledi(1, n)
1415 w_ud = pp
1416 g_ud = filledi(0, n)
1418 w_tn = cb*pp
1419 g_tn = filledi(6, n)
1420 w_te = sb*pp
1421 g_te = filledi(6, n)
1423 w_pp = pp
1424 g_pp = filledi(7, n)
1426 w_dvn = cb*pp
1427 g_dvn = filledi(9, n)
1428 w_dve = sb*pp
1429 g_dve = filledi(9, n)
1430 w_dvd = pp
1431 g_dvd = filledi(8, n)
1433 return (
1434 ('displacement.n', w_un, g_un),
1435 ('displacement.e', w_ue, g_ue),
1436 ('displacement.d', w_ud, g_ud),
1437 ('vertical_tilt.n', w_tn, g_tn),
1438 ('vertical_tilt.e', w_te, g_te),
1439 ('pore_pressure', w_pp, g_pp),
1440 ('darcy_velocity.n', w_dvn, g_dvn),
1441 ('darcy_velocity.e', w_dve, g_dve),
1442 ('darcy_velocity.d', w_dvd, g_dvd),
1443 )
1445 def moments(self):
1446 return self.pp
1448 def centroid(self):
1449 from pyrocko.gf.seismosizer import PorePressurePointSource
1450 time, lat, lon, north_shift, east_shift, depth = \
1451 self.centroid_position()
1453 return PorePressurePointSource(
1454 time=time,
1455 lat=lat,
1456 lon=lon,
1457 north_shift=north_shift,
1458 east_shift=east_shift,
1459 depth=depth,
1460 pp=float(num.sum(self.pp)))
1462 @classmethod
1463 def combine(cls, sources, **kwargs):
1464 '''
1465 Combine several discretized source models.
1467 Concatenenates all point sources in the given discretized ``sources``.
1468 Care must be taken when using this function that the external amplitude
1469 factors and reference times of the parameterized (undiscretized)
1470 sources match or are accounted for.
1471 '''
1473 if 'pp' not in kwargs:
1474 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1476 return super(DiscretizedPorePressureSource, cls).combine(sources,
1477 **kwargs)
1480class Region(Object):
1481 name = String.T(optional=True)
1484class RectangularRegion(Region):
1485 lat_min = Float.T()
1486 lat_max = Float.T()
1487 lon_min = Float.T()
1488 lon_max = Float.T()
1491class CircularRegion(Region):
1492 lat = Float.T()
1493 lon = Float.T()
1494 radius = Float.T()
1497class Config(Object):
1498 '''
1499 Green's function store meta information.
1501 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1502 configuration types are:
1504 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1505 spherical symmetry, 1D earth model, single receiver depth
1507 * Problem is invariant to horizontal translations and rotations around
1508 vertical axis.
1509 * All receivers must be at the same depth (e.g. at the surface)
1510 * High level index variables: ``(source_depth, receiver_distance,
1511 component)``
1513 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1514 spherical symmetry, 1D earth model, variable receiver depth
1516 * Symmetries like in Type A but has additional index for receiver depth
1517 * High level index variables: ``(source_depth, receiver_distance,
1518 receiver_depth, component)``
1520 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1521 constraints but fixed receiver positions
1523 * Cartesian source volume around a reference point
1524 * High level index variables: ``(ireceiver, source_depth,
1525 source_east_shift, source_north_shift, component)``
1526 '''
1528 id = StringID.T(
1529 help='Name of the store. May consist of upper and lower-case letters, '
1530 'digits, dots and underscores. The name must start with a '
1531 'letter.')
1533 derived_from_id = StringID.T(
1534 optional=True,
1535 help='Name of the original store, if this store has been derived from '
1536 'another one (e.g. extracted subset).')
1538 version = String.T(
1539 default='1.0',
1540 optional=True,
1541 help='User-defined version string. Use <major>.<minor> format.')
1543 modelling_code_id = StringID.T(
1544 optional=True,
1545 help='Identifier of the backend used to compute the store.')
1547 author = Unicode.T(
1548 optional=True,
1549 help='Comma-separated list of author names.')
1551 author_email = String.T(
1552 optional=True,
1553 help="Author's contact email address.")
1555 created_time = Timestamp.T(
1556 optional=True,
1557 help='Time of creation of the store.')
1559 regions = List.T(
1560 Region.T(),
1561 help='Geographical regions for which the store is representative.')
1563 scope_type = ScopeType.T(
1564 optional=True,
1565 help='Distance range scope of the store '
1566 '(%s).' % fmt_choices(ScopeType))
1568 waveform_type = WaveType.T(
1569 optional=True,
1570 help='Wave type stored (%s).' % fmt_choices(WaveType))
1572 nearfield_terms = NearfieldTermsType.T(
1573 optional=True,
1574 help='Information about the inclusion of near-field terms in the '
1575 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1577 description = String.T(
1578 optional=True,
1579 help='Free form textual description of the GF store.')
1581 references = List.T(
1582 Reference.T(),
1583 help='Reference list to cite the modelling code, earth model or '
1584 'related work.')
1586 earthmodel_1d = Earthmodel1D.T(
1587 optional=True,
1588 help='Layered earth model in ND (named discontinuity) format.')
1590 earthmodel_receiver_1d = Earthmodel1D.T(
1591 optional=True,
1592 help='Receiver-side layered earth model in ND format.')
1594 can_interpolate_source = Bool.T(
1595 optional=True,
1596 help='Hint to indicate if the spatial sampling of the store is dense '
1597 'enough for multi-linear interpolation at the source.')
1599 can_interpolate_receiver = Bool.T(
1600 optional=True,
1601 help='Hint to indicate if the spatial sampling of the store is dense '
1602 'enough for multi-linear interpolation at the receiver.')
1604 frequency_min = Float.T(
1605 optional=True,
1606 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1608 frequency_max = Float.T(
1609 optional=True,
1610 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1612 sample_rate = Float.T(
1613 optional=True,
1614 help='Sample rate of the GF store [Hz].')
1616 factor = Float.T(
1617 default=1.0,
1618 help='Gain value, factored out of the stored GF samples. '
1619 '(may not work properly, keep at 1.0).',
1620 optional=True)
1622 component_scheme = ComponentScheme.T(
1623 default='elastic10',
1624 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1626 stored_quantity = QuantityType.T(
1627 optional=True,
1628 help='Physical quantity of stored values (%s). If not given, a '
1629 'default is used based on the GF component scheme. The default '
1630 'for the ``"elastic*"`` family of component schemes is '
1631 '``"displacement"``.' % fmt_choices(QuantityType))
1633 tabulated_phases = List.T(
1634 TPDef.T(),
1635 help='Mapping of phase names to phase definitions, for which travel '
1636 'time tables are available in the GF store.')
1638 ncomponents = Int.T(
1639 optional=True,
1640 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1642 uuid = String.T(
1643 optional=True,
1644 help='Heuristic hash value which can be used to uniquely identify the '
1645 'GF store for practical purposes.')
1647 reference = String.T(
1648 optional=True,
1649 help='Store reference name composed of the store\'s :gattr:`id` and '
1650 'the first six letters of its :gattr:`uuid`.')
1652 def __init__(self, **kwargs):
1653 self._do_auto_updates = False
1654 Object.__init__(self, **kwargs)
1655 self._index_function = None
1656 self._indices_function = None
1657 self._vicinity_function = None
1658 self.validate(regularize=True, depth=1)
1659 self._do_auto_updates = True
1660 self.update()
1662 def check_ncomponents(self):
1663 ncomponents = component_scheme_to_description[
1664 self.component_scheme].ncomponents
1666 if self.ncomponents is None:
1667 self.ncomponents = ncomponents
1668 elif ncomponents != self.ncomponents:
1669 raise InvalidNComponents(
1670 'ncomponents=%i incompatible with component_scheme="%s"' % (
1671 self.ncomponents, self.component_scheme))
1673 def __setattr__(self, name, value):
1674 Object.__setattr__(self, name, value)
1675 try:
1676 self.T.get_property(name)
1677 if self._do_auto_updates:
1678 self.update()
1680 except ValueError:
1681 pass
1683 def update(self):
1684 self.check_ncomponents()
1685 self._update()
1686 self._make_index_functions()
1688 def irecord(self, *args):
1689 return self._index_function(*args)
1691 def irecords(self, *args):
1692 return self._indices_function(*args)
1694 def vicinity(self, *args):
1695 return self._vicinity_function(*args)
1697 def vicinities(self, *args):
1698 return self._vicinities_function(*args)
1700 def grid_interpolation_coefficients(self, *args):
1701 return self._grid_interpolation_coefficients(*args)
1703 def nodes(self, level=None, minlevel=None):
1704 return nodes(self.coords[minlevel:level])
1706 def iter_nodes(self, level=None, minlevel=None):
1707 return nditer_outer(self.coords[minlevel:level])
1709 def iter_extraction(self, gdef, level=None):
1710 i = 0
1711 arrs = []
1712 ntotal = 1
1713 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1714 if gdef and len(gdef) > i:
1715 sssn = gdef[i]
1716 else:
1717 sssn = (None,)*4
1719 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1720 ntotal *= len(arr)
1722 arrs.append(arr)
1723 i += 1
1725 arrs.append(self.coords[-1])
1726 return nditer_outer(arrs[:level])
1728 def make_sum_params(self, source, receiver, implementation='c',
1729 nthreads=0):
1730 assert implementation in ['c', 'python']
1732 out = []
1733 delays = source.times
1734 for comp, weights, icomponents in source.make_weights(
1735 receiver,
1736 self.component_scheme):
1738 weights *= self.factor
1740 args = self.make_indexing_args(source, receiver, icomponents)
1741 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1742 out.append((comp, args, delays_expanded, weights))
1744 return out
1746 def short_info(self):
1747 raise NotImplementedError('should be implemented in subclass')
1749 def get_shear_moduli(self, lat, lon, points,
1750 interpolation=None):
1751 '''
1752 Get shear moduli at given points from contained velocity model.
1754 :param lat: surface origin for coordinate system of ``points``
1755 :param points: NumPy array of shape ``(N, 3)``, where each row is
1756 a point ``(north, east, depth)``, relative to origin at
1757 ``(lat, lon)``
1758 :param interpolation: interpolation method. Choose from
1759 ``('nearest_neighbor', 'multilinear')``
1760 :returns: NumPy array of length N with extracted shear moduli at each
1761 point
1763 The default implementation retrieves and interpolates the shear moduli
1764 from the contained 1D velocity profile.
1765 '''
1766 return self.get_material_property(lat, lon, points,
1767 parameter='shear_moduli',
1768 interpolation=interpolation)
1770 def get_lambda_moduli(self, lat, lon, points,
1771 interpolation=None):
1772 '''
1773 Get lambda moduli at given points from contained velocity model.
1775 :param lat: surface origin for coordinate system of ``points``
1776 :param points: NumPy array of shape ``(N, 3)``, where each row is
1777 a point ``(north, east, depth)``, relative to origin at
1778 ``(lat, lon)``
1779 :param interpolation: interpolation method. Choose from
1780 ``('nearest_neighbor', 'multilinear')``
1781 :returns: NumPy array of length N with extracted shear moduli at each
1782 point
1784 The default implementation retrieves and interpolates the lambda moduli
1785 from the contained 1D velocity profile.
1786 '''
1787 return self.get_material_property(lat, lon, points,
1788 parameter='lambda_moduli',
1789 interpolation=interpolation)
1791 def get_bulk_moduli(self, lat, lon, points,
1792 interpolation=None):
1793 '''
1794 Get bulk moduli at given points from contained velocity model.
1796 :param lat: surface origin for coordinate system of ``points``
1797 :param points: NumPy array of shape ``(N, 3)``, where each row is
1798 a point ``(north, east, depth)``, relative to origin at
1799 ``(lat, lon)``
1800 :param interpolation: interpolation method. Choose from
1801 ``('nearest_neighbor', 'multilinear')``
1802 :returns: NumPy array of length N with extracted shear moduli at each
1803 point
1805 The default implementation retrieves and interpolates the lambda moduli
1806 from the contained 1D velocity profile.
1807 '''
1808 lambda_moduli = self.get_material_property(
1809 lat, lon, points, parameter='lambda_moduli',
1810 interpolation=interpolation)
1811 shear_moduli = self.get_material_property(
1812 lat, lon, points, parameter='shear_moduli',
1813 interpolation=interpolation)
1814 return lambda_moduli + (2 / 3) * shear_moduli
1816 def get_vs(self, lat, lon, points, interpolation=None):
1817 '''
1818 Get Vs at given points from contained velocity model.
1820 :param lat: surface origin for coordinate system of ``points``
1821 :param points: NumPy array of shape ``(N, 3)``, where each row is
1822 a point ``(north, east, depth)``, relative to origin at
1823 ``(lat, lon)``
1824 :param interpolation: interpolation method. Choose from
1825 ``('nearest_neighbor', 'multilinear')``
1826 :returns: NumPy array of length N with extracted shear moduli at each
1827 point
1829 The default implementation retrieves and interpolates Vs
1830 from the contained 1D velocity profile.
1831 '''
1832 return self.get_material_property(lat, lon, points,
1833 parameter='vs',
1834 interpolation=interpolation)
1836 def get_vp(self, lat, lon, points, interpolation=None):
1837 '''
1838 Get Vp at given points from contained velocity model.
1840 :param lat: surface origin for coordinate system of ``points``
1841 :param points: NumPy array of shape ``(N, 3)``, where each row is
1842 a point ``(north, east, depth)``, relative to origin at
1843 ``(lat, lon)``
1844 :param interpolation: interpolation method. Choose from
1845 ``('nearest_neighbor', 'multilinear')``
1846 :returns: NumPy array of length N with extracted shear moduli at each
1847 point
1849 The default implementation retrieves and interpolates Vp
1850 from the contained 1D velocity profile.
1851 '''
1852 return self.get_material_property(lat, lon, points,
1853 parameter='vp',
1854 interpolation=interpolation)
1856 def get_rho(self, lat, lon, points, interpolation=None):
1857 '''
1858 Get rho at given points from contained velocity model.
1860 :param lat: surface origin for coordinate system of ``points``
1861 :param points: NumPy array of shape ``(N, 3)``, where each row is
1862 a point ``(north, east, depth)``, relative to origin at
1863 ``(lat, lon)``
1864 :param interpolation: interpolation method. Choose from
1865 ``('nearest_neighbor', 'multilinear')``
1866 :returns: NumPy array of length N with extracted shear moduli at each
1867 point
1869 The default implementation retrieves and interpolates rho
1870 from the contained 1D velocity profile.
1871 '''
1872 return self.get_material_property(lat, lon, points,
1873 parameter='rho',
1874 interpolation=interpolation)
1876 def get_material_property(self, lat, lon, points, parameter='vs',
1877 interpolation=None):
1879 if interpolation is None:
1880 raise TypeError('Interpolation method not defined! available: '
1881 "multilinear", "nearest_neighbor")
1883 earthmod = self.earthmodel_1d
1884 store_depth_profile = self.get_source_depths()
1885 z_profile = earthmod.profile('z')
1887 if parameter == 'vs':
1888 vs_profile = earthmod.profile('vs')
1889 profile = num.interp(
1890 store_depth_profile, z_profile, vs_profile)
1892 elif parameter == 'vp':
1893 vp_profile = earthmod.profile('vp')
1894 profile = num.interp(
1895 store_depth_profile, z_profile, vp_profile)
1897 elif parameter == 'rho':
1898 rho_profile = earthmod.profile('rho')
1900 profile = num.interp(
1901 store_depth_profile, z_profile, rho_profile)
1903 elif parameter == 'shear_moduli':
1904 vs_profile = earthmod.profile('vs')
1905 rho_profile = earthmod.profile('rho')
1907 store_vs_profile = num.interp(
1908 store_depth_profile, z_profile, vs_profile)
1909 store_rho_profile = num.interp(
1910 store_depth_profile, z_profile, rho_profile)
1912 profile = num.power(store_vs_profile, 2) * store_rho_profile
1914 elif parameter == 'lambda_moduli':
1915 vs_profile = earthmod.profile('vs')
1916 vp_profile = earthmod.profile('vp')
1917 rho_profile = earthmod.profile('rho')
1919 store_vs_profile = num.interp(
1920 store_depth_profile, z_profile, vs_profile)
1921 store_vp_profile = num.interp(
1922 store_depth_profile, z_profile, vp_profile)
1923 store_rho_profile = num.interp(
1924 store_depth_profile, z_profile, rho_profile)
1926 profile = store_rho_profile * (
1927 num.power(store_vp_profile, 2) -
1928 num.power(store_vs_profile, 2) * 2)
1929 else:
1930 raise TypeError(
1931 'parameter %s not available' % parameter)
1933 if interpolation == 'multilinear':
1934 kind = 'linear'
1935 elif interpolation == 'nearest_neighbor':
1936 kind = 'nearest'
1937 else:
1938 raise TypeError(
1939 'Interpolation method %s not available' % interpolation)
1941 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1943 try:
1944 return interpolator(points[:, 2])
1945 except ValueError:
1946 raise OutOfBounds()
1948 def is_static(self):
1949 for code in ('psgrn_pscmp', 'poel'):
1950 if self.modelling_code_id.startswith(code):
1951 return True
1952 return False
1954 def is_dynamic(self):
1955 return not self.is_static()
1957 def get_source_depths(self):
1958 raise NotImplementedError('must be implemented in subclass')
1960 def get_tabulated_phase(self, phase_id):
1961 '''
1962 Get tabulated phase definition.
1963 '''
1965 for pdef in self.tabulated_phases:
1966 if pdef.id == phase_id:
1967 return pdef
1969 raise StoreError('No such phase: %s' % phase_id)
1971 def fix_ttt_holes(self, sptree, mode):
1972 raise StoreError(
1973 'Cannot fix travel time table holes in GF stores of type %s.'
1974 % self.short_type)
1977class ConfigTypeA(Config):
1978 '''
1979 Cylindrical symmetry, 1D earth model, single receiver depth
1981 * Problem is invariant to horizontal translations and rotations around
1982 vertical axis.
1984 * All receivers must be at the same depth (e.g. at the surface)
1985 High level index variables: ``(source_depth, distance,
1986 component)``
1988 * The ``distance`` is the surface distance between source and receiver
1989 points.
1990 '''
1992 receiver_depth = Float.T(
1993 default=0.0,
1994 help='Fixed receiver depth [m].')
1996 source_depth_min = Float.T(
1997 help='Minimum source depth [m].')
1999 source_depth_max = Float.T(
2000 help='Maximum source depth [m].')
2002 source_depth_delta = Float.T(
2003 help='Grid spacing of source depths [m]')
2005 distance_min = Float.T(
2006 help='Minimum source-receiver surface distance [m].')
2008 distance_max = Float.T(
2009 help='Maximum source-receiver surface distance [m].')
2011 distance_delta = Float.T(
2012 help='Grid spacing of source-receiver surface distance [m].')
2014 short_type = 'A'
2016 provided_schemes = [
2017 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2019 def get_surface_distance(self, args):
2020 return args[1]
2022 def get_distance(self, args):
2023 return math.sqrt(args[0]**2 + args[1]**2)
2025 def get_source_depth(self, args):
2026 return args[0]
2028 def get_source_depths(self):
2029 return self.coords[0]
2031 def get_receiver_depth(self, args):
2032 return self.receiver_depth
2034 def _update(self):
2035 self.mins = num.array(
2036 [self.source_depth_min, self.distance_min], dtype=float)
2037 self.maxs = num.array(
2038 [self.source_depth_max, self.distance_max], dtype=float)
2039 self.deltas = num.array(
2040 [self.source_depth_delta, self.distance_delta],
2041 dtype=float)
2042 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2043 vicinity_eps).astype(int) + 1
2044 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2045 self.deltat = 1.0/self.sample_rate
2046 self.nrecords = num.product(self.ns) * self.ncomponents
2047 self.coords = tuple(num.linspace(mi, ma, n) for
2048 (mi, ma, n) in
2049 zip(self.mins, self.effective_maxs, self.ns)) + \
2050 (num.arange(self.ncomponents),)
2052 self.nsource_depths, self.ndistances = self.ns
2054 def _make_index_functions(self):
2056 amin, bmin = self.mins
2057 da, db = self.deltas
2058 na, nb = self.ns
2060 ng = self.ncomponents
2062 def index_function(a, b, ig):
2063 ia = int(round((a - amin) / da))
2064 ib = int(round((b - bmin) / db))
2065 try:
2066 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2067 except ValueError:
2068 raise OutOfBounds()
2070 def indices_function(a, b, ig):
2071 ia = num.round((a - amin) / da).astype(int)
2072 ib = num.round((b - bmin) / db).astype(int)
2073 try:
2074 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2075 except ValueError:
2076 for ia_, ib_, ig_ in zip(ia, ib, ig):
2077 try:
2078 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2079 except ValueError:
2080 raise OutOfBounds()
2082 def grid_interpolation_coefficients(a, b):
2083 ias = indi12((a - amin) / da, na)
2084 ibs = indi12((b - bmin) / db, nb)
2085 return ias, ibs
2087 def vicinity_function(a, b, ig):
2088 ias, ibs = grid_interpolation_coefficients(a, b)
2090 if not (0 <= ig < ng):
2091 raise OutOfBounds()
2093 indis = []
2094 weights = []
2095 for ia, va in ias:
2096 iia = ia*nb*ng
2097 for ib, vb in ibs:
2098 indis.append(iia + ib*ng + ig)
2099 weights.append(va*vb)
2101 return num.array(indis), num.array(weights)
2103 def vicinities_function(a, b, ig):
2105 xa = (a - amin) / da
2106 xb = (b - bmin) / db
2108 xa_fl = num.floor(xa)
2109 xa_ce = num.ceil(xa)
2110 xb_fl = num.floor(xb)
2111 xb_ce = num.ceil(xb)
2112 va_fl = 1.0 - (xa - xa_fl)
2113 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2114 vb_fl = 1.0 - (xb - xb_fl)
2115 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2117 ia_fl = xa_fl.astype(int)
2118 ia_ce = xa_ce.astype(int)
2119 ib_fl = xb_fl.astype(int)
2120 ib_ce = xb_ce.astype(int)
2122 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2123 raise OutOfBounds()
2125 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2126 raise OutOfBounds()
2128 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2129 raise OutOfBounds()
2131 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2132 raise OutOfBounds()
2134 irecords = num.empty(a.size*4, dtype=int)
2135 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2136 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2137 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2138 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2140 weights = num.empty(a.size*4, dtype=float)
2141 weights[0::4] = va_fl * vb_fl
2142 weights[1::4] = va_ce * vb_fl
2143 weights[2::4] = va_fl * vb_ce
2144 weights[3::4] = va_ce * vb_ce
2146 return irecords, weights
2148 self._index_function = index_function
2149 self._indices_function = indices_function
2150 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2151 self._vicinity_function = vicinity_function
2152 self._vicinities_function = vicinities_function
2154 def make_indexing_args(self, source, receiver, icomponents):
2155 nc = icomponents.size
2156 dists = source.distances_to(receiver)
2157 n = dists.size
2158 return (num.tile(source.depths, nc//n),
2159 num.tile(dists, nc//n),
2160 icomponents)
2162 def make_indexing_args1(self, source, receiver):
2163 return (source.depth, source.distance_to(receiver))
2165 @property
2166 def short_extent(self):
2167 return '%g:%g:%g x %g:%g:%g' % (
2168 self.source_depth_min/km,
2169 self.source_depth_max/km,
2170 self.source_depth_delta/km,
2171 self.distance_min/km,
2172 self.distance_max/km,
2173 self.distance_delta/km)
2175 def fix_ttt_holes(self, sptree, mode):
2176 from pyrocko import eikonal_ext, spit
2178 nodes = self.nodes(level=-1)
2180 delta = self.deltas[-1]
2181 assert num.all(delta == self.deltas)
2183 nsources, ndistances = self.ns
2185 points = num.zeros((nodes.shape[0], 3))
2186 points[:, 0] = nodes[:, 1]
2187 points[:, 2] = nodes[:, 0]
2189 speeds = self.get_material_property(
2190 0., 0., points,
2191 parameter='vp' if mode == cake.P else 'vs',
2192 interpolation='multilinear')
2194 speeds = speeds.reshape((nsources, ndistances))
2196 times = sptree.interpolate_many(nodes)
2198 times[num.isnan(times)] = -1.
2199 times = times.reshape(speeds.shape)
2201 try:
2202 eikonal_ext.eikonal_solver_fmm_cartesian(
2203 speeds, times, delta)
2204 except eikonal_ext.EikonalExtError as e:
2205 if str(e).endswith('please check results'):
2206 logger.debug(
2207 'Got a warning from eikonal solver '
2208 '- may be ok...')
2209 else:
2210 raise
2212 def func(x):
2213 ibs, ics = \
2214 self.grid_interpolation_coefficients(*x)
2216 t = 0
2217 for ib, vb in ibs:
2218 for ic, vc in ics:
2219 t += times[ib, ic] * vb * vc
2221 return t
2223 return spit.SPTree(
2224 f=func,
2225 ftol=sptree.ftol,
2226 xbounds=sptree.xbounds,
2227 xtols=sptree.xtols)
2230class ConfigTypeB(Config):
2231 '''
2232 Cylindrical symmetry, 1D earth model, variable receiver depth
2234 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2235 receiver depth
2237 * High level index variables: ``(receiver_depth, source_depth,
2238 receiver_distance, component)``
2239 '''
2241 receiver_depth_min = Float.T(
2242 help='Minimum receiver depth [m].')
2244 receiver_depth_max = Float.T(
2245 help='Maximum receiver depth [m].')
2247 receiver_depth_delta = Float.T(
2248 help='Grid spacing of receiver depths [m]')
2250 source_depth_min = Float.T(
2251 help='Minimum source depth [m].')
2253 source_depth_max = Float.T(
2254 help='Maximum source depth [m].')
2256 source_depth_delta = Float.T(
2257 help='Grid spacing of source depths [m]')
2259 distance_min = Float.T(
2260 help='Minimum source-receiver surface distance [m].')
2262 distance_max = Float.T(
2263 help='Maximum source-receiver surface distance [m].')
2265 distance_delta = Float.T(
2266 help='Grid spacing of source-receiver surface distances [m].')
2268 short_type = 'B'
2270 provided_schemes = [
2271 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2273 def get_distance(self, args):
2274 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2276 def get_surface_distance(self, args):
2277 return args[2]
2279 def get_source_depth(self, args):
2280 return args[1]
2282 def get_receiver_depth(self, args):
2283 return args[0]
2285 def get_source_depths(self):
2286 return self.coords[1]
2288 def _update(self):
2289 self.mins = num.array([
2290 self.receiver_depth_min,
2291 self.source_depth_min,
2292 self.distance_min],
2293 dtype=float)
2295 self.maxs = num.array([
2296 self.receiver_depth_max,
2297 self.source_depth_max,
2298 self.distance_max],
2299 dtype=float)
2301 self.deltas = num.array([
2302 self.receiver_depth_delta,
2303 self.source_depth_delta,
2304 self.distance_delta],
2305 dtype=float)
2307 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2308 vicinity_eps).astype(int) + 1
2309 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2310 self.deltat = 1.0/self.sample_rate
2311 self.nrecords = num.product(self.ns) * self.ncomponents
2312 self.coords = tuple(num.linspace(mi, ma, n) for
2313 (mi, ma, n) in
2314 zip(self.mins, self.effective_maxs, self.ns)) + \
2315 (num.arange(self.ncomponents),)
2316 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2318 def _make_index_functions(self):
2320 amin, bmin, cmin = self.mins
2321 da, db, dc = self.deltas
2322 na, nb, nc = self.ns
2323 ng = self.ncomponents
2325 def index_function(a, b, c, ig):
2326 ia = int(round((a - amin) / da))
2327 ib = int(round((b - bmin) / db))
2328 ic = int(round((c - cmin) / dc))
2329 try:
2330 return num.ravel_multi_index((ia, ib, ic, ig),
2331 (na, nb, nc, ng))
2332 except ValueError:
2333 raise OutOfBounds()
2335 def indices_function(a, b, c, ig):
2336 ia = num.round((a - amin) / da).astype(int)
2337 ib = num.round((b - bmin) / db).astype(int)
2338 ic = num.round((c - cmin) / dc).astype(int)
2339 try:
2340 return num.ravel_multi_index((ia, ib, ic, ig),
2341 (na, nb, nc, ng))
2342 except ValueError:
2343 raise OutOfBounds()
2345 def grid_interpolation_coefficients(a, b, c):
2346 ias = indi12((a - amin) / da, na)
2347 ibs = indi12((b - bmin) / db, nb)
2348 ics = indi12((c - cmin) / dc, nc)
2349 return ias, ibs, ics
2351 def vicinity_function(a, b, c, ig):
2352 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2354 if not (0 <= ig < ng):
2355 raise OutOfBounds()
2357 indis = []
2358 weights = []
2359 for ia, va in ias:
2360 iia = ia*nb*nc*ng
2361 for ib, vb in ibs:
2362 iib = ib*nc*ng
2363 for ic, vc in ics:
2364 indis.append(iia + iib + ic*ng + ig)
2365 weights.append(va*vb*vc)
2367 return num.array(indis), num.array(weights)
2369 def vicinities_function(a, b, c, ig):
2371 xa = (a - amin) / da
2372 xb = (b - bmin) / db
2373 xc = (c - cmin) / dc
2375 xa_fl = num.floor(xa)
2376 xa_ce = num.ceil(xa)
2377 xb_fl = num.floor(xb)
2378 xb_ce = num.ceil(xb)
2379 xc_fl = num.floor(xc)
2380 xc_ce = num.ceil(xc)
2381 va_fl = 1.0 - (xa - xa_fl)
2382 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2383 vb_fl = 1.0 - (xb - xb_fl)
2384 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2385 vc_fl = 1.0 - (xc - xc_fl)
2386 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2388 ia_fl = xa_fl.astype(int)
2389 ia_ce = xa_ce.astype(int)
2390 ib_fl = xb_fl.astype(int)
2391 ib_ce = xb_ce.astype(int)
2392 ic_fl = xc_fl.astype(int)
2393 ic_ce = xc_ce.astype(int)
2395 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2396 raise OutOfBounds()
2398 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2399 raise OutOfBounds()
2401 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2402 raise OutOfBounds()
2404 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2405 raise OutOfBounds()
2407 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2408 raise OutOfBounds()
2410 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2411 raise OutOfBounds()
2413 irecords = num.empty(a.size*8, dtype=int)
2414 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2415 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2416 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2417 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2418 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2419 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2420 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2421 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2423 weights = num.empty(a.size*8, dtype=float)
2424 weights[0::8] = va_fl * vb_fl * vc_fl
2425 weights[1::8] = va_ce * vb_fl * vc_fl
2426 weights[2::8] = va_fl * vb_ce * vc_fl
2427 weights[3::8] = va_ce * vb_ce * vc_fl
2428 weights[4::8] = va_fl * vb_fl * vc_ce
2429 weights[5::8] = va_ce * vb_fl * vc_ce
2430 weights[6::8] = va_fl * vb_ce * vc_ce
2431 weights[7::8] = va_ce * vb_ce * vc_ce
2433 return irecords, weights
2435 self._index_function = index_function
2436 self._indices_function = indices_function
2437 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2438 self._vicinity_function = vicinity_function
2439 self._vicinities_function = vicinities_function
2441 def make_indexing_args(self, source, receiver, icomponents):
2442 nc = icomponents.size
2443 dists = source.distances_to(receiver)
2444 n = dists.size
2445 receiver_depths = num.empty(nc)
2446 receiver_depths.fill(receiver.depth)
2447 return (receiver_depths,
2448 num.tile(source.depths, nc//n),
2449 num.tile(dists, nc//n),
2450 icomponents)
2452 def make_indexing_args1(self, source, receiver):
2453 return (receiver.depth,
2454 source.depth,
2455 source.distance_to(receiver))
2457 @property
2458 def short_extent(self):
2459 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2460 self.receiver_depth_min/km,
2461 self.receiver_depth_max/km,
2462 self.receiver_depth_delta/km,
2463 self.source_depth_min/km,
2464 self.source_depth_max/km,
2465 self.source_depth_delta/km,
2466 self.distance_min/km,
2467 self.distance_max/km,
2468 self.distance_delta/km)
2470 def fix_ttt_holes(self, sptree, mode):
2471 from pyrocko import eikonal_ext, spit
2473 nodes_sr = self.nodes(minlevel=1, level=-1)
2475 delta = self.deltas[-1]
2476 assert num.all(delta == self.deltas[1:])
2478 nreceivers, nsources, ndistances = self.ns
2480 points = num.zeros((nodes_sr.shape[0], 3))
2481 points[:, 0] = nodes_sr[:, 1]
2482 points[:, 2] = nodes_sr[:, 0]
2484 speeds = self.get_material_property(
2485 0., 0., points,
2486 parameter='vp' if mode == cake.P else 'vs',
2487 interpolation='multilinear')
2489 speeds = speeds.reshape((nsources, ndistances))
2491 receiver_times = []
2492 for ireceiver in range(nreceivers):
2493 nodes = num.hstack([
2494 num_full(
2495 (nodes_sr.shape[0], 1),
2496 self.coords[0][ireceiver]),
2497 nodes_sr])
2499 times = sptree.interpolate_many(nodes)
2501 times[num.isnan(times)] = -1.
2503 times = times.reshape(speeds.shape)
2505 try:
2506 eikonal_ext.eikonal_solver_fmm_cartesian(
2507 speeds, times, delta)
2508 except eikonal_ext.EikonalExtError as e:
2509 if str(e).endswith('please check results'):
2510 logger.debug(
2511 'Got a warning from eikonal solver '
2512 '- may be ok...')
2513 else:
2514 raise
2516 receiver_times.append(times)
2518 def func(x):
2519 ias, ibs, ics = \
2520 self.grid_interpolation_coefficients(*x)
2522 t = 0
2523 for ia, va in ias:
2524 times = receiver_times[ia]
2525 for ib, vb in ibs:
2526 for ic, vc in ics:
2527 t += times[ib, ic] * va * vb * vc
2529 return t
2531 return spit.SPTree(
2532 f=func,
2533 ftol=sptree.ftol,
2534 xbounds=sptree.xbounds,
2535 xtols=sptree.xtols)
2538class ConfigTypeC(Config):
2539 '''
2540 No symmetrical constraints, one fixed receiver position.
2542 * Cartesian 3D source volume around a reference point
2544 * High level index variables: ``(source_depth,
2545 source_east_shift, source_north_shift, component)``
2546 '''
2548 receiver = Receiver.T(
2549 help='Receiver location')
2551 source_origin = Location.T(
2552 help='Origin of the source volume grid.')
2554 source_depth_min = Float.T(
2555 help='Minimum source depth [m].')
2557 source_depth_max = Float.T(
2558 help='Maximum source depth [m].')
2560 source_depth_delta = Float.T(
2561 help='Source depth grid spacing [m].')
2563 source_east_shift_min = Float.T(
2564 help='Minimum easting of source grid [m].')
2566 source_east_shift_max = Float.T(
2567 help='Maximum easting of source grid [m].')
2569 source_east_shift_delta = Float.T(
2570 help='Source volume grid spacing in east direction [m].')
2572 source_north_shift_min = Float.T(
2573 help='Minimum northing of source grid [m].')
2575 source_north_shift_max = Float.T(
2576 help='Maximum northing of source grid [m].')
2578 source_north_shift_delta = Float.T(
2579 help='Source volume grid spacing in north direction [m].')
2581 short_type = 'C'
2583 provided_schemes = ['elastic18']
2585 def get_surface_distance(self, args):
2586 _, source_east_shift, source_north_shift, _ = args
2587 sorig = self.source_origin
2588 sloc = Location(
2589 lat=sorig.lat,
2590 lon=sorig.lon,
2591 north_shift=sorig.north_shift + source_north_shift,
2592 east_shift=sorig.east_shift + source_east_shift)
2594 return self.receiver.distance_to(sloc)
2596 def get_distance(self, args):
2597 sdepth, source_east_shift, source_north_shift, _ = args
2598 sorig = self.source_origin
2599 sloc = Location(
2600 lat=sorig.lat,
2601 lon=sorig.lon,
2602 north_shift=sorig.north_shift + source_north_shift,
2603 east_shift=sorig.east_shift + source_east_shift)
2605 return self.receiver.distance_3d_to(sloc)
2607 def get_source_depth(self, args):
2608 return args[0]
2610 def get_receiver_depth(self, args):
2611 return self.receiver.depth
2613 def get_source_depths(self):
2614 return self.coords[0]
2616 def _update(self):
2617 self.mins = num.array([
2618 self.source_depth_min,
2619 self.source_east_shift_min,
2620 self.source_north_shift_min],
2621 dtype=float)
2623 self.maxs = num.array([
2624 self.source_depth_max,
2625 self.source_east_shift_max,
2626 self.source_north_shift_max],
2627 dtype=float)
2629 self.deltas = num.array([
2630 self.source_depth_delta,
2631 self.source_east_shift_delta,
2632 self.source_north_shift_delta],
2633 dtype=float)
2635 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2636 vicinity_eps).astype(int) + 1
2637 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2638 self.deltat = 1.0/self.sample_rate
2639 self.nrecords = num.product(self.ns) * self.ncomponents
2641 self.coords = tuple(
2642 num.linspace(mi, ma, n) for (mi, ma, n) in
2643 zip(self.mins, self.effective_maxs, self.ns)) + \
2644 (num.arange(self.ncomponents),)
2646 def _make_index_functions(self):
2648 amin, bmin, cmin = self.mins
2649 da, db, dc = self.deltas
2650 na, nb, nc = self.ns
2651 ng = self.ncomponents
2653 def index_function(a, b, c, ig):
2654 ia = int(round((a - amin) / da))
2655 ib = int(round((b - bmin) / db))
2656 ic = int(round((c - cmin) / dc))
2657 try:
2658 return num.ravel_multi_index((ia, ib, ic, ig),
2659 (na, nb, nc, ng))
2660 except ValueError:
2661 raise OutOfBounds(values=(a, b, c, ig))
2663 def indices_function(a, b, c, ig):
2664 ia = num.round((a - amin) / da).astype(int)
2665 ib = num.round((b - bmin) / db).astype(int)
2666 ic = num.round((c - cmin) / dc).astype(int)
2668 try:
2669 return num.ravel_multi_index((ia, ib, ic, ig),
2670 (na, nb, nc, ng))
2671 except ValueError:
2672 raise OutOfBounds()
2674 def vicinity_function(a, b, c, ig):
2675 ias = indi12((a - amin) / da, na)
2676 ibs = indi12((b - bmin) / db, nb)
2677 ics = indi12((c - cmin) / dc, nc)
2679 if not (0 <= ig < ng):
2680 raise OutOfBounds()
2682 indis = []
2683 weights = []
2684 for ia, va in ias:
2685 iia = ia*nb*nc*ng
2686 for ib, vb in ibs:
2687 iib = ib*nc*ng
2688 for ic, vc in ics:
2689 indis.append(iia + iib + ic*ng + ig)
2690 weights.append(va*vb*vc)
2692 return num.array(indis), num.array(weights)
2694 def vicinities_function(a, b, c, ig):
2696 xa = (a-amin) / da
2697 xb = (b-bmin) / db
2698 xc = (c-cmin) / dc
2700 xa_fl = num.floor(xa)
2701 xa_ce = num.ceil(xa)
2702 xb_fl = num.floor(xb)
2703 xb_ce = num.ceil(xb)
2704 xc_fl = num.floor(xc)
2705 xc_ce = num.ceil(xc)
2706 va_fl = 1.0 - (xa - xa_fl)
2707 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2708 vb_fl = 1.0 - (xb - xb_fl)
2709 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2710 vc_fl = 1.0 - (xc - xc_fl)
2711 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2713 ia_fl = xa_fl.astype(int)
2714 ia_ce = xa_ce.astype(int)
2715 ib_fl = xb_fl.astype(int)
2716 ib_ce = xb_ce.astype(int)
2717 ic_fl = xc_fl.astype(int)
2718 ic_ce = xc_ce.astype(int)
2720 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2721 raise OutOfBounds()
2723 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2724 raise OutOfBounds()
2726 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2727 raise OutOfBounds()
2729 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2730 raise OutOfBounds()
2732 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2733 raise OutOfBounds()
2735 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2736 raise OutOfBounds()
2738 irecords = num.empty(a.size*8, dtype=int)
2739 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2740 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2741 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2742 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2743 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2744 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2745 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2746 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2748 weights = num.empty(a.size*8, dtype=float)
2749 weights[0::8] = va_fl * vb_fl * vc_fl
2750 weights[1::8] = va_ce * vb_fl * vc_fl
2751 weights[2::8] = va_fl * vb_ce * vc_fl
2752 weights[3::8] = va_ce * vb_ce * vc_fl
2753 weights[4::8] = va_fl * vb_fl * vc_ce
2754 weights[5::8] = va_ce * vb_fl * vc_ce
2755 weights[6::8] = va_fl * vb_ce * vc_ce
2756 weights[7::8] = va_ce * vb_ce * vc_ce
2758 return irecords, weights
2760 self._index_function = index_function
2761 self._indices_function = indices_function
2762 self._vicinity_function = vicinity_function
2763 self._vicinities_function = vicinities_function
2765 def make_indexing_args(self, source, receiver, icomponents):
2766 nc = icomponents.size
2768 dists = source.distances_to(self.source_origin)
2769 azis, _ = source.azibazis_to(self.source_origin)
2771 source_north_shifts = - num.cos(d2r*azis) * dists
2772 source_east_shifts = - num.sin(d2r*azis) * dists
2773 source_depths = source.depths - self.source_origin.depth
2775 n = dists.size
2777 return (num.tile(source_depths, nc//n),
2778 num.tile(source_east_shifts, nc//n),
2779 num.tile(source_north_shifts, nc//n),
2780 icomponents)
2782 def make_indexing_args1(self, source, receiver):
2783 dist = source.distance_to(self.source_origin)
2784 azi, _ = source.azibazi_to(self.source_origin)
2786 source_north_shift = - num.cos(d2r*azi) * dist
2787 source_east_shift = - num.sin(d2r*azi) * dist
2788 source_depth = source.depth - self.source_origin.depth
2790 return (source_depth,
2791 source_east_shift,
2792 source_north_shift)
2794 @property
2795 def short_extent(self):
2796 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2797 self.source_depth_min/km,
2798 self.source_depth_max/km,
2799 self.source_depth_delta/km,
2800 self.source_east_shift_min/km,
2801 self.source_east_shift_max/km,
2802 self.source_east_shift_delta/km,
2803 self.source_north_shift_min/km,
2804 self.source_north_shift_max/km,
2805 self.source_north_shift_delta/km)
2808class Weighting(Object):
2809 factor = Float.T(default=1.0)
2812class Taper(Object):
2813 tmin = Timing.T()
2814 tmax = Timing.T()
2815 tfade = Float.T(default=0.0)
2816 shape = StringChoice.T(
2817 choices=['cos', 'linear'],
2818 default='cos',
2819 optional=True)
2822class SimplePattern(SObject):
2824 _pool = {}
2826 def __init__(self, pattern):
2827 self._pattern = pattern
2828 SObject.__init__(self)
2830 def __str__(self):
2831 return self._pattern
2833 @property
2834 def regex(self):
2835 pool = SimplePattern._pool
2836 if self.pattern not in pool:
2837 rpat = '|'.join(fnmatch.translate(x) for
2838 x in self.pattern.split('|'))
2839 pool[self.pattern] = re.compile(rpat, re.I)
2841 return pool[self.pattern]
2843 def match(self, s):
2844 return self.regex.match(s)
2847class WaveformType(StringChoice):
2848 choices = ['dis', 'vel', 'acc',
2849 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2850 'envelope_dis', 'envelope_vel', 'envelope_acc']
2853class ChannelSelection(Object):
2854 pattern = SimplePattern.T(optional=True)
2855 min_sample_rate = Float.T(optional=True)
2856 max_sample_rate = Float.T(optional=True)
2859class StationSelection(Object):
2860 includes = SimplePattern.T()
2861 excludes = SimplePattern.T()
2862 distance_min = Float.T(optional=True)
2863 distance_max = Float.T(optional=True)
2864 azimuth_min = Float.T(optional=True)
2865 azimuth_max = Float.T(optional=True)
2868class WaveformSelection(Object):
2869 channel_selection = ChannelSelection.T(optional=True)
2870 station_selection = StationSelection.T(optional=True)
2871 taper = Taper.T()
2872 # filter = FrequencyResponse.T()
2873 waveform_type = WaveformType.T(default='dis')
2874 weighting = Weighting.T(optional=True)
2875 sample_rate = Float.T(optional=True)
2876 gf_store_id = StringID.T(optional=True)
2879def indi12(x, n):
2880 '''
2881 Get linear interpolation index and weight.
2882 '''
2884 r = round(x)
2885 if abs(r - x) < vicinity_eps:
2886 i = int(r)
2887 if not (0 <= i < n):
2888 raise OutOfBounds()
2890 return ((int(r), 1.),)
2891 else:
2892 f = math.floor(x)
2893 i = int(f)
2894 if not (0 <= i < n-1):
2895 raise OutOfBounds()
2897 v = x-f
2898 return ((i, 1.-v), (i + 1, v))
2901def float_or_none(s):
2902 units = {
2903 'k': 1e3,
2904 'M': 1e6,
2905 }
2907 factor = 1.0
2908 if s and s[-1] in units:
2909 factor = units[s[-1]]
2910 s = s[:-1]
2911 if not s:
2912 raise ValueError('unit without a number: \'%s\'' % s)
2914 if s:
2915 return float(s) * factor
2916 else:
2917 return None
2920class GridSpecError(Exception):
2921 def __init__(self, s):
2922 Exception.__init__(self, 'invalid grid specification: %s' % s)
2925def parse_grid_spec(spec):
2926 try:
2927 result = []
2928 for dspec in spec.split(','):
2929 t = dspec.split('@')
2930 num = start = stop = step = None
2931 if len(t) == 2:
2932 num = int(t[1])
2933 if num <= 0:
2934 raise GridSpecError(spec)
2936 elif len(t) > 2:
2937 raise GridSpecError(spec)
2939 s = t[0]
2940 v = [float_or_none(x) for x in s.split(':')]
2941 if len(v) == 1:
2942 start = stop = v[0]
2943 if len(v) >= 2:
2944 start, stop = v[0:2]
2945 if len(v) == 3:
2946 step = v[2]
2948 if len(v) > 3 or (len(v) > 2 and num is not None):
2949 raise GridSpecError(spec)
2951 if step == 0.0:
2952 raise GridSpecError(spec)
2954 result.append((start, stop, step, num))
2956 except ValueError:
2957 raise GridSpecError(spec)
2959 return result
2962def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2963 swap = step is not None and step < 0.
2964 if start is None:
2965 start = ma if swap else mi
2966 if stop is None:
2967 stop = mi if swap else ma
2968 if step is None:
2969 step = -inc if ma < mi else inc
2970 if num is None:
2971 if (step < 0) != (stop-start < 0):
2972 raise GridSpecError()
2974 num = int(round((stop-start)/step))+1
2975 stop2 = start + (num-1)*step
2976 if abs(stop-stop2) > eps:
2977 num = int(math.floor((stop-start)/step))+1
2978 stop = start + (num-1)*step
2979 else:
2980 stop = stop2
2982 if start == stop:
2983 num = 1
2985 return start, stop, num
2988def nditer_outer(x):
2989 return num.nditer(
2990 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
2993def nodes(xs):
2994 ns = [x.size for x in xs]
2995 nnodes = num.prod(ns)
2996 ndim = len(xs)
2997 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
2998 for idim in range(ndim-1, -1, -1):
2999 x = xs[idim]
3000 nrepeat = num.prod(ns[idim+1:], dtype=int)
3001 ntile = num.prod(ns[:idim], dtype=int)
3002 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
3004 return nodes
3007def filledi(x, n):
3008 a = num.empty(n, dtype=int)
3009 a.fill(x)
3010 return a
3013config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3015discretized_source_classes = [
3016 DiscretizedExplosionSource,
3017 DiscretizedSFSource,
3018 DiscretizedMTSource,
3019 DiscretizedPorePressureSource]
3022__all__ = '''
3023Earthmodel1D
3024StringID
3025ScopeType
3026WaveformType
3027QuantityType
3028NearfieldTermsType
3029Reference
3030Region
3031CircularRegion
3032RectangularRegion
3033PhaseSelect
3034InvalidTimingSpecification
3035Timing
3036TPDef
3037OutOfBounds
3038Location
3039Receiver
3040'''.split() + [
3041 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3042ComponentScheme
3043component_scheme_to_description
3044component_schemes
3045Config
3046GridSpecError
3047Weighting
3048Taper
3049SimplePattern
3050WaveformType
3051ChannelSelection
3052StationSelection
3053WaveformSelection
3054nditer_outer
3055dump
3056load
3057discretized_source_classes
3058config_type_classes
3059UnavailableScheme
3060InterpolationMethod
3061SeismosizerTrace
3062SeismosizerResult
3063Result
3064StaticResult
3065'''.split()