Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/meta.py: 74%

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

4# ---|P------/S----------~Lg----------

6'''

7Data model for Green's function store meta-data.

8'''

10import math

11import re

12import fnmatch

13import logging

15import numpy as num

16from scipy.interpolate import interp1d

18from pyrocko.guts import (Object, SObject, String, StringChoice,

19 StringPattern, Unicode, Float, Bool, Int, TBase,

20 List, ValidationError, Timestamp, Tuple, Dict)

21from pyrocko.guts import dump, load # noqa

22from pyrocko.guts_array import literal, Array

23from pyrocko.model import Location, gnss

25from pyrocko import cake, orthodrome, spit, moment_tensor, trace

26from pyrocko.util import num_full

28from .error import StoreError

30try:

31 new_str = unicode

32except NameError:

33 new_str = str

35guts_prefix = 'pf'

37d2r = math.pi / 180.

38r2d = 1.0 / d2r

39km = 1000.

40vicinity_eps = 1e-5

42logger = logging.getLogger('pyrocko.gf.meta')

45def fmt_choices(cls):

46 return 'choices: %s' % ', '.join(

47 "``'%s'``" % choice for choice in cls.choices)

50class InterpolationMethod(StringChoice):

51 choices = ['nearest_neighbor', 'multilinear']

54class SeismosizerTrace(Object):

56 codes = Tuple.T(

57 4, String.T(),

58 default=('', 'STA', '', 'Z'),

59 help='network, station, location and channel codes')

61 data = Array.T(

62 shape=(None,),

63 dtype=num.float32,

64 serialize_as='base64',

65 serialize_dtype=num.dtype('<f4'),

66 help='numpy array with data samples')

68 deltat = Float.T(

69 default=1.0,

70 help='sampling interval [s]')

72 tmin = Timestamp.T(

73 default=Timestamp.D('1970-01-01 00:00:00'),

74 help='time of first sample as a system timestamp [s]')

76 def pyrocko_trace(self):

77 c = self.codes

78 return trace.Trace(c[0], c[1], c[2], c[3],

79 ydata=self.data,

80 deltat=self.deltat,

81 tmin=self.tmin)

83 @classmethod

84 def from_pyrocko_trace(cls, tr, **kwargs):

85 d = dict(

86 codes=tr.nslc_id,

87 tmin=tr.tmin,

88 deltat=tr.deltat,

89 data=num.asarray(tr.get_ydata(), dtype=num.float32))

90 d.update(kwargs)

91 return cls(**d)

94class SeismosizerResult(Object):

95 n_records_stacked = Int.T(optional=True, default=1)

96 t_stack = Float.T(optional=True, default=0.)

99class Result(SeismosizerResult):

100 trace = SeismosizerTrace.T(optional=True)

101 n_shared_stacking = Int.T(optional=True, default=1)

102 t_optimize = Float.T(optional=True, default=0.)

103

104

105class StaticResult(SeismosizerResult):

106 result = Dict.T(

107 String.T(),

108 Array.T(shape=(None,), dtype=float, serialize_as='base64'))

109

110

111class GNSSCampaignResult(StaticResult):

112 campaign = gnss.GNSSCampaign.T(

113 optional=True)

114

115

116class SatelliteResult(StaticResult):

117

118 theta = Array.T(

119 optional=True,

120 shape=(None,), dtype=float, serialize_as='base64')

121

122 phi = Array.T(

123 optional=True,

124 shape=(None,), dtype=float, serialize_as='base64')

125

126

127class KiteSceneResult(SatelliteResult):

128

129 shape = Tuple.T()

130

131 def get_scene(self, component='displacement.los'):

132 try:

133 from kite import Scene

134 except ImportError:

135 raise ImportError('Kite not installed')

136

137 def reshape(arr):

138 return arr.reshape(self.shape)

139

140 displacement = self.result[component]

141

142 displacement, theta, phi = map(

143 reshape, (displacement, self.theta, self.phi))

144

145 sc = Scene(

146 displacement=displacement,

147 theta=theta, phi=phi,

148 config=self.config)

149

150 return sc

151

152

153class ComponentSchemeDescription(Object):

154 name = String.T()

155 source_terms = List.T(String.T())

156 ncomponents = Int.T()

157 default_stored_quantity = String.T(optional=True)

158 provided_components = List.T(String.T())

159

160

161component_scheme_descriptions = [

162 ComponentSchemeDescription(

163 name='elastic2',

164 source_terms=['m0'],

165 ncomponents=2,

166 default_stored_quantity='displacement',

167 provided_components=[

168 'n', 'e', 'd']),

169 ComponentSchemeDescription(

170 name='elastic8',

171 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],

172 ncomponents=8,

173 default_stored_quantity='displacement',

174 provided_components=[

175 'n', 'e', 'd']),

176 ComponentSchemeDescription(

177 name='elastic10',

178 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],

179 ncomponents=10,

180 default_stored_quantity='displacement',

181 provided_components=[

182 'n', 'e', 'd']),

183 ComponentSchemeDescription(

184 name='elastic18',

185 source_terms=['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med'],

186 ncomponents=18,

187 default_stored_quantity='displacement',

188 provided_components=[

189 'n', 'e', 'd']),

190 ComponentSchemeDescription(

191 name='elastic5',

192 source_terms=['fn', 'fe', 'fd'],

193 ncomponents=5,

194 default_stored_quantity='displacement',

195 provided_components=[

196 'n', 'e', 'd']),

197 ComponentSchemeDescription(

198 name='poroelastic10',

199 source_terms=['pore_pressure'],

200 ncomponents=10,

201 default_stored_quantity=None,

202 provided_components=[

203 'displacement.n', 'displacement.e', 'displacement.d',

204 'vertical_tilt.n', 'vertical_tilt.e',

205 'pore_pressure',

206 'darcy_velocity.n', 'darcy_velocity.e', 'darcy_velocity.d'])]

207

208

209# new names?

210# 'mt_to_displacement_1d'

211# 'mt_to_displacement_1d_ff_only'

212# 'mt_to_gravity_1d'

213# 'mt_to_stress_1d'

214# 'explosion_to_displacement_1d'

215# 'sf_to_displacement_1d'

216# 'mt_to_displacement_3d'

217# 'mt_to_gravity_3d'

218# 'mt_to_stress_3d'

219# 'pore_pressure_to_displacement_1d'

220# 'pore_pressure_to_vertical_tilt_1d'

221# 'pore_pressure_to_pore_pressure_1d'

222# 'pore_pressure_to_darcy_velocity_1d'

223

224

225component_schemes = [c.name for c in component_scheme_descriptions]

226component_scheme_to_description = dict(

227 (c.name, c) for c in component_scheme_descriptions)

228

229

230class ComponentScheme(StringChoice):

231 '''

232 Different Green's Function component schemes are available:

233

234 ================= =========================================================

235 Name Description

236 ================= =========================================================

237 ``elastic10`` Elastodynamic for

238 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and

239 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, MT

240 sources only

241 ``elastic8`` Elastodynamic for far-field only

242 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and

243 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores,

244 MT sources only

245 ``elastic2`` Elastodynamic for

246 :py:class:`~pyrocko.gf.meta.ConfigTypeA` and

247 :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, purely

248 isotropic sources only

249 ``elastic5`` Elastodynamic for

250 :py:class:`~pyrocko.gf.meta.ConfigTypeA`

251 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores, SF

252 sources only

253 ``elastic18`` Elastodynamic for

254 :py:class:`~pyrocko.gf.meta.ConfigTypeC` stores, MT

255 sources only

256 ``poroelastic10`` Poroelastic for :py:class:`~pyrocko.gf.meta.ConfigTypeA`

257 and :py:class:`~pyrocko.gf.meta.ConfigTypeB` stores

258 ================= =========================================================

259 '''

260

261 choices = component_schemes

262

263

264class Earthmodel1D(Object):

265 dummy_for = cake.LayeredModel

266 dummy_for_description = 'pyrocko.cake.LayeredModel'

267

268 class __T(TBase):

269 def regularize_extra(self, val):

270 if isinstance(val, str):

271 val = cake.LayeredModel.from_scanlines(

272 cake.read_nd_model_str(val))

273

274 return val

275

276 def to_save(self, val):

277 return literal(cake.write_nd_model_str(val))

278

279

280class StringID(StringPattern):

281 pattern = r'^[A-Za-z][A-Za-z0-9._]{0,64}$'

282

283

284class ScopeType(StringChoice):

285 choices = [

286 'global',

287 'regional',

288 'local',

289 ]

290

291

292class WaveType(StringChoice):

293 choices = [

294 'full waveform',

295 'bodywave',

296 'P wave',

297 'S wave',

298 'surface wave',

299 ]

300

301

302class NearfieldTermsType(StringChoice):

303 choices = [

304 'complete',

305 'incomplete',

306 'missing',

307 ]

308

309

310class QuantityType(StringChoice):

311 choices = [

312 'displacement',

313 'rotation',

314 'velocity',

315 'acceleration',

316 'pressure',

317 'tilt',

318 'pore_pressure',

319 'darcy_velocity',

320 'vertical_tilt']

321

322

323class Reference(Object):

324 id = StringID.T()

325 type = String.T()

326 title = Unicode.T()

327 journal = Unicode.T(optional=True)

328 volume = Unicode.T(optional=True)

329 number = Unicode.T(optional=True)

330 pages = Unicode.T(optional=True)

331 year = String.T()

332 note = Unicode.T(optional=True)

333 issn = String.T(optional=True)

334 doi = String.T(optional=True)

335 url = String.T(optional=True)

336 eprint = String.T(optional=True)

337 authors = List.T(Unicode.T())

338 publisher = Unicode.T(optional=True)

339 keywords = Unicode.T(optional=True)

340 note = Unicode.T(optional=True)

341 abstract = Unicode.T(optional=True)

342

343 @classmethod

344 def from_bibtex(cls, filename=None, stream=None):

345

346 from pybtex.database.input import bibtex

347

348 parser = bibtex.Parser()

349

350 if filename is not None:

351 bib_data = parser.parse_file(filename)

352 elif stream is not None:

353 bib_data = parser.parse_stream(stream)

354

355 references = []

356

357 for id_, entry in bib_data.entries.items():

358 d = {}

359 avail = entry.fields.keys()

360 for prop in cls.T.properties:

361 if prop.name == 'authors' or prop.name not in avail:

362 continue

363

364 d[prop.name] = entry.fields[prop.name]

365

366 if 'author' in entry.persons:

367 d['authors'] = []

368 for person in entry.persons['author']:

369 d['authors'].append(new_str(person))

370

371 c = Reference(id=id_, type=entry.type, **d)

372 references.append(c)

373

374 return references

375

376

377_fpat = r'[+-]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?'

378_spat = StringID.pattern[1:-1]

379_cake_pat = cake.PhaseDef.allowed_characters_pattern

380_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic

381

382_ppat = r'(' \

383 r'cake:' + _cake_pat + \

384 r'|iaspei:' + _iaspei_pat + \

385 r'|vel_surface:' + _fpat + \

386 r'|vel:' + _fpat + \

387 r'|stored:' + _spat + \

388 r')'

389

390

391timing_regex_old = re.compile(

392 r'^((first|last)?$(' + _spat + r'(\|' + _spat + r')*)$|(' +

393 _spat + r'))?(' + _fpat + ')?$')

394

395

396timing_regex = re.compile(

397 r'^((first|last)?\{(' + _ppat + r'(\|' + _ppat + r')*)\}|(' +

398 _ppat + r'))?(' + _fpat + '(S|%)?)?$')

399

400

401class PhaseSelect(StringChoice):

402 choices = ['', 'first', 'last']

403

404

405class InvalidTimingSpecification(ValidationError):

406 pass

407

408

409class TimingAttributeError(Exception):

410 pass

411

412

413class Timing(SObject):

