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