shifter.py 2.79 KB
Newer Older
1
2
3
from pyrocko.guts import Object, Float, String, List
from pyrocko import cake, gf, spit
from pyrocko.gf import Earthmodel1D, meta
4
5
6
import numpy as num
import logging

Sebastian Heimann's avatar
Sebastian Heimann committed
7
8

guts_prefix = 'lassie'
9
logger = logging.getLogger('lassie.config')
Sebastian Heimann's avatar
Sebastian Heimann committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


class Shifter(Object):
    pass


class VelocityShifter(Shifter):
    velocity = Float.T()

    def get_table(self, grid, receivers):
        distances = grid.distances(receivers)
        return distances / self.velocity

    def get_vmin(self):
        return self.velocity


27
class CakePhaseShifter(Shifter):
28
29
    earthmodel_1d = Earthmodel1D.T()
    timing = meta.Timing.T()
30

Marius Kriegerowski's avatar
Marius Kriegerowski committed
31
32
33
    def get_vmin(self):
        return min(filter(lambda x: x!=0., self.earthmodel_1d.profile('vs')))

34
    def get_table(self, grid, receivers):
Marius Kriegerowski's avatar
Marius Kriegerowski committed
35
36
        phases = [p.split(':')[1] for p in self.timing.phase_defs]
        tabulated_phases = [gf.meta.TPDef(id=p, definition=p) for p in phases]
37
38
        distances = grid.lateral_distances(receivers)
        s_depths = grid._coords[2]
39
40
41
        x_bounds = num.array(((grid.zmin, grid.zmax),
                              (min(num.ravel(distances)),
                               max(num.ravel(distances)))))
42

Marius Kriegerowski's avatar
Marius Kriegerowski committed
43
44
        x_tolerance = num.array((grid.dz/10., grid.dx/10.))
        t_tolerance = grid.max_delta()/(self.get_vmin()*20.)
45
        interpolated_tts = {}
Marius Kriegerowski's avatar
Marius Kriegerowski committed
46
        for phase_def in tabulated_phases:
47
48
49
50
51
52
53
54

            def evaluate(args):
                source_depth, x = args
                t = []
                rays = self.earthmodel_1d.arrivals(
                    phases=phase_def.phases,
                    distances=[x*cake.m2d],
                    zstart=source_depth,
Marius Kriegerowski's avatar
Marius Kriegerowski committed
55
                    zstop=0.)
56
57
58
59
60
61
62
63
64

                for ray in rays:
                    t.append(ray.t)

                if t:
                    return min(t)
                else:
                    return None

65
            logger.info('prepare tabulated phases: %s' % phase_def.id)
66
67
            sptree = spit.SPTree(
                f=evaluate,
Marius Kriegerowski's avatar
Marius Kriegerowski committed
68
                ftol=t_tolerance,
69
70
71
                xbounds=x_bounds,
                xtols=x_tolerance)

72
            interpolated_tts["stored:"+str(phase_def.id)] = sptree
73
74

            #sptree.dump(filename='sptree_%s.yaml' % phase_def.id)
75
76

        arrivals = num.zeros((len(distances)*len(s_depths), len(receivers)))
77
78
79
80

        def interpolate(phase_id):
            return interpolated_tts[phase_id].interpolate_many

81
        for i_r, r in enumerate(receivers):
82
            distances = grid.lateral_distances([r])
Marius Kriegerowski's avatar
Marius Kriegerowski committed
83
84
85
            a = num.tile(num.ravel(distances), len(s_depths))
            b = num.repeat(s_depths, len(distances))
            coords = num.vstack((b, a)).T
86
            arr = self.timing.evaluate(interpolate, coords)
Marius Kriegerowski's avatar
Marius Kriegerowski committed
87
            arrivals[:, i_r] = num.ravel(arr)
88

89
        return arrivals
90
91


Sebastian Heimann's avatar
Sebastian Heimann committed
92
93
94
__all__ = [
    'Shifter',
    'VelocityShifter',
95
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
96
]