414 '''

415 Definition of a time instant relative to one or more named phase arrivals.

416

417 Instances of this class can be used e.g. in cutting and tapering

418 operations. They can hold an absolute time or an offset to a named phase

419 arrival or group of such arrivals.

420

421 Timings can be instantiated from a simple string defintion i.e. with

422 ``Timing(str)`` where ``str`` is something like

423 ``'SELECT{PHASE_DEFS}[+-]OFFSET[S|%]'`` where ``'SELECT'`` is ``'first'``,

424 ``'last'`` or empty, ``'PHASE_DEFS'`` is a ``'|'``-separated list of phase

425 definitions, and ``'OFFSET'`` is the time offset in seconds. If a ``'%'``

426 is appended, it is interpreted as percent. If the an ``'S'`` is appended

427 to ``'OFFSET'``, it is interpreted as a surface slowness in [s/km].

428

429 Phase definitions can be specified in either of the following ways:

430

431 * ``'stored:PHASE_ID'`` - retrieves value from stored travel time table

432 * ``'cake:CAKE_PHASE_DEF'`` - evaluates first arrival of phase with cake

433 (see :py:class:`pyrocko.cake.PhaseDef`)

434 * ``'vel_surface:VELOCITY'`` - arrival according to surface distance /

435 velocity [km/s]

436 * ``'vel:VELOCITY'`` - arrival according to 3D-distance / velocity [km/s]

437

438 **Examples:**

439

440 * ``'100'`` : absolute time; 100 s

441 * ``'{stored:P}-100'`` : 100 s before arrival of P phase according to

442 stored travel time table named ``'P'``

443 * ``'{stored:P}-5.1S'`` : 10% before arrival of P phase according to

444 stored travel time table named ``'P'``

445 * ``'{stored:P}-10%'`` : 10% before arrival of P phase according to

446 stored travel time table named ``'P'``

447 * ``'{stored:A|stored:B}'`` : time instant of phase arrival A, or B if A is

448 undefined for a given geometry

449 * ``'first{stored:A|stored:B}'`` : as above, but the earlier arrival of A

450 and B is chosen, if both phases are defined for a given geometry

451 * ``'last{stored:A|stored:B}'`` : as above but the later arrival is chosen

452 * ``'first{stored:A|stored:B|stored:C}-100'`` : 100 s before first out of

453 arrivals A, B, and C

454 '''

455

456 def __init__(self, s=None, **kwargs):

457

458 if s is not None:

459 offset_is = None

460 s = re.sub(r'\s+', '', s)

461 try:

462 offset = float(s.rstrip('S'))

463

464 if s.endswith('S'):

465 offset_is = 'slowness'

466

467 phase_defs = []

468 select = ''

469

470 except ValueError:

471 matched = False

472 m = timing_regex.match(s)

473 if m:

474 if m.group(3):

475 phase_defs = m.group(3).split('|')

476 elif m.group(19):

477 phase_defs = [m.group(19)]

478 else:

479 phase_defs = []

480

481 select = m.group(2) or ''

482

483 offset = 0.0

484 soff = m.group(27)

485 if soff:

486 offset = float(soff.rstrip('S%'))

487 if soff.endswith('S'):

488 offset_is = 'slowness'

489 elif soff.endswith('%'):

490 offset_is = 'percent'

491

492 matched = True

493

494 else:

495 m = timing_regex_old.match(s)

496 if m:

497 if m.group(3):

498 phase_defs = m.group(3).split('|')

499 elif m.group(5):

500 phase_defs = [m.group(5)]

501 else:

502 phase_defs = []

503

504 select = m.group(2) or ''

505

506 offset = 0.0

507 if m.group(6):

508 offset = float(m.group(6))

509

510 phase_defs = [

511 'stored:' + phase_def for phase_def in phase_defs]

512

513 matched = True

514

515 if not matched:

516 raise InvalidTimingSpecification(s)

517

518 kwargs = dict(

519 phase_defs=phase_defs,

520 select=select,

521 offset=offset,

522 offset_is=offset_is)

523

524 SObject.__init__(self, **kwargs)

525

526 def __str__(self):

527 s = []

528 if self.phase_defs:

529 sphases = '|'.join(self.phase_defs)

530 # if len(self.phase_defs) > 1 or self.select:

531 sphases = '{%s}' % sphases

532

533 if self.select:

534 sphases = self.select + sphases

535

536 s.append(sphases)

537

538 if self.offset != 0.0 or not self.phase_defs:

539 s.append('%+g' % self.offset)

540 if self.offset_is == 'slowness':

541 s.append('S')

542 elif self.offset_is == 'percent':

543 s.append('%')

544

545 return ''.join(s)

546

547 def evaluate(self, get_phase, args, attributes=None):

548 try:

549 if self.offset_is == 'slowness' and self.offset != 0.0:

550 phase_offset = get_phase(

551 'vel_surface:%g' % (1.0/self.offset))

552 offset = phase_offset(args)

553 else:

554 offset = self.offset

555

556 if self.phase_defs:

557 extra_args = ()

558 if attributes:

559 extra_args = (attributes,)

560

561 phases = [

562 get_phase(phase_def, *extra_args)

563 for phase_def in self.phase_defs]

564

565 results = [phase(args) for phase in phases]

566 results = [

567 result if isinstance(result, tuple) else (result,)

568 for result in results]

569

570 results = [

571 result for result in results

572 if result[0] is not None]

573

574 if self.select == 'first':

575 results.sort(key=lambda result: result[0])

576 elif self.select == 'last':

577 results.sort(key=lambda result: -result[0])

578

579 if not results:

580 if attributes is not None:

581 return (None,) * (len(attributes) + 1)

582 else:

583 return None

584

585 else:

586 t, values = results[0][0], results[0][1:]

587 if self.offset_is == 'percent':

588 t = t*(1.+offset/100.)

589 else:

590 t = t + offset

591

592 if attributes is not None:

593 return (t, ) + values

594 else:

595 return t

596

597 else:

598 if attributes is not None:

599 raise TimingAttributeError(

600 'Cannot get phase attributes from offset-only '

601 'definition.')

602 return offset

603

604 except spit.OutOfBounds:

605 raise OutOfBounds(args)

606

607 phase_defs = List.T(String.T())

608 offset = Float.T(default=0.0)

609 offset_is = String.T(optional=True)

610 select = PhaseSelect.T(

611 default='',

612 help=("Can be either ``'%s'``, ``'%s'``, or ``'%s'``. " %

613 tuple(PhaseSelect.choices)))

614

615

616def mkdefs(s):

617 defs = []

618 for x in s.split(','):

619 try:

620 defs.append(float(x))

621 except ValueError:

622 if x.startswith('!'):

623 defs.extend(cake.PhaseDef.classic(x[1:]))

624 else:

625 defs.append(cake.PhaseDef(x))

626

627 return defs

628

629

630class TPDef(Object):

631 '''

632 Maps an arrival phase identifier to an arrival phase definition.

633 '''

634

635 id = StringID.T(

636 help='name used to identify the phase')

637 definition = String.T(

638 help='definition of the phase in either cake syntax as defined in '

639 ':py:class:`pyrocko.cake.PhaseDef`, or, if prepended with an '

640 '``!``, as a *classic phase name*, or, if it is a simple '

641 'number, as a constant horizontal velocity.')

642

643 @property

644 def phases(self):

645 return [x for x in mkdefs(self.definition)

646 if isinstance(x, cake.PhaseDef)]

647

648 @property

649 def horizontal_velocities(self):

650 return [x for x in mkdefs(self.definition) if isinstance(x, float)]

651

652

653class OutOfBounds(Exception):

654 def __init__(self, values=None, reason=None):

655 Exception.__init__(self)

656 self.values = values

657 self.reason = reason

658 self.context = None

659

660 def __str__(self):

661 scontext = ''

662 if self.context:

663 scontext = '\n%s' % str(self.context)

664

665 if self.reason:

666 scontext += '\n%s' % self.reason

667

668 if self.values:

669 return 'out of bounds: (%s)%s' % (

670 ', '.join('%g' % x for x in self.values), scontext)

671 else:

672 return 'out of bounds%s' % scontext

673

674

675class MultiLocation(Object):

676 '''

677 Unstructured grid of point coordinates.

678 '''

679

680 lats = Array.T(

681 optional=True, shape=(None,), dtype=float,

682 help='Latitudes of targets.')

683 lons = Array.T(

684 optional=True, shape=(None,), dtype=float,

685 help='Longitude of targets.')

686 north_shifts = Array.T(

687 optional=True, shape=(None,), dtype=float,

688 help='North shifts of targets.')

689 east_shifts = Array.T(

690 optional=True, shape=(None,), dtype=float,

691 help='East shifts of targets.')

692 elevation = Array.T(

693 optional=True, shape=(None,), dtype=float,

694 help='Elevations of targets.')

695

696 def __init__(self, *args, **kwargs):

697 self._coords5 = None

698 Object.__init__(self, *args, **kwargs)

699

700 @property

701 def coords5(self):

702 if self._coords5 is not None:

703 return self._coords5

704 props = [self.lats, self.lons, self.north_shifts, self.east_shifts,

705 self.elevation]

706 sizes = [p.size for p in props if p is not None]

707 if not sizes:

708 raise AttributeError('Empty StaticTarget')

709 if num.unique(sizes).size != 1:

710 raise AttributeError('Inconsistent coordinate sizes.')

711

712 self._coords5 = num.zeros((sizes[0], 5))

713 for idx, p in enumerate(props):

714 if p is not None:

715 self._coords5[:, idx] = p.flatten()

716

717 return self._coords5

718

719 @property

720 def ncoords(self):

721 if self.coords5.shape[0] is None:

722 return 0

723 return int(self.coords5.shape[0])

724

725 def get_latlon(self):

726 '''

727 Get all coordinates as lat/lon.

728

729 :returns: Coordinates in Latitude, Longitude

730 :rtype: :class:`numpy.ndarray`, (N, 2)

731 '''

732 latlons = num.empty((self.ncoords, 2))

733 for i in range(self.ncoords):

734 latlons[i, :] = orthodrome.ne_to_latlon(*self.coords5[i, :4])

735 return latlons

736

737 def get_corner_coords(self):

738 '''

739 Returns the corner coordinates of the multi-location object.

740

741 :returns: In lat/lon: lower left, upper left, upper right, lower right

742 :rtype: tuple

743 '''

744 latlon = self.get_latlon()

745 ll = (latlon[:, 0].min(), latlon[:, 1].min())

746 ur = (latlon[:, 0].max(), latlon[:, 1].max())

747 return (ll, (ll[0], ur[1]), ur, (ur[0], ll[1]))

748

749

750class Receiver(Location):

751 codes = Tuple.T(

752 3, String.T(),

753 optional=True,

754 help='network, station, and location codes')

755

756 def pyrocko_station(self):

757 from pyrocko import model

758

759 lat, lon = self.effective_latlon

760

761 return model.Station(

762 network=self.codes[0],

763 station=self.codes[1],

764 location=self.codes[2],

765 lat=lat,

766 lon=lon,

767 depth=self.depth)

768

769

770def g(x, d):

771 if x is None:

772 return d

773 else:

774 return x

775

776

777class UnavailableScheme(Exception):

778 '''

779 Raised when it is not possible to use the requested component scheme.

780 '''

781 pass

782

783

784class InvalidNComponents(Exception):

785 '''

786 Raised when ``ncomponents`` is incompatible with ``component_scheme``.

787 '''

788 pass

789

790

791class DiscretizedSource(Object):

792 '''

793 Base class for discretized sources.

794

795 To compute synthetic seismograms, the parameterized source models (see of

796 :py:class:`~pyrocko.gf.seismosizer.Source` derived classes) are first

797 discretized into a number of point sources. These spacio-temporal point

798 source distributions are represented by subclasses of the

799 :py:class:`DiscretizedSource`. For elastodynamic problems there is the

800 :py:class:`DiscretizedMTSource` for moment tensor point source

801 distributions and the :py:class:`DiscretizedExplosionSource` for pure

802 explosion/implosion type source distributions. The geometry calculations

803 are implemented in the base class. How Green's function components have to

804 be weighted and sumed is defined in the derived classes.

805

806 Like in the :py:class:`pyrocko.model.location.Location` class, the

807 positions of the point sources contained in the discretized source are

808 defined by a common reference point (:py:attr:`lat`, :py:attr:`lon`) and

809 cartesian offsets to this (:py:attr:`north_shifts`, :py:attr:`east_shifts`,

810 :py:attr:`depths`). Alternatively latitude and longitude of each contained

811 point source can be specified directly (:py:attr:`lats`, :py:attr:`lons`).

812 '''

813

814 times = Array.T(shape=(None,), dtype=float)

815 lats = Array.T(shape=(None,), dtype=float, optional=True)

816 lons = Array.T(shape=(None,), dtype=float, optional=True)

817 lat = Float.T(optional=True)

818 lon = Float.T(optional=True)

819 north_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)

820 east_shifts = Array.T(shape=(None,), dtype=num.float64, optional=True)

821 depths = Array.T(shape=(None,), dtype=num.float64)

822 dl = Float.T(optional=True)

823 dw = Float.T(optional=True)

824 nl = Float.T(optional=True)

825 nw = Float.T(optional=True)

826

827 @classmethod

828 def check_scheme(cls, scheme):

829 '''

830 Check if given GF component scheme is supported by the class.

831

832 Raises :py:class:`UnavailableScheme` if the given scheme is not

833 supported by this discretized source class.

834 '''

835

836 if scheme not in cls.provided_schemes:

837 raise UnavailableScheme(

838 'source type "%s" does not support GF component scheme "%s"' %

839 (cls.__name__, scheme))

840

841 def __init__(self, **kwargs):

842 Object.__init__(self, **kwargs)

843 self._latlons = None

844

