Browse Source

improved cake shifter

deluxe_memory
Sebastian Heimann 4 years ago
parent
commit
d057c605ab
  1. 13
      src/common.py
  2. 14
      src/config.py
  3. 7
      src/core.py
  4. 11
      src/grid.py
  5. 4
      src/ifc.py
  6. 149
      src/shifter.py

13
src/common.py

@ -1,4 +1,6 @@
from collections import defaultdict
from pyrocko.guts import Object, String
from pyrocko.gf import Earthmodel1D
guts_prefix = 'lassie'
@ -7,6 +9,14 @@ class LassieError(Exception):
pass
class Earthmodel(Object):
id = String.T()
class CakeEarthmodel(Earthmodel):
earthmodel_1d = Earthmodel1D.T()
def grouped_by(l, key):
d = defaultdict(list)
for x in l:
@ -14,6 +24,9 @@ def grouped_by(l, key):
return d
__all__ = [
'LassieError',
'Earthmodel',
'CakeEarthmodel',
]

14
src/config.py

@ -3,8 +3,10 @@ import logging
from pyrocko.guts import Object, String, Float, Timestamp, List, Bool, load
from pyrocko import model
from pyrocko.fdsn import station as fs
from pyrocko.gf import TPDef
from lassie import receiver, ifc, grid, geo
from lassie.common import Earthmodel
guts_prefix = 'lassie'
@ -96,6 +98,18 @@ class Config(Object):
help='fill incomplete trace time windows with zeros '
'(and let edge effects ruin your day)')
earthmodels = List.T(
Earthmodel.T(),
help='list of earthmodels usable in shifters')
tabulated_phases = List.T(
TPDef.T(),
help='list of tabulated phase definitions usable shifters')
cache_path = String.T(
default='lassie_cache',
help='directory where lassie stores tabulated phases etc.')
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._receivers = None

7
src/core.py

