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='elastic2',
160 source_terms=['m0'],
161 ncomponents=2,
162 default_stored_quantity='displacement',
163 provided_components=[
164 'n', 'e', 'd']),
165 ComponentSchemeDescription(
166 name='elastic8',
167 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
168 ncomponents=8,
169 default_stored_quantity='displacement',
170 provided_components=[
171 'n', 'e', 'd']),
172 ComponentSchemeDescription(
173 name='elastic10',
174 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
175 ncomponents=10,
176 default_stored_quantity='displacement',
177 provided_components=[
178 'n', 'e', 'd']),
179 ComponentSchemeDescription(
180 name='elastic18',
181 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
182 ncomponents=18,
183 default_stored_quantity='displacement',
184 provided_components=[
185 'n', 'e', 'd']),
186 ComponentSchemeDescription(
187 name='elastic5',
188 source_terms=['fn', 'fe', 'fd'],
189 ncomponents=5,
190 default_stored_quantity='displacement',
191 provided_components=[
192 'n', 'e', 'd']),
193 ComponentSchemeDescription(
194 name='poroelastic10',
195 source_terms=['pore_pressure'],
196 ncomponents=10,
197 default_stored_quantity=None,
198 provided_components=[
199 'displacement.n', 'displacement.e', 'displacement.d',
200 'vertical_tilt.n', 'vertical_tilt.e',
201 'pore_pressure',
202 'darcy_velocity.n', 'darcy_velocity.e', 'darcy_velocity.d'])]
205# new names?
206# 'mt_to_displacement_1d'
207# 'mt_to_displacement_1d_ff_only'
208# 'mt_to_gravity_1d'
209# 'mt_to_stress_1d'
210# 'explosion_to_displacement_1d'
211# 'sf_to_displacement_1d'
212# 'mt_to_displacement_3d'
213# 'mt_to_gravity_3d'
214# 'mt_to_stress_3d'
215# 'pore_pressure_to_displacement_1d'
216# 'pore_pressure_to_vertical_tilt_1d'
217# 'pore_pressure_to_pore_pressure_1d'
218# 'pore_pressure_to_darcy_velocity_1d'
221component_schemes = [c.name for c in component_scheme_descriptions]
222component_scheme_to_description = dict(
223 (c.name, c) for c in component_scheme_descriptions)
226class ComponentScheme(StringChoice):
227 '''
228 Different Green's Function component schemes are available:
230 ================= =========================================================
231 Name Description
232 ================= =========================================================
233 ``elastic10`` Elastodynamic for
234 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
235 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT
236 sources only
237 ``elastic8`` Elastodynamic for far-field only
238 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
239 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores,
240 MT sources only
241 ``elastic2`` Elastodynamic for
242 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
243 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, purely
244 isotropic sources only
245 ``elastic5`` Elastodynamic for
246 :py:class:`~pyrocko.gf.meta.ConfigTypeA`
247 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, SF
248 sources only
249 ``elastic18`` Elastodynamic for
250 :py:class:`~pyrocko.gf.meta.ConfigTypeC` stores, MT
251 sources only
252 ``poroelastic10`` Poroelastic for :py:class:`~pyrocko.gf.meta.ConfigTypeA`
253 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores
254 ================= =========================================================
255 '''
257 choices = component_schemes
260class Earthmodel1D(Object):
261 dummy_for = cake.LayeredModel
263 class __T(TBase):
264 def regularize_extra(self, val):
265 if isinstance(val, str):
266 val = cake.LayeredModel.from_scanlines(
267 cake.read_nd_model_str(val))
269 return val
271 def to_save(self, val):
272 return literal(cake.write_nd_model_str(val))
275class StringID(StringPattern):
276 pattern = r'^[A-Za-z][A-Za-z0-9._]{0,64}$'
279class ScopeType(StringChoice):
280 choices = [
281 'global',
282 'regional',
283 'local',
284 ]
287class WaveType(StringChoice):
288 choices = [
289 'full waveform',
290 'bodywave',
291 'P wave',
292 'S wave',
293 'surface wave',
294 ]
297class NearfieldTermsType(StringChoice):
298 choices = [
299 'complete',
300 'incomplete',
301 'missing',
302 ]
305class QuantityType(StringChoice):
306 choices = [
307 'displacement',
308 'rotation',
309 'velocity',
310 'acceleration',
311 'pressure',
312 'tilt',
313 'pore_pressure',
314 'darcy_velocity',
315 'vertical_tilt']
318class Reference(Object):
319 id = StringID.T()
320 type = String.T()
321 title = Unicode.T()
322 journal = Unicode.T(optional=True)
323 volume = Unicode.T(optional=True)
324 number = Unicode.T(optional=True)
325 pages = Unicode.T(optional=True)
326 year = String.T()
327 note = Unicode.T(optional=True)
328 issn = String.T(optional=True)
329 doi = String.T(optional=True)
330 url = String.T(optional=True)
331 eprint = String.T(optional=True)
332 authors = List.T(Unicode.T())
333 publisher = Unicode.T(optional=True)
334 keywords = Unicode.T(optional=True)
335 note = Unicode.T(optional=True)
336 abstract = Unicode.T(optional=True)
338 @classmethod
339 def from_bibtex(cls, filename=None, stream=None):
341 from pybtex.database.input import bibtex
343 parser = bibtex.Parser()
345 if filename is not None:
346 bib_data = parser.parse_file(filename)
347 elif stream is not None:
348 bib_data = parser.parse_stream(stream)
350 references = []
352 for id_, entry in bib_data.entries.items():
353 d = {}
354 avail = entry.fields.keys()
355 for prop in cls.T.properties:
356 if prop.name == 'authors' or prop.name not in avail:
357 continue
359 d[prop.name] = entry.fields[prop.name]
361 if 'author' in entry.persons:
362 d['authors'] = []
363 for person in entry.persons['author']:
364 d['authors'].append(new_str(person))
366 c = Reference(id=id_, type=entry.type, **d)
367 references.append(c)
369 return references
372_fpat = r'[+-]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?'
373_spat = StringID.pattern[1:-1]
374_cake_pat = cake.PhaseDef.allowed_characters_pattern
375_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic
377_ppat = r'(' \
378 r'cake:' + _cake_pat + \
379 r'|iaspei:' + _iaspei_pat + \
380 r'|vel_surface:' + _fpat + \
381 r'|vel:' + _fpat + \
382 r'|stored:' + _spat + \
383 r')'
386timing_regex_old = re.compile(
387 r'^((first|last)?\((' + _spat + r'(\|' + _spat + r')*)\)|(' +
388 _spat + r'))?(' + _fpat + ')?$')
391timing_regex = re.compile(
392 r'^((first|last)?\{(' + _ppat + r'(\|' + _ppat + r')*)\}|(' +
393 _ppat + r'))?(' + _fpat + '(S|%)?)?$')
396class PhaseSelect(StringChoice):
397 choices = ['', 'first', 'last']
400class InvalidTimingSpecification(ValidationError):
401 pass
404class Timing(SObject):
405 '''
406 Definition of a time instant relative to one or more named phase arrivals.
408 Instances of this class can be used e.g. in cutting and tapering
409 operations. They can hold an absolute time or an offset to a named phase
410 arrival or group of such arrivals.
412 Timings can be instantiated from a simple string defintion i.e. with
413 ``Timing(str)`` where ``str`` is something like
414 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``,
415 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase
416 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'``
417 is appended, it is interpreted as percent. If the an ``'S'`` is appended
418 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km].
420 Phase definitions can be specified in either of the following ways:
422 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table
423 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake
424 (see :py:class:`pyrocko.cake.PhaseDef`)
425 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance /
426 velocity [km/s]
427 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s]
429 **Examples:**
431 * ``'100'`` : absolute time; 100 s
432 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to
433 stored travel time table named ``'P'``
434 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to
435 stored travel time table named ``'P'``
436 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to
437 stored travel time table named ``'P'``
438 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is
439 undefined for a given geometry
440 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A
441 and B is chosen, if both phases are defined for a given geometry
442 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen
443 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of
444 arrivals A, B, and C
445 '''
447 def __init__(self, s=None, **kwargs):
449 if s is not None:
450 offset_is = None
451 s = re.sub(r'\s+', '', s)
452 try:
453 offset = float(s.rstrip('S'))
455 if s.endswith('S'):
456 offset_is = 'slowness'
458 phase_defs = []
459 select = ''
461 except ValueError:
462 matched = False
463 m = timing_regex.match(s)
464 if m:
465 if m.group(3):
466 phase_defs = m.group(3).split('|')
467 elif m.group(19):
468 phase_defs = [m.group(19)]
469 else:
470 phase_defs = []
472 select = m.group(2) or ''
474 offset = 0.0
475 soff = m.group(27)
476 if soff:
477 offset = float(soff.rstrip('S%'))
478 if soff.endswith('S'):
479 offset_is = 'slowness'
480 elif soff.endswith('%'):
481 offset_is = 'percent'
483 matched = True
485 else:
486 m = timing_regex_old.match(s)
487 if m:
488 if m.group(3):
489 phase_defs = m.group(3).split('|')
490 elif m.group(5):
491 phase_defs = [m.group(5)]
492 else:
493 phase_defs = []
495 select = m.group(2) or ''
497 offset = 0.0
498 if m.group(6):
499 offset = float(m.group(6))
501 phase_defs = [
502 'stored:' + phase_def for phase_def in phase_defs]
504 matched = True
506 if not matched:
507 raise InvalidTimingSpecification(s)
509 kwargs = dict(
510 phase_defs=phase_defs,
511 select=select,
512 offset=offset,
513 offset_is=offset_is)
515 SObject.__init__(self, **kwargs)
517 def __str__(self):
518 s = []
519 if self.phase_defs:
520 sphases = '|'.join(self.phase_defs)
521 # if len(self.phase_defs) > 1 or self.select:
522 sphases = '{%s}' % sphases
524 if self.select:
525 sphases = self.select + sphases
527 s.append(sphases)
529 if self.offset != 0.0 or not self.phase_defs:
530 s.append('%+g' % self.offset)
531 if self.offset_is == 'slowness':
532 s.append('S')
533 elif self.offset_is == 'percent':
534 s.append('%')
536 return ''.join(s)
538 def evaluate(self, get_phase, args):
539 try:
540 if self.offset_is == 'slowness' and self.offset != 0.0:
541 phase_offset = get_phase(
542 'vel_surface:%g' % (1.0/self.offset))
543 offset = phase_offset(args)
544 else:
545 offset = self.offset
547 if self.phase_defs:
548 phases = [
549 get_phase(phase_def) for phase_def in self.phase_defs]
550 times = [phase(args) for phase in phases]
551 if self.offset_is == 'percent':
552 times = [t*(1.+offset/100.)
553 for t in times if t is not None]
554 else:
555 times = [t+offset for t in times if t is not None]
557 if not times:
558 return None
559 elif self.select == 'first':
560 return min(times)
561 elif self.select == 'last':
562 return max(times)
563 else:
564 return times[0]
565 else:
566 return offset
568 except spit.OutOfBounds:
569 raise OutOfBounds(args)
571 phase_defs = List.T(String.T())
572 offset = Float.T(default=0.0)
573 offset_is = String.T(optional=True)
574 select = PhaseSelect.T(
575 default='',
576 help=('Can be either ``\'%s\'``, ``\'%s\'``, or ``\'%s\'``. ' %
577 tuple(PhaseSelect.choices)))
580def mkdefs(s):
581 defs = []
582 for x in s.split(','):
583 try:
584 defs.append(float(x))
585 except ValueError:
586 if x.startswith('!'):
587 defs.extend(cake.PhaseDef.classic(x[1:]))
588 else:
589 defs.append(cake.PhaseDef(x))
591 return defs
594class TPDef(Object):
595 '''
596 Maps an arrival phase identifier to an arrival phase definition.
597 '''
599 id = StringID.T(
600 help='name used to identify the phase')
601 definition = String.T(
602 help='definition of the phase in either cake syntax as defined in '
603 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an '
604 '``!``, as a *classic phase name*, or, if it is a simple '
605 'number, as a constant horizontal velocity.')
607 @property
608 def phases(self):
609 return [x for x in mkdefs(self.definition)
610 if isinstance(x, cake.PhaseDef)]
612 @property
613 def horizontal_velocities(self):
614 return [x for x in mkdefs(self.definition) if isinstance(x, float)]
617class OutOfBounds(Exception):
618 def __init__(self, values=None, reason=None):
619 Exception.__init__(self)
620 self.values = values
621 self.reason = reason
622 self.context = None
624 def __str__(self):
625 scontext = ''
626 if self.context:
627 scontext = '\n%s' % str(self.context)
629 if self.reason:
630 scontext += '\n%s' % self.reason
632 if self.values:
633 return 'out of bounds: (%s)%s' % (
634 ', '.join('%g' % x for x in self.values), scontext)
635 else:
636 return 'out of bounds%s' % scontext
639class MultiLocation(Object):
641 lats = Array.T(
642 optional=True, shape=(None,), dtype=float,
643 help='Latitudes of targets.')
644 lons = Array.T(
645 optional=True, shape=(None,), dtype=float,
646 help='Longitude of targets.')
647 north_shifts = Array.T(
648 optional=True, shape=(None,), dtype=float,
649 help='North shifts of targets.')
650 east_shifts = Array.T(
651 optional=True, shape=(None,), dtype=float,
652 help='East shifts of targets.')
653 elevation = Array.T(
654 optional=True, shape=(None,), dtype=float,
655 help='Elevations of targets.')
657 def __init__(self, *args, **kwargs):
658 self._coords5 = None
659 Object.__init__(self, *args, **kwargs)
661 @property
662 def coords5(self):
663 if self._coords5 is not None:
664 return self._coords5
665 props = [self.lats, self.lons, self.north_shifts, self.east_shifts,
666 self.elevation]
667 sizes = [p.size for p in props if p is not None]
668 if not sizes:
669 raise AttributeError('Empty StaticTarget')
670 if num.unique(sizes).size != 1:
671 raise AttributeError('Inconsistent coordinate sizes.')
673 self._coords5 = num.zeros((sizes[0], 5))
674 for idx, p in enumerate(props):
675 if p is not None:
676 self._coords5[:, idx] = p.flatten()
678 return self._coords5
680 @property
681 def ncoords(self):
682 if self.coords5.shape[0] is None:
683 return 0
684 return int(self.coords5.shape[0])
686 def get_latlon(self):
687 '''
688 Get all coordinates as lat/lon.
690 :returns: Coordinates in Latitude, Longitude
691 :rtype: :class:`numpy.ndarray`, (N, 2)
692 '''
693 latlons = num.empty((self.ncoords, 2))
694 for i in range(self.ncoords):
695 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4])
696 return latlons
698 def get_corner_coords(self):
699 '''
700 Returns the corner coordinates of the multi-location object.
702 :returns: In lat/lon: lower left, upper left, upper right, lower right
703 :rtype: tuple
704 '''
705 latlon = self.get_latlon()
706 ll = (latlon[:, 0].min(), latlon[:, 1].min())
707 ur = (latlon[:, 0].max(), latlon[:, 1].max())
708 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1]))
711class Receiver(Location):
712 codes = Tuple.T(
713 3, String.T(),
714 optional=True,
715 help='network, station, and location codes')
717 def pyrocko_station(self):
718 from pyrocko import model
720 lat, lon = self.effective_latlon
722 return model.Station(
723 network=self.codes[0],
724 station=self.codes[1],
725 location=self.codes[2],
726 lat=lat,
727 lon=lon,
728 depth=self.depth)
731def g(x, d):
732 if x is None:
733 return d
734 else:
735 return x
738class UnavailableScheme(Exception):
739 pass
742class InvalidNComponents(Exception):
743 pass
746class DiscretizedSource(Object):
747 '''
748 Base class for discretized sources.
750 To compute synthetic seismograms, the parameterized source models (see of
751 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first
752 discretized into a number of point sources. These spacio-temporal point
753 source distributions are represented by subclasses of the
754 :py:class:`DiscretizedSource`. For elastodynamic problems there is the
755 :py:class:`DiscretizedMTSource` for moment tensor point source
756 distributions and the :py:class:`DiscretizedExplosionSource` for pure
757 explosion/implosion type source distributions. The geometry calculations
758 are implemented in the base class. How Green's function components have to
759 be weighted and sumed is defined in the derived classes.
761 Like in the :py:class:`Location` class, the positions of the point sources
762 contained in the discretized source are defined by a common reference point
763 (:py:attr:`lat`, :py:attr:`lon`) and cartesian offsets to this
764 (:py:attr:`north_shifts`, :py:attr:`east_shifts`, :py:attr:`depths`).
765 Alternatively latitude and longitude of each contained point source can be
766 specified directly (:py:attr:`lats`, :py:attr:`lons`).
767 '''
769 times = Array.T(shape=(None,), dtype=float)
770 lats = Array.T(shape=(None,), dtype=float, optional=True)
771 lons = Array.T(shape=(None,), dtype=float, optional=True)
772 lat = Float.T(optional=True)
773 lon = Float.T(optional=True)
774 north_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
775 east_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
776 depths = Array.T(shape=(None,), dtype=num.float64)
777 dl = Float.T(optional=True)
778 dw = Float.T(optional=True)
779 nl = Float.T(optional=True)
780 nw = Float.T(optional=True)
782 @classmethod
783 def check_scheme(cls, scheme):
784 '''
785 Check if given GF component scheme is supported by the class.
787 Raises :py:class:`UnavailableScheme` if the given scheme is not
788 supported by this discretized source class.
789 '''
791 if scheme not in cls.provided_schemes:
792 raise UnavailableScheme(
793 'source type "%s" does not support GF component scheme "%s"' %
794 (cls.__name__, scheme))
796 def __init__(self, **kwargs):
797 Object.__init__(self, **kwargs)
798 self._latlons = None
800 def __setattr__(self, name, value):
801 if name in ('lat', 'lon', 'north_shifts', 'east_shifts',
802 'lats', 'lons'):
803 self.__dict__['_latlons'] = None
805 Object.__setattr__(self, name, value)
807 def get_source_terms(self, scheme):
808 raise NotImplementedError()
810 def make_weights(self, receiver, scheme):
811 raise NotImplementedError()
813 @property
814 def effective_latlons(self):
815 '''
816 Property holding the offest-corrected lats and lons of all points.
817 '''
819 if self._latlons is None:
820 if self.lats is not None and self.lons is not None:
821 if (self.north_shifts is not None and
822 self.east_shifts is not None):
823 self._latlons = orthodrome.ne_to_latlon(
824 self.lats, self.lons,
825 self.north_shifts, self.east_shifts)
826 else:
827 self._latlons = self.lats, self.lons
828 else:
829 lat = g(self.lat, 0.0)
830 lon = g(self.lon, 0.0)
831 self._latlons = orthodrome.ne_to_latlon(
832 lat, lon, self.north_shifts, self.east_shifts)
834 return self._latlons
836 @property
837 def effective_north_shifts(self):
838 if self.north_shifts is not None:
839 return self.north_shifts
840 else:
841 return num.zeros(self.times.size)
843 @property
844 def effective_east_shifts(self):
845 if self.east_shifts is not None:
846 return self.east_shifts
847 else:
848 return num.zeros(self.times.size)
850 def same_origin(self, receiver):
851 '''
852 Check if receiver has same reference point.
853 '''
855 return (g(self.lat, 0.0) == receiver.lat and
856 g(self.lon, 0.0) == receiver.lon and
857 self.lats is None and self.lons is None)
859 def azibazis_to(self, receiver):
860 '''
861 Compute azimuths and backaziumuths to/from receiver, for all contained
862 points.
863 '''
865 if self.same_origin(receiver):
866 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts,
867 receiver.north_shift - self.north_shifts)
868 bazis = azis + 180.
869 else:
870 slats, slons = self.effective_latlons
871 rlat, rlon = receiver.effective_latlon
872 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon)
873 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons)
875 return azis, bazis
877 def distances_to(self, receiver):
878 '''
879 Compute distances to receiver for all contained points.
880 '''
881 if self.same_origin(receiver):
882 return num.sqrt((self.north_shifts - receiver.north_shift)**2 +
883 (self.east_shifts - receiver.east_shift)**2)
885 else:
886 slats, slons = self.effective_latlons
887 rlat, rlon = receiver.effective_latlon
888 return orthodrome.distance_accurate50m_numpy(slats, slons,
889 rlat, rlon)
891 def element_coords(self, i):
892 if self.lats is not None and self.lons is not None:
893 lat = float(self.lats[i])
894 lon = float(self.lons[i])
895 else:
896 lat = self.lat
897 lon = self.lon
899 if self.north_shifts is not None and self.east_shifts is not None:
900 north_shift = float(self.north_shifts[i])
901 east_shift = float(self.east_shifts[i])
903 else:
904 north_shift = east_shift = 0.0
906 return lat, lon, north_shift, east_shift
908 def coords5(self):
909 xs = num.zeros((self.nelements, 5))
911 if self.lats is not None and self.lons is not None:
912 xs[:, 0] = self.lats
913 xs[:, 1] = self.lons
914 else:
915 xs[:, 0] = g(self.lat, 0.0)
916 xs[:, 1] = g(self.lon, 0.0)
918 if self.north_shifts is not None and self.east_shifts is not None:
919 xs[:, 2] = self.north_shifts
920 xs[:, 3] = self.east_shifts
922 xs[:, 4] = self.depths
924 return xs
926 @property
927 def nelements(self):
928 return self.times.size
930 @classmethod
931 def combine(cls, sources, **kwargs):
932 '''
933 Combine several discretized source models.
935 Concatenenates all point sources in the given discretized ``sources``.
936 Care must be taken when using this function that the external amplitude
937 factors and reference times of the parameterized (undiscretized)
938 sources match or are accounted for.
939 '''
941 first = sources[0]
943 if not all(type(s) == type(first) for s in sources):
944 raise Exception('DiscretizedSource.combine must be called with '
945 'sources of same type.')
947 latlons = []
948 for s in sources:
949 latlons.append(s.effective_latlons)
951 lats, lons = num.hstack(latlons)
953 if all((s.lats is None and s.lons is None) for s in sources):
954 rlats = num.array([s.lat for s in sources], dtype=float)
955 rlons = num.array([s.lon for s in sources], dtype=float)
956 same_ref = num.all(
957 rlats == rlats[0]) and num.all(rlons == rlons[0])
958 else:
959 same_ref = False
961 cat = num.concatenate
962 times = cat([s.times for s in sources])
963 depths = cat([s.depths for s in sources])
965 if same_ref:
966 lat = first.lat
967 lon = first.lon
968 north_shifts = cat([s.effective_north_shifts for s in sources])
969 east_shifts = cat([s.effective_east_shifts for s in sources])
970 lats = None
971 lons = None
972 else:
973 lat = None
974 lon = None
975 north_shifts = None
976 east_shifts = None
978 return cls(
979 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
980 north_shifts=north_shifts, east_shifts=east_shifts,
981 depths=depths, **kwargs)
983 def centroid_position(self):
984 moments = self.moments()
985 norm = num.sum(moments)
986 if norm != 0.0:
987 w = moments / num.sum(moments)
988 else:
989 w = num.ones(moments.size)
991 if self.lats is not None and self.lons is not None:
992 lats, lons = self.effective_latlons
993 rlat, rlon = num.mean(lats), num.mean(lons)
994 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
995 else:
996 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
997 n, e = self.north_shifts, self.east_shifts
999 cn = num.sum(n*w)
1000 ce = num.sum(e*w)
1001 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1003 if self.lats is not None and self.lons is not None:
1004 lat = clat
1005 lon = clon
1006 north_shift = 0.
1007 east_shift = 0.
1008 else:
1009 lat = g(self.lat, 0.0)
1010 lon = g(self.lon, 0.0)
1011 north_shift = cn
1012 east_shift = ce
1014 depth = num.sum(self.depths*w)
1015 time = num.sum(self.times*w)
1016 return tuple(float(x) for x in
1017 (time, lat, lon, north_shift, east_shift, depth))
1020class DiscretizedExplosionSource(DiscretizedSource):
1021 m0s = Array.T(shape=(None,), dtype=float)
1023 provided_schemes = (
1024 'elastic2',
1025 'elastic8',
1026 'elastic10',
1027 )
1029 def get_source_terms(self, scheme):
1030 self.check_scheme(scheme)
1032 if scheme == 'elastic2':
1033 return self.m0s[:, num.newaxis].copy()
1035 elif scheme in ('elastic8', 'elastic10'):
1036 m6s = num.zeros((self.m0s.size, 6))
1037 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1038 return m6s
1039 else:
1040 assert False
1042 def make_weights(self, receiver, scheme):
1043 self.check_scheme(scheme)
1045 azis, bazis = self.azibazis_to(receiver)
1047 sb = num.sin(bazis*d2r-num.pi)
1048 cb = num.cos(bazis*d2r-num.pi)
1050 m0s = self.m0s
1051 n = azis.size
1053 cat = num.concatenate
1054 rep = num.repeat
1056 if scheme == 'elastic2':
1057 w_n = cb*m0s
1058 g_n = filledi(0, n)
1059 w_e = sb*m0s
1060 g_e = filledi(0, n)
1061 w_d = m0s
1062 g_d = filledi(1, n)
1064 elif scheme == 'elastic8':
1065 w_n = cat((cb*m0s, cb*m0s))
1066 g_n = rep((0, 2), n)
1067 w_e = cat((sb*m0s, sb*m0s))
1068 g_e = rep((0, 2), n)
1069 w_d = cat((m0s, m0s))
1070 g_d = rep((5, 7), n)
1072 elif scheme == 'elastic10':
1073 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1074 g_n = rep((0, 2, 8), n)
1075 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1076 g_e = rep((0, 2, 8), n)
1077 w_d = cat((m0s, m0s, m0s))
1078 g_d = rep((5, 7, 9), n)
1080 else:
1081 assert False
1083 return (
1084 ('displacement.n', w_n, g_n),
1085 ('displacement.e', w_e, g_e),
1086 ('displacement.d', w_d, g_d),
1087 )
1089 def split(self):
1090 from pyrocko.gf.seismosizer import ExplosionSource
1091 sources = []
1092 for i in range(self.nelements):
1093 lat, lon, north_shift, east_shift = self.element_coords(i)
1094 sources.append(ExplosionSource(
1095 time=float(self.times[i]),
1096 lat=lat,
1097 lon=lon,
1098 north_shift=north_shift,
1099 east_shift=east_shift,
1100 depth=float(self.depths[i]),
1101 moment=float(self.m0s[i])))
1103 return sources
1105 def moments(self):
1106 return self.m0s
1108 def centroid(self):
1109 from pyrocko.gf.seismosizer import ExplosionSource
1110 time, lat, lon, north_shift, east_shift, depth = \
1111 self.centroid_position()
1113 return ExplosionSource(
1114 time=time,
1115 lat=lat,
1116 lon=lon,
1117 north_shift=north_shift,
1118 east_shift=east_shift,
1119 depth=depth,
1120 moment=float(num.sum(self.m0s)))
1122 @classmethod
1123 def combine(cls, sources, **kwargs):
1124 '''
1125 Combine several discretized source models.
1127 Concatenenates all point sources in the given discretized ``sources``.
1128 Care must be taken when using this function that the external amplitude
1129 factors and reference times of the parameterized (undiscretized)
1130 sources match or are accounted for.
1131 '''
1133 if 'm0s' not in kwargs:
1134 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1136 return super(DiscretizedExplosionSource, cls).combine(sources,
1137 **kwargs)
1140class DiscretizedSFSource(DiscretizedSource):
1141 forces = Array.T(shape=(None, 3), dtype=float)
1143 provided_schemes = (
1144 'elastic5',
1145 )
1147 def get_source_terms(self, scheme):
1148 self.check_scheme(scheme)
1150 return self.forces
1152 def make_weights(self, receiver, scheme):
1153 self.check_scheme(scheme)
1155 azis, bazis = self.azibazis_to(receiver)
1157 sa = num.sin(azis*d2r)
1158 ca = num.cos(azis*d2r)
1159 sb = num.sin(bazis*d2r-num.pi)
1160 cb = num.cos(bazis*d2r-num.pi)
1162 forces = self.forces
1163 fn = forces[:, 0]
1164 fe = forces[:, 1]
1165 fd = forces[:, 2]
1167 f0 = fd
1168 f1 = ca * fn + sa * fe
1169 f2 = ca * fe - sa * fn
1171 n = azis.size
1173 if scheme == 'elastic5':
1174 ioff = 0
1176 cat = num.concatenate
1177 rep = num.repeat
1179 w_n = cat((cb*f0, cb*f1, -sb*f2))
1180 g_n = ioff + rep((0, 1, 2), n)
1181 w_e = cat((sb*f0, sb*f1, cb*f2))
1182 g_e = ioff + rep((0, 1, 2), n)
1183 w_d = cat((f0, f1))
1184 g_d = ioff + rep((3, 4), n)
1186 return (
1187 ('displacement.n', w_n, g_n),
1188 ('displacement.e', w_e, g_e),
1189 ('displacement.d', w_d, g_d),
1190 )
1192 @classmethod
1193 def combine(cls, sources, **kwargs):
1194 '''
1195 Combine several discretized source models.
1197 Concatenenates all point sources in the given discretized ``sources``.
1198 Care must be taken when using this function that the external amplitude
1199 factors and reference times of the parameterized (undiscretized)
1200 sources match or are accounted for.
1201 '''
1203 if 'forces' not in kwargs:
1204 kwargs['forces'] = num.vstack([s.forces for s in sources])
1206 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1208 def moments(self):
1209 return num.sum(self.forces**2, axis=1)
1211 def centroid(self):
1212 from pyrocko.gf.seismosizer import SFSource
1213 time, lat, lon, north_shift, east_shift, depth = \
1214 self.centroid_position()
1216 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1217 return SFSource(
1218 time=time,
1219 lat=lat,
1220 lon=lon,
1221 north_shift=north_shift,
1222 east_shift=east_shift,
1223 depth=depth,
1224 fn=fn,
1225 fe=fe,
1226 fd=fd)
1229class DiscretizedMTSource(DiscretizedSource):
1230 m6s = Array.T(
1231 shape=(None, 6), dtype=float,
1232 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1234 provided_schemes = (
1235 'elastic8',
1236 'elastic10',
1237 'elastic18',
1238 )
1240 def get_source_terms(self, scheme):
1241 self.check_scheme(scheme)
1242 return self.m6s
1244 def make_weights(self, receiver, scheme):
1245 self.check_scheme(scheme)
1247 m6s = self.m6s
1248 n = m6s.shape[0]
1249 rep = num.repeat
1251 if scheme == 'elastic18':
1252 w_n = m6s.flatten()
1253 g_n = num.tile(num.arange(0, 6), n)
1254 w_e = m6s.flatten()
1255 g_e = num.tile(num.arange(6, 12), n)
1256 w_d = m6s.flatten()
1257 g_d = num.tile(num.arange(12, 18), n)
1259 else:
1260 azis, bazis = self.azibazis_to(receiver)
1262 sa = num.sin(azis*d2r)
1263 ca = num.cos(azis*d2r)
1264 s2a = num.sin(2.*azis*d2r)
1265 c2a = num.cos(2.*azis*d2r)
1266 sb = num.sin(bazis*d2r-num.pi)
1267 cb = num.cos(bazis*d2r-num.pi)
1269 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1270 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1271 f2 = m6s[:, 2]
1272 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1273 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1274 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1276 cat = num.concatenate
1278 if scheme == 'elastic8':
1279 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1280 g_n = rep((0, 1, 2, 3, 4), n)
1281 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1282 g_e = rep((0, 1, 2, 3, 4), n)
1283 w_d = cat((f0, f1, f2))
1284 g_d = rep((5, 6, 7), n)
1286 elif scheme == 'elastic10':
1287 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1288 g_n = rep((0, 1, 2, 8, 3, 4), n)
1289 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1290 g_e = rep((0, 1, 2, 8, 3, 4), n)
1291 w_d = cat((f0, f1, f2, f5))
1292 g_d = rep((5, 6, 7, 9), n)
1294 else:
1295 assert False
1297 return (
1298 ('displacement.n', w_n, g_n),
1299 ('displacement.e', w_e, g_e),
1300 ('displacement.d', w_d, g_d),
1301 )
1303 def split(self):
1304 from pyrocko.gf.seismosizer import MTSource
1305 sources = []
1306 for i in range(self.nelements):
1307 lat, lon, north_shift, east_shift = self.element_coords(i)
1308 sources.append(MTSource(
1309 time=float(self.times[i]),
1310 lat=lat,
1311 lon=lon,
1312 north_shift=north_shift,
1313 east_shift=east_shift,
1314 depth=float(self.depths[i]),
1315 m6=self.m6s[i]))
1317 return sources
1319 def moments(self):
1320 moments = num.array(
1321 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))
1322 for m6 in self.m6s])
1323 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)
1325 def get_moment_rate(self, deltat=None):
1326 moments = self.moments()
1327 times = self.times
1328 times -= times.min()
1330 t_max = times.max()
1331 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat
1332 mom_times[mom_times > t_max] = t_max
1334 # Right open histrogram bins
1335 mom, _ = num.histogram(
1336 times,
1337 bins=mom_times,
1338 weights=moments)
1340 deltat = num.diff(mom_times)
1341 mom_rate = mom / deltat
1343 return mom_rate, mom_times[1:]
1345 def centroid(self):
1346 from pyrocko.gf.seismosizer import MTSource
1347 time, lat, lon, north_shift, east_shift, depth = \
1348 self.centroid_position()
1350 return MTSource(
1351 time=time,
1352 lat=lat,
1353 lon=lon,
1354 north_shift=north_shift,
1355 east_shift=east_shift,
1356 depth=depth,
1357 m6=num.sum(self.m6s, axis=0))
1359 @classmethod
1360 def combine(cls, sources, **kwargs):
1361 '''
1362 Combine several discretized source models.
1364 Concatenenates all point sources in the given discretized ``sources``.
1365 Care must be taken when using this function that the external amplitude
1366 factors and reference times of the parameterized (undiscretized)
1367 sources match or are accounted for.
1368 '''
1370 if 'm6s' not in kwargs:
1371 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1373 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1376class DiscretizedPorePressureSource(DiscretizedSource):
1377 pp = Array.T(shape=(None,), dtype=float)
1379 provided_schemes = (
1380 'poroelastic10',
1381 )
1383 def get_source_terms(self, scheme):
1384 self.check_scheme(scheme)
1385 return self.pp[:, num.newaxis].copy()
1387 def make_weights(self, receiver, scheme):
1388 self.check_scheme(scheme)
1390 azis, bazis = self.azibazis_to(receiver)
1392 sb = num.sin(bazis*d2r-num.pi)
1393 cb = num.cos(bazis*d2r-num.pi)
1395 pp = self.pp
1396 n = bazis.size
1398 w_un = cb*pp
1399 g_un = filledi(1, n)
1400 w_ue = sb*pp
1401 g_ue = filledi(1, n)
1402 w_ud = pp
1403 g_ud = filledi(0, n)
1405 w_tn = cb*pp
1406 g_tn = filledi(6, n)
1407 w_te = sb*pp
1408 g_te = filledi(6, n)
1410 w_pp = pp
1411 g_pp = filledi(7, n)
1413 w_dvn = cb*pp
1414 g_dvn = filledi(9, n)
1415 w_dve = sb*pp
1416 g_dve = filledi(9, n)
1417 w_dvd = pp
1418 g_dvd = filledi(8, n)
1420 return (
1421 ('displacement.n', w_un, g_un),
1422 ('displacement.e', w_ue, g_ue),
1423 ('displacement.d', w_ud, g_ud),
1424 ('vertical_tilt.n', w_tn, g_tn),
1425 ('vertical_tilt.e', w_te, g_te),
1426 ('pore_pressure', w_pp, g_pp),
1427 ('darcy_velocity.n', w_dvn, g_dvn),
1428 ('darcy_velocity.e', w_dve, g_dve),
1429 ('darcy_velocity.d', w_dvd, g_dvd),
1430 )
1432 def moments(self):
1433 return self.pp
1435 def centroid(self):
1436 from pyrocko.gf.seismosizer import PorePressurePointSource
1437 time, lat, lon, north_shift, east_shift, depth = \
1438 self.centroid_position()
1440 return PorePressurePointSource(
1441 time=time,
1442 lat=lat,
1443 lon=lon,
1444 north_shift=north_shift,
1445 east_shift=east_shift,
1446 depth=depth,
1447 pp=float(num.sum(self.pp)))
1449 @classmethod
1450 def combine(cls, sources, **kwargs):
1451 '''
1452 Combine several discretized source models.
1454 Concatenenates all point sources in the given discretized ``sources``.
1455 Care must be taken when using this function that the external amplitude
1456 factors and reference times of the parameterized (undiscretized)
1457 sources match or are accounted for.
1458 '''
1460 if 'pp' not in kwargs:
1461 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1463 return super(DiscretizedPorePressureSource, cls).combine(sources,
1464 **kwargs)
1467class Region(Object):
1468 name = String.T(optional=True)
1471class RectangularRegion(Region):
1472 lat_min = Float.T()
1473 lat_max = Float.T()
1474 lon_min = Float.T()
1475 lon_max = Float.T()
1478class CircularRegion(Region):
1479 lat = Float.T()
1480 lon = Float.T()
1481 radius = Float.T()
1484class Config(Object):
1485 '''
1486 Green's function store meta information.
1488 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1489 configuration types are:
1491 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1492 spherical symmetry, 1D earth model, single receiver depth
1494 * Problem is invariant to horizontal translations and rotations around
1495 vertical axis.
1496 * All receivers must be at the same depth (e.g. at the surface)
1497 * High level index variables: ``(source_depth, receiver_distance,
1498 component)``
1500 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1501 spherical symmetry, 1D earth model, variable receiver depth
1503 * Symmetries like in Type A but has additional index for receiver depth
1504 * High level index variables: ``(source_depth, receiver_distance,
1505 receiver_depth, component)``
1507 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1508 constraints but fixed receiver positions
1510 * Cartesian source volume around a reference point
1511 * High level index variables: ``(ireceiver, source_depth,
1512 source_east_shift, source_north_shift, component)``
1513 '''
1515 id = StringID.T(
1516 help='Name of the store. May consist of upper and lower-case letters, '
1517 'digits, dots and underscores. The name must start with a '
1518 'letter.')
1520 derived_from_id = StringID.T(
1521 optional=True,
1522 help='Name of the original store, if this store has been derived from '
1523 'another one (e.g. extracted subset).')
1525 version = String.T(
1526 default='1.0',
1527 optional=True,
1528 help='User-defined version string. Use <major>.<minor> format.')
1530 modelling_code_id = StringID.T(
1531 optional=True,
1532 help='Identifier of the backend used to compute the store.')
1534 author = Unicode.T(
1535 optional=True,
1536 help='Comma-separated list of author names.')
1538 author_email = String.T(
1539 optional=True,
1540 help="Author's contact email address.")
1542 created_time = Timestamp.T(
1543 optional=True,
1544 help='Time of creation of the store.')
1546 regions = List.T(
1547 Region.T(),
1548 help='Geographical regions for which the store is representative.')
1550 scope_type = ScopeType.T(
1551 optional=True,
1552 help='Distance range scope of the store '
1553 '(%s).' % fmt_choices(ScopeType))
1555 waveform_type = WaveType.T(
1556 optional=True,
1557 help='Wave type stored (%s).' % fmt_choices(WaveType))
1559 nearfield_terms = NearfieldTermsType.T(
1560 optional=True,
1561 help='Information about the inclusion of near-field terms in the '
1562 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1564 description = String.T(
1565 optional=True,
1566 help='Free form textual description of the GF store.')
1568 references = List.T(
1569 Reference.T(),
1570 help='Reference list to cite the modelling code, earth model or '
1571 'related work.')
1573 earthmodel_1d = Earthmodel1D.T(
1574 optional=True,
1575 help='Layered earth model in ND (named discontinuity) format.')
1577 earthmodel_receiver_1d = Earthmodel1D.T(
1578 optional=True,
1579 help='Receiver-side layered earth model in ND format.')
1581 can_interpolate_source = Bool.T(
1582 optional=True,
1583 help='Hint to indicate if the spatial sampling of the store is dense '
1584 'enough for multi-linear interpolation at the source.')
1586 can_interpolate_receiver = Bool.T(
1587 optional=True,
1588 help='Hint to indicate if the spatial sampling of the store is dense '
1589 'enough for multi-linear interpolation at the receiver.')
1591 frequency_min = Float.T(
1592 optional=True,
1593 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1595 frequency_max = Float.T(
1596 optional=True,
1597 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1599 sample_rate = Float.T(
1600 optional=True,
1601 help='Sample rate of the GF store [Hz].')
1603 factor = Float.T(
1604 default=1.0,
1605 help='Gain value, factored out of the stored GF samples. '
1606 '(may not work properly, keep at 1.0).',
1607 optional=True)
1609 component_scheme = ComponentScheme.T(
1610 default='elastic10',
1611 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1613 stored_quantity = QuantityType.T(
1614 optional=True,
1615 help='Physical quantity of stored values (%s). If not given, a '
1616 'default is used based on the GF component scheme. The default '
1617 'for the ``"elastic*"`` family of component schemes is '
1618 '``"displacement"``.' % fmt_choices(QuantityType))
1620 tabulated_phases = List.T(
1621 TPDef.T(),
1622 help='Mapping of phase names to phase definitions, for which travel '
1623 'time tables are available in the GF store.')
1625 ncomponents = Int.T(
1626 optional=True,
1627 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1629 uuid = String.T(
1630 optional=True,
1631 help='Heuristic hash value which can be used to uniquely identify the '
1632 'GF store for practical purposes.')
1634 reference = String.T(
1635 optional=True,
1636 help='Store reference name composed of the store\'s :gattr:`id` and '
1637 'the first six letters of its :gattr:`uuid`.')
1639 def __init__(self, **kwargs):
1640 self._do_auto_updates = False
1641 Object.__init__(self, **kwargs)
1642 self._index_function = None
1643 self._indices_function = None
1644 self._vicinity_function = None
1645 self.validate(regularize=True, depth=1)
1646 self._do_auto_updates = True
1647 self.update()
1649 def check_ncomponents(self):
1650 ncomponents = component_scheme_to_description[
1651 self.component_scheme].ncomponents
1653 if self.ncomponents is None:
1654 self.ncomponents = ncomponents
1655 elif ncomponents != self.ncomponents:
1656 raise InvalidNComponents(
1657 'ncomponents=%i incompatible with component_scheme="%s"' % (
1658 self.ncomponents, self.component_scheme))
1660 def __setattr__(self, name, value):
1661 Object.__setattr__(self, name, value)
1662 try:
1663 self.T.get_property(name)
1664 if self._do_auto_updates:
1665 self.update()
1667 except ValueError:
1668 pass
1670 def update(self):
1671 self.check_ncomponents()
1672 self._update()
1673 self._make_index_functions()
1675 def irecord(self, *args):
1676 return self._index_function(*args)
1678 def irecords(self, *args):
1679 return self._indices_function(*args)
1681 def vicinity(self, *args):
1682 return self._vicinity_function(*args)
1684 def vicinities(self, *args):
1685 return self._vicinities_function(*args)
1687 def grid_interpolation_coefficients(self, *args):
1688 return self._grid_interpolation_coefficients(*args)
1690 def nodes(self, level=None, minlevel=None):
1691 return nodes(self.coords[minlevel:level])
1693 def iter_nodes(self, level=None, minlevel=None):
1694 return nditer_outer(self.coords[minlevel:level])
1696 def iter_extraction(self, gdef, level=None):
1697 i = 0
1698 arrs = []
1699 ntotal = 1
1700 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1701 if gdef and len(gdef) > i:
1702 sssn = gdef[i]
1703 else:
1704 sssn = (None,)*4
1706 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1707 ntotal *= len(arr)
1709 arrs.append(arr)
1710 i += 1
1712 arrs.append(self.coords[-1])
1713 return nditer_outer(arrs[:level])
1715 def make_sum_params(self, source, receiver, implementation='c',
1716 nthreads=0):
1717 assert implementation in ['c', 'python']
1719 out = []
1720 delays = source.times
1721 for comp, weights, icomponents in source.make_weights(
1722 receiver,
1723 self.component_scheme):
1725 weights *= self.factor
1727 args = self.make_indexing_args(source, receiver, icomponents)
1728 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1729 out.append((comp, args, delays_expanded, weights))
1731 return out
1733 def short_info(self):
1734 raise NotImplementedError('should be implemented in subclass')
1736 def get_shear_moduli(self, lat, lon, points,
1737 interpolation=None):
1738 '''
1739 Get shear moduli at given points from contained velocity model.
1741 :param lat: surface origin for coordinate system of ``points``
1742 :param points: NumPy array of shape ``(N, 3)``, where each row is
1743 a point ``(north, east, depth)``, relative to origin at
1744 ``(lat, lon)``
1745 :param interpolation: interpolation method. Choose from
1746 ``('nearest_neighbor', 'multilinear')``
1747 :returns: NumPy array of length N with extracted shear moduli at each
1748 point
1750 The default implementation retrieves and interpolates the shear moduli
1751 from the contained 1D velocity profile.
1752 '''
1753 return self.get_material_property(lat, lon, points,
1754 parameter='shear_moduli',
1755 interpolation=interpolation)
1757 def get_lambda_moduli(self, lat, lon, points,
1758 interpolation=None):
1759 '''
1760 Get lambda moduli at given points from contained velocity model.
1762 :param lat: surface origin for coordinate system of ``points``
1763 :param points: NumPy array of shape ``(N, 3)``, where each row is
1764 a point ``(north, east, depth)``, relative to origin at
1765 ``(lat, lon)``
1766 :param interpolation: interpolation method. Choose from
1767 ``('nearest_neighbor', 'multilinear')``
1768 :returns: NumPy array of length N with extracted shear moduli at each
1769 point
1771 The default implementation retrieves and interpolates the lambda moduli
1772 from the contained 1D velocity profile.
1773 '''
1774 return self.get_material_property(lat, lon, points,
1775 parameter='lambda_moduli',
1776 interpolation=interpolation)
1778 def get_bulk_moduli(self, lat, lon, points,
1779 interpolation=None):
1780 '''
1781 Get bulk moduli at given points from contained velocity model.
1783 :param lat: surface origin for coordinate system of ``points``
1784 :param points: NumPy array of shape ``(N, 3)``, where each row is
1785 a point ``(north, east, depth)``, relative to origin at
1786 ``(lat, lon)``
1787 :param interpolation: interpolation method. Choose from
1788 ``('nearest_neighbor', 'multilinear')``
1789 :returns: NumPy array of length N with extracted shear moduli at each
1790 point
1792 The default implementation retrieves and interpolates the lambda moduli
1793 from the contained 1D velocity profile.
1794 '''
1795 lambda_moduli = self.get_material_property(
1796 lat, lon, points, parameter='lambda_moduli',
1797 interpolation=interpolation)
1798 shear_moduli = self.get_material_property(
1799 lat, lon, points, parameter='shear_moduli',
1800 interpolation=interpolation)
1801 return lambda_moduli + (2 / 3) * shear_moduli
1803 def get_vs(self, lat, lon, points, interpolation=None):
1804 '''
1805 Get Vs at given points from contained velocity model.
1807 :param lat: surface origin for coordinate system of ``points``
1808 :param points: NumPy array of shape ``(N, 3)``, where each row is
1809 a point ``(north, east, depth)``, relative to origin at
1810 ``(lat, lon)``
1811 :param interpolation: interpolation method. Choose from
1812 ``('nearest_neighbor', 'multilinear')``
1813 :returns: NumPy array of length N with extracted shear moduli at each
1814 point
1816 The default implementation retrieves and interpolates Vs
1817 from the contained 1D velocity profile.
1818 '''
1819 return self.get_material_property(lat, lon, points,
1820 parameter='vs',
1821 interpolation=interpolation)
1823 def get_vp(self, lat, lon, points, interpolation=None):
1824 '''
1825 Get Vp at given points from contained velocity model.
1827 :param lat: surface origin for coordinate system of ``points``
1828 :param points: NumPy array of shape ``(N, 3)``, where each row is
1829 a point ``(north, east, depth)``, relative to origin at
1830 ``(lat, lon)``
1831 :param interpolation: interpolation method. Choose from
1832 ``('nearest_neighbor', 'multilinear')``
1833 :returns: NumPy array of length N with extracted shear moduli at each
1834 point
1836 The default implementation retrieves and interpolates Vp
1837 from the contained 1D velocity profile.
1838 '''
1839 return self.get_material_property(lat, lon, points,
1840 parameter='vp',
1841 interpolation=interpolation)
1843 def get_rho(self, lat, lon, points, interpolation=None):
1844 '''
1845 Get rho at given points from contained velocity model.
1847 :param lat: surface origin for coordinate system of ``points``
1848 :param points: NumPy array of shape ``(N, 3)``, where each row is
1849 a point ``(north, east, depth)``, relative to origin at
1850 ``(lat, lon)``
1851 :param interpolation: interpolation method. Choose from
1852 ``('nearest_neighbor', 'multilinear')``
1853 :returns: NumPy array of length N with extracted shear moduli at each
1854 point
1856 The default implementation retrieves and interpolates rho
1857 from the contained 1D velocity profile.
1858 '''
1859 return self.get_material_property(lat, lon, points,
1860 parameter='rho',
1861 interpolation=interpolation)
1863 def get_material_property(self, lat, lon, points, parameter='vs',
1864 interpolation=None):
1866 if interpolation is None:
1867 raise TypeError('Interpolation method not defined! available: '
1868 "multilinear", "nearest_neighbor")
1870 earthmod = self.earthmodel_1d
1871 store_depth_profile = self.get_source_depths()
1872 z_profile = earthmod.profile('z')
1874 if parameter == 'vs':
1875 vs_profile = earthmod.profile('vs')
1876 profile = num.interp(
1877 store_depth_profile, z_profile, vs_profile)
1879 elif parameter == 'vp':
1880 vp_profile = earthmod.profile('vp')
1881 profile = num.interp(
1882 store_depth_profile, z_profile, vp_profile)
1884 elif parameter == 'rho':
1885 rho_profile = earthmod.profile('rho')
1887 profile = num.interp(
1888 store_depth_profile, z_profile, rho_profile)
1890 elif parameter == 'shear_moduli':
1891 vs_profile = earthmod.profile('vs')
1892 rho_profile = earthmod.profile('rho')
1894 store_vs_profile = num.interp(
1895 store_depth_profile, z_profile, vs_profile)
1896 store_rho_profile = num.interp(
1897 store_depth_profile, z_profile, rho_profile)
1899 profile = num.power(store_vs_profile, 2) * store_rho_profile
1901 elif parameter == 'lambda_moduli':
1902 vs_profile = earthmod.profile('vs')
1903 vp_profile = earthmod.profile('vp')
1904 rho_profile = earthmod.profile('rho')
1906 store_vs_profile = num.interp(
1907 store_depth_profile, z_profile, vs_profile)
1908 store_vp_profile = num.interp(
1909 store_depth_profile, z_profile, vp_profile)
1910 store_rho_profile = num.interp(
1911 store_depth_profile, z_profile, rho_profile)
1913 profile = store_rho_profile * (
1914 num.power(store_vp_profile, 2) -
1915 num.power(store_vs_profile, 2) * 2)
1916 else:
1917 raise TypeError(
1918 'parameter %s not available' % parameter)
1920 if interpolation == 'multilinear':
1921 kind = 'linear'
1922 elif interpolation == 'nearest_neighbor':
1923 kind = 'nearest'
1924 else:
1925 raise TypeError(
1926 'Interpolation method %s not available' % interpolation)
1928 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1930 try:
1931 return interpolator(points[:, 2])
1932 except ValueError:
1933 raise OutOfBounds()
1935 def is_static(self):
1936 for code in ('psgrn_pscmp', 'poel'):
1937 if self.modelling_code_id.startswith(code):
1938 return True
1939 return False
1941 def is_dynamic(self):
1942 return not self.is_static()
1944 def get_source_depths(self):
1945 raise NotImplementedError('must be implemented in subclass')
1947 def get_tabulated_phase(self, phase_id):
1948 '''
1949 Get tabulated phase definition.
1950 '''
1952 for pdef in self.tabulated_phases:
1953 if pdef.id == phase_id:
1954 return pdef
1956 raise StoreError('No such phase: %s' % phase_id)
1958 def fix_ttt_holes(self, sptree, mode):
1959 raise StoreError(
1960 'Cannot fix travel time table holes in GF stores of type %s.'
1961 % self.short_type)
1964class ConfigTypeA(Config):
1965 '''
1966 Cylindrical symmetry, 1D earth model, single receiver depth
1968 * Problem is invariant to horizontal translations and rotations around
1969 vertical axis.
1971 * All receivers must be at the same depth (e.g. at the surface)
1972 High level index variables: ``(source_depth, distance,
1973 component)``
1975 * The ``distance`` is the surface distance between source and receiver
1976 points.
1977 '''
1979 receiver_depth = Float.T(
1980 default=0.0,
1981 help='Fixed receiver depth [m].')
1983 source_depth_min = Float.T(
1984 help='Minimum source depth [m].')
1986 source_depth_max = Float.T(
1987 help='Maximum source depth [m].')
1989 source_depth_delta = Float.T(
1990 help='Grid spacing of source depths [m]')
1992 distance_min = Float.T(
1993 help='Minimum source-receiver surface distance [m].')
1995 distance_max = Float.T(
1996 help='Maximum source-receiver surface distance [m].')
1998 distance_delta = Float.T(
1999 help='Grid spacing of source-receiver surface distance [m].')
2001 short_type = 'A'
2003 provided_schemes = [
2004 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2006 def get_surface_distance(self, args):
2007 return args[1]
2009 def get_distance(self, args):
2010 return math.sqrt(args[0]**2 + args[1]**2)
2012 def get_source_depth(self, args):
2013 return args[0]
2015 def get_source_depths(self):
2016 return self.coords[0]
2018 def get_receiver_depth(self, args):
2019 return self.receiver_depth
2021 def _update(self):
2022 self.mins = num.array(
2023 [self.source_depth_min, self.distance_min], dtype=float)
2024 self.maxs = num.array(
2025 [self.source_depth_max, self.distance_max], dtype=float)
2026 self.deltas = num.array(
2027 [self.source_depth_delta, self.distance_delta],
2028 dtype=float)
2029 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2030 vicinity_eps).astype(int) + 1
2031 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2032 self.deltat = 1.0/self.sample_rate
2033 self.nrecords = num.product(self.ns) * self.ncomponents
2034 self.coords = tuple(num.linspace(mi, ma, n) for
2035 (mi, ma, n) in
2036 zip(self.mins, self.effective_maxs, self.ns)) + \
2037 (num.arange(self.ncomponents),)
2039 self.nsource_depths, self.ndistances = self.ns
2041 def _make_index_functions(self):
2043 amin, bmin = self.mins
2044 da, db = self.deltas
2045 na, nb = self.ns
2047 ng = self.ncomponents
2049 def index_function(a, b, ig):
2050 ia = int(round((a - amin) / da))
2051 ib = int(round((b - bmin) / db))
2052 try:
2053 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2054 except ValueError:
2055 raise OutOfBounds()
2057 def indices_function(a, b, ig):
2058 ia = num.round((a - amin) / da).astype(int)
2059 ib = num.round((b - bmin) / db).astype(int)
2060 try:
2061 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2062 except ValueError:
2063 for ia_, ib_, ig_ in zip(ia, ib, ig):
2064 try:
2065 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2066 except ValueError:
2067 raise OutOfBounds()
2069 def grid_interpolation_coefficients(a, b):
2070 ias = indi12((a - amin) / da, na)
2071 ibs = indi12((b - bmin) / db, nb)
2072 return ias, ibs
2074 def vicinity_function(a, b, ig):
2075 ias, ibs = grid_interpolation_coefficients(a, b)
2077 if not (0 <= ig < ng):
2078 raise OutOfBounds()
2080 indis = []
2081 weights = []
2082 for ia, va in ias:
2083 iia = ia*nb*ng
2084 for ib, vb in ibs:
2085 indis.append(iia + ib*ng + ig)
2086 weights.append(va*vb)
2088 return num.array(indis), num.array(weights)
2090 def vicinities_function(a, b, ig):
2092 xa = (a - amin) / da
2093 xb = (b - bmin) / db
2095 xa_fl = num.floor(xa)
2096 xa_ce = num.ceil(xa)
2097 xb_fl = num.floor(xb)
2098 xb_ce = num.ceil(xb)
2099 va_fl = 1.0 - (xa - xa_fl)
2100 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2101 vb_fl = 1.0 - (xb - xb_fl)
2102 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2104 ia_fl = xa_fl.astype(int)
2105 ia_ce = xa_ce.astype(int)
2106 ib_fl = xb_fl.astype(int)
2107 ib_ce = xb_ce.astype(int)
2109 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2110 raise OutOfBounds()
2112 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2113 raise OutOfBounds()
2115 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2116 raise OutOfBounds()
2118 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2119 raise OutOfBounds()
2121 irecords = num.empty(a.size*4, dtype=int)
2122 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2123 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2124 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2125 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2127 weights = num.empty(a.size*4, dtype=float)
2128 weights[0::4] = va_fl * vb_fl
2129 weights[1::4] = va_ce * vb_fl
2130 weights[2::4] = va_fl * vb_ce
2131 weights[3::4] = va_ce * vb_ce
2133 return irecords, weights
2135 self._index_function = index_function
2136 self._indices_function = indices_function
2137 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2138 self._vicinity_function = vicinity_function
2139 self._vicinities_function = vicinities_function
2141 def make_indexing_args(self, source, receiver, icomponents):
2142 nc = icomponents.size
2143 dists = source.distances_to(receiver)
2144 n = dists.size
2145 return (num.tile(source.depths, nc//n),
2146 num.tile(dists, nc//n),
2147 icomponents)
2149 def make_indexing_args1(self, source, receiver):
2150 return (source.depth, source.distance_to(receiver))
2152 @property
2153 def short_extent(self):
2154 return '%g:%g:%g x %g:%g:%g' % (
2155 self.source_depth_min/km,
2156 self.source_depth_max/km,
2157 self.source_depth_delta/km,
2158 self.distance_min/km,
2159 self.distance_max/km,
2160 self.distance_delta/km)
2162 def fix_ttt_holes(self, sptree, mode):
2163 from pyrocko import eikonal_ext, spit
2165 nodes = self.nodes(level=-1)
2167 delta = self.deltas[-1]
2168 assert num.all(delta == self.deltas)
2170 nsources, ndistances = self.ns
2172 points = num.zeros((nodes.shape[0], 3))
2173 points[:, 0] = nodes[:, 1]
2174 points[:, 2] = nodes[:, 0]
2176 speeds = self.get_material_property(
2177 0., 0., points,
2178 parameter='vp' if mode == cake.P else 'vs',
2179 interpolation='multilinear')
2181 speeds = speeds.reshape((nsources, ndistances))
2183 times = sptree.interpolate_many(nodes)
2185 times[num.isnan(times)] = -1.
2186 times = times.reshape(speeds.shape)
2188 try:
2189 eikonal_ext.eikonal_solver_fmm_cartesian(
2190 speeds, times, delta)
2191 except eikonal_ext.EikonalExtError as e:
2192 if str(e).endswith('please check results'):
2193 logger.debug(
2194 'Got a warning from eikonal solver '
2195 '- may be ok...')
2196 else:
2197 raise
2199 def func(x):
2200 ibs, ics = \
2201 self.grid_interpolation_coefficients(*x)
2203 t = 0
2204 for ib, vb in ibs:
2205 for ic, vc in ics:
2206 t += times[ib, ic] * vb * vc
2208 return t
2210 return spit.SPTree(
2211 f=func,
2212 ftol=sptree.ftol,
2213 xbounds=sptree.xbounds,
2214 xtols=sptree.xtols)
2217class ConfigTypeB(Config):
2218 '''
2219 Cylindrical symmetry, 1D earth model, variable receiver depth
2221 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2222 receiver depth
2224 * High level index variables: ``(receiver_depth, source_depth,
2225 receiver_distance, component)``
2226 '''
2228 receiver_depth_min = Float.T(
2229 help='Minimum receiver depth [m].')
2231 receiver_depth_max = Float.T(
2232 help='Maximum receiver depth [m].')
2234 receiver_depth_delta = Float.T(
2235 help='Grid spacing of receiver depths [m]')
2237 source_depth_min = Float.T(
2238 help='Minimum source depth [m].')
2240 source_depth_max = Float.T(
2241 help='Maximum source depth [m].')
2243 source_depth_delta = Float.T(
2244 help='Grid spacing of source depths [m]')
2246 distance_min = Float.T(
2247 help='Minimum source-receiver surface distance [m].')
2249 distance_max = Float.T(
2250 help='Maximum source-receiver surface distance [m].')
2252 distance_delta = Float.T(
2253 help='Grid spacing of source-receiver surface distances [m].')
2255 short_type = 'B'
2257 provided_schemes = [
2258 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2260 def get_distance(self, args):
2261 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2263 def get_surface_distance(self, args):
2264 return args[2]
2266 def get_source_depth(self, args):
2267 return args[1]
2269 def get_receiver_depth(self, args):
2270 return args[0]
2272 def get_source_depths(self):
2273 return self.coords[1]
2275 def _update(self):
2276 self.mins = num.array([
2277 self.receiver_depth_min,
2278 self.source_depth_min,
2279 self.distance_min],
2280 dtype=float)
2282 self.maxs = num.array([
2283 self.receiver_depth_max,
2284 self.source_depth_max,
2285 self.distance_max],
2286 dtype=float)
2288 self.deltas = num.array([
2289 self.receiver_depth_delta,
2290 self.source_depth_delta,
2291 self.distance_delta],
2292 dtype=float)
2294 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2295 vicinity_eps).astype(int) + 1
2296 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2297 self.deltat = 1.0/self.sample_rate
2298 self.nrecords = num.product(self.ns) * self.ncomponents
2299 self.coords = tuple(num.linspace(mi, ma, n) for
2300 (mi, ma, n) in
2301 zip(self.mins, self.effective_maxs, self.ns)) + \
2302 (num.arange(self.ncomponents),)
2303 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2305 def _make_index_functions(self):
2307 amin, bmin, cmin = self.mins
2308 da, db, dc = self.deltas
2309 na, nb, nc = self.ns
2310 ng = self.ncomponents
2312 def index_function(a, b, c, ig):
2313 ia = int(round((a - amin) / da))
2314 ib = int(round((b - bmin) / db))
2315 ic = int(round((c - cmin) / dc))
2316 try:
2317 return num.ravel_multi_index((ia, ib, ic, ig),
2318 (na, nb, nc, ng))
2319 except ValueError:
2320 raise OutOfBounds()
2322 def indices_function(a, b, c, ig):
2323 ia = num.round((a - amin) / da).astype(int)
2324 ib = num.round((b - bmin) / db).astype(int)
2325 ic = num.round((c - cmin) / dc).astype(int)
2326 try:
2327 return num.ravel_multi_index((ia, ib, ic, ig),
2328 (na, nb, nc, ng))
2329 except ValueError:
2330 raise OutOfBounds()
2332 def grid_interpolation_coefficients(a, b, c):
2333 ias = indi12((a - amin) / da, na)
2334 ibs = indi12((b - bmin) / db, nb)
2335 ics = indi12((c - cmin) / dc, nc)
2336 return ias, ibs, ics
2338 def vicinity_function(a, b, c, ig):
2339 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2341 if not (0 <= ig < ng):
2342 raise OutOfBounds()
2344 indis = []
2345 weights = []
2346 for ia, va in ias:
2347 iia = ia*nb*nc*ng
2348 for ib, vb in ibs:
2349 iib = ib*nc*ng
2350 for ic, vc in ics:
2351 indis.append(iia + iib + ic*ng + ig)
2352 weights.append(va*vb*vc)
2354 return num.array(indis), num.array(weights)
2356 def vicinities_function(a, b, c, ig):
2358 xa = (a - amin) / da
2359 xb = (b - bmin) / db
2360 xc = (c - cmin) / dc
2362 xa_fl = num.floor(xa)
2363 xa_ce = num.ceil(xa)
2364 xb_fl = num.floor(xb)
2365 xb_ce = num.ceil(xb)
2366 xc_fl = num.floor(xc)
2367 xc_ce = num.ceil(xc)
2368 va_fl = 1.0 - (xa - xa_fl)
2369 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2370 vb_fl = 1.0 - (xb - xb_fl)
2371 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2372 vc_fl = 1.0 - (xc - xc_fl)
2373 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2375 ia_fl = xa_fl.astype(int)
2376 ia_ce = xa_ce.astype(int)
2377 ib_fl = xb_fl.astype(int)
2378 ib_ce = xb_ce.astype(int)
2379 ic_fl = xc_fl.astype(int)
2380 ic_ce = xc_ce.astype(int)
2382 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2383 raise OutOfBounds()
2385 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2386 raise OutOfBounds()
2388 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2389 raise OutOfBounds()
2391 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2392 raise OutOfBounds()
2394 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2395 raise OutOfBounds()
2397 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2398 raise OutOfBounds()
2400 irecords = num.empty(a.size*8, dtype=int)
2401 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2402 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2403 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2404 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2405 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2406 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2407 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2408 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2410 weights = num.empty(a.size*8, dtype=float)
2411 weights[0::8] = va_fl * vb_fl * vc_fl
2412 weights[1::8] = va_ce * vb_fl * vc_fl
2413 weights[2::8] = va_fl * vb_ce * vc_fl
2414 weights[3::8] = va_ce * vb_ce * vc_fl
2415 weights[4::8] = va_fl * vb_fl * vc_ce
2416 weights[5::8] = va_ce * vb_fl * vc_ce
2417 weights[6::8] = va_fl * vb_ce * vc_ce
2418 weights[7::8] = va_ce * vb_ce * vc_ce
2420 return irecords, weights
2422 self._index_function = index_function
2423 self._indices_function = indices_function
2424 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2425 self._vicinity_function = vicinity_function
2426 self._vicinities_function = vicinities_function
2428 def make_indexing_args(self, source, receiver, icomponents):
2429 nc = icomponents.size
2430 dists = source.distances_to(receiver)
2431 n = dists.size
2432 receiver_depths = num.empty(nc)
2433 receiver_depths.fill(receiver.depth)
2434 return (receiver_depths,
2435 num.tile(source.depths, nc//n),
2436 num.tile(dists, nc//n),
2437 icomponents)
2439 def make_indexing_args1(self, source, receiver):
2440 return (receiver.depth,
2441 source.depth,
2442 source.distance_to(receiver))
2444 @property
2445 def short_extent(self):
2446 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2447 self.receiver_depth_min/km,
2448 self.receiver_depth_max/km,
2449 self.receiver_depth_delta/km,
2450 self.source_depth_min/km,
2451 self.source_depth_max/km,
2452 self.source_depth_delta/km,
2453 self.distance_min/km,
2454 self.distance_max/km,
2455 self.distance_delta/km)
2457 def fix_ttt_holes(self, sptree, mode):
2458 from pyrocko import eikonal_ext, spit
2460 nodes_sr = self.nodes(minlevel=1, level=-1)
2462 delta = self.deltas[-1]
2463 assert num.all(delta == self.deltas[1:])
2465 nreceivers, nsources, ndistances = self.ns
2467 points = num.zeros((nodes_sr.shape[0], 3))
2468 points[:, 0] = nodes_sr[:, 1]
2469 points[:, 2] = nodes_sr[:, 0]
2471 speeds = self.get_material_property(
2472 0., 0., points,
2473 parameter='vp' if mode == cake.P else 'vs',
2474 interpolation='multilinear')
2476 speeds = speeds.reshape((nsources, ndistances))
2478 receiver_times = []
2479 for ireceiver in range(nreceivers):
2480 nodes = num.hstack([
2481 num_full(
2482 (nodes_sr.shape[0], 1),
2483 self.coords[0][ireceiver]),
2484 nodes_sr])
2486 times = sptree.interpolate_many(nodes)
2488 times[num.isnan(times)] = -1.
2490 times = times.reshape(speeds.shape)
2492 try:
2493 eikonal_ext.eikonal_solver_fmm_cartesian(
2494 speeds, times, delta)
2495 except eikonal_ext.EikonalExtError as e:
2496 if str(e).endswith('please check results'):
2497 logger.debug(
2498 'Got a warning from eikonal solver '
2499 '- may be ok...')
2500 else:
2501 raise
2503 receiver_times.append(times)
2505 def func(x):
2506 ias, ibs, ics = \
2507 self.grid_interpolation_coefficients(*x)
2509 t = 0
2510 for ia, va in ias:
2511 times = receiver_times[ia]
2512 for ib, vb in ibs:
2513 for ic, vc in ics:
2514 t += times[ib, ic] * va * vb * vc
2516 return t
2518 return spit.SPTree(
2519 f=func,
2520 ftol=sptree.ftol,
2521 xbounds=sptree.xbounds,
2522 xtols=sptree.xtols)
2525class ConfigTypeC(Config):
2526 '''
2527 No symmetrical constraints, one fixed receiver position.
2529 * Cartesian 3D source volume around a reference point
2531 * High level index variables: ``(source_depth,
2532 source_east_shift, source_north_shift, component)``
2533 '''
2535 receiver = Receiver.T(
2536 help='Receiver location')
2538 source_origin = Location.T(
2539 help='Origin of the source volume grid.')
2541 source_depth_min = Float.T(
2542 help='Minimum source depth [m].')
2544 source_depth_max = Float.T(
2545 help='Maximum source depth [m].')
2547 source_depth_delta = Float.T(
2548 help='Source depth grid spacing [m].')
2550 source_east_shift_min = Float.T(
2551 help='Minimum easting of source grid [m].')
2553 source_east_shift_max = Float.T(
2554 help='Maximum easting of source grid [m].')
2556 source_east_shift_delta = Float.T(
2557 help='Source volume grid spacing in east direction [m].')
2559 source_north_shift_min = Float.T(
2560 help='Minimum northing of source grid [m].')
2562 source_north_shift_max = Float.T(
2563 help='Maximum northing of source grid [m].')
2565 source_north_shift_delta = Float.T(
2566 help='Source volume grid spacing in north direction [m].')
2568 short_type = 'C'
2570 provided_schemes = ['elastic18']
2572 def get_surface_distance(self, args):
2573 _, source_east_shift, source_north_shift, _ = args
2574 sorig = self.source_origin
2575 sloc = Location(
2576 lat=sorig.lat,
2577 lon=sorig.lon,
2578 north_shift=sorig.north_shift + source_north_shift,
2579 east_shift=sorig.east_shift + source_east_shift)
2581 return self.receiver.distance_to(sloc)
2583 def get_distance(self, args):
2584 sdepth, source_east_shift, source_north_shift, _ = args
2585 sorig = self.source_origin
2586 sloc = Location(
2587 lat=sorig.lat,
2588 lon=sorig.lon,
2589 north_shift=sorig.north_shift + source_north_shift,
2590 east_shift=sorig.east_shift + source_east_shift)
2592 return self.receiver.distance_3d_to(sloc)
2594 def get_source_depth(self, args):
2595 return args[0]
2597 def get_receiver_depth(self, args):
2598 return self.receiver.depth
2600 def get_source_depths(self):
2601 return self.coords[0]
2603 def _update(self):
2604 self.mins = num.array([
2605 self.source_depth_min,
2606 self.source_east_shift_min,
2607 self.source_north_shift_min],
2608 dtype=float)
2610 self.maxs = num.array([
2611 self.source_depth_max,
2612 self.source_east_shift_max,
2613 self.source_north_shift_max],
2614 dtype=float)
2616 self.deltas = num.array([
2617 self.source_depth_delta,
2618 self.source_east_shift_delta,
2619 self.source_north_shift_delta],
2620 dtype=float)
2622 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2623 vicinity_eps).astype(int) + 1
2624 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2625 self.deltat = 1.0/self.sample_rate
2626 self.nrecords = num.product(self.ns) * self.ncomponents
2628 self.coords = tuple(
2629 num.linspace(mi, ma, n) for (mi, ma, n) in
2630 zip(self.mins, self.effective_maxs, self.ns)) + \
2631 (num.arange(self.ncomponents),)
2633 def _make_index_functions(self):
2635 amin, bmin, cmin = self.mins
2636 da, db, dc = self.deltas
2637 na, nb, nc = self.ns
2638 ng = self.ncomponents
2640 def index_function(a, b, c, ig):
2641 ia = int(round((a - amin) / da))
2642 ib = int(round((b - bmin) / db))
2643 ic = int(round((c - cmin) / dc))
2644 try:
2645 return num.ravel_multi_index((ia, ib, ic, ig),
2646 (na, nb, nc, ng))
2647 except ValueError:
2648 raise OutOfBounds(values=(a, b, c, ig))
2650 def indices_function(a, b, c, ig):
2651 ia = num.round((a - amin) / da).astype(int)
2652 ib = num.round((b - bmin) / db).astype(int)
2653 ic = num.round((c - cmin) / dc).astype(int)
2655 try:
2656 return num.ravel_multi_index((ia, ib, ic, ig),
2657 (na, nb, nc, ng))
2658 except ValueError:
2659 raise OutOfBounds()
2661 def vicinity_function(a, b, c, ig):
2662 ias = indi12((a - amin) / da, na)
2663 ibs = indi12((b - bmin) / db, nb)
2664 ics = indi12((c - cmin) / dc, nc)
2666 if not (0 <= ig < ng):
2667 raise OutOfBounds()
2669 indis = []
2670 weights = []
2671 for ia, va in ias:
2672 iia = ia*nb*nc*ng
2673 for ib, vb in ibs:
2674 iib = ib*nc*ng
2675 for ic, vc in ics:
2676 indis.append(iia + iib + ic*ng + ig)
2677 weights.append(va*vb*vc)
2679 return num.array(indis), num.array(weights)
2681 def vicinities_function(a, b, c, ig):
2683 xa = (a-amin) / da
2684 xb = (b-bmin) / db
2685 xc = (c-cmin) / dc
2687 xa_fl = num.floor(xa)
2688 xa_ce = num.ceil(xa)
2689 xb_fl = num.floor(xb)
2690 xb_ce = num.ceil(xb)
2691 xc_fl = num.floor(xc)
2692 xc_ce = num.ceil(xc)
2693 va_fl = 1.0 - (xa - xa_fl)
2694 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2695 vb_fl = 1.0 - (xb - xb_fl)
2696 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2697 vc_fl = 1.0 - (xc - xc_fl)
2698 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2700 ia_fl = xa_fl.astype(int)
2701 ia_ce = xa_ce.astype(int)
2702 ib_fl = xb_fl.astype(int)
2703 ib_ce = xb_ce.astype(int)
2704 ic_fl = xc_fl.astype(int)
2705 ic_ce = xc_ce.astype(int)
2707 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2708 raise OutOfBounds()
2710 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2711 raise OutOfBounds()
2713 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2714 raise OutOfBounds()
2716 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2717 raise OutOfBounds()
2719 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2720 raise OutOfBounds()
2722 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2723 raise OutOfBounds()
2725 irecords = num.empty(a.size*8, dtype=int)
2726 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2727 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2728 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2729 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2730 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2731 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2732 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2733 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2735 weights = num.empty(a.size*8, dtype=float)
2736 weights[0::8] = va_fl * vb_fl * vc_fl
2737 weights[1::8] = va_ce * vb_fl * vc_fl
2738 weights[2::8] = va_fl * vb_ce * vc_fl
2739 weights[3::8] = va_ce * vb_ce * vc_fl
2740 weights[4::8] = va_fl * vb_fl * vc_ce
2741 weights[5::8] = va_ce * vb_fl * vc_ce
2742 weights[6::8] = va_fl * vb_ce * vc_ce
2743 weights[7::8] = va_ce * vb_ce * vc_ce
2745 return irecords, weights
2747 self._index_function = index_function
2748 self._indices_function = indices_function
2749 self._vicinity_function = vicinity_function
2750 self._vicinities_function = vicinities_function
2752 def make_indexing_args(self, source, receiver, icomponents):
2753 nc = icomponents.size
2755 dists = source.distances_to(self.source_origin)
2756 azis, _ = source.azibazis_to(self.source_origin)
2758 source_north_shifts = - num.cos(d2r*azis) * dists
2759 source_east_shifts = - num.sin(d2r*azis) * dists
2760 source_depths = source.depths - self.source_origin.depth
2762 n = dists.size
2764 return (num.tile(source_depths, nc//n),
2765 num.tile(source_east_shifts, nc//n),
2766 num.tile(source_north_shifts, nc//n),
2767 icomponents)
2769 def make_indexing_args1(self, source, receiver):
2770 dist = source.distance_to(self.source_origin)
2771 azi, _ = source.azibazi_to(self.source_origin)
2773 source_north_shift = - num.cos(d2r*azi) * dist
2774 source_east_shift = - num.sin(d2r*azi) * dist
2775 source_depth = source.depth - self.source_origin.depth
2777 return (source_depth,
2778 source_east_shift,
2779 source_north_shift)
2781 @property
2782 def short_extent(self):
2783 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2784 self.source_depth_min/km,
2785 self.source_depth_max/km,
2786 self.source_depth_delta/km,
2787 self.source_east_shift_min/km,
2788 self.source_east_shift_max/km,
2789 self.source_east_shift_delta/km,
2790 self.source_north_shift_min/km,
2791 self.source_north_shift_max/km,
2792 self.source_north_shift_delta/km)
2795class Weighting(Object):
2796 factor = Float.T(default=1.0)
2799class Taper(Object):
2800 tmin = Timing.T()
2801 tmax = Timing.T()
2802 tfade = Float.T(default=0.0)
2803 shape = StringChoice.T(
2804 choices=['cos', 'linear'],
2805 default='cos',
2806 optional=True)
2809class SimplePattern(SObject):
2811 _pool = {}
2813 def __init__(self, pattern):
2814 self._pattern = pattern
2815 SObject.__init__(self)
2817 def __str__(self):
2818 return self._pattern
2820 @property
2821 def regex(self):
2822 pool = SimplePattern._pool
2823 if self.pattern not in pool:
2824 rpat = '|'.join(fnmatch.translate(x) for
2825 x in self.pattern.split('|'))
2826 pool[self.pattern] = re.compile(rpat, re.I)
2828 return pool[self.pattern]
2830 def match(self, s):
2831 return self.regex.match(s)
2834class WaveformType(StringChoice):
2835 choices = ['dis', 'vel', 'acc',
2836 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2837 'envelope_dis', 'envelope_vel', 'envelope_acc']
2840class ChannelSelection(Object):
2841 pattern = SimplePattern.T(optional=True)
2842 min_sample_rate = Float.T(optional=True)
2843 max_sample_rate = Float.T(optional=True)
2846class StationSelection(Object):
2847 includes = SimplePattern.T()
2848 excludes = SimplePattern.T()
2849 distance_min = Float.T(optional=True)
2850 distance_max = Float.T(optional=True)
2851 azimuth_min = Float.T(optional=True)
2852 azimuth_max = Float.T(optional=True)
2855class WaveformSelection(Object):
2856 channel_selection = ChannelSelection.T(optional=True)
2857 station_selection = StationSelection.T(optional=True)
2858 taper = Taper.T()
2859 # filter = FrequencyResponse.T()
2860 waveform_type = WaveformType.T(default='dis')
2861 weighting = Weighting.T(optional=True)
2862 sample_rate = Float.T(optional=True)
2863 gf_store_id = StringID.T(optional=True)
2866def indi12(x, n):
2867 '''
2868 Get linear interpolation index and weight.
2869 '''
2871 r = round(x)
2872 if abs(r - x) < vicinity_eps:
2873 i = int(r)
2874 if not (0 <= i < n):
2875 raise OutOfBounds()
2877 return ((int(r), 1.),)
2878 else:
2879 f = math.floor(x)
2880 i = int(f)
2881 if not (0 <= i < n-1):
2882 raise OutOfBounds()
2884 v = x-f
2885 return ((i, 1.-v), (i + 1, v))
2888def float_or_none(s):
2889 units = {
2890 'k': 1e3,
2891 'M': 1e6,
2892 }
2894 factor = 1.0
2895 if s and s[-1] in units:
2896 factor = units[s[-1]]
2897 s = s[:-1]
2898 if not s:
2899 raise ValueError('unit without a number: \'%s\'' % s)
2901 if s:
2902 return float(s) * factor
2903 else:
2904 return None
2907class GridSpecError(Exception):
2908 def __init__(self, s):
2909 Exception.__init__(self, 'invalid grid specification: %s' % s)
2912def parse_grid_spec(spec):
2913 try:
2914 result = []
2915 for dspec in spec.split(','):
2916 t = dspec.split('@')
2917 num = start = stop = step = None
2918 if len(t) == 2:
2919 num = int(t[1])
2920 if num <= 0:
2921 raise GridSpecError(spec)
2923 elif len(t) > 2:
2924 raise GridSpecError(spec)
2926 s = t[0]
2927 v = [float_or_none(x) for x in s.split(':')]
2928 if len(v) == 1:
2929 start = stop = v[0]
2930 if len(v) >= 2:
2931 start, stop = v[0:2]
2932 if len(v) == 3:
2933 step = v[2]
2935 if len(v) > 3 or (len(v) > 2 and num is not None):
2936 raise GridSpecError(spec)
2938 if step == 0.0:
2939 raise GridSpecError(spec)
2941 result.append((start, stop, step, num))
2943 except ValueError:
2944 raise GridSpecError(spec)
2946 return result
2949def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2950 swap = step is not None and step < 0.
2951 if start is None:
2952 start = ma if swap else mi
2953 if stop is None:
2954 stop = mi if swap else ma
2955 if step is None:
2956 step = -inc if ma < mi else inc
2957 if num is None:
2958 if (step < 0) != (stop-start < 0):
2959 raise GridSpecError()
2961 num = int(round((stop-start)/step))+1
2962 stop2 = start + (num-1)*step
2963 if abs(stop-stop2) > eps:
2964 num = int(math.floor((stop-start)/step))+1
2965 stop = start + (num-1)*step
2966 else:
2967 stop = stop2
2969 if start == stop:
2970 num = 1
2972 return start, stop, num
2975def nditer_outer(x):
2976 return num.nditer(
2977 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
2980def nodes(xs):
2981 ns = [x.size for x in xs]
2982 nnodes = num.prod(ns)
2983 ndim = len(xs)
2984 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
2985 for idim in range(ndim-1, -1, -1):
2986 x = xs[idim]
2987 nrepeat = num.prod(ns[idim+1:], dtype=int)
2988 ntile = num.prod(ns[:idim], dtype=int)
2989 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
2991 return nodes
2994def filledi(x, n):
2995 a = num.empty(n, dtype=int)
2996 a.fill(x)
2997 return a
3000config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3002discretized_source_classes = [
3003 DiscretizedExplosionSource,
3004 DiscretizedSFSource,
3005 DiscretizedMTSource,
3006 DiscretizedPorePressureSource]
3009__all__ = '''
3010Earthmodel1D
3011StringID
3012ScopeType
3013WaveformType
3014QuantityType
3015NearfieldTermsType
3016Reference
3017Region
3018CircularRegion
3019RectangularRegion
3020PhaseSelect
3021InvalidTimingSpecification
3022Timing
3023TPDef
3024OutOfBounds
3025Location
3026Receiver
3027'''.split() + [
3028 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3029ComponentScheme
3030component_scheme_to_description
3031component_schemes
3032Config
3033GridSpecError
3034Weighting
3035Taper
3036SimplePattern
3037WaveformType
3038ChannelSelection
3039StationSelection
3040WaveformSelection
3041nditer_outer
3042dump
3043load
3044discretized_source_classes
3045config_type_classes
3046UnavailableScheme
3047InterpolationMethod
3048SeismosizerTrace
3049SeismosizerResult
3050Result
3051StaticResult
3052'''.split()