845 def __setattr__(self, name, value):

846 if name in ('lat', 'lon', 'north_shifts', 'east_shifts',

847 'lats', 'lons'):

848 self.__dict__['_latlons'] = None

849

850 Object.__setattr__(self, name, value)

851

852 def get_source_terms(self, scheme):

853 raise NotImplementedError()

854

855 def make_weights(self, receiver, scheme):

856 raise NotImplementedError()

857

858 @property

859 def effective_latlons(self):

860 '''

861 Property holding the offest-corrected lats and lons of all points.

862 '''

863

864 if self._latlons is None:

865 if self.lats is not None and self.lons is not None:

866 if (self.north_shifts is not None and

867 self.east_shifts is not None):

868 self._latlons = orthodrome.ne_to_latlon(

869 self.lats, self.lons,

870 self.north_shifts, self.east_shifts)

871 else:

872 self._latlons = self.lats, self.lons

873 else:

874 lat = g(self.lat, 0.0)

875 lon = g(self.lon, 0.0)

876 self._latlons = orthodrome.ne_to_latlon(

877 lat, lon, self.north_shifts, self.east_shifts)

878

879 return self._latlons

880

881 @property

882 def effective_north_shifts(self):

883 if self.north_shifts is not None:

884 return self.north_shifts

885 else:

886 return num.zeros(self.times.size)

887

888 @property

889 def effective_east_shifts(self):

890 if self.east_shifts is not None:

891 return self.east_shifts

892 else:

893 return num.zeros(self.times.size)

894

895 def same_origin(self, receiver):

896 '''

897 Check if receiver has same reference point.

898 '''

899

900 return (g(self.lat, 0.0) == receiver.lat and

901 g(self.lon, 0.0) == receiver.lon and

902 self.lats is None and self.lons is None)

903

904 def azibazis_to(self, receiver):

905 '''

906 Compute azimuths and backaziumuths to/from receiver, for all contained

907 points.

908 '''

909

910 if self.same_origin(receiver):

911 azis = r2d * num.arctan2(receiver.east_shift - self.east_shifts,

912 receiver.north_shift - self.north_shifts)

913 bazis = azis + 180.

914 else:

915 slats, slons = self.effective_latlons

916 rlat, rlon = receiver.effective_latlon

917 azis = orthodrome.azimuth_numpy(slats, slons, rlat, rlon)

918 bazis = orthodrome.azimuth_numpy(rlat, rlon, slats, slons)

919

920 return azis, bazis

921

922 def distances_to(self, receiver):

923 '''

924 Compute distances to receiver for all contained points.

925 '''

926 if self.same_origin(receiver):

927 return num.sqrt((self.north_shifts - receiver.north_shift)**2 +

928 (self.east_shifts - receiver.east_shift)**2)

929

930 else:

931 slats, slons = self.effective_latlons

932 rlat, rlon = receiver.effective_latlon

933 return orthodrome.distance_accurate50m_numpy(slats, slons,

934 rlat, rlon)

935

936 def element_coords(self, i):

937 if self.lats is not None and self.lons is not None:

938 lat = float(self.lats[i])

939 lon = float(self.lons[i])

940 else:

941 lat = self.lat

942 lon = self.lon

943

944 if self.north_shifts is not None and self.east_shifts is not None:

945 north_shift = float(self.north_shifts[i])

946 east_shift = float(self.east_shifts[i])

947

948 else:

949 north_shift = east_shift = 0.0

950

951 return lat, lon, north_shift, east_shift

952

953 def coords5(self):

954 xs = num.zeros((self.nelements, 5))

955

956 if self.lats is not None and self.lons is not None:

957 xs[:, 0] = self.lats

958 xs[:, 1] = self.lons

959 else:

960 xs[:, 0] = g(self.lat, 0.0)

961 xs[:, 1] = g(self.lon, 0.0)

962

963 if self.north_shifts is not None and self.east_shifts is not None:

964 xs[:, 2] = self.north_shifts

965 xs[:, 3] = self.east_shifts

966

967 xs[:, 4] = self.depths

968

969 return xs

970

971 @property

972 def nelements(self):

973 return self.times.size

974

975 @classmethod

976 def combine(cls, sources, **kwargs):

977 '''

978 Combine several discretized source models.

979

980 Concatenenates all point sources in the given discretized ``sources``.

981 Care must be taken when using this function that the external amplitude

982 factors and reference times of the parameterized (undiscretized)

983 sources match or are accounted for.

984 '''

985

986 first = sources[0]

987

988 if not all(type(s) == type(first) for s in sources):

989 raise Exception('DiscretizedSource.combine must be called with '

990 'sources of same type.')

991

992 latlons = []

993 for s in sources:

994 latlons.append(s.effective_latlons)

995

996 lats, lons = num.hstack(latlons)

997

998 if all((s.lats is None and s.lons is None) for s in sources):

999 rlats = num.array([s.lat for s in sources], dtype=float)

1000 rlons = num.array([s.lon for s in sources], dtype=float)

1001 same_ref = num.all(

1002 rlats == rlats[0]) and num.all(rlons == rlons[0])

1003 else:

1004 same_ref = False

1005

1006 cat = num.concatenate

1007 times = cat([s.times for s in sources])

1008 print(times)

1009 depths = cat([s.depths for s in sources])

1010

1011 if same_ref:

1012 lat = first.lat

1013 lon = first.lon

1014 north_shifts = cat([s.effective_north_shifts for s in sources])

1015 east_shifts = cat([s.effective_east_shifts for s in sources])

1016 lats = None

1017 lons = None

1018 else:

1019 lat = None

1020 lon = None

1021 north_shifts = None

1022 east_shifts = None

1023

1024 return cls(

1025 times=times, lat=lat, lon=lon, lats=lats, lons=lons,

1026 north_shifts=north_shifts, east_shifts=east_shifts,

1027 depths=depths, **kwargs)

1028

1029 def centroid_position(self):

1030 moments = self.moments()

1031 norm = num.sum(moments)

1032 if norm != 0.0:

1033 w = moments / num.sum(moments)

1034 else:

1035 w = num.ones(moments.size)

1036

1037 if self.lats is not None and self.lons is not None:

1038 lats, lons = self.effective_latlons

1039 rlat, rlon = num.mean(lats), num.mean(lons)

1040 n, e = orthodrome.latlon_to_ne_numpy(rlat, rlon, lats, lons)

1041 else:

1042 rlat, rlon = g(self.lat, 0.0), g(self.lon, 0.0)

1043 n, e = self.north_shifts, self.east_shifts

1044

1045 cn = num.sum(n*w)

1046 ce = num.sum(e*w)

1047 clat, clon = orthodrome.ne_to_latlon(rlat, rlon, cn, ce)

1048

1049 if self.lats is not None and self.lons is not None:

1050 lat = clat

1051 lon = clon

1052 north_shift = 0.

1053 east_shift = 0.

1054 else:

1055 lat = g(self.lat, 0.0)

1056 lon = g(self.lon, 0.0)

1057 north_shift = cn

1058 east_shift = ce

1059

1060 depth = num.sum(self.depths*w)

1061 time = num.sum(self.times*w)

1062 return tuple(float(x) for x in

1063 (time, lat, lon, north_shift, east_shift, depth))

1064

1065

1066class DiscretizedExplosionSource(DiscretizedSource):

1067 m0s = Array.T(shape=(None,), dtype=float)

1068

1069 provided_schemes = (

1070 'elastic2',

1071 'elastic8',

1072 'elastic10',

1073 )

1074

1075 def get_source_terms(self, scheme):

1076 self.check_scheme(scheme)

1077

1078 if scheme == 'elastic2':

1079 return self.m0s[:, num.newaxis].copy()

1080

1081 elif scheme in ('elastic8', 'elastic10'):

1082 m6s = num.zeros((self.m0s.size, 6))

1083 m6s[:, 0:3] = self.m0s[:, num.newaxis]

1084 return m6s

1085 else:

1086 assert False

1087

1088 def make_weights(self, receiver, scheme):

1089 self.check_scheme(scheme)

1090

1091 azis, bazis = self.azibazis_to(receiver)

1092

1093 sb = num.sin(bazis*d2r-num.pi)

1094 cb = num.cos(bazis*d2r-num.pi)

1095

1096 m0s = self.m0s

1097 n = azis.size

1098

1099 cat = num.concatenate

1100 rep = num.repeat

1101

1102 if scheme == 'elastic2':

1103 w_n = cb*m0s

1104 g_n = filledi(0, n)

1105 w_e = sb*m0s

1106 g_e = filledi(0, n)

1107 w_d = m0s

1108 g_d = filledi(1, n)

1109

1110 elif scheme == 'elastic8':

1111 w_n = cat((cb*m0s, cb*m0s))

1112 g_n = rep((0, 2), n)

1113 w_e = cat((sb*m0s, sb*m0s))

1114 g_e = rep((0, 2), n)

1115 w_d = cat((m0s, m0s))

1116 g_d = rep((5, 7), n)

1117

1118 elif scheme == 'elastic10':

1119 w_n = cat((cb*m0s, cb*m0s, cb*m0s))

1120 g_n = rep((0, 2, 8), n)

1121 w_e = cat((sb*m0s, sb*m0s, sb*m0s))

1122 g_e = rep((0, 2, 8), n)

1123 w_d = cat((m0s, m0s, m0s))

1124 g_d = rep((5, 7, 9), n)

1125

1126 else:

1127 assert False

1128

1129 return (

1130 ('displacement.n', w_n, g_n),

1131 ('displacement.e', w_e, g_e),

1132 ('displacement.d', w_d, g_d),

1133 )

1134

1135 def split(self):

1136 from pyrocko.gf.seismosizer import ExplosionSource

1137 sources = []

1138 for i in range(self.nelements):

1139 lat, lon, north_shift, east_shift = self.element_coords(i)

1140 sources.append(ExplosionSource(

1141 time=float(self.times[i]),

1142 lat=lat,

1143 lon=lon,

1144 north_shift=north_shift,

1145 east_shift=east_shift,

1146 depth=float(self.depths[i]),

1147 moment=float(self.m0s[i])))

1148

1149 return sources

1150

1151 def moments(self):

1152 return self.m0s

1153

1154 def centroid(self):

1155 from pyrocko.gf.seismosizer import ExplosionSource

1156 time, lat, lon, north_shift, east_shift, depth = \

1157 self.centroid_position()

1158

1159 return ExplosionSource(

1160 time=time,

1161 lat=lat,

1162 lon=lon,

1163 north_shift=north_shift,

1164 east_shift=east_shift,

1165 depth=depth,

1166 moment=float(num.sum(self.m0s)))

1167

1168 @classmethod

1169 def combine(cls, sources, **kwargs):

1170 '''

1171 Combine several discretized source models.

1172

1173 Concatenenates all point sources in the given discretized ``sources``.

1174 Care must be taken when using this function that the external amplitude

1175 factors and reference times of the parameterized (undiscretized)

1176 sources match or are accounted for.

1177 '''

1178

1179 if 'm0s' not in kwargs:

1180 kwargs['m0s'] = num.concatenate([s.m0s for s in sources])

1181

1182 return super(DiscretizedExplosionSource, cls).combine(sources,

1183 **kwargs)

1184

1185

1186class DiscretizedSFSource(DiscretizedSource):

1187 forces = Array.T(shape=(None, 3), dtype=float)

1188

1189 provided_schemes = (

1190 'elastic5',

1191 )

1192

1193 def get_source_terms(self, scheme):

1194 self.check_scheme(scheme)

1195

1196 return self.forces

1197

1198 def make_weights(self, receiver, scheme):

1199 self.check_scheme(scheme)

1200

1201 azis, bazis = self.azibazis_to(receiver)

1202

1203 sa = num.sin(azis*d2r)

1204 ca = num.cos(azis*d2r)

1205 sb = num.sin(bazis*d2r-num.pi)

1206 cb = num.cos(bazis*d2r-num.pi)

1207

1208 forces = self.forces

1209 fn = forces[:, 0]

1210 fe = forces[:, 1]

1211 fd = forces[:, 2]

1212

1213 f0 = fd

1214 f1 = ca * fn + sa * fe

1215 f2 = ca * fe - sa * fn

1216

1217 n = azis.size

1218

1219 if scheme == 'elastic5':

1220 ioff = 0

1221

1222 cat = num.concatenate

1223 rep = num.repeat

1224

1225 w_n = cat((cb*f0, cb*f1, -sb*f2))

1226 g_n = ioff + rep((0, 1, 2), n)

1227 w_e = cat((sb*f0, sb*f1, cb*f2))

1228 g_e = ioff + rep((0, 1, 2), n)

1229 w_d = cat((f0, f1))

1230 g_d = ioff + rep((3, 4), n)

1231

1232 return (

1233 ('displacement.n', w_n, g_n),

1234 ('displacement.e', w_e, g_e),

1235 ('displacement.d', w_d, g_d),

1236 )

