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

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 depths = cat([s.depths for s in sources])

1028

1029 if same_ref:

1030 lat = first.lat

1031 lon = first.lon

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

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

1034 lats = None

1035 lons = None

1036 else:

1037 lat = None

1038 lon = None

1039 north_shifts = None

1040 east_shifts = None

1041

1042 return cls(

1043 times=times, lat=lat, lon=lon, lats=lats, lons=lons,

1044 north_shifts=north_shifts, east_shifts=east_shifts,

1045 depths=depths, **kwargs)

1046

1047 def centroid_position(self):

1048 moments = self.moments()

1049 norm = num.sum(moments)

1050 if norm != 0.0:

1051 w = moments / num.sum(moments)

1052 else:

1053 w = num.ones(moments.size) / moments.size

1054

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

1056 lats, lons = self.effective_latlons

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

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

1059 else:

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

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

1062

1063 cn = num.sum(n*w)

1064 ce = num.sum(e*w)

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

1066

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

1068 lat = clat

1069 lon = clon

1070 north_shift = 0.

1071 east_shift = 0.

1072 else:

1073 lat = g(self.lat, 0.0)

1074 lon = g(self.lon, 0.0)

1075 north_shift = cn

1076 east_shift = ce

1077

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

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

1080 return tuple(float(x) for x in

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

1082

1083

1084class DiscretizedExplosionSource(DiscretizedSource):

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

1086

1087 provided_schemes = (

1088 'elastic2',

1089 'elastic8',

1090 'elastic10',

1091 'rotational8',

1092 )

1093

1094 def get_source_terms(self, scheme):

1095 self.check_scheme(scheme)

1096

1097 if scheme == 'elastic2':

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

1099

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

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

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

1103 return m6s

1104 else:

1105 assert False

1106

1107 def make_weights(self, receiver, scheme):

1108 self.check_scheme(scheme)

1109

1110 azis, bazis = self.azibazis_to(receiver)

1111

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

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

1114

1115 m0s = self.m0s

1116 n = azis.size

1117

1118 cat = num.concatenate

1119 rep = num.repeat

1120

1121 if scheme == 'elastic2':

1122 w_n = cb*m0s

1123 g_n = filledi(0, n)

1124 w_e = sb*m0s

1125 g_e = filledi(0, n)

1126 w_d = m0s

1127 g_d = filledi(1, n)

1128

1129 elif scheme == 'elastic8':

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

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

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

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

1134 w_d = cat((m0s, m0s))

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

1136

1137 elif scheme == 'elastic10':

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

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

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

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

1142 w_d = cat((m0s, m0s, m0s))

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

1144

1145 elif scheme == 'rotational8':

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

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

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

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

1150 w_d = num.zeros(0)

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

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

1153

1154 else:

1155 assert False

1156

1157 return (

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

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

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

1161 )

1162

1163 def split(self):

1164 from pyrocko.gf.seismosizer import ExplosionSource

1165 sources = []

1166 for i in range(self.nelements):

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

1168 sources.append(ExplosionSource(

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

1170 lat=lat,

1171 lon=lon,

1172 north_shift=north_shift,

1173 east_shift=east_shift,

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

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

1176

1177 return sources

1178

1179 def moments(self):

1180 return self.m0s

1181

1182 def centroid(self):

1183 from pyrocko.gf.seismosizer import ExplosionSource

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

1185 self.centroid_position()

1186

1187 return ExplosionSource(

1188 time=time,

1189 lat=lat,

1190 lon=lon,

1191 north_shift=north_shift,

1192 east_shift=east_shift,

1193 depth=depth,

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

1195

1196 @classmethod

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

1198 '''

1199 Combine several discretized source models.

1200

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

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

1203 factors and reference times of the parameterized (undiscretized)

1204 sources match or are accounted for.

1205 '''

1206

1207 if 'm0s' not in kwargs:

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

1209

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

1211 **kwargs)

1212

1213

1214class DiscretizedSFSource(DiscretizedSource):

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

1216

1217 provided_schemes = (

1218 'elastic5',

1219 )

1220

1221 def get_source_terms(self, scheme):

1222 self.check_scheme(scheme)

1223

1224 return self.forces

1225

1226 def make_weights(self, receiver, scheme):

1227 self.check_scheme(scheme)

1228

1229 azis, bazis = self.azibazis_to(receiver)

1230

1231 sa = num.sin(azis*d2r)

1232 ca = num.cos(azis*d2r)

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

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

1235

1236 forces = self.forces

1237 fn = forces[:, 0]

1238 fe = forces[:, 1]

1239 fd = forces[:, 2]

1240

1241 f0 = fd

1242 f1 = ca * fn + sa * fe

1243 f2 = ca * fe - sa * fn

1244

1245 n = azis.size

1246

1247 if scheme == 'elastic5':

1248 ioff = 0

1249

1250 cat = num.concatenate

1251 rep = num.repeat

1252

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

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

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

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

1257 w_d = cat((f0, f1))

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

1259

1260 return (

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

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

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

1264 )

1265

1266 @classmethod

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

1268 '''

1269 Combine several discretized source models.

1270

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

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

1273 factors and reference times of the parameterized (undiscretized)

1274 sources match or are accounted for.

1275 '''

1276

1277 if 'forces' not in kwargs:

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

1279

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

1281

1282 def moments(self):

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

1284

1285 def centroid(self):

1286 from pyrocko.gf.seismosizer import SFSource

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

1288 self.centroid_position()

1289

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

1291 return SFSource(

1292 time=time,

1293 lat=lat,

1294 lon=lon,

1295 north_shift=north_shift,

1296 east_shift=east_shift,

1297 depth=depth,

1298 fn=fn,

1299 fe=fe,

1300 fd=fd)

1301

1302

1303class DiscretizedMTSource(DiscretizedSource):

1304 m6s = Array.T(

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

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

1307

1308 provided_schemes = (

1309 'elastic8',

1310 'elastic10',

1311 'elastic18',

1312 'rotational8',

1313 )

1314

1315 def get_source_terms(self, scheme):

1316 self.check_scheme(scheme)

1317 return self.m6s

1318

1319 def make_weights(self, receiver, scheme):

1320 self.check_scheme(scheme)

1321

1322 m6s = self.m6s

1323 n = m6s.shape[0]

1324 rep = num.repeat

1325

1326 if scheme == 'elastic18':

1327 w_n = m6s.flatten()

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

1329 w_e = m6s.flatten()

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

1331 w_d = m6s.flatten()

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

1333

1334 else:

1335 azis, bazis = self.azibazis_to(receiver)

1336

1337 sa = num.sin(azis*d2r)

1338 ca = num.cos(azis*d2r)

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

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

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

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

1343

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

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

1346 f2 = m6s[:, 2]

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

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

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

1350

1351 cat = num.concatenate

1352

1353 if scheme == 'elastic8':

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

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

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

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

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

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

1360

1361 elif scheme == 'elastic10':

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

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

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

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

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

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

1368

1369 elif scheme == 'rotational18':

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

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

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

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

1374 w_d = cat((f3, f4))

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

1376

1377 else:

1378 assert False

1379

1380 return (

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

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

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

1384 )

1385

1386 def split(self):

1387 from pyrocko.gf.seismosizer import MTSource

1388 sources = []

1389 for i in range(self.nelements):

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

1391 sources.append(MTSource(

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

1393 lat=lat,

1394 lon=lon,

1395 north_shift=north_shift,

1396 east_shift=east_shift,

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

1398 m6=self.m6s[i]))

1399

1400 return sources

1401

1402 def moments(self):

1403 moments = num.array(

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

1405 for m6 in self.m6s])

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

1407

1408 def get_moment_rate(self, deltat=None):

1409 moments = self.moments()

1410 times = self.times

1411 times -= times.min()

1412

1413 t_max = times.max()

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

1415 mom_times[mom_times > t_max] = t_max

1416

1417 # Right open histrogram bins

1418 mom, _ = num.histogram(

1419 times,

1420 bins=mom_times,

1421 weights=moments)

1422

1423 deltat = num.diff(mom_times)

1424 mom_rate = mom / deltat

1425

1426 return mom_rate, mom_times[1:]

1427

1428 def centroid(self):

1429 from pyrocko.gf.seismosizer import MTSource

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

1431 self.centroid_position()

1432

1433 return MTSource(

1434 time=time,

1435 lat=lat,

1436 lon=lon,

1437 north_shift=north_shift,

1438 east_shift=east_shift,

1439 depth=depth,

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

1441

1442 @classmethod

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

1444 '''

1445 Combine several discretized source models.

1446

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

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

1449 factors and reference times of the parameterized (undiscretized)

1450 sources match or are accounted for.

1451 '''

1452

1453 if 'm6s' not in kwargs:

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

1455

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

1457

1458

1459class DiscretizedPorePressureSource(DiscretizedSource):

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

1461

1462 provided_schemes = (

1463 'poroelastic10',

1464 )

1465

1466 def get_source_terms(self, scheme):

1467 self.check_scheme(scheme)

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

1469

1470 def make_weights(self, receiver, scheme):

1471 self.check_scheme(scheme)

1472

1473 azis, bazis = self.azibazis_to(receiver)

1474

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

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

1477

1478 pp = self.pp

1479 n = bazis.size

1480

1481 w_un = cb*pp

1482 g_un = filledi(1, n)

1483 w_ue = sb*pp

1484 g_ue = filledi(1, n)

1485 w_ud = pp

1486 g_ud = filledi(0, n)

1487

1488 w_tn = cb*pp

1489 g_tn = filledi(6, n)

1490 w_te = sb*pp

1491 g_te = filledi(6, n)

1492

1493 w_pp = pp

1494 g_pp = filledi(7, n)

1495

1496 w_dvn = cb*pp

1497 g_dvn = filledi(9, n)

1498 w_dve = sb*pp

1499 g_dve = filledi(9, n)

1500 w_dvd = pp

1501 g_dvd = filledi(8, n)

1502

1503 return (

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

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

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

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

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

1509 ('pore_pressure', w_pp, g_pp),

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

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

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

1513 )

1514

1515 def moments(self):

1516 return self.pp

1517

1518 def centroid(self):

1519 from pyrocko.gf.seismosizer import PorePressurePointSource

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

1521 self.centroid_position()

1522

1523 return PorePressurePointSource(

1524 time=time,

1525 lat=lat,

1526 lon=lon,

1527 north_shift=north_shift,

1528 east_shift=east_shift,

1529 depth=depth,

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

1531

1532 @classmethod

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

1534 '''

1535 Combine several discretized source models.

1536

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

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

1539 factors and reference times of the parameterized (undiscretized)

1540 sources match or are accounted for.

1541 '''

1542

1543 if 'pp' not in kwargs:

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

1545

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

1547 **kwargs)

1548

1549

1550class Region(Object):

1551 name = String.T(optional=True)

1552

1553

1554class RectangularRegion(Region):

1555 lat_min = Float.T()

1556 lat_max = Float.T()

1557 lon_min = Float.T()

1558 lon_max = Float.T()

1559

1560

1561class CircularRegion(Region):

1562 lat = Float.T()

1563 lon = Float.T()

1564 radius = Float.T()

1565

1566

1567class Config(Object):

1568 '''

1569 Green's function store meta information.

1570

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

1572 configuration types are:

1573

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

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

1576

1577 * Problem is invariant to horizontal translations and rotations around

1578 vertical axis.

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

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

1581 component)``

1582

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

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

1585

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

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

1588 receiver_depth, component)``

1589

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

1591 constraints but fixed receiver positions

1592

1593 * Cartesian source volume around a reference point

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

1595 source_east_shift, source_north_shift, component)``

1596 '''

1597

1598 id = StringID.T(

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

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

1601 'letter.')

1602

1603 derived_from_id = StringID.T(

1604 optional=True,

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

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

1607

1608 version = String.T(

1609 default='1.0',

1610 optional=True,

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

1612

1613 modelling_code_id = StringID.T(

1614 optional=True,

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

1616

1617 author = Unicode.T(

1618 optional=True,

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

1620

1621 author_email = String.T(

1622 optional=True,

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

1624

1625 created_time = Timestamp.T(

1626 optional=True,

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

1628

1629 regions = List.T(

1630 Region.T(),

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

1632

1633 scope_type = ScopeType.T(

1634 optional=True,

1635 help='Distance range scope of the store '

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

1637

1638 waveform_type = WaveType.T(

1639 optional=True,

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

1641

1642 nearfield_terms = NearfieldTermsType.T(

1643 optional=True,

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

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

1646

1647 description = String.T(

1648 optional=True,

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

1650

1651 references = List.T(

1652 Reference.T(),

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

1654 'related work.')

1655

1656 earthmodel_1d = Earthmodel1D.T(

1657 optional=True,

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

1659

1660 earthmodel_receiver_1d = Earthmodel1D.T(

1661 optional=True,

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

1663

1664 can_interpolate_source = Bool.T(

1665 optional=True,

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

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

1668

1669 can_interpolate_receiver = Bool.T(

1670 optional=True,

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

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

1673

1674 frequency_min = Float.T(

1675 optional=True,

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

1677

1678 frequency_max = Float.T(

1679 optional=True,

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

1681

1682 sample_rate = Float.T(

1683 optional=True,

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

1685

1686 factor = Float.T(

1687 default=1.0,

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

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

1690 optional=True)

1691

1692 component_scheme = ComponentScheme.T(

1693 default='elastic10',

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

1695

1696 stored_quantity = QuantityType.T(

1697 optional=True,

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

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

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

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

1702

1703 tabulated_phases = List.T(

1704 TPDef.T(),

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

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

1707

1708 ncomponents = Int.T(

1709 optional=True,

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

1711

1712 uuid = String.T(

1713 optional=True,

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

1715 'GF store for practical purposes.')

1716

1717 reference = String.T(

1718 optional=True,

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

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

1721

1722 def __init__(self, **kwargs):

1723 self._do_auto_updates = False

1724 Object.__init__(self, **kwargs)

1725 self._index_function = None

1726 self._indices_function = None

1727 self._vicinity_function = None

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

1729 self._do_auto_updates = True

1730 self.update()

1731

1732 def check_ncomponents(self):

1733 ncomponents = component_scheme_to_description[

1734 self.component_scheme].ncomponents

1735

1736 if self.ncomponents is None:

1737 self.ncomponents = ncomponents

1738 elif ncomponents != self.ncomponents:

1739 raise InvalidNComponents(

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

1741 self.ncomponents, self.component_scheme))

1742

1743 def __setattr__(self, name, value):

1744 Object.__setattr__(self, name, value)

1745 try:

1746 self.T.get_property(name)

1747 if self._do_auto_updates:

1748 self.update()

1749

1750 except ValueError:

1751 pass

1752

1753 def update(self):

1754 self.check_ncomponents()

1755 self._update()

1756 self._make_index_functions()

1757

1758 def irecord(self, *args):

1759 return self._index_function(*args)

1760

1761 def irecords(self, *args):

1762 return self._indices_function(*args)

1763

1764 def vicinity(self, *args):

1765 return self._vicinity_function(*args)

1766

1767 def vicinities(self, *args):

1768 return self._vicinities_function(*args)

1769

1770 def grid_interpolation_coefficients(self, *args):

1771 return self._grid_interpolation_coefficients(*args)

1772

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

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

1775

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

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

1778

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

1780 i = 0

1781 arrs = []

1782 ntotal = 1

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

1784 if gdef and len(gdef) > i:

1785 sssn = gdef[i]

1786 else:

1787 sssn = (None,)*4

1788

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

1790 ntotal *= len(arr)

1791

1792 arrs.append(arr)

1793 i += 1

1794

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

1796 return nditer_outer(arrs[:level])

1797

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

1799 nthreads=0):

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

1801

1802 out = []

1803 delays = source.times

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

1805 receiver,

1806 self.component_scheme):

1807

1808 weights *= self.factor

1809

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

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

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

1813

1814 return out

1815

1816 def short_info(self):

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

1818

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

1820 interpolation=None):

1821 '''

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

1823

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

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

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

1827 ``(lat, lon)``

1828 :param interpolation: interpolation method. Choose from

1829 ``('nearest_neighbor', 'multilinear')``

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

1831 point

1832

1833 The default implementation retrieves and interpolates the shear moduli

1834 from the contained 1D velocity profile.

1835 '''

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

1837 parameter='shear_moduli',

1838 interpolation=interpolation)

1839

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

1841 interpolation=None):

1842 '''

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

1844

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

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

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

1848 ``(lat, lon)``

1849 :param interpolation: interpolation method. Choose from

1850 ``('nearest_neighbor', 'multilinear')``

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

1852 point

1853

1854 The default implementation retrieves and interpolates the lambda moduli

1855 from the contained 1D velocity profile.

1856 '''

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

1858 parameter='lambda_moduli',

1859 interpolation=interpolation)

1860

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

1862 interpolation=None):

1863 '''

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

1865

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

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

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

1869 ``(lat, lon)``

1870 :param interpolation: interpolation method. Choose from

1871 ``('nearest_neighbor', 'multilinear')``

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

1873 point

1874

1875 The default implementation retrieves and interpolates the lambda moduli

1876 from the contained 1D velocity profile.

1877 '''

1878 lambda_moduli = self.get_material_property(

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

1880 interpolation=interpolation)

1881 shear_moduli = self.get_material_property(

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

1883 interpolation=interpolation)

1884 return lambda_moduli + (2 / 3) * shear_moduli

1885

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

1887 '''

1888 Get Vs at given points from contained velocity model.

1889

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

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

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

1893 ``(lat, lon)``

1894 :param interpolation: interpolation method. Choose from

1895 ``('nearest_neighbor', 'multilinear')``

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

1897 point

1898

1899 The default implementation retrieves and interpolates Vs

1900 from the contained 1D velocity profile.

1901 '''

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

1903 parameter='vs',

1904 interpolation=interpolation)

1905

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

1907 '''

1908 Get Vp at given points from contained velocity model.

1909

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

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

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

1913 ``(lat, lon)``

1914 :param interpolation: interpolation method. Choose from

1915 ``('nearest_neighbor', 'multilinear')``

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

1917 point

1918

1919 The default implementation retrieves and interpolates Vp

1920 from the contained 1D velocity profile.

1921 '''

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

1923 parameter='vp',

1924 interpolation=interpolation)

1925

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

1927 '''

1928 Get rho at given points from contained velocity model.

1929

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

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

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

1933 ``(lat, lon)``

1934 :param interpolation: interpolation method. Choose from

1935 ``('nearest_neighbor', 'multilinear')``

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

1937 point

1938

1939 The default implementation retrieves and interpolates rho

1940 from the contained 1D velocity profile.

1941 '''

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

1943 parameter='rho',

1944 interpolation=interpolation)

1945

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

1947 interpolation=None):

1948

1949 if interpolation is None:

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

1951 'multilinear', 'nearest_neighbor')

1952

1953 earthmod = self.earthmodel_1d

1954 store_depth_profile = self.get_source_depths()

1955 z_profile = earthmod.profile('z')

1956

1957 if parameter == 'vs':

1958 vs_profile = earthmod.profile('vs')

1959 profile = num.interp(

1960 store_depth_profile, z_profile, vs_profile)

1961

1962 elif parameter == 'vp':

1963 vp_profile = earthmod.profile('vp')

1964 profile = num.interp(

1965 store_depth_profile, z_profile, vp_profile)

1966

1967 elif parameter == 'rho':

1968 rho_profile = earthmod.profile('rho')

1969

1970 profile = num.interp(

1971 store_depth_profile, z_profile, rho_profile)

1972

1973 elif parameter == 'shear_moduli':

1974 vs_profile = earthmod.profile('vs')

1975 rho_profile = earthmod.profile('rho')

1976

1977 store_vs_profile = num.interp(

1978 store_depth_profile, z_profile, vs_profile)

1979 store_rho_profile = num.interp(

1980 store_depth_profile, z_profile, rho_profile)

1981

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

1983

1984 elif parameter == 'lambda_moduli':

1985 vs_profile = earthmod.profile('vs')

1986 vp_profile = earthmod.profile('vp')

1987 rho_profile = earthmod.profile('rho')

1988

1989 store_vs_profile = num.interp(

1990 store_depth_profile, z_profile, vs_profile)

1991 store_vp_profile = num.interp(

1992 store_depth_profile, z_profile, vp_profile)

1993 store_rho_profile = num.interp(

1994 store_depth_profile, z_profile, rho_profile)

1995

1996 profile = store_rho_profile * (

1997 num.power(store_vp_profile, 2) -

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

1999 else:

2000 raise TypeError(

2001 'parameter %s not available' % parameter)

2002

2003 if interpolation == 'multilinear':

2004 kind = 'linear'

2005 elif interpolation == 'nearest_neighbor':

2006 kind = 'nearest'

2007 else:

2008 raise TypeError(

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

2010

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

2012

2013 try:

2014 return interpolator(points[:, 2])

2015 except ValueError:

2016 raise OutOfBounds()

2017

2018 def is_static(self):

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

2020 if self.modelling_code_id.startswith(code):

2021 return True

2022 return False

2023

2024 def is_dynamic(self):

2025 return not self.is_static()

2026

2027 def get_source_depths(self):

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

2029

2030 def get_tabulated_phase(self, phase_id):

2031 '''

2032 Get tabulated phase definition.

2033 '''

2034

2035 for pdef in self.tabulated_phases:

2036 if pdef.id == phase_id:

2037 return pdef

2038

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

2040

2041 def fix_ttt_holes(self, sptree, mode):

2042 raise StoreError(

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

2044 % self.short_type)

2045

2046 @property

2047 def effective_stored_quantity(self):

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

2049 component_scheme_to_description[self.component_scheme] \

2050 .default_stored_quantity

2051

2052

2053class ConfigTypeA(Config):

2054 '''

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

2056

2057 * Problem is invariant to horizontal translations and rotations around

2058 vertical axis.

2059

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

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

2062 component)``

2063

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

2065 points.

2066 '''

2067

2068 receiver_depth = Float.T(

2069 default=0.0,

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

2071

2072 source_depth_min = Float.T(

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

2074

2075 source_depth_max = Float.T(

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

2077

2078 source_depth_delta = Float.T(

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

2080

2081 distance_min = Float.T(

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

2083

2084 distance_max = Float.T(

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

2086

2087 distance_delta = Float.T(

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

2089

2090 fd_distance_delta = Float.T(

2091 optional=True,

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

2093

2094 short_type = 'A'

2095

2096 provided_schemes = [

2097 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10',

2098 'rotational8']

2099

2100 def get_surface_distance(self, args):

2101 return args[1]

2102

2103 def get_distance(self, args):

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

2105

2106 def get_source_depth(self, args):

2107 return args[0]

2108

2109 def get_source_depths(self):

2110 return self.coords[0]

2111

2112 def get_receiver_depth(self, args):

2113 return self.receiver_depth

2114

2115 def _update(self):

2116 self.mins = num.array(

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

2118 self.maxs = num.array(

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

2120 self.deltas = num.array(

2121 [self.source_depth_delta, self.distance_delta],

2122 dtype=float)

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

2124 vicinity_eps).astype(int) + 1

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

2126 self.deltat = 1.0/self.sample_rate

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

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

2129 (mi, ma, n) in

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

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

2132

2133 self.nsource_depths, self.ndistances = self.ns

2134

2135 def _make_index_functions(self):

2136

2137 amin, bmin = self.mins

2138 da, db = self.deltas

2139 na, nb = self.ns

2140

2141 ng = self.ncomponents

2142

2143 def index_function(a, b, ig):

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

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

2146 try:

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

2148 except ValueError:

2149 raise OutOfBounds()

2150

2151 def indices_function(a, b, ig):

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

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

2154 try:

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

2156 except ValueError:

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

2158 try:

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

2160 except ValueError:

2161 raise OutOfBounds()

2162

2163 def grid_interpolation_coefficients(a, b):

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

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

2166 return ias, ibs

2167

2168 def vicinity_function(a, b, ig):

2169 ias, ibs = grid_interpolation_coefficients(a, b)

2170

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

2172 raise OutOfBounds()

2173

2174 indis = []

2175 weights = []

2176 for ia, va in ias:

2177 iia = ia*nb*ng

2178 for ib, vb in ibs:

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

2180 weights.append(va*vb)

2181

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

2183

2184 def vicinities_function(a, b, ig):

2185

2186 xa = (a - amin) / da

2187 xb = (b - bmin) / db

2188

2189 xa_fl = num.floor(xa)

2190 xa_ce = num.ceil(xa)

2191 xb_fl = num.floor(xb)

2192 xb_ce = num.ceil(xb)

2193 va_fl = 1.0 - (xa - xa_fl)

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

2195 vb_fl = 1.0 - (xb - xb_fl)

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

2197

2198 ia_fl = xa_fl.astype(int)

2199 ia_ce = xa_ce.astype(int)

2200 ib_fl = xb_fl.astype(int)

2201 ib_ce = xb_ce.astype(int)

2202

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

2204 raise OutOfBounds()

2205

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

2207 raise OutOfBounds()

2208

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

2210 raise OutOfBounds()

2211

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

2213 raise OutOfBounds()

2214

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

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

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

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

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

2220

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

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

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

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

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

2226

2227 return irecords, weights

2228

2229 self._index_function = index_function

2230 self._indices_function = indices_function

2231 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2232 self._vicinity_function = vicinity_function

2233 self._vicinities_function = vicinities_function

2234

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

2236 nc = icomponents.size

2237 dists = source.distances_to(receiver)

2238 n = dists.size

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

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

2241 icomponents)

2242

2243 def make_indexing_args1(self, source, receiver):

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

2245

2246 @property

2247 def short_extent(self):

2248 return '%g:%g:%g x %g:%g:%g' % (

2249 self.source_depth_min/km,

2250 self.source_depth_max/km,

2251 self.source_depth_delta/km,

2252 self.distance_min/km,

2253 self.distance_max/km,

2254 self.distance_delta/km)

2255

2256 def fix_ttt_holes(self, sptree, mode):

2257 from pyrocko import eikonal_ext, spit

2258

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

2260

2261 delta = self.deltas[-1]

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

2263

2264 nsources, ndistances = self.ns

2265

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

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

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

2269

2270 speeds = self.get_material_property(

2271 0., 0., points,

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

2273 interpolation='multilinear')

2274

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

2276

2277 times = sptree.interpolate_many(nodes)

2278

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

2280 times = times.reshape(speeds.shape)

2281

2282 try:

2283 eikonal_ext.eikonal_solver_fmm_cartesian(

2284 speeds, times, delta)

2285 except eikonal_ext.EikonalExtError as e:

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

2287 logger.debug(

2288 'Got a warning from eikonal solver '

2289 '- may be ok...')

2290 else:

2291 raise

2292

2293 def func(x):

2294 ibs, ics = \

2295 self.grid_interpolation_coefficients(*x)

2296

2297 t = 0

2298 for ib, vb in ibs:

2299 for ic, vc in ics:

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

2301

2302 return t

2303

2304 return spit.SPTree(

2305 f=func,

2306 ftol=sptree.ftol,

2307 xbounds=sptree.xbounds,

2308 xtols=sptree.xtols)

2309

2310

2311class ConfigTypeB(Config):

2312 '''

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

2314

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

2316 receiver depth

2317

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

2319 receiver_distance, component)``

2320 '''

2321

2322 receiver_depth_min = Float.T(

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

2324

2325 receiver_depth_max = Float.T(

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

2327

2328 receiver_depth_delta = Float.T(

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

2330

2331 source_depth_min = Float.T(

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

2333

2334 source_depth_max = Float.T(

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

2336

2337 source_depth_delta = Float.T(

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

2339

2340 distance_min = Float.T(

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

2342

2343 distance_max = Float.T(

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

2345

2346 distance_delta = Float.T(

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

2348

2349 fd_distance_delta = Float.T(

2350 optional=True,

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

2352

2353 fd_receiver_depth_delta = Float.T(

2354 optional=True,

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

2356

2357 short_type = 'B'

2358

2359 provided_schemes = [

2360 'elastic2', 'elastic5', 'elastic8', 'elastic10', 'poroelastic10',

2361 'rotational8']

2362

2363 def get_distance(self, args):

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

2365

2366 def get_surface_distance(self, args):

2367 return args[2]

2368

2369 def get_source_depth(self, args):

2370 return args[1]

2371

2372 def get_receiver_depth(self, args):

2373 return args[0]

2374

2375 def get_source_depths(self):

2376 return self.coords[1]

2377

2378 def _update(self):

2379 self.mins = num.array([

2380 self.receiver_depth_min,

2381 self.source_depth_min,

2382 self.distance_min],

2383 dtype=float)

2384

2385 self.maxs = num.array([

2386 self.receiver_depth_max,

2387 self.source_depth_max,

2388 self.distance_max],

2389 dtype=float)

2390

2391 self.deltas = num.array([

2392 self.receiver_depth_delta,

2393 self.source_depth_delta,

2394 self.distance_delta],

2395 dtype=float)

2396

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

2398 vicinity_eps).astype(int) + 1

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

2400 self.deltat = 1.0/self.sample_rate

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

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

2403 (mi, ma, n) in

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

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

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

2407

2408 def _make_index_functions(self):

2409

2410 amin, bmin, cmin = self.mins

2411 da, db, dc = self.deltas

2412 na, nb, nc = self.ns

2413 ng = self.ncomponents

2414

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

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

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

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

2419 try:

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

2421 (na, nb, nc, ng))

2422 except ValueError:

2423 raise OutOfBounds()

2424

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

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

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

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

2429 try:

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

2431 (na, nb, nc, ng))

2432 except ValueError:

2433 raise OutOfBounds()

2434

2435 def grid_interpolation_coefficients(a, b, c):

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

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

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

2439 return ias, ibs, ics

2440

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

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

2443

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

2445 raise OutOfBounds()

2446

2447 indis = []

2448 weights = []

2449 for ia, va in ias:

2450 iia = ia*nb*nc*ng

2451 for ib, vb in ibs:

2452 iib = ib*nc*ng

2453 for ic, vc in ics:

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

2455 weights.append(va*vb*vc)

2456

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

2458

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

2460

2461 xa = (a - amin) / da

2462 xb = (b - bmin) / db

2463 xc = (c - cmin) / dc

2464

2465 xa_fl = num.floor(xa)

2466 xa_ce = num.ceil(xa)

2467 xb_fl = num.floor(xb)

2468 xb_ce = num.ceil(xb)

2469 xc_fl = num.floor(xc)

2470 xc_ce = num.ceil(xc)

2471 va_fl = 1.0 - (xa - xa_fl)

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

2473 vb_fl = 1.0 - (xb - xb_fl)

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

2475 vc_fl = 1.0 - (xc - xc_fl)

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

2477

2478 ia_fl = xa_fl.astype(int)

2479 ia_ce = xa_ce.astype(int)

2480 ib_fl = xb_fl.astype(int)

2481 ib_ce = xb_ce.astype(int)

2482 ic_fl = xc_fl.astype(int)

2483 ic_ce = xc_ce.astype(int)

2484

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

2486 raise OutOfBounds()

2487

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

2489 raise OutOfBounds()

2490

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

2492 raise OutOfBounds()

2493

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

2495 raise OutOfBounds()

2496

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

2498 raise OutOfBounds()

2499

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

2501 raise OutOfBounds()

2502

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

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

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

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

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

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

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

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

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

2512

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

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

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

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

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

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

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

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

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

2522

2523 return irecords, weights

2524

2525 self._index_function = index_function

2526 self._indices_function = indices_function

2527 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2528 self._vicinity_function = vicinity_function

2529 self._vicinities_function = vicinities_function

2530

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

2532 nc = icomponents.size

2533 dists = source.distances_to(receiver)

2534 n = dists.size

2535 receiver_depths = num.empty(nc)

2536 receiver_depths.fill(receiver.depth)

2537 return (receiver_depths,

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

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

2540 icomponents)

2541

2542 def make_indexing_args1(self, source, receiver):

2543 return (receiver.depth,

2544 source.depth,

2545 source.distance_to(receiver))

2546

2547 @property

2548 def short_extent(self):

2549 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2550 self.receiver_depth_min/km,

2551 self.receiver_depth_max/km,

2552 self.receiver_depth_delta/km,

2553 self.source_depth_min/km,

2554 self.source_depth_max/km,

2555 self.source_depth_delta/km,

2556 self.distance_min/km,

2557 self.distance_max/km,

2558 self.distance_delta/km)

2559

2560 def fix_ttt_holes(self, sptree, mode):

2561 from pyrocko import eikonal_ext, spit

2562

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

2564

2565 delta = self.deltas[-1]

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

2567

2568 nreceivers, nsources, ndistances = self.ns

2569

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

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

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

2573

2574 speeds = self.get_material_property(

2575 0., 0., points,

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

2577 interpolation='multilinear')

2578

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

2580

2581 receiver_times = []

2582 for ireceiver in range(nreceivers):

2583 nodes = num.hstack([

2584 num.full(

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

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

2587 nodes_sr])

2588

2589 times = sptree.interpolate_many(nodes)

2590

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

2592

2593 times = times.reshape(speeds.shape)

2594

2595 try:

2596 eikonal_ext.eikonal_solver_fmm_cartesian(

2597 speeds, times, delta)

2598 except eikonal_ext.EikonalExtError as e:

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

2600 logger.debug(

2601 'Got a warning from eikonal solver '

2602 '- may be ok...')

2603 else:

2604 raise

2605

2606 receiver_times.append(times)

2607

2608 def func(x):

2609 ias, ibs, ics = \

2610 self.grid_interpolation_coefficients(*x)

2611

2612 t = 0

2613 for ia, va in ias:

2614 times = receiver_times[ia]

2615 for ib, vb in ibs:

2616 for ic, vc in ics:

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

2618

2619 return t

2620

2621 return spit.SPTree(

2622 f=func,

2623 ftol=sptree.ftol,

2624 xbounds=sptree.xbounds,

2625 xtols=sptree.xtols)

2626

2627

2628class ConfigTypeC(Config):

2629 '''

2630 No symmetrical constraints, one fixed receiver position.

2631

2632 * Cartesian 3D source volume around a reference point

2633

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

2635 source_east_shift, source_north_shift, component)``

2636 '''

2637

2638 receiver = Receiver.T(

2639 help='Receiver location')

2640

2641 source_origin = Location.T(

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

2643

2644 source_depth_min = Float.T(

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

2646

2647 source_depth_max = Float.T(

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

2649

2650 source_depth_delta = Float.T(

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

2652

2653 source_east_shift_min = Float.T(

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

2655

2656 source_east_shift_max = Float.T(

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

2658

2659 source_east_shift_delta = Float.T(

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

2661

2662 source_north_shift_min = Float.T(

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

2664

2665 source_north_shift_max = Float.T(

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

2667

2668 source_north_shift_delta = Float.T(

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

2670

2671 short_type = 'C'

2672

2673 provided_schemes = ['elastic18']

2674

2675 def get_surface_distance(self, args):

2676 _, source_east_shift, source_north_shift, _ = args

2677 sorig = self.source_origin

2678 sloc = Location(

2679 lat=sorig.lat,

2680 lon=sorig.lon,

2681 north_shift=sorig.north_shift + source_north_shift,

2682 east_shift=sorig.east_shift + source_east_shift)

2683

2684 return self.receiver.distance_to(sloc)

2685

2686 def get_distance(self, args):

2687 sdepth, source_east_shift, source_north_shift, _ = args

2688 sorig = self.source_origin

2689 sloc = Location(

2690 lat=sorig.lat,

2691 lon=sorig.lon,

2692 north_shift=sorig.north_shift + source_north_shift,

2693 east_shift=sorig.east_shift + source_east_shift)

2694

2695 return self.receiver.distance_3d_to(sloc)

2696

2697 def get_source_depth(self, args):

2698 return args[0]

2699

2700 def get_receiver_depth(self, args):

2701 return self.receiver.depth

2702

2703 def get_source_depths(self):

2704 return self.coords[0]

2705

2706 def _update(self):

2707 self.mins = num.array([

2708 self.source_depth_min,

2709 self.source_east_shift_min,

2710 self.source_north_shift_min],

2711 dtype=float)

2712

2713 self.maxs = num.array([

2714 self.source_depth_max,

2715 self.source_east_shift_max,

2716 self.source_north_shift_max],

2717 dtype=float)

2718

2719 self.deltas = num.array([

2720 self.source_depth_delta,

2721 self.source_east_shift_delta,

2722 self.source_north_shift_delta],

2723 dtype=float)

2724

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

2726 vicinity_eps).astype(int) + 1

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

2728 self.deltat = 1.0/self.sample_rate

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

2730

2731 self.coords = tuple(

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

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

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

2735

2736 def _make_index_functions(self):

2737

2738 amin, bmin, cmin = self.mins

2739 da, db, dc = self.deltas

2740 na, nb, nc = self.ns

2741 ng = self.ncomponents

2742

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

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

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

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

2747 try:

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

2749 (na, nb, nc, ng))

2750 except ValueError:

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

2752

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

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

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

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

2757

2758 try:

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

2760 (na, nb, nc, ng))

2761 except ValueError:

2762 raise OutOfBounds()

2763

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

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

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

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

2768

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

2770 raise OutOfBounds()

2771

2772 indis = []

2773 weights = []

2774 for ia, va in ias:

2775 iia = ia*nb*nc*ng

2776 for ib, vb in ibs:

2777 iib = ib*nc*ng

2778 for ic, vc in ics:

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

2780 weights.append(va*vb*vc)

2781

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

2783

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

2785

2786 xa = (a-amin) / da

2787 xb = (b-bmin) / db

2788 xc = (c-cmin) / dc

2789

2790 xa_fl = num.floor(xa)

2791 xa_ce = num.ceil(xa)

2792 xb_fl = num.floor(xb)

2793 xb_ce = num.ceil(xb)

2794 xc_fl = num.floor(xc)

2795 xc_ce = num.ceil(xc)

2796 va_fl = 1.0 - (xa - xa_fl)

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

2798 vb_fl = 1.0 - (xb - xb_fl)

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

2800 vc_fl = 1.0 - (xc - xc_fl)

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

2802

2803 ia_fl = xa_fl.astype(int)

2804 ia_ce = xa_ce.astype(int)

2805 ib_fl = xb_fl.astype(int)

2806 ib_ce = xb_ce.astype(int)

2807 ic_fl = xc_fl.astype(int)

2808 ic_ce = xc_ce.astype(int)

2809

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

2811 raise OutOfBounds()

2812

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

2814 raise OutOfBounds()

2815

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

2817 raise OutOfBounds()

2818

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

2820 raise OutOfBounds()

2821

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

2823 raise OutOfBounds()

2824

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

2826 raise OutOfBounds()

2827

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

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

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

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

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

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

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

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

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

2837

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

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

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

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

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

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

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

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

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

2847

2848 return irecords, weights

2849

2850 self._index_function = index_function

2851 self._indices_function = indices_function

2852 self._vicinity_function = vicinity_function

2853 self._vicinities_function = vicinities_function

2854

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

2856 nc = icomponents.size

2857

2858 dists = source.distances_to(self.source_origin)

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

2860

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

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

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

2864

2865 n = dists.size

2866

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

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

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

2870 icomponents)

2871

2872 def make_indexing_args1(self, source, receiver):

2873 dist = source.distance_to(self.source_origin)

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

2875

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

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

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

2879

2880 return (source_depth,

2881 source_east_shift,

2882 source_north_shift)

2883

2884 @property

2885 def short_extent(self):

2886 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2887 self.source_depth_min/km,

2888 self.source_depth_max/km,

2889 self.source_depth_delta/km,

2890 self.source_east_shift_min/km,

2891 self.source_east_shift_max/km,

2892 self.source_east_shift_delta/km,

2893 self.source_north_shift_min/km,

2894 self.source_north_shift_max/km,

2895 self.source_north_shift_delta/km)

2896

2897

2898class Weighting(Object):

2899 factor = Float.T(default=1.0)

2900

2901

2902class Taper(Object):

2903 tmin = Timing.T()

2904 tmax = Timing.T()

2905 tfade = Float.T(default=0.0)

2906 shape = StringChoice.T(

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

2908 default='cos',

2909 optional=True)

2910

2911

2912class SimplePattern(SObject):

2913

2914 _pool = {}

2915

2916 def __init__(self, pattern):

2917 self._pattern = pattern

2918 SObject.__init__(self)

2919

2920 def __str__(self):

2921 return self._pattern

2922

2923 @property

2924 def regex(self):

2925 pool = SimplePattern._pool

2926 if self.pattern not in pool:

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

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

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

2930

2931 return pool[self.pattern]

2932

2933 def match(self, s):

2934 return self.regex.match(s)

2935

2936

2937class WaveformType(StringChoice):

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

2939 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',

2940 'envelope_dis', 'envelope_vel', 'envelope_acc']

2941

2942

2943class ChannelSelection(Object):

2944 pattern = SimplePattern.T(optional=True)

2945 min_sample_rate = Float.T(optional=True)

2946 max_sample_rate = Float.T(optional=True)

2947

2948

2949class StationSelection(Object):

2950 includes = SimplePattern.T()

2951 excludes = SimplePattern.T()

2952 distance_min = Float.T(optional=True)

2953 distance_max = Float.T(optional=True)

2954 azimuth_min = Float.T(optional=True)

2955 azimuth_max = Float.T(optional=True)

2956

2957

2958class WaveformSelection(Object):

2959 channel_selection = ChannelSelection.T(optional=True)

2960 station_selection = StationSelection.T(optional=True)

2961 taper = Taper.T()

2962 # filter = FrequencyResponse.T()

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

2964 weighting = Weighting.T(optional=True)

2965 sample_rate = Float.T(optional=True)

2966 gf_store_id = StringID.T(optional=True)

2967

2968

2969def indi12(x, n):

2970 '''

2971 Get linear interpolation index and weight.

2972 '''

2973

2974 r = round(x)

2975 if abs(r - x) < vicinity_eps:

2976 i = int(r)

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

2978 raise OutOfBounds()

2979

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

2981 else:

2982 f = math.floor(x)

2983 i = int(f)

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

2985 raise OutOfBounds()

2986

2987 v = x-f

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

2989

2990

2991def float_or_none(s):

2992 units = {

2993 'k': 1e3,

2994 'M': 1e6,

2995 }

2996

2997 factor = 1.0

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

2999 factor = units[s[-1]]

3000 s = s[:-1]

3001 if not s:

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

3003

3004 if s:

3005 return float(s) * factor

3006 else:

3007 return None

3008

3009

3010class GridSpecError(Exception):

3011 def __init__(self, s):

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

3013

3014

3015def parse_grid_spec(spec):

3016 try:

3017 result = []

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

3019 t = dspec.split('@')

3020 num = start = stop = step = None

3021 if len(t) == 2:

3022 num = int(t[1])

3023 if num <= 0:

3024 raise GridSpecError(spec)

3025

3026 elif len(t) > 2:

3027 raise GridSpecError(spec)

3028

3029 s = t[0]

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

3031 if len(v) == 1:

3032 start = stop = v[0]

3033 if len(v) >= 2:

3034 start, stop = v[0:2]

3035 if len(v) == 3:

3036 step = v[2]

3037

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

3039 raise GridSpecError(spec)

3040

3041 if step == 0.0:

3042 raise GridSpecError(spec)

3043

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

3045

3046 except ValueError:

3047 raise GridSpecError(spec)

3048

3049 return result

3050

3051

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

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

3054 if start is None:

3055 start = ma if swap else mi

3056 if stop is None:

3057 stop = mi if swap else ma

3058 if step is None:

3059 step = -inc if ma < mi else inc

3060 if num is None:

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

3062 raise GridSpecError()

3063

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

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

3066 if abs(stop-stop2) > eps:

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

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

3069 else:

3070 stop = stop2

3071

3072 if start == stop:

3073 num = 1

3074

3075 return start, stop, num

3076

3077

3078def nditer_outer(x):

3079 return num.nditer(

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

3081

3082

3083def nodes(xs):

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

3085 nnodes = num.prod(ns)

3086 ndim = len(xs)

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

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

3089 x = xs[idim]

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

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

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

3093

3094 return nodes

3095

3096

3097def filledi(x, n):

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

3099 a.fill(x)

3100 return a

3101

3102

3103config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]

3104

3105discretized_source_classes = [

3106 DiscretizedExplosionSource,

3107 DiscretizedSFSource,

3108 DiscretizedMTSource,

3109 DiscretizedPorePressureSource]

3110

3111

3112__all__ = '''

3113Earthmodel1D

3114StringID

3115ScopeType

3116WaveformType

3117QuantityType

3118NearfieldTermsType

3119Reference

3120Region

3121CircularRegion

3122RectangularRegion

3123PhaseSelect

3124InvalidTimingSpecification

3125Timing

3126TPDef

3127OutOfBounds

3128Location

3129Receiver

3130'''.split() + [

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

3132ComponentScheme

3133component_scheme_to_description

3134component_schemes

3135Config

3136GridSpecError

3137Weighting

3138Taper

3139SimplePattern

3140WaveformType

3141ChannelSelection

3142StationSelection

3143WaveformSelection

3144nditer_outer

3145dump

3146load

3147discretized_source_classes

3148config_type_classes

3149UnavailableScheme

3150InterpolationMethod

3151SeismosizerTrace

3152SeismosizerResult

3153Result

3154StaticResult

3155TimingAttributeError

3156'''.split()