@ -125,6 +125,7 @@ def scan(
ifcs = config.image_function_contributions
for ifc in ifcs:
ifc.setup(config)
ifc.prescan(p)
shift_tables = []
@ -286,10 +287,16 @@ def scan(
shift_table = shift_tables[iifc]
ok = num.isfinite(shift_table[:, istations_selected])
bad = num.logical_not(ok)
shifts = -num.round(
shift_table[:, istations_selected] /
deltat_cf).astype(num.int32)
weights[bad] = 0.0
shifts[bad] = num.max(shifts[ok])
pdata.append((list(trs_selected), shift_table, ifc))
iwmin = int(round((wmin-tpeaksearch-t0) / deltat_cf))

11
src/grid.py

@ -58,6 +58,11 @@ class Carthesian3DGrid(Grid):
return self._coords
def depths(self):
nx, ny, _ = self._shape()
_, _, z = self._get_coords()
return num.repeat(z, nx*ny)
def index_to_location(self, i):
nx, ny, nz = self._shape()
iz, iy, ix = num.unravel_index(i, (nz, ny, nx))
@ -68,7 +73,7 @@ class Carthesian3DGrid(Grid):
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
rx, ry = receiver.receivers_coords(
rx, ry = geo.points_coords(
receivers, system=('ne', self.lat, self.lon))
na = num.newaxis
@ -77,9 +82,9 @@ class Carthesian3DGrid(Grid):
distances = num.sqrt(
(x[na, na, :, na] - rx[na, na, na, :])**2 +
(y[na, :, na, na] - ry[na, na, na, :])**2).reshape((nx*ny, nr))
(y[na, :, na, na] - ry[na, na, na, :])**2).reshape((nx*ny*nr))
return distances
return num.tile(distances, nz).reshape((nx*ny*nz, nr))
def distances(self, receivers):
nx, ny, nz = self._shape()

4
src/ifc.py

@ -26,6 +26,10 @@ class IFC(Object):
Object.__init__(self, *args, **kwargs)
self.shifter.t_tolerance = 1./(self.fmax * 2.)
def setup(self, config):
if self.shifter:
self.shifter.setup(config)
def get_table(self, grid, receivers):
return self.shifter.get_table(grid, receivers)

149
src/shifter.py

@ -1,16 +1,22 @@
import logging
import hashlib
from StringIO import StringIO
import os.path as op
import numpy as num
from pyrocko.guts import Object, Float, String, List
from pyrocko import cake, gf, spit
from pyrocko.gf import Earthmodel1D, meta
from pyrocko.guts import Object, Float, String
from pyrocko import cake, spit, util
from pyrocko.gf import meta
from lassie import geo
from lassie.common import LassieError, CakeEarthmodel
guts_prefix = 'lassie'
logger = logging.getLogger('lassie.config')
class Shifter(Object):
pass
def setup(self, config):
pass
class VelocityShifter(Shifter):
@ -46,66 +52,115 @@ class VelocityShifter(Shifter):
class CakePhaseShifter(Shifter):
earthmodel_1d = Earthmodel1D.T()
earthmodel_id = String.T()
timing = meta.Timing.T()
def setup(self, config):
Shifter.setup(self, config)
self._earthmodels = config.earthmodels
self._tabulated_phases = config.tabulated_phases
self._cache_path = config.cache_path
def get_earthmodel(self):
for earthmodel in self._earthmodels:
if isinstance(earthmodel, CakeEarthmodel):
if earthmodel.id == self.earthmodel_id:
return earthmodel.earthmodel_1d
raise LassieError(
'no cake earthmodel with id "%s" found' % self.earthmodel_id)
def get_vmin(self):
return min(filter(lambda x: x!=0., self.earthmodel_1d.profile('vs')))
vs = self.get_earthmodel().profile('vs')
return num.min(vs[vs != 0.0])
def ttt_path(self, ehash):
return op.join(self._cache_path, ehash + '.spit')
def ttt_hash(self, earthmodel, phases, x_bounds, x_tolerance, t_tolerance):
f = StringIO()
earthmodel.profile('vp').dump(f)
earthmodel.profile('vs').dump(f)
earthmodel.profile('rho').dump(f)
f.write(','.join(phase.definition() for phase in phases))
x_bounds.dump(f)
x_tolerance.dump(f)
f.write('%f' % t_tolerance)
s = f.getvalue()
h = hashlib.sha1(s).hexdigest()
f.close()
return h
def get_table(self, grid, receivers):
phases = [p.split(':')[1] for p in self.timing.phase_defs]
tabulated_phases = [gf.meta.TPDef(id=p, definition=p) for p in phases]
distances = grid.lateral_distances(receivers)
s_depths = grid._coords[2]
x_bounds = num.array(((grid.zmin, grid.zmax),
(min(num.ravel(distances)),
max(num.ravel(distances)))))
r_depths = num.array([r.z for r in receivers], dtype=num.float)
s_depths = grid.depths()
x_bounds = num.array(
[[num.min(r_depths), num.max(r_depths)],
[num.min(s_depths), num.max(s_depths)],
[num.min(distances), num.max(distances)]], dtype=num.float)
x_tolerance = num.array((grid.dz/2., grid.dz/2., grid.dx/2.))
t_tolerance = grid.max_delta()/(self.get_vmin()*5.)
earthmodel = self.get_earthmodel()
x_tolerance = num.array((grid.dz/10., grid.dx/10.))
t_tolerance = grid.max_delta()/(self.get_vmin()*20.)
interpolated_tts = {}
for phase_def in tabulated_phases:
def evaluate(args):
source_depth, x = args
t = []
rays = self.earthmodel_1d.arrivals(
phases=phase_def.phases,
distances=[x*cake.m2d],
zstart=source_depth,
zstop=0.)
for ray in rays:
t.append(ray.t)
if t:
return min(t)
else:
return None
logger.info('prepare tabulated phases: %s' % phase_def.id)
sptree = spit.SPTree(
f=evaluate,
ftol=t_tolerance,
xbounds=x_bounds,
xtols=x_tolerance)
interpolated_tts["stored:"+str(phase_def.id)] = sptree
for phase_def in self._tabulated_phases:
ttt_hash = self.ttt_hash(
earthmodel, phase_def.phases, x_bounds, x_tolerance,
t_tolerance)
fpath = self.ttt_path(ttt_hash)
if not op.exists(fpath):
def evaluate(args):
receiver_depth, source_depth, x = args
t = []
rays = earthmodel.arrivals(
phases=phase_def.phases,
distances=[x*cake.m2d],
zstart=source_depth,
zstop=receiver_depth)
for ray in rays:
t.append(ray.t)
if t:
return min(t)
else:
return None
logger.info(
'prepare tabulated phases: %s [%s]' % (
phase_def.id, ttt_hash))
sptree = spit.SPTree(
f=evaluate,
ftol=t_tolerance,
xbounds=x_bounds,
xtols=x_tolerance)
util.ensuredirs(fpath)
sptree.dump(filename=fpath)
else:
sptree = spit.SPTree(filename=fpath)
#sptree.dump(filename='sptree_%s.yaml' % phase_def.id)
interpolated_tts["stored:"+str(phase_def.id)] = sptree
arrivals = num.zeros((len(distances)*len(s_depths), len(receivers)))
arrivals = num.zeros(distances.shape)
def interpolate(phase_id):
return interpolated_tts[phase_id].interpolate_many
for i_r, r in enumerate(receivers):
distances = grid.lateral_distances([r])
a = num.tile(num.ravel(distances), len(s_depths))
b = num.repeat(s_depths, len(distances))
coords = num.vstack((b, a)).T
r_depths = num.zeros(distances.shape[0]) + r.z
coords = num.zeros((distances.shape[0], 3))
coords[:, 0] = r_depths
coords[:, 1] = s_depths
coords[:, 2] = distances[:, i_r]
arr = self.timing.evaluate(interpolate, coords)
arrivals[:, i_r] = num.ravel(arr)
arrivals[:, i_r] = arr
return arrivals

Loading…
Cancel
Save