1237

1238 @classmethod

1239 def combine(cls, sources, **kwargs):

1240 '''

1241 Combine several discretized source models.

1242

1243 Concatenenates all point sources in the given discretized ``sources``.

1244 Care must be taken when using this function that the external amplitude

1245 factors and reference times of the parameterized (undiscretized)

1246 sources match or are accounted for.

1247 '''

1248

1249 if 'forces' not in kwargs:

1250 kwargs['forces'] = num.vstack([s.forces for s in sources])

1251

1252 return super(DiscretizedSFSource, cls).combine(sources, **kwargs)

1253

1254 def moments(self):

1255 return num.sum(self.forces**2, axis=1)

1256

1257 def centroid(self):

1258 from pyrocko.gf.seismosizer import SFSource

1259 time, lat, lon, north_shift, east_shift, depth = \

1260 self.centroid_position()

1261

1262 fn, fe, fd = map(float, num.sum(self.forces, axis=0))

1263 return SFSource(

1264 time=time,

1265 lat=lat,

1266 lon=lon,

1267 north_shift=north_shift,

1268 east_shift=east_shift,

1269 depth=depth,

1270 fn=fn,

1271 fe=fe,

1272 fd=fd)

1273

1274

1275class DiscretizedMTSource(DiscretizedSource):

1276 m6s = Array.T(

1277 shape=(None, 6), dtype=float,

1278 help='rows with (m_nn, m_ee, m_dd, m_ne, m_nd, m_ed)')

1279

1280 provided_schemes = (

1281 'elastic8',

1282 'elastic10',

1283 'elastic18',

1284 )

1285

1286 def get_source_terms(self, scheme):

1287 self.check_scheme(scheme)

1288 return self.m6s

1289

1290 def make_weights(self, receiver, scheme):

1291 self.check_scheme(scheme)

1292

1293 m6s = self.m6s

1294 n = m6s.shape[0]

1295 rep = num.repeat

1296

1297 if scheme == 'elastic18':

1298 w_n = m6s.flatten()

1299 g_n = num.tile(num.arange(0, 6), n)

1300 w_e = m6s.flatten()

1301 g_e = num.tile(num.arange(6, 12), n)

1302 w_d = m6s.flatten()

1303 g_d = num.tile(num.arange(12, 18), n)

1304

1305 else:

1306 azis, bazis = self.azibazis_to(receiver)

1307

1308 sa = num.sin(azis*d2r)

1309 ca = num.cos(azis*d2r)

1310 s2a = num.sin(2.*azis*d2r)

1311 c2a = num.cos(2.*azis*d2r)

1312 sb = num.sin(bazis*d2r-num.pi)

1313 cb = num.cos(bazis*d2r-num.pi)

1314

1315 f0 = m6s[:, 0]*ca**2 + m6s[:, 1]*sa**2 + m6s[:, 3]*s2a

1316 f1 = m6s[:, 4]*ca + m6s[:, 5]*sa

1317 f2 = m6s[:, 2]

1318 f3 = 0.5*(m6s[:, 1]-m6s[:, 0])*s2a + m6s[:, 3]*c2a

1319 f4 = m6s[:, 5]*ca - m6s[:, 4]*sa

1320 f5 = m6s[:, 0]*sa**2 + m6s[:, 1]*ca**2 - m6s[:, 3]*s2a

1321

1322 cat = num.concatenate

1323

1324 if scheme == 'elastic8':

1325 w_n = cat((cb*f0, cb*f1, cb*f2, -sb*f3, -sb*f4))

1326 g_n = rep((0, 1, 2, 3, 4), n)

1327 w_e = cat((sb*f0, sb*f1, sb*f2, cb*f3, cb*f4))

1328 g_e = rep((0, 1, 2, 3, 4), n)

1329 w_d = cat((f0, f1, f2))

1330 g_d = rep((5, 6, 7), n)

1331

1332 elif scheme == 'elastic10':

1333 w_n = cat((cb*f0, cb*f1, cb*f2, cb*f5, -sb*f3, -sb*f4))

1334 g_n = rep((0, 1, 2, 8, 3, 4), n)

1335 w_e = cat((sb*f0, sb*f1, sb*f2, sb*f5, cb*f3, cb*f4))

1336 g_e = rep((0, 1, 2, 8, 3, 4), n)

1337 w_d = cat((f0, f1, f2, f5))

1338 g_d = rep((5, 6, 7, 9), n)

1339

1340 else:

1341 assert False

1342

1343 return (

1344 ('displacement.n', w_n, g_n),

1345 ('displacement.e', w_e, g_e),

1346 ('displacement.d', w_d, g_d),

1347 )

1348

1349 def split(self):

1350 from pyrocko.gf.seismosizer import MTSource

1351 sources = []

1352 for i in range(self.nelements):

1353 lat, lon, north_shift, east_shift = self.element_coords(i)

1354 sources.append(MTSource(

1355 time=float(self.times[i]),

1356 lat=lat,

1357 lon=lon,

1358 north_shift=north_shift,

1359 east_shift=east_shift,

1360 depth=float(self.depths[i]),

1361 m6=self.m6s[i]))

1362

1363 return sources

1364

1365 def moments(self):

1366 moments = num.array(

1367 [num.linalg.eigvalsh(moment_tensor.symmat6(*m6))

1368 for m6 in self.m6s])

1369 return num.linalg.norm(moments, axis=1) / num.sqrt(2.)

1370

1371 def get_moment_rate(self, deltat=None):

1372 moments = self.moments()

1373 times = self.times

1374 times -= times.min()

1375

1376 t_max = times.max()

1377 mom_times = num.arange(0, t_max + 2 * deltat, deltat) - deltat

1378 mom_times[mom_times > t_max] = t_max

1379

1380 # Right open histrogram bins

1381 mom, _ = num.histogram(

1382 times,

1383 bins=mom_times,

1384 weights=moments)

1385

1386 deltat = num.diff(mom_times)

1387 mom_rate = mom / deltat

1388

1389 return mom_rate, mom_times[1:]

1390

1391 def centroid(self):

1392 from pyrocko.gf.seismosizer import MTSource

1393 time, lat, lon, north_shift, east_shift, depth = \

1394 self.centroid_position()

1395

1396 return MTSource(

1397 time=time,

1398 lat=lat,

1399 lon=lon,

1400 north_shift=north_shift,

1401 east_shift=east_shift,

1402 depth=depth,

1403 m6=num.sum(self.m6s, axis=0))

1404

1405 @classmethod

1406 def combine(cls, sources, **kwargs):

1407 '''

1408 Combine several discretized source models.

1409

1410 Concatenenates all point sources in the given discretized ``sources``.

1411 Care must be taken when using this function that the external amplitude

1412 factors and reference times of the parameterized (undiscretized)

1413 sources match or are accounted for.

1414 '''

1415

1416 if 'm6s' not in kwargs:

1417 kwargs['m6s'] = num.vstack([s.m6s for s in sources])

1418

1419 return super(DiscretizedMTSource, cls).combine(sources, **kwargs)

1420

1421

1422class DiscretizedPorePressureSource(DiscretizedSource):

1423 pp = Array.T(shape=(None,), dtype=float)

1424

1425 provided_schemes = (

1426 'poroelastic10',

1427 )

1428

1429 def get_source_terms(self, scheme):

1430 self.check_scheme(scheme)

1431 return self.pp[:, num.newaxis].copy()

1432

1433 def make_weights(self, receiver, scheme):

1434 self.check_scheme(scheme)

1435

1436 azis, bazis = self.azibazis_to(receiver)

1437

1438 sb = num.sin(bazis*d2r-num.pi)

1439 cb = num.cos(bazis*d2r-num.pi)

1440

1441 pp = self.pp

1442 n = bazis.size

1443

1444 w_un = cb*pp

1445 g_un = filledi(1, n)

1446 w_ue = sb*pp

1447 g_ue = filledi(1, n)

1448 w_ud = pp

1449 g_ud = filledi(0, n)

1450

1451 w_tn = cb*pp

1452 g_tn = filledi(6, n)

1453 w_te = sb*pp

1454 g_te = filledi(6, n)

1455

1456 w_pp = pp

1457 g_pp = filledi(7, n)

1458

1459 w_dvn = cb*pp

1460 g_dvn = filledi(9, n)

1461 w_dve = sb*pp

1462 g_dve = filledi(9, n)

1463 w_dvd = pp

1464 g_dvd = filledi(8, n)

1465

1466 return (

1467 ('displacement.n', w_un, g_un),

1468 ('displacement.e', w_ue, g_ue),

1469 ('displacement.d', w_ud, g_ud),

1470 ('vertical_tilt.n', w_tn, g_tn),

1471 ('vertical_tilt.e', w_te, g_te),

1472 ('pore_pressure', w_pp, g_pp),

1473 ('darcy_velocity.n', w_dvn, g_dvn),

1474 ('darcy_velocity.e', w_dve, g_dve),

1475 ('darcy_velocity.d', w_dvd, g_dvd),

1476 )

1477

1478 def moments(self):

1479 return self.pp

1480

1481 def centroid(self):

1482 from pyrocko.gf.seismosizer import PorePressurePointSource

1483 time, lat, lon, north_shift, east_shift, depth = \

1484 self.centroid_position()

1485

1486 return PorePressurePointSource(

1487 time=time,

1488 lat=lat,

1489 lon=lon,

1490 north_shift=north_shift,

1491 east_shift=east_shift,

1492 depth=depth,

1493 pp=float(num.sum(self.pp)))

1494

1495 @classmethod

1496 def combine(cls, sources, **kwargs):

1497 '''

1498 Combine several discretized source models.

1499

1500 Concatenenates all point sources in the given discretized ``sources``.

1501 Care must be taken when using this function that the external amplitude

1502 factors and reference times of the parameterized (undiscretized)

1503 sources match or are accounted for.

1504 '''

1505

1506 if 'pp' not in kwargs:

1507 kwargs['pp'] = num.concatenate([s.pp for s in sources])

1508

1509 return super(DiscretizedPorePressureSource, cls).combine(sources,

1510 **kwargs)

1511

1512

1513class Region(Object):

1514 name = String.T(optional=True)

1515

1516

1517class RectangularRegion(Region):

1518 lat_min = Float.T()

1519 lat_max = Float.T()

1520 lon_min = Float.T()

1521 lon_max = Float.T()

1522

1523

1524class CircularRegion(Region):

1525 lat = Float.T()

1526 lon = Float.T()

1527 radius = Float.T()

1528

1529

1530class Config(Object):

1531 '''

1532 Green's function store meta information.

1533

1534 Currently implemented :py:class:`~pyrocko.gf.store.Store`

1535 configuration types are:

1536

1537 * :py:class:`~pyrocko.gf.meta.ConfigTypeA` - cylindrical or

1538 spherical symmetry, 1D earth model, single receiver depth

1539

1540 * Problem is invariant to horizontal translations and rotations around

1541 vertical axis.

1542 * All receivers must be at the same depth (e.g. at the surface)

1543 * High level index variables: ``(source_depth, receiver_distance,

1544 component)``

1545

1546 * :py:class:`~pyrocko.gf.meta.ConfigTypeB` - cylindrical or

1547 spherical symmetry, 1D earth model, variable receiver depth

1548

1549 * Symmetries like in Type A but has additional index for receiver depth

1550 * High level index variables: ``(source_depth, receiver_distance,

1551 receiver_depth, component)``

1552

1553 * :py:class:`~pyrocko.gf.meta.ConfigTypeC` - no symmetrical

1554 constraints but fixed receiver positions

1555

1556 * Cartesian source volume around a reference point

1557 * High level index variables: ``(ireceiver, source_depth,

1558 source_east_shift, source_north_shift, component)``

1559 '''

1560

1561 id = StringID.T(

1562 help='Name of the store. May consist of upper and lower-case letters, '

1563 'digits, dots and underscores. The name must start with a '

1564 'letter.')

1565

1566 derived_from_id = StringID.T(

1567 optional=True,

1568 help='Name of the original store, if this store has been derived from '

1569 'another one (e.g. extracted subset).')

1570

1571 version = String.T(

1572 default='1.0',

1573 optional=True,

1574 help='User-defined version string. Use <major>.<minor> format.')

1575

1576 modelling_code_id = StringID.T(

1577 optional=True,

1578 help='Identifier of the backend used to compute the store.')

1579

1580 author = Unicode.T(

1581 optional=True,

1582 help='Comma-separated list of author names.')

1583

1584 author_email = String.T(

1585 optional=True,

1586 help="Author's contact email address.")

1587

1588 created_time = Timestamp.T(

1589 optional=True,

1590 help='Time of creation of the store.')

1591

1592 regions = List.T(

1593 Region.T(),

1594 help='Geographical regions for which the store is representative.')

