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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

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

6'''

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

8'''

10import math

11import re

12import fnmatch

13import logging

15import numpy as num

16from scipy.interpolate import interp1d

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

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

20 List, ValidationError, Timestamp, Tuple, Dict)

21from pyrocko.guts import dump, load # noqa

22from pyrocko.guts_array import literal, Array

23from pyrocko.model import Location, gnss

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

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='elastic18',

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

185 ncomponents=18,

186 default_stored_quantity='displacement',

187 provided_components=[

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

189 ComponentSchemeDescription(

190 name='elastic5',

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

192 ncomponents=5,

193 default_stored_quantity='displacement',

194 provided_components=[

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

196 ComponentSchemeDescription(

197 name='poroelastic10',

198 source_terms=['pore_pressure'],

199 ncomponents=10,

200 default_stored_quantity=None,

201 provided_components=[

202 'displacement.n', 'displacement.e', 'displacement.d',

203 'vertical_tilt.n', 'vertical_tilt.e',

204 'pore_pressure',

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

206

207

208# new names?

209# 'mt_to_displacement_1d'

210# 'mt_to_displacement_1d_ff_only'

211# 'mt_to_gravity_1d'

212# 'mt_to_stress_1d'

213# 'explosion_to_displacement_1d'

214# 'sf_to_displacement_1d'

215# 'mt_to_displacement_3d'

216# 'mt_to_gravity_3d'

217# 'mt_to_stress_3d'

218# 'pore_pressure_to_displacement_1d'

219# 'pore_pressure_to_vertical_tilt_1d'

220# 'pore_pressure_to_pore_pressure_1d'

221# 'pore_pressure_to_darcy_velocity_1d'

222

223

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

225component_scheme_to_description = dict(

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

227

228

229class ComponentScheme(StringChoice):

230 '''

231 Different Green's Function component schemes are available:

232

233 ================= =========================================================

234 Name Description

235 ================= =========================================================

236 ``elastic10`` Elastodynamic for

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

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

239 sources only

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

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

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

243 MT sources only

244 ``elastic2`` Elastodynamic for

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

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

247 isotropic sources only

248 ``elastic5`` Elastodynamic for

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

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

251 sources only

252 ``elastic18`` Elastodynamic for

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

254 sources only

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

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

257 ================= =========================================================

258 '''

259

260 choices = component_schemes

261

262

263class Earthmodel1D(Object):

264 dummy_for = cake.LayeredModel

265 dummy_for_description = 'pyrocko.cake.LayeredModel'

266

267 class __T(TBase):

268 def regularize_extra(self, val):

269 if isinstance(val, str):

270 val = cake.LayeredModel.from_scanlines(

271 cake.read_nd_model_str(val))

272

273 return val

274

275 def to_save(self, val):

276 return literal(cake.write_nd_model_str(val))

277

278

279class StringID(StringPattern):

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

281

282

283class ScopeType(StringChoice):

284 choices = [

285 'global',

286 'regional',

287 'local',

288 ]

289

290

291class WaveType(StringChoice):

292 choices = [

293 'full waveform',

294 'bodywave',

295 'P wave',

296 'S wave',

297 'surface wave',

298 ]

299

300

301class NearfieldTermsType(StringChoice):

302 choices = [

303 'complete',

304 'incomplete',

305 'missing',

306 ]

307

308

309class QuantityType(StringChoice):

310 choices = [

311 'displacement',

312 'rotation',

313 'velocity',

314 'acceleration',

315 'pressure',

316 'tilt',

317 'pore_pressure',

318 'darcy_velocity',

319 'vertical_tilt']

320

321

322class Reference(Object):

323 id = StringID.T()

324 type = String.T()

325 title = Unicode.T()

326 journal = Unicode.T(optional=True)

327 volume = Unicode.T(optional=True)

328 number = Unicode.T(optional=True)

329 pages = Unicode.T(optional=True)

330 year = String.T()

331 note = Unicode.T(optional=True)

332 issn = String.T(optional=True)

333 doi = String.T(optional=True)

334 url = String.T(optional=True)

335 eprint = String.T(optional=True)

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

337 publisher = Unicode.T(optional=True)

338 keywords = Unicode.T(optional=True)

339 note = Unicode.T(optional=True)

340 abstract = Unicode.T(optional=True)

341

342 @classmethod

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

344

345 from pybtex.database.input import bibtex

346

347 parser = bibtex.Parser()

348

349 if filename is not None:

350 bib_data = parser.parse_file(filename)

351 elif stream is not None:

352 bib_data = parser.parse_stream(stream)

353

354 references = []

355

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

357 d = {}

358 avail = entry.fields.keys()

359 for prop in cls.T.properties:

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

361 continue

362

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

364

365 if 'author' in entry.persons:

366 d['authors'] = []

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

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

369

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

371 references.append(c)

372

373 return references

374

375

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

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

378_cake_pat = cake.PhaseDef.allowed_characters_pattern

379_iaspei_pat = cake.PhaseDef.allowed_characters_pattern_classic

380

381_ppat = r'(' \

382 r'cake:' + _cake_pat + \

383 r'|iaspei:' + _iaspei_pat + \

384 r'|vel_surface:' + _fpat + \

385 r'|vel:' + _fpat + \

386 r'|stored:' + _spat + \

387 r')'

388

389

390timing_regex_old = re.compile(

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

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

393

394

395timing_regex = re.compile(

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

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

398

399

400class PhaseSelect(StringChoice):

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

402

403

404class InvalidTimingSpecification(ValidationError):

405 pass

406

407

408class TimingAttributeError(Exception):

409 pass

410

411

412class Timing(SObject):

413 '''

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

415

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

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

418 arrival or group of such arrivals.

419

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

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

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

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

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

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

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

427

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

429

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

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

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

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

434 velocity [km/s]

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

436

437 **Examples:**

438

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

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

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

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

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

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

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

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

447 undefined for a given geometry

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

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

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

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

452 arrivals A, B, and C

453 '''

454

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

456

457 if s is not None:

458 offset_is = None

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

460 try:

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

462

463 if s.endswith('S'):

464 offset_is = 'slowness'

465

466 phase_defs = []

467 select = ''

468

469 except ValueError:

470 matched = False

471 m = timing_regex.match(s)

472 if m:

473 if m.group(3):

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

475 elif m.group(19):

476 phase_defs = [m.group(19)]

477 else:

478 phase_defs = []

479

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

481

482 offset = 0.0

483 soff = m.group(27)

484 if soff:

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

486 if soff.endswith('S'):

487 offset_is = 'slowness'

488 elif soff.endswith('%'):

489 offset_is = 'percent'

490

491 matched = True

492

493 else:

494 m = timing_regex_old.match(s)

495 if m:

496 if m.group(3):

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

498 elif m.group(5):

499 phase_defs = [m.group(5)]

500 else:

501 phase_defs = []

502

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

504

505 offset = 0.0

506 if m.group(6):

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

508

509 phase_defs = [

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

511

512 matched = True

513

514 if not matched:

515 raise InvalidTimingSpecification(s)

516

517 kwargs = dict(

518 phase_defs=phase_defs,

519 select=select,

520 offset=offset,

521 offset_is=offset_is)

522

523 SObject.__init__(self, **kwargs)

524

525 def __str__(self):

526 s = []

527 if self.phase_defs:

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

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

530 sphases = '{%s}' % sphases

531

532 if self.select:

533 sphases = self.select + sphases

534

535 s.append(sphases)

536

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

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

539 if self.offset_is == 'slowness':

540 s.append('S')

541 elif self.offset_is == 'percent':

542 s.append('%')

543

544 return ''.join(s)

545

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

547 try:

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

549 phase_offset = get_phase(

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

551 offset = phase_offset(args)

552 else:

553 offset = self.offset

554

555 if self.phase_defs:

556 extra_args = ()

557 if attributes:

558 extra_args = (attributes,)

559

560 phases = [

561 get_phase(phase_def, *extra_args)

562 for phase_def in self.phase_defs]

563

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

565 results = [

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

567 for result in results]

568

569 results = [

570 result for result in results

571 if result[0] is not None]

572

573 if self.select == 'first':

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

575 elif self.select == 'last':

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

577

578 if not results:

579 if attributes is not None:

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

581 else:

582 return None

583

584 else:

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

586 if self.offset_is == 'percent':

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

588 else:

589 t = t + offset

590

591 if attributes is not None:

592 return (t, ) + values

593 else:

594 return t

595

596 else:

597 if attributes is not None:

598 raise TimingAttributeError(

599 'Cannot get phase attributes from offset-only '

600 'definition.')

601 return offset

602

603 except spit.OutOfBounds:

604 raise OutOfBounds(args)

605

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

607 offset = Float.T(default=0.0)

608 offset_is = String.T(optional=True)

609 select = PhaseSelect.T(

610 default='',

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

612 tuple(PhaseSelect.choices)))

613

614

615def mkdefs(s):

616 defs = []

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

618 try:

619 defs.append(float(x))

620 except ValueError:

621 if x.startswith('!'):

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

623 else:

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

625

626 return defs

627

628

629class TPDef(Object):

630 '''

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

632 '''

633

634 id = StringID.T(

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

636 definition = String.T(

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

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

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

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

641

642 @property

643 def phases(self):

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

645 if isinstance(x, cake.PhaseDef)]

646

647 @property

648 def horizontal_velocities(self):

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

650

651

652class OutOfBounds(Exception):

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

654 Exception.__init__(self)

655 self.values = values

656 self.reason = reason

657 self.context = None

658

659 def __str__(self):

660 scontext = ''

661 if self.context:

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

663

664 if self.reason:

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

666

667 if self.values:

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

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

670 else:

671 return 'out of bounds%s' % scontext

672

673

674class MultiLocation(Object):

675 '''

676 Unstructured grid of point coordinates.

677 '''

678

679 lats = Array.T(

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

681 help='Latitudes of targets.')

682 lons = Array.T(

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

684 help='Longitude of targets.')

685 north_shifts = Array.T(

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

687 help='North shifts of targets.')

688 east_shifts = Array.T(

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

690 help='East shifts of targets.')

691 elevation = Array.T(

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

693 help='Elevations of targets.')

694

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

696 self._coords5 = None

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

698

699 @property

700 def coords5(self):

701 if self._coords5 is not None:

702 return self._coords5

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

704 self.elevation]

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

706 if not sizes:

707 raise AttributeError('Empty StaticTarget')

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

709 raise AttributeError('Inconsistent coordinate sizes.')

710

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

712 for idx, p in enumerate(props):

713 if p is not None:

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

715

716 return self._coords5

717

718 @property

719 def ncoords(self):

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

721 return 0

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

723

724 def get_latlon(self):

725 '''

726 Get all coordinates as lat/lon.

727

728 :returns: Coordinates in Latitude, Longitude

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

730 '''

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

732 for i in range(self.ncoords):

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

734 return latlons

735

736 def get_corner_coords(self):

737 '''

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

739

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

741 :rtype: tuple

742 '''

743 latlon = self.get_latlon()

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

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

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

747

748

749class Receiver(Location):

750 codes = Tuple.T(

751 3, String.T(),

752 optional=True,

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

754

755 def pyrocko_station(self):

756 from pyrocko import model

757

758 lat, lon = self.effective_latlon

759

760 return model.Station(

761 network=self.codes[0],

762 station=self.codes[1],

763 location=self.codes[2],

764 lat=lat,

765 lon=lon,

766 depth=self.depth)

767

768

769def g(x, d):

770 if x is None:

771 return d

772 else:

773 return x

774

775

776class UnavailableScheme(Exception):

777 '''

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

779 '''

780 pass

781

782

783class InvalidNComponents(Exception):

784 '''

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

786 '''

787 pass

788

789

790class DiscretizedSource(Object):

791 '''

792 Base class for discretized sources.

793

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

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

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

797 source distributions are represented by subclasses of the

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

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

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

801 explosion/implosion type source distributions. The geometry calculations

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

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

804

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

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

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

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

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

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

811 '''

812

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

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

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

816 lat = Float.T(optional=True)

817 lon = Float.T(optional=True)

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

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

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

821 dl = Float.T(optional=True)

822 dw = Float.T(optional=True)

823 nl = Float.T(optional=True)

824 nw = Float.T(optional=True)

825

826 @classmethod

827 def check_scheme(cls, scheme):

828 '''

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

830

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

832 supported by this discretized source class.

833 '''

834

835 if scheme not in cls.provided_schemes:

836 raise UnavailableScheme(

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

838 (cls.__name__, scheme))

839

840 def __init__(self, **kwargs):

841 Object.__init__(self, **kwargs)

842 self._latlons = None

843

844 def __setattr__(self, name, value):

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

846 'lats', 'lons'):

847 self.__dict__['_latlons'] = None

848

849 Object.__setattr__(self, name, value)

850

851 def get_source_terms(self, scheme):

852 raise NotImplementedError()

853

854 def make_weights(self, receiver, scheme):

855 raise NotImplementedError()

856

857 @property

858 def effective_latlons(self):

859 '''

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

861 '''

862

863 if self._latlons is None:

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

865 if (self.north_shifts is not None and

866 self.east_shifts is not None):

867 self._latlons = orthodrome.ne_to_latlon(

868 self.lats, self.lons,

869 self.north_shifts, self.east_shifts)

870 else:

871 self._latlons = self.lats, self.lons

872 else:

873 lat = g(self.lat, 0.0)

874 lon = g(self.lon, 0.0)

875 self._latlons = orthodrome.ne_to_latlon(

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

877

878 return self._latlons

879

880 @property

881 def effective_north_shifts(self):

882 if self.north_shifts is not None:

883 return self.north_shifts

884 else:

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

886

887 @property

888 def effective_east_shifts(self):

889 if self.east_shifts is not None:

890 return self.east_shifts

891 else:

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

893

894 def same_origin(self, receiver):

895 '''

896 Check if receiver has same reference point.

897 '''

898

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

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

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

902

903 def azibazis_to(self, receiver):

904 '''

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

906 points.

907 '''

908

909 if self.same_origin(receiver):

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

911 receiver.north_shift - self.north_shifts)

912 bazis = azis + 180.

913 else:

914 slats, slons = self.effective_latlons

915 rlat, rlon = receiver.effective_latlon

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

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

918

919 return azis, bazis

920

921 def distances_to(self, receiver):

922 '''

923 Compute distances to receiver for all contained points.

924 '''

925 if self.same_origin(receiver):

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

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

928

929 else:

930 slats, slons = self.effective_latlons

931 rlat, rlon = receiver.effective_latlon

932 return orthodrome.distance_accurate50m_numpy(slats, slons,

933 rlat, rlon)

934

935 def element_coords(self, i):

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

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

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

939 else:

940 lat = self.lat

941 lon = self.lon

942

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

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

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

946

947 else:

948 north_shift = east_shift = 0.0

949

950 return lat, lon, north_shift, east_shift

951

952 def coords5(self):

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

954

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

956 xs[:, 0] = self.lats

957 xs[:, 1] = self.lons

958 else:

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

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

961

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

963 xs[:, 2] = self.north_shifts

964 xs[:, 3] = self.east_shifts

965

966 xs[:, 4] = self.depths

967

968 return xs

969

970 @property

971 def nelements(self):

972 return self.times.size

973

974 @classmethod

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

976 '''

977 Combine several discretized source models.

978

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

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

981 factors and reference times of the parameterized (undiscretized)

982 sources match or are accounted for.

983 '''

984

985 first = sources[0]

986

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

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

989 'sources of same type.')

990

991 latlons = []

992 for s in sources:

993 latlons.append(s.effective_latlons)

994

995 lats, lons = num.hstack(latlons)

996

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

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

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

1000 same_ref = num.all(

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

1002 else:

1003 same_ref = False

1004

1005 cat = num.concatenate

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

1007 print(times)

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

1009

1010 if same_ref:

1011 lat = first.lat

1012 lon = first.lon

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

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

1015 lats = None

1016 lons = None

1017 else:

1018 lat = None

1019 lon = None

1020 north_shifts = None

1021 east_shifts = None

1022

1023 return cls(

1024 times=times, lat=lat, lon=lon, lats=lats, lons=lons,

1025 north_shifts=north_shifts, east_shifts=east_shifts,

1026 depths=depths, **kwargs)

1027

1028 def centroid_position(self):

1029 moments = self.moments()

1030 norm = num.sum(moments)

1031 if norm != 0.0:

1032 w = moments / num.sum(moments)

1033 else:

1034 w = num.ones(moments.size)

1035

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

1037 lats, lons = self.effective_latlons

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

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

1040 else:

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

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

1043

1044 cn = num.sum(n*w)

1045 ce = num.sum(e*w)

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

1047

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

1049 lat = clat

1050 lon = clon

1051 north_shift = 0.

1052 east_shift = 0.

1053 else:

1054 lat = g(self.lat, 0.0)

1055 lon = g(self.lon, 0.0)

1056 north_shift = cn

1057 east_shift = ce

1058

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

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

1061 return tuple(float(x) for x in

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

1063

1064

1065class DiscretizedExplosionSource(DiscretizedSource):

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

1067

1068 provided_schemes = (

1069 'elastic2',

1070 'elastic8',

1071 'elastic10',

1072 )

1073

1074 def get_source_terms(self, scheme):

1075 self.check_scheme(scheme)

1076

1077 if scheme == 'elastic2':

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

1079

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

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

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

1083 return m6s

1084 else:

1085 assert False

1086

1087 def make_weights(self, receiver, scheme):

1088 self.check_scheme(scheme)

1089

1090 azis, bazis = self.azibazis_to(receiver)

1091

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

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

1094

1095 m0s = self.m0s

1096 n = azis.size

1097

1098 cat = num.concatenate

1099 rep = num.repeat

1100

1101 if scheme == 'elastic2':

1102 w_n = cb*m0s

1103 g_n = filledi(0, n)

1104 w_e = sb*m0s

1105 g_e = filledi(0, n)

1106 w_d = m0s

1107 g_d = filledi(1, n)

1108

1109 elif scheme == 'elastic8':

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

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

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

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

1114 w_d = cat((m0s, m0s))

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

1116

1117 elif scheme == 'elastic10':

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

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

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

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

1122 w_d = cat((m0s, m0s, m0s))

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

1124

1125 else:

1126 assert False

1127

1128 return (

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

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

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

1132 )

1133

1134 def split(self):

1135 from pyrocko.gf.seismosizer import ExplosionSource

1136 sources = []

1137 for i in range(self.nelements):

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

1139 sources.append(ExplosionSource(

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

1141 lat=lat,

1142 lon=lon,

1143 north_shift=north_shift,

1144 east_shift=east_shift,

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

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

1147

1148 return sources

1149

1150 def moments(self):

1151 return self.m0s

1152

1153 def centroid(self):

1154 from pyrocko.gf.seismosizer import ExplosionSource

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

1156 self.centroid_position()

1157

1158 return ExplosionSource(

1159 time=time,

1160 lat=lat,

1161 lon=lon,

1162 north_shift=north_shift,

1163 east_shift=east_shift,

1164 depth=depth,

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

1166

1167 @classmethod

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

1169 '''

1170 Combine several discretized source models.

1171

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

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

1174 factors and reference times of the parameterized (undiscretized)

1175 sources match or are accounted for.

1176 '''

1177

1178 if 'm0s' not in kwargs:

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

1180

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

1182 **kwargs)

1183

1184

1185class DiscretizedSFSource(DiscretizedSource):

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

1187

1188 provided_schemes = (

1189 'elastic5',

1190 )

1191

1192 def get_source_terms(self, scheme):

1193 self.check_scheme(scheme)

1194

1195 return self.forces

1196

1197 def make_weights(self, receiver, scheme):

1198 self.check_scheme(scheme)

1199

1200 azis, bazis = self.azibazis_to(receiver)

1201

1202 sa = num.sin(azis*d2r)

1203 ca = num.cos(azis*d2r)

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

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

1206

1207 forces = self.forces

1208 fn = forces[:, 0]

1209 fe = forces[:, 1]

1210 fd = forces[:, 2]

1211

1212 f0 = fd

1213 f1 = ca * fn + sa * fe

1214 f2 = ca * fe - sa * fn

1215

1216 n = azis.size

1217

1218 if scheme == 'elastic5':

1219 ioff = 0

1220

1221 cat = num.concatenate

1222 rep = num.repeat

1223

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

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

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

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

1228 w_d = cat((f0, f1))

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

1230

1231 return (

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

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

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

1235 )

1236

1237 @classmethod

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

1239 '''

1240 Combine several discretized source models.

1241

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

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

1244 factors and reference times of the parameterized (undiscretized)

1245 sources match or are accounted for.

1246 '''

1247

1248 if 'forces' not in kwargs:

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

1250

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

1252

1253 def moments(self):

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

1255

1256 def centroid(self):

1257 from pyrocko.gf.seismosizer import SFSource

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

1259 self.centroid_position()

1260

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

1262 return SFSource(

1263 time=time,

1264 lat=lat,

1265 lon=lon,

1266 north_shift=north_shift,

1267 east_shift=east_shift,

1268 depth=depth,

1269 fn=fn,

1270 fe=fe,

1271 fd=fd)

1272

1273

1274class DiscretizedMTSource(DiscretizedSource):

1275 m6s = Array.T(

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

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

1278

1279 provided_schemes = (

1280 'elastic8',

1281 'elastic10',

1282 'elastic18',

1283 )

1284

1285 def get_source_terms(self, scheme):

1286 self.check_scheme(scheme)

1287 return self.m6s

1288

1289 def make_weights(self, receiver, scheme):

1290 self.check_scheme(scheme)

1291

1292 m6s = self.m6s

1293 n = m6s.shape[0]

1294 rep = num.repeat

1295

1296 if scheme == 'elastic18':

1297 w_n = m6s.flatten()

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

1299 w_e = m6s.flatten()

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

1301 w_d = m6s.flatten()

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

1303

1304 else:

1305 azis, bazis = self.azibazis_to(receiver)

1306

1307 sa = num.sin(azis*d2r)

1308 ca = num.cos(azis*d2r)

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

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

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

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

1313

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

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

1316 f2 = m6s[:, 2]

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

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

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

1320

1321 cat = num.concatenate

1322

1323 if scheme == 'elastic8':

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

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

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

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

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

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

1330

1331 elif scheme == 'elastic10':

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

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

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

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

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

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

1338

1339 else:

1340 assert False

1341

1342 return (

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

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

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

1346 )

1347

1348 def split(self):

1349 from pyrocko.gf.seismosizer import MTSource

1350 sources = []

1351 for i in range(self.nelements):

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

1353 sources.append(MTSource(

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

1355 lat=lat,

1356 lon=lon,

1357 north_shift=north_shift,

1358 east_shift=east_shift,

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

1360 m6=self.m6s[i]))

1361

1362 return sources

1363

1364 def moments(self):

1365 moments = num.array(

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

1367 for m6 in self.m6s])

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

1369

1370 def get_moment_rate(self, deltat=None):

1371 moments = self.moments()

1372 times = self.times

1373 times -= times.min()

1374

1375 t_max = times.max()

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

1377 mom_times[mom_times > t_max] = t_max

1378

1379 # Right open histrogram bins

1380 mom, _ = num.histogram(

1381 times,

1382 bins=mom_times,

1383 weights=moments)

1384

1385 deltat = num.diff(mom_times)

1386 mom_rate = mom / deltat

1387

1388 return mom_rate, mom_times[1:]

1389

1390 def centroid(self):

1391 from pyrocko.gf.seismosizer import MTSource

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

1393 self.centroid_position()

1394

1395 return MTSource(

1396 time=time,

1397 lat=lat,

1398 lon=lon,

1399 north_shift=north_shift,

1400 east_shift=east_shift,

1401 depth=depth,

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

1403

1404 @classmethod

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

1406 '''

1407 Combine several discretized source models.

1408

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

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

1411 factors and reference times of the parameterized (undiscretized)

1412 sources match or are accounted for.

1413 '''

1414

1415 if 'm6s' not in kwargs:

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

1417

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

1419

1420

1421class DiscretizedPorePressureSource(DiscretizedSource):

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

1423

1424 provided_schemes = (

1425 'poroelastic10',

1426 )

1427

1428 def get_source_terms(self, scheme):

1429 self.check_scheme(scheme)

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

1431

1432 def make_weights(self, receiver, scheme):

1433 self.check_scheme(scheme)

1434

1435 azis, bazis = self.azibazis_to(receiver)

1436

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

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

1439

1440 pp = self.pp

1441 n = bazis.size

1442

1443 w_un = cb*pp

1444 g_un = filledi(1, n)

1445 w_ue = sb*pp

1446 g_ue = filledi(1, n)

1447 w_ud = pp

1448 g_ud = filledi(0, n)

1449

1450 w_tn = cb*pp

1451 g_tn = filledi(6, n)

1452 w_te = sb*pp

1453 g_te = filledi(6, n)

1454

1455 w_pp = pp

1456 g_pp = filledi(7, n)

1457

1458 w_dvn = cb*pp

1459 g_dvn = filledi(9, n)

1460 w_dve = sb*pp

1461 g_dve = filledi(9, n)

1462 w_dvd = pp

1463 g_dvd = filledi(8, n)

1464

1465 return (

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

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

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

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

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

1471 ('pore_pressure', w_pp, g_pp),

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

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

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

1475 )

1476

1477 def moments(self):

1478 return self.pp

1479

1480 def centroid(self):

1481 from pyrocko.gf.seismosizer import PorePressurePointSource

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

1483 self.centroid_position()

1484

1485 return PorePressurePointSource(

1486 time=time,

1487 lat=lat,

1488 lon=lon,

1489 north_shift=north_shift,

1490 east_shift=east_shift,

1491 depth=depth,

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

1493

1494 @classmethod

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

1496 '''

1497 Combine several discretized source models.

1498

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

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

1501 factors and reference times of the parameterized (undiscretized)

1502 sources match or are accounted for.

1503 '''

1504

1505 if 'pp' not in kwargs:

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

1507

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

1509 **kwargs)

1510

1511

1512class Region(Object):

1513 name = String.T(optional=True)

1514

1515

1516class RectangularRegion(Region):

1517 lat_min = Float.T()

1518 lat_max = Float.T()

1519 lon_min = Float.T()

1520 lon_max = Float.T()

1521

1522

1523class CircularRegion(Region):

1524 lat = Float.T()

1525 lon = Float.T()

1526 radius = Float.T()

1527

1528

1529class Config(Object):

1530 '''

1531 Green's function store meta information.

1532

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

1534 configuration types are:

1535

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

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

1538

1539 * Problem is invariant to horizontal translations and rotations around

1540 vertical axis.

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

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

1543 component)``

1544

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

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

1547

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

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

1550 receiver_depth, component)``

1551

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

1553 constraints but fixed receiver positions

1554

1555 * Cartesian source volume around a reference point

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

1557 source_east_shift, source_north_shift, component)``

1558 '''

1559

1560 id = StringID.T(

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

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

1563 'letter.')

1564

1565 derived_from_id = StringID.T(

1566 optional=True,

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

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

1569

1570 version = String.T(

1571 default='1.0',

1572 optional=True,

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

1574

1575 modelling_code_id = StringID.T(

1576 optional=True,

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

1578

1579 author = Unicode.T(

1580 optional=True,

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

1582

1583 author_email = String.T(

1584 optional=True,

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

1586

1587 created_time = Timestamp.T(

1588 optional=True,

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

1590

1591 regions = List.T(

1592 Region.T(),

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

1594

1595 scope_type = ScopeType.T(

1596 optional=True,

1597 help='Distance range scope of the store '

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

1599

1600 waveform_type = WaveType.T(

1601 optional=True,

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

1603

1604 nearfield_terms = NearfieldTermsType.T(

1605 optional=True,

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

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

1608

1609 description = String.T(

1610 optional=True,

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

1612

1613 references = List.T(

1614 Reference.T(),

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

1616 'related work.')

1617

1618 earthmodel_1d = Earthmodel1D.T(

1619 optional=True,

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

1621

1622 earthmodel_receiver_1d = Earthmodel1D.T(

1623 optional=True,

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

1625

1626 can_interpolate_source = Bool.T(

1627 optional=True,

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

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

1630

1631 can_interpolate_receiver = Bool.T(

1632 optional=True,

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

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

1635

1636 frequency_min = Float.T(

1637 optional=True,

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

1639

1640 frequency_max = Float.T(

1641 optional=True,

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

1643

1644 sample_rate = Float.T(

1645 optional=True,

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

1647

1648 factor = Float.T(

1649 default=1.0,

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

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

1652 optional=True)

1653

1654 component_scheme = ComponentScheme.T(

1655 default='elastic10',

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

1657

1658 stored_quantity = QuantityType.T(

1659 optional=True,

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

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

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

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

1664

1665 tabulated_phases = List.T(

1666 TPDef.T(),

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

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

1669

1670 ncomponents = Int.T(

1671 optional=True,

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

1673

1674 uuid = String.T(

1675 optional=True,

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

1677 'GF store for practical purposes.')

1678

1679 reference = String.T(

1680 optional=True,

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

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

1683

1684 def __init__(self, **kwargs):

1685 self._do_auto_updates = False

1686 Object.__init__(self, **kwargs)

1687 self._index_function = None

1688 self._indices_function = None

1689 self._vicinity_function = None

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

1691 self._do_auto_updates = True

1692 self.update()

1693

1694 def check_ncomponents(self):

1695 ncomponents = component_scheme_to_description[

1696 self.component_scheme].ncomponents

1697

1698 if self.ncomponents is None:

1699 self.ncomponents = ncomponents

1700 elif ncomponents != self.ncomponents:

1701 raise InvalidNComponents(

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

1703 self.ncomponents, self.component_scheme))

1704

1705 def __setattr__(self, name, value):

1706 Object.__setattr__(self, name, value)

1707 try:

1708 self.T.get_property(name)

1709 if self._do_auto_updates:

1710 self.update()

1711

1712 except ValueError:

1713 pass

1714

1715 def update(self):

1716 self.check_ncomponents()

1717 self._update()

1718 self._make_index_functions()

1719

1720 def irecord(self, *args):

1721 return self._index_function(*args)

1722

1723 def irecords(self, *args):

1724 return self._indices_function(*args)

1725

1726 def vicinity(self, *args):

1727 return self._vicinity_function(*args)

1728

1729 def vicinities(self, *args):

1730 return self._vicinities_function(*args)

1731

1732 def grid_interpolation_coefficients(self, *args):

1733 return self._grid_interpolation_coefficients(*args)

1734

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

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

1737

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

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

1740

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

1742 i = 0

1743 arrs = []

1744 ntotal = 1

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

1746 if gdef and len(gdef) > i:

1747 sssn = gdef[i]

1748 else:

1749 sssn = (None,)*4

1750

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

1752 ntotal *= len(arr)

1753

1754 arrs.append(arr)

1755 i += 1

1756

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

1758 return nditer_outer(arrs[:level])

1759

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

1761 nthreads=0):

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

1763

1764 out = []

1765 delays = source.times

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

1767 receiver,

1768 self.component_scheme):

1769

1770 weights *= self.factor

1771

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

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

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

1775

1776 return out

1777

1778 def short_info(self):

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

1780

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

1782 interpolation=None):

1783 '''

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

1785

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

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

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

1789 ``(lat, lon)``

1790 :param interpolation: interpolation method. Choose from

1791 ``('nearest_neighbor', 'multilinear')``

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

1793 point

1794

1795 The default implementation retrieves and interpolates the shear moduli

1796 from the contained 1D velocity profile.

1797 '''

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

1799 parameter='shear_moduli',

1800 interpolation=interpolation)

1801

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

1803 interpolation=None):

1804 '''

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

1806

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

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

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

1810 ``(lat, lon)``

1811 :param interpolation: interpolation method. Choose from

1812 ``('nearest_neighbor', 'multilinear')``

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

1814 point

1815

1816 The default implementation retrieves and interpolates the lambda moduli

1817 from the contained 1D velocity profile.

1818 '''

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

1820 parameter='lambda_moduli',

1821 interpolation=interpolation)

1822

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

1824 interpolation=None):

1825 '''

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

1827

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

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

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

1831 ``(lat, lon)``

1832 :param interpolation: interpolation method. Choose from

1833 ``('nearest_neighbor', 'multilinear')``

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

1835 point

1836

1837 The default implementation retrieves and interpolates the lambda moduli

1838 from the contained 1D velocity profile.

1839 '''

1840 lambda_moduli = self.get_material_property(

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

1842 interpolation=interpolation)

1843 shear_moduli = self.get_material_property(

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

1845 interpolation=interpolation)

1846 return lambda_moduli + (2 / 3) * shear_moduli

1847

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

1849 '''

1850 Get Vs at given points from contained velocity model.

1851

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

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

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

1855 ``(lat, lon)``

1856 :param interpolation: interpolation method. Choose from

1857 ``('nearest_neighbor', 'multilinear')``

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

1859 point

1860

1861 The default implementation retrieves and interpolates Vs

1862 from the contained 1D velocity profile.

1863 '''

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

1865 parameter='vs',

1866 interpolation=interpolation)

1867

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

1869 '''

1870 Get Vp at given points from contained velocity model.

1871

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

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

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

1875 ``(lat, lon)``

1876 :param interpolation: interpolation method. Choose from

1877 ``('nearest_neighbor', 'multilinear')``

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

1879 point

1880

1881 The default implementation retrieves and interpolates Vp

1882 from the contained 1D velocity profile.

1883 '''

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

1885 parameter='vp',

1886 interpolation=interpolation)

1887

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

1889 '''

1890 Get rho at given points from contained velocity model.

1891

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

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

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

1895 ``(lat, lon)``

1896 :param interpolation: interpolation method. Choose from

1897 ``('nearest_neighbor', 'multilinear')``

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

1899 point

1900

1901 The default implementation retrieves and interpolates rho

1902 from the contained 1D velocity profile.

1903 '''

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

1905 parameter='rho',

1906 interpolation=interpolation)

1907

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

1909 interpolation=None):

1910

1911 if interpolation is None:

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

1913 'multilinear', 'nearest_neighbor')

1914

1915 earthmod = self.earthmodel_1d

1916 store_depth_profile = self.get_source_depths()

1917 z_profile = earthmod.profile('z')

1918

1919 if parameter == 'vs':

1920 vs_profile = earthmod.profile('vs')

1921 profile = num.interp(

1922 store_depth_profile, z_profile, vs_profile)

1923

1924 elif parameter == 'vp':

1925 vp_profile = earthmod.profile('vp')

1926 profile = num.interp(

1927 store_depth_profile, z_profile, vp_profile)

1928

1929 elif parameter == 'rho':

1930 rho_profile = earthmod.profile('rho')

1931

1932 profile = num.interp(

1933 store_depth_profile, z_profile, rho_profile)

1934

1935 elif parameter == 'shear_moduli':

1936 vs_profile = earthmod.profile('vs')

1937 rho_profile = earthmod.profile('rho')

1938

1939 store_vs_profile = num.interp(

1940 store_depth_profile, z_profile, vs_profile)

1941 store_rho_profile = num.interp(

1942 store_depth_profile, z_profile, rho_profile)

1943

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

1945

1946 elif parameter == 'lambda_moduli':

1947 vs_profile = earthmod.profile('vs')

1948 vp_profile = earthmod.profile('vp')

1949 rho_profile = earthmod.profile('rho')

1950

1951 store_vs_profile = num.interp(

1952 store_depth_profile, z_profile, vs_profile)

1953 store_vp_profile = num.interp(

1954 store_depth_profile, z_profile, vp_profile)

1955 store_rho_profile = num.interp(

1956 store_depth_profile, z_profile, rho_profile)

1957

1958 profile = store_rho_profile * (

1959 num.power(store_vp_profile, 2) -

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

1961 else:

1962 raise TypeError(

1963 'parameter %s not available' % parameter)

1964

1965 if interpolation == 'multilinear':

1966 kind = 'linear'

1967 elif interpolation == 'nearest_neighbor':

1968 kind = 'nearest'

1969 else:

1970 raise TypeError(

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

1972

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

1974

1975 try:

1976 return interpolator(points[:, 2])

1977 except ValueError:

1978 raise OutOfBounds()

1979

1980 def is_static(self):

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

1982 if self.modelling_code_id.startswith(code):

1983 return True

1984 return False

1985

1986 def is_dynamic(self):

1987 return not self.is_static()

1988

1989 def get_source_depths(self):

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

1991

1992 def get_tabulated_phase(self, phase_id):

1993 '''

1994 Get tabulated phase definition.

1995 '''

1996

1997 for pdef in self.tabulated_phases:

1998 if pdef.id == phase_id:

1999 return pdef

2000

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

2002

2003 def fix_ttt_holes(self, sptree, mode):

2004 raise StoreError(

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

2006 % self.short_type)

2007

2008

2009class ConfigTypeA(Config):

2010 '''

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

2012

2013 * Problem is invariant to horizontal translations and rotations around

2014 vertical axis.

2015

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

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

2018 component)``

2019

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

2021 points.

2022 '''

2023

2024 receiver_depth = Float.T(

2025 default=0.0,

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

2027

2028 source_depth_min = Float.T(

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

2030

2031 source_depth_max = Float.T(

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

2033

2034 source_depth_delta = Float.T(

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

2036

2037 distance_min = Float.T(

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

2039

2040 distance_max = Float.T(

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

2042

2043 distance_delta = Float.T(

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

2045

2046 short_type = 'A'

2047

2048 provided_schemes = [

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

2050

2051 def get_surface_distance(self, args):

2052 return args[1]

2053

2054 def get_distance(self, args):

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

2056

2057 def get_source_depth(self, args):

2058 return args[0]

2059

2060 def get_source_depths(self):

2061 return self.coords[0]

2062

2063 def get_receiver_depth(self, args):

2064 return self.receiver_depth

2065

2066 def _update(self):

2067 self.mins = num.array(

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

2069 self.maxs = num.array(

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

2071 self.deltas = num.array(

2072 [self.source_depth_delta, self.distance_delta],

2073 dtype=float)

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

2075 vicinity_eps).astype(int) + 1

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

2077 self.deltat = 1.0/self.sample_rate

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

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

2080 (mi, ma, n) in

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

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

2083

2084 self.nsource_depths, self.ndistances = self.ns

2085

2086 def _make_index_functions(self):

2087

2088 amin, bmin = self.mins

2089 da, db = self.deltas

2090 na, nb = self.ns

2091

2092 ng = self.ncomponents

2093

2094 def index_function(a, b, ig):

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

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

2097 try:

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

2099 except ValueError:

2100 raise OutOfBounds()

2101

2102 def indices_function(a, b, ig):

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

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

2105 try:

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

2107 except ValueError:

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

2109 try:

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

2111 except ValueError:

2112 raise OutOfBounds()

2113

2114 def grid_interpolation_coefficients(a, b):

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

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

2117 return ias, ibs

2118

2119 def vicinity_function(a, b, ig):

2120 ias, ibs = grid_interpolation_coefficients(a, b)

2121

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

2123 raise OutOfBounds()

2124

2125 indis = []

2126 weights = []

2127 for ia, va in ias:

2128 iia = ia*nb*ng

2129 for ib, vb in ibs:

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

2131 weights.append(va*vb)

2132

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

2134

2135 def vicinities_function(a, b, ig):

2136

2137 xa = (a - amin) / da

2138 xb = (b - bmin) / db

2139

2140 xa_fl = num.floor(xa)

2141 xa_ce = num.ceil(xa)

2142 xb_fl = num.floor(xb)

2143 xb_ce = num.ceil(xb)

2144 va_fl = 1.0 - (xa - xa_fl)

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

2146 vb_fl = 1.0 - (xb - xb_fl)

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

2148

2149 ia_fl = xa_fl.astype(int)

2150 ia_ce = xa_ce.astype(int)

2151 ib_fl = xb_fl.astype(int)

2152 ib_ce = xb_ce.astype(int)

2153

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

2155 raise OutOfBounds()

2156

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

2158 raise OutOfBounds()

2159

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

2161 raise OutOfBounds()

2162

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

2164 raise OutOfBounds()

2165

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

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

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

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

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

2171

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

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

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

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

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

2177

2178 return irecords, weights

2179

2180 self._index_function = index_function

2181 self._indices_function = indices_function

2182 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2183 self._vicinity_function = vicinity_function

2184 self._vicinities_function = vicinities_function

2185

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

2187 nc = icomponents.size

2188 dists = source.distances_to(receiver)

2189 n = dists.size

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

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

2192 icomponents)

2193

2194 def make_indexing_args1(self, source, receiver):

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

2196

2197 @property

2198 def short_extent(self):

2199 return '%g:%g:%g x %g:%g:%g' % (

2200 self.source_depth_min/km,

2201 self.source_depth_max/km,

2202 self.source_depth_delta/km,

2203 self.distance_min/km,

2204 self.distance_max/km,

2205 self.distance_delta/km)

2206

2207 def fix_ttt_holes(self, sptree, mode):

2208 from pyrocko import eikonal_ext, spit

2209

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

2211

2212 delta = self.deltas[-1]

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

2214

2215 nsources, ndistances = self.ns

2216

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

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

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

2220

2221 speeds = self.get_material_property(

2222 0., 0., points,

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

2224 interpolation='multilinear')

2225

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

2227

2228 times = sptree.interpolate_many(nodes)

2229

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

2231 times = times.reshape(speeds.shape)

2232

2233 try:

2234 eikonal_ext.eikonal_solver_fmm_cartesian(

2235 speeds, times, delta)

2236 except eikonal_ext.EikonalExtError as e:

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

2238 logger.debug(

2239 'Got a warning from eikonal solver '

2240 '- may be ok...')

2241 else:

2242 raise

2243

2244 def func(x):

2245 ibs, ics = \

2246 self.grid_interpolation_coefficients(*x)

2247

2248 t = 0

2249 for ib, vb in ibs:

2250 for ic, vc in ics:

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

2252

2253 return t

2254

2255 return spit.SPTree(

2256 f=func,

2257 ftol=sptree.ftol,

2258 xbounds=sptree.xbounds,

2259 xtols=sptree.xtols)

2260

2261

2262class ConfigTypeB(Config):

2263 '''

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

2265

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

2267 receiver depth

2268

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

2270 receiver_distance, component)``

2271 '''

2272

2273 receiver_depth_min = Float.T(

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

2275

2276 receiver_depth_max = Float.T(

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

2278

2279 receiver_depth_delta = Float.T(

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

2281

2282 source_depth_min = Float.T(

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

2284

2285 source_depth_max = Float.T(

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

2287

2288 source_depth_delta = Float.T(

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

2290

2291 distance_min = Float.T(

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

2293

2294 distance_max = Float.T(

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

2296

2297 distance_delta = Float.T(

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

2299

2300 short_type = 'B'

2301

2302 provided_schemes = [

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

2304

2305 def get_distance(self, args):

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

2307

2308 def get_surface_distance(self, args):

2309 return args[2]

2310

2311 def get_source_depth(self, args):

2312 return args[1]

2313

2314 def get_receiver_depth(self, args):

2315 return args[0]

2316

2317 def get_source_depths(self):

2318 return self.coords[1]

2319

2320 def _update(self):

2321 self.mins = num.array([

2322 self.receiver_depth_min,

2323 self.source_depth_min,

2324 self.distance_min],

2325 dtype=float)

2326

2327 self.maxs = num.array([

2328 self.receiver_depth_max,

2329 self.source_depth_max,

2330 self.distance_max],

2331 dtype=float)

2332

2333 self.deltas = num.array([

2334 self.receiver_depth_delta,

2335 self.source_depth_delta,

2336 self.distance_delta],

2337 dtype=float)

2338

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

2340 vicinity_eps).astype(int) + 1

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

2342 self.deltat = 1.0/self.sample_rate

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

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

2345 (mi, ma, n) in

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

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

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

2349

2350 def _make_index_functions(self):

2351

2352 amin, bmin, cmin = self.mins

2353 da, db, dc = self.deltas

2354 na, nb, nc = self.ns

2355 ng = self.ncomponents

2356

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

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

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

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

2361 try:

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

2363 (na, nb, nc, ng))

2364 except ValueError:

2365 raise OutOfBounds()

2366

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

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

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

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

2371 try:

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

2373 (na, nb, nc, ng))

2374 except ValueError:

2375 raise OutOfBounds()

2376

2377 def grid_interpolation_coefficients(a, b, c):

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

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

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

2381 return ias, ibs, ics

2382

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

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

2385

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

2387 raise OutOfBounds()

2388

2389 indis = []

2390 weights = []

2391 for ia, va in ias:

2392 iia = ia*nb*nc*ng

2393 for ib, vb in ibs:

2394 iib = ib*nc*ng

2395 for ic, vc in ics:

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

2397 weights.append(va*vb*vc)

2398

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

2400

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

2402

2403 xa = (a - amin) / da

2404 xb = (b - bmin) / db

2405 xc = (c - cmin) / dc

2406

2407 xa_fl = num.floor(xa)

2408 xa_ce = num.ceil(xa)

2409 xb_fl = num.floor(xb)

2410 xb_ce = num.ceil(xb)

2411 xc_fl = num.floor(xc)

2412 xc_ce = num.ceil(xc)

2413 va_fl = 1.0 - (xa - xa_fl)

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

2415 vb_fl = 1.0 - (xb - xb_fl)

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

2417 vc_fl = 1.0 - (xc - xc_fl)

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

2419

2420 ia_fl = xa_fl.astype(int)

2421 ia_ce = xa_ce.astype(int)

2422 ib_fl = xb_fl.astype(int)

2423 ib_ce = xb_ce.astype(int)

2424 ic_fl = xc_fl.astype(int)

2425 ic_ce = xc_ce.astype(int)

2426

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

2428 raise OutOfBounds()

2429

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

2431 raise OutOfBounds()

2432

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

2434 raise OutOfBounds()

2435

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

2437 raise OutOfBounds()

2438

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

2440 raise OutOfBounds()

2441

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

2443 raise OutOfBounds()

2444

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

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

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

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

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

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

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

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

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

2454

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

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

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

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

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

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

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

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

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

2464

2465 return irecords, weights

2466

2467 self._index_function = index_function

2468 self._indices_function = indices_function

2469 self._grid_interpolation_coefficients = grid_interpolation_coefficients

2470 self._vicinity_function = vicinity_function

2471 self._vicinities_function = vicinities_function

2472

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

2474 nc = icomponents.size

2475 dists = source.distances_to(receiver)

2476 n = dists.size

2477 receiver_depths = num.empty(nc)

2478 receiver_depths.fill(receiver.depth)

2479 return (receiver_depths,

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

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

2482 icomponents)

2483

2484 def make_indexing_args1(self, source, receiver):

2485 return (receiver.depth,

2486 source.depth,

2487 source.distance_to(receiver))

2488

2489 @property

2490 def short_extent(self):

2491 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2492 self.receiver_depth_min/km,

2493 self.receiver_depth_max/km,

2494 self.receiver_depth_delta/km,

2495 self.source_depth_min/km,

2496 self.source_depth_max/km,

2497 self.source_depth_delta/km,

2498 self.distance_min/km,

2499 self.distance_max/km,

2500 self.distance_delta/km)

2501

2502 def fix_ttt_holes(self, sptree, mode):

2503 from pyrocko import eikonal_ext, spit

2504

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

2506

2507 delta = self.deltas[-1]

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

2509

2510 nreceivers, nsources, ndistances = self.ns

2511

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

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

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

2515

2516 speeds = self.get_material_property(

2517 0., 0., points,

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

2519 interpolation='multilinear')

2520

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

2522

2523 receiver_times = []

2524 for ireceiver in range(nreceivers):

2525 nodes = num.hstack([

2526 num.full(

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

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

2529 nodes_sr])

2530

2531 times = sptree.interpolate_many(nodes)

2532

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

2534

2535 times = times.reshape(speeds.shape)

2536

2537 try:

2538 eikonal_ext.eikonal_solver_fmm_cartesian(

2539 speeds, times, delta)

2540 except eikonal_ext.EikonalExtError as e:

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

2542 logger.debug(

2543 'Got a warning from eikonal solver '

2544 '- may be ok...')

2545 else:

2546 raise

2547

2548 receiver_times.append(times)

2549

2550 def func(x):

2551 ias, ibs, ics = \

2552 self.grid_interpolation_coefficients(*x)

2553

2554 t = 0

2555 for ia, va in ias:

2556 times = receiver_times[ia]

2557 for ib, vb in ibs:

2558 for ic, vc in ics:

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

2560

2561 return t

2562

2563 return spit.SPTree(

2564 f=func,

2565 ftol=sptree.ftol,

2566 xbounds=sptree.xbounds,

2567 xtols=sptree.xtols)

2568

2569

2570class ConfigTypeC(Config):

2571 '''

2572 No symmetrical constraints, one fixed receiver position.

2573

2574 * Cartesian 3D source volume around a reference point

2575

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

2577 source_east_shift, source_north_shift, component)``

2578 '''

2579

2580 receiver = Receiver.T(

2581 help='Receiver location')

2582

2583 source_origin = Location.T(

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

2585

2586 source_depth_min = Float.T(

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

2588

2589 source_depth_max = Float.T(

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

2591

2592 source_depth_delta = Float.T(

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

2594

2595 source_east_shift_min = Float.T(

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

2597

2598 source_east_shift_max = Float.T(

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

2600

2601 source_east_shift_delta = Float.T(

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

2603

2604 source_north_shift_min = Float.T(

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

2606

2607 source_north_shift_max = Float.T(

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

2609

2610 source_north_shift_delta = Float.T(

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

2612

2613 short_type = 'C'

2614

2615 provided_schemes = ['elastic18']

2616

2617 def get_surface_distance(self, args):

2618 _, source_east_shift, source_north_shift, _ = args

2619 sorig = self.source_origin

2620 sloc = Location(

2621 lat=sorig.lat,

2622 lon=sorig.lon,

2623 north_shift=sorig.north_shift + source_north_shift,

2624 east_shift=sorig.east_shift + source_east_shift)

2625

2626 return self.receiver.distance_to(sloc)

2627

2628 def get_distance(self, args):

2629 sdepth, source_east_shift, source_north_shift, _ = args

2630 sorig = self.source_origin

2631 sloc = Location(

2632 lat=sorig.lat,

2633 lon=sorig.lon,

2634 north_shift=sorig.north_shift + source_north_shift,

2635 east_shift=sorig.east_shift + source_east_shift)

2636

2637 return self.receiver.distance_3d_to(sloc)

2638

2639 def get_source_depth(self, args):

2640 return args[0]

2641

2642 def get_receiver_depth(self, args):

2643 return self.receiver.depth

2644

2645 def get_source_depths(self):

2646 return self.coords[0]

2647

2648 def _update(self):

2649 self.mins = num.array([

2650 self.source_depth_min,

2651 self.source_east_shift_min,

2652 self.source_north_shift_min],

2653 dtype=float)

2654

2655 self.maxs = num.array([

2656 self.source_depth_max,

2657 self.source_east_shift_max,

2658 self.source_north_shift_max],

2659 dtype=float)

2660

2661 self.deltas = num.array([

2662 self.source_depth_delta,

2663 self.source_east_shift_delta,

2664 self.source_north_shift_delta],

2665 dtype=float)

2666

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

2668 vicinity_eps).astype(int) + 1

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

2670 self.deltat = 1.0/self.sample_rate

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

2672

2673 self.coords = tuple(

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

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

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

2677

2678 def _make_index_functions(self):

2679

2680 amin, bmin, cmin = self.mins

2681 da, db, dc = self.deltas

2682 na, nb, nc = self.ns

2683 ng = self.ncomponents

2684

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

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

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

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

2689 try:

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

2691 (na, nb, nc, ng))

2692 except ValueError:

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

2694

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

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

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

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

2699

2700 try:

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

2702 (na, nb, nc, ng))

2703 except ValueError:

2704 raise OutOfBounds()

2705

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

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

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

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

2710

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

2712 raise OutOfBounds()

2713

2714 indis = []

2715 weights = []

2716 for ia, va in ias:

2717 iia = ia*nb*nc*ng

2718 for ib, vb in ibs:

2719 iib = ib*nc*ng

2720 for ic, vc in ics:

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

2722 weights.append(va*vb*vc)

2723

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

2725

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

2727

2728 xa = (a-amin) / da

2729 xb = (b-bmin) / db

2730 xc = (c-cmin) / dc

2731

2732 xa_fl = num.floor(xa)

2733 xa_ce = num.ceil(xa)

2734 xb_fl = num.floor(xb)

2735 xb_ce = num.ceil(xb)

2736 xc_fl = num.floor(xc)

2737 xc_ce = num.ceil(xc)

2738 va_fl = 1.0 - (xa - xa_fl)

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

2740 vb_fl = 1.0 - (xb - xb_fl)

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

2742 vc_fl = 1.0 - (xc - xc_fl)

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

2744

2745 ia_fl = xa_fl.astype(int)

2746 ia_ce = xa_ce.astype(int)

2747 ib_fl = xb_fl.astype(int)

2748 ib_ce = xb_ce.astype(int)

2749 ic_fl = xc_fl.astype(int)

2750 ic_ce = xc_ce.astype(int)

2751

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

2753 raise OutOfBounds()

2754

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

2756 raise OutOfBounds()

2757

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

2759 raise OutOfBounds()

2760

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

2762 raise OutOfBounds()

2763

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

2765 raise OutOfBounds()

2766

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

2768 raise OutOfBounds()

2769

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

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

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

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

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

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

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

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

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

2779

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

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

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

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

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

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

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

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

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

2789

2790 return irecords, weights

2791

2792 self._index_function = index_function

2793 self._indices_function = indices_function

2794 self._vicinity_function = vicinity_function

2795 self._vicinities_function = vicinities_function

2796

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

2798 nc = icomponents.size

2799

2800 dists = source.distances_to(self.source_origin)

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

2802

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

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

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

2806

2807 n = dists.size

2808

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

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

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

2812 icomponents)

2813

2814 def make_indexing_args1(self, source, receiver):

2815 dist = source.distance_to(self.source_origin)

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

2817

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

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

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

2821

2822 return (source_depth,

2823 source_east_shift,

2824 source_north_shift)

2825

2826 @property

2827 def short_extent(self):

2828 return '%g:%g:%g x %g:%g:%g x %g:%g:%g' % (

2829 self.source_depth_min/km,

2830 self.source_depth_max/km,

2831 self.source_depth_delta/km,

2832 self.source_east_shift_min/km,

2833 self.source_east_shift_max/km,

2834 self.source_east_shift_delta/km,

2835 self.source_north_shift_min/km,

2836 self.source_north_shift_max/km,

2837 self.source_north_shift_delta/km)

2838

2839

2840class Weighting(Object):

2841 factor = Float.T(default=1.0)

2842

2843

2844class Taper(Object):

2845 tmin = Timing.T()

2846 tmax = Timing.T()

2847 tfade = Float.T(default=0.0)

2848 shape = StringChoice.T(

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

2850 default='cos',

2851 optional=True)

2852

2853

2854class SimplePattern(SObject):

2855

2856 _pool = {}

2857

2858 def __init__(self, pattern):

2859 self._pattern = pattern

2860 SObject.__init__(self)

2861

2862 def __str__(self):

2863 return self._pattern

2864

2865 @property

2866 def regex(self):

2867 pool = SimplePattern._pool

2868 if self.pattern not in pool:

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

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

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

2872

2873 return pool[self.pattern]

2874

2875 def match(self, s):

2876 return self.regex.match(s)

2877

2878

2879class WaveformType(StringChoice):

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

2881 'amp_spec_dis', 'amp_spec_vel', 'amp_spec_acc',

2882 'envelope_dis', 'envelope_vel', 'envelope_acc']

2883

2884

2885class ChannelSelection(Object):

2886 pattern = SimplePattern.T(optional=True)

2887 min_sample_rate = Float.T(optional=True)

2888 max_sample_rate = Float.T(optional=True)

2889

2890

2891class StationSelection(Object):

2892 includes = SimplePattern.T()

2893 excludes = SimplePattern.T()

2894 distance_min = Float.T(optional=True)

2895 distance_max = Float.T(optional=True)

2896 azimuth_min = Float.T(optional=True)

2897 azimuth_max = Float.T(optional=True)

2898

2899

2900class WaveformSelection(Object):

2901 channel_selection = ChannelSelection.T(optional=True)

2902 station_selection = StationSelection.T(optional=True)

2903 taper = Taper.T()

2904 # filter = FrequencyResponse.T()

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

2906 weighting = Weighting.T(optional=True)

2907 sample_rate = Float.T(optional=True)

2908 gf_store_id = StringID.T(optional=True)

2909

2910

2911def indi12(x, n):

2912 '''

2913 Get linear interpolation index and weight.

2914 '''

2915

2916 r = round(x)

2917 if abs(r - x) < vicinity_eps:

2918 i = int(r)

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

2920 raise OutOfBounds()

2921

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

2923 else:

2924 f = math.floor(x)

2925 i = int(f)

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

2927 raise OutOfBounds()

2928

2929 v = x-f

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

2931

2932

2933def float_or_none(s):

2934 units = {

2935 'k': 1e3,

2936 'M': 1e6,

2937 }

2938

2939 factor = 1.0

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

2941 factor = units[s[-1]]

2942 s = s[:-1]

2943 if not s:

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

2945

2946 if s:

2947 return float(s) * factor

2948 else:

2949 return None

2950

2951

2952class GridSpecError(Exception):

2953 def __init__(self, s):

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

2955

2956

2957def parse_grid_spec(spec):

2958 try:

2959 result = []

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

2961 t = dspec.split('@')

2962 num = start = stop = step = None

2963 if len(t) == 2:

2964 num = int(t[1])

2965 if num <= 0:

2966 raise GridSpecError(spec)

2967

2968 elif len(t) > 2:

2969 raise GridSpecError(spec)

2970

2971 s = t[0]

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

2973 if len(v) == 1:

2974 start = stop = v[0]

2975 if len(v) >= 2:

2976 start, stop = v[0:2]

2977 if len(v) == 3:

2978 step = v[2]

2979

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

2981 raise GridSpecError(spec)

2982

2983 if step == 0.0:

2984 raise GridSpecError(spec)

2985

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

2987

2988 except ValueError:

2989 raise GridSpecError(spec)

2990

2991 return result

2992

2993

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

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

2996 if start is None:

2997 start = ma if swap else mi

2998 if stop is None:

2999 stop = mi if swap else ma

3000 if step is None:

3001 step = -inc if ma < mi else inc

3002 if num is None:

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

3004 raise GridSpecError()

3005

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

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

3008 if abs(stop-stop2) > eps:

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

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

3011 else:

3012 stop = stop2

3013

3014 if start == stop:

3015 num = 1

3016

3017 return start, stop, num

3018

3019

3020def nditer_outer(x):

3021 return num.nditer(

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

3023

3024

3025def nodes(xs):

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

3027 nnodes = num.prod(ns)

3028 ndim = len(xs)

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

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

3031 x = xs[idim]

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

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

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

3035

3036 return nodes

3037

3038

3039def filledi(x, n):

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

3041 a.fill(x)

3042 return a

3043

3044

3045config_type_classes = [ConfigTypeA, ConfigTypeB, ConfigTypeC]

3046

3047discretized_source_classes = [

3048 DiscretizedExplosionSource,

3049 DiscretizedSFSource,

3050 DiscretizedMTSource,

3051 DiscretizedPorePressureSource]

3052

3053

3054__all__ = '''

3055Earthmodel1D

3056StringID

3057ScopeType

3058WaveformType

3059QuantityType

3060NearfieldTermsType

3061Reference

3062Region

3063CircularRegion

3064RectangularRegion

3065PhaseSelect

3066InvalidTimingSpecification

3067Timing

3068TPDef

3069OutOfBounds

3070Location

3071Receiver

3072'''.split() + [

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

3074ComponentScheme

3075component_scheme_to_description

3076component_schemes

3077Config

3078GridSpecError

3079Weighting

3080Taper

3081SimplePattern

3082WaveformType

3083ChannelSelection

3084StationSelection

3085WaveformSelection

3086nditer_outer

3087dump

3088load

3089discretized_source_classes

3090config_type_classes

3091UnavailableScheme

3092InterpolationMethod

3093SeismosizerTrace

3094SeismosizerResult

3095Result

3096StaticResult

3097TimingAttributeError

3098'''.split()