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 print(times)
964 depths = cat([s.depths for s in sources])
966 if same_ref:
967 lat = first.lat
968 lon = first.lon
969 north_shifts = cat([s.effective_north_shifts for s in sources])
970 east_shifts = cat([s.effective_east_shifts for s in sources])
971 lats = None
972 lons = None
973 else:
974 lat = None
975 lon = None
976 north_shifts = None
977 east_shifts = None
979 return cls(
980 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
981 north_shifts=north_shifts, east_shifts=east_shifts,
982 depths=depths, **kwargs)
984 def centroid_position(self):
985 moments = self.moments()
986 norm = num.sum(moments)
987 if norm != 0.0:
988 w = moments / num.sum(moments)
989 else:
990 w = num.ones(moments.size)
992 if self.lats is not None and self.lons is not None:
993 lats, lons = self.effective_latlons
994 rlat, rlon = num.mean(lats), num.mean(lons)
995 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
996 else:
997 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
998 n, e = self.north_shifts, self.east_shifts
1000 cn = num.sum(n*w)
1001 ce = num.sum(e*w)
1002 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1004 if self.lats is not None and self.lons is not None:
1005 lat = clat
1006 lon = clon
1007 north_shift = 0.
1008 east_shift = 0.
1009 else:
1010 lat = g(self.lat, 0.0)
1011 lon = g(self.lon, 0.0)
1012 north_shift = cn
1013 east_shift = ce
1015 depth = num.sum(self.depths*w)
1016 time = num.sum(self.times*w)
1017 return tuple(float(x) for x in
1018 (time, lat, lon, north_shift, east_shift, depth))
1021class DiscretizedExplosionSource(DiscretizedSource):
1022 m0s = Array.T(shape=(None,), dtype=float)
1024 provided_schemes = (
1025 'elastic2',
1026 'elastic8',
1027 'elastic10',
1028 )
1030 def get_source_terms(self, scheme):
1031 self.check_scheme(scheme)
1033 if scheme == 'elastic2':
1034 return self.m0s[:, num.newaxis].copy()
1036 elif scheme in ('elastic8', 'elastic10'):
1037 m6s = num.zeros((self.m0s.size, 6))
1038 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1039 return m6s
1040 else:
1041 assert False
1043 def make_weights(self, receiver, scheme):
1044 self.check_scheme(scheme)
1046 azis, bazis = self.azibazis_to(receiver)
1048 sb = num.sin(bazis*d2r-num.pi)
1049 cb = num.cos(bazis*d2r-num.pi)
1051 m0s = self.m0s
1052 n = azis.size
1054 cat = num.concatenate
1055 rep = num.repeat
1057 if scheme == 'elastic2':
1058 w_n = cb*m0s
1059 g_n = filledi(0, n)
1060 w_e = sb*m0s
1061 g_e = filledi(0, n)
1062 w_d = m0s
1063 g_d = filledi(1, n)
1065 elif scheme == 'elastic8':
1066 w_n = cat((cb*m0s, cb*m0s))
1067 g_n = rep((0, 2), n)
1068 w_e = cat((sb*m0s, sb*m0s))
1069 g_e = rep((0, 2), n)
1070 w_d = cat((m0s, m0s))
1071 g_d = rep((5, 7), n)
1073 elif scheme == 'elastic10':
1074 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1075 g_n = rep((0, 2, 8), n)
1076 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1077 g_e = rep((0, 2, 8), n)
1078 w_d = cat((m0s, m0s, m0s))
1079 g_d = rep((5, 7, 9), n)
1081 else:
1082 assert False
1084 return (
1085 ('displacement.n', w_n, g_n),
1086 ('displacement.e', w_e, g_e),
1087 ('displacement.d', w_d, g_d),
1088 )
1090 def split(self):
1091 from pyrocko.gf.seismosizer import ExplosionSource
1092 sources = []
1093 for i in range(self.nelements):
1094 lat, lon, north_shift, east_shift = self.element_coords(i)
1095 sources.append(ExplosionSource(
1096 time=float(self.times[i]),
1097 lat=lat,
1098 lon=lon,
1099 north_shift=north_shift,
1100 east_shift=east_shift,
1101 depth=float(self.depths[i]),
1102 moment=float(self.m0s[i])))
1104 return sources
1106 def moments(self):
1107 return self.m0s
1109 def centroid(self):
1110 from pyrocko.gf.seismosizer import ExplosionSource
1111 time, lat, lon, north_shift, east_shift, depth = \
1112 self.centroid_position()
1114 return ExplosionSource(
1115 time=time,
1116 lat=lat,
1117 lon=lon,
1118 north_shift=north_shift,
1119 east_shift=east_shift,
1120 depth=depth,
1121 moment=float(num.sum(self.m0s)))
1123 @classmethod
1124 def combine(cls, sources, **kwargs):
1125 '''
1126 Combine several discretized source models.
1128 Concatenenates all point sources in the given discretized ``sources``.
1129 Care must be taken when using this function that the external amplitude
1130 factors and reference times of the parameterized (undiscretized)
1131 sources match or are accounted for.
1132 '''
1134 if 'm0s' not in kwargs:
1135 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1137 return super(DiscretizedExplosionSource, cls).combine(sources,
1138 **kwargs)
1141class DiscretizedSFSource(DiscretizedSource):
1142 forces = Array.T(shape=(None, 3), dtype=float)
1144 provided_schemes = (
1145 'elastic5',
1146 )
1148 def get_source_terms(self, scheme):
1149 self.check_scheme(scheme)
1151 return self.forces
1153 def make_weights(self, receiver, scheme):
1154 self.check_scheme(scheme)
1156 azis, bazis = self.azibazis_to(receiver)
1158 sa = num.sin(azis*d2r)
1159 ca = num.cos(azis*d2r)
1160 sb = num.sin(bazis*d2r-num.pi)
1161 cb = num.cos(bazis*d2r-num.pi)
1163 forces = self.forces
1164 fn = forces[:, 0]
1165 fe = forces[:, 1]
1166 fd = forces[:, 2]
1168 f0 = fd
1169 f1 = ca * fn + sa * fe
1170 f2 = ca * fe - sa * fn
1172 n = azis.size
1174 if scheme == 'elastic5':
1175 ioff = 0
1177 cat = num.concatenate
1178 rep = num.repeat
1180 w_n = cat((cb*f0, cb*f1, -sb*f2))
1181 g_n = ioff + rep((0, 1, 2), n)
1182 w_e = cat((sb*f0, sb*f1, cb*f2))
1183 g_e = ioff + rep((0, 1, 2), n)
1184 w_d = cat((f0, f1))
1185 g_d = ioff + rep((3, 4), n)
1187 return (
1188 ('displacement.n', w_n, g_n),
1189 ('displacement.e', w_e, g_e),
1190 ('displacement.d', w_d, g_d),
1191 )
1193 @classmethod
1194 def combine(cls, sources, **kwargs):
1195 '''
1196 Combine several discretized source models.
1198 Concatenenates all point sources in the given discretized ``sources``.
1199 Care must be taken when using this function that the external amplitude
1200 factors and reference times of the parameterized (undiscretized)
1201 sources match or are accounted for.
1202 '''
1204 if 'forces' not in kwargs:
1205 kwargs['forces'] = num.vstack([s.forces for s in sources])
1207 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1209 def moments(self):
1210 return num.sum(self.forces**2, axis=1)
1212 def centroid(self):
1213 from pyrocko.gf.seismosizer import SFSource
1214 time, lat, lon, north_shift, east_shift, depth = \
1215 self.centroid_position()
1217 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1218 return SFSource(
1219 time=time,
1220 lat=lat,
1221 lon=lon,
1222 north_shift=north_shift,
1223 east_shift=east_shift,
1224 depth=depth,
1225 fn=fn,
1226 fe=fe,
1227 fd=fd)
1230class DiscretizedMTSource(DiscretizedSource):
1231 m6s = Array.T(
1232 shape=(None, 6), dtype=float,
1233 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1235 provided_schemes = (
1236 'elastic8',
1237 'elastic10',
1238 'elastic18',
1239 )
1241 def get_source_terms(self, scheme):
1242 self.check_scheme(scheme)
1243 return self.m6s
1245 def make_weights(self, receiver, scheme):
1246 self.check_scheme(scheme)
1248 m6s = self.m6s
1249 n = m6s.shape[0]
1250 rep = num.repeat
1252 if scheme == 'elastic18':
1253 w_n = m6s.flatten()
1254 g_n = num.tile(num.arange(0, 6), n)
1255 w_e = m6s.flatten()
1256 g_e = num.tile(num.arange(6, 12), n)
1257 w_d = m6s.flatten()
1258 g_d = num.tile(num.arange(12, 18), n)
1260 else:
1261 azis, bazis = self.azibazis_to(receiver)
1263 sa = num.sin(azis*d2r)
1264 ca = num.cos(azis*d2r)
1265 s2a = num.sin(2.*azis*d2r)
1266 c2a = num.cos(2.*azis*d2r)
1267 sb = num.sin(bazis*d2r-num.pi)
1268 cb = num.cos(bazis*d2r-num.pi)
1270 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1271 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1272 f2 = m6s[:, 2]
1273 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1274 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1275 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1277 cat = num.concatenate
1279 if scheme == 'elastic8':
1280 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1281 g_n = rep((0, 1, 2, 3, 4), n)
1282 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1283 g_e = rep((0, 1, 2, 3, 4), n)
1284 w_d = cat((f0, f1, f2))
1285 g_d = rep((5, 6, 7), n)
1287 elif scheme == 'elastic10':
1288 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1289 g_n = rep((0, 1, 2, 8, 3, 4), n)
1290 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1291 g_e = rep((0, 1, 2, 8, 3, 4), n)
1292 w_d = cat((f0, f1, f2, f5))
1293 g_d = rep((5, 6, 7, 9), n)
1295 else:
1296 assert False
1298 return (
1299 ('displacement.n', w_n, g_n),
1300 ('displacement.e', w_e, g_e),
1301 ('displacement.d', w_d, g_d),
1302 )
1304 def split(self):
1305 from pyrocko.gf.seismosizer import MTSource
1306 sources = []
1307 for i in range(self.nelements):
1308 lat, lon, north_shift, east_shift = self.element_coords(i)
1309 sources.append(MTSource(
1310 time=float(self.times[i]),
1311 lat=lat,
1312 lon=lon,
1313 north_shift=north_shift,
1314 east_shift=east_shift,
1315 depth=float(self.depths[i]),
1316 m6=self.m6s[i]))
1318 return sources
1320 def moments(self):
1321 moments = num.array(
1322 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))
1323 for m6 in self.m6s])
1324 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)
1326 def get_moment_rate(self, deltat=None):
1327 moments = self.moments()
1328 times = self.times
1329 times -= times.min()
1331 t_max = times.max()
1332 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat
1333 mom_times[mom_times > t_max] = t_max
1335 # Right open histrogram bins
1336 mom, _ = num.histogram(
1337 times,
1338 bins=mom_times,
1339 weights=moments)
1341 deltat = num.diff(mom_times)
1342 mom_rate = mom / deltat
1344 return mom_rate, mom_times[1:]
1346 def centroid(self):
1347 from pyrocko.gf.seismosizer import MTSource
1348 time, lat, lon, north_shift, east_shift, depth = \
1349 self.centroid_position()
1351 return MTSource(
1352 time=time,
1353 lat=lat,
1354 lon=lon,
1355 north_shift=north_shift,
1356 east_shift=east_shift,
1357 depth=depth,
1358 m6=num.sum(self.m6s, axis=0))
1360 @classmethod
1361 def combine(cls, sources, **kwargs):
1362 '''
1363 Combine several discretized source models.
1365 Concatenenates all point sources in the given discretized ``sources``.
1366 Care must be taken when using this function that the external amplitude
1367 factors and reference times of the parameterized (undiscretized)
1368 sources match or are accounted for.
1369 '''
1371 if 'm6s' not in kwargs:
1372 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1374 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1377class DiscretizedPorePressureSource(DiscretizedSource):
1378 pp = Array.T(shape=(None,), dtype=float)
1380 provided_schemes = (
1381 'poroelastic10',
1382 )
1384 def get_source_terms(self, scheme):
1385 self.check_scheme(scheme)
1386 return self.pp[:, num.newaxis].copy()
1388 def make_weights(self, receiver, scheme):
1389 self.check_scheme(scheme)
1391 azis, bazis = self.azibazis_to(receiver)
1393 sb = num.sin(bazis*d2r-num.pi)
1394 cb = num.cos(bazis*d2r-num.pi)
1396 pp = self.pp
1397 n = bazis.size
1399 w_un = cb*pp
1400 g_un = filledi(1, n)
1401 w_ue = sb*pp
1402 g_ue = filledi(1, n)
1403 w_ud = pp
1404 g_ud = filledi(0, n)
1406 w_tn = cb*pp
1407 g_tn = filledi(6, n)
1408 w_te = sb*pp
1409 g_te = filledi(6, n)
1411 w_pp = pp
1412 g_pp = filledi(7, n)
1414 w_dvn = cb*pp
1415 g_dvn = filledi(9, n)
1416 w_dve = sb*pp
1417 g_dve = filledi(9, n)
1418 w_dvd = pp
1419 g_dvd = filledi(8, n)
1421 return (
1422 ('displacement.n', w_un, g_un),
1423 ('displacement.e', w_ue, g_ue),
1424 ('displacement.d', w_ud, g_ud),
1425 ('vertical_tilt.n', w_tn, g_tn),
1426 ('vertical_tilt.e', w_te, g_te),
1427 ('pore_pressure', w_pp, g_pp),
1428 ('darcy_velocity.n', w_dvn, g_dvn),
1429 ('darcy_velocity.e', w_dve, g_dve),
1430 ('darcy_velocity.d', w_dvd, g_dvd),
1431 )
1433 def moments(self):
1434 return self.pp
1436 def centroid(self):
1437 from pyrocko.gf.seismosizer import PorePressurePointSource
1438 time, lat, lon, north_shift, east_shift, depth = \
1439 self.centroid_position()
1441 return PorePressurePointSource(
1442 time=time,
1443 lat=lat,
1444 lon=lon,
1445 north_shift=north_shift,
1446 east_shift=east_shift,
1447 depth=depth,
1448 pp=float(num.sum(self.pp)))
1450 @classmethod
1451 def combine(cls, sources, **kwargs):
1452 '''
1453 Combine several discretized source models.
1455 Concatenenates all point sources in the given discretized ``sources``.
1456 Care must be taken when using this function that the external amplitude
1457 factors and reference times of the parameterized (undiscretized)
1458 sources match or are accounted for.
1459 '''
1461 if 'pp' not in kwargs:
1462 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1464 return super(DiscretizedPorePressureSource, cls).combine(sources,
1465 **kwargs)
1468class Region(Object):
1469 name = String.T(optional=True)
1472class RectangularRegion(Region):
1473 lat_min = Float.T()
1474 lat_max = Float.T()
1475 lon_min = Float.T()
1476 lon_max = Float.T()
1479class CircularRegion(Region):
1480 lat = Float.T()
1481 lon = Float.T()
1482 radius = Float.T()
1485class Config(Object):
1486 '''
1487 Green's function store meta information.
1489 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1490 configuration types are:
1492 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1493 spherical symmetry, 1D earth model, single receiver depth
1495 * Problem is invariant to horizontal translations and rotations around
1496 vertical axis.
1497 * All receivers must be at the same depth (e.g. at the surface)
1498 * High level index variables: ``(source_depth, receiver_distance,
1499 component)``
1501 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1502 spherical symmetry, 1D earth model, variable receiver depth
1504 * Symmetries like in Type A but has additional index for receiver depth
1505 * High level index variables: ``(source_depth, receiver_distance,
1506 receiver_depth, component)``
1508 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1509 constraints but fixed receiver positions
1511 * Cartesian source volume around a reference point
1512 * High level index variables: ``(ireceiver, source_depth,
1513 source_east_shift, source_north_shift, component)``
1514 '''
1516 id = StringID.T(
1517 help='Name of the store. May consist of upper and lower-case letters, '
1518 'digits, dots and underscores. The name must start with a '
1519 'letter.')
1521 derived_from_id = StringID.T(
1522 optional=True,
1523 help='Name of the original store, if this store has been derived from '
1524 'another one (e.g. extracted subset).')
1526 version = String.T(
1527 default='1.0',
1528 optional=True,
1529 help='User-defined version string. Use <major>.<minor> format.')
1531 modelling_code_id = StringID.T(
1532 optional=True,
1533 help='Identifier of the backend used to compute the store.')
1535 author = Unicode.T(
1536 optional=True,
1537 help='Comma-separated list of author names.')
1539 author_email = String.T(
1540 optional=True,
1541 help="Author's contact email address.")
1543 created_time = Timestamp.T(
1544 optional=True,
1545 help='Time of creation of the store.')
1547 regions = List.T(
1548 Region.T(),
1549 help='Geographical regions for which the store is representative.')
1551 scope_type = ScopeType.T(
1552 optional=True,
1553 help='Distance range scope of the store '
1554 '(%s).' % fmt_choices(ScopeType))
1556 waveform_type = WaveType.T(
1557 optional=True,
1558 help='Wave type stored (%s).' % fmt_choices(WaveType))
1560 nearfield_terms = NearfieldTermsType.T(
1561 optional=True,
1562 help='Information about the inclusion of near-field terms in the '
1563 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1565 description = String.T(
1566 optional=True,
1567 help='Free form textual description of the GF store.')
1569 references = List.T(
1570 Reference.T(),
1571 help='Reference list to cite the modelling code, earth model or '
1572 'related work.')
1574 earthmodel_1d = Earthmodel1D.T(
1575 optional=True,
1576 help='Layered earth model in ND (named discontinuity) format.')
1578 earthmodel_receiver_1d = Earthmodel1D.T(
1579 optional=True,
1580 help='Receiver-side layered earth model in ND format.')
1582 can_interpolate_source = Bool.T(
1583 optional=True,
1584 help='Hint to indicate if the spatial sampling of the store is dense '
1585 'enough for multi-linear interpolation at the source.')
1587 can_interpolate_receiver = Bool.T(
1588 optional=True,
1589 help='Hint to indicate if the spatial sampling of the store is dense '
1590 'enough for multi-linear interpolation at the receiver.')
1592 frequency_min = Float.T(
1593 optional=True,
1594 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1596 frequency_max = Float.T(
1597 optional=True,
1598 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1600 sample_rate = Float.T(
1601 optional=True,
1602 help='Sample rate of the GF store [Hz].')
1604 factor = Float.T(
1605 default=1.0,
1606 help='Gain value, factored out of the stored GF samples. '
1607 '(may not work properly, keep at 1.0).',
1608 optional=True)
1610 component_scheme = ComponentScheme.T(
1611 default='elastic10',
1612 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1614 stored_quantity = QuantityType.T(
1615 optional=True,
1616 help='Physical quantity of stored values (%s). If not given, a '
1617 'default is used based on the GF component scheme. The default '
1618 'for the ``"elastic*"`` family of component schemes is '
1619 '``"displacement"``.' % fmt_choices(QuantityType))
1621 tabulated_phases = List.T(
1622 TPDef.T(),
1623 help='Mapping of phase names to phase definitions, for which travel '
1624 'time tables are available in the GF store.')
1626 ncomponents = Int.T(
1627 optional=True,
1628 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1630 uuid = String.T(
1631 optional=True,
1632 help='Heuristic hash value which can be used to uniquely identify the '
1633 'GF store for practical purposes.')
1635 reference = String.T(
1636 optional=True,
1637 help="Store reference name composed of the store's :gattr:`id` and "
1638 'the first six letters of its :gattr:`uuid`.')
1640 def __init__(self, **kwargs):
1641 self._do_auto_updates = False
1642 Object.__init__(self, **kwargs)
1643 self._index_function = None
1644 self._indices_function = None
1645 self._vicinity_function = None
1646 self.validate(regularize=True, depth=1)
1647 self._do_auto_updates = True
1648 self.update()
1650 def check_ncomponents(self):
1651 ncomponents = component_scheme_to_description[
1652 self.component_scheme].ncomponents
1654 if self.ncomponents is None:
1655 self.ncomponents = ncomponents
1656 elif ncomponents != self.ncomponents:
1657 raise InvalidNComponents(
1658 'ncomponents=%i incompatible with component_scheme="%s"' % (
1659 self.ncomponents, self.component_scheme))
1661 def __setattr__(self, name, value):
1662 Object.__setattr__(self, name, value)
1663 try:
1664 self.T.get_property(name)
1665 if self._do_auto_updates:
1666 self.update()
1668 except ValueError:
1669 pass
1671 def update(self):
1672 self.check_ncomponents()
1673 self._update()
1674 self._make_index_functions()
1676 def irecord(self, *args):
1677 return self._index_function(*args)
1679 def irecords(self, *args):
1680 return self._indices_function(*args)
1682 def vicinity(self, *args):
1683 return self._vicinity_function(*args)
1685 def vicinities(self, *args):
1686 return self._vicinities_function(*args)
1688 def grid_interpolation_coefficients(self, *args):
1689 return self._grid_interpolation_coefficients(*args)
1691 def nodes(self, level=None, minlevel=None):
1692 return nodes(self.coords[minlevel:level])
1694 def iter_nodes(self, level=None, minlevel=None):
1695 return nditer_outer(self.coords[minlevel:level])
1697 def iter_extraction(self, gdef, level=None):
1698 i = 0
1699 arrs = []
1700 ntotal = 1
1701 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1702 if gdef and len(gdef) > i:
1703 sssn = gdef[i]
1704 else:
1705 sssn = (None,)*4
1707 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1708 ntotal *= len(arr)
1710 arrs.append(arr)
1711 i += 1
1713 arrs.append(self.coords[-1])
1714 return nditer_outer(arrs[:level])
1716 def make_sum_params(self, source, receiver, implementation='c',
1717 nthreads=0):
1718 assert implementation in ['c', 'python']
1720 out = []
1721 delays = source.times
1722 for comp, weights, icomponents in source.make_weights(
1723 receiver,
1724 self.component_scheme):
1726 weights *= self.factor
1728 args = self.make_indexing_args(source, receiver, icomponents)
1729 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1730 out.append((comp, args, delays_expanded, weights))
1732 return out
1734 def short_info(self):
1735 raise NotImplementedError('should be implemented in subclass')
1737 def get_shear_moduli(self, lat, lon, points,
1738 interpolation=None):
1739 '''
1740 Get shear moduli at given points from contained velocity model.
1742 :param lat: surface origin for coordinate system of ``points``
1743 :param points: NumPy array of shape ``(N, 3)``, where each row is
1744 a point ``(north, east, depth)``, relative to origin at
1745 ``(lat, lon)``
1746 :param interpolation: interpolation method. Choose from
1747 ``('nearest_neighbor', 'multilinear')``
1748 :returns: NumPy array of length N with extracted shear moduli at each
1749 point
1751 The default implementation retrieves and interpolates the shear moduli
1752 from the contained 1D velocity profile.
1753 '''
1754 return self.get_material_property(lat, lon, points,
1755 parameter='shear_moduli',
1756 interpolation=interpolation)
1758 def get_lambda_moduli(self, lat, lon, points,
1759 interpolation=None):
1760 '''
1761 Get lambda moduli at given points from contained velocity model.
1763 :param lat: surface origin for coordinate system of ``points``
1764 :param points: NumPy array of shape ``(N, 3)``, where each row is
1765 a point ``(north, east, depth)``, relative to origin at
1766 ``(lat, lon)``
1767 :param interpolation: interpolation method. Choose from
1768 ``('nearest_neighbor', 'multilinear')``
1769 :returns: NumPy array of length N with extracted shear moduli at each
1770 point
1772 The default implementation retrieves and interpolates the lambda moduli
1773 from the contained 1D velocity profile.
1774 '''
1775 return self.get_material_property(lat, lon, points,
1776 parameter='lambda_moduli',
1777 interpolation=interpolation)
1779 def get_bulk_moduli(self, lat, lon, points,
1780 interpolation=None):
1781 '''
1782 Get bulk moduli at given points from contained velocity model.
1784 :param lat: surface origin for coordinate system of ``points``
1785 :param points: NumPy array of shape ``(N, 3)``, where each row is
1786 a point ``(north, east, depth)``, relative to origin at
1787 ``(lat, lon)``
1788 :param interpolation: interpolation method. Choose from
1789 ``('nearest_neighbor', 'multilinear')``
1790 :returns: NumPy array of length N with extracted shear moduli at each
1791 point
1793 The default implementation retrieves and interpolates the lambda moduli
1794 from the contained 1D velocity profile.
1795 '''
1796 lambda_moduli = self.get_material_property(
1797 lat, lon, points, parameter='lambda_moduli',
1798 interpolation=interpolation)
1799 shear_moduli = self.get_material_property(
1800 lat, lon, points, parameter='shear_moduli',
1801 interpolation=interpolation)
1802 return lambda_moduli + (2 / 3) * shear_moduli
1804 def get_vs(self, lat, lon, points, interpolation=None):
1805 '''
1806 Get Vs at given points from contained velocity model.
1808 :param lat: surface origin for coordinate system of ``points``
1809 :param points: NumPy array of shape ``(N, 3)``, where each row is
1810 a point ``(north, east, depth)``, relative to origin at
1811 ``(lat, lon)``
1812 :param interpolation: interpolation method. Choose from
1813 ``('nearest_neighbor', 'multilinear')``
1814 :returns: NumPy array of length N with extracted shear moduli at each
1815 point
1817 The default implementation retrieves and interpolates Vs
1818 from the contained 1D velocity profile.
1819 '''
1820 return self.get_material_property(lat, lon, points,
1821 parameter='vs',
1822 interpolation=interpolation)
1824 def get_vp(self, lat, lon, points, interpolation=None):
1825 '''
1826 Get Vp at given points from contained velocity model.
1828 :param lat: surface origin for coordinate system of ``points``
1829 :param points: NumPy array of shape ``(N, 3)``, where each row is
1830 a point ``(north, east, depth)``, relative to origin at
1831 ``(lat, lon)``
1832 :param interpolation: interpolation method. Choose from
1833 ``('nearest_neighbor', 'multilinear')``
1834 :returns: NumPy array of length N with extracted shear moduli at each
1835 point
1837 The default implementation retrieves and interpolates Vp
1838 from the contained 1D velocity profile.
1839 '''
1840 return self.get_material_property(lat, lon, points,
1841 parameter='vp',
1842 interpolation=interpolation)
1844 def get_rho(self, lat, lon, points, interpolation=None):
1845 '''
1846 Get rho at given points from contained velocity model.
1848 :param lat: surface origin for coordinate system of ``points``
1849 :param points: NumPy array of shape ``(N, 3)``, where each row is
1850 a point ``(north, east, depth)``, relative to origin at
1851 ``(lat, lon)``
1852 :param interpolation: interpolation method. Choose from
1853 ``('nearest_neighbor', 'multilinear')``
1854 :returns: NumPy array of length N with extracted shear moduli at each
1855 point
1857 The default implementation retrieves and interpolates rho
1858 from the contained 1D velocity profile.
1859 '''
1860 return self.get_material_property(lat, lon, points,
1861 parameter='rho',
1862 interpolation=interpolation)
1864 def get_material_property(self, lat, lon, points, parameter='vs',
1865 interpolation=None):
1867 if interpolation is None:
1868 raise TypeError('Interpolation method not defined! available: '
1869 'multilinear', 'nearest_neighbor')
1871 earthmod = self.earthmodel_1d
1872 store_depth_profile = self.get_source_depths()
1873 z_profile = earthmod.profile('z')
1875 if parameter == 'vs':
1876 vs_profile = earthmod.profile('vs')
1877 profile = num.interp(
1878 store_depth_profile, z_profile, vs_profile)
1880 elif parameter == 'vp':
1881 vp_profile = earthmod.profile('vp')
1882 profile = num.interp(
1883 store_depth_profile, z_profile, vp_profile)
1885 elif parameter == 'rho':
1886 rho_profile = earthmod.profile('rho')
1888 profile = num.interp(
1889 store_depth_profile, z_profile, rho_profile)
1891 elif parameter == 'shear_moduli':
1892 vs_profile = earthmod.profile('vs')
1893 rho_profile = earthmod.profile('rho')
1895 store_vs_profile = num.interp(
1896 store_depth_profile, z_profile, vs_profile)
1897 store_rho_profile = num.interp(
1898 store_depth_profile, z_profile, rho_profile)
1900 profile = num.power(store_vs_profile, 2) * store_rho_profile
1902 elif parameter == 'lambda_moduli':
1903 vs_profile = earthmod.profile('vs')
1904 vp_profile = earthmod.profile('vp')
1905 rho_profile = earthmod.profile('rho')
1907 store_vs_profile = num.interp(
1908 store_depth_profile, z_profile, vs_profile)
1909 store_vp_profile = num.interp(
1910 store_depth_profile, z_profile, vp_profile)
1911 store_rho_profile = num.interp(
1912 store_depth_profile, z_profile, rho_profile)
1914 profile = store_rho_profile * (
1915 num.power(store_vp_profile, 2) -
1916 num.power(store_vs_profile, 2) * 2)
1917 else:
1918 raise TypeError(
1919 'parameter %s not available' % parameter)
1921 if interpolation == 'multilinear':
1922 kind = 'linear'
1923 elif interpolation == 'nearest_neighbor':
1924 kind = 'nearest'
1925 else:
1926 raise TypeError(
1927 'Interpolation method %s not available' % interpolation)
1929 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1931 try:
1932 return interpolator(points[:, 2])
1933 except ValueError:
1934 raise OutOfBounds()
1936 def is_static(self):
1937 for code in ('psgrn_pscmp', 'poel'):
1938 if self.modelling_code_id.startswith(code):
1939 return True
1940 return False
1942 def is_dynamic(self):
1943 return not self.is_static()
1945 def get_source_depths(self):
1946 raise NotImplementedError('must be implemented in subclass')
1948 def get_tabulated_phase(self, phase_id):
1949 '''
1950 Get tabulated phase definition.
1951 '''
1953 for pdef in self.tabulated_phases:
1954 if pdef.id == phase_id:
1955 return pdef
1957 raise StoreError('No such phase: %s' % phase_id)
1959 def fix_ttt_holes(self, sptree, mode):
1960 raise StoreError(
1961 'Cannot fix travel time table holes in GF stores of type %s.'
1962 % self.short_type)
1965class ConfigTypeA(Config):
1966 '''
1967 Cylindrical symmetry, 1D earth model, single receiver depth
1969 * Problem is invariant to horizontal translations and rotations around
1970 vertical axis.
1972 * All receivers must be at the same depth (e.g. at the surface)
1973 High level index variables: ``(source_depth, distance,
1974 component)``
1976 * The ``distance`` is the surface distance between source and receiver
1977 points.
1978 '''
1980 receiver_depth = Float.T(
1981 default=0.0,
1982 help='Fixed receiver depth [m].')
1984 source_depth_min = Float.T(
1985 help='Minimum source depth [m].')
1987 source_depth_max = Float.T(
1988 help='Maximum source depth [m].')
1990 source_depth_delta = Float.T(
1991 help='Grid spacing of source depths [m]')
1993 distance_min = Float.T(
1994 help='Minimum source-receiver surface distance [m].')
1996 distance_max = Float.T(
1997 help='Maximum source-receiver surface distance [m].')
1999 distance_delta = Float.T(
2000 help='Grid spacing of source-receiver surface distance [m].')
2002 short_type = 'A'
2004 provided_schemes = [
2005 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2007 def get_surface_distance(self, args):
2008 return args[1]
2010 def get_distance(self, args):
2011 return math.sqrt(args[0]**2 + args[1]**2)
2013 def get_source_depth(self, args):
2014 return args[0]
2016 def get_source_depths(self):
2017 return self.coords[0]
2019 def get_receiver_depth(self, args):
2020 return self.receiver_depth
2022 def _update(self):
2023 self.mins = num.array(
2024 [self.source_depth_min, self.distance_min], dtype=float)
2025 self.maxs = num.array(
2026 [self.source_depth_max, self.distance_max], dtype=float)
2027 self.deltas = num.array(
2028 [self.source_depth_delta, self.distance_delta],
2029 dtype=float)
2030 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2031 vicinity_eps).astype(int) + 1
2032 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2033 self.deltat = 1.0/self.sample_rate
2034 self.nrecords = num.product(self.ns) * self.ncomponents
2035 self.coords = tuple(num.linspace(mi, ma, n) for
2036 (mi, ma, n) in
2037 zip(self.mins, self.effective_maxs, self.ns)) + \
2038 (num.arange(self.ncomponents),)
2040 self.nsource_depths, self.ndistances = self.ns
2042 def _make_index_functions(self):
2044 amin, bmin = self.mins
2045 da, db = self.deltas
2046 na, nb = self.ns
2048 ng = self.ncomponents
2050 def index_function(a, b, ig):
2051 ia = int(round((a - amin) / da))
2052 ib = int(round((b - bmin) / db))
2053 try:
2054 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2055 except ValueError:
2056 raise OutOfBounds()
2058 def indices_function(a, b, ig):
2059 ia = num.round((a - amin) / da).astype(int)
2060 ib = num.round((b - bmin) / db).astype(int)
2061 try:
2062 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2063 except ValueError:
2064 for ia_, ib_, ig_ in zip(ia, ib, ig):
2065 try:
2066 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2067 except ValueError:
2068 raise OutOfBounds()
2070 def grid_interpolation_coefficients(a, b):
2071 ias = indi12((a - amin) / da, na)
2072 ibs = indi12((b - bmin) / db, nb)
2073 return ias, ibs
2075 def vicinity_function(a, b, ig):
2076 ias, ibs = grid_interpolation_coefficients(a, b)
2078 if not (0 <= ig < ng):
2079 raise OutOfBounds()
2081 indis = []
2082 weights = []
2083 for ia, va in ias:
2084 iia = ia*nb*ng
2085 for ib, vb in ibs:
2086 indis.append(iia + ib*ng + ig)
2087 weights.append(va*vb)
2089 return num.array(indis), num.array(weights)
2091 def vicinities_function(a, b, ig):
2093 xa = (a - amin) / da
2094 xb = (b - bmin) / db
2096 xa_fl = num.floor(xa)
2097 xa_ce = num.ceil(xa)
2098 xb_fl = num.floor(xb)
2099 xb_ce = num.ceil(xb)
2100 va_fl = 1.0 - (xa - xa_fl)
2101 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2102 vb_fl = 1.0 - (xb - xb_fl)
2103 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2105 ia_fl = xa_fl.astype(int)
2106 ia_ce = xa_ce.astype(int)
2107 ib_fl = xb_fl.astype(int)
2108 ib_ce = xb_ce.astype(int)
2110 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2111 raise OutOfBounds()
2113 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2114 raise OutOfBounds()
2116 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2117 raise OutOfBounds()
2119 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2120 raise OutOfBounds()
2122 irecords = num.empty(a.size*4, dtype=int)
2123 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2124 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2125 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2126 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2128 weights = num.empty(a.size*4, dtype=float)
2129 weights[0::4] = va_fl * vb_fl
2130 weights[1::4] = va_ce * vb_fl
2131 weights[2::4] = va_fl * vb_ce
2132 weights[3::4] = va_ce * vb_ce
2134 return irecords, weights
2136 self._index_function = index_function
2137 self._indices_function = indices_function
2138 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2139 self._vicinity_function = vicinity_function
2140 self._vicinities_function = vicinities_function
2142 def make_indexing_args(self, source, receiver, icomponents):
2143 nc = icomponents.size
2144 dists = source.distances_to(receiver)
2145 n = dists.size
2146 return (num.tile(source.depths, nc//n),
2147 num.tile(dists, nc//n),
2148 icomponents)
2150 def make_indexing_args1(self, source, receiver):
2151 return (source.depth, source.distance_to(receiver))
2153 @property
2154 def short_extent(self):
2155 return '%g:%g:%g x %g:%g:%g' % (
2156 self.source_depth_min/km,
2157 self.source_depth_max/km,
2158 self.source_depth_delta/km,
2159 self.distance_min/km,
2160 self.distance_max/km,
2161 self.distance_delta/km)
2163 def fix_ttt_holes(self, sptree, mode):
2164 from pyrocko import eikonal_ext, spit
2166 nodes = self.nodes(level=-1)
2168 delta = self.deltas[-1]
2169 assert num.all(delta == self.deltas)
2171 nsources, ndistances = self.ns
2173 points = num.zeros((nodes.shape[0], 3))
2174 points[:, 0] = nodes[:, 1]
2175 points[:, 2] = nodes[:, 0]
2177 speeds = self.get_material_property(
2178 0., 0., points,
2179 parameter='vp' if mode == cake.P else 'vs',
2180 interpolation='multilinear')
2182 speeds = speeds.reshape((nsources, ndistances))
2184 times = sptree.interpolate_many(nodes)
2186 times[num.isnan(times)] = -1.
2187 times = times.reshape(speeds.shape)
2189 try:
2190 eikonal_ext.eikonal_solver_fmm_cartesian(
2191 speeds, times, delta)
2192 except eikonal_ext.EikonalExtError as e:
2193 if str(e).endswith('please check results'):
2194 logger.debug(
2195 'Got a warning from eikonal solver '
2196 '- may be ok...')
2197 else:
2198 raise
2200 def func(x):
2201 ibs, ics = \
2202 self.grid_interpolation_coefficients(*x)
2204 t = 0
2205 for ib, vb in ibs:
2206 for ic, vc in ics:
2207 t += times[ib, ic] * vb * vc
2209 return t
2211 return spit.SPTree(
2212 f=func,
2213 ftol=sptree.ftol,
2214 xbounds=sptree.xbounds,
2215 xtols=sptree.xtols)
2218class ConfigTypeB(Config):
2219 '''
2220 Cylindrical symmetry, 1D earth model, variable receiver depth
2222 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2223 receiver depth
2225 * High level index variables: ``(receiver_depth, source_depth,
2226 receiver_distance, component)``
2227 '''
2229 receiver_depth_min = Float.T(
2230 help='Minimum receiver depth [m].')
2232 receiver_depth_max = Float.T(
2233 help='Maximum receiver depth [m].')
2235 receiver_depth_delta = Float.T(
2236 help='Grid spacing of receiver depths [m]')
2238 source_depth_min = Float.T(
2239 help='Minimum source depth [m].')
2241 source_depth_max = Float.T(
2242 help='Maximum source depth [m].')
2244 source_depth_delta = Float.T(
2245 help='Grid spacing of source depths [m]')
2247 distance_min = Float.T(
2248 help='Minimum source-receiver surface distance [m].')
2250 distance_max = Float.T(
2251 help='Maximum source-receiver surface distance [m].')
2253 distance_delta = Float.T(
2254 help='Grid spacing of source-receiver surface distances [m].')
2256 short_type = 'B'
2258 provided_schemes = [
2259 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2261 def get_distance(self, args):
2262 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2264 def get_surface_distance(self, args):
2265 return args[2]
2267 def get_source_depth(self, args):
2268 return args[1]
2270 def get_receiver_depth(self, args):
2271 return args[0]
2273 def get_source_depths(self):
2274 return self.coords[1]
2276 def _update(self):
2277 self.mins = num.array([
2278 self.receiver_depth_min,
2279 self.source_depth_min,
2280 self.distance_min],
2281 dtype=float)
2283 self.maxs = num.array([
2284 self.receiver_depth_max,
2285 self.source_depth_max,
2286 self.distance_max],
2287 dtype=float)
2289 self.deltas = num.array([
2290 self.receiver_depth_delta,
2291 self.source_depth_delta,
2292 self.distance_delta],
2293 dtype=float)
2295 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2296 vicinity_eps).astype(int) + 1
2297 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2298 self.deltat = 1.0/self.sample_rate
2299 self.nrecords = num.product(self.ns) * self.ncomponents
2300 self.coords = tuple(num.linspace(mi, ma, n) for
2301 (mi, ma, n) in
2302 zip(self.mins, self.effective_maxs, self.ns)) + \
2303 (num.arange(self.ncomponents),)
2304 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2306 def _make_index_functions(self):
2308 amin, bmin, cmin = self.mins
2309 da, db, dc = self.deltas
2310 na, nb, nc = self.ns
2311 ng = self.ncomponents
2313 def index_function(a, b, c, ig):
2314 ia = int(round((a - amin) / da))
2315 ib = int(round((b - bmin) / db))
2316 ic = int(round((c - cmin) / dc))
2317 try:
2318 return num.ravel_multi_index((ia, ib, ic, ig),
2319 (na, nb, nc, ng))
2320 except ValueError:
2321 raise OutOfBounds()
2323 def indices_function(a, b, c, ig):
2324 ia = num.round((a - amin) / da).astype(int)
2325 ib = num.round((b - bmin) / db).astype(int)
2326 ic = num.round((c - cmin) / dc).astype(int)
2327 try:
2328 return num.ravel_multi_index((ia, ib, ic, ig),
2329 (na, nb, nc, ng))
2330 except ValueError:
2331 raise OutOfBounds()
2333 def grid_interpolation_coefficients(a, b, c):
2334 ias = indi12((a - amin) / da, na)
2335 ibs = indi12((b - bmin) / db, nb)
2336 ics = indi12((c - cmin) / dc, nc)
2337 return ias, ibs, ics
2339 def vicinity_function(a, b, c, ig):
2340 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2342 if not (0 <= ig < ng):
2343 raise OutOfBounds()
2345 indis = []
2346 weights = []
2347 for ia, va in ias:
2348 iia = ia*nb*nc*ng
2349 for ib, vb in ibs:
2350 iib = ib*nc*ng
2351 for ic, vc in ics:
2352 indis.append(iia + iib + ic*ng + ig)
2353 weights.append(va*vb*vc)
2355 return num.array(indis), num.array(weights)
2357 def vicinities_function(a, b, c, ig):
2359 xa = (a - amin) / da
2360 xb = (b - bmin) / db
2361 xc = (c - cmin) / dc
2363 xa_fl = num.floor(xa)
2364 xa_ce = num.ceil(xa)
2365 xb_fl = num.floor(xb)
2366 xb_ce = num.ceil(xb)
2367 xc_fl = num.floor(xc)
2368 xc_ce = num.ceil(xc)
2369 va_fl = 1.0 - (xa - xa_fl)
2370 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2371 vb_fl = 1.0 - (xb - xb_fl)
2372 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2373 vc_fl = 1.0 - (xc - xc_fl)
2374 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2376 ia_fl = xa_fl.astype(int)
2377 ia_ce = xa_ce.astype(int)
2378 ib_fl = xb_fl.astype(int)
2379 ib_ce = xb_ce.astype(int)
2380 ic_fl = xc_fl.astype(int)
2381 ic_ce = xc_ce.astype(int)
2383 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2384 raise OutOfBounds()
2386 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2387 raise OutOfBounds()
2389 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2390 raise OutOfBounds()
2392 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2393 raise OutOfBounds()
2395 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2396 raise OutOfBounds()
2398 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2399 raise OutOfBounds()
2401 irecords = num.empty(a.size*8, dtype=int)
2402 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2403 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2404 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2405 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2406 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2407 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2408 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2409 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2411 weights = num.empty(a.size*8, dtype=float)
2412 weights[0::8] = va_fl * vb_fl * vc_fl
2413 weights[1::8] = va_ce * vb_fl * vc_fl
2414 weights[2::8] = va_fl * vb_ce * vc_fl
2415 weights[3::8] = va_ce * vb_ce * vc_fl
2416 weights[4::8] = va_fl * vb_fl * vc_ce
2417 weights[5::8] = va_ce * vb_fl * vc_ce
2418 weights[6::8] = va_fl * vb_ce * vc_ce
2419 weights[7::8] = va_ce * vb_ce * vc_ce
2421 return irecords, weights
2423 self._index_function = index_function
2424 self._indices_function = indices_function
2425 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2426 self._vicinity_function = vicinity_function
2427 self._vicinities_function = vicinities_function
2429 def make_indexing_args(self, source, receiver, icomponents):
2430 nc = icomponents.size
2431 dists = source.distances_to(receiver)
2432 n = dists.size
2433 receiver_depths = num.empty(nc)
2434 receiver_depths.fill(receiver.depth)
2435 return (receiver_depths,
2436 num.tile(source.depths, nc//n),
2437 num.tile(dists, nc//n),
2438 icomponents)
2440 def make_indexing_args1(self, source, receiver):
2441 return (receiver.depth,
2442 source.depth,
2443 source.distance_to(receiver))
2445 @property
2446 def short_extent(self):
2447 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2448 self.receiver_depth_min/km,
2449 self.receiver_depth_max/km,
2450 self.receiver_depth_delta/km,
2451 self.source_depth_min/km,
2452 self.source_depth_max/km,
2453 self.source_depth_delta/km,
2454 self.distance_min/km,
2455 self.distance_max/km,
2456 self.distance_delta/km)
2458 def fix_ttt_holes(self, sptree, mode):
2459 from pyrocko import eikonal_ext, spit
2461 nodes_sr = self.nodes(minlevel=1, level=-1)
2463 delta = self.deltas[-1]
2464 assert num.all(delta == self.deltas[1:])
2466 nreceivers, nsources, ndistances = self.ns
2468 points = num.zeros((nodes_sr.shape[0], 3))
2469 points[:, 0] = nodes_sr[:, 1]
2470 points[:, 2] = nodes_sr[:, 0]
2472 speeds = self.get_material_property(
2473 0., 0., points,
2474 parameter='vp' if mode == cake.P else 'vs',
2475 interpolation='multilinear')
2477 speeds = speeds.reshape((nsources, ndistances))
2479 receiver_times = []
2480 for ireceiver in range(nreceivers):
2481 nodes = num.hstack([
2482 num_full(
2483 (nodes_sr.shape[0], 1),
2484 self.coords[0][ireceiver]),
2485 nodes_sr])
2487 times = sptree.interpolate_many(nodes)
2489 times[num.isnan(times)] = -1.
2491 times = times.reshape(speeds.shape)
2493 try:
2494 eikonal_ext.eikonal_solver_fmm_cartesian(
2495 speeds, times, delta)
2496 except eikonal_ext.EikonalExtError as e:
2497 if str(e).endswith('please check results'):
2498 logger.debug(
2499 'Got a warning from eikonal solver '
2500 '- may be ok...')
2501 else:
2502 raise
2504 receiver_times.append(times)
2506 def func(x):
2507 ias, ibs, ics = \
2508 self.grid_interpolation_coefficients(*x)
2510 t = 0
2511 for ia, va in ias:
2512 times = receiver_times[ia]
2513 for ib, vb in ibs:
2514 for ic, vc in ics:
2515 t += times[ib, ic] * va * vb * vc
2517 return t
2519 return spit.SPTree(
2520 f=func,
2521 ftol=sptree.ftol,
2522 xbounds=sptree.xbounds,
2523 xtols=sptree.xtols)
2526class ConfigTypeC(Config):
2527 '''
2528 No symmetrical constraints, one fixed receiver position.
2530 * Cartesian 3D source volume around a reference point
2532 * High level index variables: ``(source_depth,
2533 source_east_shift, source_north_shift, component)``
2534 '''
2536 receiver = Receiver.T(
2537 help='Receiver location')
2539 source_origin = Location.T(
2540 help='Origin of the source volume grid.')
2542 source_depth_min = Float.T(
2543 help='Minimum source depth [m].')
2545 source_depth_max = Float.T(
2546 help='Maximum source depth [m].')
2548 source_depth_delta = Float.T(
2549 help='Source depth grid spacing [m].')
2551 source_east_shift_min = Float.T(
2552 help='Minimum easting of source grid [m].')
2554 source_east_shift_max = Float.T(
2555 help='Maximum easting of source grid [m].')
2557 source_east_shift_delta = Float.T(
2558 help='Source volume grid spacing in east direction [m].')
2560 source_north_shift_min = Float.T(
2561 help='Minimum northing of source grid [m].')
2563 source_north_shift_max = Float.T(
2564 help='Maximum northing of source grid [m].')
2566 source_north_shift_delta = Float.T(
2567 help='Source volume grid spacing in north direction [m].')
2569 short_type = 'C'
2571 provided_schemes = ['elastic18']
2573 def get_surface_distance(self, args):
2574 _, source_east_shift, source_north_shift, _ = args
2575 sorig = self.source_origin
2576 sloc = Location(
2577 lat=sorig.lat,
2578 lon=sorig.lon,
2579 north_shift=sorig.north_shift + source_north_shift,
2580 east_shift=sorig.east_shift + source_east_shift)
2582 return self.receiver.distance_to(sloc)
2584 def get_distance(self, args):
2585 sdepth, source_east_shift, source_north_shift, _ = args
2586 sorig = self.source_origin
2587 sloc = Location(
2588 lat=sorig.lat,
2589 lon=sorig.lon,
2590 north_shift=sorig.north_shift + source_north_shift,
2591 east_shift=sorig.east_shift + source_east_shift)
2593 return self.receiver.distance_3d_to(sloc)
2595 def get_source_depth(self, args):
2596 return args[0]
2598 def get_receiver_depth(self, args):
2599 return self.receiver.depth
2601 def get_source_depths(self):
2602 return self.coords[0]
2604 def _update(self):
2605 self.mins = num.array([
2606 self.source_depth_min,
2607 self.source_east_shift_min,
2608 self.source_north_shift_min],
2609 dtype=float)
2611 self.maxs = num.array([
2612 self.source_depth_max,
2613 self.source_east_shift_max,
2614 self.source_north_shift_max],
2615 dtype=float)
2617 self.deltas = num.array([
2618 self.source_depth_delta,
2619 self.source_east_shift_delta,
2620 self.source_north_shift_delta],
2621 dtype=float)
2623 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2624 vicinity_eps).astype(int) + 1
2625 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2626 self.deltat = 1.0/self.sample_rate
2627 self.nrecords = num.product(self.ns) * self.ncomponents
2629 self.coords = tuple(
2630 num.linspace(mi, ma, n) for (mi, ma, n) in
2631 zip(self.mins, self.effective_maxs, self.ns)) + \
2632 (num.arange(self.ncomponents),)
2634 def _make_index_functions(self):
2636 amin, bmin, cmin = self.mins
2637 da, db, dc = self.deltas
2638 na, nb, nc = self.ns
2639 ng = self.ncomponents
2641 def index_function(a, b, c, ig):
2642 ia = int(round((a - amin) / da))
2643 ib = int(round((b - bmin) / db))
2644 ic = int(round((c - cmin) / dc))
2645 try:
2646 return num.ravel_multi_index((ia, ib, ic, ig),
2647 (na, nb, nc, ng))
2648 except ValueError:
2649 raise OutOfBounds(values=(a, b, c, ig))
2651 def indices_function(a, b, c, ig):
2652 ia = num.round((a - amin) / da).astype(int)
2653 ib = num.round((b - bmin) / db).astype(int)
2654 ic = num.round((c - cmin) / dc).astype(int)
2656 try:
2657 return num.ravel_multi_index((ia, ib, ic, ig),
2658 (na, nb, nc, ng))
2659 except ValueError:
2660 raise OutOfBounds()
2662 def vicinity_function(a, b, c, ig):
2663 ias = indi12((a - amin) / da, na)
2664 ibs = indi12((b - bmin) / db, nb)
2665 ics = indi12((c - cmin) / dc, nc)
2667 if not (0 <= ig < ng):
2668 raise OutOfBounds()
2670 indis = []
2671 weights = []
2672 for ia, va in ias:
2673 iia = ia*nb*nc*ng
2674 for ib, vb in ibs:
2675 iib = ib*nc*ng
2676 for ic, vc in ics:
2677 indis.append(iia + iib + ic*ng + ig)
2678 weights.append(va*vb*vc)
2680 return num.array(indis), num.array(weights)
2682 def vicinities_function(a, b, c, ig):
2684 xa = (a-amin) / da
2685 xb = (b-bmin) / db
2686 xc = (c-cmin) / dc
2688 xa_fl = num.floor(xa)
2689 xa_ce = num.ceil(xa)
2690 xb_fl = num.floor(xb)
2691 xb_ce = num.ceil(xb)
2692 xc_fl = num.floor(xc)
2693 xc_ce = num.ceil(xc)
2694 va_fl = 1.0 - (xa - xa_fl)
2695 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2696 vb_fl = 1.0 - (xb - xb_fl)
2697 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2698 vc_fl = 1.0 - (xc - xc_fl)
2699 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2701 ia_fl = xa_fl.astype(int)
2702 ia_ce = xa_ce.astype(int)
2703 ib_fl = xb_fl.astype(int)
2704 ib_ce = xb_ce.astype(int)
2705 ic_fl = xc_fl.astype(int)
2706 ic_ce = xc_ce.astype(int)
2708 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2709 raise OutOfBounds()
2711 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2712 raise OutOfBounds()
2714 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2715 raise OutOfBounds()
2717 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2718 raise OutOfBounds()
2720 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2721 raise OutOfBounds()
2723 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2724 raise OutOfBounds()
2726 irecords = num.empty(a.size*8, dtype=int)
2727 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2728 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2729 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2730 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2731 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2732 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2733 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2734 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2736 weights = num.empty(a.size*8, dtype=float)
2737 weights[0::8] = va_fl * vb_fl * vc_fl
2738 weights[1::8] = va_ce * vb_fl * vc_fl
2739 weights[2::8] = va_fl * vb_ce * vc_fl
2740 weights[3::8] = va_ce * vb_ce * vc_fl
2741 weights[4::8] = va_fl * vb_fl * vc_ce
2742 weights[5::8] = va_ce * vb_fl * vc_ce
2743 weights[6::8] = va_fl * vb_ce * vc_ce
2744 weights[7::8] = va_ce * vb_ce * vc_ce
2746 return irecords, weights
2748 self._index_function = index_function
2749 self._indices_function = indices_function
2750 self._vicinity_function = vicinity_function
2751 self._vicinities_function = vicinities_function
2753 def make_indexing_args(self, source, receiver, icomponents):
2754 nc = icomponents.size
2756 dists = source.distances_to(self.source_origin)
2757 azis, _ = source.azibazis_to(self.source_origin)
2759 source_north_shifts = - num.cos(d2r*azis) * dists
2760 source_east_shifts = - num.sin(d2r*azis) * dists
2761 source_depths = source.depths - self.source_origin.depth
2763 n = dists.size
2765 return (num.tile(source_depths, nc//n),
2766 num.tile(source_east_shifts, nc//n),
2767 num.tile(source_north_shifts, nc//n),
2768 icomponents)
2770 def make_indexing_args1(self, source, receiver):
2771 dist = source.distance_to(self.source_origin)
2772 azi, _ = source.azibazi_to(self.source_origin)
2774 source_north_shift = - num.cos(d2r*azi) * dist
2775 source_east_shift = - num.sin(d2r*azi) * dist
2776 source_depth = source.depth - self.source_origin.depth
2778 return (source_depth,
2779 source_east_shift,
2780 source_north_shift)
2782 @property
2783 def short_extent(self):
2784 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2785 self.source_depth_min/km,
2786 self.source_depth_max/km,
2787 self.source_depth_delta/km,
2788 self.source_east_shift_min/km,
2789 self.source_east_shift_max/km,
2790 self.source_east_shift_delta/km,
2791 self.source_north_shift_min/km,
2792 self.source_north_shift_max/km,
2793 self.source_north_shift_delta/km)
2796class Weighting(Object):
2797 factor = Float.T(default=1.0)
2800class Taper(Object):
2801 tmin = Timing.T()
2802 tmax = Timing.T()
2803 tfade = Float.T(default=0.0)
2804 shape = StringChoice.T(
2805 choices=['cos', 'linear'],
2806 default='cos',
2807 optional=True)
2810class SimplePattern(SObject):
2812 _pool = {}
2814 def __init__(self, pattern):
2815 self._pattern = pattern
2816 SObject.__init__(self)
2818 def __str__(self):
2819 return self._pattern
2821 @property
2822 def regex(self):
2823 pool = SimplePattern._pool
2824 if self.pattern not in pool:
2825 rpat = '|'.join(fnmatch.translate(x) for
2826 x in self.pattern.split('|'))
2827 pool[self.pattern] = re.compile(rpat, re.I)
2829 return pool[self.pattern]
2831 def match(self, s):
2832 return self.regex.match(s)
2835class WaveformType(StringChoice):
2836 choices = ['dis', 'vel', 'acc',
2837 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2838 'envelope_dis', 'envelope_vel', 'envelope_acc']
2841class ChannelSelection(Object):
2842 pattern = SimplePattern.T(optional=True)
2843 min_sample_rate = Float.T(optional=True)
2844 max_sample_rate = Float.T(optional=True)
2847class StationSelection(Object):
2848 includes = SimplePattern.T()
2849 excludes = SimplePattern.T()
2850 distance_min = Float.T(optional=True)
2851 distance_max = Float.T(optional=True)
2852 azimuth_min = Float.T(optional=True)
2853 azimuth_max = Float.T(optional=True)
2856class WaveformSelection(Object):
2857 channel_selection = ChannelSelection.T(optional=True)
2858 station_selection = StationSelection.T(optional=True)
2859 taper = Taper.T()
2860 # filter = FrequencyResponse.T()
2861 waveform_type = WaveformType.T(default='dis')
2862 weighting = Weighting.T(optional=True)
2863 sample_rate = Float.T(optional=True)
2864 gf_store_id = StringID.T(optional=True)
2867def indi12(x, n):
2868 '''
2869 Get linear interpolation index and weight.
2870 '''
2872 r = round(x)
2873 if abs(r - x) < vicinity_eps:
2874 i = int(r)
2875 if not (0 <= i < n):
2876 raise OutOfBounds()
2878 return ((int(r), 1.),)
2879 else:
2880 f = math.floor(x)
2881 i = int(f)
2882 if not (0 <= i < n-1):
2883 raise OutOfBounds()
2885 v = x-f
2886 return ((i, 1.-v), (i + 1, v))
2889def float_or_none(s):
2890 units = {
2891 'k': 1e3,
2892 'M': 1e6,
2893 }
2895 factor = 1.0
2896 if s and s[-1] in units:
2897 factor = units[s[-1]]
2898 s = s[:-1]
2899 if not s:
2900 raise ValueError("unit without a number: '%s'" % s)
2902 if s:
2903 return float(s) * factor
2904 else:
2905 return None
2908class GridSpecError(Exception):
2909 def __init__(self, s):
2910 Exception.__init__(self, 'invalid grid specification: %s' % s)
2913def parse_grid_spec(spec):
2914 try:
2915 result = []
2916 for dspec in spec.split(','):
2917 t = dspec.split('@')
2918 num = start = stop = step = None
2919 if len(t) == 2:
2920 num = int(t[1])
2921 if num <= 0:
2922 raise GridSpecError(spec)
2924 elif len(t) > 2:
2925 raise GridSpecError(spec)
2927 s = t[0]
2928 v = [float_or_none(x) for x in s.split(':')]
2929 if len(v) == 1:
2930 start = stop = v[0]
2931 if len(v) >= 2:
2932 start, stop = v[0:2]
2933 if len(v) == 3:
2934 step = v[2]
2936 if len(v) > 3 or (len(v) > 2 and num is not None):
2937 raise GridSpecError(spec)
2939 if step == 0.0:
2940 raise GridSpecError(spec)
2942 result.append((start, stop, step, num))
2944 except ValueError:
2945 raise GridSpecError(spec)
2947 return result
2950def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2951 swap = step is not None and step < 0.
2952 if start is None:
2953 start = ma if swap else mi
2954 if stop is None:
2955 stop = mi if swap else ma
2956 if step is None:
2957 step = -inc if ma < mi else inc
2958 if num is None:
2959 if (step < 0) != (stop-start < 0):
2960 raise GridSpecError()
2962 num = int(round((stop-start)/step))+1
2963 stop2 = start + (num-1)*step
2964 if abs(stop-stop2) > eps:
2965 num = int(math.floor((stop-start)/step))+1
2966 stop = start + (num-1)*step
2967 else:
2968 stop = stop2
2970 if start == stop:
2971 num = 1
2973 return start, stop, num
2976def nditer_outer(x):
2977 return num.nditer(
2978 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
2981def nodes(xs):
2982 ns = [x.size for x in xs]
2983 nnodes = num.prod(ns)
2984 ndim = len(xs)
2985 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
2986 for idim in range(ndim-1, -1, -1):
2987 x = xs[idim]
2988 nrepeat = num.prod(ns[idim+1:], dtype=int)
2989 ntile = num.prod(ns[:idim], dtype=int)
2990 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
2992 return nodes
2995def filledi(x, n):
2996 a = num.empty(n, dtype=int)
2997 a.fill(x)
2998 return a
3001config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3003discretized_source_classes = [
3004 DiscretizedExplosionSource,
3005 DiscretizedSFSource,
3006 DiscretizedMTSource,
3007 DiscretizedPorePressureSource]
3010__all__ = '''
3011Earthmodel1D
3012StringID
3013ScopeType
3014WaveformType
3015QuantityType
3016NearfieldTermsType
3017Reference
3018Region
3019CircularRegion
3020RectangularRegion
3021PhaseSelect
3022InvalidTimingSpecification
3023Timing
3024TPDef
3025OutOfBounds
3026Location
3027Receiver
3028'''.split() + [
3029 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3030ComponentScheme
3031component_scheme_to_description
3032component_schemes
3033Config
3034GridSpecError
3035Weighting
3036Taper
3037SimplePattern
3038WaveformType
3039ChannelSelection
3040StationSelection
3041WaveformSelection
3042nditer_outer
3043dump
3044load
3045discretized_source_classes
3046config_type_classes
3047UnavailableScheme
3048InterpolationMethod
3049SeismosizerTrace
3050SeismosizerResult
3051Result
3052StaticResult
3053'''.split()