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

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

27from .error import StoreError

29try:

30 new_str = unicode

31except NameError:

32 new_str = str

34guts_prefix = 'pf'

36d2r = math.pi / 180.

37r2d = 1.0 / d2r

38km = 1000.

39vicinity_eps = 1e-5

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

44def fmt_choices(cls):

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

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

49class InterpolationMethod(StringChoice):

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

53class SeismosizerTrace(Object):

55 codes = Tuple.T(

56 4, String.T(),

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

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

60 data = Array.T(

61 shape=(None,),

62 dtype=num.float32,

63 serialize_as='base64',

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

65 help='numpy array with data samples')

67 deltat = Float.T(

68 default=1.0,

69 help='sampling interval [s]')

71 tmin = Timestamp.T(

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

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

75 def pyrocko_trace(self):

76 c = self.codes

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

78 ydata=self.data,

79 deltat=self.deltat,

80 tmin=self.tmin)

82 @classmethod

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

84 d = dict(

85 codes=tr.nslc_id,

86 tmin=tr.tmin,

87 deltat=tr.deltat,

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

89 d.update(kwargs)

90 return cls(**d)

93class SeismosizerResult(Object):

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

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

98class Result(SeismosizerResult):

99 trace = SeismosizerTrace.T(optional=True)

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

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

102

103

104class StaticResult(SeismosizerResult):

105 result = Dict.T(

106 String.T(),

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

108

109

110class GNSSCampaignResult(StaticResult):

111 campaign = gnss.GNSSCampaign.T(

112 optional=True)

113

114

115class SatelliteResult(StaticResult):

116

117 theta = Array.T(

118 optional=True,

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

120

121 phi = Array.T(

122 optional=True,

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

124

125

126class KiteSceneResult(SatelliteResult):

127

128 shape = Tuple.T()

129

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

131 try:

132 from kite import Scene

133 except ImportError:

134 raise ImportError('Kite not installed')

135

136 def reshape(arr):

137 return arr.reshape(self.shape)

138

139 displacement = self.result[component]

140

141 displacement, theta, phi = map(

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

143

144 sc = Scene(

145 displacement=displacement,

146 theta=theta, phi=phi,

147 config=self.config)

148

149 return sc

150

151

152class ComponentSchemeDescription(Object):

153 name = String.T()

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

155 ncomponents = Int.T()

156 default_stored_quantity = String.T(optional=True)

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

158

159

160component_scheme_descriptions = [

161 ComponentSchemeDescription(

162 name='elastic2',

163 source_terms=['m0'],

164 ncomponents=2,

165 default_stored_quantity='displacement',

166 provided_components=[

167 'n', 'e', 'd']),

168 ComponentSchemeDescription(

169 name='elastic8',

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

171 ncomponents=8,

172 default_stored_quantity='displacement',

173 provided_components=[

174 'n', 'e', 'd']),

175 ComponentSchemeDescription(

176 name='elastic10',

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

178 ncomponents=10,

179 default_stored_quantity='displacement',

180 provided_components=[

181 'n', 'e', 'd']),

182 ComponentSchemeDescription(

183 name='elastic10_fd',

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

185 ncomponents=30,

186 default_stored_quantity='displacement',

187 provided_components=[

188 'n', 'e', 'd']),

189 ComponentSchemeDescription(

190 name='elastic18',

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

192 ncomponents=18,

193 default_stored_quantity='displacement',

194 provided_components=[

195 'n', 'e', 'd']),

196 ComponentSchemeDescription(

197 name='elastic5',

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

199 ncomponents=5,

200 default_stored_quantity='displacement',

201 provided_components=[

202 'n', 'e', 'd']),

203 ComponentSchemeDescription(

204 name='rotational8',

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

206 ncomponents=8,

207 default_stored_quantity='rotation_displacement',

208 provided_components=[

209 'n', 'e', 'd']),

210 ComponentSchemeDescription(

211 name='poroelastic10',

212 source_terms=['pore_pressure'],

213 ncomponents=10,

214 default_stored_quantity=None,

215 provided_components=[

216 'displacement.n', 'displacement.e', 'displacement.d',

217 'vertical_tilt.n', 'vertical_tilt.e',

218 'pore_pressure',

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

220

221

222# new names?

223# 'mt_to_displacement_1d'

224# 'mt_to_displacement_1d_ff_only'

225# 'mt_to_gravity_1d'

226# 'mt_to_stress_1d'

227# 'explosion_to_displacement_1d'

228# 'sf_to_displacement_1d'

229# 'mt_to_displacement_3d'

230# 'mt_to_gravity_3d'

231# 'mt_to_stress_3d'

232# 'pore_pressure_to_displacement_1d'

233# 'pore_pressure_to_vertical_tilt_1d'

234# 'pore_pressure_to_pore_pressure_1d'

235# 'pore_pressure_to_darcy_velocity_1d'

236

237

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

239component_scheme_to_description = dict(

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

241

242

243class ComponentScheme(StringChoice):

244 '''

245 Different Green's Function component schemes are available:

246

247 ================= =========================================================

248 Name Description

249 ================= =========================================================

250 ``elastic10`` Elastodynamic for

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

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

253 sources only

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

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

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

257 MT sources only

258 ``elastic2`` Elastodynamic for

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

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

261 isotropic sources only

262 ``elastic5`` Elastodynamic for

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

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

265 sources only

266 ``elastic18`` Elastodynamic for

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

268 sources only

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

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

271 ``rotational8`` Elastodynamic rotational motions for

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

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

274 sources only

275 ================= =========================================================

276 '''

277

278 choices = component_schemes

279

280

281class Earthmodel1D(Object):

282 dummy_for = cake.LayeredModel

283 dummy_for_description = 'pyrocko.cake.LayeredModel'

284

285 class __T(TBase):

286 def regularize_extra(self, val):

287 if isinstance(val, str):

288 val = cake.LayeredModel.from_scanlines(

289 cake.read_nd_model_str(val))

290

291 return val

292

293 def to_save(self, val):

294 return literal(cake.write_nd_model_str(val))

295

296

297class StringID(StringPattern):

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

299

300

301class ScopeType(StringChoice):

302 choices = [

303 'global',

304 'regional',

305 'local',

306 ]

307

308

309class WaveType(StringChoice):

310 choices = [

311 'full waveform',

312 'bodywave',

313 'P wave',

314 'S wave',

315 'surface wave',

316 ]

317

318

319class NearfieldTermsType(StringChoice):

320 choices = [

321 'complete',

322 'incomplete',

323 'missing',

324 ]

325

326

327class QuantityType(StringChoice):

328 choices = [

329 'displacement',

330 'velocity',

331 'acceleration',

332 'rotation_displacement',

333 'rotation_velocity',

334 'rotation_acceleration',

335 'pressure',

336 'tilt',

337 'pore_pressure',

338 'darcy_velocity',

339 'vertical_tilt']

340

341

342class Reference(Object):

343 id = StringID.T()

344 type = String.T()

345 title = Unicode.T()

346 journal = Unicode.T(optional=True)

347 volume = Unicode.T(optional=True)

348 number = Unicode.T(optional=True)

349 pages = Unicode.T(optional=True)

350 year = String.T()

351 note = Unicode.T(optional=True)

352 issn = String.T(optional=True)

353 doi = String.T(optional=True)

354 url = String.T(optional=True)

355 eprint = String.T(optional=True)

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

357 publisher = Unicode.T(optional=True)

358 keywords = Unicode.T(optional=True)

359 note = Unicode.T(optional=True)

360 abstract = Unicode.T(optional=True)

361

362 @classmethod

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

364

365 from pybtex.database.input import bibtex

366

367 parser = bibtex.Parser()

368

369 if filename is not None:

370 bib_data = parser.parse_file(filename)

371 elif stream is not None:

372 bib_data = parser.parse_stream(stream)

373

374 references = []

375

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

377 d = {}

378 avail = entry.fields.keys()

379 for prop in cls.T.properties:

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

381 continue

382

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

384

385 if 'author' in entry.persons:

386 d['authors'] = []

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

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

389

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

391 references.append(c)

392

393 return references

394

395

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

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

398_cake_pat = cake.PhaseDef.allowed_characters_pattern

399_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic

400

401_ppat = r'(' \

402 r'cake:' + _cake_pat + \

403 r'|iaspei:' + _iaspei_pat + \

404 r'|vel_surface:' + _fpat + \

405 r'|vel:' + _fpat + \

406 r'|stored:' + _spat + \

407 r')'

408

409

410timing_regex_old = re.compile(

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

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

413

414

415timing_regex = re.compile(

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

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

418

419

420class PhaseSelect(StringChoice):

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

422

423

424class InvalidTimingSpecification(ValidationError):

425 pass

426

427

428class TimingAttributeError(Exception):

429 pass

430

431

432class Timing(SObject):

433 '''

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

435

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

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

438 arrival or group of such arrivals.

439

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

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

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

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

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

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

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

447

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

449

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

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

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

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

454 velocity [km/s]

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

456

457 **Examples:**

458

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

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

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

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

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

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

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

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

467 undefined for a given geometry

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

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

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

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

472 arrivals A, B, and C

473 '''

474

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

476

477 if s is not None:

478 offset_is = None

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

480 try:

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

482

483 if s.endswith('S'):

484 offset_is = 'slowness'

485

486 phase_defs = []

487 select = ''

488

489 except ValueError:

490 matched = False

491 m = timing_regex.match(s)

492 if m:

493 if m.group(3):

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

495 elif m.group(19):

496 phase_defs = [m.group(19)]

497 else:

498 phase_defs = []

499

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

501

502 offset = 0.0

503 soff = m.group(27)

504 if soff:

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

506 if soff.endswith('S'):

507 offset_is = 'slowness'

508 elif soff.endswith('%'):

509 offset_is = 'percent'

510

511 matched = True

512

513 else:

514 m = timing_regex_old.match(s)

515 if m:

516 if m.group(3):

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

518 elif m.group(5):

519 phase_defs = [m.group(5)]

520 else:

521 phase_defs = []

522

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

524

525 offset = 0.0

526 if m.group(6):

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

528

529 phase_defs = [

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

531

532 matched = True

533

534 if not matched:

535 raise InvalidTimingSpecification(s)

536

537 kwargs = dict(

538 phase_defs=phase_defs,

539 select=select,

540 offset=offset,

541 offset_is=offset_is)

542

543 SObject.__init__(self, **kwargs)

544

545 def __str__(self):

546 s = []

547 if self.phase_defs:

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

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

550 sphases = '{%s}' % sphases

551

552 if self.select:

553 sphases = self.select + sphases

554

555 s.append(sphases)

556

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

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

559 if self.offset_is == 'slowness':

560 s.append('S')

561 elif self.offset_is == 'percent':

562 s.append('%')

563

564 return ''.join(s)

565

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

567 try:

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

569 phase_offset = get_phase(

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

571 offset = phase_offset(args)

572 else:

573 offset = self.offset

574

575 if self.phase_defs:

576 extra_args = ()

577 if attributes:

578 extra_args = (attributes,)

579

580 phases = [

581 get_phase(phase_def, *extra_args)

582 for phase_def in self.phase_defs]

583

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

585 results = [

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

587 for result in results]

588

589 results = [

590 result for result in results

591 if result[0] is not None]

592

593 if self.select == 'first':

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

595 elif self.select == 'last':

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

597

598 if not results:

599 if attributes is not None:

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

601 else:

602 return None

603

604 else:

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

606 if self.offset_is == 'percent':

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

608 else:

609 t = t + offset

610

611 if attributes is not None:

612 return (t, ) + values

613 else:

614 return t

615

616 else:

617 if attributes is not None:

618 raise TimingAttributeError(

619 'Cannot get phase attributes from offset-only '

620 'definition.')

621 return offset

622

623 except spit.OutOfBounds:

624 raise OutOfBounds(args)

625

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

627 offset = Float.T(default=0.0)

628 offset_is = String.T(optional=True)

629 select = PhaseSelect.T(

630 default='',

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

632 tuple(PhaseSelect.choices)))

633

634

635def mkdefs(s):

636 defs = []

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

638 try:

639 defs.append(float(x))

640 except ValueError:

641 if x.startswith('!'):

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

643 else:

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

645

646 return defs

647

648

649class TPDef(Object):

650 '''

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

652 '''

653

654 id = StringID.T(

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

656 definition = String.T(

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

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

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

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

661

662 @property

663 def phases(self):

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

665 if isinstance(x, cake.PhaseDef)]

666

667 @property

668 def horizontal_velocities(self):

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

670

671

672class OutOfBounds(Exception):

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

674 Exception.__init__(self)

675 self.values = values

676 self.reason = reason

677 self.context = None

678

679 def __str__(self):

680 scontext = ''

681 if self.context:

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

683

684 if self.reason:

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

686

687 if self.values:

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

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

690 else:

691 return 'out of bounds%s' % scontext

692

693

694class MultiLocation(Object):

695 '''

696 Unstructured grid of point coordinates.

697 '''

698

699 lats = Array.T(

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

701 help='Latitudes of targets.')

702 lons = Array.T(

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

704 help='Longitude of targets.')

705 north_shifts = Array.T(

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

707 help='North shifts of targets.')

708 east_shifts = Array.T(

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

710 help='East shifts of targets.')

711 elevation = Array.T(

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

713 help='Elevations of targets.')

714

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

716 self._coords5 = None

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

718

719 @property

720 def coords5(self):

721 if self._coords5 is not None:

722 return self._coords5

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

724 self.elevation]

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

726 if not sizes:

727 raise AttributeError('Empty StaticTarget')

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

729 raise AttributeError('Inconsistent coordinate sizes.')

730

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

732 for idx, p in enumerate(props):

733 if p is not None:

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

735

736 return self._coords5

737

738 @property

739 def ncoords(self):

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

741 return 0

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

743

744 def get_latlon(self):

745 '''

746 Get all coordinates as lat/lon.

747

748 :returns: Coordinates in Latitude, Longitude

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

750 '''

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

752 for i in range(self.ncoords):

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

754 return latlons

755

756 def get_corner_coords(self):

757 '''

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

759

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

761 :rtype: tuple

762 '''

763 latlon = self.get_latlon()

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

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

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

767

768

769class Receiver(Location):

770 codes = Tuple.T(

771 3, String.T(),

772 optional=True,

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

774

775 def pyrocko_station(self):

776 from pyrocko import model

777

778 lat, lon = self.effective_latlon

779

780 return model.Station(

781 network=self.codes[0],

782 station=self.codes[1],

783 location=self.codes[2],

784 lat=lat,

785 lon=lon,

786 depth=self.depth)

787

788

789def g(x, d):

790 if x is None:

791 return d

792 else:

793 return x

794

795

796class UnavailableScheme(Exception):

797 '''

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

799 '''

800 pass

801

802

803class InvalidNComponents(Exception):

804 '''

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

806 '''

807 pass

808

809

810class DiscretizedSource(Object):

811 '''

812 Base class for discretized sources.

813

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

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

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

817 source distributions are represented by subclasses of the

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

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

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

821 explosion/implosion type source distributions. The geometry calculations

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

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

824

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

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

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

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

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

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

831 '''

832

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

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

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

836 lat = Float.T(optional=True)

837 lon = Float.T(optional=True)

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

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

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

841 dl = Float.T(optional=True)

842 dw = Float.T(optional=True)

843 nl = Float.T(optional=True)

844 nw = Float.T(optional=True)

845

846 @classmethod

847 def check_scheme(cls, scheme):

848 '''

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

850

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

852 supported by this discretized source class.

853 '''

854

855 if scheme not in cls.provided_schemes:

856 raise UnavailableScheme(

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

858 (cls.__name__, scheme))

859

860 def __init__(self, **kwargs):

861 Object.__init__(self, **kwargs)

862 self._latlons = None

863

864 def __setattr__(self, name, value):

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

866 'lats', 'lons'):

867 self.__dict__['_latlons'] = None

868

869 Object.__setattr__(self, name, value)

870

871 def get_source_terms(self, scheme):

872 raise NotImplementedError()

873

874 def make_weights(self, receiver, scheme):

875 raise NotImplementedError()

876

877 @property

878 def effective_latlons(self):

879 '''

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

881 '''

882

883 if self._latlons is None:

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

885 if (self.north_shifts is not None and

886 self.east_shifts is not None):

887 self._latlons = orthodrome.ne_to_latlon(

888 self.lats, self.lons,

889 self.north_shifts, self.east_shifts)

890 else:

891 self._latlons = self.lats, self.lons

892 else:

893 lat = g(self.lat, 0.0)

894 lon = g(self.lon, 0.0)

895 self._latlons = orthodrome.ne_to_latlon(

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

897

898 return self._latlons

899

900 @property

901 def effective_north_shifts(self):

902 if self.north_shifts is not None:

903 return self.north_shifts

904 else:

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

906

907 @property

908 def effective_east_shifts(self):

909 if self.east_shifts is not None:

910 return self.east_shifts

911 else:

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

913

914 def same_origin(self, receiver):

915 '''

916 Check if receiver has same reference point.

917 '''

918

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

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

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

922

923 def azibazis_to(self, receiver):

924 '''

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

926 points.

927 '''

928

929 if self.same_origin(receiver):

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

931 receiver.north_shift - self.north_shifts)

932 bazis = azis + 180.

933 else:

934 slats, slons = self.effective_latlons

935 rlat, rlon = receiver.effective_latlon

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

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

938

939 return azis, bazis

940

941 def distances_to(self, receiver):

942 '''

943 Compute distances to receiver for all contained points.

944 '''

945 if self.same_origin(receiver):

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

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

948

949 else:

950 slats, slons = self.effective_latlons

951 rlat, rlon = receiver.effective_latlon

952 return orthodrome.distance_accurate50m_numpy(slats, slons,

953 rlat, rlon)

954

955 def element_coords(self, i):

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

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

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

959 else:

960 lat = self.lat

961 lon = self.lon

962

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

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

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

966

967 else:

968 north_shift = east_shift = 0.0

969

970 return lat, lon, north_shift, east_shift

971

972 def coords5(self):

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

974

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

976 xs[:, 0] = self.lats

977 xs[:, 1] = self.lons

978 else:

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

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

981

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

983 xs[:, 2] = self.north_shifts

984 xs[:, 3] = self.east_shifts

985

986 xs[:, 4] = self.depths

987

988 return xs

989

990 @property

991 def nelements(self):

992 return self.times.size

993

994 @classmethod

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

996 '''

997 Combine several discretized source models.

998

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

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

1001 factors and reference times of the parameterized (undiscretized)

1002 sources match or are accounted for.

1003 '''

1004

1005 first = sources[0]

1006

1007 if not all(type(s) is type(first) for s in sources):

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

1009 'sources of same type.')

1010

1011 latlons = []

1012 for s in sources:

1013 latlons.append(s.effective_latlons)

1014

1015 lats, lons = num.hstack(latlons)

1016

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

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

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

1020 same_ref = num.all(

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

1022 else:

1023 same_ref = False

1024

1025 cat = num.concatenate

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

1027 print(times)

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

1029

1030 if same_ref:

1031 lat = first.lat

1032 lon = first.lon

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

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

1035 lats = None

1036 lons = None

1037 else:

1038 lat = None

1039 lon = None

1040 north_shifts = None

1041 east_shifts = None

1042

1043 return cls(

1044 times=times, lat=lat, lon=lon, lats=lats, lons=lons,

1045 north_shifts=north_shifts, east_shifts=east_shifts,

1046 depths=depths, **kwargs)

1047

1048 def centroid_position(self):

1049 moments = self.moments()

1050 norm = num.sum(moments)

1051 if norm != 0.0:

1052 w = moments / num.sum(moments)

1053 else:

1054 w = num.ones(moments.size)

1055

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

1057 lats, lons = self.effective_latlons

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

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

1060 else:

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

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

1063

1064 cn = num.sum(n*w)

1065 ce = num.sum(e*w)

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

1067

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

1069 lat = clat

1070 lon = clon

1071 north_shift = 0.

1072 east_shift = 0.

1073 else:

1074 lat = g(self.lat, 0.0)

1075 lon = g(self.lon, 0.0)

1076 north_shift = cn

1077 east_shift = ce

1078

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

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

1081 return tuple(float(x) for x in

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

1083

1084

1085class DiscretizedExplosionSource(DiscretizedSource):

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

1087

1088 provided_schemes = (

1089 'elastic2',

1090 'elastic8',

1091 'elastic10',

1092 'rotational8',

1093 )

1094

1095 def get_source_terms(self, scheme):

1096 self.check_scheme(scheme)

1097

1098 if scheme == 'elastic2':

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

1100

1101 elif scheme in ('elastic8', 'elastic10', 'rotational8'):

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

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

1104 return m6s

1105 else:

1106 assert False

1107

1108 def make_weights(self, receiver, scheme):

1109 self.check_scheme(scheme)

1110

1111 azis, bazis = self.azibazis_to(receiver)

1112

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

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

1115

1116 m0s = self.m0s

1117 n = azis.size

1118

1119 cat = num.concatenate

1120 rep = num.repeat

1121

1122 if scheme == 'elastic2':

1123 w_n = cb*m0s

1124 g_n = filledi(0, n)

1125 w_e = sb*m0s

1126 g_e = filledi(0, n)

1127 w_d = m0s

1128 g_d = filledi(1, n)

1129

1130 elif scheme == 'elastic8':

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

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

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

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

1135 w_d = cat((m0s, m0s))

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

1137

1138 elif scheme == 'elastic10':

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

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

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

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

1143 w_d = cat((m0s, m0s, m0s))

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

1145

1146 elif scheme == 'rotational8':

1147 w_n = cat((sb*m0s, sb*m0s, sb*m0s))

1148 g_n = rep((2, 4, 7), n)

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

1150 g_e = rep((2, 4, 7), n)

1151 w_d = num.zeros(0)

1152 g_d = num.zeros(0, dtype=int)

1153 logger.warning('untested code executed, check for sign errors')

1154

1155 else:

1156 assert False

1157

1158 return (

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

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

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

1162 )

1163

1164 def split(self):

1165 from pyrocko.gf.seismosizer import ExplosionSource

1166 sources = []

1167 for i in range(self.nelements):

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

1169 sources.append(ExplosionSource(

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

1171 lat=lat,

1172 lon=lon,

1173 north_shift=north_shift,

1174 east_shift=east_shift,

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

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

1177

1178 return sources

1179

1180 def moments(self):

1181 return self.m0s

1182

1183 def centroid(self):

1184 from pyrocko.gf.seismosizer import ExplosionSource

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

1186 self.centroid_position()

1187

1188 return ExplosionSource(

1189 time=time,

1190 lat=lat,

1191 lon=lon,

1192 north_shift=north_shift,

1193 east_shift=east_shift,

1194 depth=depth,

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

1196

1197 @classmethod

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

1199 '''

1200 Combine several discretized source models.

1201

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

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

1204 factors and reference times of the parameterized (undiscretized)

1205 sources match or are accounted for.

1206 '''

1207

1208 if 'm0s' not in kwargs:

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

1210

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

1212 **kwargs)

1213

1214

1215class DiscretizedSFSource(DiscretizedSource):

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

1217

1218 provided_schemes = (

1219 'elastic5',

1220 )

1221

1222 def get_source_terms(self, scheme):

1223 self.check_scheme(scheme)

1224

1225 return self.forces

1226

1227 def make_weights(self, receiver, scheme):

1228 self.check_scheme(scheme)

1229

1230 azis, bazis = self.azibazis_to(receiver)

1231

1232 sa = num.sin(azis*d2r)

1233 ca = num.cos(azis*d2r)

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

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

1236

1237 forces = self.forces

1238 fn = forces[:, 0]

1239 fe = forces[:, 1]

1240 fd = forces[:, 2]

1241

1242 f0 = fd

1243 f1 = ca * fn + sa * fe

1244 f2 = ca * fe - sa * fn

1245

1246 n = azis.size

1247

1248 if scheme == 'elastic5':

1249 ioff = 0

1250

1251 cat = num.concatenate

1252 rep = num.repeat

1253

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

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

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

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

1258 w_d = cat((f0, f1))

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

1260

1261 return (

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

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

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

1265 )

1266

1267 @classmethod

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

1269 '''

1270 Combine several discretized source models.

1271

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

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

1274 factors and reference times of the parameterized (undiscretized)

1275 sources match or are accounted for.

1276 '''

1277

1278 if 'forces' not in kwargs:

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

1280

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

1282

1283 def moments(self):

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

1285

1286 def centroid(self):

1287 from pyrocko.gf.seismosizer import SFSource

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

1289 self.centroid_position()

1290

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

1292 return SFSource(

1293 time=time,

1294 lat=lat,

1295 lon=lon,

1296 north_shift=north_shift,

1297 east_shift=east_shift,

1298 depth=depth,

1299 fn=fn,

1300 fe=fe,

1301 fd=fd)

1302

1303

1304class DiscretizedMTSource(DiscretizedSource):

1305 m6s = Array.T(

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

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

1308

1309 provided_schemes = (

1310 'elastic8',

1311 'elastic10',

1312 'elastic18',

1313 'rotational8',

1314 )

1315

1316 def get_source_terms(self, scheme):

1317 self.check_scheme(scheme)

1318 return self.m6s

1319

1320 def make_weights(self, receiver, scheme):

1321 self.check_scheme(scheme)

1322

1323 m6s = self.m6s

1324 n = m6s.shape[0]

1325 rep = num.repeat

1326

1327 if scheme == 'elastic18':

1328 w_n = m6s.flatten()

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

1330 w_e = m6s.flatten()

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

1332 w_d = m6s.flatten()

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

1334

1335 else:

1336 azis, bazis = self.azibazis_to(receiver)

1337

1338 sa = num.sin(azis*d2r)

1339 ca = num.cos(azis*d2r)

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

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

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

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

1344

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

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

1347 f2 = m6s[:, 2]

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

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

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

1351

1352 cat = num.concatenate

1353

1354 if scheme == 'elastic8':

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

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

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

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

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

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

1361

1362 elif scheme == 'elastic10':

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

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

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

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

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

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

1369

1370 elif scheme == 'rotational18':

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

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

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

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

1375 w_d = cat((f3, f4))

1376 g_d = rep((5, 6), n)

1377

1378 else:

1379 assert False

1380

1381 return (

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

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

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

1385 )

1386

1387 def split(self):

1388 from pyrocko.gf.seismosizer import MTSource

1389 sources = []

1390 for i in range(self.nelements):

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

1392 sources.append(MTSource(

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

1394 lat=lat,

1395 lon=lon,

1396 north_shift=north_shift,

1397 east_shift=east_shift,

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

1399 m6=self.m6s[i]))

1400

1401 return sources

1402

1403 def moments(self):

1404 moments = num.array(

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

1406 for m6 in self.m6s])

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

1408

1409 def get_moment_rate(self, deltat=None):

1410 moments = self.moments()

1411 times = self.times

1412 times -= times.min()

1413

1414 t_max = times.max()

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

1416 mom_times[mom_times > t_max] = t_max

1417

1418 # Right open histrogram bins

1419 mom, _ = num.histogram(

1420 times,

1421 bins=mom_times,

1422 weights=moments)

1423

1424 deltat = num.diff(mom_times)

1425 mom_rate = mom / deltat

1426

1427 return mom_rate, mom_times[1:]

1428

1429 def centroid(self):

1430 from pyrocko.gf.seismosizer import MTSource

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

1432 self.centroid_position()

1433

1434 return MTSource(

1435 time=time,

1436 lat=lat,

1437 lon=lon,

1438 north_shift=north_shift,

1439 east_shift=east_shift,

1440 depth=depth,

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

1442

1443 @classmethod

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

1445 '''

1446 Combine several discretized source models.

1447

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

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

1450 factors and reference times of the parameterized (undiscretized)

1451 sources match or are accounted for.

1452 '''

1453

1454 if 'm6s' not in kwargs:

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

1456

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

1458

1459

1460class DiscretizedPorePressureSource(DiscretizedSource):

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

1462

1463 provided_schemes = (

1464 'poroelastic10',

1465 )

1466

1467 def get_source_terms(self, scheme):

1468 self.check_scheme(scheme)

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

1470

1471 def make_weights(self, receiver, scheme):

1472 self.check_scheme(scheme)

1473

1474 azis, bazis = self.azibazis_to(receiver)

1475

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

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

1478

1479 pp = self.pp

1480 n = bazis.size

1481

1482 w_un = cb*pp

1483 g_un = filledi(1, n)

1484 w_ue = sb*pp

1485 g_ue = filledi(1, n)

1486 w_ud = pp

1487 g_ud = filledi(0, n)

1488

1489 w_tn = cb*pp

1490 g_tn = filledi(6, n)

1491 w_te = sb*pp

1492 g_te = filledi(6, n)

1493

1494 w_pp = pp

1495 g_pp = filledi(7, n)

1496

1497 w_dvn = cb*pp

1498 g_dvn = filledi(9, n)

1499 w_dve = sb*pp

1500 g_dve = filledi(9, n)

1501 w_dvd = pp

1502 g_dvd = filledi(8, n)

1503

1504 return (

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

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

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

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

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

1510 ('pore_pressure', w_pp, g_pp),

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

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

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

1514 )

1515

1516 def moments(self):

1517 return self.pp

1518

1519 def centroid(self):

1520 from pyrocko.gf.seismosizer import PorePressurePointSource

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

1522 self.centroid_position()

1523

1524 return PorePressurePointSource(

1525 time=time,

1526 lat=lat,

1527 lon=lon,

1528 north_shift=north_shift,

1529 east_shift=east_shift,

1530 depth=depth,

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

1532

1533 @classmethod

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

1535 '''

1536 Combine several discretized source models.

1537

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

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

1540 factors and reference times of the parameterized (undiscretized)

1541 sources match or are accounted for.

1542 '''

1543

1544 if 'pp' not in kwargs:

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

1546

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

1548 **kwargs)

1549

1550

1551class Region(Object):

1552 name = String.T(optional=True)

1553

1554

1555class RectangularRegion(Region):

1556 lat_min = Float.T()

1557 lat_max = Float.T()

1558 lon_min = Float.T()

1559 lon_max = Float.T()

1560

1561

1562class CircularRegion(Region):

1563 lat = Float.T()

1564 lon = Float.T()

1565 radius = Float.T()

1566

1567

1568class Config(Object):

1569 '''

1570 Green's function store meta information.

1571

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

1573 configuration types are:

1574

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

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

1577

1578 * Problem is invariant to horizontal translations and rotations around

1579 vertical axis.

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

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

1582 component)``

1583

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

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

1586

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

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

1589 receiver_depth, component)``

1590

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

1592 constraints but fixed receiver positions

1593

1594 * Cartesian source volume around a reference point

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

1596 source_east_shift, source_north_shift, component)``

1597 '''

1598

1599 id = StringID.T(

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

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

1602 'letter.')

1603

1604 derived_from_id = StringID.T(

1605 optional=True,

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

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

1608

1609 version = String.T(

1610 default='1.0',

1611 optional=True,

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

1613

1614 modelling_code_id = StringID.T(

1615 optional=True,

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

1617

1618 author = Unicode.T(

1619 optional=True,

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

1621

1622 author_email = String.T(

1623 optional=True,

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

1625

1626 created_time = Timestamp.T(

1627 optional=True,

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

1629

1630 regions = List.T(

1631 Region.T(),

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

1633

1634 scope_type = ScopeType.T(

1635 optional=True,

1636 help='Distance range scope of the store '

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

1638

1639 waveform_type = WaveType.T(

1640 optional=True,

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

1642

1643 nearfield_terms = NearfieldTermsType.T(

1644 optional=True,

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

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

1647

1648 description = String.T(

1649 optional=True,

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

1651

1652 references = List.T(

1653 Reference.T(),

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

1655 'related work.')

1656

1657 earthmodel_1d = Earthmodel1D.T(

1658 optional=True,

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

1660

1661 earthmodel_receiver_1d = Earthmodel1D.T(

1662 optional=True,

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

1664

1665 can_interpolate_source = Bool.T(

1666 optional=True,

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

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

1669

1670 can_interpolate_receiver = Bool.T(

1671 optional=True,

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

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

1674

1675 frequency_min = Float.T(

1676 optional=True,

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

1678

1679 frequency_max = Float.T(

1680 optional=True,

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

1682

1683 sample_rate = Float.T(

1684 optional=True,

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

1686

1687 factor = Float.T(

1688 default=1.0,

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

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

1691 optional=True)

1692

1693 component_scheme = ComponentScheme.T(

1694 default='elastic10',

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

1696

1697 stored_quantity = QuantityType.T(

1698 optional=True,

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

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

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

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

1703

1704 tabulated_phases = List.T(

1705 TPDef.T(),

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

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

1708

1709 ncomponents = Int.T(

1710 optional=True,

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

1712

1713 uuid = String.T(

1714 optional=True,

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

1716 'GF store for practical purposes.')

1717

1718 reference = String.T(

1719 optional=True,

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

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

1722

1723 def __init__(self, **kwargs):

1724 self._do_auto_updates = False

1725 Object.__init__(self, **kwargs)

1726 self._index_function = None

1727 self._indices_function = None

1728 self._vicinity_function = None

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

1730 self._do_auto_updates = True

1731 self.update()

1732

1733 def check_ncomponents(self):

1734 ncomponents = component_scheme_to_description[

1735 self.component_scheme].ncomponents

1736

1737 if self.ncomponents is None:

1738 self.ncomponents = ncomponents

1739 elif ncomponents != self.ncomponents:

1740 raise InvalidNComponents(

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

1742 self.ncomponents, self.component_scheme))

1743

1744 def __setattr__(self, name, value):

1745 Object.__setattr__(self, name, value)

1746 try:

1747 self.T.get_property(name)

1748 if self._do_auto_updates:

1749 self.update()

1750

1751 except ValueError:

1752 pass

1753

1754 def update(self):

1755 self.check_ncomponents()

1756 self._update()

1757 self._make_index_functions()

1758

1759 def irecord(self, *args):

1760 return self._index_function(*args)

1761

1762 def irecords(self, *args):

1763 return self._indices_function(*args)

1764

1765 def vicinity(self, *args):

1766 return self._vicinity_function(*args)

1767

1768 def vicinities(self, *args):

1769 return self._vicinities_function(*args)

1770

1771 def grid_interpolation_coefficients(self, *args):

1772 return self._grid_interpolation_coefficients(*args)

1773

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

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

1776

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

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

1779

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

1781 i = 0

1782 arrs = []

1783 ntotal = 1

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

1785 if gdef and len(gdef) > i:

1786 sssn = gdef[i]

1787 else:

1788 sssn = (None,)*4

1789

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

1791 ntotal *= len(arr)

1792

1793 arrs.append(arr)

1794 i += 1

1795

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

1797 return nditer_outer(arrs[:level])

1798

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

1800 nthreads=0):

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

1802

1803 out = []

1804 delays = source.times

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

1806 receiver,

1807 self.component_scheme):

1808

1809 weights *= self.factor

1810

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

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

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

1814

1815 return out

1816

1817 def short_info(self):

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

1819

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

1821 interpolation=None):

1822 '''

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

1824

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

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

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

1828 ``(lat, lon)``

1829 :param interpolation: interpolation method. Choose from

1830 ``('nearest_neighbor', 'multilinear')``

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

1832 point

1833

1834 The default implementation retrieves and interpolates the shear moduli

1835 from the contained 1D velocity profile.

1836 '''

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

1838 parameter='shear_moduli',

1839 interpolation=interpolation)

1840

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

1842 interpolation=None):

1843 '''

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

1845

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

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

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

1849 ``(lat, lon)``

1850 :param interpolation: interpolation method. Choose from

1851 ``('nearest_neighbor', 'multilinear')``

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

1853 point

1854

1855 The default implementation retrieves and interpolates the lambda moduli

1856 from the contained 1D velocity profile.

1857 '''

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

1859 parameter='lambda_moduli',

1860 interpolation=interpolation)

1861

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

1863 interpolation=None):

1864 '''

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

1866

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

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

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

1870 ``(lat, lon)``

1871 :param interpolation: interpolation method. Choose from

1872 ``('nearest_neighbor', 'multilinear')``

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

1874 point

1875

1876 The default implementation retrieves and interpolates the lambda moduli

1877 from the contained 1D velocity profile.

1878 '''

1879 lambda_moduli = self.get_material_property(

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

1881 interpolation=interpolation)

1882 shear_moduli = self.get_material_property(

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

1884 interpolation=interpolation)

1885 return lambda_moduli + (2 / 3) * shear_moduli

1886

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

1888 '''

1889 Get Vs at given points from contained velocity model.

1890

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

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

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

1894 ``(lat, lon)``

1895 :param interpolation: interpolation method. Choose from

1896 ``('nearest_neighbor', 'multilinear')``

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

1898 point

1899

1900 The default implementation retrieves and interpolates Vs

1901 from the contained 1D velocity profile.

1902 '''

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

1904 parameter='vs',

1905 interpolation=interpolation)

1906

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

1908 '''

1909 Get Vp at given points from contained velocity model.

1910

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

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

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

1914 ``(lat, lon)``

1915 :param interpolation: interpolation method. Choose from

1916 ``('nearest_neighbor', 'multilinear')``

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

1918 point

1919

1920 The default implementation retrieves and interpolates Vp

1921 from the contained 1D velocity profile.

1922 '''

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

1924 parameter='vp',

1925 interpolation=interpolation)

1926

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

1928 '''

1929 Get rho at given points from contained velocity model.

1930

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

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

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

1934 ``(lat, lon)``

1935 :param interpolation: interpolation method. Choose from

1936 ``('nearest_neighbor', 'multilinear')``

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

1938 point

1939

1940 The default implementation retrieves and interpolates rho

1941 from the contained 1D velocity profile.

1942 '''

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

1944 parameter='rho',

1945 interpolation=interpolation)

1946

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

1948 interpolation=None):

1949

1950 if interpolation is None:

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

1952 'multilinear', 'nearest_neighbor')

1953

1954 earthmod = self.earthmodel_1d

1955 store_depth_profile = self.get_source_depths()

1956 z_profile = earthmod.profile('z')

1957

1958 if parameter == 'vs':

1959 vs_profile = earthmod.profile('vs')

1960 profile = num.interp(

1961 store_depth_profile, z_profile, vs_profile)

1962

1963 elif parameter == 'vp':

1964 vp_profile = earthmod.profile('vp')

1965 profile = num.interp(

1966 store_depth_profile, z_profile, vp_profile)

1967

1968 elif parameter == 'rho':

1969 rho_profile = earthmod.profile('rho')

1970

1971 profile = num.interp(

1972 store_depth_profile, z_profile, rho_profile)

1973

1974 elif parameter == 'shear_moduli':

1975 vs_profile = earthmod.profile('vs')

1976 rho_profile = earthmod.profile('rho')

1977

1978 store_vs_profile = num.interp(

1979 store_depth_profile, z_profile, vs_profile)

1980 store_rho_profile = num.interp(

1981 store_depth_profile, z_profile, rho_profile)

1982

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

1984

1985 elif parameter == 'lambda_moduli':

1986 vs_profile = earthmod.profile('vs')

1987 vp_profile = earthmod.profile('vp')

1988 rho_profile = earthmod.profile('rho')

1989

1990 store_vs_profile = num.interp(

1991 store_depth_profile, z_profile, vs_profile)

1992 store_vp_profile = num.interp(

1993 store_depth_profile, z_profile, vp_profile)

1994 store_rho_profile = num.interp(

1995 store_depth_profile, z_profile, rho_profile)

1996

1997 profile = store_rho_profile * (

1998 num.power(store_vp_profile, 2) -

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

2000 else:

2001 raise TypeError(

2002 'parameter %s not available' % parameter)

2003

2004 if interpolation == 'multilinear':

2005 kind = 'linear'

2006 elif interpolation == 'nearest_neighbor':

2007 kind = 'nearest'

2008 else:

2009 raise TypeError(

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

2011

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

2013

2014 try:

2015 return interpolator(points[:, 2])

2016 except ValueError:

2017 raise OutOfBounds()

2018

2019 def is_static(self):

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

2021 if self.modelling_code_id.startswith(code):

2022 return True

2023 return False

2024

2025 def is_dynamic(self):

2026 return not self.is_static()

2027

2028 def get_source_depths(self):

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

2030

2031 def get_tabulated_phase(self, phase_id):

2032 '''

2033 Get tabulated phase definition.

2034 '''

2035

2036 for pdef in self.tabulated_phases:

2037 if pdef.id == phase_id:

2038 return pdef

2039

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

2041

2042 def fix_ttt_holes(self, sptree, mode):

2043 raise StoreError(

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

2045 % self.short_type)

2046

2047 @property

2048 def effective_stored_quantity(self):

2049 return self.stored_quantity if self.stored_quantity is not None else \

2050 component_scheme_to_description[self.component_scheme] \

2051 .default_stored_quantity

2052

2053

2054class ConfigTypeA(Config):

2055 '''

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

2057

2058 * Problem is invariant to horizontal translations and rotations around

2059 vertical axis.

2060

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

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

2063 component)``

2064

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

2066 points.

2067 '''

2068

2069 receiver_depth = Float.T(

2070 default=0.0,

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

2072

2073 source_depth_min = Float.T(

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

2075

2076 source_depth_max = Float.T(

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

2078

2079 source_depth_delta = Float.T(

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

2081

2082 distance_min = Float.T(

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

2084

2085 distance_max = Float.T(

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

2087

2088 distance_delta = Float.T(

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

2090

2091 fd_distance_delta = Float.T(

2092 optional=True,

2093 help='Finite differences interval for FD stores [m].')

2094

2095 short_type = 'A'

2096

2097 provided_schemes = [

2098 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10',

2099 'rotational8']

2100

2101 def get_surface_distance(self, args):

2102 return args[1]

2103

2104 def get_distance(self, args):

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

2106

2107 def get_source_depth(self, args):

2108 return args[0]

2109

2110 def get_source_depths(self):

2111 return self.coords[0]

2112

2113 def get_receiver_depth(self, args):

2114 return self.receiver_depth

2115

2116 def _update(self):

2117 self.mins = num.array(

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

2119 self.maxs = num.array(

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

2121 self.deltas = num.array(

2122 [self.source_depth_delta, self.distance_delta],

2123 dtype=float)

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

2125 vicinity_eps).astype(int) + 1

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

2127 self.deltat = 1.0/self.sample_rate

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

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

2130 (mi, ma, n) in

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

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

2133

2134 self.nsource_depths, self.ndistances = self.ns

2135

2136 def _make_index_functions(self):

2137

2138 amin, bmin = self.mins

2139 da, db = self.deltas

2140 na, nb = self.ns

2141

2142 ng = self.ncomponents

2143

2144 def index_function(a, b, ig):

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

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

2147 try:

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

2149 except ValueError:

2150 raise OutOfBounds()

2151

2152 def indices_function(a, b, ig):

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

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

2155 try:

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

2157 except ValueError:

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

2159 try:

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

2161 except ValueError:

2162 raise OutOfBounds()

2163

2164 def grid_interpolation_coefficients(a, b):

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

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

2167 return ias, ibs

2168

2169 def vicinity_function(a, b, ig):

2170 ias, ibs = grid_interpolation_coefficients(a, b)

2171

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

2173 raise OutOfBounds()

2174

2175 indis = []

2176 weights = []

2177 for ia, va in ias:

2178 iia = ia*nb*ng

2179 for ib, vb in ibs:

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

2181 weights.append(va*vb)

2182

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

2184

2185 def vicinities_function(a, b, ig):

2186

2187 xa = (a - amin) / da

2188 xb = (b - bmin) / db

2189

2190 xa_fl = num.floor(xa)

2191 xa_ce = num.ceil(xa)

2192 xb_fl = num.floor(xb)

2193 xb_ce = num.ceil(xb)

2194 va_fl = 1.0 - (xa - xa_fl)

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

2196 vb_fl = 1.0 - (xb - xb_fl)

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

2198

2199 ia_fl = xa_fl.astype(int)

2200 ia_ce = xa_ce.astype(int)

2201 ib_fl = xb_fl.astype(int)

2202 ib_ce = xb_ce.astype(int)

2203

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

2205 raise OutOfBounds()

2206

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

2208 raise OutOfBounds()

2209

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

2211 raise OutOfBounds()

2212

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

2214 raise OutOfBounds()

2215

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

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

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

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

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

2221

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

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

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

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

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

2227

2228 return irecords, weights

2229

2230 self._index_function = index_function

2231 self._indices_function = indices_function

2232 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2233 self._vicinity_function = vicinity_function

2234 self._vicinities_function = vicinities_function

2235

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

2237 nc = icomponents.size

2238 dists = source.distances_to(receiver)

2239 n = dists.size

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

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

2242 icomponents)

2243

2244 def make_indexing_args1(self, source, receiver):

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

2246

2247 @property

2248 def short_extent(self):

2249 return '%g:%g:%g x %g:%g:%g' % (

2250 self.source_depth_min/km,

2251 self.source_depth_max/km,

2252 self.source_depth_delta/km,

2253 self.distance_min/km,

2254 self.distance_max/km,

2255 self.distance_delta/km)

2256

2257 def fix_ttt_holes(self, sptree, mode):

2258 from pyrocko import eikonal_ext, spit

2259

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

2261

2262 delta = self.deltas[-1]

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

2264

2265 nsources, ndistances = self.ns

2266

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

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

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

2270

2271 speeds = self.get_material_property(

2272 0., 0., points,

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

2274 interpolation='multilinear')

2275

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

2277

2278 times = sptree.interpolate_many(nodes)

2279

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

2281 times = times.reshape(speeds.shape)

2282

2283 try:

2284 eikonal_ext.eikonal_solver_fmm_cartesian(

2285 speeds, times, delta)

2286 except eikonal_ext.EikonalExtError as e:

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

2288 logger.debug(

2289 'Got a warning from eikonal solver '

2290 '- may be ok...')

2291 else:

2292 raise

2293

2294 def func(x):

2295 ibs, ics = \

2296 self.grid_interpolation_coefficients(*x)

2297

2298 t = 0

2299 for ib, vb in ibs:

2300 for ic, vc in ics:

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

2302

2303 return t

2304

2305 return spit.SPTree(

2306 f=func,

2307 ftol=sptree.ftol,

2308 xbounds=sptree.xbounds,

2309 xtols=sptree.xtols)

2310

2311

2312class ConfigTypeB(Config):

2313 '''

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

2315

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

2317 receiver depth

2318

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

2320 receiver_distance, component)``

2321 '''

2322

2323 receiver_depth_min = Float.T(

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

2325

2326 receiver_depth_max = Float.T(

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

2328

2329 receiver_depth_delta = Float.T(

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

2331

2332 source_depth_min = Float.T(

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

2334

2335 source_depth_max = Float.T(

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

2337

2338 source_depth_delta = Float.T(

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

2340

2341 distance_min = Float.T(

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

2343

2344 distance_max = Float.T(

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

2346

2347 distance_delta = Float.T(

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

2349

2350 fd_distance_delta = Float.T(

2351 optional=True,

2352 help='Finite differences interval for FD stores [m].')

2353

2354 fd_receiver_depth_delta = Float.T(

2355 optional=True,

2356 help='Finite differences interval for FD stores [m].')

2357

2358 short_type = 'B'

2359

2360 provided_schemes = [

2361 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10',

2362 'rotational8']

2363

2364 def get_distance(self, args):

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

2366

2367 def get_surface_distance(self, args):

2368 return args[2]

2369

2370 def get_source_depth(self, args):

2371 return args[1]

2372

2373 def get_receiver_depth(self, args):

2374 return args[0]

2375

2376 def get_source_depths(self):

2377 return self.coords[1]

2378

2379 def _update(self):

2380 self.mins = num.array([

2381 self.receiver_depth_min,

2382 self.source_depth_min,

2383 self.distance_min],

2384 dtype=float)

2385

2386 self.maxs = num.array([

2387 self.receiver_depth_max,

2388 self.source_depth_max,

2389 self.distance_max],

2390 dtype=float)

2391

2392 self.deltas = num.array([

2393 self.receiver_depth_delta,

2394 self.source_depth_delta,

2395 self.distance_delta],

2396 dtype=float)

2397

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

2399 vicinity_eps).astype(int) + 1

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

2401 self.deltat = 1.0/self.sample_rate

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

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

2404 (mi, ma, n) in

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

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

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

2408

2409 def _make_index_functions(self):

2410

2411 amin, bmin, cmin = self.mins

2412 da, db, dc = self.deltas

2413 na, nb, nc = self.ns

2414 ng = self.ncomponents

2415

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

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

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

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

2420 try:

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

2422 (na, nb, nc, ng))

2423 except ValueError:

2424 raise OutOfBounds()

2425

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

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

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

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

2430 try:

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

2432 (na, nb, nc, ng))

2433 except ValueError:

2434 raise OutOfBounds()

2435

2436 def grid_interpolation_coefficients(a, b, c):

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

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

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

2440 return ias, ibs, ics

2441

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

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

2444

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

2446 raise OutOfBounds()

2447

2448 indis = []

2449 weights = []

2450 for ia, va in ias:

2451 iia = ia*nb*nc*ng

2452 for ib, vb in ibs:

2453 iib = ib*nc*ng

2454 for ic, vc in ics:

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

2456 weights.append(va*vb*vc)

2457

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

2459

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

2461

2462 xa = (a - amin) / da

2463 xb = (b - bmin) / db

2464 xc = (c - cmin) / dc

2465

2466 xa_fl = num.floor(xa)

2467 xa_ce = num.ceil(xa)

2468 xb_fl = num.floor(xb)

2469 xb_ce = num.ceil(xb)

2470 xc_fl = num.floor(xc)

2471 xc_ce = num.ceil(xc)

2472 va_fl = 1.0 - (xa - xa_fl)

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

2474 vb_fl = 1.0 - (xb - xb_fl)

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

2476 vc_fl = 1.0 - (xc - xc_fl)

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

2478

2479 ia_fl = xa_fl.astype(int)

2480 ia_ce = xa_ce.astype(int)

2481 ib_fl = xb_fl.astype(int)

2482 ib_ce = xb_ce.astype(int)

2483 ic_fl = xc_fl.astype(int)

2484 ic_ce = xc_ce.astype(int)

2485

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

2487 raise OutOfBounds()

2488

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

2490 raise OutOfBounds()

2491

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

2493 raise OutOfBounds()

2494

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

2496 raise OutOfBounds()

2497

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

2499 raise OutOfBounds()

2500

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

2502 raise OutOfBounds()

2503

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

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

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

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

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

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

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

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

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

2513

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

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

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

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

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

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

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

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

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

2523

2524 return irecords, weights

2525

2526 self._index_function = index_function

2527 self._indices_function = indices_function

2528 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2529 self._vicinity_function = vicinity_function

2530 self._vicinities_function = vicinities_function

2531

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

2533 nc = icomponents.size

2534 dists = source.distances_to(receiver)

2535 n = dists.size

2536 receiver_depths = num.empty(nc)

2537 receiver_depths.fill(receiver.depth)

2538 return (receiver_depths,

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

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

2541 icomponents)

2542

2543 def make_indexing_args1(self, source, receiver):

2544 return (receiver.depth,

2545 source.depth,

2546 source.distance_to(receiver))

2547

2548 @property

2549 def short_extent(self):

2550 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2551 self.receiver_depth_min/km,

2552 self.receiver_depth_max/km,

2553 self.receiver_depth_delta/km,

2554 self.source_depth_min/km,

2555 self.source_depth_max/km,

2556 self.source_depth_delta/km,

2557 self.distance_min/km,

2558 self.distance_max/km,

2559 self.distance_delta/km)

2560

2561 def fix_ttt_holes(self, sptree, mode):

2562 from pyrocko import eikonal_ext, spit

2563

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

2565

2566 delta = self.deltas[-1]

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

2568

2569 nreceivers, nsources, ndistances = self.ns

2570

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

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

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

2574

2575 speeds = self.get_material_property(

2576 0., 0., points,

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

2578 interpolation='multilinear')

2579

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

2581

2582 receiver_times = []

2583 for ireceiver in range(nreceivers):

2584 nodes = num.hstack([

2585 num.full(

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

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

2588 nodes_sr])

2589

2590 times = sptree.interpolate_many(nodes)

2591

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

2593

2594 times = times.reshape(speeds.shape)

2595

2596 try:

2597 eikonal_ext.eikonal_solver_fmm_cartesian(

2598 speeds, times, delta)

2599 except eikonal_ext.EikonalExtError as e:

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

2601 logger.debug(

2602 'Got a warning from eikonal solver '

2603 '- may be ok...')

2604 else:

2605 raise

2606

2607 receiver_times.append(times)

2608

2609 def func(x):

2610 ias, ibs, ics = \

2611 self.grid_interpolation_coefficients(*x)

2612

2613 t = 0

2614 for ia, va in ias:

2615 times = receiver_times[ia]

2616 for ib, vb in ibs:

2617 for ic, vc in ics:

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

2619

2620 return t

2621

2622 return spit.SPTree(

2623 f=func,

2624 ftol=sptree.ftol,

2625 xbounds=sptree.xbounds,

2626 xtols=sptree.xtols)

2627

2628

2629class ConfigTypeC(Config):

2630 '''

2631 No symmetrical constraints, one fixed receiver position.

2632

2633 * Cartesian 3D source volume around a reference point

2634

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

2636 source_east_shift, source_north_shift, component)``

2637 '''

2638

2639 receiver = Receiver.T(

2640 help='Receiver location')

2641

2642 source_origin = Location.T(

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

2644

2645 source_depth_min = Float.T(

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

2647

2648 source_depth_max = Float.T(

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

2650

2651 source_depth_delta = Float.T(

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

2653

2654 source_east_shift_min = Float.T(

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

2656

2657 source_east_shift_max = Float.T(

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

2659

2660 source_east_shift_delta = Float.T(

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

2662

2663 source_north_shift_min = Float.T(

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

2665

2666 source_north_shift_max = Float.T(

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

2668

2669 source_north_shift_delta = Float.T(

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

2671

2672 short_type = 'C'

2673

2674 provided_schemes = ['elastic18']

2675

2676 def get_surface_distance(self, args):

2677 _, source_east_shift, source_north_shift, _ = args

2678 sorig = self.source_origin

2679 sloc = Location(

2680 lat=sorig.lat,

2681 lon=sorig.lon,

2682 north_shift=sorig.north_shift + source_north_shift,

2683 east_shift=sorig.east_shift + source_east_shift)

2684

2685 return self.receiver.distance_to(sloc)

2686

2687 def get_distance(self, args):

2688 sdepth, source_east_shift, source_north_shift, _ = args

2689 sorig = self.source_origin

2690 sloc = Location(

2691 lat=sorig.lat,

2692 lon=sorig.lon,

2693 north_shift=sorig.north_shift + source_north_shift,

2694 east_shift=sorig.east_shift + source_east_shift)

2695

2696 return self.receiver.distance_3d_to(sloc)

2697

2698 def get_source_depth(self, args):

2699 return args[0]

2700

2701 def get_receiver_depth(self, args):

2702 return self.receiver.depth

2703

2704 def get_source_depths(self):

2705 return self.coords[0]

2706

2707 def _update(self):

2708 self.mins = num.array([

2709 self.source_depth_min,

2710 self.source_east_shift_min,

2711 self.source_north_shift_min],

2712 dtype=float)

2713

2714 self.maxs = num.array([

2715 self.source_depth_max,

2716 self.source_east_shift_max,

2717 self.source_north_shift_max],

2718 dtype=float)

2719

2720 self.deltas = num.array([

2721 self.source_depth_delta,

2722 self.source_east_shift_delta,

2723 self.source_north_shift_delta],

2724 dtype=float)

2725

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

2727 vicinity_eps).astype(int) + 1

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

2729 self.deltat = 1.0/self.sample_rate

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

2731

2732 self.coords = tuple(

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

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

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

2736

2737 def _make_index_functions(self):

2738

2739 amin, bmin, cmin = self.mins

2740 da, db, dc = self.deltas

2741 na, nb, nc = self.ns

2742 ng = self.ncomponents

2743

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

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

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

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

2748 try:

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

2750 (na, nb, nc, ng))

2751 except ValueError:

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

2753

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

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

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

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

2758

2759 try:

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

2761 (na, nb, nc, ng))

2762 except ValueError:

2763 raise OutOfBounds()

2764

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

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

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

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

2769

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

2771 raise OutOfBounds()

2772

2773 indis = []

2774 weights = []

2775 for ia, va in ias:

2776 iia = ia*nb*nc*ng

2777 for ib, vb in ibs:

2778 iib = ib*nc*ng

2779 for ic, vc in ics:

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

2781 weights.append(va*vb*vc)

2782

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

2784

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

2786

2787 xa = (a-amin) / da

2788 xb = (b-bmin) / db

2789 xc = (c-cmin) / dc

2790

2791 xa_fl = num.floor(xa)

2792 xa_ce = num.ceil(xa)

2793 xb_fl = num.floor(xb)

2794 xb_ce = num.ceil(xb)

2795 xc_fl = num.floor(xc)

2796 xc_ce = num.ceil(xc)

2797 va_fl = 1.0 - (xa - xa_fl)

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

2799 vb_fl = 1.0 - (xb - xb_fl)

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

2801 vc_fl = 1.0 - (xc - xc_fl)

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

2803

2804 ia_fl = xa_fl.astype(int)

2805 ia_ce = xa_ce.astype(int)

2806 ib_fl = xb_fl.astype(int)

2807 ib_ce = xb_ce.astype(int)

2808 ic_fl = xc_fl.astype(int)

2809 ic_ce = xc_ce.astype(int)

2810

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

2812 raise OutOfBounds()

2813

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

2815 raise OutOfBounds()

2816

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

2818 raise OutOfBounds()

2819

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

2821 raise OutOfBounds()

2822

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

2824 raise OutOfBounds()

2825

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

2827 raise OutOfBounds()

2828

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

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

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

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

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

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

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

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

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

2838

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

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

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

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

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

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

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

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

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

2848

2849 return irecords, weights

2850

2851 self._index_function = index_function

2852 self._indices_function = indices_function

2853 self._vicinity_function = vicinity_function

2854 self._vicinities_function = vicinities_function

2855

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

2857 nc = icomponents.size

2858

2859 dists = source.distances_to(self.source_origin)

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

2861

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

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

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

2865

2866 n = dists.size

2867

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

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

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

2871 icomponents)

2872

2873 def make_indexing_args1(self, source, receiver):

2874 dist = source.distance_to(self.source_origin)

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

2876

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

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

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

2880

2881 return (source_depth,

2882 source_east_shift,

2883 source_north_shift)

2884

2885 @property

2886 def short_extent(self):

2887 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2888 self.source_depth_min/km,

2889 self.source_depth_max/km,

2890 self.source_depth_delta/km,

2891 self.source_east_shift_min/km,

2892 self.source_east_shift_max/km,

2893 self.source_east_shift_delta/km,

2894 self.source_north_shift_min/km,

2895 self.source_north_shift_max/km,

2896 self.source_north_shift_delta/km)

2897

2898

2899class Weighting(Object):

2900 factor = Float.T(default=1.0)

2901

2902

2903class Taper(Object):

2904 tmin = Timing.T()

2905 tmax = Timing.T()

2906 tfade = Float.T(default=0.0)

2907 shape = StringChoice.T(

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

2909 default='cos',

2910 optional=True)

2911

2912

2913class SimplePattern(SObject):

2914

2915 _pool = {}

2916

2917 def __init__(self, pattern):

2918 self._pattern = pattern

2919 SObject.__init__(self)

2920

2921 def __str__(self):

2922 return self._pattern

2923

2924 @property

2925 def regex(self):

2926 pool = SimplePattern._pool

2927 if self.pattern not in pool:

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

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

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

2931

2932 return pool[self.pattern]

2933

2934 def match(self, s):

2935 return self.regex.match(s)

2936

2937

2938class WaveformType(StringChoice):

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

2940 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',

2941 'envelope_dis', 'envelope_vel', 'envelope_acc']

2942

2943

2944class ChannelSelection(Object):

2945 pattern = SimplePattern.T(optional=True)

2946 min_sample_rate = Float.T(optional=True)

2947 max_sample_rate = Float.T(optional=True)

2948

2949

2950class StationSelection(Object):

2951 includes = SimplePattern.T()

2952 excludes = SimplePattern.T()

2953 distance_min = Float.T(optional=True)

2954 distance_max = Float.T(optional=True)

2955 azimuth_min = Float.T(optional=True)

2956 azimuth_max = Float.T(optional=True)

2957

2958

2959class WaveformSelection(Object):

2960 channel_selection = ChannelSelection.T(optional=True)

2961 station_selection = StationSelection.T(optional=True)

2962 taper = Taper.T()

2963 # filter = FrequencyResponse.T()

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

2965 weighting = Weighting.T(optional=True)

2966 sample_rate = Float.T(optional=True)

2967 gf_store_id = StringID.T(optional=True)

2968

2969

2970def indi12(x, n):

2971 '''

2972 Get linear interpolation index and weight.

2973 '''

2974

2975 r = round(x)

2976 if abs(r - x) < vicinity_eps:

2977 i = int(r)

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

2979 raise OutOfBounds()

2980

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

2982 else:

2983 f = math.floor(x)

2984 i = int(f)

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

2986 raise OutOfBounds()

2987

2988 v = x-f

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

2990

2991

2992def float_or_none(s):

2993 units = {

2994 'k': 1e3,

2995 'M': 1e6,

2996 }

2997

2998 factor = 1.0

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

3000 factor = units[s[-1]]

3001 s = s[:-1]

3002 if not s:

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

3004

3005 if s:

3006 return float(s) * factor

3007 else:

3008 return None

3009

3010

3011class GridSpecError(Exception):

3012 def __init__(self, s):

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

3014

3015

3016def parse_grid_spec(spec):

3017 try:

3018 result = []

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

3020 t = dspec.split('@')

3021 num = start = stop = step = None

3022 if len(t) == 2:

3023 num = int(t[1])

3024 if num <= 0:

3025 raise GridSpecError(spec)

3026

3027 elif len(t) > 2:

3028 raise GridSpecError(spec)

3029

3030 s = t[0]

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

3032 if len(v) == 1:

3033 start = stop = v[0]

3034 if len(v) >= 2:

3035 start, stop = v[0:2]

3036 if len(v) == 3:

3037 step = v[2]

3038

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

3040 raise GridSpecError(spec)

3041

3042 if step == 0.0:

3043 raise GridSpecError(spec)

3044

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

3046

3047 except ValueError:

3048 raise GridSpecError(spec)

3049

3050 return result

3051

3052

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

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

3055 if start is None:

3056 start = ma if swap else mi

3057 if stop is None:

3058 stop = mi if swap else ma

3059 if step is None:

3060 step = -inc if ma < mi else inc

3061 if num is None:

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

3063 raise GridSpecError()

3064

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

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

3067 if abs(stop-stop2) > eps:

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

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

3070 else:

3071 stop = stop2

3072

3073 if start == stop:

3074 num = 1

3075

3076 return start, stop, num

3077

3078

3079def nditer_outer(x):

3080 return num.nditer(

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

3082

3083

3084def nodes(xs):

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

3086 nnodes = num.prod(ns)

3087 ndim = len(xs)

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

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

3090 x = xs[idim]

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

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

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

3094

3095 return nodes

3096

3097

3098def filledi(x, n):

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

3100 a.fill(x)

3101 return a

3102

3103

3104config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]

3105

3106discretized_source_classes = [

3107 DiscretizedExplosionSource,

3108 DiscretizedSFSource,

3109 DiscretizedMTSource,

3110 DiscretizedPorePressureSource]

3111

3112

3113__all__ = '''

3114Earthmodel1D

3115StringID

3116ScopeType

3117WaveformType

3118QuantityType

3119NearfieldTermsType

3120Reference

3121Region

3122CircularRegion

3123RectangularRegion

3124PhaseSelect

3125InvalidTimingSpecification

3126Timing

3127TPDef

3128OutOfBounds

3129Location

3130Receiver

3131'''.split() + [

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

3133ComponentScheme

3134component_scheme_to_description

3135component_schemes

3136Config

3137GridSpecError

3138Weighting

3139Taper

3140SimplePattern

3141WaveformType

3142ChannelSelection

3143StationSelection

3144WaveformSelection

3145nditer_outer

3146dump

3147load

3148discretized_source_classes

3149config_type_classes

3150UnavailableScheme

3151InterpolationMethod

3152SeismosizerTrace

3153SeismosizerResult

3154Result

3155StaticResult

3156TimingAttributeError

3157'''.split()