1595

1596 scope_type = ScopeType.T(

1597 optional=True,

1598 help='Distance range scope of the store '

1599 '(%s).' % fmt_choices(ScopeType))

1600

1601 waveform_type = WaveType.T(

1602 optional=True,

1603 help='Wave type stored (%s).' % fmt_choices(WaveType))

1604

1605 nearfield_terms = NearfieldTermsType.T(

1606 optional=True,

1607 help='Information about the inclusion of near-field terms in the '

1608 'modelling (%s).' % fmt_choices(NearfieldTermsType))

1609

1610 description = String.T(

1611 optional=True,

1612 help='Free form textual description of the GF store.')

1613

1614 references = List.T(

1615 Reference.T(),

1616 help='Reference list to cite the modelling code, earth model or '

1617 'related work.')

1618

1619 earthmodel_1d = Earthmodel1D.T(

1620 optional=True,

1621 help='Layered earth model in ND (named discontinuity) format.')

1622

1623 earthmodel_receiver_1d = Earthmodel1D.T(

1624 optional=True,

1625 help='Receiver-side layered earth model in ND format.')

1626

1627 can_interpolate_source = Bool.T(

1628 optional=True,

1629 help='Hint to indicate if the spatial sampling of the store is dense '

1630 'enough for multi-linear interpolation at the source.')

1631

1632 can_interpolate_receiver = Bool.T(

1633 optional=True,

1634 help='Hint to indicate if the spatial sampling of the store is dense '

1635 'enough for multi-linear interpolation at the receiver.')

1636

1637 frequency_min = Float.T(

1638 optional=True,

1639 help='Hint to indicate the lower bound of valid frequencies [Hz].')

1640

1641 frequency_max = Float.T(

1642 optional=True,

1643 help='Hint to indicate the upper bound of valid frequencies [Hz].')

1644

1645 sample_rate = Float.T(

1646 optional=True,

1647 help='Sample rate of the GF store [Hz].')

1648

1649 factor = Float.T(

1650 default=1.0,

1651 help='Gain value, factored out of the stored GF samples. '

1652 '(may not work properly, keep at 1.0).',

1653 optional=True)

1654

1655 component_scheme = ComponentScheme.T(

1656 default='elastic10',

1657 help='GF component scheme (%s).' % fmt_choices(ComponentScheme))

1658

1659 stored_quantity = QuantityType.T(

1660 optional=True,

1661 help='Physical quantity of stored values (%s). If not given, a '

1662 'default is used based on the GF component scheme. The default '

1663 'for the ``"elastic*"`` family of component schemes is '

1664 '``"displacement"``.' % fmt_choices(QuantityType))

1665

1666 tabulated_phases = List.T(

1667 TPDef.T(),

1668 help='Mapping of phase names to phase definitions, for which travel '

1669 'time tables are available in the GF store.')

1670

1671 ncomponents = Int.T(

1672 optional=True,

1673 help='Number of GF components. Use :gattr:`component_scheme` instead.')

1674

1675 uuid = String.T(

1676 optional=True,

1677 help='Heuristic hash value which can be used to uniquely identify the '

1678 'GF store for practical purposes.')

1679

1680 reference = String.T(

1681 optional=True,

1682 help="Store reference name composed of the store's :gattr:`id` and "

1683 'the first six letters of its :gattr:`uuid`.')

1684

1685 def __init__(self, **kwargs):

1686 self._do_auto_updates = False

1687 Object.__init__(self, **kwargs)

1688 self._index_function = None

1689 self._indices_function = None

1690 self._vicinity_function = None

1691 self.validate(regularize=True, depth=1)

1692 self._do_auto_updates = True

1693 self.update()

1694

1695 def check_ncomponents(self):

1696 ncomponents = component_scheme_to_description[

1697 self.component_scheme].ncomponents

1698

1699 if self.ncomponents is None:

1700 self.ncomponents = ncomponents

1701 elif ncomponents != self.ncomponents:

1702 raise InvalidNComponents(

1703 'ncomponents=%i incompatible with component_scheme="%s"' % (

1704 self.ncomponents, self.component_scheme))

1705

1706 def __setattr__(self, name, value):

1707 Object.__setattr__(self, name, value)

1708 try:

1709 self.T.get_property(name)

1710 if self._do_auto_updates:

1711 self.update()

1712

1713 except ValueError:

1714 pass

1715

1716 def update(self):

1717 self.check_ncomponents()

1718 self._update()

1719 self._make_index_functions()

1720

1721 def irecord(self, *args):

1722 return self._index_function(*args)

1723

1724 def irecords(self, *args):

1725 return self._indices_function(*args)

1726

1727 def vicinity(self, *args):

1728 return self._vicinity_function(*args)

1729

1730 def vicinities(self, *args):

1731 return self._vicinities_function(*args)

1732

1733 def grid_interpolation_coefficients(self, *args):

1734 return self._grid_interpolation_coefficients(*args)

1735

1736 def nodes(self, level=None, minlevel=None):

1737 return nodes(self.coords[minlevel:level])

1738

1739 def iter_nodes(self, level=None, minlevel=None):

1740 return nditer_outer(self.coords[minlevel:level])

1741

1742 def iter_extraction(self, gdef, level=None):

1743 i = 0

1744 arrs = []

1745 ntotal = 1

1746 for mi, ma, inc in zip(self.mins, self.effective_maxs, self.deltas):

1747 if gdef and len(gdef) > i:

1748 sssn = gdef[i]

1749 else:

1750 sssn = (None,)*4

1751

1752 arr = num.linspace(*start_stop_num(*(sssn + (mi, ma, inc))))

1753 ntotal *= len(arr)

1754

1755 arrs.append(arr)

1756 i += 1

1757

1758 arrs.append(self.coords[-1])

1759 return nditer_outer(arrs[:level])

1760

1761 def make_sum_params(self, source, receiver, implementation='c',

1762 nthreads=0):

1763 assert implementation in ['c', 'python']

1764

1765 out = []

1766 delays = source.times

1767 for comp, weights, icomponents in source.make_weights(

1768 receiver,

1769 self.component_scheme):

1770

1771 weights *= self.factor

1772

1773 args = self.make_indexing_args(source, receiver, icomponents)

1774 delays_expanded = num.tile(delays, icomponents.size//delays.size)

1775 out.append((comp, args, delays_expanded, weights))

1776

1777 return out

1778

1779 def short_info(self):

1780 raise NotImplementedError('should be implemented in subclass')

1781

1782 def get_shear_moduli(self, lat, lon, points,

1783 interpolation=None):

1784 '''

1785 Get shear moduli at given points from contained velocity model.

1786

1787 :param lat: surface origin for coordinate system of ``points``

1788 :param points: NumPy array of shape ``(N, 3)``, where each row is

1789 a point ``(north, east, depth)``, relative to origin at

1790 ``(lat, lon)``

1791 :param interpolation: interpolation method. Choose from

1792 ``('nearest_neighbor', 'multilinear')``

1793 :returns: NumPy array of length N with extracted shear moduli at each

1794 point

1795

1796 The default implementation retrieves and interpolates the shear moduli

1797 from the contained 1D velocity profile.

1798 '''

1799 return self.get_material_property(lat, lon, points,

1800 parameter='shear_moduli',

1801 interpolation=interpolation)

1802

1803 def get_lambda_moduli(self, lat, lon, points,

1804 interpolation=None):

1805 '''

1806 Get lambda moduli at given points from contained velocity model.

1807

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

1816

1817 The default implementation retrieves and interpolates the lambda moduli

1818 from the contained 1D velocity profile.

1819 '''

1820 return self.get_material_property(lat, lon, points,

1821 parameter='lambda_moduli',

1822 interpolation=interpolation)

1823

1824 def get_bulk_moduli(self, lat, lon, points,

1825 interpolation=None):

1826 '''

1827 Get bulk moduli at given points from contained velocity model.

1828

1829 :param lat: surface origin for coordinate system of ``points``

1830 :param points: NumPy array of shape ``(N, 3)``, where each row is

1831 a point ``(north, east, depth)``, relative to origin at

1832 ``(lat, lon)``

1833 :param interpolation: interpolation method. Choose from

1834 ``('nearest_neighbor', 'multilinear')``

1835 :returns: NumPy array of length N with extracted shear moduli at each

1836 point

1837

1838 The default implementation retrieves and interpolates the lambda moduli

1839 from the contained 1D velocity profile.

1840 '''

1841 lambda_moduli = self.get_material_property(

1842 lat, lon, points, parameter='lambda_moduli',

1843 interpolation=interpolation)

1844 shear_moduli = self.get_material_property(

1845 lat, lon, points, parameter='shear_moduli',

1846 interpolation=interpolation)

1847 return lambda_moduli + (2 / 3) * shear_moduli

1848

1849 def get_vs(self, lat, lon, points, interpolation=None):

1850 '''

1851 Get Vs at given points from contained velocity model.

1852

1853 :param lat: surface origin for coordinate system of ``points``

1854 :param points: NumPy array of shape ``(N, 3)``, where each row is

1855 a point ``(north, east, depth)``, relative to origin at

1856 ``(lat, lon)``

1857 :param interpolation: interpolation method. Choose from

1858 ``('nearest_neighbor', 'multilinear')``

1859 :returns: NumPy array of length N with extracted shear moduli at each

1860 point

1861

1862 The default implementation retrieves and interpolates Vs

1863 from the contained 1D velocity profile.

1864 '''

1865 return self.get_material_property(lat, lon, points,

1866 parameter='vs',

1867 interpolation=interpolation)

1868

1869 def get_vp(self, lat, lon, points, interpolation=None):

1870 '''

1871 Get Vp at given points from contained velocity model.

1872

1873 :param lat: surface origin for coordinate system of ``points``

1874 :param points: NumPy array of shape ``(N, 3)``, where each row is

1875 a point ``(north, east, depth)``, relative to origin at

1876 ``(lat, lon)``

1877 :param interpolation: interpolation method. Choose from

1878 ``('nearest_neighbor', 'multilinear')``

1879 :returns: NumPy array of length N with extracted shear moduli at each

1880 point

1881

1882 The default implementation retrieves and interpolates Vp

1883 from the contained 1D velocity profile.

1884 '''

1885 return self.get_material_property(lat, lon, points,

1886 parameter='vp',

1887 interpolation=interpolation)

1888

1889 def get_rho(self, lat, lon, points, interpolation=None):

1890 '''

1891 Get rho at given points from contained velocity model.

1892

1893 :param lat: surface origin for coordinate system of ``points``

1894 :param points: NumPy array of shape ``(N, 3)``, where each row is

1895 a point ``(north, east, depth)``, relative to origin at

1896 ``(lat, lon)``

1897 :param interpolation: interpolation method. Choose from

1898 ``('nearest_neighbor', 'multilinear')``

1899 :returns: NumPy array of length N with extracted shear moduli at each

1900 point

1901

1902 The default implementation retrieves and interpolates rho

1903 from the contained 1D velocity profile.

1904 '''

1905 return self.get_material_property(lat, lon, points,

1906 parameter='rho',

1907 interpolation=interpolation)

1908

1909 def get_material_property(self, lat, lon, points, parameter='vs',

1910 interpolation=None):

1911

1912 if interpolation is None:

1913 raise TypeError('Interpolation method not defined! available: '

1914 'multilinear', 'nearest_neighbor')

1915

1916 earthmod = self.earthmodel_1d

1917 store_depth_profile = self.get_source_depths()

1918 z_profile = earthmod.profile('z')

1919

1920 if parameter == 'vs':

1921 vs_profile = earthmod.profile('vs')

1922 profile = num.interp(

1923 store_depth_profile, z_profile, vs_profile)

1924

1925 elif parameter == 'vp':

1926 vp_profile = earthmod.profile('vp')

1927 profile = num.interp(

1928 store_depth_profile, z_profile, vp_profile)

1929

1930 elif parameter == 'rho':

1931 rho_profile = earthmod.profile('rho')

1932

1933 profile = num.interp(

1934 store_depth_profile, z_profile, rho_profile)

1935

1936 elif parameter == 'shear_moduli':

1937 vs_profile = earthmod.profile('vs')

1938 rho_profile = earthmod.profile('rho')

1939

1940 store_vs_profile = num.interp(

1941 store_depth_profile, z_profile, vs_profile)

1942 store_rho_profile = num.interp(

1943 store_depth_profile, z_profile, rho_profile)

1944

1945 profile = num.power(store_vs_profile, 2) * store_rho_profile

1946

1947 elif parameter == 'lambda_moduli':

1948 vs_profile = earthmod.profile('vs')

1949 vp_profile = earthmod.profile('vp')

1950 rho_profile = earthmod.profile('rho')

1951

1952 store_vs_profile = num.interp(

1953 store_depth_profile, z_profile, vs_profile)

1954 store_vp_profile = num.interp(

1955 store_depth_profile, z_profile, vp_profile)

1956 store_rho_profile = num.interp(

1957 store_depth_profile, z_profile, rho_profile)

1958

1959 profile = store_rho_profile * (

1960 num.power(store_vp_profile, 2) -

1961 num.power(store_vs_profile, 2) * 2)

1962 else:

1963 raise TypeError(

1964 'parameter %s not available' % parameter)

1965

1966 if interpolation == 'multilinear':

1967 kind = 'linear'

1968 elif interpolation == 'nearest_neighbor':

1969 kind = 'nearest'

1970 else:

1971 raise TypeError(

1972 'Interpolation method %s not available' % interpolation)

1973

1974 interpolator = interp1d(store_depth_profile, profile, kind=kind)

1975

1976 try:

1977 return interpolator(points[:, 2])

1978 except ValueError:

1979 raise OutOfBounds()

1980

1981 def is_static(self):

1982 for code in ('psgrn_pscmp', 'poel'):

1983 if self.modelling_code_id.startswith(code):

1984 return True

1985 return False

1986

1987 def is_dynamic(self):

1988 return not self.is_static()

1989

1990 def get_source_depths(self):

1991 raise NotImplementedError('must be implemented in subclass')

1992

1993 def get_tabulated_phase(self, phase_id):

1994 '''

1995 Get tabulated phase definition.

1996 '''

1997

1998 for pdef in self.tabulated_phases:

1999 if pdef.id == phase_id:

2000 return pdef

2001

2002 raise StoreError('No such phase: %s' % phase_id)

2003

2004 def fix_ttt_holes(self, sptree, mode):

2005 raise StoreError(

2006 'Cannot fix travel time table holes in GF stores of type %s.'

2007 % self.short_type)

2008

2009

2010class ConfigTypeA(Config):

2011 '''

2012 Cylindrical symmetry, 1D earth model, single receiver depth

2013

2014 * Problem is invariant to horizontal translations and rotations around

2015 vertical axis.

2016

2017 * All receivers must be at the same depth (e.g. at the surface)

2018 High level index variables: ``(source_depth, distance,

2019 component)``

2020

2021 * The ``distance`` is the surface distance between source and receiver

2022 points.

2023 '''

2024

2025 receiver_depth = Float.T(

2026 default=0.0,

2027 help='Fixed receiver depth [m].')

2028

2029 source_depth_min = Float.T(

2030 help='Minimum source depth [m].')

2031

2032 source_depth_max = Float.T(

2033 help='Maximum source depth [m].')

2034

2035 source_depth_delta = Float.T(

2036 help='Grid spacing of source depths [m]')

2037

2038 distance_min = Float.T(

2039 help='Minimum source-receiver surface distance [m].')

2040

2041 distance_max = Float.T(

2042 help='Maximum source-receiver surface distance [m].')

2043

2044 distance_delta = Float.T(

2045 help='Grid spacing of source-receiver surface distance [m].')

2046

2047 short_type = 'A'

2048

2049 provided_schemes = [

2050 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']

2051

2052 def get_surface_distance(self, args):

2053 return args[1]

2054

2055 def get_distance(self, args):

2056 return math.sqrt(args[0]**2 + args[1]**2)

2057

2058 def get_source_depth(self, args):

2059 return args[0]

2060

2061 def get_source_depths(self):

2062 return self.coords[0]

2063

2064 def get_receiver_depth(self, args):

2065 return self.receiver_depth

2066

2067 def _update(self):

2068 self.mins = num.array(

2069 [self.source_depth_min, self.distance_min], dtype=float)

2070 self.maxs = num.array(

2071 [self.source_depth_max, self.distance_max], dtype=float)

2072 self.deltas = num.array(

2073 [self.source_depth_delta, self.distance_delta],

2074 dtype=float)

2075 self.ns = num.floor((self.maxs - self.mins) / self.deltas +

2076 vicinity_eps).astype(int) + 1

2077 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)

2078 self.deltat = 1.0/self.sample_rate

2079 self.nrecords = num.prod(self.ns) * self.ncomponents

2080 self.coords = tuple(num.linspace(mi, ma, n) for

2081 (mi, ma, n) in

2082 zip(self.mins, self.effective_maxs, self.ns)) + \

2083 (num.arange(self.ncomponents),)

2084

2085 self.nsource_depths, self.ndistances = self.ns

2086

2087 def _make_index_functions(self):

2088

2089 amin, bmin = self.mins

2090 da, db = self.deltas

2091 na, nb = self.ns

2092

2093 ng = self.ncomponents

2094

2095 def index_function(a, b, ig):

2096 ia = int(round((a - amin) / da))

2097 ib = int(round((b - bmin) / db))

2098 try:

2099 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))

