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 TimingAttributeError(Exception):
405 pass
408class Timing(SObject):
409 '''
410 Definition of a time instant relative to one or more named phase arrivals.
412 Instances of this class can be used e.g. in cutting and tapering
413 operations. They can hold an absolute time or an offset to a named phase
414 arrival or group of such arrivals.
416 Timings can be instantiated from a simple string defintion i.e. with
417 ``Timing(str)`` where ``str`` is something like
418 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``,
419 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase
420 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'``
421 is appended, it is interpreted as percent. If the an ``'S'`` is appended
422 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km].
424 Phase definitions can be specified in either of the following ways:
426 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table
427 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake
428 (see :py:class:`pyrocko.cake.PhaseDef`)
429 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance /
430 velocity [km/s]
431 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s]
433 **Examples:**
435 * ``'100'`` : absolute time; 100 s
436 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to
437 stored travel time table named ``'P'``
438 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to
439 stored travel time table named ``'P'``
440 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to
441 stored travel time table named ``'P'``
442 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is
443 undefined for a given geometry
444 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A
445 and B is chosen, if both phases are defined for a given geometry
446 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen
447 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of
448 arrivals A, B, and C
449 '''
451 def __init__(self, s=None, **kwargs):
453 if s is not None:
454 offset_is = None
455 s = re.sub(r'\s+', '', s)
456 try:
457 offset = float(s.rstrip('S'))
459 if s.endswith('S'):
460 offset_is = 'slowness'
462 phase_defs = []
463 select = ''
465 except ValueError:
466 matched = False
467 m = timing_regex.match(s)
468 if m:
469 if m.group(3):
470 phase_defs = m.group(3).split('|')
471 elif m.group(19):
472 phase_defs = [m.group(19)]
473 else:
474 phase_defs = []
476 select = m.group(2) or ''
478 offset = 0.0
479 soff = m.group(27)
480 if soff:
481 offset = float(soff.rstrip('S%'))
482 if soff.endswith('S'):
483 offset_is = 'slowness'
484 elif soff.endswith('%'):
485 offset_is = 'percent'
487 matched = True
489 else:
490 m = timing_regex_old.match(s)
491 if m:
492 if m.group(3):
493 phase_defs = m.group(3).split('|')
494 elif m.group(5):
495 phase_defs = [m.group(5)]
496 else:
497 phase_defs = []
499 select = m.group(2) or ''
501 offset = 0.0
502 if m.group(6):
503 offset = float(m.group(6))
505 phase_defs = [
506 'stored:' + phase_def for phase_def in phase_defs]
508 matched = True
510 if not matched:
511 raise InvalidTimingSpecification(s)
513 kwargs = dict(
514 phase_defs=phase_defs,
515 select=select,
516 offset=offset,
517 offset_is=offset_is)
519 SObject.__init__(self, **kwargs)
521 def __str__(self):
522 s = []
523 if self.phase_defs:
524 sphases = '|'.join(self.phase_defs)
525 # if len(self.phase_defs) > 1 or self.select:
526 sphases = '{%s}' % sphases
528 if self.select:
529 sphases = self.select + sphases
531 s.append(sphases)
533 if self.offset != 0.0 or not self.phase_defs:
534 s.append('%+g' % self.offset)
535 if self.offset_is == 'slowness':
536 s.append('S')
537 elif self.offset_is == 'percent':
538 s.append('%')
540 return ''.join(s)
542 def evaluate(self, get_phase, args, attributes=None):
543 try:
544 if self.offset_is == 'slowness' and self.offset != 0.0:
545 phase_offset = get_phase(
546 'vel_surface:%g' % (1.0/self.offset))
547 offset = phase_offset(args)
548 else:
549 offset = self.offset
551 if self.phase_defs:
552 extra_args = ()
553 if attributes:
554 extra_args = (attributes,)
556 phases = [
557 get_phase(phase_def, *extra_args)
558 for phase_def in self.phase_defs]
560 results = [phase(args) for phase in phases]
561 results = [
562 result if isinstance(result, tuple) else (result,)
563 for result in results]
565 results = [
566 result for result in results
567 if result[0] is not None]
569 if self.select == 'first':
570 results.sort(key=lambda result: result[0])
571 elif self.select == 'last':
572 results.sort(key=lambda result: -result[0])
574 if not results:
575 if attributes is not None:
576 return (None,) * (len(attributes) + 1)
577 else:
578 return None
580 else:
581 t, values = results[0][0], results[0][1:]
582 if self.offset_is == 'percent':
583 t = t*(1.+offset/100.)
584 else:
585 t = t + offset
587 if attributes is not None:
588 return (t, ) + values
589 else:
590 return t
592 else:
593 if attributes is not None:
594 raise TimingAttributeError(
595 'Cannot get phase attributes from offset-only '
596 'definition.')
597 return offset
599 except spit.OutOfBounds:
600 raise OutOfBounds(args)
602 phase_defs = List.T(String.T())
603 offset = Float.T(default=0.0)
604 offset_is = String.T(optional=True)
605 select = PhaseSelect.T(
606 default='',
607 help=("Can be either ``'%s'``, ``'%s'``, or ``'%s'``. " %
608 tuple(PhaseSelect.choices)))
611def mkdefs(s):
612 defs = []
613 for x in s.split(','):
614 try:
615 defs.append(float(x))
616 except ValueError:
617 if x.startswith('!'):
618 defs.extend(cake.PhaseDef.classic(x[1:]))
619 else:
620 defs.append(cake.PhaseDef(x))
622 return defs
625class TPDef(Object):
626 '''
627 Maps an arrival phase identifier to an arrival phase definition.
628 '''
630 id = StringID.T(
631 help='name used to identify the phase')
632 definition = String.T(
633 help='definition of the phase in either cake syntax as defined in '
634 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an '
635 '``!``, as a *classic phase name*, or, if it is a simple '
636 'number, as a constant horizontal velocity.')
638 @property
639 def phases(self):
640 return [x for x in mkdefs(self.definition)
641 if isinstance(x, cake.PhaseDef)]
643 @property
644 def horizontal_velocities(self):
645 return [x for x in mkdefs(self.definition) if isinstance(x, float)]
648class OutOfBounds(Exception):
649 def __init__(self, values=None, reason=None):
650 Exception.__init__(self)
651 self.values = values
652 self.reason = reason
653 self.context = None
655 def __str__(self):
656 scontext = ''
657 if self.context:
658 scontext = '\n%s' % str(self.context)
660 if self.reason:
661 scontext += '\n%s' % self.reason
663 if self.values:
664 return 'out of bounds: (%s)%s' % (
665 ', '.join('%g' % x for x in self.values), scontext)
666 else:
667 return 'out of bounds%s' % scontext
670class MultiLocation(Object):
672 lats = Array.T(
673 optional=True, shape=(None,), dtype=float,
674 help='Latitudes of targets.')
675 lons = Array.T(
676 optional=True, shape=(None,), dtype=float,
677 help='Longitude of targets.')
678 north_shifts = Array.T(
679 optional=True, shape=(None,), dtype=float,
680 help='North shifts of targets.')
681 east_shifts = Array.T(
682 optional=True, shape=(None,), dtype=float,
683 help='East shifts of targets.')
684 elevation = Array.T(
685 optional=True, shape=(None,), dtype=float,
686 help='Elevations of targets.')
688 def __init__(self, *args, **kwargs):
689 self._coords5 = None
690 Object.__init__(self, *args, **kwargs)
692 @property
693 def coords5(self):
694 if self._coords5 is not None:
695 return self._coords5
696 props = [self.lats, self.lons, self.north_shifts, self.east_shifts,
697 self.elevation]
698 sizes = [p.size for p in props if p is not None]
699 if not sizes:
700 raise AttributeError('Empty StaticTarget')
701 if num.unique(sizes).size != 1:
702 raise AttributeError('Inconsistent coordinate sizes.')
704 self._coords5 = num.zeros((sizes[0], 5))
705 for idx, p in enumerate(props):
706 if p is not None:
707 self._coords5[:, idx] = p.flatten()
709 return self._coords5
711 @property
712 def ncoords(self):
713 if self.coords5.shape[0] is None:
714 return 0
715 return int(self.coords5.shape[0])
717 def get_latlon(self):
718 '''
719 Get all coordinates as lat/lon.
721 :returns: Coordinates in Latitude, Longitude
722 :rtype: :class:`numpy.ndarray`, (N, 2)
723 '''
724 latlons = num.empty((self.ncoords, 2))
725 for i in range(self.ncoords):
726 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4])
727 return latlons
729 def get_corner_coords(self):
730 '''
731 Returns the corner coordinates of the multi-location object.
733 :returns: In lat/lon: lower left, upper left, upper right, lower right
734 :rtype: tuple
735 '''
736 latlon = self.get_latlon()
737 ll = (latlon[:, 0].min(), latlon[:, 1].min())
738 ur = (latlon[:, 0].max(), latlon[:, 1].max())
739 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1]))
742class Receiver(Location):
743 codes = Tuple.T(
744 3, String.T(),
745 optional=True,
746 help='network, station, and location codes')
748 def pyrocko_station(self):
749 from pyrocko import model
751 lat, lon = self.effective_latlon
753 return model.Station(
754 network=self.codes[0],
755 station=self.codes[1],
756 location=self.codes[2],
757 lat=lat,
758 lon=lon,
759 depth=self.depth)
762def g(x, d):
763 if x is None:
764 return d
765 else:
766 return x
769class UnavailableScheme(Exception):
770 pass
773class InvalidNComponents(Exception):
774 pass
777class DiscretizedSource(Object):
778 '''
779 Base class for discretized sources.
781 To compute synthetic seismograms, the parameterized source models (see of
782 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first
783 discretized into a number of point sources. These spacio-temporal point
784 source distributions are represented by subclasses of the
785 :py:class:`DiscretizedSource`. For elastodynamic problems there is the
786 :py:class:`DiscretizedMTSource` for moment tensor point source
787 distributions and the :py:class:`DiscretizedExplosionSource` for pure
788 explosion/implosion type source distributions. The geometry calculations
789 are implemented in the base class. How Green's function components have to
790 be weighted and sumed is defined in the derived classes.
792 Like in the :py:class:`Location` class, the positions of the point sources
793 contained in the discretized source are defined by a common reference point
794 (:py:attr:`lat`, :py:attr:`lon`) and cartesian offsets to this
795 (:py:attr:`north_shifts`, :py:attr:`east_shifts`, :py:attr:`depths`).
796 Alternatively latitude and longitude of each contained point source can be
797 specified directly (:py:attr:`lats`, :py:attr:`lons`).
798 '''
800 times = Array.T(shape=(None,), dtype=float)
801 lats = Array.T(shape=(None,), dtype=float, optional=True)
802 lons = Array.T(shape=(None,), dtype=float, optional=True)
803 lat = Float.T(optional=True)
804 lon = Float.T(optional=True)
805 north_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
806 east_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)
807 depths = Array.T(shape=(None,), dtype=num.float64)
808 dl = Float.T(optional=True)
809 dw = Float.T(optional=True)
810 nl = Float.T(optional=True)
811 nw = Float.T(optional=True)
813 @classmethod
814 def check_scheme(cls, scheme):
815 '''
816 Check if given GF component scheme is supported by the class.
818 Raises :py:class:`UnavailableScheme` if the given scheme is not
819 supported by this discretized source class.
820 '''
822 if scheme not in cls.provided_schemes:
823 raise UnavailableScheme(
824 'source type "%s" does not support GF component scheme "%s"' %
825 (cls.__name__, scheme))
827 def __init__(self, **kwargs):
828 Object.__init__(self, **kwargs)
829 self._latlons = None
831 def __setattr__(self, name, value):
832 if name in ('lat', 'lon', 'north_shifts', 'east_shifts',
833 'lats', 'lons'):
834 self.__dict__['_latlons'] = None
836 Object.__setattr__(self, name, value)
838 def get_source_terms(self, scheme):
839 raise NotImplementedError()
841 def make_weights(self, receiver, scheme):
842 raise NotImplementedError()
844 @property
845 def effective_latlons(self):
846 '''
847 Property holding the offest-corrected lats and lons of all points.
848 '''
850 if self._latlons is None:
851 if self.lats is not None and self.lons is not None:
852 if (self.north_shifts is not None and
853 self.east_shifts is not None):
854 self._latlons = orthodrome.ne_to_latlon(
855 self.lats, self.lons,
856 self.north_shifts, self.east_shifts)
857 else:
858 self._latlons = self.lats, self.lons
859 else:
860 lat = g(self.lat, 0.0)
861 lon = g(self.lon, 0.0)
862 self._latlons = orthodrome.ne_to_latlon(
863 lat, lon, self.north_shifts, self.east_shifts)
865 return self._latlons
867 @property
868 def effective_north_shifts(self):
869 if self.north_shifts is not None:
870 return self.north_shifts
871 else:
872 return num.zeros(self.times.size)
874 @property
875 def effective_east_shifts(self):
876 if self.east_shifts is not None:
877 return self.east_shifts
878 else:
879 return num.zeros(self.times.size)
881 def same_origin(self, receiver):
882 '''
883 Check if receiver has same reference point.
884 '''
886 return (g(self.lat, 0.0) == receiver.lat and
887 g(self.lon, 0.0) == receiver.lon and
888 self.lats is None and self.lons is None)
890 def azibazis_to(self, receiver):
891 '''
892 Compute azimuths and backaziumuths to/from receiver, for all contained
893 points.
894 '''
896 if self.same_origin(receiver):
897 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts,
898 receiver.north_shift - self.north_shifts)
899 bazis = azis + 180.
900 else:
901 slats, slons = self.effective_latlons
902 rlat, rlon = receiver.effective_latlon
903 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon)
904 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons)
906 return azis, bazis
908 def distances_to(self, receiver):
909 '''
910 Compute distances to receiver for all contained points.
911 '''
912 if self.same_origin(receiver):
913 return num.sqrt((self.north_shifts - receiver.north_shift)**2 +
914 (self.east_shifts - receiver.east_shift)**2)
916 else:
917 slats, slons = self.effective_latlons
918 rlat, rlon = receiver.effective_latlon
919 return orthodrome.distance_accurate50m_numpy(slats, slons,
920 rlat, rlon)
922 def element_coords(self, i):
923 if self.lats is not None and self.lons is not None:
924 lat = float(self.lats[i])
925 lon = float(self.lons[i])
926 else:
927 lat = self.lat
928 lon = self.lon
930 if self.north_shifts is not None and self.east_shifts is not None:
931 north_shift = float(self.north_shifts[i])
932 east_shift = float(self.east_shifts[i])
934 else:
935 north_shift = east_shift = 0.0
937 return lat, lon, north_shift, east_shift
939 def coords5(self):
940 xs = num.zeros((self.nelements, 5))
942 if self.lats is not None and self.lons is not None:
943 xs[:, 0] = self.lats
944 xs[:, 1] = self.lons
945 else:
946 xs[:, 0] = g(self.lat, 0.0)
947 xs[:, 1] = g(self.lon, 0.0)
949 if self.north_shifts is not None and self.east_shifts is not None:
950 xs[:, 2] = self.north_shifts
951 xs[:, 3] = self.east_shifts
953 xs[:, 4] = self.depths
955 return xs
957 @property
958 def nelements(self):
959 return self.times.size
961 @classmethod
962 def combine(cls, sources, **kwargs):
963 '''
964 Combine several discretized source models.
966 Concatenenates all point sources in the given discretized ``sources``.
967 Care must be taken when using this function that the external amplitude
968 factors and reference times of the parameterized (undiscretized)
969 sources match or are accounted for.
970 '''
972 first = sources[0]
974 if not all(type(s) == type(first) for s in sources):
975 raise Exception('DiscretizedSource.combine must be called with '
976 'sources of same type.')
978 latlons = []
979 for s in sources:
980 latlons.append(s.effective_latlons)
982 lats, lons = num.hstack(latlons)
984 if all((s.lats is None and s.lons is None) for s in sources):
985 rlats = num.array([s.lat for s in sources], dtype=float)
986 rlons = num.array([s.lon for s in sources], dtype=float)
987 same_ref = num.all(
988 rlats == rlats[0]) and num.all(rlons == rlons[0])
989 else:
990 same_ref = False
992 cat = num.concatenate
993 times = cat([s.times for s in sources])
994 print(times)
995 depths = cat([s.depths for s in sources])
997 if same_ref:
998 lat = first.lat
999 lon = first.lon
1000 north_shifts = cat([s.effective_north_shifts for s in sources])
1001 east_shifts = cat([s.effective_east_shifts for s in sources])
1002 lats = None
1003 lons = None
1004 else:
1005 lat = None
1006 lon = None
1007 north_shifts = None
1008 east_shifts = None
1010 return cls(
1011 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
1012 north_shifts=north_shifts, east_shifts=east_shifts,
1013 depths=depths, **kwargs)
1015 def centroid_position(self):
1016 moments = self.moments()
1017 norm = num.sum(moments)
1018 if norm != 0.0:
1019 w = moments / num.sum(moments)
1020 else:
1021 w = num.ones(moments.size)
1023 if self.lats is not None and self.lons is not None:
1024 lats, lons = self.effective_latlons
1025 rlat, rlon = num.mean(lats), num.mean(lons)
1026 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
1027 else:
1028 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
1029 n, e = self.north_shifts, self.east_shifts
1031 cn = num.sum(n*w)
1032 ce = num.sum(e*w)
1033 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1035 if self.lats is not None and self.lons is not None:
1036 lat = clat
1037 lon = clon
1038 north_shift = 0.
1039 east_shift = 0.
1040 else:
1041 lat = g(self.lat, 0.0)
1042 lon = g(self.lon, 0.0)
1043 north_shift = cn
1044 east_shift = ce
1046 depth = num.sum(self.depths*w)
1047 time = num.sum(self.times*w)
1048 return tuple(float(x) for x in
1049 (time, lat, lon, north_shift, east_shift, depth))
1052class DiscretizedExplosionSource(DiscretizedSource):
1053 m0s = Array.T(shape=(None,), dtype=float)
1055 provided_schemes = (
1056 'elastic2',
1057 'elastic8',
1058 'elastic10',
1059 )
1061 def get_source_terms(self, scheme):
1062 self.check_scheme(scheme)
1064 if scheme == 'elastic2':
1065 return self.m0s[:, num.newaxis].copy()
1067 elif scheme in ('elastic8', 'elastic10'):
1068 m6s = num.zeros((self.m0s.size, 6))
1069 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1070 return m6s
1071 else:
1072 assert False
1074 def make_weights(self, receiver, scheme):
1075 self.check_scheme(scheme)
1077 azis, bazis = self.azibazis_to(receiver)
1079 sb = num.sin(bazis*d2r-num.pi)
1080 cb = num.cos(bazis*d2r-num.pi)
1082 m0s = self.m0s
1083 n = azis.size
1085 cat = num.concatenate
1086 rep = num.repeat
1088 if scheme == 'elastic2':
1089 w_n = cb*m0s
1090 g_n = filledi(0, n)
1091 w_e = sb*m0s
1092 g_e = filledi(0, n)
1093 w_d = m0s
1094 g_d = filledi(1, n)
1096 elif scheme == 'elastic8':
1097 w_n = cat((cb*m0s, cb*m0s))
1098 g_n = rep((0, 2), n)
1099 w_e = cat((sb*m0s, sb*m0s))
1100 g_e = rep((0, 2), n)
1101 w_d = cat((m0s, m0s))
1102 g_d = rep((5, 7), n)
1104 elif scheme == 'elastic10':
1105 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1106 g_n = rep((0, 2, 8), n)
1107 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1108 g_e = rep((0, 2, 8), n)
1109 w_d = cat((m0s, m0s, m0s))
1110 g_d = rep((5, 7, 9), n)
1112 else:
1113 assert False
1115 return (
1116 ('displacement.n', w_n, g_n),
1117 ('displacement.e', w_e, g_e),
1118 ('displacement.d', w_d, g_d),
1119 )
1121 def split(self):
1122 from pyrocko.gf.seismosizer import ExplosionSource
1123 sources = []
1124 for i in range(self.nelements):
1125 lat, lon, north_shift, east_shift = self.element_coords(i)
1126 sources.append(ExplosionSource(
1127 time=float(self.times[i]),
1128 lat=lat,
1129 lon=lon,
1130 north_shift=north_shift,
1131 east_shift=east_shift,
1132 depth=float(self.depths[i]),
1133 moment=float(self.m0s[i])))
1135 return sources
1137 def moments(self):
1138 return self.m0s
1140 def centroid(self):
1141 from pyrocko.gf.seismosizer import ExplosionSource
1142 time, lat, lon, north_shift, east_shift, depth = \
1143 self.centroid_position()
1145 return ExplosionSource(
1146 time=time,
1147 lat=lat,
1148 lon=lon,
1149 north_shift=north_shift,
1150 east_shift=east_shift,
1151 depth=depth,
1152 moment=float(num.sum(self.m0s)))
1154 @classmethod
1155 def combine(cls, sources, **kwargs):
1156 '''
1157 Combine several discretized source models.
1159 Concatenenates all point sources in the given discretized ``sources``.
1160 Care must be taken when using this function that the external amplitude
1161 factors and reference times of the parameterized (undiscretized)
1162 sources match or are accounted for.
1163 '''
1165 if 'm0s' not in kwargs:
1166 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1168 return super(DiscretizedExplosionSource, cls).combine(sources,
1169 **kwargs)
1172class DiscretizedSFSource(DiscretizedSource):
1173 forces = Array.T(shape=(None, 3), dtype=float)
1175 provided_schemes = (
1176 'elastic5',
1177 )
1179 def get_source_terms(self, scheme):
1180 self.check_scheme(scheme)
1182 return self.forces
1184 def make_weights(self, receiver, scheme):
1185 self.check_scheme(scheme)
1187 azis, bazis = self.azibazis_to(receiver)
1189 sa = num.sin(azis*d2r)
1190 ca = num.cos(azis*d2r)
1191 sb = num.sin(bazis*d2r-num.pi)
1192 cb = num.cos(bazis*d2r-num.pi)
1194 forces = self.forces
1195 fn = forces[:, 0]
1196 fe = forces[:, 1]
1197 fd = forces[:, 2]
1199 f0 = fd
1200 f1 = ca * fn + sa * fe
1201 f2 = ca * fe - sa * fn
1203 n = azis.size
1205 if scheme == 'elastic5':
1206 ioff = 0
1208 cat = num.concatenate
1209 rep = num.repeat
1211 w_n = cat((cb*f0, cb*f1, -sb*f2))
1212 g_n = ioff + rep((0, 1, 2), n)
1213 w_e = cat((sb*f0, sb*f1, cb*f2))
1214 g_e = ioff + rep((0, 1, 2), n)
1215 w_d = cat((f0, f1))
1216 g_d = ioff + rep((3, 4), n)
1218 return (
1219 ('displacement.n', w_n, g_n),
1220 ('displacement.e', w_e, g_e),
1221 ('displacement.d', w_d, g_d),
1222 )
1224 @classmethod
1225 def combine(cls, sources, **kwargs):
1226 '''
1227 Combine several discretized source models.
1229 Concatenenates all point sources in the given discretized ``sources``.
1230 Care must be taken when using this function that the external amplitude
1231 factors and reference times of the parameterized (undiscretized)
1232 sources match or are accounted for.
1233 '''
1235 if 'forces' not in kwargs:
1236 kwargs['forces'] = num.vstack([s.forces for s in sources])
1238 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1240 def moments(self):
1241 return num.sum(self.forces**2, axis=1)
1243 def centroid(self):
1244 from pyrocko.gf.seismosizer import SFSource
1245 time, lat, lon, north_shift, east_shift, depth = \
1246 self.centroid_position()
1248 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1249 return SFSource(
1250 time=time,
1251 lat=lat,
1252 lon=lon,
1253 north_shift=north_shift,
1254 east_shift=east_shift,
1255 depth=depth,
1256 fn=fn,
1257 fe=fe,
1258 fd=fd)
1261class DiscretizedMTSource(DiscretizedSource):
1262 m6s = Array.T(
1263 shape=(None, 6), dtype=float,
1264 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1266 provided_schemes = (
1267 'elastic8',
1268 'elastic10',
1269 'elastic18',
1270 )
1272 def get_source_terms(self, scheme):
1273 self.check_scheme(scheme)
1274 return self.m6s
1276 def make_weights(self, receiver, scheme):
1277 self.check_scheme(scheme)
1279 m6s = self.m6s
1280 n = m6s.shape[0]
1281 rep = num.repeat
1283 if scheme == 'elastic18':
1284 w_n = m6s.flatten()
1285 g_n = num.tile(num.arange(0, 6), n)
1286 w_e = m6s.flatten()
1287 g_e = num.tile(num.arange(6, 12), n)
1288 w_d = m6s.flatten()
1289 g_d = num.tile(num.arange(12, 18), n)
1291 else:
1292 azis, bazis = self.azibazis_to(receiver)
1294 sa = num.sin(azis*d2r)
1295 ca = num.cos(azis*d2r)
1296 s2a = num.sin(2.*azis*d2r)
1297 c2a = num.cos(2.*azis*d2r)
1298 sb = num.sin(bazis*d2r-num.pi)
1299 cb = num.cos(bazis*d2r-num.pi)
1301 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1302 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1303 f2 = m6s[:, 2]
1304 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1305 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1306 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1308 cat = num.concatenate
1310 if scheme == 'elastic8':
1311 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1312 g_n = rep((0, 1, 2, 3, 4), n)
1313 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1314 g_e = rep((0, 1, 2, 3, 4), n)
1315 w_d = cat((f0, f1, f2))
1316 g_d = rep((5, 6, 7), n)
1318 elif scheme == 'elastic10':
1319 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1320 g_n = rep((0, 1, 2, 8, 3, 4), n)
1321 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1322 g_e = rep((0, 1, 2, 8, 3, 4), n)
1323 w_d = cat((f0, f1, f2, f5))
1324 g_d = rep((5, 6, 7, 9), n)
1326 else:
1327 assert False
1329 return (
1330 ('displacement.n', w_n, g_n),
1331 ('displacement.e', w_e, g_e),
1332 ('displacement.d', w_d, g_d),
1333 )
1335 def split(self):
1336 from pyrocko.gf.seismosizer import MTSource
1337 sources = []
1338 for i in range(self.nelements):
1339 lat, lon, north_shift, east_shift = self.element_coords(i)
1340 sources.append(MTSource(
1341 time=float(self.times[i]),
1342 lat=lat,
1343 lon=lon,
1344 north_shift=north_shift,
1345 east_shift=east_shift,
1346 depth=float(self.depths[i]),
1347 m6=self.m6s[i]))
1349 return sources
1351 def moments(self):
1352 moments = num.array(
1353 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))
1354 for m6 in self.m6s])
1355 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)
1357 def get_moment_rate(self, deltat=None):
1358 moments = self.moments()
1359 times = self.times
1360 times -= times.min()
1362 t_max = times.max()
1363 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat
1364 mom_times[mom_times > t_max] = t_max
1366 # Right open histrogram bins
1367 mom, _ = num.histogram(
1368 times,
1369 bins=mom_times,
1370 weights=moments)
1372 deltat = num.diff(mom_times)
1373 mom_rate = mom / deltat
1375 return mom_rate, mom_times[1:]
1377 def centroid(self):
1378 from pyrocko.gf.seismosizer import MTSource
1379 time, lat, lon, north_shift, east_shift, depth = \
1380 self.centroid_position()
1382 return MTSource(
1383 time=time,
1384 lat=lat,
1385 lon=lon,
1386 north_shift=north_shift,
1387 east_shift=east_shift,
1388 depth=depth,
1389 m6=num.sum(self.m6s, axis=0))
1391 @classmethod
1392 def combine(cls, sources, **kwargs):
1393 '''
1394 Combine several discretized source models.
1396 Concatenenates all point sources in the given discretized ``sources``.
1397 Care must be taken when using this function that the external amplitude
1398 factors and reference times of the parameterized (undiscretized)
1399 sources match or are accounted for.
1400 '''
1402 if 'm6s' not in kwargs:
1403 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1405 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1408class DiscretizedPorePressureSource(DiscretizedSource):
1409 pp = Array.T(shape=(None,), dtype=float)
1411 provided_schemes = (
1412 'poroelastic10',
1413 )
1415 def get_source_terms(self, scheme):
1416 self.check_scheme(scheme)
1417 return self.pp[:, num.newaxis].copy()
1419 def make_weights(self, receiver, scheme):
1420 self.check_scheme(scheme)
1422 azis, bazis = self.azibazis_to(receiver)
1424 sb = num.sin(bazis*d2r-num.pi)
1425 cb = num.cos(bazis*d2r-num.pi)
1427 pp = self.pp
1428 n = bazis.size
1430 w_un = cb*pp
1431 g_un = filledi(1, n)
1432 w_ue = sb*pp
1433 g_ue = filledi(1, n)
1434 w_ud = pp
1435 g_ud = filledi(0, n)
1437 w_tn = cb*pp
1438 g_tn = filledi(6, n)
1439 w_te = sb*pp
1440 g_te = filledi(6, n)
1442 w_pp = pp
1443 g_pp = filledi(7, n)
1445 w_dvn = cb*pp
1446 g_dvn = filledi(9, n)
1447 w_dve = sb*pp
1448 g_dve = filledi(9, n)
1449 w_dvd = pp
1450 g_dvd = filledi(8, n)
1452 return (
1453 ('displacement.n', w_un, g_un),
1454 ('displacement.e', w_ue, g_ue),
1455 ('displacement.d', w_ud, g_ud),
1456 ('vertical_tilt.n', w_tn, g_tn),
1457 ('vertical_tilt.e', w_te, g_te),
1458 ('pore_pressure', w_pp, g_pp),
1459 ('darcy_velocity.n', w_dvn, g_dvn),
1460 ('darcy_velocity.e', w_dve, g_dve),
1461 ('darcy_velocity.d', w_dvd, g_dvd),
1462 )
1464 def moments(self):
1465 return self.pp
1467 def centroid(self):
1468 from pyrocko.gf.seismosizer import PorePressurePointSource
1469 time, lat, lon, north_shift, east_shift, depth = \
1470 self.centroid_position()
1472 return PorePressurePointSource(
1473 time=time,
1474 lat=lat,
1475 lon=lon,
1476 north_shift=north_shift,
1477 east_shift=east_shift,
1478 depth=depth,
1479 pp=float(num.sum(self.pp)))
1481 @classmethod
1482 def combine(cls, sources, **kwargs):
1483 '''
1484 Combine several discretized source models.
1486 Concatenenates all point sources in the given discretized ``sources``.
1487 Care must be taken when using this function that the external amplitude
1488 factors and reference times of the parameterized (undiscretized)
1489 sources match or are accounted for.
1490 '''
1492 if 'pp' not in kwargs:
1493 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1495 return super(DiscretizedPorePressureSource, cls).combine(sources,
1496 **kwargs)
1499class Region(Object):
1500 name = String.T(optional=True)
1503class RectangularRegion(Region):
1504 lat_min = Float.T()
1505 lat_max = Float.T()
1506 lon_min = Float.T()
1507 lon_max = Float.T()
1510class CircularRegion(Region):
1511 lat = Float.T()
1512 lon = Float.T()
1513 radius = Float.T()
1516class Config(Object):
1517 '''
1518 Green's function store meta information.
1520 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1521 configuration types are:
1523 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1524 spherical symmetry, 1D earth model, single receiver depth
1526 * Problem is invariant to horizontal translations and rotations around
1527 vertical axis.
1528 * All receivers must be at the same depth (e.g. at the surface)
1529 * High level index variables: ``(source_depth, receiver_distance,
1530 component)``
1532 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1533 spherical symmetry, 1D earth model, variable receiver depth
1535 * Symmetries like in Type A but has additional index for receiver depth
1536 * High level index variables: ``(source_depth, receiver_distance,
1537 receiver_depth, component)``
1539 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1540 constraints but fixed receiver positions
1542 * Cartesian source volume around a reference point
1543 * High level index variables: ``(ireceiver, source_depth,
1544 source_east_shift, source_north_shift, component)``
1545 '''
1547 id = StringID.T(
1548 help='Name of the store. May consist of upper and lower-case letters, '
1549 'digits, dots and underscores. The name must start with a '
1550 'letter.')
1552 derived_from_id = StringID.T(
1553 optional=True,
1554 help='Name of the original store, if this store has been derived from '
1555 'another one (e.g. extracted subset).')
1557 version = String.T(
1558 default='1.0',
1559 optional=True,
1560 help='User-defined version string. Use <major>.<minor> format.')
1562 modelling_code_id = StringID.T(
1563 optional=True,
1564 help='Identifier of the backend used to compute the store.')
1566 author = Unicode.T(
1567 optional=True,
1568 help='Comma-separated list of author names.')
1570 author_email = String.T(
1571 optional=True,
1572 help="Author's contact email address.")
1574 created_time = Timestamp.T(
1575 optional=True,
1576 help='Time of creation of the store.')
1578 regions = List.T(
1579 Region.T(),
1580 help='Geographical regions for which the store is representative.')
1582 scope_type = ScopeType.T(
1583 optional=True,
1584 help='Distance range scope of the store '
1585 '(%s).' % fmt_choices(ScopeType))
1587 waveform_type = WaveType.T(
1588 optional=True,
1589 help='Wave type stored (%s).' % fmt_choices(WaveType))
1591 nearfield_terms = NearfieldTermsType.T(
1592 optional=True,
1593 help='Information about the inclusion of near-field terms in the '
1594 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1596 description = String.T(
1597 optional=True,
1598 help='Free form textual description of the GF store.')
1600 references = List.T(
1601 Reference.T(),
1602 help='Reference list to cite the modelling code, earth model or '
1603 'related work.')
1605 earthmodel_1d = Earthmodel1D.T(
1606 optional=True,
1607 help='Layered earth model in ND (named discontinuity) format.')
1609 earthmodel_receiver_1d = Earthmodel1D.T(
1610 optional=True,
1611 help='Receiver-side layered earth model in ND format.')
1613 can_interpolate_source = Bool.T(
1614 optional=True,
1615 help='Hint to indicate if the spatial sampling of the store is dense '
1616 'enough for multi-linear interpolation at the source.')
1618 can_interpolate_receiver = Bool.T(
1619 optional=True,
1620 help='Hint to indicate if the spatial sampling of the store is dense '
1621 'enough for multi-linear interpolation at the receiver.')
1623 frequency_min = Float.T(
1624 optional=True,
1625 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1627 frequency_max = Float.T(
1628 optional=True,
1629 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1631 sample_rate = Float.T(
1632 optional=True,
1633 help='Sample rate of the GF store [Hz].')
1635 factor = Float.T(
1636 default=1.0,
1637 help='Gain value, factored out of the stored GF samples. '
1638 '(may not work properly, keep at 1.0).',
1639 optional=True)
1641 component_scheme = ComponentScheme.T(
1642 default='elastic10',
1643 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1645 stored_quantity = QuantityType.T(
1646 optional=True,
1647 help='Physical quantity of stored values (%s). If not given, a '
1648 'default is used based on the GF component scheme. The default '
1649 'for the ``"elastic*"`` family of component schemes is '
1650 '``"displacement"``.' % fmt_choices(QuantityType))
1652 tabulated_phases = List.T(
1653 TPDef.T(),
1654 help='Mapping of phase names to phase definitions, for which travel '
1655 'time tables are available in the GF store.')
1657 ncomponents = Int.T(
1658 optional=True,
1659 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1661 uuid = String.T(
1662 optional=True,
1663 help='Heuristic hash value which can be used to uniquely identify the '
1664 'GF store for practical purposes.')
1666 reference = String.T(
1667 optional=True,
1668 help="Store reference name composed of the store's :gattr:`id` and "
1669 'the first six letters of its :gattr:`uuid`.')
1671 def __init__(self, **kwargs):
1672 self._do_auto_updates = False
1673 Object.__init__(self, **kwargs)
1674 self._index_function = None
1675 self._indices_function = None
1676 self._vicinity_function = None
1677 self.validate(regularize=True, depth=1)
1678 self._do_auto_updates = True
1679 self.update()
1681 def check_ncomponents(self):
1682 ncomponents = component_scheme_to_description[
1683 self.component_scheme].ncomponents
1685 if self.ncomponents is None:
1686 self.ncomponents = ncomponents
1687 elif ncomponents != self.ncomponents:
1688 raise InvalidNComponents(
1689 'ncomponents=%i incompatible with component_scheme="%s"' % (
1690 self.ncomponents, self.component_scheme))
1692 def __setattr__(self, name, value):
1693 Object.__setattr__(self, name, value)
1694 try:
1695 self.T.get_property(name)
1696 if self._do_auto_updates:
1697 self.update()
1699 except ValueError:
1700 pass
1702 def update(self):
1703 self.check_ncomponents()
1704 self._update()
1705 self._make_index_functions()
1707 def irecord(self, *args):
1708 return self._index_function(*args)
1710 def irecords(self, *args):
1711 return self._indices_function(*args)
1713 def vicinity(self, *args):
1714 return self._vicinity_function(*args)
1716 def vicinities(self, *args):
1717 return self._vicinities_function(*args)
1719 def grid_interpolation_coefficients(self, *args):
1720 return self._grid_interpolation_coefficients(*args)
1722 def nodes(self, level=None, minlevel=None):
1723 return nodes(self.coords[minlevel:level])
1725 def iter_nodes(self, level=None, minlevel=None):
1726 return nditer_outer(self.coords[minlevel:level])
1728 def iter_extraction(self, gdef, level=None):
1729 i = 0
1730 arrs = []
1731 ntotal = 1
1732 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1733 if gdef and len(gdef) > i:
1734 sssn = gdef[i]
1735 else:
1736 sssn = (None,)*4
1738 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1739 ntotal *= len(arr)
1741 arrs.append(arr)
1742 i += 1
1744 arrs.append(self.coords[-1])
1745 return nditer_outer(arrs[:level])
1747 def make_sum_params(self, source, receiver, implementation='c',
1748 nthreads=0):
1749 assert implementation in ['c', 'python']
1751 out = []
1752 delays = source.times
1753 for comp, weights, icomponents in source.make_weights(
1754 receiver,
1755 self.component_scheme):
1757 weights *= self.factor
1759 args = self.make_indexing_args(source, receiver, icomponents)
1760 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1761 out.append((comp, args, delays_expanded, weights))
1763 return out
1765 def short_info(self):
1766 raise NotImplementedError('should be implemented in subclass')
1768 def get_shear_moduli(self, lat, lon, points,
1769 interpolation=None):
1770 '''
1771 Get shear moduli at given points from contained velocity model.
1773 :param lat: surface origin for coordinate system of ``points``
1774 :param points: NumPy array of shape ``(N, 3)``, where each row is
1775 a point ``(north, east, depth)``, relative to origin at
1776 ``(lat, lon)``
1777 :param interpolation: interpolation method. Choose from
1778 ``('nearest_neighbor', 'multilinear')``
1779 :returns: NumPy array of length N with extracted shear moduli at each
1780 point
1782 The default implementation retrieves and interpolates the shear moduli
1783 from the contained 1D velocity profile.
1784 '''
1785 return self.get_material_property(lat, lon, points,
1786 parameter='shear_moduli',
1787 interpolation=interpolation)
1789 def get_lambda_moduli(self, lat, lon, points,
1790 interpolation=None):
1791 '''
1792 Get lambda moduli at given points from contained velocity model.
1794 :param lat: surface origin for coordinate system of ``points``
1795 :param points: NumPy array of shape ``(N, 3)``, where each row is
1796 a point ``(north, east, depth)``, relative to origin at
1797 ``(lat, lon)``
1798 :param interpolation: interpolation method. Choose from
1799 ``('nearest_neighbor', 'multilinear')``
1800 :returns: NumPy array of length N with extracted shear moduli at each
1801 point
1803 The default implementation retrieves and interpolates the lambda moduli
1804 from the contained 1D velocity profile.
1805 '''
1806 return self.get_material_property(lat, lon, points,
1807 parameter='lambda_moduli',
1808 interpolation=interpolation)
1810 def get_bulk_moduli(self, lat, lon, points,
1811 interpolation=None):
1812 '''
1813 Get bulk moduli at given points from contained velocity model.
1815 :param lat: surface origin for coordinate system of ``points``
1816 :param points: NumPy array of shape ``(N, 3)``, where each row is
1817 a point ``(north, east, depth)``, relative to origin at
1818 ``(lat, lon)``
1819 :param interpolation: interpolation method. Choose from
1820 ``('nearest_neighbor', 'multilinear')``
1821 :returns: NumPy array of length N with extracted shear moduli at each
1822 point
1824 The default implementation retrieves and interpolates the lambda moduli
1825 from the contained 1D velocity profile.
1826 '''
1827 lambda_moduli = self.get_material_property(
1828 lat, lon, points, parameter='lambda_moduli',
1829 interpolation=interpolation)
1830 shear_moduli = self.get_material_property(
1831 lat, lon, points, parameter='shear_moduli',
1832 interpolation=interpolation)
1833 return lambda_moduli + (2 / 3) * shear_moduli
1835 def get_vs(self, lat, lon, points, interpolation=None):
1836 '''
1837 Get Vs at given points from contained velocity model.
1839 :param lat: surface origin for coordinate system of ``points``
1840 :param points: NumPy array of shape ``(N, 3)``, where each row is
1841 a point ``(north, east, depth)``, relative to origin at
1842 ``(lat, lon)``
1843 :param interpolation: interpolation method. Choose from
1844 ``('nearest_neighbor', 'multilinear')``
1845 :returns: NumPy array of length N with extracted shear moduli at each
1846 point
1848 The default implementation retrieves and interpolates Vs
1849 from the contained 1D velocity profile.
1850 '''
1851 return self.get_material_property(lat, lon, points,
1852 parameter='vs',
1853 interpolation=interpolation)
1855 def get_vp(self, lat, lon, points, interpolation=None):
1856 '''
1857 Get Vp at given points from contained velocity model.
1859 :param lat: surface origin for coordinate system of ``points``
1860 :param points: NumPy array of shape ``(N, 3)``, where each row is
1861 a point ``(north, east, depth)``, relative to origin at
1862 ``(lat, lon)``
1863 :param interpolation: interpolation method. Choose from
1864 ``('nearest_neighbor', 'multilinear')``
1865 :returns: NumPy array of length N with extracted shear moduli at each
1866 point
1868 The default implementation retrieves and interpolates Vp
1869 from the contained 1D velocity profile.
1870 '''
1871 return self.get_material_property(lat, lon, points,
1872 parameter='vp',
1873 interpolation=interpolation)
1875 def get_rho(self, lat, lon, points, interpolation=None):
1876 '''
1877 Get rho at given points from contained velocity model.
1879 :param lat: surface origin for coordinate system of ``points``
1880 :param points: NumPy array of shape ``(N, 3)``, where each row is
1881 a point ``(north, east, depth)``, relative to origin at
1882 ``(lat, lon)``
1883 :param interpolation: interpolation method. Choose from
1884 ``('nearest_neighbor', 'multilinear')``
1885 :returns: NumPy array of length N with extracted shear moduli at each
1886 point
1888 The default implementation retrieves and interpolates rho
1889 from the contained 1D velocity profile.
1890 '''
1891 return self.get_material_property(lat, lon, points,
1892 parameter='rho',
1893 interpolation=interpolation)
1895 def get_material_property(self, lat, lon, points, parameter='vs',
1896 interpolation=None):
1898 if interpolation is None:
1899 raise TypeError('Interpolation method not defined! available: '
1900 'multilinear', 'nearest_neighbor')
1902 earthmod = self.earthmodel_1d
1903 store_depth_profile = self.get_source_depths()
1904 z_profile = earthmod.profile('z')
1906 if parameter == 'vs':
1907 vs_profile = earthmod.profile('vs')
1908 profile = num.interp(
1909 store_depth_profile, z_profile, vs_profile)
1911 elif parameter == 'vp':
1912 vp_profile = earthmod.profile('vp')
1913 profile = num.interp(
1914 store_depth_profile, z_profile, vp_profile)
1916 elif parameter == 'rho':
1917 rho_profile = earthmod.profile('rho')
1919 profile = num.interp(
1920 store_depth_profile, z_profile, rho_profile)
1922 elif parameter == 'shear_moduli':
1923 vs_profile = earthmod.profile('vs')
1924 rho_profile = earthmod.profile('rho')
1926 store_vs_profile = num.interp(
1927 store_depth_profile, z_profile, vs_profile)
1928 store_rho_profile = num.interp(
1929 store_depth_profile, z_profile, rho_profile)
1931 profile = num.power(store_vs_profile, 2) * store_rho_profile
1933 elif parameter == 'lambda_moduli':
1934 vs_profile = earthmod.profile('vs')
1935 vp_profile = earthmod.profile('vp')
1936 rho_profile = earthmod.profile('rho')
1938 store_vs_profile = num.interp(
1939 store_depth_profile, z_profile, vs_profile)
1940 store_vp_profile = num.interp(
1941 store_depth_profile, z_profile, vp_profile)
1942 store_rho_profile = num.interp(
1943 store_depth_profile, z_profile, rho_profile)
1945 profile = store_rho_profile * (
1946 num.power(store_vp_profile, 2) -
1947 num.power(store_vs_profile, 2) * 2)
1948 else:
1949 raise TypeError(
1950 'parameter %s not available' % parameter)
1952 if interpolation == 'multilinear':
1953 kind = 'linear'
1954 elif interpolation == 'nearest_neighbor':
1955 kind = 'nearest'
1956 else:
1957 raise TypeError(
1958 'Interpolation method %s not available' % interpolation)
1960 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1962 try:
1963 return interpolator(points[:, 2])
1964 except ValueError:
1965 raise OutOfBounds()
1967 def is_static(self):
1968 for code in ('psgrn_pscmp', 'poel'):
1969 if self.modelling_code_id.startswith(code):
1970 return True
1971 return False
1973 def is_dynamic(self):
1974 return not self.is_static()
1976 def get_source_depths(self):
1977 raise NotImplementedError('must be implemented in subclass')
1979 def get_tabulated_phase(self, phase_id):
1980 '''
1981 Get tabulated phase definition.
1982 '''
1984 for pdef in self.tabulated_phases:
1985 if pdef.id == phase_id:
1986 return pdef
1988 raise StoreError('No such phase: %s' % phase_id)
1990 def fix_ttt_holes(self, sptree, mode):
1991 raise StoreError(
1992 'Cannot fix travel time table holes in GF stores of type %s.'
1993 % self.short_type)
1996class ConfigTypeA(Config):
1997 '''
1998 Cylindrical symmetry, 1D earth model, single receiver depth
2000 * Problem is invariant to horizontal translations and rotations around
2001 vertical axis.
2003 * All receivers must be at the same depth (e.g. at the surface)
2004 High level index variables: ``(source_depth, distance,
2005 component)``
2007 * The ``distance`` is the surface distance between source and receiver
2008 points.
2009 '''
2011 receiver_depth = Float.T(
2012 default=0.0,
2013 help='Fixed receiver depth [m].')
2015 source_depth_min = Float.T(
2016 help='Minimum source depth [m].')
2018 source_depth_max = Float.T(
2019 help='Maximum source depth [m].')
2021 source_depth_delta = Float.T(
2022 help='Grid spacing of source depths [m]')
2024 distance_min = Float.T(
2025 help='Minimum source-receiver surface distance [m].')
2027 distance_max = Float.T(
2028 help='Maximum source-receiver surface distance [m].')
2030 distance_delta = Float.T(
2031 help='Grid spacing of source-receiver surface distance [m].')
2033 short_type = 'A'
2035 provided_schemes = [
2036 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2038 def get_surface_distance(self, args):
2039 return args[1]
2041 def get_distance(self, args):
2042 return math.sqrt(args[0]**2 + args[1]**2)
2044 def get_source_depth(self, args):
2045 return args[0]
2047 def get_source_depths(self):
2048 return self.coords[0]
2050 def get_receiver_depth(self, args):
2051 return self.receiver_depth
2053 def _update(self):
2054 self.mins = num.array(
2055 [self.source_depth_min, self.distance_min], dtype=float)
2056 self.maxs = num.array(
2057 [self.source_depth_max, self.distance_max], dtype=float)
2058 self.deltas = num.array(
2059 [self.source_depth_delta, self.distance_delta],
2060 dtype=float)
2061 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2062 vicinity_eps).astype(int) + 1
2063 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2064 self.deltat = 1.0/self.sample_rate
2065 self.nrecords = num.prod(self.ns) * self.ncomponents
2066 self.coords = tuple(num.linspace(mi, ma, n) for
2067 (mi, ma, n) in
2068 zip(self.mins, self.effective_maxs, self.ns)) + \
2069 (num.arange(self.ncomponents),)
2071 self.nsource_depths, self.ndistances = self.ns
2073 def _make_index_functions(self):
2075 amin, bmin = self.mins
2076 da, db = self.deltas
2077 na, nb = self.ns
2079 ng = self.ncomponents
2081 def index_function(a, b, ig):
2082 ia = int(round((a - amin) / da))
2083 ib = int(round((b - bmin) / db))
2084 try:
2085 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2086 except ValueError:
2087 raise OutOfBounds()
2089 def indices_function(a, b, ig):
2090 ia = num.round((a - amin) / da).astype(int)
2091 ib = num.round((b - bmin) / db).astype(int)
2092 try:
2093 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2094 except ValueError:
2095 for ia_, ib_, ig_ in zip(ia, ib, ig):
2096 try:
2097 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2098 except ValueError:
2099 raise OutOfBounds()
2101 def grid_interpolation_coefficients(a, b):
2102 ias = indi12((a - amin) / da, na)
2103 ibs = indi12((b - bmin) / db, nb)
2104 return ias, ibs
2106 def vicinity_function(a, b, ig):
2107 ias, ibs = grid_interpolation_coefficients(a, b)
2109 if not (0 <= ig < ng):
2110 raise OutOfBounds()
2112 indis = []
2113 weights = []
2114 for ia, va in ias:
2115 iia = ia*nb*ng
2116 for ib, vb in ibs:
2117 indis.append(iia + ib*ng + ig)
2118 weights.append(va*vb)
2120 return num.array(indis), num.array(weights)
2122 def vicinities_function(a, b, ig):
2124 xa = (a - amin) / da
2125 xb = (b - bmin) / db
2127 xa_fl = num.floor(xa)
2128 xa_ce = num.ceil(xa)
2129 xb_fl = num.floor(xb)
2130 xb_ce = num.ceil(xb)
2131 va_fl = 1.0 - (xa - xa_fl)
2132 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2133 vb_fl = 1.0 - (xb - xb_fl)
2134 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2136 ia_fl = xa_fl.astype(int)
2137 ia_ce = xa_ce.astype(int)
2138 ib_fl = xb_fl.astype(int)
2139 ib_ce = xb_ce.astype(int)
2141 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2142 raise OutOfBounds()
2144 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2145 raise OutOfBounds()
2147 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2148 raise OutOfBounds()
2150 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2151 raise OutOfBounds()
2153 irecords = num.empty(a.size*4, dtype=int)
2154 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2155 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2156 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2157 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2159 weights = num.empty(a.size*4, dtype=float)
2160 weights[0::4] = va_fl * vb_fl
2161 weights[1::4] = va_ce * vb_fl
2162 weights[2::4] = va_fl * vb_ce
2163 weights[3::4] = va_ce * vb_ce
2165 return irecords, weights
2167 self._index_function = index_function
2168 self._indices_function = indices_function
2169 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2170 self._vicinity_function = vicinity_function
2171 self._vicinities_function = vicinities_function
2173 def make_indexing_args(self, source, receiver, icomponents):
2174 nc = icomponents.size
2175 dists = source.distances_to(receiver)
2176 n = dists.size
2177 return (num.tile(source.depths, nc//n),
2178 num.tile(dists, nc//n),
2179 icomponents)
2181 def make_indexing_args1(self, source, receiver):
2182 return (source.depth, source.distance_to(receiver))
2184 @property
2185 def short_extent(self):
2186 return '%g:%g:%g x %g:%g:%g' % (
2187 self.source_depth_min/km,
2188 self.source_depth_max/km,
2189 self.source_depth_delta/km,
2190 self.distance_min/km,
2191 self.distance_max/km,
2192 self.distance_delta/km)
2194 def fix_ttt_holes(self, sptree, mode):
2195 from pyrocko import eikonal_ext, spit
2197 nodes = self.nodes(level=-1)
2199 delta = self.deltas[-1]
2200 assert num.all(delta == self.deltas)
2202 nsources, ndistances = self.ns
2204 points = num.zeros((nodes.shape[0], 3))
2205 points[:, 0] = nodes[:, 1]
2206 points[:, 2] = nodes[:, 0]
2208 speeds = self.get_material_property(
2209 0., 0., points,
2210 parameter='vp' if mode == cake.P else 'vs',
2211 interpolation='multilinear')
2213 speeds = speeds.reshape((nsources, ndistances))
2215 times = sptree.interpolate_many(nodes)
2217 times[num.isnan(times)] = -1.
2218 times = times.reshape(speeds.shape)
2220 try:
2221 eikonal_ext.eikonal_solver_fmm_cartesian(
2222 speeds, times, delta)
2223 except eikonal_ext.EikonalExtError as e:
2224 if str(e).endswith('please check results'):
2225 logger.debug(
2226 'Got a warning from eikonal solver '
2227 '- may be ok...')
2228 else:
2229 raise
2231 def func(x):
2232 ibs, ics = \
2233 self.grid_interpolation_coefficients(*x)
2235 t = 0
2236 for ib, vb in ibs:
2237 for ic, vc in ics:
2238 t += times[ib, ic] * vb * vc
2240 return t
2242 return spit.SPTree(
2243 f=func,
2244 ftol=sptree.ftol,
2245 xbounds=sptree.xbounds,
2246 xtols=sptree.xtols)
2249class ConfigTypeB(Config):
2250 '''
2251 Cylindrical symmetry, 1D earth model, variable receiver depth
2253 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2254 receiver depth
2256 * High level index variables: ``(receiver_depth, source_depth,
2257 receiver_distance, component)``
2258 '''
2260 receiver_depth_min = Float.T(
2261 help='Minimum receiver depth [m].')
2263 receiver_depth_max = Float.T(
2264 help='Maximum receiver depth [m].')
2266 receiver_depth_delta = Float.T(
2267 help='Grid spacing of receiver depths [m]')
2269 source_depth_min = Float.T(
2270 help='Minimum source depth [m].')
2272 source_depth_max = Float.T(
2273 help='Maximum source depth [m].')
2275 source_depth_delta = Float.T(
2276 help='Grid spacing of source depths [m]')
2278 distance_min = Float.T(
2279 help='Minimum source-receiver surface distance [m].')
2281 distance_max = Float.T(
2282 help='Maximum source-receiver surface distance [m].')
2284 distance_delta = Float.T(
2285 help='Grid spacing of source-receiver surface distances [m].')
2287 short_type = 'B'
2289 provided_schemes = [
2290 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2292 def get_distance(self, args):
2293 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2295 def get_surface_distance(self, args):
2296 return args[2]
2298 def get_source_depth(self, args):
2299 return args[1]
2301 def get_receiver_depth(self, args):
2302 return args[0]
2304 def get_source_depths(self):
2305 return self.coords[1]
2307 def _update(self):
2308 self.mins = num.array([
2309 self.receiver_depth_min,
2310 self.source_depth_min,
2311 self.distance_min],
2312 dtype=float)
2314 self.maxs = num.array([
2315 self.receiver_depth_max,
2316 self.source_depth_max,
2317 self.distance_max],
2318 dtype=float)
2320 self.deltas = num.array([
2321 self.receiver_depth_delta,
2322 self.source_depth_delta,
2323 self.distance_delta],
2324 dtype=float)
2326 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2327 vicinity_eps).astype(int) + 1
2328 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2329 self.deltat = 1.0/self.sample_rate
2330 self.nrecords = num.prod(self.ns) * self.ncomponents
2331 self.coords = tuple(num.linspace(mi, ma, n) for
2332 (mi, ma, n) in
2333 zip(self.mins, self.effective_maxs, self.ns)) + \
2334 (num.arange(self.ncomponents),)
2335 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2337 def _make_index_functions(self):
2339 amin, bmin, cmin = self.mins
2340 da, db, dc = self.deltas
2341 na, nb, nc = self.ns
2342 ng = self.ncomponents
2344 def index_function(a, b, c, ig):
2345 ia = int(round((a - amin) / da))
2346 ib = int(round((b - bmin) / db))
2347 ic = int(round((c - cmin) / dc))
2348 try:
2349 return num.ravel_multi_index((ia, ib, ic, ig),
2350 (na, nb, nc, ng))
2351 except ValueError:
2352 raise OutOfBounds()
2354 def indices_function(a, b, c, ig):
2355 ia = num.round((a - amin) / da).astype(int)
2356 ib = num.round((b - bmin) / db).astype(int)
2357 ic = num.round((c - cmin) / dc).astype(int)
2358 try:
2359 return num.ravel_multi_index((ia, ib, ic, ig),
2360 (na, nb, nc, ng))
2361 except ValueError:
2362 raise OutOfBounds()
2364 def grid_interpolation_coefficients(a, b, c):
2365 ias = indi12((a - amin) / da, na)
2366 ibs = indi12((b - bmin) / db, nb)
2367 ics = indi12((c - cmin) / dc, nc)
2368 return ias, ibs, ics
2370 def vicinity_function(a, b, c, ig):
2371 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2373 if not (0 <= ig < ng):
2374 raise OutOfBounds()
2376 indis = []
2377 weights = []
2378 for ia, va in ias:
2379 iia = ia*nb*nc*ng
2380 for ib, vb in ibs:
2381 iib = ib*nc*ng
2382 for ic, vc in ics:
2383 indis.append(iia + iib + ic*ng + ig)
2384 weights.append(va*vb*vc)
2386 return num.array(indis), num.array(weights)
2388 def vicinities_function(a, b, c, ig):
2390 xa = (a - amin) / da
2391 xb = (b - bmin) / db
2392 xc = (c - cmin) / dc
2394 xa_fl = num.floor(xa)
2395 xa_ce = num.ceil(xa)
2396 xb_fl = num.floor(xb)
2397 xb_ce = num.ceil(xb)
2398 xc_fl = num.floor(xc)
2399 xc_ce = num.ceil(xc)
2400 va_fl = 1.0 - (xa - xa_fl)
2401 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2402 vb_fl = 1.0 - (xb - xb_fl)
2403 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2404 vc_fl = 1.0 - (xc - xc_fl)
2405 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2407 ia_fl = xa_fl.astype(int)
2408 ia_ce = xa_ce.astype(int)
2409 ib_fl = xb_fl.astype(int)
2410 ib_ce = xb_ce.astype(int)
2411 ic_fl = xc_fl.astype(int)
2412 ic_ce = xc_ce.astype(int)
2414 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2415 raise OutOfBounds()
2417 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2418 raise OutOfBounds()
2420 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2421 raise OutOfBounds()
2423 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2424 raise OutOfBounds()
2426 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2427 raise OutOfBounds()
2429 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2430 raise OutOfBounds()
2432 irecords = num.empty(a.size*8, dtype=int)
2433 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2434 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2435 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2436 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2437 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2438 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2439 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2440 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2442 weights = num.empty(a.size*8, dtype=float)
2443 weights[0::8] = va_fl * vb_fl * vc_fl
2444 weights[1::8] = va_ce * vb_fl * vc_fl
2445 weights[2::8] = va_fl * vb_ce * vc_fl
2446 weights[3::8] = va_ce * vb_ce * vc_fl
2447 weights[4::8] = va_fl * vb_fl * vc_ce
2448 weights[5::8] = va_ce * vb_fl * vc_ce
2449 weights[6::8] = va_fl * vb_ce * vc_ce
2450 weights[7::8] = va_ce * vb_ce * vc_ce
2452 return irecords, weights
2454 self._index_function = index_function
2455 self._indices_function = indices_function
2456 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2457 self._vicinity_function = vicinity_function
2458 self._vicinities_function = vicinities_function
2460 def make_indexing_args(self, source, receiver, icomponents):
2461 nc = icomponents.size
2462 dists = source.distances_to(receiver)
2463 n = dists.size
2464 receiver_depths = num.empty(nc)
2465 receiver_depths.fill(receiver.depth)
2466 return (receiver_depths,
2467 num.tile(source.depths, nc//n),
2468 num.tile(dists, nc//n),
2469 icomponents)
2471 def make_indexing_args1(self, source, receiver):
2472 return (receiver.depth,
2473 source.depth,
2474 source.distance_to(receiver))
2476 @property
2477 def short_extent(self):
2478 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2479 self.receiver_depth_min/km,
2480 self.receiver_depth_max/km,
2481 self.receiver_depth_delta/km,
2482 self.source_depth_min/km,
2483 self.source_depth_max/km,
2484 self.source_depth_delta/km,
2485 self.distance_min/km,
2486 self.distance_max/km,
2487 self.distance_delta/km)
2489 def fix_ttt_holes(self, sptree, mode):
2490 from pyrocko import eikonal_ext, spit
2492 nodes_sr = self.nodes(minlevel=1, level=-1)
2494 delta = self.deltas[-1]
2495 assert num.all(delta == self.deltas[1:])
2497 nreceivers, nsources, ndistances = self.ns
2499 points = num.zeros((nodes_sr.shape[0], 3))
2500 points[:, 0] = nodes_sr[:, 1]
2501 points[:, 2] = nodes_sr[:, 0]
2503 speeds = self.get_material_property(
2504 0., 0., points,
2505 parameter='vp' if mode == cake.P else 'vs',
2506 interpolation='multilinear')
2508 speeds = speeds.reshape((nsources, ndistances))
2510 receiver_times = []
2511 for ireceiver in range(nreceivers):
2512 nodes = num.hstack([
2513 num_full(
2514 (nodes_sr.shape[0], 1),
2515 self.coords[0][ireceiver]),
2516 nodes_sr])
2518 times = sptree.interpolate_many(nodes)
2520 times[num.isnan(times)] = -1.
2522 times = times.reshape(speeds.shape)
2524 try:
2525 eikonal_ext.eikonal_solver_fmm_cartesian(
2526 speeds, times, delta)
2527 except eikonal_ext.EikonalExtError as e:
2528 if str(e).endswith('please check results'):
2529 logger.debug(
2530 'Got a warning from eikonal solver '
2531 '- may be ok...')
2532 else:
2533 raise
2535 receiver_times.append(times)
2537 def func(x):
2538 ias, ibs, ics = \
2539 self.grid_interpolation_coefficients(*x)
2541 t = 0
2542 for ia, va in ias:
2543 times = receiver_times[ia]
2544 for ib, vb in ibs:
2545 for ic, vc in ics:
2546 t += times[ib, ic] * va * vb * vc
2548 return t
2550 return spit.SPTree(
2551 f=func,
2552 ftol=sptree.ftol,
2553 xbounds=sptree.xbounds,
2554 xtols=sptree.xtols)
2557class ConfigTypeC(Config):
2558 '''
2559 No symmetrical constraints, one fixed receiver position.
2561 * Cartesian 3D source volume around a reference point
2563 * High level index variables: ``(source_depth,
2564 source_east_shift, source_north_shift, component)``
2565 '''
2567 receiver = Receiver.T(
2568 help='Receiver location')
2570 source_origin = Location.T(
2571 help='Origin of the source volume grid.')
2573 source_depth_min = Float.T(
2574 help='Minimum source depth [m].')
2576 source_depth_max = Float.T(
2577 help='Maximum source depth [m].')
2579 source_depth_delta = Float.T(
2580 help='Source depth grid spacing [m].')
2582 source_east_shift_min = Float.T(
2583 help='Minimum easting of source grid [m].')
2585 source_east_shift_max = Float.T(
2586 help='Maximum easting of source grid [m].')
2588 source_east_shift_delta = Float.T(
2589 help='Source volume grid spacing in east direction [m].')
2591 source_north_shift_min = Float.T(
2592 help='Minimum northing of source grid [m].')
2594 source_north_shift_max = Float.T(
2595 help='Maximum northing of source grid [m].')
2597 source_north_shift_delta = Float.T(
2598 help='Source volume grid spacing in north direction [m].')
2600 short_type = 'C'
2602 provided_schemes = ['elastic18']
2604 def get_surface_distance(self, args):
2605 _, source_east_shift, source_north_shift, _ = args
2606 sorig = self.source_origin
2607 sloc = Location(
2608 lat=sorig.lat,
2609 lon=sorig.lon,
2610 north_shift=sorig.north_shift + source_north_shift,
2611 east_shift=sorig.east_shift + source_east_shift)
2613 return self.receiver.distance_to(sloc)
2615 def get_distance(self, args):
2616 sdepth, source_east_shift, source_north_shift, _ = args
2617 sorig = self.source_origin
2618 sloc = Location(
2619 lat=sorig.lat,
2620 lon=sorig.lon,
2621 north_shift=sorig.north_shift + source_north_shift,
2622 east_shift=sorig.east_shift + source_east_shift)
2624 return self.receiver.distance_3d_to(sloc)
2626 def get_source_depth(self, args):
2627 return args[0]
2629 def get_receiver_depth(self, args):
2630 return self.receiver.depth
2632 def get_source_depths(self):
2633 return self.coords[0]
2635 def _update(self):
2636 self.mins = num.array([
2637 self.source_depth_min,
2638 self.source_east_shift_min,
2639 self.source_north_shift_min],
2640 dtype=float)
2642 self.maxs = num.array([
2643 self.source_depth_max,
2644 self.source_east_shift_max,
2645 self.source_north_shift_max],
2646 dtype=float)
2648 self.deltas = num.array([
2649 self.source_depth_delta,
2650 self.source_east_shift_delta,
2651 self.source_north_shift_delta],
2652 dtype=float)
2654 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2655 vicinity_eps).astype(int) + 1
2656 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2657 self.deltat = 1.0/self.sample_rate
2658 self.nrecords = num.prod(self.ns) * self.ncomponents
2660 self.coords = tuple(
2661 num.linspace(mi, ma, n) for (mi, ma, n) in
2662 zip(self.mins, self.effective_maxs, self.ns)) + \
2663 (num.arange(self.ncomponents),)
2665 def _make_index_functions(self):
2667 amin, bmin, cmin = self.mins
2668 da, db, dc = self.deltas
2669 na, nb, nc = self.ns
2670 ng = self.ncomponents
2672 def index_function(a, b, c, ig):
2673 ia = int(round((a - amin) / da))
2674 ib = int(round((b - bmin) / db))
2675 ic = int(round((c - cmin) / dc))
2676 try:
2677 return num.ravel_multi_index((ia, ib, ic, ig),
2678 (na, nb, nc, ng))
2679 except ValueError:
2680 raise OutOfBounds(values=(a, b, c, ig))
2682 def indices_function(a, b, c, ig):
2683 ia = num.round((a - amin) / da).astype(int)
2684 ib = num.round((b - bmin) / db).astype(int)
2685 ic = num.round((c - cmin) / dc).astype(int)
2687 try:
2688 return num.ravel_multi_index((ia, ib, ic, ig),
2689 (na, nb, nc, ng))
2690 except ValueError:
2691 raise OutOfBounds()
2693 def vicinity_function(a, b, c, ig):
2694 ias = indi12((a - amin) / da, na)
2695 ibs = indi12((b - bmin) / db, nb)
2696 ics = indi12((c - cmin) / dc, nc)
2698 if not (0 <= ig < ng):
2699 raise OutOfBounds()
2701 indis = []
2702 weights = []
2703 for ia, va in ias:
2704 iia = ia*nb*nc*ng
2705 for ib, vb in ibs:
2706 iib = ib*nc*ng
2707 for ic, vc in ics:
2708 indis.append(iia + iib + ic*ng + ig)
2709 weights.append(va*vb*vc)
2711 return num.array(indis), num.array(weights)
2713 def vicinities_function(a, b, c, ig):
2715 xa = (a-amin) / da
2716 xb = (b-bmin) / db
2717 xc = (c-cmin) / dc
2719 xa_fl = num.floor(xa)
2720 xa_ce = num.ceil(xa)
2721 xb_fl = num.floor(xb)
2722 xb_ce = num.ceil(xb)
2723 xc_fl = num.floor(xc)
2724 xc_ce = num.ceil(xc)
2725 va_fl = 1.0 - (xa - xa_fl)
2726 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2727 vb_fl = 1.0 - (xb - xb_fl)
2728 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2729 vc_fl = 1.0 - (xc - xc_fl)
2730 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2732 ia_fl = xa_fl.astype(int)
2733 ia_ce = xa_ce.astype(int)
2734 ib_fl = xb_fl.astype(int)
2735 ib_ce = xb_ce.astype(int)
2736 ic_fl = xc_fl.astype(int)
2737 ic_ce = xc_ce.astype(int)
2739 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2740 raise OutOfBounds()
2742 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2743 raise OutOfBounds()
2745 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2746 raise OutOfBounds()
2748 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2749 raise OutOfBounds()
2751 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2752 raise OutOfBounds()
2754 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2755 raise OutOfBounds()
2757 irecords = num.empty(a.size*8, dtype=int)
2758 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2759 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2760 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2761 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2762 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2763 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2764 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2765 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2767 weights = num.empty(a.size*8, dtype=float)
2768 weights[0::8] = va_fl * vb_fl * vc_fl
2769 weights[1::8] = va_ce * vb_fl * vc_fl
2770 weights[2::8] = va_fl * vb_ce * vc_fl
2771 weights[3::8] = va_ce * vb_ce * vc_fl
2772 weights[4::8] = va_fl * vb_fl * vc_ce
2773 weights[5::8] = va_ce * vb_fl * vc_ce
2774 weights[6::8] = va_fl * vb_ce * vc_ce
2775 weights[7::8] = va_ce * vb_ce * vc_ce
2777 return irecords, weights
2779 self._index_function = index_function
2780 self._indices_function = indices_function
2781 self._vicinity_function = vicinity_function
2782 self._vicinities_function = vicinities_function
2784 def make_indexing_args(self, source, receiver, icomponents):
2785 nc = icomponents.size
2787 dists = source.distances_to(self.source_origin)
2788 azis, _ = source.azibazis_to(self.source_origin)
2790 source_north_shifts = - num.cos(d2r*azis) * dists
2791 source_east_shifts = - num.sin(d2r*azis) * dists
2792 source_depths = source.depths - self.source_origin.depth
2794 n = dists.size
2796 return (num.tile(source_depths, nc//n),
2797 num.tile(source_east_shifts, nc//n),
2798 num.tile(source_north_shifts, nc//n),
2799 icomponents)
2801 def make_indexing_args1(self, source, receiver):
2802 dist = source.distance_to(self.source_origin)
2803 azi, _ = source.azibazi_to(self.source_origin)
2805 source_north_shift = - num.cos(d2r*azi) * dist
2806 source_east_shift = - num.sin(d2r*azi) * dist
2807 source_depth = source.depth - self.source_origin.depth
2809 return (source_depth,
2810 source_east_shift,
2811 source_north_shift)
2813 @property
2814 def short_extent(self):
2815 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2816 self.source_depth_min/km,
2817 self.source_depth_max/km,
2818 self.source_depth_delta/km,
2819 self.source_east_shift_min/km,
2820 self.source_east_shift_max/km,
2821 self.source_east_shift_delta/km,
2822 self.source_north_shift_min/km,
2823 self.source_north_shift_max/km,
2824 self.source_north_shift_delta/km)
2827class Weighting(Object):
2828 factor = Float.T(default=1.0)
2831class Taper(Object):
2832 tmin = Timing.T()
2833 tmax = Timing.T()
2834 tfade = Float.T(default=0.0)
2835 shape = StringChoice.T(
2836 choices=['cos', 'linear'],
2837 default='cos',
2838 optional=True)
2841class SimplePattern(SObject):
2843 _pool = {}
2845 def __init__(self, pattern):
2846 self._pattern = pattern
2847 SObject.__init__(self)
2849 def __str__(self):
2850 return self._pattern
2852 @property
2853 def regex(self):
2854 pool = SimplePattern._pool
2855 if self.pattern not in pool:
2856 rpat = '|'.join(fnmatch.translate(x) for
2857 x in self.pattern.split('|'))
2858 pool[self.pattern] = re.compile(rpat, re.I)
2860 return pool[self.pattern]
2862 def match(self, s):
2863 return self.regex.match(s)
2866class WaveformType(StringChoice):
2867 choices = ['dis', 'vel', 'acc',
2868 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2869 'envelope_dis', 'envelope_vel', 'envelope_acc']
2872class ChannelSelection(Object):
2873 pattern = SimplePattern.T(optional=True)
2874 min_sample_rate = Float.T(optional=True)
2875 max_sample_rate = Float.T(optional=True)
2878class StationSelection(Object):
2879 includes = SimplePattern.T()
2880 excludes = SimplePattern.T()
2881 distance_min = Float.T(optional=True)
2882 distance_max = Float.T(optional=True)
2883 azimuth_min = Float.T(optional=True)
2884 azimuth_max = Float.T(optional=True)
2887class WaveformSelection(Object):
2888 channel_selection = ChannelSelection.T(optional=True)
2889 station_selection = StationSelection.T(optional=True)
2890 taper = Taper.T()
2891 # filter = FrequencyResponse.T()
2892 waveform_type = WaveformType.T(default='dis')
2893 weighting = Weighting.T(optional=True)
2894 sample_rate = Float.T(optional=True)
2895 gf_store_id = StringID.T(optional=True)
2898def indi12(x, n):
2899 '''
2900 Get linear interpolation index and weight.
2901 '''
2903 r = round(x)
2904 if abs(r - x) < vicinity_eps:
2905 i = int(r)
2906 if not (0 <= i < n):
2907 raise OutOfBounds()
2909 return ((int(r), 1.),)
2910 else:
2911 f = math.floor(x)
2912 i = int(f)
2913 if not (0 <= i < n-1):
2914 raise OutOfBounds()
2916 v = x-f
2917 return ((i, 1.-v), (i + 1, v))
2920def float_or_none(s):
2921 units = {
2922 'k': 1e3,
2923 'M': 1e6,
2924 }
2926 factor = 1.0
2927 if s and s[-1] in units:
2928 factor = units[s[-1]]
2929 s = s[:-1]
2930 if not s:
2931 raise ValueError("unit without a number: '%s'" % s)
2933 if s:
2934 return float(s) * factor
2935 else:
2936 return None
2939class GridSpecError(Exception):
2940 def __init__(self, s):
2941 Exception.__init__(self, 'invalid grid specification: %s' % s)
2944def parse_grid_spec(spec):
2945 try:
2946 result = []
2947 for dspec in spec.split(','):
2948 t = dspec.split('@')
2949 num = start = stop = step = None
2950 if len(t) == 2:
2951 num = int(t[1])
2952 if num <= 0:
2953 raise GridSpecError(spec)
2955 elif len(t) > 2:
2956 raise GridSpecError(spec)
2958 s = t[0]
2959 v = [float_or_none(x) for x in s.split(':')]
2960 if len(v) == 1:
2961 start = stop = v[0]
2962 if len(v) >= 2:
2963 start, stop = v[0:2]
2964 if len(v) == 3:
2965 step = v[2]
2967 if len(v) > 3 or (len(v) > 2 and num is not None):
2968 raise GridSpecError(spec)
2970 if step == 0.0:
2971 raise GridSpecError(spec)
2973 result.append((start, stop, step, num))
2975 except ValueError:
2976 raise GridSpecError(spec)
2978 return result
2981def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2982 swap = step is not None and step < 0.
2983 if start is None:
2984 start = ma if swap else mi
2985 if stop is None:
2986 stop = mi if swap else ma
2987 if step is None:
2988 step = -inc if ma < mi else inc
2989 if num is None:
2990 if (step < 0) != (stop-start < 0):
2991 raise GridSpecError()
2993 num = int(round((stop-start)/step))+1
2994 stop2 = start + (num-1)*step
2995 if abs(stop-stop2) > eps:
2996 num = int(math.floor((stop-start)/step))+1
2997 stop = start + (num-1)*step
2998 else:
2999 stop = stop2
3001 if start == stop:
3002 num = 1
3004 return start, stop, num
3007def nditer_outer(x):
3008 return num.nditer(
3009 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
3012def nodes(xs):
3013 ns = [x.size for x in xs]
3014 nnodes = num.prod(ns)
3015 ndim = len(xs)
3016 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
3017 for idim in range(ndim-1, -1, -1):
3018 x = xs[idim]
3019 nrepeat = num.prod(ns[idim+1:], dtype=int)
3020 ntile = num.prod(ns[:idim], dtype=int)
3021 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
3023 return nodes
3026def filledi(x, n):
3027 a = num.empty(n, dtype=int)
3028 a.fill(x)
3029 return a
3032config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3034discretized_source_classes = [
3035 DiscretizedExplosionSource,
3036 DiscretizedSFSource,
3037 DiscretizedMTSource,
3038 DiscretizedPorePressureSource]
3041__all__ = '''
3042Earthmodel1D
3043StringID
3044ScopeType
3045WaveformType
3046QuantityType
3047NearfieldTermsType
3048Reference
3049Region
3050CircularRegion
3051RectangularRegion
3052PhaseSelect
3053InvalidTimingSpecification
3054Timing
3055TPDef
3056OutOfBounds
3057Location
3058Receiver
3059'''.split() + [
3060 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3061ComponentScheme
3062component_scheme_to_description
3063component_schemes
3064Config
3065GridSpecError
3066Weighting
3067Taper
3068SimplePattern
3069WaveformType
3070ChannelSelection
3071StationSelection
3072WaveformSelection
3073nditer_outer
3074dump
3075load
3076discretized_source_classes
3077config_type_classes
3078UnavailableScheme
3079InterpolationMethod
3080SeismosizerTrace
3081SeismosizerResult
3082Result
3083StaticResult
3084TimingAttributeError
3085'''.split()