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 depths = cat([s.depths for s in sources])
996 if same_ref:
997 lat = first.lat
998 lon = first.lon
999 north_shifts = cat([s.effective_north_shifts for s in sources])
1000 east_shifts = cat([s.effective_east_shifts for s in sources])
1001 lats = None
1002 lons = None
1003 else:
1004 lat = None
1005 lon = None
1006 north_shifts = None
1007 east_shifts = None
1009 return cls(
1010 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
1011 north_shifts=north_shifts, east_shifts=east_shifts,
1012 depths=depths, **kwargs)
1014 def centroid_position(self):
1015 moments = self.moments()
1016 norm = num.sum(moments)
1017 if norm != 0.0:
1018 w = moments / num.sum(moments)
1019 else:
1020 w = num.ones(moments.size)
1022 if self.lats is not None and self.lons is not None:
1023 lats, lons = self.effective_latlons
1024 rlat, rlon = num.mean(lats), num.mean(lons)
1025 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
1026 else:
1027 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
1028 n, e = self.north_shifts, self.east_shifts
1030 cn = num.sum(n*w)
1031 ce = num.sum(e*w)
1032 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1034 if self.lats is not None and self.lons is not None:
1035 lat = clat
1036 lon = clon
1037 north_shift = 0.
1038 east_shift = 0.
1039 else:
1040 lat = g(self.lat, 0.0)
1041 lon = g(self.lon, 0.0)
1042 north_shift = cn
1043 east_shift = ce
1045 depth = num.sum(self.depths*w)
1046 time = num.sum(self.times*w)
1047 return tuple(float(x) for x in
1048 (time, lat, lon, north_shift, east_shift, depth))
1051class DiscretizedExplosionSource(DiscretizedSource):
1052 m0s = Array.T(shape=(None,), dtype=float)
1054 provided_schemes = (
1055 'elastic2',
1056 'elastic8',
1057 'elastic10',
1058 )
1060 def get_source_terms(self, scheme):
1061 self.check_scheme(scheme)
1063 if scheme == 'elastic2':
1064 return self.m0s[:, num.newaxis].copy()
1066 elif scheme in ('elastic8', 'elastic10'):
1067 m6s = num.zeros((self.m0s.size, 6))
1068 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1069 return m6s
1070 else:
1071 assert False
1073 def make_weights(self, receiver, scheme):
1074 self.check_scheme(scheme)
1076 azis, bazis = self.azibazis_to(receiver)
1078 sb = num.sin(bazis*d2r-num.pi)
1079 cb = num.cos(bazis*d2r-num.pi)
1081 m0s = self.m0s
1082 n = azis.size
1084 cat = num.concatenate
1085 rep = num.repeat
1087 if scheme == 'elastic2':
1088 w_n = cb*m0s
1089 g_n = filledi(0, n)
1090 w_e = sb*m0s
1091 g_e = filledi(0, n)
1092 w_d = m0s
1093 g_d = filledi(1, n)
1095 elif scheme == 'elastic8':
1096 w_n = cat((cb*m0s, cb*m0s))
1097 g_n = rep((0, 2), n)
1098 w_e = cat((sb*m0s, sb*m0s))
1099 g_e = rep((0, 2), n)
1100 w_d = cat((m0s, m0s))
1101 g_d = rep((5, 7), n)
1103 elif scheme == 'elastic10':
1104 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1105 g_n = rep((0, 2, 8), n)
1106 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1107 g_e = rep((0, 2, 8), n)
1108 w_d = cat((m0s, m0s, m0s))
1109 g_d = rep((5, 7, 9), n)
1111 else:
1112 assert False
1114 return (
1115 ('displacement.n', w_n, g_n),
1116 ('displacement.e', w_e, g_e),
1117 ('displacement.d', w_d, g_d),
1118 )
1120 def split(self):
1121 from pyrocko.gf.seismosizer import ExplosionSource
1122 sources = []
1123 for i in range(self.nelements):
1124 lat, lon, north_shift, east_shift = self.element_coords(i)
1125 sources.append(ExplosionSource(
1126 time=float(self.times[i]),
1127 lat=lat,
1128 lon=lon,
1129 north_shift=north_shift,
1130 east_shift=east_shift,
1131 depth=float(self.depths[i]),
1132 moment=float(self.m0s[i])))
1134 return sources
1136 def moments(self):
1137 return self.m0s
1139 def centroid(self):
1140 from pyrocko.gf.seismosizer import ExplosionSource
1141 time, lat, lon, north_shift, east_shift, depth = \
1142 self.centroid_position()
1144 return ExplosionSource(
1145 time=time,
1146 lat=lat,
1147 lon=lon,
1148 north_shift=north_shift,
1149 east_shift=east_shift,
1150 depth=depth,
1151 moment=float(num.sum(self.m0s)))
1153 @classmethod
1154 def combine(cls, sources, **kwargs):
1155 '''
1156 Combine several discretized source models.
1158 Concatenenates all point sources in the given discretized ``sources``.
1159 Care must be taken when using this function that the external amplitude
1160 factors and reference times of the parameterized (undiscretized)
1161 sources match or are accounted for.
1162 '''
1164 if 'm0s' not in kwargs:
1165 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1167 return super(DiscretizedExplosionSource, cls).combine(sources,
1168 **kwargs)
1171class DiscretizedSFSource(DiscretizedSource):
1172 forces = Array.T(shape=(None, 3), dtype=float)
1174 provided_schemes = (
1175 'elastic5',
1176 )
1178 def get_source_terms(self, scheme):
1179 self.check_scheme(scheme)
1181 return self.forces
1183 def make_weights(self, receiver, scheme):
1184 self.check_scheme(scheme)
1186 azis, bazis = self.azibazis_to(receiver)
1188 sa = num.sin(azis*d2r)
1189 ca = num.cos(azis*d2r)
1190 sb = num.sin(bazis*d2r-num.pi)
1191 cb = num.cos(bazis*d2r-num.pi)
1193 forces = self.forces
1194 fn = forces[:, 0]
1195 fe = forces[:, 1]
1196 fd = forces[:, 2]
1198 f0 = fd
1199 f1 = ca * fn + sa * fe
1200 f2 = ca * fe - sa * fn
1202 n = azis.size
1204 if scheme == 'elastic5':
1205 ioff = 0
1207 cat = num.concatenate
1208 rep = num.repeat
1210 w_n = cat((cb*f0, cb*f1, -sb*f2))
1211 g_n = ioff + rep((0, 1, 2), n)
1212 w_e = cat((sb*f0, sb*f1, cb*f2))
1213 g_e = ioff + rep((0, 1, 2), n)
1214 w_d = cat((f0, f1))
1215 g_d = ioff + rep((3, 4), n)
1217 return (
1218 ('displacement.n', w_n, g_n),
1219 ('displacement.e', w_e, g_e),
1220 ('displacement.d', w_d, g_d),
1221 )
1223 @classmethod
1224 def combine(cls, sources, **kwargs):
1225 '''
1226 Combine several discretized source models.
1228 Concatenenates all point sources in the given discretized ``sources``.
1229 Care must be taken when using this function that the external amplitude
1230 factors and reference times of the parameterized (undiscretized)
1231 sources match or are accounted for.
1232 '''
1234 if 'forces' not in kwargs:
1235 kwargs['forces'] = num.vstack([s.forces for s in sources])
1237 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1239 def moments(self):
1240 return num.sum(self.forces**2, axis=1)
1242 def centroid(self):
1243 from pyrocko.gf.seismosizer import SFSource
1244 time, lat, lon, north_shift, east_shift, depth = \
1245 self.centroid_position()
1247 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1248 return SFSource(
1249 time=time,
1250 lat=lat,
1251 lon=lon,
1252 north_shift=north_shift,
1253 east_shift=east_shift,
1254 depth=depth,
1255 fn=fn,
1256 fe=fe,
1257 fd=fd)
1260class DiscretizedMTSource(DiscretizedSource):
1261 m6s = Array.T(
1262 shape=(None, 6), dtype=float,
1263 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1265 provided_schemes = (
1266 'elastic8',
1267 'elastic10',
1268 'elastic18',
1269 )
1271 def get_source_terms(self, scheme):
1272 self.check_scheme(scheme)
1273 return self.m6s
1275 def make_weights(self, receiver, scheme):
1276 self.check_scheme(scheme)
1278 m6s = self.m6s
1279 n = m6s.shape[0]
1280 rep = num.repeat
1282 if scheme == 'elastic18':
1283 w_n = m6s.flatten()
1284 g_n = num.tile(num.arange(0, 6), n)
1285 w_e = m6s.flatten()
1286 g_e = num.tile(num.arange(6, 12), n)
1287 w_d = m6s.flatten()
1288 g_d = num.tile(num.arange(12, 18), n)
1290 else:
1291 azis, bazis = self.azibazis_to(receiver)
1293 sa = num.sin(azis*d2r)
1294 ca = num.cos(azis*d2r)
1295 s2a = num.sin(2.*azis*d2r)
1296 c2a = num.cos(2.*azis*d2r)
1297 sb = num.sin(bazis*d2r-num.pi)
1298 cb = num.cos(bazis*d2r-num.pi)
1300 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1301 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1302 f2 = m6s[:, 2]
1303 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1304 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1305 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1307 cat = num.concatenate
1309 if scheme == 'elastic8':
1310 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1311 g_n = rep((0, 1, 2, 3, 4), n)
1312 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1313 g_e = rep((0, 1, 2, 3, 4), n)
1314 w_d = cat((f0, f1, f2))
1315 g_d = rep((5, 6, 7), n)
1317 elif scheme == 'elastic10':
1318 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1319 g_n = rep((0, 1, 2, 8, 3, 4), n)
1320 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1321 g_e = rep((0, 1, 2, 8, 3, 4), n)
1322 w_d = cat((f0, f1, f2, f5))
1323 g_d = rep((5, 6, 7, 9), n)
1325 else:
1326 assert False
1328 return (
1329 ('displacement.n', w_n, g_n),
1330 ('displacement.e', w_e, g_e),
1331 ('displacement.d', w_d, g_d),
1332 )
1334 def split(self):
1335 from pyrocko.gf.seismosizer import MTSource
1336 sources = []
1337 for i in range(self.nelements):
1338 lat, lon, north_shift, east_shift = self.element_coords(i)
1339 sources.append(MTSource(
1340 time=float(self.times[i]),
1341 lat=lat,
1342 lon=lon,
1343 north_shift=north_shift,
1344 east_shift=east_shift,
1345 depth=float(self.depths[i]),
1346 m6=self.m6s[i]))
1348 return sources
1350 def moments(self):
1351 moments = num.array(
1352 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))
1353 for m6 in self.m6s])
1354 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)
1356 def get_moment_rate(self, deltat=None):
1357 moments = self.moments()
1358 times = self.times
1359 times -= times.min()
1361 t_max = times.max()
1362 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat
1363 mom_times[mom_times > t_max] = t_max
1365 # Right open histrogram bins
1366 mom, _ = num.histogram(
1367 times,
1368 bins=mom_times,
1369 weights=moments)
1371 deltat = num.diff(mom_times)
1372 mom_rate = mom / deltat
1374 return mom_rate, mom_times[1:]
1376 def centroid(self):
1377 from pyrocko.gf.seismosizer import MTSource
1378 time, lat, lon, north_shift, east_shift, depth = \
1379 self.centroid_position()
1381 return MTSource(
1382 time=time,
1383 lat=lat,
1384 lon=lon,
1385 north_shift=north_shift,
1386 east_shift=east_shift,
1387 depth=depth,
1388 m6=num.sum(self.m6s, axis=0))
1390 @classmethod
1391 def combine(cls, sources, **kwargs):
1392 '''
1393 Combine several discretized source models.
1395 Concatenenates all point sources in the given discretized ``sources``.
1396 Care must be taken when using this function that the external amplitude
1397 factors and reference times of the parameterized (undiscretized)
1398 sources match or are accounted for.
1399 '''
1401 if 'm6s' not in kwargs:
1402 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1404 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1407class DiscretizedPorePressureSource(DiscretizedSource):
1408 pp = Array.T(shape=(None,), dtype=float)
1410 provided_schemes = (
1411 'poroelastic10',
1412 )
1414 def get_source_terms(self, scheme):
1415 self.check_scheme(scheme)
1416 return self.pp[:, num.newaxis].copy()
1418 def make_weights(self, receiver, scheme):
1419 self.check_scheme(scheme)
1421 azis, bazis = self.azibazis_to(receiver)
1423 sb = num.sin(bazis*d2r-num.pi)
1424 cb = num.cos(bazis*d2r-num.pi)
1426 pp = self.pp
1427 n = bazis.size
1429 w_un = cb*pp
1430 g_un = filledi(1, n)
1431 w_ue = sb*pp
1432 g_ue = filledi(1, n)
1433 w_ud = pp
1434 g_ud = filledi(0, n)
1436 w_tn = cb*pp
1437 g_tn = filledi(6, n)
1438 w_te = sb*pp
1439 g_te = filledi(6, n)
1441 w_pp = pp
1442 g_pp = filledi(7, n)
1444 w_dvn = cb*pp
1445 g_dvn = filledi(9, n)
1446 w_dve = sb*pp
1447 g_dve = filledi(9, n)
1448 w_dvd = pp
1449 g_dvd = filledi(8, n)
1451 return (
1452 ('displacement.n', w_un, g_un),
1453 ('displacement.e', w_ue, g_ue),
1454 ('displacement.d', w_ud, g_ud),
1455 ('vertical_tilt.n', w_tn, g_tn),
1456 ('vertical_tilt.e', w_te, g_te),
1457 ('pore_pressure', w_pp, g_pp),
1458 ('darcy_velocity.n', w_dvn, g_dvn),
1459 ('darcy_velocity.e', w_dve, g_dve),
1460 ('darcy_velocity.d', w_dvd, g_dvd),
1461 )
1463 def moments(self):
1464 return self.pp
1466 def centroid(self):
1467 from pyrocko.gf.seismosizer import PorePressurePointSource
1468 time, lat, lon, north_shift, east_shift, depth = \
1469 self.centroid_position()
1471 return PorePressurePointSource(
1472 time=time,
1473 lat=lat,
1474 lon=lon,
1475 north_shift=north_shift,
1476 east_shift=east_shift,
1477 depth=depth,
1478 pp=float(num.sum(self.pp)))
1480 @classmethod
1481 def combine(cls, sources, **kwargs):
1482 '''
1483 Combine several discretized source models.
1485 Concatenenates all point sources in the given discretized ``sources``.
1486 Care must be taken when using this function that the external amplitude
1487 factors and reference times of the parameterized (undiscretized)
1488 sources match or are accounted for.
1489 '''
1491 if 'pp' not in kwargs:
1492 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1494 return super(DiscretizedPorePressureSource, cls).combine(sources,
1495 **kwargs)
1498class Region(Object):
1499 name = String.T(optional=True)
1502class RectangularRegion(Region):
1503 lat_min = Float.T()
1504 lat_max = Float.T()
1505 lon_min = Float.T()
1506 lon_max = Float.T()
1509class CircularRegion(Region):
1510 lat = Float.T()
1511 lon = Float.T()
1512 radius = Float.T()
1515class Config(Object):
1516 '''
1517 Green's function store meta information.
1519 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1520 configuration types are:
1522 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1523 spherical symmetry, 1D earth model, single receiver depth
1525 * Problem is invariant to horizontal translations and rotations around
1526 vertical axis.
1527 * All receivers must be at the same depth (e.g. at the surface)
1528 * High level index variables: ``(source_depth, receiver_distance,
1529 component)``
1531 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1532 spherical symmetry, 1D earth model, variable receiver depth
1534 * Symmetries like in Type A but has additional index for receiver depth
1535 * High level index variables: ``(source_depth, receiver_distance,
1536 receiver_depth, component)``
1538 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1539 constraints but fixed receiver positions
1541 * Cartesian source volume around a reference point
1542 * High level index variables: ``(ireceiver, source_depth,
1543 source_east_shift, source_north_shift, component)``
1544 '''
1546 id = StringID.T(
1547 help='Name of the store. May consist of upper and lower-case letters, '
1548 'digits, dots and underscores. The name must start with a '
1549 'letter.')
1551 derived_from_id = StringID.T(
1552 optional=True,
1553 help='Name of the original store, if this store has been derived from '
1554 'another one (e.g. extracted subset).')
1556 version = String.T(
1557 default='1.0',
1558 optional=True,
1559 help='User-defined version string. Use <major>.<minor> format.')
1561 modelling_code_id = StringID.T(
1562 optional=True,
1563 help='Identifier of the backend used to compute the store.')
1565 author = Unicode.T(
1566 optional=True,
1567 help='Comma-separated list of author names.')
1569 author_email = String.T(
1570 optional=True,
1571 help="Author's contact email address.")
1573 created_time = Timestamp.T(
1574 optional=True,
1575 help='Time of creation of the store.')
1577 regions = List.T(
1578 Region.T(),
1579 help='Geographical regions for which the store is representative.')
1581 scope_type = ScopeType.T(
1582 optional=True,
1583 help='Distance range scope of the store '
1584 '(%s).' % fmt_choices(ScopeType))
1586 waveform_type = WaveType.T(
1587 optional=True,
1588 help='Wave type stored (%s).' % fmt_choices(WaveType))
1590 nearfield_terms = NearfieldTermsType.T(
1591 optional=True,
1592 help='Information about the inclusion of near-field terms in the '
1593 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1595 description = String.T(
1596 optional=True,
1597 help='Free form textual description of the GF store.')
1599 references = List.T(
1600 Reference.T(),
1601 help='Reference list to cite the modelling code, earth model or '
1602 'related work.')
1604 earthmodel_1d = Earthmodel1D.T(
1605 optional=True,
1606 help='Layered earth model in ND (named discontinuity) format.')
1608 earthmodel_receiver_1d = Earthmodel1D.T(
1609 optional=True,
1610 help='Receiver-side layered earth model in ND format.')
1612 can_interpolate_source = Bool.T(
1613 optional=True,
1614 help='Hint to indicate if the spatial sampling of the store is dense '
1615 'enough for multi-linear interpolation at the source.')
1617 can_interpolate_receiver = Bool.T(
1618 optional=True,
1619 help='Hint to indicate if the spatial sampling of the store is dense '
1620 'enough for multi-linear interpolation at the receiver.')
1622 frequency_min = Float.T(
1623 optional=True,
1624 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1626 frequency_max = Float.T(
1627 optional=True,
1628 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1630 sample_rate = Float.T(
1631 optional=True,
1632 help='Sample rate of the GF store [Hz].')
1634 factor = Float.T(
1635 default=1.0,
1636 help='Gain value, factored out of the stored GF samples. '
1637 '(may not work properly, keep at 1.0).',
1638 optional=True)
1640 component_scheme = ComponentScheme.T(
1641 default='elastic10',
1642 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1644 stored_quantity = QuantityType.T(
1645 optional=True,
1646 help='Physical quantity of stored values (%s). If not given, a '
1647 'default is used based on the GF component scheme. The default '
1648 'for the ``"elastic*"`` family of component schemes is '
1649 '``"displacement"``.' % fmt_choices(QuantityType))
1651 tabulated_phases = List.T(
1652 TPDef.T(),
1653 help='Mapping of phase names to phase definitions, for which travel '
1654 'time tables are available in the GF store.')
1656 ncomponents = Int.T(
1657 optional=True,
1658 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1660 uuid = String.T(
1661 optional=True,
1662 help='Heuristic hash value which can be used to uniquely identify the '
1663 'GF store for practical purposes.')
1665 reference = String.T(
1666 optional=True,
1667 help="Store reference name composed of the store's :gattr:`id` and "
1668 'the first six letters of its :gattr:`uuid`.')
1670 def __init__(self, **kwargs):
1671 self._do_auto_updates = False
1672 Object.__init__(self, **kwargs)
1673 self._index_function = None
1674 self._indices_function = None
1675 self._vicinity_function = None
1676 self.validate(regularize=True, depth=1)
1677 self._do_auto_updates = True
1678 self.update()
1680 def check_ncomponents(self):
1681 ncomponents = component_scheme_to_description[
1682 self.component_scheme].ncomponents
1684 if self.ncomponents is None:
1685 self.ncomponents = ncomponents
1686 elif ncomponents != self.ncomponents:
1687 raise InvalidNComponents(
1688 'ncomponents=%i incompatible with component_scheme="%s"' % (
1689 self.ncomponents, self.component_scheme))
1691 def __setattr__(self, name, value):
1692 Object.__setattr__(self, name, value)
1693 try:
1694 self.T.get_property(name)
1695 if self._do_auto_updates:
1696 self.update()
1698 except ValueError:
1699 pass
1701 def update(self):
1702 self.check_ncomponents()
1703 self._update()
1704 self._make_index_functions()
1706 def irecord(self, *args):
1707 return self._index_function(*args)
1709 def irecords(self, *args):
1710 return self._indices_function(*args)
1712 def vicinity(self, *args):
1713 return self._vicinity_function(*args)
1715 def vicinities(self, *args):
1716 return self._vicinities_function(*args)
1718 def grid_interpolation_coefficients(self, *args):
1719 return self._grid_interpolation_coefficients(*args)
1721 def nodes(self, level=None, minlevel=None):
1722 return nodes(self.coords[minlevel:level])
1724 def iter_nodes(self, level=None, minlevel=None):
1725 return nditer_outer(self.coords[minlevel:level])
1727 def iter_extraction(self, gdef, level=None):
1728 i = 0
1729 arrs = []
1730 ntotal = 1
1731 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1732 if gdef and len(gdef) > i:
1733 sssn = gdef[i]
1734 else:
1735 sssn = (None,)*4
1737 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1738 ntotal *= len(arr)
1740 arrs.append(arr)
1741 i += 1
1743 arrs.append(self.coords[-1])
1744 return nditer_outer(arrs[:level])
1746 def make_sum_params(self, source, receiver, implementation='c',
1747 nthreads=0):
1748 assert implementation in ['c', 'python']
1750 out = []
1751 delays = source.times
1752 for comp, weights, icomponents in source.make_weights(
1753 receiver,
1754 self.component_scheme):
1756 weights *= self.factor
1758 args = self.make_indexing_args(source, receiver, icomponents)
1759 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1760 out.append((comp, args, delays_expanded, weights))
1762 return out
1764 def short_info(self):
1765 raise NotImplementedError('should be implemented in subclass')
1767 def get_shear_moduli(self, lat, lon, points,
1768 interpolation=None):
1769 '''
1770 Get shear moduli at given points from contained velocity model.
1772 :param lat: surface origin for coordinate system of ``points``
1773 :param points: NumPy array of shape ``(N, 3)``, where each row is
1774 a point ``(north, east, depth)``, relative to origin at
1775 ``(lat, lon)``
1776 :param interpolation: interpolation method. Choose from
1777 ``('nearest_neighbor', 'multilinear')``
1778 :returns: NumPy array of length N with extracted shear moduli at each
1779 point
1781 The default implementation retrieves and interpolates the shear moduli
1782 from the contained 1D velocity profile.
1783 '''
1784 return self.get_material_property(lat, lon, points,
1785 parameter='shear_moduli',
1786 interpolation=interpolation)
1788 def get_lambda_moduli(self, lat, lon, points,
1789 interpolation=None):
1790 '''
1791 Get lambda moduli at given points from contained velocity model.
1793 :param lat: surface origin for coordinate system of ``points``
1794 :param points: NumPy array of shape ``(N, 3)``, where each row is
1795 a point ``(north, east, depth)``, relative to origin at
1796 ``(lat, lon)``
1797 :param interpolation: interpolation method. Choose from
1798 ``('nearest_neighbor', 'multilinear')``
1799 :returns: NumPy array of length N with extracted shear moduli at each
1800 point
1802 The default implementation retrieves and interpolates the lambda moduli
1803 from the contained 1D velocity profile.
1804 '''
1805 return self.get_material_property(lat, lon, points,
1806 parameter='lambda_moduli',
1807 interpolation=interpolation)
1809 def get_bulk_moduli(self, lat, lon, points,
1810 interpolation=None):
1811 '''
1812 Get bulk moduli at given points from contained velocity model.
1814 :param lat: surface origin for coordinate system of ``points``
1815 :param points: NumPy array of shape ``(N, 3)``, where each row is
1816 a point ``(north, east, depth)``, relative to origin at
1817 ``(lat, lon)``
1818 :param interpolation: interpolation method. Choose from
1819 ``('nearest_neighbor', 'multilinear')``
1820 :returns: NumPy array of length N with extracted shear moduli at each
1821 point
1823 The default implementation retrieves and interpolates the lambda moduli
1824 from the contained 1D velocity profile.
1825 '''
1826 lambda_moduli = self.get_material_property(
1827 lat, lon, points, parameter='lambda_moduli',
1828 interpolation=interpolation)
1829 shear_moduli = self.get_material_property(
1830 lat, lon, points, parameter='shear_moduli',
1831 interpolation=interpolation)
1832 return lambda_moduli + (2 / 3) * shear_moduli
1834 def get_vs(self, lat, lon, points, interpolation=None):
1835 '''
1836 Get Vs at given points from contained velocity model.
1838 :param lat: surface origin for coordinate system of ``points``
1839 :param points: NumPy array of shape ``(N, 3)``, where each row is
1840 a point ``(north, east, depth)``, relative to origin at
1841 ``(lat, lon)``
1842 :param interpolation: interpolation method. Choose from
1843 ``('nearest_neighbor', 'multilinear')``
1844 :returns: NumPy array of length N with extracted shear moduli at each
1845 point
1847 The default implementation retrieves and interpolates Vs
1848 from the contained 1D velocity profile.
1849 '''
1850 return self.get_material_property(lat, lon, points,
1851 parameter='vs',
1852 interpolation=interpolation)
1854 def get_vp(self, lat, lon, points, interpolation=None):
1855 '''
1856 Get Vp at given points from contained velocity model.
1858 :param lat: surface origin for coordinate system of ``points``
1859 :param points: NumPy array of shape ``(N, 3)``, where each row is
1860 a point ``(north, east, depth)``, relative to origin at
1861 ``(lat, lon)``
1862 :param interpolation: interpolation method. Choose from
1863 ``('nearest_neighbor', 'multilinear')``
1864 :returns: NumPy array of length N with extracted shear moduli at each
1865 point
1867 The default implementation retrieves and interpolates Vp
1868 from the contained 1D velocity profile.
1869 '''
1870 return self.get_material_property(lat, lon, points,
1871 parameter='vp',
1872 interpolation=interpolation)
1874 def get_rho(self, lat, lon, points, interpolation=None):
1875 '''
1876 Get rho at given points from contained velocity model.
1878 :param lat: surface origin for coordinate system of ``points``
1879 :param points: NumPy array of shape ``(N, 3)``, where each row is
1880 a point ``(north, east, depth)``, relative to origin at
1881 ``(lat, lon)``
1882 :param interpolation: interpolation method. Choose from
1883 ``('nearest_neighbor', 'multilinear')``
1884 :returns: NumPy array of length N with extracted shear moduli at each
1885 point
1887 The default implementation retrieves and interpolates rho
1888 from the contained 1D velocity profile.
1889 '''
1890 return self.get_material_property(lat, lon, points,
1891 parameter='rho',
1892 interpolation=interpolation)
1894 def get_material_property(self, lat, lon, points, parameter='vs',
1895 interpolation=None):
1897 if interpolation is None:
1898 raise TypeError('Interpolation method not defined! available: '
1899 'multilinear', 'nearest_neighbor')
1901 earthmod = self.earthmodel_1d
1902 store_depth_profile = self.get_source_depths()
1903 z_profile = earthmod.profile('z')
1905 if parameter == 'vs':
1906 vs_profile = earthmod.profile('vs')
1907 profile = num.interp(
1908 store_depth_profile, z_profile, vs_profile)
1910 elif parameter == 'vp':
1911 vp_profile = earthmod.profile('vp')
1912 profile = num.interp(
1913 store_depth_profile, z_profile, vp_profile)
1915 elif parameter == 'rho':
1916 rho_profile = earthmod.profile('rho')
1918 profile = num.interp(
1919 store_depth_profile, z_profile, rho_profile)
1921 elif parameter == 'shear_moduli':
1922 vs_profile = earthmod.profile('vs')
1923 rho_profile = earthmod.profile('rho')
1925 store_vs_profile = num.interp(
1926 store_depth_profile, z_profile, vs_profile)
1927 store_rho_profile = num.interp(
1928 store_depth_profile, z_profile, rho_profile)
1930 profile = num.power(store_vs_profile, 2) * store_rho_profile
1932 elif parameter == 'lambda_moduli':
1933 vs_profile = earthmod.profile('vs')
1934 vp_profile = earthmod.profile('vp')
1935 rho_profile = earthmod.profile('rho')
1937 store_vs_profile = num.interp(
1938 store_depth_profile, z_profile, vs_profile)
1939 store_vp_profile = num.interp(
1940 store_depth_profile, z_profile, vp_profile)
1941 store_rho_profile = num.interp(
1942 store_depth_profile, z_profile, rho_profile)
1944 profile = store_rho_profile * (
1945 num.power(store_vp_profile, 2) -
1946 num.power(store_vs_profile, 2) * 2)
1947 else:
1948 raise TypeError(
1949 'parameter %s not available' % parameter)
1951 if interpolation == 'multilinear':
1952 kind = 'linear'
1953 elif interpolation == 'nearest_neighbor':
1954 kind = 'nearest'
1955 else:
1956 raise TypeError(
1957 'Interpolation method %s not available' % interpolation)
1959 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1961 try:
1962 return interpolator(points[:, 2])
1963 except ValueError:
1964 raise OutOfBounds()
1966 def is_static(self):
1967 for code in ('psgrn_pscmp', 'poel'):
1968 if self.modelling_code_id.startswith(code):
1969 return True
1970 return False
1972 def is_dynamic(self):
1973 return not self.is_static()
1975 def get_source_depths(self):
1976 raise NotImplementedError('must be implemented in subclass')
1978 def get_tabulated_phase(self, phase_id):
1979 '''
1980 Get tabulated phase definition.
1981 '''
1983 for pdef in self.tabulated_phases:
1984 if pdef.id == phase_id:
1985 return pdef
1987 raise StoreError('No such phase: %s' % phase_id)
1989 def fix_ttt_holes(self, sptree, mode):
1990 raise StoreError(
1991 'Cannot fix travel time table holes in GF stores of type %s.'
1992 % self.short_type)
1995class ConfigTypeA(Config):
1996 '''
1997 Cylindrical symmetry, 1D earth model, single receiver depth
1999 * Problem is invariant to horizontal translations and rotations around
2000 vertical axis.
2002 * All receivers must be at the same depth (e.g. at the surface)
2003 High level index variables: ``(source_depth, distance,
2004 component)``
2006 * The ``distance`` is the surface distance between source and receiver
2007 points.
2008 '''
2010 receiver_depth = Float.T(
2011 default=0.0,
2012 help='Fixed receiver depth [m].')
2014 source_depth_min = Float.T(
2015 help='Minimum source depth [m].')
2017 source_depth_max = Float.T(
2018 help='Maximum source depth [m].')
2020 source_depth_delta = Float.T(
2021 help='Grid spacing of source depths [m]')
2023 distance_min = Float.T(
2024 help='Minimum source-receiver surface distance [m].')
2026 distance_max = Float.T(
2027 help='Maximum source-receiver surface distance [m].')
2029 distance_delta = Float.T(
2030 help='Grid spacing of source-receiver surface distance [m].')
2032 short_type = 'A'
2034 provided_schemes = [
2035 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2037 def get_surface_distance(self, args):
2038 return args[1]
2040 def get_distance(self, args):
2041 return math.sqrt(args[0]**2 + args[1]**2)
2043 def get_source_depth(self, args):
2044 return args[0]
2046 def get_source_depths(self):
2047 return self.coords[0]
2049 def get_receiver_depth(self, args):
2050 return self.receiver_depth
2052 def _update(self):
2053 self.mins = num.array(
2054 [self.source_depth_min, self.distance_min], dtype=float)
2055 self.maxs = num.array(
2056 [self.source_depth_max, self.distance_max], dtype=float)
2057 self.deltas = num.array(
2058 [self.source_depth_delta, self.distance_delta],
2059 dtype=float)
2060 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2061 vicinity_eps).astype(int) + 1
2062 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2063 self.deltat = 1.0/self.sample_rate
2064 self.nrecords = num.product(self.ns) * self.ncomponents
2065 self.coords = tuple(num.linspace(mi, ma, n) for
2066 (mi, ma, n) in
2067 zip(self.mins, self.effective_maxs, self.ns)) + \
2068 (num.arange(self.ncomponents),)
2070 self.nsource_depths, self.ndistances = self.ns
2072 def _make_index_functions(self):
2074 amin, bmin = self.mins
2075 da, db = self.deltas
2076 na, nb = self.ns
2078 ng = self.ncomponents
2080 def index_function(a, b, ig):
2081 ia = int(round((a - amin) / da))
2082 ib = int(round((b - bmin) / db))
2083 try:
2084 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2085 except ValueError:
2086 raise OutOfBounds()
2088 def indices_function(a, b, ig):
2089 ia = num.round((a - amin) / da).astype(int)
2090 ib = num.round((b - bmin) / db).astype(int)
2091 try:
2092 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2093 except ValueError:
2094 for ia_, ib_, ig_ in zip(ia, ib, ig):
2095 try:
2096 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2097 except ValueError:
2098 raise OutOfBounds()
2100 def grid_interpolation_coefficients(a, b):
2101 ias = indi12((a - amin) / da, na)
2102 ibs = indi12((b - bmin) / db, nb)
2103 return ias, ibs
2105 def vicinity_function(a, b, ig):
2106 ias, ibs = grid_interpolation_coefficients(a, b)
2108 if not (0 <= ig < ng):
2109 raise OutOfBounds()
2111 indis = []
2112 weights = []
2113 for ia, va in ias:
2114 iia = ia*nb*ng
2115 for ib, vb in ibs:
2116 indis.append(iia + ib*ng + ig)
2117 weights.append(va*vb)
2119 return num.array(indis), num.array(weights)
2121 def vicinities_function(a, b, ig):
2123 xa = (a - amin) / da
2124 xb = (b - bmin) / db
2126 xa_fl = num.floor(xa)
2127 xa_ce = num.ceil(xa)
2128 xb_fl = num.floor(xb)
2129 xb_ce = num.ceil(xb)
2130 va_fl = 1.0 - (xa - xa_fl)
2131 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2132 vb_fl = 1.0 - (xb - xb_fl)
2133 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2135 ia_fl = xa_fl.astype(int)
2136 ia_ce = xa_ce.astype(int)
2137 ib_fl = xb_fl.astype(int)
2138 ib_ce = xb_ce.astype(int)
2140 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2141 raise OutOfBounds()
2143 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2144 raise OutOfBounds()
2146 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2147 raise OutOfBounds()
2149 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2150 raise OutOfBounds()
2152 irecords = num.empty(a.size*4, dtype=int)
2153 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2154 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2155 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2156 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2158 weights = num.empty(a.size*4, dtype=float)
2159 weights[0::4] = va_fl * vb_fl
2160 weights[1::4] = va_ce * vb_fl
2161 weights[2::4] = va_fl * vb_ce
2162 weights[3::4] = va_ce * vb_ce
2164 return irecords, weights
2166 self._index_function = index_function
2167 self._indices_function = indices_function
2168 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2169 self._vicinity_function = vicinity_function
2170 self._vicinities_function = vicinities_function
2172 def make_indexing_args(self, source, receiver, icomponents):
2173 nc = icomponents.size
2174 dists = source.distances_to(receiver)
2175 n = dists.size
2176 return (num.tile(source.depths, nc//n),
2177 num.tile(dists, nc//n),
2178 icomponents)
2180 def make_indexing_args1(self, source, receiver):
2181 return (source.depth, source.distance_to(receiver))
2183 @property
2184 def short_extent(self):
2185 return '%g:%g:%g x %g:%g:%g' % (
2186 self.source_depth_min/km,
2187 self.source_depth_max/km,
2188 self.source_depth_delta/km,
2189 self.distance_min/km,
2190 self.distance_max/km,
2191 self.distance_delta/km)
2193 def fix_ttt_holes(self, sptree, mode):
2194 from pyrocko import eikonal_ext, spit
2196 nodes = self.nodes(level=-1)
2198 delta = self.deltas[-1]
2199 assert num.all(delta == self.deltas)
2201 nsources, ndistances = self.ns
2203 points = num.zeros((nodes.shape[0], 3))
2204 points[:, 0] = nodes[:, 1]
2205 points[:, 2] = nodes[:, 0]
2207 speeds = self.get_material_property(
2208 0., 0., points,
2209 parameter='vp' if mode == cake.P else 'vs',
2210 interpolation='multilinear')
2212 speeds = speeds.reshape((nsources, ndistances))
2214 times = sptree.interpolate_many(nodes)
2216 times[num.isnan(times)] = -1.
2217 times = times.reshape(speeds.shape)
2219 try:
2220 eikonal_ext.eikonal_solver_fmm_cartesian(
2221 speeds, times, delta)
2222 except eikonal_ext.EikonalExtError as e:
2223 if str(e).endswith('please check results'):
2224 logger.debug(
2225 'Got a warning from eikonal solver '
2226 '- may be ok...')
2227 else:
2228 raise
2230 def func(x):
2231 ibs, ics = \
2232 self.grid_interpolation_coefficients(*x)
2234 t = 0
2235 for ib, vb in ibs:
2236 for ic, vc in ics:
2237 t += times[ib, ic] * vb * vc
2239 return t
2241 return spit.SPTree(
2242 f=func,
2243 ftol=sptree.ftol,
2244 xbounds=sptree.xbounds,
2245 xtols=sptree.xtols)
2248class ConfigTypeB(Config):
2249 '''
2250 Cylindrical symmetry, 1D earth model, variable receiver depth
2252 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2253 receiver depth
2255 * High level index variables: ``(receiver_depth, source_depth,
2256 receiver_distance, component)``
2257 '''
2259 receiver_depth_min = Float.T(
2260 help='Minimum receiver depth [m].')
2262 receiver_depth_max = Float.T(
2263 help='Maximum receiver depth [m].')
2265 receiver_depth_delta = Float.T(
2266 help='Grid spacing of receiver depths [m]')
2268 source_depth_min = Float.T(
2269 help='Minimum source depth [m].')
2271 source_depth_max = Float.T(
2272 help='Maximum source depth [m].')
2274 source_depth_delta = Float.T(
2275 help='Grid spacing of source depths [m]')
2277 distance_min = Float.T(
2278 help='Minimum source-receiver surface distance [m].')
2280 distance_max = Float.T(
2281 help='Maximum source-receiver surface distance [m].')
2283 distance_delta = Float.T(
2284 help='Grid spacing of source-receiver surface distances [m].')
2286 short_type = 'B'
2288 provided_schemes = [
2289 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2291 def get_distance(self, args):
2292 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2294 def get_surface_distance(self, args):
2295 return args[2]
2297 def get_source_depth(self, args):
2298 return args[1]
2300 def get_receiver_depth(self, args):
2301 return args[0]
2303 def get_source_depths(self):
2304 return self.coords[1]
2306 def _update(self):
2307 self.mins = num.array([
2308 self.receiver_depth_min,
2309 self.source_depth_min,
2310 self.distance_min],
2311 dtype=float)
2313 self.maxs = num.array([
2314 self.receiver_depth_max,
2315 self.source_depth_max,
2316 self.distance_max],
2317 dtype=float)
2319 self.deltas = num.array([
2320 self.receiver_depth_delta,
2321 self.source_depth_delta,
2322 self.distance_delta],
2323 dtype=float)
2325 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2326 vicinity_eps).astype(int) + 1
2327 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2328 self.deltat = 1.0/self.sample_rate
2329 self.nrecords = num.product(self.ns) * self.ncomponents
2330 self.coords = tuple(num.linspace(mi, ma, n) for
2331 (mi, ma, n) in
2332 zip(self.mins, self.effective_maxs, self.ns)) + \
2333 (num.arange(self.ncomponents),)
2334 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2336 def _make_index_functions(self):
2338 amin, bmin, cmin = self.mins
2339 da, db, dc = self.deltas
2340 na, nb, nc = self.ns
2341 ng = self.ncomponents
2343 def index_function(a, b, c, ig):
2344 ia = int(round((a - amin) / da))
2345 ib = int(round((b - bmin) / db))
2346 ic = int(round((c - cmin) / dc))
2347 try:
2348 return num.ravel_multi_index((ia, ib, ic, ig),
2349 (na, nb, nc, ng))
2350 except ValueError:
2351 raise OutOfBounds()
2353 def indices_function(a, b, c, ig):
2354 ia = num.round((a - amin) / da).astype(int)
2355 ib = num.round((b - bmin) / db).astype(int)
2356 ic = num.round((c - cmin) / dc).astype(int)
2357 try:
2358 return num.ravel_multi_index((ia, ib, ic, ig),
2359 (na, nb, nc, ng))
2360 except ValueError:
2361 raise OutOfBounds()
2363 def grid_interpolation_coefficients(a, b, c):
2364 ias = indi12((a - amin) / da, na)
2365 ibs = indi12((b - bmin) / db, nb)
2366 ics = indi12((c - cmin) / dc, nc)
2367 return ias, ibs, ics
2369 def vicinity_function(a, b, c, ig):
2370 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2372 if not (0 <= ig < ng):
2373 raise OutOfBounds()
2375 indis = []
2376 weights = []
2377 for ia, va in ias:
2378 iia = ia*nb*nc*ng
2379 for ib, vb in ibs:
2380 iib = ib*nc*ng
2381 for ic, vc in ics:
2382 indis.append(iia + iib + ic*ng + ig)
2383 weights.append(va*vb*vc)
2385 return num.array(indis), num.array(weights)
2387 def vicinities_function(a, b, c, ig):
2389 xa = (a - amin) / da
2390 xb = (b - bmin) / db
2391 xc = (c - cmin) / dc
2393 xa_fl = num.floor(xa)
2394 xa_ce = num.ceil(xa)
2395 xb_fl = num.floor(xb)
2396 xb_ce = num.ceil(xb)
2397 xc_fl = num.floor(xc)
2398 xc_ce = num.ceil(xc)
2399 va_fl = 1.0 - (xa - xa_fl)
2400 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2401 vb_fl = 1.0 - (xb - xb_fl)
2402 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2403 vc_fl = 1.0 - (xc - xc_fl)
2404 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2406 ia_fl = xa_fl.astype(int)
2407 ia_ce = xa_ce.astype(int)
2408 ib_fl = xb_fl.astype(int)
2409 ib_ce = xb_ce.astype(int)
2410 ic_fl = xc_fl.astype(int)
2411 ic_ce = xc_ce.astype(int)
2413 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2414 raise OutOfBounds()
2416 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2417 raise OutOfBounds()
2419 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2420 raise OutOfBounds()
2422 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2423 raise OutOfBounds()
2425 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2426 raise OutOfBounds()
2428 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2429 raise OutOfBounds()
2431 irecords = num.empty(a.size*8, dtype=int)
2432 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2433 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2434 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2435 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2436 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2437 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2438 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2439 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2441 weights = num.empty(a.size*8, dtype=float)
2442 weights[0::8] = va_fl * vb_fl * vc_fl
2443 weights[1::8] = va_ce * vb_fl * vc_fl
2444 weights[2::8] = va_fl * vb_ce * vc_fl
2445 weights[3::8] = va_ce * vb_ce * vc_fl
2446 weights[4::8] = va_fl * vb_fl * vc_ce
2447 weights[5::8] = va_ce * vb_fl * vc_ce
2448 weights[6::8] = va_fl * vb_ce * vc_ce
2449 weights[7::8] = va_ce * vb_ce * vc_ce
2451 return irecords, weights
2453 self._index_function = index_function
2454 self._indices_function = indices_function
2455 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2456 self._vicinity_function = vicinity_function
2457 self._vicinities_function = vicinities_function
2459 def make_indexing_args(self, source, receiver, icomponents):
2460 nc = icomponents.size
2461 dists = source.distances_to(receiver)
2462 n = dists.size
2463 receiver_depths = num.empty(nc)
2464 receiver_depths.fill(receiver.depth)
2465 return (receiver_depths,
2466 num.tile(source.depths, nc//n),
2467 num.tile(dists, nc//n),
2468 icomponents)
2470 def make_indexing_args1(self, source, receiver):
2471 return (receiver.depth,
2472 source.depth,
2473 source.distance_to(receiver))
2475 @property
2476 def short_extent(self):
2477 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2478 self.receiver_depth_min/km,
2479 self.receiver_depth_max/km,
2480 self.receiver_depth_delta/km,
2481 self.source_depth_min/km,
2482 self.source_depth_max/km,
2483 self.source_depth_delta/km,
2484 self.distance_min/km,
2485 self.distance_max/km,
2486 self.distance_delta/km)
2488 def fix_ttt_holes(self, sptree, mode):
2489 from pyrocko import eikonal_ext, spit
2491 nodes_sr = self.nodes(minlevel=1, level=-1)
2493 delta = self.deltas[-1]
2494 assert num.all(delta == self.deltas[1:])
2496 nreceivers, nsources, ndistances = self.ns
2498 points = num.zeros((nodes_sr.shape[0], 3))
2499 points[:, 0] = nodes_sr[:, 1]
2500 points[:, 2] = nodes_sr[:, 0]
2502 speeds = self.get_material_property(
2503 0., 0., points,
2504 parameter='vp' if mode == cake.P else 'vs',
2505 interpolation='multilinear')
2507 speeds = speeds.reshape((nsources, ndistances))
2509 receiver_times = []
2510 for ireceiver in range(nreceivers):
2511 nodes = num.hstack([
2512 num_full(
2513 (nodes_sr.shape[0], 1),
2514 self.coords[0][ireceiver]),
2515 nodes_sr])
2517 times = sptree.interpolate_many(nodes)
2519 times[num.isnan(times)] = -1.
2521 times = times.reshape(speeds.shape)
2523 try:
2524 eikonal_ext.eikonal_solver_fmm_cartesian(
2525 speeds, times, delta)
2526 except eikonal_ext.EikonalExtError as e:
2527 if str(e).endswith('please check results'):
2528 logger.debug(
2529 'Got a warning from eikonal solver '
2530 '- may be ok...')
2531 else:
2532 raise
2534 receiver_times.append(times)
2536 def func(x):
2537 ias, ibs, ics = \
2538 self.grid_interpolation_coefficients(*x)
2540 t = 0
2541 for ia, va in ias:
2542 times = receiver_times[ia]
2543 for ib, vb in ibs:
2544 for ic, vc in ics:
2545 t += times[ib, ic] * va * vb * vc
2547 return t
2549 return spit.SPTree(
2550 f=func,
2551 ftol=sptree.ftol,
2552 xbounds=sptree.xbounds,
2553 xtols=sptree.xtols)
2556class ConfigTypeC(Config):
2557 '''
2558 No symmetrical constraints, one fixed receiver position.
2560 * Cartesian 3D source volume around a reference point
2562 * High level index variables: ``(source_depth,
2563 source_east_shift, source_north_shift, component)``
2564 '''
2566 receiver = Receiver.T(
2567 help='Receiver location')
2569 source_origin = Location.T(
2570 help='Origin of the source volume grid.')
2572 source_depth_min = Float.T(
2573 help='Minimum source depth [m].')
2575 source_depth_max = Float.T(
2576 help='Maximum source depth [m].')
2578 source_depth_delta = Float.T(
2579 help='Source depth grid spacing [m].')
2581 source_east_shift_min = Float.T(
2582 help='Minimum easting of source grid [m].')
2584 source_east_shift_max = Float.T(
2585 help='Maximum easting of source grid [m].')
2587 source_east_shift_delta = Float.T(
2588 help='Source volume grid spacing in east direction [m].')
2590 source_north_shift_min = Float.T(
2591 help='Minimum northing of source grid [m].')
2593 source_north_shift_max = Float.T(
2594 help='Maximum northing of source grid [m].')
2596 source_north_shift_delta = Float.T(
2597 help='Source volume grid spacing in north direction [m].')
2599 short_type = 'C'
2601 provided_schemes = ['elastic18']
2603 def get_surface_distance(self, args):
2604 _, source_east_shift, source_north_shift, _ = args
2605 sorig = self.source_origin
2606 sloc = Location(
2607 lat=sorig.lat,
2608 lon=sorig.lon,
2609 north_shift=sorig.north_shift + source_north_shift,
2610 east_shift=sorig.east_shift + source_east_shift)
2612 return self.receiver.distance_to(sloc)
2614 def get_distance(self, args):
2615 sdepth, source_east_shift, source_north_shift, _ = args
2616 sorig = self.source_origin
2617 sloc = Location(
2618 lat=sorig.lat,
2619 lon=sorig.lon,
2620 north_shift=sorig.north_shift + source_north_shift,
2621 east_shift=sorig.east_shift + source_east_shift)
2623 return self.receiver.distance_3d_to(sloc)
2625 def get_source_depth(self, args):
2626 return args[0]
2628 def get_receiver_depth(self, args):
2629 return self.receiver.depth
2631 def get_source_depths(self):
2632 return self.coords[0]
2634 def _update(self):
2635 self.mins = num.array([
2636 self.source_depth_min,
2637 self.source_east_shift_min,
2638 self.source_north_shift_min],
2639 dtype=float)
2641 self.maxs = num.array([
2642 self.source_depth_max,
2643 self.source_east_shift_max,
2644 self.source_north_shift_max],
2645 dtype=float)
2647 self.deltas = num.array([
2648 self.source_depth_delta,
2649 self.source_east_shift_delta,
2650 self.source_north_shift_delta],
2651 dtype=float)
2653 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2654 vicinity_eps).astype(int) + 1
2655 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2656 self.deltat = 1.0/self.sample_rate
2657 self.nrecords = num.product(self.ns) * self.ncomponents
2659 self.coords = tuple(
2660 num.linspace(mi, ma, n) for (mi, ma, n) in
2661 zip(self.mins, self.effective_maxs, self.ns)) + \
2662 (num.arange(self.ncomponents),)
2664 def _make_index_functions(self):
2666 amin, bmin, cmin = self.mins
2667 da, db, dc = self.deltas
2668 na, nb, nc = self.ns
2669 ng = self.ncomponents
2671 def index_function(a, b, c, ig):
2672 ia = int(round((a - amin) / da))
2673 ib = int(round((b - bmin) / db))
2674 ic = int(round((c - cmin) / dc))
2675 try:
2676 return num.ravel_multi_index((ia, ib, ic, ig),
2677 (na, nb, nc, ng))
2678 except ValueError:
2679 raise OutOfBounds(values=(a, b, c, ig))
2681 def indices_function(a, b, c, ig):
2682 ia = num.round((a - amin) / da).astype(int)
2683 ib = num.round((b - bmin) / db).astype(int)
2684 ic = num.round((c - cmin) / dc).astype(int)
2686 try:
2687 return num.ravel_multi_index((ia, ib, ic, ig),
2688 (na, nb, nc, ng))
2689 except ValueError:
2690 raise OutOfBounds()
2692 def vicinity_function(a, b, c, ig):
2693 ias = indi12((a - amin) / da, na)
2694 ibs = indi12((b - bmin) / db, nb)
2695 ics = indi12((c - cmin) / dc, nc)
2697 if not (0 <= ig < ng):
2698 raise OutOfBounds()
2700 indis = []
2701 weights = []
2702 for ia, va in ias:
2703 iia = ia*nb*nc*ng
2704 for ib, vb in ibs:
2705 iib = ib*nc*ng
2706 for ic, vc in ics:
2707 indis.append(iia + iib + ic*ng + ig)
2708 weights.append(va*vb*vc)
2710 return num.array(indis), num.array(weights)
2712 def vicinities_function(a, b, c, ig):
2714 xa = (a-amin) / da
2715 xb = (b-bmin) / db
2716 xc = (c-cmin) / dc
2718 xa_fl = num.floor(xa)
2719 xa_ce = num.ceil(xa)
2720 xb_fl = num.floor(xb)
2721 xb_ce = num.ceil(xb)
2722 xc_fl = num.floor(xc)
2723 xc_ce = num.ceil(xc)
2724 va_fl = 1.0 - (xa - xa_fl)
2725 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2726 vb_fl = 1.0 - (xb - xb_fl)
2727 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2728 vc_fl = 1.0 - (xc - xc_fl)
2729 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2731 ia_fl = xa_fl.astype(int)
2732 ia_ce = xa_ce.astype(int)
2733 ib_fl = xb_fl.astype(int)
2734 ib_ce = xb_ce.astype(int)
2735 ic_fl = xc_fl.astype(int)
2736 ic_ce = xc_ce.astype(int)
2738 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2739 raise OutOfBounds()
2741 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2742 raise OutOfBounds()
2744 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2745 raise OutOfBounds()
2747 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2748 raise OutOfBounds()
2750 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2751 raise OutOfBounds()
2753 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2754 raise OutOfBounds()
2756 irecords = num.empty(a.size*8, dtype=int)
2757 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2758 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2759 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2760 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2761 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2762 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2763 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2764 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2766 weights = num.empty(a.size*8, dtype=float)
2767 weights[0::8] = va_fl * vb_fl * vc_fl
2768 weights[1::8] = va_ce * vb_fl * vc_fl
2769 weights[2::8] = va_fl * vb_ce * vc_fl
2770 weights[3::8] = va_ce * vb_ce * vc_fl
2771 weights[4::8] = va_fl * vb_fl * vc_ce
2772 weights[5::8] = va_ce * vb_fl * vc_ce
2773 weights[6::8] = va_fl * vb_ce * vc_ce
2774 weights[7::8] = va_ce * vb_ce * vc_ce
2776 return irecords, weights
2778 self._index_function = index_function
2779 self._indices_function = indices_function
2780 self._vicinity_function = vicinity_function
2781 self._vicinities_function = vicinities_function
2783 def make_indexing_args(self, source, receiver, icomponents):
2784 nc = icomponents.size
2786 dists = source.distances_to(self.source_origin)
2787 azis, _ = source.azibazis_to(self.source_origin)
2789 source_north_shifts = - num.cos(d2r*azis) * dists
2790 source_east_shifts = - num.sin(d2r*azis) * dists
2791 source_depths = source.depths - self.source_origin.depth
2793 n = dists.size
2795 return (num.tile(source_depths, nc//n),
2796 num.tile(source_east_shifts, nc//n),
2797 num.tile(source_north_shifts, nc//n),
2798 icomponents)
2800 def make_indexing_args1(self, source, receiver):
2801 dist = source.distance_to(self.source_origin)
2802 azi, _ = source.azibazi_to(self.source_origin)
2804 source_north_shift = - num.cos(d2r*azi) * dist
2805 source_east_shift = - num.sin(d2r*azi) * dist
2806 source_depth = source.depth - self.source_origin.depth
2808 return (source_depth,
2809 source_east_shift,
2810 source_north_shift)
2812 @property
2813 def short_extent(self):
2814 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2815 self.source_depth_min/km,
2816 self.source_depth_max/km,
2817 self.source_depth_delta/km,
2818 self.source_east_shift_min/km,
2819 self.source_east_shift_max/km,
2820 self.source_east_shift_delta/km,
2821 self.source_north_shift_min/km,
2822 self.source_north_shift_max/km,
2823 self.source_north_shift_delta/km)
2826class Weighting(Object):
2827 factor = Float.T(default=1.0)
2830class Taper(Object):
2831 tmin = Timing.T()
2832 tmax = Timing.T()
2833 tfade = Float.T(default=0.0)
2834 shape = StringChoice.T(
2835 choices=['cos', 'linear'],
2836 default='cos',
2837 optional=True)
2840class SimplePattern(SObject):
2842 _pool = {}
2844 def __init__(self, pattern):
2845 self._pattern = pattern
2846 SObject.__init__(self)
2848 def __str__(self):
2849 return self._pattern
2851 @property
2852 def regex(self):
2853 pool = SimplePattern._pool
2854 if self.pattern not in pool:
2855 rpat = '|'.join(fnmatch.translate(x) for
2856 x in self.pattern.split('|'))
2857 pool[self.pattern] = re.compile(rpat, re.I)
2859 return pool[self.pattern]
2861 def match(self, s):
2862 return self.regex.match(s)
2865class WaveformType(StringChoice):
2866 choices = ['dis', 'vel', 'acc',
2867 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2868 'envelope_dis', 'envelope_vel', 'envelope_acc']
2871class ChannelSelection(Object):
2872 pattern = SimplePattern.T(optional=True)
2873 min_sample_rate = Float.T(optional=True)
2874 max_sample_rate = Float.T(optional=True)
2877class StationSelection(Object):
2878 includes = SimplePattern.T()
2879 excludes = SimplePattern.T()
2880 distance_min = Float.T(optional=True)
2881 distance_max = Float.T(optional=True)
2882 azimuth_min = Float.T(optional=True)
2883 azimuth_max = Float.T(optional=True)
2886class WaveformSelection(Object):
2887 channel_selection = ChannelSelection.T(optional=True)
2888 station_selection = StationSelection.T(optional=True)
2889 taper = Taper.T()
2890 # filter = FrequencyResponse.T()
2891 waveform_type = WaveformType.T(default='dis')
2892 weighting = Weighting.T(optional=True)
2893 sample_rate = Float.T(optional=True)
2894 gf_store_id = StringID.T(optional=True)
2897def indi12(x, n):
2898 '''
2899 Get linear interpolation index and weight.
2900 '''
2902 r = round(x)
2903 if abs(r - x) < vicinity_eps:
2904 i = int(r)
2905 if not (0 <= i < n):
2906 raise OutOfBounds()
2908 return ((int(r), 1.),)
2909 else:
2910 f = math.floor(x)
2911 i = int(f)
2912 if not (0 <= i < n-1):
2913 raise OutOfBounds()
2915 v = x-f
2916 return ((i, 1.-v), (i + 1, v))
2919def float_or_none(s):
2920 units = {
2921 'k': 1e3,
2922 'M': 1e6,
2923 }
2925 factor = 1.0
2926 if s and s[-1] in units:
2927 factor = units[s[-1]]
2928 s = s[:-1]
2929 if not s:
2930 raise ValueError("unit without a number: '%s'" % s)
2932 if s:
2933 return float(s) * factor
2934 else:
2935 return None
2938class GridSpecError(Exception):
2939 def __init__(self, s):
2940 Exception.__init__(self, 'invalid grid specification: %s' % s)
2943def parse_grid_spec(spec):
2944 try:
2945 result = []
2946 for dspec in spec.split(','):
2947 t = dspec.split('@')
2948 num = start = stop = step = None
2949 if len(t) == 2:
2950 num = int(t[1])
2951 if num <= 0:
2952 raise GridSpecError(spec)
2954 elif len(t) > 2:
2955 raise GridSpecError(spec)
2957 s = t[0]
2958 v = [float_or_none(x) for x in s.split(':')]
2959 if len(v) == 1:
2960 start = stop = v[0]
2961 if len(v) >= 2:
2962 start, stop = v[0:2]
2963 if len(v) == 3:
2964 step = v[2]
2966 if len(v) > 3 or (len(v) > 2 and num is not None):
2967 raise GridSpecError(spec)
2969 if step == 0.0:
2970 raise GridSpecError(spec)
2972 result.append((start, stop, step, num))
2974 except ValueError:
2975 raise GridSpecError(spec)
2977 return result
2980def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2981 swap = step is not None and step < 0.
2982 if start is None:
2983 start = ma if swap else mi
2984 if stop is None:
2985 stop = mi if swap else ma
2986 if step is None:
2987 step = -inc if ma < mi else inc
2988 if num is None:
2989 if (step < 0) != (stop-start < 0):
2990 raise GridSpecError()
2992 num = int(round((stop-start)/step))+1
2993 stop2 = start + (num-1)*step
2994 if abs(stop-stop2) > eps:
2995 num = int(math.floor((stop-start)/step))+1
2996 stop = start + (num-1)*step
2997 else:
2998 stop = stop2
3000 if start == stop:
3001 num = 1
3003 return start, stop, num
3006def nditer_outer(x):
3007 return num.nditer(
3008 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
3011def nodes(xs):
3012 ns = [x.size for x in xs]
3013 nnodes = num.prod(ns)
3014 ndim = len(xs)
3015 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
3016 for idim in range(ndim-1, -1, -1):
3017 x = xs[idim]
3018 nrepeat = num.prod(ns[idim+1:], dtype=int)
3019 ntile = num.prod(ns[:idim], dtype=int)
3020 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
3022 return nodes
3025def filledi(x, n):
3026 a = num.empty(n, dtype=int)
3027 a.fill(x)
3028 return a
3031config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3033discretized_source_classes = [
3034 DiscretizedExplosionSource,
3035 DiscretizedSFSource,
3036 DiscretizedMTSource,
3037 DiscretizedPorePressureSource]
3040__all__ = '''
3041Earthmodel1D
3042StringID
3043ScopeType
3044WaveformType
3045QuantityType
3046NearfieldTermsType
3047Reference
3048Region
3049CircularRegion
3050RectangularRegion
3051PhaseSelect
3052InvalidTimingSpecification
3053Timing
3054TPDef
3055OutOfBounds
3056Location
3057Receiver
3058'''.split() + [
3059 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3060ComponentScheme
3061component_scheme_to_description
3062component_schemes
3063Config
3064GridSpecError
3065Weighting
3066Taper
3067SimplePattern
3068WaveformType
3069ChannelSelection
3070StationSelection
3071WaveformSelection
3072nditer_outer
3073dump
3074load
3075discretized_source_classes
3076config_type_classes
3077UnavailableScheme
3078InterpolationMethod
3079SeismosizerTrace
3080SeismosizerResult
3081Result
3082StaticResult
3083TimingAttributeError
3084'''.split()