2100 except ValueError:

2101 raise OutOfBounds()

2102

2103 def indices_function(a, b, ig):

2104 ia = num.round((a - amin) / da).astype(int)

2105 ib = num.round((b - bmin) / db).astype(int)

2106 try:

2107 return num.ravel_multi_index((ia, ib, ig), (na, nb, ng))

2108 except ValueError:

2109 for ia_, ib_, ig_ in zip(ia, ib, ig):

2110 try:

2111 num.ravel_multi_index((ia_, ib_, ig_), (na, nb, ng))

2112 except ValueError:

2113 raise OutOfBounds()

2114

2115 def grid_interpolation_coefficients(a, b):

2116 ias = indi12((a - amin) / da, na)

2117 ibs = indi12((b - bmin) / db, nb)

2118 return ias, ibs

2119

2120 def vicinity_function(a, b, ig):

2121 ias, ibs = grid_interpolation_coefficients(a, b)

2122

2123 if not (0 <= ig < ng):

2124 raise OutOfBounds()

2125

2126 indis = []

2127 weights = []

2128 for ia, va in ias:

2129 iia = ia*nb*ng

2130 for ib, vb in ibs:

2131 indis.append(iia + ib*ng + ig)

2132 weights.append(va*vb)

2133

2134 return num.array(indis), num.array(weights)

2135

2136 def vicinities_function(a, b, ig):

2137

2138 xa = (a - amin) / da

2139 xb = (b - bmin) / db

2140

2141 xa_fl = num.floor(xa)

2142 xa_ce = num.ceil(xa)

2143 xb_fl = num.floor(xb)

2144 xb_ce = num.ceil(xb)

2145 va_fl = 1.0 - (xa - xa_fl)

2146 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)

2147 vb_fl = 1.0 - (xb - xb_fl)

2148 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)

2149

2150 ia_fl = xa_fl.astype(int)

2151 ia_ce = xa_ce.astype(int)

2152 ib_fl = xb_fl.astype(int)

2153 ib_ce = xb_ce.astype(int)

2154

2155 if num.any(ia_fl < 0) or num.any(ia_fl >= na):

2156 raise OutOfBounds()

2157

2158 if num.any(ia_ce < 0) or num.any(ia_ce >= na):

2159 raise OutOfBounds()

2160

2161 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):

2162 raise OutOfBounds()

2163

2164 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):

2165 raise OutOfBounds()

2166

2167 irecords = num.empty(a.size*4, dtype=int)

2168 irecords[0::4] = ia_fl*nb*ng + ib_fl*ng + ig

2169 irecords[1::4] = ia_ce*nb*ng + ib_fl*ng + ig

2170 irecords[2::4] = ia_fl*nb*ng + ib_ce*ng + ig

2171 irecords[3::4] = ia_ce*nb*ng + ib_ce*ng + ig

2172

2173 weights = num.empty(a.size*4, dtype=float)

2174 weights[0::4] = va_fl * vb_fl

2175 weights[1::4] = va_ce * vb_fl

2176 weights[2::4] = va_fl * vb_ce

2177 weights[3::4] = va_ce * vb_ce

2178

2179 return irecords, weights

2180

2181 self._index_function = index_function

2182 self._indices_function = indices_function

2183 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2184 self._vicinity_function = vicinity_function

2185 self._vicinities_function = vicinities_function

2186

2187 def make_indexing_args(self, source, receiver, icomponents):

2188 nc = icomponents.size

2189 dists = source.distances_to(receiver)

2190 n = dists.size

2191 return (num.tile(source.depths, nc//n),

2192 num.tile(dists, nc//n),

2193 icomponents)

2194

2195 def make_indexing_args1(self, source, receiver):

2196 return (source.depth, source.distance_to(receiver))

2197

2198 @property

2199 def short_extent(self):

2200 return '%g:%g:%g x %g:%g:%g' % (

2201 self.source_depth_min/km,

2202 self.source_depth_max/km,

2203 self.source_depth_delta/km,

2204 self.distance_min/km,

2205 self.distance_max/km,

2206 self.distance_delta/km)

2207

2208 def fix_ttt_holes(self, sptree, mode):

2209 from pyrocko import eikonal_ext, spit

2210

2211 nodes = self.nodes(level=-1)

2212

2213 delta = self.deltas[-1]

2214 assert num.all(delta == self.deltas)

2215

2216 nsources, ndistances = self.ns

2217

2218 points = num.zeros((nodes.shape[0], 3))

2219 points[:, 0] = nodes[:, 1]

2220 points[:, 2] = nodes[:, 0]

2221

2222 speeds = self.get_material_property(

2223 0., 0., points,

2224 parameter='vp' if mode == cake.P else 'vs',

2225 interpolation='multilinear')

2226

2227 speeds = speeds.reshape((nsources, ndistances))

2228

2229 times = sptree.interpolate_many(nodes)

2230

2231 times[num.isnan(times)] = -1.

2232 times = times.reshape(speeds.shape)

2233

2234 try:

2235 eikonal_ext.eikonal_solver_fmm_cartesian(

2236 speeds, times, delta)

2237 except eikonal_ext.EikonalExtError as e:

2238 if str(e).endswith('please check results'):

2239 logger.debug(

2240 'Got a warning from eikonal solver '

2241 '- may be ok...')

2242 else:

2243 raise

2244

2245 def func(x):

2246 ibs, ics = \

2247 self.grid_interpolation_coefficients(*x)

2248

2249 t = 0

2250 for ib, vb in ibs:

2251 for ic, vc in ics:

2252 t += times[ib, ic] * vb * vc

2253

2254 return t

2255

2256 return spit.SPTree(

2257 f=func,

2258 ftol=sptree.ftol,

2259 xbounds=sptree.xbounds,

2260 xtols=sptree.xtols)

2261

2262

2263class ConfigTypeB(Config):

2264 '''

2265 Cylindrical symmetry, 1D earth model, variable receiver depth

2266

2267 * Symmetries like in :py:class:`ConfigTypeA` but has additional index for

2268 receiver depth

2269

2270 * High level index variables: ``(receiver_depth, source_depth,

2271 receiver_distance, component)``

2272 '''

2273

2274 receiver_depth_min = Float.T(

2275 help='Minimum receiver depth [m].')

2276

2277 receiver_depth_max = Float.T(

2278 help='Maximum receiver depth [m].')

2279

2280 receiver_depth_delta = Float.T(

2281 help='Grid spacing of receiver depths [m]')

2282

2283 source_depth_min = Float.T(

2284 help='Minimum source depth [m].')

2285

2286 source_depth_max = Float.T(

2287 help='Maximum source depth [m].')

2288

2289 source_depth_delta = Float.T(

2290 help='Grid spacing of source depths [m]')

2291

2292 distance_min = Float.T(

2293 help='Minimum source-receiver surface distance [m].')

2294

2295 distance_max = Float.T(

2296 help='Maximum source-receiver surface distance [m].')

2297

2298 distance_delta = Float.T(

2299 help='Grid spacing of source-receiver surface distances [m].')

2300

2301 short_type = 'B'

2302

2303 provided_schemes = [

2304 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10']

2305

2306 def get_distance(self, args):

2307 return math.sqrt((args[1] - args[0])**2 + args[2]**2)

2308

2309 def get_surface_distance(self, args):

2310 return args[2]

2311

2312 def get_source_depth(self, args):

2313 return args[1]

2314

2315 def get_receiver_depth(self, args):

2316 return args[0]

2317

2318 def get_source_depths(self):

2319 return self.coords[1]

2320

2321 def _update(self):

2322 self.mins = num.array([

2323 self.receiver_depth_min,

2324 self.source_depth_min,

2325 self.distance_min],

2326 dtype=float)

2327

2328 self.maxs = num.array([

2329 self.receiver_depth_max,

2330 self.source_depth_max,

2331 self.distance_max],

2332 dtype=float)

2333

2334 self.deltas = num.array([

2335 self.receiver_depth_delta,

2336 self.source_depth_delta,

2337 self.distance_delta],

2338 dtype=float)

2339

2340 self.ns = num.floor((self.maxs - self.mins) / self.deltas +

2341 vicinity_eps).astype(int) + 1

2342 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)

2343 self.deltat = 1.0/self.sample_rate

2344 self.nrecords = num.prod(self.ns) * self.ncomponents

2345 self.coords = tuple(num.linspace(mi, ma, n) for

2346 (mi, ma, n) in

2347 zip(self.mins, self.effective_maxs, self.ns)) + \

2348 (num.arange(self.ncomponents),)

2349 self.nreceiver_depths, self.nsource_depths, self.ndistances = self.ns

2350

2351 def _make_index_functions(self):

2352

2353 amin, bmin, cmin = self.mins

2354 da, db, dc = self.deltas

2355 na, nb, nc = self.ns

2356 ng = self.ncomponents

2357

2358 def index_function(a, b, c, ig):

2359 ia = int(round((a - amin) / da))

2360 ib = int(round((b - bmin) / db))

2361 ic = int(round((c - cmin) / dc))

2362 try:

2363 return num.ravel_multi_index((ia, ib, ic, ig),

2364 (na, nb, nc, ng))

2365 except ValueError:

2366 raise OutOfBounds()

2367

2368 def indices_function(a, b, c, ig):

2369 ia = num.round((a - amin) / da).astype(int)

2370 ib = num.round((b - bmin) / db).astype(int)

2371 ic = num.round((c - cmin) / dc).astype(int)

2372 try:

2373 return num.ravel_multi_index((ia, ib, ic, ig),

2374 (na, nb, nc, ng))

2375 except ValueError:

2376 raise OutOfBounds()

2377

2378 def grid_interpolation_coefficients(a, b, c):

2379 ias = indi12((a - amin) / da, na)

2380 ibs = indi12((b - bmin) / db, nb)

2381 ics = indi12((c - cmin) / dc, nc)

2382 return ias, ibs, ics

2383

2384 def vicinity_function(a, b, c, ig):

2385 ias, ibs, ics = grid_interpolation_coefficients(a, b, c)

2386

2387 if not (0 <= ig < ng):

2388 raise OutOfBounds()

2389

2390 indis = []

2391 weights = []

2392 for ia, va in ias:

2393 iia = ia*nb*nc*ng

2394 for ib, vb in ibs:

2395 iib = ib*nc*ng

2396 for ic, vc in ics:

2397 indis.append(iia + iib + ic*ng + ig)

2398 weights.append(va*vb*vc)

2399

2400 return num.array(indis), num.array(weights)

2401

2402 def vicinities_function(a, b, c, ig):

2403

2404 xa = (a - amin) / da

2405 xb = (b - bmin) / db

2406 xc = (c - cmin) / dc

2407

2408 xa_fl = num.floor(xa)

2409 xa_ce = num.ceil(xa)

2410 xb_fl = num.floor(xb)

2411 xb_ce = num.ceil(xb)

2412 xc_fl = num.floor(xc)

2413 xc_ce = num.ceil(xc)

2414 va_fl = 1.0 - (xa - xa_fl)

2415 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)

2416 vb_fl = 1.0 - (xb - xb_fl)

2417 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)

