1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import math
8import re
9import fnmatch
10import logging
12import numpy as num
13from scipy.interpolate import interp1d
15from pyrocko.guts import (Object, SObject, String, StringChoice,
16 StringPattern, Unicode, Float, Bool, Int, TBase,
17 List, ValidationError, Timestamp, Tuple, Dict)
18from pyrocko.guts import dump, load # noqa
19from pyrocko.guts_array import literal, Array
20from pyrocko.model import Location, gnss
22from pyrocko import cake, orthodrome, spit, moment_tensor, trace
23from pyrocko.util import num_full
25from .error import StoreError
27try:
28 new_str = unicode
29except NameError:
30 new_str = str
32guts_prefix = 'pf'
34d2r = math.pi / 180.
35r2d = 1.0 / d2r
36km = 1000.
37vicinity_eps = 1e-5
39logger = logging.getLogger('pyrocko.gf.meta')
42def fmt_choices(cls):
43 return 'choices: %s' % ', '.join(
44 "``'%s'``" % choice for choice in cls.choices)
47class InterpolationMethod(StringChoice):
48 choices = ['nearest_neighbor', 'multilinear']
51class SeismosizerTrace(Object):
53 codes = Tuple.T(
54 4, String.T(),
55 default=('', 'STA', '', 'Z'),
56 help='network, station, location and channel codes')
58 data = Array.T(
59 shape=(None,),
60 dtype=num.float32,
61 serialize_as='base64',
62 serialize_dtype=num.dtype('<f4'),
63 help='numpy array with data samples')
65 deltat = Float.T(
66 default=1.0,
67 help='sampling interval [s]')
69 tmin = Timestamp.T(
70 default=Timestamp.D('1970-01-01 00:00:00'),
71 help='time of first sample as a system timestamp [s]')
73 def pyrocko_trace(self):
74 c = self.codes
75 return trace.Trace(c[0], c[1], c[2], c[3],
76 ydata=self.data,
77 deltat=self.deltat,
78 tmin=self.tmin)
80 @classmethod
81 def from_pyrocko_trace(cls, tr, **kwargs):
82 d = dict(
83 codes=tr.nslc_id,
84 tmin=tr.tmin,
85 deltat=tr.deltat,
86 data=num.asarray(tr.get_ydata(), dtype=num.float32))
87 d.update(kwargs)
88 return cls(**d)
91class SeismosizerResult(Object):
92 n_records_stacked = Int.T(optional=True, default=1)
93 t_stack = Float.T(optional=True, default=0.)
96class Result(SeismosizerResult):
97 trace = SeismosizerTrace.T(optional=True)
98 n_shared_stacking = Int.T(optional=True, default=1)
99 t_optimize = Float.T(optional=True, default=0.)
102class StaticResult(SeismosizerResult):
103 result = Dict.T(
104 String.T(),
105 Array.T(shape=(None,), dtype=float, serialize_as='base64'))
108class GNSSCampaignResult(StaticResult):
109 campaign = gnss.GNSSCampaign.T(
110 optional=True)
113class SatelliteResult(StaticResult):
115 theta = Array.T(
116 optional=True,
117 shape=(None,), dtype=float, serialize_as='base64')
119 phi = Array.T(
120 optional=True,
121 shape=(None,), dtype=float, serialize_as='base64')
124class KiteSceneResult(SatelliteResult):
126 shape = Tuple.T()
128 def get_scene(self, component='displacement.los'):
129 try:
130 from kite import Scene
131 except ImportError:
132 raise ImportError('Kite not installed')
134 def reshape(arr):
135 return arr.reshape(self.shape)
137 displacement = self.result[component]
139 displacement, theta, phi = map(
140 reshape, (displacement, self.theta, self.phi))
142 sc = Scene(
143 displacement=displacement,
144 theta=theta, phi=phi,
145 config=self.config)
147 return sc
150class ComponentSchemeDescription(Object):
151 name = String.T()
152 source_terms = List.T(String.T())
153 ncomponents = Int.T()
154 default_stored_quantity = String.T(optional=True)
155 provided_components = List.T(String.T())
158component_scheme_descriptions = [
159 ComponentSchemeDescription(
160 name='elastic2',
161 source_terms=['m0'],
162 ncomponents=2,
163 default_stored_quantity='displacement',
164 provided_components=[
165 'n', 'e', 'd']),
166 ComponentSchemeDescription(
167 name='elastic8',
168 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
169 ncomponents=8,
170 default_stored_quantity='displacement',
171 provided_components=[
172 'n', 'e', 'd']),
173 ComponentSchemeDescription(
174 name='elastic10',
175 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
176 ncomponents=10,
177 default_stored_quantity='displacement',
178 provided_components=[
179 'n', 'e', 'd']),
180 ComponentSchemeDescription(
181 name='elastic18',
182 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],
183 ncomponents=18,
184 default_stored_quantity='displacement',
185 provided_components=[
186 'n', 'e', 'd']),
187 ComponentSchemeDescription(
188 name='elastic5',
189 source_terms=['fn', 'fe', 'fd'],
190 ncomponents=5,
191 default_stored_quantity='displacement',
192 provided_components=[
193 'n', 'e', 'd']),
194 ComponentSchemeDescription(
195 name='poroelastic10',
196 source_terms=['pore_pressure'],
197 ncomponents=10,
198 default_stored_quantity=None,
199 provided_components=[
200 'displacement.n', 'displacement.e', 'displacement.d',
201 'vertical_tilt.n', 'vertical_tilt.e',
202 'pore_pressure',
203 'darcy_velocity.n', 'darcy_velocity.e', 'darcy_velocity.d'])]
206# new names?
207# 'mt_to_displacement_1d'
208# 'mt_to_displacement_1d_ff_only'
209# 'mt_to_gravity_1d'
210# 'mt_to_stress_1d'
211# 'explosion_to_displacement_1d'
212# 'sf_to_displacement_1d'
213# 'mt_to_displacement_3d'
214# 'mt_to_gravity_3d'
215# 'mt_to_stress_3d'
216# 'pore_pressure_to_displacement_1d'
217# 'pore_pressure_to_vertical_tilt_1d'
218# 'pore_pressure_to_pore_pressure_1d'
219# 'pore_pressure_to_darcy_velocity_1d'
222component_schemes = [c.name for c in component_scheme_descriptions]
223component_scheme_to_description = dict(
224 (c.name, c) for c in component_scheme_descriptions)
227class ComponentScheme(StringChoice):
228 '''
229 Different Green's Function component schemes are available:
231 ================= =========================================================
232 Name Description
233 ================= =========================================================
234 ``elastic10`` Elastodynamic for
235 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
236 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT
237 sources only
238 ``elastic8`` Elastodynamic for far-field only
239 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
240 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores,
241 MT sources only
242 ``elastic2`` Elastodynamic for
243 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and
244 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, purely
245 isotropic sources only
246 ``elastic5`` Elastodynamic for
247 :py:class:`~pyrocko.gf.meta.ConfigTypeA`
248 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, SF
249 sources only
250 ``elastic18`` Elastodynamic for
251 :py:class:`~pyrocko.gf.meta.ConfigTypeC` stores, MT
252 sources only
253 ``poroelastic10`` Poroelastic for :py:class:`~pyrocko.gf.meta.ConfigTypeA`
254 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores
255 ================= =========================================================
256 '''
258 choices = component_schemes
261class Earthmodel1D(Object):
262 dummy_for = cake.LayeredModel
264 class __T(TBase):
265 def regularize_extra(self, val):
266 if isinstance(val, str):
267 val = cake.LayeredModel.from_scanlines(
268 cake.read_nd_model_str(val))
270 return val
272 def to_save(self, val):
273 return literal(cake.write_nd_model_str(val))
276class StringID(StringPattern):
277 pattern = r'^[A-Za-z][A-Za-z0-9._]{0,64}$'
280class ScopeType(StringChoice):
281 choices = [
282 'global',
283 'regional',
284 'local',
285 ]
288class WaveType(StringChoice):
289 choices = [
290 'full waveform',
291 'bodywave',
292 'P wave',
293 'S wave',
294 'surface wave',
295 ]
298class NearfieldTermsType(StringChoice):
299 choices = [
300 'complete',
301 'incomplete',
302 'missing',
303 ]
306class QuantityType(StringChoice):
307 choices = [
308 'displacement',
309 'rotation',
310 'velocity',
311 'acceleration',
312 'pressure',
313 'tilt',
314 'pore_pressure',
315 'darcy_velocity',
316 'vertical_tilt']
319class Reference(Object):
320 id = StringID.T()
321 type = String.T()
322 title = Unicode.T()
323 journal = Unicode.T(optional=True)
324 volume = Unicode.T(optional=True)
325 number = Unicode.T(optional=True)
326 pages = Unicode.T(optional=True)
327 year = String.T()
328 note = Unicode.T(optional=True)
329 issn = String.T(optional=True)
330 doi = String.T(optional=True)
331 url = String.T(optional=True)
332 eprint = String.T(optional=True)
333 authors = List.T(Unicode.T())
334 publisher = Unicode.T(optional=True)
335 keywords = Unicode.T(optional=True)
336 note = Unicode.T(optional=True)
337 abstract = Unicode.T(optional=True)
339 @classmethod
340 def from_bibtex(cls, filename=None, stream=None):
342 from pybtex.database.input import bibtex
344 parser = bibtex.Parser()
346 if filename is not None:
347 bib_data = parser.parse_file(filename)
348 elif stream is not None:
349 bib_data = parser.parse_stream(stream)
351 references = []
353 for id_, entry in bib_data.entries.items():
354 d = {}
355 avail = entry.fields.keys()
356 for prop in cls.T.properties:
357 if prop.name == 'authors' or prop.name not in avail:
358 continue
360 d[prop.name] = entry.fields[prop.name]
362 if 'author' in entry.persons:
363 d['authors'] = []
364 for person in entry.persons['author']:
365 d['authors'].append(new_str(person))
367 c = Reference(id=id_, type=entry.type, **d)
368 references.append(c)
370 return references
373_fpat = r'[+-]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?'
374_spat = StringID.pattern[1:-1]
375_cake_pat = cake.PhaseDef.allowed_characters_pattern
376_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic
378_ppat = r'(' \
379 r'cake:' + _cake_pat + \
380 r'|iaspei:' + _iaspei_pat + \
381 r'|vel_surface:' + _fpat + \
382 r'|vel:' + _fpat + \
383 r'|stored:' + _spat + \
384 r')'
387timing_regex_old = re.compile(
388 r'^((first|last)?\((' + _spat + r'(\|' + _spat + r')*)\)|(' +
389 _spat + r'))?(' + _fpat + ')?$')
392timing_regex = re.compile(
393 r'^((first|last)?\{(' + _ppat + r'(\|' + _ppat + r')*)\}|(' +
394 _ppat + r'))?(' + _fpat + '(S|%)?)?$')
397class PhaseSelect(StringChoice):
398 choices = ['', 'first', 'last']
401class InvalidTimingSpecification(ValidationError):
402 pass
405class Timing(SObject):
406 '''
407 Definition of a time instant relative to one or more named phase arrivals.
409 Instances of this class can be used e.g. in cutting and tapering
410 operations. They can hold an absolute time or an offset to a named phase
411 arrival or group of such arrivals.
413 Timings can be instantiated from a simple string defintion i.e. with
414 ``Timing(str)`` where ``str`` is something like
415 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``,
416 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase
417 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'``
418 is appended, it is interpreted as percent. If the an ``'S'`` is appended
419 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km].
421 Phase definitions can be specified in either of the following ways:
423 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table
424 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake
425 (see :py:class:`pyrocko.cake.PhaseDef`)
426 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance /
427 velocity [km/s]
428 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s]
430 **Examples:**
432 * ``'100'`` : absolute time; 100 s
433 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to
434 stored travel time table named ``'P'``
435 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to
436 stored travel time table named ``'P'``
437 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to
438 stored travel time table named ``'P'``
439 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is
440 undefined for a given geometry
441 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A
442 and B is chosen, if both phases are defined for a given geometry
443 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen
444 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of
445 arrivals A, B, and C
446 '''
448 def __init__(self, s=None, **kwargs):
450 if s is not None:
451 offset_is = None
452 s = re.sub(r'\s+', '', s)
453 try:
454 offset = float(s.rstrip('S'))
456 if s.endswith('S'):
457 offset_is = 'slowness'
459 phase_defs = []
460 select = ''
462 except ValueError:
463 matched = False
464 m = timing_regex.match(s)
465 if m:
466 if m.group(3):
467 phase_defs = m.group(3).split('|')
468 elif m.group(19):
469 phase_defs = [m.group(19)]
470 else:
471 phase_defs = []
473 select = m.group(2) or ''
475 offset = 0.0
476 soff = m.group(27)
477 if soff:
478 offset = float(soff.rstrip('S%'))
479 if soff.endswith('S'):
480 offset_is = 'slowness'
481 elif soff.endswith('%'):
482 offset_is = 'percent'
484 matched = True
486 else:
487 m = timing_regex_old.match(s)
488 if m:
489 if m.group(3):
490 phase_defs = m.group(3).split('|')
491 elif m.group(5):
492 phase_defs = [m.group(5)]
493 else:
494 phase_defs = []
496 select = m.group(2) or ''
498 offset = 0.0
499 if m.group(6):
500 offset = float(m.group(6))
502 phase_defs = [
503 'stored:' + phase_def for phase_def in phase_defs]
505 matched = True
507 if not matched:
508 raise InvalidTimingSpecification(s)
510 kwargs = dict(
511 phase_defs=phase_defs,
512 select=select,
513 offset=offset,
514 offset_is=offset_is)
516 SObject.__init__(self, **kwargs)
518 def __str__(self):
519 s = []
520 if self.phase_defs:
521 sphases = '|'.join(self.phase_defs)
522 # if len(self.phase_defs) > 1 or self.select:
523 sphases = '{%s}' % sphases
525 if self.select:
526 sphases = self.select + sphases
528 s.append(sphases)
530 if self.offset != 0.0 or not self.phase_defs:
531 s.append('%+g' % self.offset)
532 if self.offset_is == 'slowness':
533 s.append('S')
534 elif self.offset_is == 'percent':
535 s.append('%')
537 return ''.join(s)
539 def evaluate(self, get_phase, args):
540 try:
541 if self.offset_is == 'slowness' and self.offset != 0.0:
542 phase_offset = get_phase(
543 'vel_surface:%g' % (1.0/self.offset))
544 offset = phase_offset(args)
545 else:
546 offset = self.offset
548 if self.phase_defs:
549 phases = [
550 get_phase(phase_def) for phase_def in self.phase_defs]
551 times = [phase(args) for phase in phases]
552 if self.offset_is == 'percent':
553 times = [t*(1.+offset/100.)
554 for t in times if t is not None]
555 else:
556 times = [t+offset for t in times if t is not None]
558 if not times:
559 return None
560 elif self.select == 'first':
561 return min(times)
562 elif self.select == 'last':
563 return max(times)
564 else:
565 return times[0]
566 else:
567 return offset
569 except spit.OutOfBounds:
570 raise OutOfBounds(args)
572 phase_defs = List.T(String.T())
573 offset = Float.T(default=0.0)
574 offset_is = String.T(optional=True)
575 select = PhaseSelect.T(
576 default='',
577 help=('Can be either ``\'%s\'``, ``\'%s\'``, or ``\'%s\'``. ' %
578 tuple(PhaseSelect.choices)))
581def mkdefs(s):
582 defs = []
583 for x in s.split(','):
584 try:
585 defs.append(float(x))
586 except ValueError:
587 if x.startswith('!'):
588 defs.extend(cake.PhaseDef.classic(x[1:]))
589 else:
590 defs.append(cake.PhaseDef(x))
592 return defs
595class TPDef(Object):
596 '''
597 Maps an arrival phase identifier to an arrival phase definition.
598 '''
600 id = StringID.T(
601 help='name used to identify the phase')
602 definition = String.T(
603 help='definition of the phase in either cake syntax as defined in '
604 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an '
605 '``!``, as a *classic phase name*, or, if it is a simple '
606 'number, as a constant horizontal velocity.')
608 @property
609 def phases(self):
610 return [x for x in mkdefs(self.definition)
611 if isinstance(x, cake.PhaseDef)]
613 @property
614 def horizontal_velocities(self):
615 return [x for x in mkdefs(self.definition) if isinstance(x, float)]
618class OutOfBounds(Exception):
619 def __init__(self, values=None, reason=None):
620 Exception.__init__(self)
621 self.values = values
622 self.reason = reason
623 self.context = None
625 def __str__(self):
626 scontext = ''
627 if self.context:
628 scontext = '\n%s' % str(self.context)
630 if self.reason:
631 scontext += '\n%s' % self.reason
633 if self.values:
634 return 'out of bounds: (%s)%s' % (
635 ','.join('%g' % x for x in self.values), scontext)
636 else:
637 return 'out of bounds%s' % scontext
640class MultiLocation(Object):
642 lats = Array.T(
643 optional=True, shape=(None,), dtype=float,
644 help='Latitudes of targets.')
645 lons = Array.T(
646 optional=True, shape=(None,), dtype=float,
647 help='Longitude of targets.')
648 north_shifts = Array.T(
649 optional=True, shape=(None,), dtype=float,
650 help='North shifts of targets.')
651 east_shifts = Array.T(
652 optional=True, shape=(None,), dtype=float,
653 help='East shifts of targets.')
654 elevation = Array.T(
655 optional=True, shape=(None,), dtype=float,
656 help='Elevations of targets.')
658 def __init__(self, *args, **kwargs):
659 self._coords5 = None
660 Object.__init__(self, *args, **kwargs)
662 @property
663 def coords5(self):
664 if self._coords5 is not None:
665 return self._coords5
666 props = [self.lats, self.lons, self.north_shifts, self.east_shifts,
667 self.elevation]
668 sizes = [p.size for p in props if p is not None]
669 if not sizes:
670 raise AttributeError('Empty StaticTarget')
671 if num.unique(sizes).size != 1:
672 raise AttributeError('Inconsistent coordinate sizes.')
674 self._coords5 = num.zeros((sizes[0], 5))
675 for idx, p in enumerate(props):
676 if p is not None:
677 self._coords5[:, idx] = p.flatten()
679 return self._coords5
681 @property
682 def ncoords(self):
683 if self.coords5.shape[0] is None:
684 return 0
685 return int(self.coords5.shape[0])
687 def get_latlon(self):
688 '''
689 Get all coordinates as lat/lon.
691 :returns: Coordinates in Latitude, Longitude
692 :rtype: :class:`numpy.ndarray`, (N, 2)
693 '''
694 latlons = num.empty((self.ncoords, 2))
695 for i in range(self.ncoords):
696 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4])
697 return latlons
699 def get_corner_coords(self):
700 '''
701 Returns the corner coordinates of the multi-location object.
703 :returns: In lat/lon: lower left, upper left, upper right, lower right
704 :rtype: tuple
705 '''
706 latlon = self.get_latlon()
707 ll = (latlon[:, 0].min(), latlon[:, 1].min())
708 ur = (latlon[:, 0].max(), latlon[:, 1].max())
709 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1]))
712class Receiver(Location):
713 codes = Tuple.T(
714 3, String.T(),
715 optional=True,
716 help='network, station, and location codes')
718 def pyrocko_station(self):
719 from pyrocko import model
721 lat, lon = self.effective_latlon
723 return model.Station(
724 network=self.codes[0],
725 station=self.codes[1],
726 location=self.codes[2],
727 lat=lat,
728 lon=lon,
729 depth=self.depth)
732def g(x, d):
733 if x is None:
734 return d
735 else:
736 return x
739class UnavailableScheme(Exception):
740 pass
743class InvalidNComponents(Exception):
744 pass
747class DiscretizedSource(Object):
748 '''
749 Base class for discretized sources.
751 To compute synthetic seismograms, the parameterized source models (see of
752 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first
753 discretized into a number of point sources. These spacio-temporal point
754 source distributions are represented by subclasses of the
755 :py:class:`DiscretizedSource`. For elastodynamic problems there is the
756 :py:class:`DiscretizedMTSource` for moment tensor point source
757 distributions and the :py:class:`DiscretizedExplosionSource` for pure
758 explosion/implosion type source distributions. The geometry calculations
759 are implemented in the base class. How Green's function components have to
760 be weighted and sumed is defined in the derived classes.
762 Like in the :py:class:`Location` class, the positions of the point sources
763 contained in the discretized source are defined by a common reference point
764 (:py:attr:`lat`, :py:attr:`lon`) and cartesian offsets to this
765 (:py:attr:`north_shifts`, :py:attr:`east_shifts`, :py:attr:`depths`).
766 Alternatively latitude and longitude of each contained point source can be
767 specified directly (:py:attr:`lats`, :py:attr:`lons`).
768 '''
770 times = Array.T(shape=(None,), dtype=float)
771 lats = Array.T(shape=(None,), dtype=float, optional=True)
772 lons = Array.T(shape=(None,), dtype=float, optional=True)
773 lat = Float.T(optional=True)
774 lon = Float.T(optional=True)
775 north_shifts = Array.T(shape=(None,), dtype=float, optional=True)
776 east_shifts = Array.T(shape=(None,), dtype=float, optional=True)
777 depths = Array.T(shape=(None,), dtype=float)
779 @classmethod
780 def check_scheme(cls, scheme):
781 '''
782 Check if given GF component scheme is supported by the class.
784 Raises :py:class:`UnavailableScheme` if the given scheme is not
785 supported by this discretized source class.
786 '''
788 if scheme not in cls.provided_schemes:
789 raise UnavailableScheme(
790 'source type "%s" does not support GF component scheme "%s"' %
791 (cls.__name__, scheme))
793 def __init__(self, **kwargs):
794 Object.__init__(self, **kwargs)
795 self._latlons = None
797 def __setattr__(self, name, value):
798 if name in ('lat', 'lon', 'north_shifts', 'east_shifts',
799 'lats', 'lons'):
800 self.__dict__['_latlons'] = None
802 Object.__setattr__(self, name, value)
804 def get_source_terms(self, scheme):
805 raise NotImplementedError()
807 def make_weights(self, receiver, scheme):
808 raise NotImplementedError()
810 @property
811 def effective_latlons(self):
812 '''
813 Property holding the offest-corrected lats and lons of all points.
814 '''
816 if self._latlons is None:
817 if self.lats is not None and self.lons is not None:
818 if (self.north_shifts is not None and
819 self.east_shifts is not None):
820 self._latlons = orthodrome.ne_to_latlon(
821 self.lats, self.lons,
822 self.north_shifts, self.east_shifts)
823 else:
824 self._latlons = self.lats, self.lons
825 else:
826 lat = g(self.lat, 0.0)
827 lon = g(self.lon, 0.0)
828 self._latlons = orthodrome.ne_to_latlon(
829 lat, lon, self.north_shifts, self.east_shifts)
831 return self._latlons
833 @property
834 def effective_north_shifts(self):
835 if self.north_shifts is not None:
836 return self.north_shifts
837 else:
838 return num.zeros(self.times.size)
840 @property
841 def effective_east_shifts(self):
842 if self.east_shifts is not None:
843 return self.east_shifts
844 else:
845 return num.zeros(self.times.size)
847 def same_origin(self, receiver):
848 '''
849 Check if receiver has same reference point.
850 '''
852 return (g(self.lat, 0.0) == receiver.lat and
853 g(self.lon, 0.0) == receiver.lon and
854 self.lats is None and self.lons is None)
856 def azibazis_to(self, receiver):
857 '''
858 Compute azimuths and backaziumuths to/from receiver, for all contained
859 points.
860 '''
862 if self.same_origin(receiver):
863 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts,
864 receiver.north_shift - self.north_shifts)
865 bazis = azis + 180.
866 else:
867 slats, slons = self.effective_latlons
868 rlat, rlon = receiver.effective_latlon
869 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon)
870 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons)
872 return azis, bazis
874 def distances_to(self, receiver):
875 '''
876 Compute distances to receiver for all contained points.
877 '''
878 if self.same_origin(receiver):
879 return num.sqrt((self.north_shifts - receiver.north_shift)**2 +
880 (self.east_shifts - receiver.east_shift)**2)
882 else:
883 slats, slons = self.effective_latlons
884 rlat, rlon = receiver.effective_latlon
885 return orthodrome.distance_accurate50m_numpy(slats, slons,
886 rlat, rlon)
888 def element_coords(self, i):
889 if self.lats is not None and self.lons is not None:
890 lat = float(self.lats[i])
891 lon = float(self.lons[i])
892 else:
893 lat = self.lat
894 lon = self.lon
896 if self.north_shifts is not None and self.east_shifts is not None:
897 north_shift = float(self.north_shifts[i])
898 east_shift = float(self.east_shifts[i])
900 else:
901 north_shift = east_shift = 0.0
903 return lat, lon, north_shift, east_shift
905 def coords5(self):
906 xs = num.zeros((self.nelements, 5))
908 if self.lats is not None and self.lons is not None:
909 xs[:, 0] = self.lats
910 xs[:, 1] = self.lons
911 else:
912 xs[:, 0] = g(self.lat, 0.0)
913 xs[:, 1] = g(self.lon, 0.0)
915 if self.north_shifts is not None and self.east_shifts is not None:
916 xs[:, 2] = self.north_shifts
917 xs[:, 3] = self.east_shifts
919 xs[:, 4] = self.depths
921 return xs
923 @property
924 def nelements(self):
925 return self.times.size
927 @classmethod
928 def combine(cls, sources, **kwargs):
929 '''
930 Combine several discretized source models.
932 Concatenenates all point sources in the given discretized ``sources``.
933 Care must be taken when using this function that the external amplitude
934 factors and reference times of the parameterized (undiscretized)
935 sources match or are accounted for.
936 '''
938 first = sources[0]
940 if not all(type(s) == type(first) for s in sources):
941 raise Exception('DiscretizedSource.combine must be called with '
942 'sources of same type.')
944 latlons = []
945 for s in sources:
946 latlons.append(s.effective_latlons)
948 lats, lons = num.hstack(latlons)
950 if all((s.lats is None and s.lons is None) for s in sources):
951 rlats = num.array([s.lat for s in sources], dtype=float)
952 rlons = num.array([s.lon for s in sources], dtype=float)
953 same_ref = num.all(
954 rlats == rlats[0]) and num.all(rlons == rlons[0])
955 else:
956 same_ref = False
958 cat = num.concatenate
959 times = cat([s.times for s in sources])
960 depths = cat([s.depths for s in sources])
962 if same_ref:
963 lat = first.lat
964 lon = first.lon
965 north_shifts = cat([s.effective_north_shifts for s in sources])
966 east_shifts = cat([s.effective_east_shifts for s in sources])
967 lats = None
968 lons = None
969 else:
970 lat = None
971 lon = None
972 north_shifts = None
973 east_shifts = None
975 return cls(
976 times=times, lat=lat, lon=lon, lats=lats, lons=lons,
977 north_shifts=north_shifts, east_shifts=east_shifts,
978 depths=depths, **kwargs)
980 def centroid_position(self):
981 moments = self.moments()
982 norm = num.sum(moments)
983 if norm != 0.0:
984 w = moments / num.sum(moments)
985 else:
986 w = num.ones(moments.size)
988 if self.lats is not None and self.lons is not None:
989 lats, lons = self.effective_latlons
990 rlat, rlon = num.mean(lats), num.mean(lons)
991 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)
992 else:
993 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)
994 n, e = self.north_shifts, self.east_shifts
996 cn = num.sum(n*w)
997 ce = num.sum(e*w)
998 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)
1000 if self.lats is not None and self.lons is not None:
1001 lat = clat
1002 lon = clon
1003 north_shift = 0.
1004 east_shift = 0.
1005 else:
1006 lat = g(self.lat, 0.0)
1007 lon = g(self.lon, 0.0)
1008 north_shift = cn
1009 east_shift = ce
1011 depth = num.sum(self.depths*w)
1012 time = num.sum(self.times*w)
1013 return tuple(float(x) for x in
1014 (time, lat, lon, north_shift, east_shift, depth))
1017class DiscretizedExplosionSource(DiscretizedSource):
1018 m0s = Array.T(shape=(None,), dtype=float)
1020 provided_schemes = (
1021 'elastic2',
1022 'elastic8',
1023 'elastic10',
1024 )
1026 def get_source_terms(self, scheme):
1027 self.check_scheme(scheme)
1029 if scheme == 'elastic2':
1030 return self.m0s[:, num.newaxis].copy()
1032 elif scheme in ('elastic8', 'elastic10'):
1033 m6s = num.zeros((self.m0s.size, 6))
1034 m6s[:, 0:3] = self.m0s[:, num.newaxis]
1035 return m6s
1036 else:
1037 assert False
1039 def make_weights(self, receiver, scheme):
1040 self.check_scheme(scheme)
1042 azis, bazis = self.azibazis_to(receiver)
1044 sb = num.sin(bazis*d2r-num.pi)
1045 cb = num.cos(bazis*d2r-num.pi)
1047 m0s = self.m0s
1048 n = azis.size
1050 cat = num.concatenate
1051 rep = num.repeat
1053 if scheme == 'elastic2':
1054 w_n = cb*m0s
1055 g_n = filledi(0, n)
1056 w_e = sb*m0s
1057 g_e = filledi(0, n)
1058 w_d = m0s
1059 g_d = filledi(1, n)
1061 elif scheme == 'elastic8':
1062 w_n = cat((cb*m0s, cb*m0s))
1063 g_n = rep((0, 2), n)
1064 w_e = cat((sb*m0s, sb*m0s))
1065 g_e = rep((0, 2), n)
1066 w_d = cat((m0s, m0s))
1067 g_d = rep((5, 7), n)
1069 elif scheme == 'elastic10':
1070 w_n = cat((cb*m0s, cb*m0s, cb*m0s))
1071 g_n = rep((0, 2, 8), n)
1072 w_e = cat((sb*m0s, sb*m0s, sb*m0s))
1073 g_e = rep((0, 2, 8), n)
1074 w_d = cat((m0s, m0s, m0s))
1075 g_d = rep((5, 7, 9), n)
1077 else:
1078 assert False
1080 return (
1081 ('displacement.n', w_n, g_n),
1082 ('displacement.e', w_e, g_e),
1083 ('displacement.d', w_d, g_d),
1084 )
1086 def split(self):
1087 from pyrocko.gf.seismosizer import ExplosionSource
1088 sources = []
1089 for i in range(self.nelements):
1090 lat, lon, north_shift, east_shift = self.element_coords(i)
1091 sources.append(ExplosionSource(
1092 time=float(self.times[i]),
1093 lat=lat,
1094 lon=lon,
1095 north_shift=north_shift,
1096 east_shift=east_shift,
1097 depth=float(self.depths[i]),
1098 moment=float(self.m0s[i])))
1100 return sources
1102 def moments(self):
1103 return self.m0s
1105 def centroid(self):
1106 from pyrocko.gf.seismosizer import ExplosionSource
1107 time, lat, lon, north_shift, east_shift, depth = \
1108 self.centroid_position()
1110 return ExplosionSource(
1111 time=time,
1112 lat=lat,
1113 lon=lon,
1114 north_shift=north_shift,
1115 east_shift=east_shift,
1116 depth=depth,
1117 moment=float(num.sum(self.m0s)))
1119 @classmethod
1120 def combine(cls, sources, **kwargs):
1121 '''
1122 Combine several discretized source models.
1124 Concatenenates all point sources in the given discretized ``sources``.
1125 Care must be taken when using this function that the external amplitude
1126 factors and reference times of the parameterized (undiscretized)
1127 sources match or are accounted for.
1128 '''
1130 if 'm0s' not in kwargs:
1131 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])
1133 return super(DiscretizedExplosionSource, cls).combine(sources,
1134 **kwargs)
1137class DiscretizedSFSource(DiscretizedSource):
1138 forces = Array.T(shape=(None, 3), dtype=float)
1140 provided_schemes = (
1141 'elastic5',
1142 )
1144 def get_source_terms(self, scheme):
1145 self.check_scheme(scheme)
1147 return self.forces
1149 def make_weights(self, receiver, scheme):
1150 self.check_scheme(scheme)
1152 azis, bazis = self.azibazis_to(receiver)
1154 sa = num.sin(azis*d2r)
1155 ca = num.cos(azis*d2r)
1156 sb = num.sin(bazis*d2r-num.pi)
1157 cb = num.cos(bazis*d2r-num.pi)
1159 forces = self.forces
1160 fn = forces[:, 0]
1161 fe = forces[:, 1]
1162 fd = forces[:, 2]
1164 f0 = fd
1165 f1 = ca * fn + sa * fe
1166 f2 = ca * fe - sa * fn
1168 n = azis.size
1170 if scheme == 'elastic5':
1171 ioff = 0
1173 cat = num.concatenate
1174 rep = num.repeat
1176 w_n = cat((cb*f0, cb*f1, -sb*f2))
1177 g_n = ioff + rep((0, 1, 2), n)
1178 w_e = cat((sb*f0, sb*f1, cb*f2))
1179 g_e = ioff + rep((0, 1, 2), n)
1180 w_d = cat((f0, f1))
1181 g_d = ioff + rep((3, 4), n)
1183 return (
1184 ('displacement.n', w_n, g_n),
1185 ('displacement.e', w_e, g_e),
1186 ('displacement.d', w_d, g_d),
1187 )
1189 @classmethod
1190 def combine(cls, sources, **kwargs):
1191 '''
1192 Combine several discretized source models.
1194 Concatenenates all point sources in the given discretized ``sources``.
1195 Care must be taken when using this function that the external amplitude
1196 factors and reference times of the parameterized (undiscretized)
1197 sources match or are accounted for.
1198 '''
1200 if 'forces' not in kwargs:
1201 kwargs['forces'] = num.vstack([s.forces for s in sources])
1203 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)
1205 def moments(self):
1206 return num.sum(self.forces**2, axis=1)
1208 def centroid(self):
1209 from pyrocko.gf.seismosizer import SFSource
1210 time, lat, lon, north_shift, east_shift, depth = \
1211 self.centroid_position()
1213 fn, fe, fd = map(float, num.sum(self.forces, axis=0))
1214 return SFSource(
1215 time=time,
1216 lat=lat,
1217 lon=lon,
1218 north_shift=north_shift,
1219 east_shift=east_shift,
1220 depth=depth,
1221 fn=fn,
1222 fe=fe,
1223 fd=fd)
1226class DiscretizedMTSource(DiscretizedSource):
1227 m6s = Array.T(
1228 shape=(None, 6), dtype=float,
1229 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')
1231 provided_schemes = (
1232 'elastic8',
1233 'elastic10',
1234 'elastic18',
1235 )
1237 def get_source_terms(self, scheme):
1238 self.check_scheme(scheme)
1239 return self.m6s
1241 def make_weights(self, receiver, scheme):
1242 self.check_scheme(scheme)
1244 m6s = self.m6s
1245 n = m6s.shape[0]
1246 rep = num.repeat
1248 if scheme == 'elastic18':
1249 w_n = m6s.flatten()
1250 g_n = num.tile(num.arange(0, 6), n)
1251 w_e = m6s.flatten()
1252 g_e = num.tile(num.arange(6, 12), n)
1253 w_d = m6s.flatten()
1254 g_d = num.tile(num.arange(12, 18), n)
1256 else:
1257 azis, bazis = self.azibazis_to(receiver)
1259 sa = num.sin(azis*d2r)
1260 ca = num.cos(azis*d2r)
1261 s2a = num.sin(2.*azis*d2r)
1262 c2a = num.cos(2.*azis*d2r)
1263 sb = num.sin(bazis*d2r-num.pi)
1264 cb = num.cos(bazis*d2r-num.pi)
1266 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a
1267 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa
1268 f2 = m6s[:, 2]
1269 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a
1270 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa
1271 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a
1273 cat = num.concatenate
1275 if scheme == 'elastic8':
1276 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))
1277 g_n = rep((0, 1, 2, 3, 4), n)
1278 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))
1279 g_e = rep((0, 1, 2, 3, 4), n)
1280 w_d = cat((f0, f1, f2))
1281 g_d = rep((5, 6, 7), n)
1283 elif scheme == 'elastic10':
1284 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))
1285 g_n = rep((0, 1, 2, 8, 3, 4), n)
1286 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))
1287 g_e = rep((0, 1, 2, 8, 3, 4), n)
1288 w_d = cat((f0, f1, f2, f5))
1289 g_d = rep((5, 6, 7, 9), n)
1291 else:
1292 assert False
1294 return (
1295 ('displacement.n', w_n, g_n),
1296 ('displacement.e', w_e, g_e),
1297 ('displacement.d', w_d, g_d),
1298 )
1300 def split(self):
1301 from pyrocko.gf.seismosizer import MTSource
1302 sources = []
1303 for i in range(self.nelements):
1304 lat, lon, north_shift, east_shift = self.element_coords(i)
1305 sources.append(MTSource(
1306 time=float(self.times[i]),
1307 lat=lat,
1308 lon=lon,
1309 north_shift=north_shift,
1310 east_shift=east_shift,
1311 depth=float(self.depths[i]),
1312 m6=self.m6s[i]))
1314 return sources
1316 def moments(self):
1317 n = self.nelements
1318 moments = num.zeros(n)
1319 for i in range(n):
1320 m = moment_tensor.symmat6(*self.m6s[i])
1321 m_evals = num.linalg.eigh(m)[0]
1323 # incorrect for non-dc sources: !!!!
1324 m0 = num.linalg.norm(m_evals)/math.sqrt(2.)
1325 moments[i] = m0
1327 return moments
1329 def centroid(self):
1330 from pyrocko.gf.seismosizer import MTSource
1331 time, lat, lon, north_shift, east_shift, depth = \
1332 self.centroid_position()
1334 return MTSource(
1335 time=time,
1336 lat=lat,
1337 lon=lon,
1338 north_shift=north_shift,
1339 east_shift=east_shift,
1340 depth=depth,
1341 m6=num.sum(self.m6s, axis=0))
1343 @classmethod
1344 def combine(cls, sources, **kwargs):
1345 '''
1346 Combine several discretized source models.
1348 Concatenenates all point sources in the given discretized ``sources``.
1349 Care must be taken when using this function that the external amplitude
1350 factors and reference times of the parameterized (undiscretized)
1351 sources match or are accounted for.
1352 '''
1354 if 'm6s' not in kwargs:
1355 kwargs['m6s'] = num.vstack([s.m6s for s in sources])
1357 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)
1360class DiscretizedPorePressureSource(DiscretizedSource):
1361 pp = Array.T(shape=(None,), dtype=float)
1363 provided_schemes = (
1364 'poroelastic10',
1365 )
1367 def get_source_terms(self, scheme):
1368 self.check_scheme(scheme)
1369 return self.pp[:, num.newaxis].copy()
1371 def make_weights(self, receiver, scheme):
1372 self.check_scheme(scheme)
1374 azis, bazis = self.azibazis_to(receiver)
1376 sb = num.sin(bazis*d2r-num.pi)
1377 cb = num.cos(bazis*d2r-num.pi)
1379 pp = self.pp
1380 n = bazis.size
1382 w_un = cb*pp
1383 g_un = filledi(1, n)
1384 w_ue = sb*pp
1385 g_ue = filledi(1, n)
1386 w_ud = pp
1387 g_ud = filledi(0, n)
1389 w_tn = cb*pp
1390 g_tn = filledi(6, n)
1391 w_te = sb*pp
1392 g_te = filledi(6, n)
1394 w_pp = pp
1395 g_pp = filledi(7, n)
1397 w_dvn = cb*pp
1398 g_dvn = filledi(9, n)
1399 w_dve = sb*pp
1400 g_dve = filledi(9, n)
1401 w_dvd = pp
1402 g_dvd = filledi(8, n)
1404 return (
1405 ('displacement.n', w_un, g_un),
1406 ('displacement.e', w_ue, g_ue),
1407 ('displacement.d', w_ud, g_ud),
1408 ('vertical_tilt.n', w_tn, g_tn),
1409 ('vertical_tilt.e', w_te, g_te),
1410 ('pore_pressure', w_pp, g_pp),
1411 ('darcy_velocity.n', w_dvn, g_dvn),
1412 ('darcy_velocity.e', w_dve, g_dve),
1413 ('darcy_velocity.d', w_dvd, g_dvd),
1414 )
1416 def moments(self):
1417 return self.pp
1419 def centroid(self):
1420 from pyrocko.gf.seismosizer import PorePressurePointSource
1421 time, lat, lon, north_shift, east_shift, depth = \
1422 self.centroid_position()
1424 return PorePressurePointSource(
1425 time=time,
1426 lat=lat,
1427 lon=lon,
1428 north_shift=north_shift,
1429 east_shift=east_shift,
1430 depth=depth,
1431 pp=float(num.sum(self.pp)))
1433 @classmethod
1434 def combine(cls, sources, **kwargs):
1435 '''
1436 Combine several discretized source models.
1438 Concatenenates all point sources in the given discretized ``sources``.
1439 Care must be taken when using this function that the external amplitude
1440 factors and reference times of the parameterized (undiscretized)
1441 sources match or are accounted for.
1442 '''
1444 if 'pp' not in kwargs:
1445 kwargs['pp'] = num.concatenate([s.pp for s in sources])
1447 return super(DiscretizedPorePressureSource, cls).combine(sources,
1448 **kwargs)
1451class Region(Object):
1452 name = String.T(optional=True)
1455class RectangularRegion(Region):
1456 lat_min = Float.T()
1457 lat_max = Float.T()
1458 lon_min = Float.T()
1459 lon_max = Float.T()
1462class CircularRegion(Region):
1463 lat = Float.T()
1464 lon = Float.T()
1465 radius = Float.T()
1468class Config(Object):
1469 '''
1470 Green's function store meta information.
1472 Currently implemented :py:class:`~pyrocko.gf.store.Store`
1473 configuration types are:
1475 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or
1476 spherical symmetry, 1D earth model, single receiver depth
1478 * Problem is invariant to horizontal translations and rotations around
1479 vertical axis.
1480 * All receivers must be at the same depth (e.g. at the surface)
1481 * High level index variables: ``(source_depth, receiver_distance,
1482 component)``
1484 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or
1485 spherical symmetry, 1D earth model, variable receiver depth
1487 * Symmetries like in Type A but has additional index for receiver depth
1488 * High level index variables: ``(source_depth, receiver_distance,
1489 receiver_depth, component)``
1491 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical
1492 constraints but fixed receiver positions
1494 * Cartesian source volume around a reference point
1495 * High level index variables: ``(ireceiver, source_depth,
1496 source_east_shift, source_north_shift, component)``
1497 '''
1499 id = StringID.T(
1500 help='Name of the store. May consist of upper and lower-case letters, '
1501 'digits, dots and underscores. The name must start with a '
1502 'letter.')
1504 derived_from_id = StringID.T(
1505 optional=True,
1506 help='Name of the original store, if this store has been derived from '
1507 'another one (e.g. extracted subset).')
1509 version = String.T(
1510 default='1.0',
1511 optional=True,
1512 help='User-defined version string. Use <major>.<minor> format.')
1514 modelling_code_id = StringID.T(
1515 optional=True,
1516 help='Identifier of the backend used to compute the store.')
1518 author = Unicode.T(
1519 optional=True,
1520 help='Comma-separated list of author names.')
1522 author_email = String.T(
1523 optional=True,
1524 help="Author's contact email address.")
1526 created_time = Timestamp.T(
1527 optional=True,
1528 help='Time of creation of the store.')
1530 regions = List.T(
1531 Region.T(),
1532 help='Geographical regions for which the store is representative.')
1534 scope_type = ScopeType.T(
1535 optional=True,
1536 help='Distance range scope of the store '
1537 '(%s).' % fmt_choices(ScopeType))
1539 waveform_type = WaveType.T(
1540 optional=True,
1541 help='Wave type stored (%s).' % fmt_choices(WaveType))
1543 nearfield_terms = NearfieldTermsType.T(
1544 optional=True,
1545 help='Information about the inclusion of near-field terms in the '
1546 'modelling (%s).' % fmt_choices(NearfieldTermsType))
1548 description = String.T(
1549 optional=True,
1550 help='Free form textual description of the GF store.')
1552 references = List.T(
1553 Reference.T(),
1554 help='Reference list to cite the modelling code, earth model or '
1555 'related work.')
1557 earthmodel_1d = Earthmodel1D.T(
1558 optional=True,
1559 help='Layered earth model in ND (named discontinuity) format.')
1561 earthmodel_receiver_1d = Earthmodel1D.T(
1562 optional=True,
1563 help='Receiver-side layered earth model in ND format.')
1565 can_interpolate_source = Bool.T(
1566 optional=True,
1567 help='Hint to indicate if the spatial sampling of the store is dense '
1568 'enough for multi-linear interpolation at the source.')
1570 can_interpolate_receiver = Bool.T(
1571 optional=True,
1572 help='Hint to indicate if the spatial sampling of the store is dense '
1573 'enough for multi-linear interpolation at the receiver.')
1575 frequency_min = Float.T(
1576 optional=True,
1577 help='Hint to indicate the lower bound of valid frequencies [Hz].')
1579 frequency_max = Float.T(
1580 optional=True,
1581 help='Hint to indicate the upper bound of valid frequencies [Hz].')
1583 sample_rate = Float.T(
1584 optional=True,
1585 help='Sample rate of the GF store [Hz].')
1587 factor = Float.T(
1588 default=1.0,
1589 help='Gain value, factored out of the stored GF samples. '
1590 '(may not work properly, keep at 1.0).',
1591 optional=True)
1593 component_scheme = ComponentScheme.T(
1594 default='elastic10',
1595 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))
1597 stored_quantity = QuantityType.T(
1598 optional=True,
1599 help='Physical quantity of stored values (%s). If not given, a '
1600 'default is used based on the GF component scheme. The default '
1601 'for the ``"elastic*"`` family of component schemes is '
1602 '``"displacement"``.' % fmt_choices(QuantityType))
1604 tabulated_phases = List.T(
1605 TPDef.T(),
1606 help='Mapping of phase names to phase definitions, for which travel '
1607 'time tables are available in the GF store.')
1609 ncomponents = Int.T(
1610 optional=True,
1611 help='Number of GF components. Use :gattr:`component_scheme` instead.')
1613 uuid = String.T(
1614 optional=True,
1615 help='Heuristic hash value which can be used to uniquely identify the '
1616 'GF store for practical purposes.')
1618 reference = String.T(
1619 optional=True,
1620 help='Store reference name composed of the store\'s :gattr:`id` and '
1621 'the first six letters of its :gattr:`uuid`.')
1623 def __init__(self, **kwargs):
1624 self._do_auto_updates = False
1625 Object.__init__(self, **kwargs)
1626 self._index_function = None
1627 self._indices_function = None
1628 self._vicinity_function = None
1629 self.validate(regularize=True, depth=1)
1630 self._do_auto_updates = True
1631 self.update()
1633 def check_ncomponents(self):
1634 ncomponents = component_scheme_to_description[
1635 self.component_scheme].ncomponents
1637 if self.ncomponents is None:
1638 self.ncomponents = ncomponents
1639 elif ncomponents != self.ncomponents:
1640 raise InvalidNComponents(
1641 'ncomponents=%i incompatible with component_scheme="%s"' % (
1642 self.ncomponents, self.component_scheme))
1644 def __setattr__(self, name, value):
1645 Object.__setattr__(self, name, value)
1646 try:
1647 self.T.get_property(name)
1648 if self._do_auto_updates:
1649 self.update()
1651 except ValueError:
1652 pass
1654 def update(self):
1655 self.check_ncomponents()
1656 self._update()
1657 self._make_index_functions()
1659 def irecord(self, *args):
1660 return self._index_function(*args)
1662 def irecords(self, *args):
1663 return self._indices_function(*args)
1665 def vicinity(self, *args):
1666 return self._vicinity_function(*args)
1668 def vicinities(self, *args):
1669 return self._vicinities_function(*args)
1671 def grid_interpolation_coefficients(self, *args):
1672 return self._grid_interpolation_coefficients(*args)
1674 def nodes(self, level=None, minlevel=None):
1675 return nodes(self.coords[minlevel:level])
1677 def iter_nodes(self, level=None, minlevel=None):
1678 return nditer_outer(self.coords[minlevel:level])
1680 def iter_extraction(self, gdef, level=None):
1681 i = 0
1682 arrs = []
1683 ntotal = 1
1684 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):
1685 if gdef and len(gdef) > i:
1686 sssn = gdef[i]
1687 else:
1688 sssn = (None,)*4
1690 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))
1691 ntotal *= len(arr)
1693 arrs.append(arr)
1694 i += 1
1696 arrs.append(self.coords[-1])
1697 return nditer_outer(arrs[:level])
1699 def make_sum_params(self, source, receiver, implementation='c',
1700 nthreads=0):
1701 assert implementation in ['c', 'python']
1703 out = []
1704 delays = source.times
1705 for comp, weights, icomponents in source.make_weights(
1706 receiver,
1707 self.component_scheme):
1709 weights *= self.factor
1711 args = self.make_indexing_args(source, receiver, icomponents)
1712 delays_expanded = num.tile(delays, icomponents.size//delays.size)
1713 out.append((comp, args, delays_expanded, weights))
1715 return out
1717 def short_info(self):
1718 raise NotImplementedError('should be implemented in subclass')
1720 def get_shear_moduli(self, lat, lon, points,
1721 interpolation=None):
1722 '''
1723 Get shear moduli at given points from contained velocity model.
1725 :param lat: surface origin for coordinate system of ``points``
1726 :param points: NumPy array of shape ``(N, 3)``, where each row is
1727 a point ``(north, east, depth)``, relative to origin at
1728 ``(lat, lon)``
1729 :param interpolation: interpolation method. Choose from
1730 ``('nearest_neighbor', 'multilinear')``
1731 :returns: NumPy array of length N with extracted shear moduli at each
1732 point
1734 The default implementation retrieves and interpolates the shear moduli
1735 from the contained 1D velocity profile.
1736 '''
1737 return self.get_material_property(lat, lon, points,
1738 parameter='shear_moduli',
1739 interpolation=interpolation)
1741 def get_lambda_moduli(self, lat, lon, points,
1742 interpolation=None):
1743 '''
1744 Get lambda moduli at given points from contained velocity model.
1746 :param lat: surface origin for coordinate system of ``points``
1747 :param points: NumPy array of shape ``(N, 3)``, where each row is
1748 a point ``(north, east, depth)``, relative to origin at
1749 ``(lat, lon)``
1750 :param interpolation: interpolation method. Choose from
1751 ``('nearest_neighbor', 'multilinear')``
1752 :returns: NumPy array of length N with extracted shear moduli at each
1753 point
1755 The default implementation retrieves and interpolates the lambda moduli
1756 from the contained 1D velocity profile.
1757 '''
1758 return self.get_material_property(lat, lon, points,
1759 parameter='lambda_moduli',
1760 interpolation=interpolation)
1762 def get_bulk_moduli(self, lat, lon, points,
1763 interpolation=None):
1764 '''
1765 Get bulk moduli at given points from contained velocity model.
1767 :param lat: surface origin for coordinate system of ``points``
1768 :param points: NumPy array of shape ``(N, 3)``, where each row is
1769 a point ``(north, east, depth)``, relative to origin at
1770 ``(lat, lon)``
1771 :param interpolation: interpolation method. Choose from
1772 ``('nearest_neighbor', 'multilinear')``
1773 :returns: NumPy array of length N with extracted shear moduli at each
1774 point
1776 The default implementation retrieves and interpolates the lambda moduli
1777 from the contained 1D velocity profile.
1778 '''
1779 lambda_moduli = self.get_material_property(
1780 lat, lon, points, parameter='lambda_moduli',
1781 interpolation=interpolation)
1782 shear_moduli = self.get_material_property(
1783 lat, lon, points, parameter='shear_moduli',
1784 interpolation=interpolation)
1785 return lambda_moduli + (2 / 3) * shear_moduli
1787 def get_vs(self, lat, lon, points, interpolation=None):
1788 '''
1789 Get Vs at given points from contained velocity model.
1791 :param lat: surface origin for coordinate system of ``points``
1792 :param points: NumPy array of shape ``(N, 3)``, where each row is
1793 a point ``(north, east, depth)``, relative to origin at
1794 ``(lat, lon)``
1795 :param interpolation: interpolation method. Choose from
1796 ``('nearest_neighbor', 'multilinear')``
1797 :returns: NumPy array of length N with extracted shear moduli at each
1798 point
1800 The default implementation retrieves and interpolates Vs
1801 from the contained 1D velocity profile.
1802 '''
1803 return self.get_material_property(lat, lon, points,
1804 parameter='vs',
1805 interpolation=interpolation)
1807 def get_vp(self, lat, lon, points, interpolation=None):
1808 '''
1809 Get Vp at given points from contained velocity model.
1811 :param lat: surface origin for coordinate system of ``points``
1812 :param points: NumPy array of shape ``(N, 3)``, where each row is
1813 a point ``(north, east, depth)``, relative to origin at
1814 ``(lat, lon)``
1815 :param interpolation: interpolation method. Choose from
1816 ``('nearest_neighbor', 'multilinear')``
1817 :returns: NumPy array of length N with extracted shear moduli at each
1818 point
1820 The default implementation retrieves and interpolates Vp
1821 from the contained 1D velocity profile.
1822 '''
1823 return self.get_material_property(lat, lon, points,
1824 parameter='vp',
1825 interpolation=interpolation)
1827 def get_rho(self, lat, lon, points, interpolation=None):
1828 '''
1829 Get rho at given points from contained velocity model.
1831 :param lat: surface origin for coordinate system of ``points``
1832 :param points: NumPy array of shape ``(N, 3)``, where each row is
1833 a point ``(north, east, depth)``, relative to origin at
1834 ``(lat, lon)``
1835 :param interpolation: interpolation method. Choose from
1836 ``('nearest_neighbor', 'multilinear')``
1837 :returns: NumPy array of length N with extracted shear moduli at each
1838 point
1840 The default implementation retrieves and interpolates rho
1841 from the contained 1D velocity profile.
1842 '''
1843 return self.get_material_property(lat, lon, points,
1844 parameter='rho',
1845 interpolation=interpolation)
1847 def get_material_property(self, lat, lon, points, parameter='vs',
1848 interpolation=None):
1850 if interpolation is None:
1851 raise TypeError('Interpolation method not defined! available: '
1852 "multilinear", "nearest_neighbor")
1854 earthmod = self.earthmodel_1d
1855 store_depth_profile = self.get_source_depths()
1856 z_profile = earthmod.profile('z')
1858 if parameter == 'vs':
1859 vs_profile = earthmod.profile('vs')
1860 profile = num.interp(
1861 store_depth_profile, z_profile, vs_profile)
1863 elif parameter == 'vp':
1864 vp_profile = earthmod.profile('vp')
1865 profile = num.interp(
1866 store_depth_profile, z_profile, vp_profile)
1868 elif parameter == 'rho':
1869 rho_profile = earthmod.profile('rho')
1871 profile = num.interp(
1872 store_depth_profile, z_profile, rho_profile)
1874 elif parameter == 'shear_moduli':
1875 vs_profile = earthmod.profile('vs')
1876 rho_profile = earthmod.profile('rho')
1878 store_vs_profile = num.interp(
1879 store_depth_profile, z_profile, vs_profile)
1880 store_rho_profile = num.interp(
1881 store_depth_profile, z_profile, rho_profile)
1883 profile = num.power(store_vs_profile, 2) * store_rho_profile
1885 elif parameter == 'lambda_moduli':
1886 vs_profile = earthmod.profile('vs')
1887 vp_profile = earthmod.profile('vp')
1888 rho_profile = earthmod.profile('rho')
1890 store_vs_profile = num.interp(
1891 store_depth_profile, z_profile, vs_profile)
1892 store_vp_profile = num.interp(
1893 store_depth_profile, z_profile, vp_profile)
1894 store_rho_profile = num.interp(
1895 store_depth_profile, z_profile, rho_profile)
1897 profile = store_rho_profile * (
1898 num.power(store_vp_profile, 2) -
1899 num.power(store_vs_profile, 2) * 2)
1900 else:
1901 raise TypeError(
1902 'parameter %s not available' % parameter)
1904 if interpolation == 'multilinear':
1905 kind = 'linear'
1906 elif interpolation == 'nearest_neighbor':
1907 kind = 'nearest'
1908 else:
1909 raise TypeError(
1910 'Interpolation method %s not available' % interpolation)
1912 interpolator = interp1d(store_depth_profile, profile, kind=kind)
1914 try:
1915 return interpolator(points[:, 2])
1916 except ValueError:
1917 raise OutOfBounds()
1919 def is_static(self):
1920 for code in ('psgrn_pscmp', 'poel'):
1921 if self.modelling_code_id.startswith(code):
1922 return True
1923 return False
1925 def is_dynamic(self):
1926 return not self.is_static()
1928 def get_source_depths(self):
1929 raise NotImplementedError('must be implemented in subclass')
1931 def get_tabulated_phase(self, phase_id):
1932 '''
1933 Get tabulated phase definition.
1934 '''
1936 for pdef in self.tabulated_phases:
1937 if pdef.id == phase_id:
1938 return pdef
1940 raise StoreError('No such phase: %s' % phase_id)
1942 def fix_ttt_holes(self, sptree, mode):
1943 raise StoreError(
1944 'Cannot fix travel time table holes in GF stores of type %s.'
1945 % self.short_type)
1948class ConfigTypeA(Config):
1949 '''
1950 Cylindrical symmetry, 1D earth model, single receiver depth
1952 * Problem is invariant to horizontal translations and rotations around
1953 vertical axis.
1955 * All receivers must be at the same depth (e.g. at the surface)
1956 High level index variables: ``(source_depth, distance,
1957 component)``
1959 * The ``distance`` is the surface distance between source and receiver
1960 points.
1961 '''
1963 receiver_depth = Float.T(
1964 default=0.0,
1965 help='Fixed receiver depth [m].')
1967 source_depth_min = Float.T(
1968 help='Minimum source depth [m].')
1970 source_depth_max = Float.T(
1971 help='Maximum source depth [m].')
1973 source_depth_delta = Float.T(
1974 help='Grid spacing of source depths [m]')
1976 distance_min = Float.T(
1977 help='Minimum source-receiver surface distance [m].')
1979 distance_max = Float.T(
1980 help='Maximum source-receiver surface distance [m].')
1982 distance_delta = Float.T(
1983 help='Grid spacing of source-receiver surface distance [m].')
1985 short_type = 'A'
1987 provided_schemes = [
1988 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
1990 def get_surface_distance(self, args):
1991 return args[1]
1993 def get_distance(self, args):
1994 return math.sqrt(args[0]**2 + args[1]**2)
1996 def get_source_depth(self, args):
1997 return args[0]
1999 def get_source_depths(self):
2000 return self.coords[0]
2002 def get_receiver_depth(self, args):
2003 return self.receiver_depth
2005 def _update(self):
2006 self.mins = num.array(
2007 [self.source_depth_min, self.distance_min], dtype=float)
2008 self.maxs = num.array(
2009 [self.source_depth_max, self.distance_max], dtype=float)
2010 self.deltas = num.array(
2011 [self.source_depth_delta, self.distance_delta],
2012 dtype=float)
2013 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2014 vicinity_eps).astype(int) + 1
2015 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2016 self.deltat = 1.0/self.sample_rate
2017 self.nrecords = num.product(self.ns) * self.ncomponents
2018 self.coords = tuple(num.linspace(mi, ma, n) for
2019 (mi, ma, n) in
2020 zip(self.mins, self.effective_maxs, self.ns)) + \
2021 (num.arange(self.ncomponents),)
2023 self.nsource_depths, self.ndistances = self.ns
2025 def _make_index_functions(self):
2027 amin, bmin = self.mins
2028 da, db = self.deltas
2029 na, nb = self.ns
2031 ng = self.ncomponents
2033 def index_function(a, b, ig):
2034 ia = int(round((a - amin) / da))
2035 ib = int(round((b - bmin) / db))
2036 try:
2037 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2038 except ValueError:
2039 raise OutOfBounds()
2041 def indices_function(a, b, ig):
2042 ia = num.round((a - amin) / da).astype(int)
2043 ib = num.round((b - bmin) / db).astype(int)
2044 try:
2045 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))
2046 except ValueError:
2047 for ia_, ib_, ig_ in zip(ia, ib, ig):
2048 try:
2049 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))
2050 except ValueError:
2051 raise OutOfBounds()
2053 def grid_interpolation_coefficients(a, b):
2054 ias = indi12((a - amin) / da, na)
2055 ibs = indi12((b - bmin) / db, nb)
2056 return ias, ibs
2058 def vicinity_function(a, b, ig):
2059 ias, ibs = grid_interpolation_coefficients(a, b)
2061 if not (0 <= ig < ng):
2062 raise OutOfBounds()
2064 indis = []
2065 weights = []
2066 for ia, va in ias:
2067 iia = ia*nb*ng
2068 for ib, vb in ibs:
2069 indis.append(iia + ib*ng + ig)
2070 weights.append(va*vb)
2072 return num.array(indis), num.array(weights)
2074 def vicinities_function(a, b, ig):
2076 xa = (a - amin) / da
2077 xb = (b - bmin) / db
2079 xa_fl = num.floor(xa)
2080 xa_ce = num.ceil(xa)
2081 xb_fl = num.floor(xb)
2082 xb_ce = num.ceil(xb)
2083 va_fl = 1.0 - (xa - xa_fl)
2084 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2085 vb_fl = 1.0 - (xb - xb_fl)
2086 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2088 ia_fl = xa_fl.astype(int)
2089 ia_ce = xa_ce.astype(int)
2090 ib_fl = xb_fl.astype(int)
2091 ib_ce = xb_ce.astype(int)
2093 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2094 raise OutOfBounds()
2096 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2097 raise OutOfBounds()
2099 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2100 raise OutOfBounds()
2102 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2103 raise OutOfBounds()
2105 irecords = num.empty(a.size*4, dtype=int)
2106 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig
2107 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig
2108 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig
2109 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig
2111 weights = num.empty(a.size*4, dtype=float)
2112 weights[0::4] = va_fl * vb_fl
2113 weights[1::4] = va_ce * vb_fl
2114 weights[2::4] = va_fl * vb_ce
2115 weights[3::4] = va_ce * vb_ce
2117 return irecords, weights
2119 self._index_function = index_function
2120 self._indices_function = indices_function
2121 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2122 self._vicinity_function = vicinity_function
2123 self._vicinities_function = vicinities_function
2125 def make_indexing_args(self, source, receiver, icomponents):
2126 nc = icomponents.size
2127 dists = source.distances_to(receiver)
2128 n = dists.size
2129 return (num.tile(source.depths, nc//n),
2130 num.tile(dists, nc//n),
2131 icomponents)
2133 def make_indexing_args1(self, source, receiver):
2134 return (source.depth, source.distance_to(receiver))
2136 @property
2137 def short_extent(self):
2138 return '%g:%g:%g x %g:%g:%g' % (
2139 self.source_depth_min/km,
2140 self.source_depth_max/km,
2141 self.source_depth_delta/km,
2142 self.distance_min/km,
2143 self.distance_max/km,
2144 self.distance_delta/km)
2146 def fix_ttt_holes(self, sptree, mode):
2147 from pyrocko import eikonal_ext, spit
2149 nodes = self.nodes(level=-1)
2151 delta = self.deltas[-1]
2152 assert num.all(delta == self.deltas)
2154 nsources, ndistances = self.ns
2156 points = num.zeros((nodes.shape[0], 3))
2157 points[:, 0] = nodes[:, 1]
2158 points[:, 2] = nodes[:, 0]
2160 speeds = self.get_material_property(
2161 0., 0., points,
2162 parameter='vp' if mode == cake.P else 'vs',
2163 interpolation='multilinear')
2165 speeds = speeds.reshape((nsources, ndistances))
2167 times = sptree.interpolate_many(nodes)
2169 times[num.isnan(times)] = -1.
2170 times = times.reshape(speeds.shape)
2172 try:
2173 eikonal_ext.eikonal_solver_fmm_cartesian(
2174 speeds, times, delta)
2175 except eikonal_ext.EikonalExtError as e:
2176 if str(e).endswith('please check results'):
2177 logger.debug(
2178 'Got a warning from eikonal solver '
2179 '- may be ok...')
2180 else:
2181 raise
2183 def func(x):
2184 ibs, ics = \
2185 self.grid_interpolation_coefficients(*x)
2187 t = 0
2188 for ib, vb in ibs:
2189 for ic, vc in ics:
2190 t += times[ib, ic] * vb * vc
2192 return t
2194 return spit.SPTree(
2195 f=func,
2196 ftol=sptree.ftol,
2197 xbounds=sptree.xbounds,
2198 xtols=sptree.xtols)
2201class ConfigTypeB(Config):
2202 '''
2203 Cylindrical symmetry, 1D earth model, variable receiver depth
2205 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for
2206 receiver depth
2208 * High level index variables: ``(receiver_depth, source_depth,
2209 receiver_distance, component)``
2210 '''
2212 receiver_depth_min = Float.T(
2213 help='Minimum receiver depth [m].')
2215 receiver_depth_max = Float.T(
2216 help='Maximum receiver depth [m].')
2218 receiver_depth_delta = Float.T(
2219 help='Grid spacing of receiver depths [m]')
2221 source_depth_min = Float.T(
2222 help='Minimum source depth [m].')
2224 source_depth_max = Float.T(
2225 help='Maximum source depth [m].')
2227 source_depth_delta = Float.T(
2228 help='Grid spacing of source depths [m]')
2230 distance_min = Float.T(
2231 help='Minimum source-receiver surface distance [m].')
2233 distance_max = Float.T(
2234 help='Maximum source-receiver surface distance [m].')
2236 distance_delta = Float.T(
2237 help='Grid spacing of source-receiver surface distances [m].')
2239 short_type = 'B'
2241 provided_schemes = [
2242 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']
2244 def get_distance(self, args):
2245 return math.sqrt((args[1] - args[0])**2 + args[2]**2)
2247 def get_surface_distance(self, args):
2248 return args[2]
2250 def get_source_depth(self, args):
2251 return args[1]
2253 def get_receiver_depth(self, args):
2254 return args[0]
2256 def get_source_depths(self):
2257 return self.coords[1]
2259 def _update(self):
2260 self.mins = num.array([
2261 self.receiver_depth_min,
2262 self.source_depth_min,
2263 self.distance_min],
2264 dtype=float)
2266 self.maxs = num.array([
2267 self.receiver_depth_max,
2268 self.source_depth_max,
2269 self.distance_max],
2270 dtype=float)
2272 self.deltas = num.array([
2273 self.receiver_depth_delta,
2274 self.source_depth_delta,
2275 self.distance_delta],
2276 dtype=float)
2278 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2279 vicinity_eps).astype(int) + 1
2280 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2281 self.deltat = 1.0/self.sample_rate
2282 self.nrecords = num.product(self.ns) * self.ncomponents
2283 self.coords = tuple(num.linspace(mi, ma, n) for
2284 (mi, ma, n) in
2285 zip(self.mins, self.effective_maxs, self.ns)) + \
2286 (num.arange(self.ncomponents),)
2287 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2289 def _make_index_functions(self):
2291 amin, bmin, cmin = self.mins
2292 da, db, dc = self.deltas
2293 na, nb, nc = self.ns
2294 ng = self.ncomponents
2296 def index_function(a, b, c, ig):
2297 ia = int(round((a - amin) / da))
2298 ib = int(round((b - bmin) / db))
2299 ic = int(round((c - cmin) / dc))
2300 try:
2301 return num.ravel_multi_index((ia, ib, ic, ig),
2302 (na, nb, nc, ng))
2303 except ValueError:
2304 raise OutOfBounds()
2306 def indices_function(a, b, c, ig):
2307 ia = num.round((a - amin) / da).astype(int)
2308 ib = num.round((b - bmin) / db).astype(int)
2309 ic = num.round((c - cmin) / dc).astype(int)
2310 try:
2311 return num.ravel_multi_index((ia, ib, ic, ig),
2312 (na, nb, nc, ng))
2313 except ValueError:
2314 raise OutOfBounds()
2316 def grid_interpolation_coefficients(a, b, c):
2317 ias = indi12((a - amin) / da, na)
2318 ibs = indi12((b - bmin) / db, nb)
2319 ics = indi12((c - cmin) / dc, nc)
2320 return ias, ibs, ics
2322 def vicinity_function(a, b, c, ig):
2323 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)
2325 if not (0 <= ig < ng):
2326 raise OutOfBounds()
2328 indis = []
2329 weights = []
2330 for ia, va in ias:
2331 iia = ia*nb*nc*ng
2332 for ib, vb in ibs:
2333 iib = ib*nc*ng
2334 for ic, vc in ics:
2335 indis.append(iia + iib + ic*ng + ig)
2336 weights.append(va*vb*vc)
2338 return num.array(indis), num.array(weights)
2340 def vicinities_function(a, b, c, ig):
2342 xa = (a - amin) / da
2343 xb = (b - bmin) / db
2344 xc = (c - cmin) / dc
2346 xa_fl = num.floor(xa)
2347 xa_ce = num.ceil(xa)
2348 xb_fl = num.floor(xb)
2349 xb_ce = num.ceil(xb)
2350 xc_fl = num.floor(xc)
2351 xc_ce = num.ceil(xc)
2352 va_fl = 1.0 - (xa - xa_fl)
2353 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2354 vb_fl = 1.0 - (xb - xb_fl)
2355 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2356 vc_fl = 1.0 - (xc - xc_fl)
2357 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2359 ia_fl = xa_fl.astype(int)
2360 ia_ce = xa_ce.astype(int)
2361 ib_fl = xb_fl.astype(int)
2362 ib_ce = xb_ce.astype(int)
2363 ic_fl = xc_fl.astype(int)
2364 ic_ce = xc_ce.astype(int)
2366 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2367 raise OutOfBounds()
2369 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2370 raise OutOfBounds()
2372 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2373 raise OutOfBounds()
2375 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2376 raise OutOfBounds()
2378 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2379 raise OutOfBounds()
2381 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2382 raise OutOfBounds()
2384 irecords = num.empty(a.size*8, dtype=int)
2385 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2386 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig
2387 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2388 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig
2389 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2390 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig
2391 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2392 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig
2394 weights = num.empty(a.size*8, dtype=float)
2395 weights[0::8] = va_fl * vb_fl * vc_fl
2396 weights[1::8] = va_ce * vb_fl * vc_fl
2397 weights[2::8] = va_fl * vb_ce * vc_fl
2398 weights[3::8] = va_ce * vb_ce * vc_fl
2399 weights[4::8] = va_fl * vb_fl * vc_ce
2400 weights[5::8] = va_ce * vb_fl * vc_ce
2401 weights[6::8] = va_fl * vb_ce * vc_ce
2402 weights[7::8] = va_ce * vb_ce * vc_ce
2404 return irecords, weights
2406 self._index_function = index_function
2407 self._indices_function = indices_function
2408 self._grid_interpolation_coefficients = grid_interpolation_coefficients
2409 self._vicinity_function = vicinity_function
2410 self._vicinities_function = vicinities_function
2412 def make_indexing_args(self, source, receiver, icomponents):
2413 nc = icomponents.size
2414 dists = source.distances_to(receiver)
2415 n = dists.size
2416 receiver_depths = num.empty(nc)
2417 receiver_depths.fill(receiver.depth)
2418 return (receiver_depths,
2419 num.tile(source.depths, nc//n),
2420 num.tile(dists, nc//n),
2421 icomponents)
2423 def make_indexing_args1(self, source, receiver):
2424 return (receiver.depth,
2425 source.depth,
2426 source.distance_to(receiver))
2428 @property
2429 def short_extent(self):
2430 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2431 self.receiver_depth_min/km,
2432 self.receiver_depth_max/km,
2433 self.receiver_depth_delta/km,
2434 self.source_depth_min/km,
2435 self.source_depth_max/km,
2436 self.source_depth_delta/km,
2437 self.distance_min/km,
2438 self.distance_max/km,
2439 self.distance_delta/km)
2441 def fix_ttt_holes(self, sptree, mode):
2442 from pyrocko import eikonal_ext, spit
2444 nodes_sr = self.nodes(minlevel=1, level=-1)
2446 delta = self.deltas[-1]
2447 assert num.all(delta == self.deltas[1:])
2449 nreceivers, nsources, ndistances = self.ns
2451 points = num.zeros((nodes_sr.shape[0], 3))
2452 points[:, 0] = nodes_sr[:, 1]
2453 points[:, 2] = nodes_sr[:, 0]
2455 speeds = self.get_material_property(
2456 0., 0., points,
2457 parameter='vp' if mode == cake.P else 'vs',
2458 interpolation='multilinear')
2460 speeds = speeds.reshape((nsources, ndistances))
2462 receiver_times = []
2463 for ireceiver in range(nreceivers):
2464 nodes = num.hstack([
2465 num_full(
2466 (nodes_sr.shape[0], 1),
2467 self.coords[0][ireceiver]),
2468 nodes_sr])
2470 times = sptree.interpolate_many(nodes)
2472 times[num.isnan(times)] = -1.
2474 times = times.reshape(speeds.shape)
2476 try:
2477 eikonal_ext.eikonal_solver_fmm_cartesian(
2478 speeds, times, delta)
2479 except eikonal_ext.EikonalExtError as e:
2480 if str(e).endswith('please check results'):
2481 logger.debug(
2482 'Got a warning from eikonal solver '
2483 '- may be ok...')
2484 else:
2485 raise
2487 receiver_times.append(times)
2489 def func(x):
2490 ias, ibs, ics = \
2491 self.grid_interpolation_coefficients(*x)
2493 t = 0
2494 for ia, va in ias:
2495 times = receiver_times[ia]
2496 for ib, vb in ibs:
2497 for ic, vc in ics:
2498 t += times[ib, ic] * va * vb * vc
2500 return t
2502 return spit.SPTree(
2503 f=func,
2504 ftol=sptree.ftol,
2505 xbounds=sptree.xbounds,
2506 xtols=sptree.xtols)
2509class ConfigTypeC(Config):
2510 '''
2511 No symmetrical constraints but fixed receiver positions.
2513 * Cartesian 3D source volume around a reference point
2515 * High level index variables: ``(ireceiver, source_depth,
2516 source_east_shift, source_north_shift, component)``
2517 '''
2519 receivers = List.T(
2520 Receiver.T(),
2521 help='List of fixed receivers.')
2523 source_origin = Location.T(
2524 help='Origin of the source volume grid.')
2526 source_depth_min = Float.T(
2527 help='Minimum source depth [m].')
2529 source_depth_max = Float.T(
2530 help='Maximum source depth [m].')
2532 source_depth_delta = Float.T(
2533 help='Source depth grid spacing [m].')
2535 source_east_shift_min = Float.T(
2536 help='Minimum easting of source grid [m].')
2538 source_east_shift_max = Float.T(
2539 help='Maximum easting of source grid [m].')
2541 source_east_shift_delta = Float.T(
2542 help='Source volume grid spacing in east direction [m].')
2544 source_north_shift_min = Float.T(
2545 help='Minimum northing of source grid [m].')
2547 source_north_shift_max = Float.T(
2548 help='Maximum northing of source grid [m].')
2550 source_north_shift_delta = Float.T(
2551 help='Source volume grid spacing in north direction [m].')
2553 short_type = 'C'
2555 provided_schemes = ['elastic18']
2557 def get_surface_distance(self, args):
2558 ireceiver, _, source_east_shift, source_north_shift, _ = args
2559 sorig = self.source_origin
2560 sloc = Location(
2561 lat=sorig.lat,
2562 lon=sorig.lon,
2563 north_shift=sorig.north_shift + source_north_shift,
2564 east_shift=sorig.east_shift + source_east_shift)
2566 return self.receivers[args[0]].distance_to(sloc)
2568 def get_distance(self, args):
2569 # to be improved...
2570 ireceiver, sdepth, source_east_shift, source_north_shift, _ = args
2571 sorig = self.source_origin
2572 sloc = Location(
2573 lat=sorig.lat,
2574 lon=sorig.lon,
2575 north_shift=sorig.north_shift + source_north_shift,
2576 east_shift=sorig.east_shift + source_east_shift)
2578 return math.sqrt(
2579 self.receivers[args[0]].distance_to(sloc)**2 + sdepth**2)
2581 def get_source_depth(self, args):
2582 return args[1]
2584 def get_receiver_depth(self, args):
2585 return self.receivers[args[0]].depth
2587 def get_source_depths(self):
2588 return self.coords[0]
2590 def _update(self):
2591 self.mins = num.array([
2592 self.source_depth_min,
2593 self.source_east_shift_min,
2594 self.source_north_shift_min],
2595 dtype=float)
2597 self.maxs = num.array([
2598 self.source_depth_max,
2599 self.source_east_shift_max,
2600 self.source_north_shift_max],
2601 dtype=float)
2603 self.deltas = num.array([
2604 self.source_depth_delta,
2605 self.source_east_shift_delta,
2606 self.source_north_shift_delta],
2607 dtype=float)
2609 self.ns = num.floor((self.maxs - self.mins) / self.deltas +
2610 vicinity_eps).astype(int) + 1
2611 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)
2612 self.deltat = 1.0/self.sample_rate
2613 self.nreceivers = len(self.receivers)
2614 self.nrecords = \
2615 self.nreceivers * num.product(self.ns) * self.ncomponents
2617 self.coords = (num.arange(self.nreceivers),) + \
2618 tuple(num.linspace(mi, ma, n) for (mi, ma, n) in
2619 zip(self.mins, self.effective_maxs, self.ns)) + \
2620 (num.arange(self.ncomponents),)
2621 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns
2623 self._distances_cache = {}
2625 def _make_index_functions(self):
2627 amin, bmin, cmin = self.mins
2628 da, db, dc = self.deltas
2629 na, nb, nc = self.ns
2630 ng = self.ncomponents
2631 nr = self.nreceivers
2633 def index_function(ir, a, b, c, ig):
2634 ia = int(round((a - amin) / da))
2635 ib = int(round((b - bmin) / db))
2636 ic = int(round((c - cmin) / dc))
2637 try:
2638 return num.ravel_multi_index((ir, ia, ib, ic, ig),
2639 (nr, na, nb, nc, ng))
2640 except ValueError:
2641 raise OutOfBounds()
2643 def indices_function(ir, a, b, c, ig):
2644 ia = num.round((a - amin) / da).astype(int)
2645 ib = num.round((b - bmin) / db).astype(int)
2646 ic = num.round((c - cmin) / dc).astype(int)
2648 try:
2649 return num.ravel_multi_index((ir, ia, ib, ic, ig),
2650 (nr, na, nb, nc, ng))
2651 except ValueError:
2652 raise OutOfBounds()
2654 def vicinity_function(ir, a, b, c, ig):
2655 ias = indi12((a - amin) / da, na)
2656 ibs = indi12((b - bmin) / db, nb)
2657 ics = indi12((c - cmin) / dc, nc)
2659 if not (0 <= ir < nr):
2660 raise OutOfBounds()
2662 if not (0 <= ig < ng):
2663 raise OutOfBounds()
2665 indis = []
2666 weights = []
2667 iir = ir*na*nb*nc*ng
2668 for ia, va in ias:
2669 iia = ia*nb*nc*ng
2670 for ib, vb in ibs:
2671 iib = ib*nc*ng
2672 for ic, vc in ics:
2673 indis.append(iir + iia + iib + ic*ng + ig)
2674 weights.append(va*vb*vc)
2676 return num.array(indis), num.array(weights)
2678 def vicinities_function(ir, a, b, c, ig):
2680 xa = (a-amin) / da
2681 xb = (b-bmin) / db
2682 xc = (c-cmin) / dc
2684 xa_fl = num.floor(xa)
2685 xa_ce = num.ceil(xa)
2686 xb_fl = num.floor(xb)
2687 xb_ce = num.ceil(xb)
2688 xc_fl = num.floor(xc)
2689 xc_ce = num.ceil(xc)
2690 va_fl = 1.0 - (xa - xa_fl)
2691 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)
2692 vb_fl = 1.0 - (xb - xb_fl)
2693 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)
2694 vc_fl = 1.0 - (xc - xc_fl)
2695 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)
2697 ia_fl = xa_fl.astype(int)
2698 ia_ce = xa_ce.astype(int)
2699 ib_fl = xb_fl.astype(int)
2700 ib_ce = xb_ce.astype(int)
2701 ic_fl = xc_fl.astype(int)
2702 ic_ce = xc_ce.astype(int)
2704 if num.any(ia_fl < 0) or num.any(ia_fl >= na):
2705 raise OutOfBounds()
2707 if num.any(ia_ce < 0) or num.any(ia_ce >= na):
2708 raise OutOfBounds()
2710 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):
2711 raise OutOfBounds()
2713 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):
2714 raise OutOfBounds()
2716 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):
2717 raise OutOfBounds()
2719 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):
2720 raise OutOfBounds()
2722 irig = ir*na*nb*nc*ng + ig
2724 irecords = num.empty(a.size*8, dtype=int)
2725 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + irig
2726 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + irig
2727 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + irig
2728 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + irig
2729 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + irig
2730 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + irig
2731 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + irig
2732 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + irig
2734 weights = num.empty(a.size*8, dtype=float)
2735 weights[0::8] = va_fl * vb_fl * vc_fl
2736 weights[1::8] = va_ce * vb_fl * vc_fl
2737 weights[2::8] = va_fl * vb_ce * vc_fl
2738 weights[3::8] = va_ce * vb_ce * vc_fl
2739 weights[4::8] = va_fl * vb_fl * vc_ce
2740 weights[5::8] = va_ce * vb_fl * vc_ce
2741 weights[6::8] = va_fl * vb_ce * vc_ce
2742 weights[7::8] = va_ce * vb_ce * vc_ce
2744 return irecords, weights
2746 self._index_function = index_function
2747 self._indices_function = indices_function
2748 self._vicinity_function = vicinity_function
2749 self._vicinities_function = vicinities_function
2751 def lookup_ireceiver(self, receiver):
2752 k = (receiver.lat, receiver.lon,
2753 receiver.north_shift, receiver.east_shift)
2754 dh = min(self.source_north_shift_delta, self.source_east_shift_delta)
2755 dv = self.source_depth_delta
2757 for irec, rec in enumerate(self.receivers):
2758 if (k, irec) not in self._distances_cache:
2759 self._distances_cache[k, irec] = math.sqrt(
2760 (receiver.distance_to(rec)/dh)**2 +
2761 ((rec.depth - receiver.depth)/dv)**2)
2763 if self._distances_cache[k, irec] < 0.1:
2764 return irec
2766 raise OutOfBounds(
2767 reason='No GFs available for receiver at (%g, %g).' %
2768 receiver.effective_latlon)
2770 def make_indexing_args(self, source, receiver, icomponents):
2771 nc = icomponents.size
2773 dists = source.distances_to(self.source_origin)
2774 azis, _ = source.azibazis_to(self.source_origin)
2776 source_north_shifts = - num.cos(d2r*azis) * dists
2777 source_east_shifts = - num.sin(d2r*azis) * dists
2778 source_depths = source.depths - self.source_origin.depth
2780 n = dists.size
2781 ireceivers = num.empty(nc, dtype=int)
2782 ireceivers.fill(self.lookup_ireceiver(receiver))
2784 return (ireceivers,
2785 num.tile(source_depths, nc//n),
2786 num.tile(source_east_shifts, nc//n),
2787 num.tile(source_north_shifts, nc//n),
2788 icomponents)
2790 def make_indexing_args1(self, source, receiver):
2791 dist = source.distance_to(self.source_origin)
2792 azi, _ = source.azibazi_to(self.source_origin)
2794 source_north_shift = - num.cos(d2r*azi) * dist
2795 source_east_shift = - num.sin(d2r*azi) * dist
2796 source_depth = source.depth - self.source_origin.depth
2798 return (self.lookup_ireceiver(receiver),
2799 source_depth,
2800 source_east_shift,
2801 source_north_shift)
2803 @property
2804 def short_extent(self):
2805 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (
2806 self.source_depth_min/km,
2807 self.source_depth_max/km,
2808 self.source_depth_delta/km,
2809 self.source_east_shift_min/km,
2810 self.source_east_shift_max/km,
2811 self.source_east_shift_delta/km,
2812 self.source_north_shift_min/km,
2813 self.source_north_shift_max/km,
2814 self.source_north_shift_delta/km)
2817class Weighting(Object):
2818 factor = Float.T(default=1.0)
2821class Taper(Object):
2822 tmin = Timing.T()
2823 tmax = Timing.T()
2824 tfade = Float.T(default=0.0)
2825 shape = StringChoice.T(
2826 choices=['cos', 'linear'],
2827 default='cos',
2828 optional=True)
2831class SimplePattern(SObject):
2833 _pool = {}
2835 def __init__(self, pattern):
2836 self._pattern = pattern
2837 SObject.__init__(self)
2839 def __str__(self):
2840 return self._pattern
2842 @property
2843 def regex(self):
2844 pool = SimplePattern._pool
2845 if self.pattern not in pool:
2846 rpat = '|'.join(fnmatch.translate(x) for
2847 x in self.pattern.split('|'))
2848 pool[self.pattern] = re.compile(rpat, re.I)
2850 return pool[self.pattern]
2852 def match(self, s):
2853 return self.regex.match(s)
2856class WaveformType(StringChoice):
2857 choices = ['dis', 'vel', 'acc',
2858 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',
2859 'envelope_dis', 'envelope_vel', 'envelope_acc']
2862class ChannelSelection(Object):
2863 pattern = SimplePattern.T(optional=True)
2864 min_sample_rate = Float.T(optional=True)
2865 max_sample_rate = Float.T(optional=True)
2868class StationSelection(Object):
2869 includes = SimplePattern.T()
2870 excludes = SimplePattern.T()
2871 distance_min = Float.T(optional=True)
2872 distance_max = Float.T(optional=True)
2873 azimuth_min = Float.T(optional=True)
2874 azimuth_max = Float.T(optional=True)
2877class WaveformSelection(Object):
2878 channel_selection = ChannelSelection.T(optional=True)
2879 station_selection = StationSelection.T(optional=True)
2880 taper = Taper.T()
2881 # filter = FrequencyResponse.T()
2882 waveform_type = WaveformType.T(default='dis')
2883 weighting = Weighting.T(optional=True)
2884 sample_rate = Float.T(optional=True)
2885 gf_store_id = StringID.T(optional=True)
2888def indi12(x, n):
2889 '''
2890 Get linear interpolation index and weight.
2891 '''
2893 r = round(x)
2894 if abs(r - x) < vicinity_eps:
2895 i = int(r)
2896 if not (0 <= i < n):
2897 raise OutOfBounds()
2899 return ((int(r), 1.),)
2900 else:
2901 f = math.floor(x)
2902 i = int(f)
2903 if not (0 <= i < n-1):
2904 raise OutOfBounds()
2906 v = x-f
2907 return ((i, 1.-v), (i + 1, v))
2910def float_or_none(s):
2911 units = {
2912 'k': 1e3,
2913 'M': 1e6,
2914 }
2916 factor = 1.0
2917 if s and s[-1] in units:
2918 factor = units[s[-1]]
2919 s = s[:-1]
2920 if not s:
2921 raise ValueError('unit without a number: \'%s\'' % s)
2923 if s:
2924 return float(s) * factor
2925 else:
2926 return None
2929class GridSpecError(Exception):
2930 def __init__(self, s):
2931 Exception.__init__(self, 'invalid grid specification: %s' % s)
2934def parse_grid_spec(spec):
2935 try:
2936 result = []
2937 for dspec in spec.split(','):
2938 t = dspec.split('@')
2939 num = start = stop = step = None
2940 if len(t) == 2:
2941 num = int(t[1])
2942 if num <= 0:
2943 raise GridSpecError(spec)
2945 elif len(t) > 2:
2946 raise GridSpecError(spec)
2948 s = t[0]
2949 v = [float_or_none(x) for x in s.split(':')]
2950 if len(v) == 1:
2951 start = stop = v[0]
2952 if len(v) >= 2:
2953 start, stop = v[0:2]
2954 if len(v) == 3:
2955 step = v[2]
2957 if len(v) > 3 or (len(v) > 2 and num is not None):
2958 raise GridSpecError(spec)
2960 if step == 0.0:
2961 raise GridSpecError(spec)
2963 result.append((start, stop, step, num))
2965 except ValueError:
2966 raise GridSpecError(spec)
2968 return result
2971def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):
2972 swap = step is not None and step < 0.
2973 if start is None:
2974 start = [mi, ma][swap]
2975 if stop is None:
2976 stop = [ma, mi][swap]
2977 if step is None:
2978 step = [inc, -inc][ma < mi]
2979 if num is None:
2980 if (step < 0) != (stop-start < 0):
2981 raise GridSpecError()
2983 num = int(round((stop-start)/step))+1
2984 stop2 = start + (num-1)*step
2985 if abs(stop-stop2) > eps:
2986 num = int(math.floor((stop-start)/step))+1
2987 stop = start + (num-1)*step
2988 else:
2989 stop = stop2
2991 if start == stop:
2992 num = 1
2994 return start, stop, num
2997def nditer_outer(x):
2998 return num.nditer(
2999 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())
3002def nodes(xs):
3003 ns = [x.size for x in xs]
3004 nnodes = num.prod(ns)
3005 ndim = len(xs)
3006 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)
3007 for idim in range(ndim-1, -1, -1):
3008 x = xs[idim]
3009 nrepeat = num.prod(ns[idim+1:], dtype=int)
3010 ntile = num.prod(ns[:idim], dtype=int)
3011 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)
3013 return nodes
3016def filledi(x, n):
3017 a = num.empty(n, dtype=int)
3018 a.fill(x)
3019 return a
3022config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]
3024discretized_source_classes = [
3025 DiscretizedExplosionSource,
3026 DiscretizedSFSource,
3027 DiscretizedMTSource,
3028 DiscretizedPorePressureSource]
3031__all__ = '''
3032Earthmodel1D
3033StringID
3034ScopeType
3035WaveformType
3036QuantityType
3037NearfieldTermsType
3038Reference
3039Region
3040CircularRegion
3041RectangularRegion
3042PhaseSelect
3043InvalidTimingSpecification
3044Timing
3045TPDef
3046OutOfBounds
3047Location
3048Receiver
3049'''.split() + [
3050 S.__name__ for S in discretized_source_classes + config_type_classes] + '''
3051ComponentScheme
3052component_scheme_to_description
3053component_schemes
3054Config
3055GridSpecError
3056Weighting
3057Taper
3058SimplePattern
3059WaveformType
3060ChannelSelection
3061StationSelection
3062WaveformSelection
3063nditer_outer
3064dump
3065load
3066discretized_source_classes
3067config_type_classes
3068UnavailableScheme
3069InterpolationMethod
3070SeismosizerTrace
3071SeismosizerResult
3072Result
3073StaticResult
3074'''.split()