2418 vc_fl = 1.0 - (xc - xc_fl)

2419 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)

2420

2421 ia_fl = xa_fl.astype(int)

2422 ia_ce = xa_ce.astype(int)

2423 ib_fl = xb_fl.astype(int)

2424 ib_ce = xb_ce.astype(int)

2425 ic_fl = xc_fl.astype(int)

2426 ic_ce = xc_ce.astype(int)

2427

2428 if num.any(ia_fl < 0) or num.any(ia_fl >= na):

2429 raise OutOfBounds()

2430

2431 if num.any(ia_ce < 0) or num.any(ia_ce >= na):

2432 raise OutOfBounds()

2433

2434 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):

2435 raise OutOfBounds()

2436

2437 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):

2438 raise OutOfBounds()

2439

2440 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):

2441 raise OutOfBounds()

2442

2443 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):

2444 raise OutOfBounds()

2445

2446 irecords = num.empty(a.size*8, dtype=int)

2447 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig

2448 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig

2449 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig

2450 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig

2451 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig

2452 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig

2453 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig

2454 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig

2455

2456 weights = num.empty(a.size*8, dtype=float)

2457 weights[0::8] = va_fl * vb_fl * vc_fl

2458 weights[1::8] = va_ce * vb_fl * vc_fl

2459 weights[2::8] = va_fl * vb_ce * vc_fl

2460 weights[3::8] = va_ce * vb_ce * vc_fl

2461 weights[4::8] = va_fl * vb_fl * vc_ce

2462 weights[5::8] = va_ce * vb_fl * vc_ce

2463 weights[6::8] = va_fl * vb_ce * vc_ce

2464 weights[7::8] = va_ce * vb_ce * vc_ce

2465

2466 return irecords, weights

2467

2468 self._index_function = index_function

2469 self._indices_function = indices_function

2470 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2471 self._vicinity_function = vicinity_function

2472 self._vicinities_function = vicinities_function

2473

2474 def make_indexing_args(self, source, receiver, icomponents):

2475 nc = icomponents.size

2476 dists = source.distances_to(receiver)

2477 n = dists.size

2478 receiver_depths = num.empty(nc)

2479 receiver_depths.fill(receiver.depth)

2480 return (receiver_depths,

2481 num.tile(source.depths, nc//n),

2482 num.tile(dists, nc//n),

2483 icomponents)

2484

2485 def make_indexing_args1(self, source, receiver):

2486 return (receiver.depth,

2487 source.depth,

2488 source.distance_to(receiver))

2489

2490 @property

2491 def short_extent(self):

2492 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2493 self.receiver_depth_min/km,

2494 self.receiver_depth_max/km,

2495 self.receiver_depth_delta/km,

2496 self.source_depth_min/km,

2497 self.source_depth_max/km,

2498 self.source_depth_delta/km,

2499 self.distance_min/km,

2500 self.distance_max/km,

2501 self.distance_delta/km)

2502

2503 def fix_ttt_holes(self, sptree, mode):

2504 from pyrocko import eikonal_ext, spit

2505

2506 nodes_sr = self.nodes(minlevel=1, level=-1)

2507

2508 delta = self.deltas[-1]

2509 assert num.all(delta == self.deltas[1:])

2510

2511 nreceivers, nsources, ndistances = self.ns

2512

2513 points = num.zeros((nodes_sr.shape[0], 3))

2514 points[:, 0] = nodes_sr[:, 1]

2515 points[:, 2] = nodes_sr[:, 0]

2516

2517 speeds = self.get_material_property(

2518 0., 0., points,

2519 parameter='vp' if mode == cake.P else 'vs',

2520 interpolation='multilinear')

2521

2522 speeds = speeds.reshape((nsources, ndistances))

2523

2524 receiver_times = []

2525 for ireceiver in range(nreceivers):

2526 nodes = num.hstack([

2527 num_full(

2528 (nodes_sr.shape[0], 1),

2529 self.coords[0][ireceiver]),

2530 nodes_sr])

2531

2532 times = sptree.interpolate_many(nodes)

2533

2534 times[num.isnan(times)] = -1.

2535

2536 times = times.reshape(speeds.shape)

2537

2538 try:

2539 eikonal_ext.eikonal_solver_fmm_cartesian(

2540 speeds, times, delta)

2541 except eikonal_ext.EikonalExtError as e:

2542 if str(e).endswith('please check results'):

2543 logger.debug(

2544 'Got a warning from eikonal solver '

2545 '- may be ok...')

2546 else:

2547 raise

2548

2549 receiver_times.append(times)

2550

2551 def func(x):

2552 ias, ibs, ics = \

2553 self.grid_interpolation_coefficients(*x)

2554

2555 t = 0

2556 for ia, va in ias:

2557 times = receiver_times[ia]

2558 for ib, vb in ibs:

2559 for ic, vc in ics:

2560 t += times[ib, ic] * va * vb * vc

2561

2562 return t

2563

2564 return spit.SPTree(

2565 f=func,

2566 ftol=sptree.ftol,

2567 xbounds=sptree.xbounds,

2568 xtols=sptree.xtols)

2569

2570

2571class ConfigTypeC(Config):

2572 '''

2573 No symmetrical constraints, one fixed receiver position.

2574

2575 * Cartesian 3D source volume around a reference point

2576

2577 * High level index variables: ``(source_depth,

2578 source_east_shift, source_north_shift, component)``

2579 '''

2580

2581 receiver = Receiver.T(

2582 help='Receiver location')

2583

2584 source_origin = Location.T(

2585 help='Origin of the source volume grid.')

2586

2587 source_depth_min = Float.T(

2588 help='Minimum source depth [m].')

2589

2590 source_depth_max = Float.T(

2591 help='Maximum source depth [m].')

2592

2593 source_depth_delta = Float.T(

2594 help='Source depth grid spacing [m].')

2595

2596 source_east_shift_min = Float.T(

2597 help='Minimum easting of source grid [m].')

2598

2599 source_east_shift_max = Float.T(

2600 help='Maximum easting of source grid [m].')

2601

2602 source_east_shift_delta = Float.T(

2603 help='Source volume grid spacing in east direction [m].')

2604

2605 source_north_shift_min = Float.T(

2606 help='Minimum northing of source grid [m].')

2607

2608 source_north_shift_max = Float.T(

2609 help='Maximum northing of source grid [m].')

2610

2611 source_north_shift_delta = Float.T(

2612 help='Source volume grid spacing in north direction [m].')

2613

2614 short_type = 'C'

2615

2616 provided_schemes = ['elastic18']

2617

2618 def get_surface_distance(self, args):

2619 _, source_east_shift, source_north_shift, _ = args

2620 sorig = self.source_origin

2621 sloc = Location(

2622 lat=sorig.lat,

2623 lon=sorig.lon,

2624 north_shift=sorig.north_shift + source_north_shift,

2625 east_shift=sorig.east_shift + source_east_shift)

2626

2627 return self.receiver.distance_to(sloc)

2628

2629 def get_distance(self, args):

2630 sdepth, source_east_shift, source_north_shift, _ = args

2631 sorig = self.source_origin

2632 sloc = Location(

2633 lat=sorig.lat,

2634 lon=sorig.lon,

2635 north_shift=sorig.north_shift + source_north_shift,

2636 east_shift=sorig.east_shift + source_east_shift)

2637

2638 return self.receiver.distance_3d_to(sloc)

2639

2640 def get_source_depth(self, args):

2641 return args[0]

2642

2643 def get_receiver_depth(self, args):

2644 return self.receiver.depth

2645

2646 def get_source_depths(self):

2647 return self.coords[0]

2648

2649 def _update(self):

2650 self.mins = num.array([

2651 self.source_depth_min,

2652 self.source_east_shift_min,

2653 self.source_north_shift_min],

2654 dtype=float)

2655

2656 self.maxs = num.array([

2657 self.source_depth_max,

2658 self.source_east_shift_max,

2659 self.source_north_shift_max],

2660 dtype=float)

2661

2662 self.deltas = num.array([

2663 self.source_depth_delta,

2664 self.source_east_shift_delta,

2665 self.source_north_shift_delta],

2666 dtype=float)

2667

2668 self.ns = num.floor((self.maxs - self.mins) / self.deltas +

2669 vicinity_eps).astype(int) + 1

2670 self.effective_maxs = self.mins + self.deltas * (self.ns - 1)

2671 self.deltat = 1.0/self.sample_rate

2672 self.nrecords = num.prod(self.ns) * self.ncomponents

2673

2674 self.coords = tuple(

2675 num.linspace(mi, ma, n) for (mi, ma, n) in

2676 zip(self.mins, self.effective_maxs, self.ns)) + \

2677 (num.arange(self.ncomponents),)

2678

2679 def _make_index_functions(self):

2680

2681 amin, bmin, cmin = self.mins

2682 da, db, dc = self.deltas

2683 na, nb, nc = self.ns

2684 ng = self.ncomponents

2685

2686 def index_function(a, b, c, ig):

2687 ia = int(round((a - amin) / da))

2688 ib = int(round((b - bmin) / db))

2689 ic = int(round((c - cmin) / dc))

2690 try:

2691 return num.ravel_multi_index((ia, ib, ic, ig),

2692 (na, nb, nc, ng))

2693 except ValueError:

2694 raise OutOfBounds(values=(a, b, c, ig))

2695

2696 def indices_function(a, b, c, ig):

2697 ia = num.round((a - amin) / da).astype(int)

2698 ib = num.round((b - bmin) / db).astype(int)

2699 ic = num.round((c - cmin) / dc).astype(int)

2700

2701 try:

2702 return num.ravel_multi_index((ia, ib, ic, ig),

2703 (na, nb, nc, ng))

2704 except ValueError:

2705 raise OutOfBounds()

2706

2707 def vicinity_function(a, b, c, ig):

2708 ias = indi12((a - amin) / da, na)

2709 ibs = indi12((b - bmin) / db, nb)

2710 ics = indi12((c - cmin) / dc, nc)

2711

2712 if not (0 <= ig < ng):

2713 raise OutOfBounds()

2714

2715 indis = []

2716 weights = []

2717 for ia, va in ias:

2718 iia = ia*nb*nc*ng

2719 for ib, vb in ibs:

2720 iib = ib*nc*ng

2721 for ic, vc in ics:

2722 indis.append(iia + iib + ic*ng + ig)

2723 weights.append(va*vb*vc)

2724

2725 return num.array(indis), num.array(weights)

2726

2727 def vicinities_function(a, b, c, ig):

2728

2729 xa = (a-amin) / da

2730 xb = (b-bmin) / db

2731 xc = (c-cmin) / dc

2732

2733 xa_fl = num.floor(xa)

2734 xa_ce = num.ceil(xa)

2735 xb_fl = num.floor(xb)

2736 xb_ce = num.ceil(xb)

2737 xc_fl = num.floor(xc)

2738 xc_ce = num.ceil(xc)

2739 va_fl = 1.0 - (xa - xa_fl)

2740 va_ce = (1.0 - (xa_ce - xa)) * (xa_ce - xa_fl)

2741 vb_fl = 1.0 - (xb - xb_fl)

2742 vb_ce = (1.0 - (xb_ce - xb)) * (xb_ce - xb_fl)

2743 vc_fl = 1.0 - (xc - xc_fl)

2744 vc_ce = (1.0 - (xc_ce - xc)) * (xc_ce - xc_fl)

2745

2746 ia_fl = xa_fl.astype(int)

2747 ia_ce = xa_ce.astype(int)

2748 ib_fl = xb_fl.astype(int)

2749 ib_ce = xb_ce.astype(int)

2750 ic_fl = xc_fl.astype(int)

2751 ic_ce = xc_ce.astype(int)

2752

2753 if num.any(ia_fl < 0) or num.any(ia_fl >= na):

2754 raise OutOfBounds()

2755

2756 if num.any(ia_ce < 0) or num.any(ia_ce >= na):

2757 raise OutOfBounds()

2758

2759 if num.any(ib_fl < 0) or num.any(ib_fl >= nb):

2760 raise OutOfBounds()

2761

2762 if num.any(ib_ce < 0) or num.any(ib_ce >= nb):

2763 raise OutOfBounds()

2764

2765 if num.any(ic_fl < 0) or num.any(ic_fl >= nc):

2766 raise OutOfBounds()

2767

2768 if num.any(ic_ce < 0) or num.any(ic_ce >= nc):

2769 raise OutOfBounds()

2770

2771 irecords = num.empty(a.size*8, dtype=int)

2772 irecords[0::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig

2773 irecords[1::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_fl*ng + ig

2774 irecords[2::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig

2775 irecords[3::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_fl*ng + ig

2776 irecords[4::8] = ia_fl*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig

2777 irecords[5::8] = ia_ce*nb*nc*ng + ib_fl*nc*ng + ic_ce*ng + ig

2778 irecords[6::8] = ia_fl*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig

2779 irecords[7::8] = ia_ce*nb*nc*ng + ib_ce*nc*ng + ic_ce*ng + ig

2780

2781 weights = num.empty(a.size*8, dtype=float)

2782 weights[0::8] = va_fl * vb_fl * vc_fl

2783 weights[1::8] = va_ce * vb_fl * vc_fl

2784 weights[2::8] = va_fl * vb_ce * vc_fl

2785 weights[3::8] = va_ce * vb_ce * vc_fl

2786 weights[4::8] = va_fl * vb_fl * vc_ce

2787 weights[5::8] = va_ce * vb_fl * vc_ce

2788 weights[6::8] = va_fl * vb_ce * vc_ce

2789 weights[7::8] = va_ce * vb_ce * vc_ce

2790

2791 return irecords, weights

2792

2793 self._index_function = index_function

2794 self._indices_function = indices_function

2795 self._vicinity_function = vicinity_function

2796 self._vicinities_function = vicinities_function

2797

2798 def make_indexing_args(self, source, receiver, icomponents):

2799 nc = icomponents.size

2800

2801 dists = source.distances_to(self.source_origin)

2802 azis, _ = source.azibazis_to(self.source_origin)

2803

2804 source_north_shifts = - num.cos(d2r*azis) * dists

2805 source_east_shifts = - num.sin(d2r*azis) * dists

2806 source_depths = source.depths - self.source_origin.depth

2807

2808 n = dists.size

2809

2810 return (num.tile(source_depths, nc//n),

2811 num.tile(source_east_shifts, nc//n),

2812 num.tile(source_north_shifts, nc//n),

2813 icomponents)

2814

2815 def make_indexing_args1(self, source, receiver):

2816 dist = source.distance_to(self.source_origin)

2817 azi, _ = source.azibazi_to(self.source_origin)

2818

2819 source_north_shift = - num.cos(d2r*azi) * dist

2820 source_east_shift = - num.sin(d2r*azi) * dist

2821 source_depth = source.depth - self.source_origin.depth

2822

2823 return (source_depth,

2824 source_east_shift,

2825 source_north_shift)

2826

2827 @property

2828 def short_extent(self):

2829 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2830 self.source_depth_min/km,

2831 self.source_depth_max/km,

2832 self.source_depth_delta/km,

2833 self.source_east_shift_min/km,

2834 self.source_east_shift_max/km,

2835 self.source_east_shift_delta/km,

2836 self.source_north_shift_min/km,

2837 self.source_north_shift_max/km,

2838 self.source_north_shift_delta/km)

2839

2840

2841class Weighting(Object):

2842 factor = Float.T(default=1.0)

2843

2844

2845class Taper(Object):

2846 tmin = Timing.T()

2847 tmax = Timing.T()

2848 tfade = Float.T(default=0.0)

2849 shape = StringChoice.T(

2850 choices=['cos', 'linear'],

2851 default='cos',

2852 optional=True)

2853

2854

2855class SimplePattern(SObject):

2856

2857 _pool = {}

2858

2859 def __init__(self, pattern):

2860 self._pattern = pattern

2861 SObject.__init__(self)

2862

2863 def __str__(self):

2864 return self._pattern

2865

2866 @property

2867 def regex(self):

2868 pool = SimplePattern._pool

2869 if self.pattern not in pool:

2870 rpat = '|'.join(fnmatch.translate(x) for

2871 x in self.pattern.split('|'))

2872 pool[self.pattern] = re.compile(rpat, re.I)

2873

2874 return pool[self.pattern]

2875

2876 def match(self, s):

2877 return self.regex.match(s)

2878

2879

2880class WaveformType(StringChoice):

2881 choices = ['dis', 'vel', 'acc',

2882 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',

2883 'envelope_dis', 'envelope_vel', 'envelope_acc']

2884

2885

2886class ChannelSelection(Object):

2887 pattern = SimplePattern.T(optional=True)

2888 min_sample_rate = Float.T(optional=True)

2889 max_sample_rate = Float.T(optional=True)

2890

2891

2892class StationSelection(Object):

2893 includes = SimplePattern.T()

2894 excludes = SimplePattern.T()

2895 distance_min = Float.T(optional=True)

2896 distance_max = Float.T(optional=True)

2897 azimuth_min = Float.T(optional=True)

2898 azimuth_max = Float.T(optional=True)

2899

2900

2901class WaveformSelection(Object):

2902 channel_selection = ChannelSelection.T(optional=True)

2903 station_selection = StationSelection.T(optional=True)

2904 taper = Taper.T()

2905 # filter = FrequencyResponse.T()

2906 waveform_type = WaveformType.T(default='dis')

2907 weighting = Weighting.T(optional=True)

2908 sample_rate = Float.T(optional=True)

2909 gf_store_id = StringID.T(optional=True)

2910

2911

2912def indi12(x, n):

2913 '''

2914 Get linear interpolation index and weight.

2915 '''

2916

2917 r = round(x)

2918 if abs(r - x) < vicinity_eps:

2919 i = int(r)

2920 if not (0 <= i < n):

2921 raise OutOfBounds()

2922

2923 return ((int(r), 1.),)

2924 else:

2925 f = math.floor(x)

2926 i = int(f)

2927 if not (0 <= i < n-1):

2928 raise OutOfBounds()

2929

2930 v = x-f

2931 return ((i, 1.-v), (i + 1, v))

2932

2933

2934def float_or_none(s):

2935 units = {

2936 'k': 1e3,

2937 'M': 1e6,

2938 }

2939

2940 factor = 1.0

2941 if s and s[-1] in units:

2942 factor = units[s[-1]]

2943 s = s[:-1]

2944 if not s:

2945 raise ValueError("unit without a number: '%s'" % s)

2946

2947 if s:

2948 return float(s) * factor

2949 else:

2950 return None

2951

2952

2953class GridSpecError(Exception):

2954 def __init__(self, s):

2955 Exception.__init__(self, 'invalid grid specification: %s' % s)

2956

2957

2958def parse_grid_spec(spec):

2959 try:

2960 result = []

2961 for dspec in spec.split(','):

2962 t = dspec.split('@')

2963 num = start = stop = step = None

2964 if len(t) == 2:

2965 num = int(t[1])

2966 if num <= 0:

2967 raise GridSpecError(spec)

2968

2969 elif len(t) > 2:

2970 raise GridSpecError(spec)

2971

2972 s = t[0]

2973 v = [float_or_none(x) for x in s.split(':')]

2974 if len(v) == 1:

2975 start = stop = v[0]

2976 if len(v) >= 2:

2977 start, stop = v[0:2]

2978 if len(v) == 3:

2979 step = v[2]

2980

2981 if len(v) > 3 or (len(v) > 2 and num is not None):

2982 raise GridSpecError(spec)

2983

2984 if step == 0.0:

2985 raise GridSpecError(spec)

2986

2987 result.append((start, stop, step, num))

2988

2989 except ValueError:

2990 raise GridSpecError(spec)

2991

2992 return result

2993

2994

2995def start_stop_num(start, stop, step, num, mi, ma, inc, eps=1e-5):

2996 swap = step is not None and step < 0.

2997 if start is None:

2998 start = ma if swap else mi

2999 if stop is None:

3000 stop = mi if swap else ma

3001 if step is None:

3002 step = -inc if ma < mi else inc

3003 if num is None:

3004 if (step < 0) != (stop-start < 0):

3005 raise GridSpecError()

3006

3007 num = int(round((stop-start)/step))+1

3008 stop2 = start + (num-1)*step

3009 if abs(stop-stop2) > eps:

3010 num = int(math.floor((stop-start)/step))+1

3011 stop = start + (num-1)*step

3012 else:

3013 stop = stop2

3014

3015 if start == stop:

3016 num = 1

3017

3018 return start, stop, num

3019

3020

3021def nditer_outer(x):

3022 return num.nditer(

3023 x, op_axes=(num.identity(len(x), dtype=int)-1).tolist())

3024

3025

3026def nodes(xs):

3027 ns = [x.size for x in xs]

3028 nnodes = num.prod(ns)

3029 ndim = len(xs)

3030 nodes = num.empty((nnodes, ndim), dtype=xs[0].dtype)

3031 for idim in range(ndim-1, -1, -1):

3032 x = xs[idim]

3033 nrepeat = num.prod(ns[idim+1:], dtype=int)

3034 ntile = num.prod(ns[:idim], dtype=int)

3035 nodes[:, idim] = num.repeat(num.tile(x, ntile), nrepeat)

3036

3037 return nodes

3038

3039

3040def filledi(x, n):

3041 a = num.empty(n, dtype=int)

3042 a.fill(x)

3043 return a

3044

3045

3046config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]

3047

3048discretized_source_classes = [

3049 DiscretizedExplosionSource,

3050 DiscretizedSFSource,

3051 DiscretizedMTSource,

3052 DiscretizedPorePressureSource]

3053

3054

3055__all__ = '''

3056Earthmodel1D

3057StringID

3058ScopeType

3059WaveformType

3060QuantityType

3061NearfieldTermsType

3062Reference

3063Region

3064CircularRegion

3065RectangularRegion

3066PhaseSelect

3067InvalidTimingSpecification

3068Timing

3069TPDef

3070OutOfBounds

3071Location

3072Receiver

3073'''.split() + [

3074 S.__name__ for S in discretized_source_classes + config_type_classes] + '''

3075ComponentScheme

3076component_scheme_to_description

3077component_schemes

3078Config

3079GridSpecError

3080Weighting

3081Taper

3082SimplePattern

3083WaveformType

3084ChannelSelection

3085StationSelection

3086WaveformSelection

3087nditer_outer

3088dump

3089load

3090discretized_source_classes

3091config_type_classes

3092UnavailableScheme

3093InterpolationMethod

3094SeismosizerTrace

3095SeismosizerResult

3096Result

3097StaticResult

3098TimingAttributeError

3099'''.split()