shifter.py 2.71 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
30
    earthmodel_1d = Earthmodel1D.T()
    tabulated_phases = List.T(gf.meta.TPDef.T())
    timing = meta.Timing.T()
31
32
33
34

    def get_table(self, grid, receivers):
        distances = grid.lateral_distances(receivers)
        s_depths = grid._coords[2]
35
36
37
        x_bounds = num.array(((grid.zmin, grid.zmax),
                              (min(num.ravel(distances)),
                               max(num.ravel(distances)))))
38
39

        x_tolerance = num.array((grid.dz/2., grid.dx/2.))
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        interpolated_tts = {}
        for phase_def in self.tabulated_phases:
            v_horizontal = phase_def.horizontal_velocities

            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.)    # TODO: receiver depths

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

                for v in v_horizontal:
                    t.append(x/(v*1000.))
                if t:
                    return min(t)
                else:
                    return None

63
            logger.info('prepare tabulated phases: %s' % phase_def.id)
64
65
66
67
68
69
            sptree = spit.SPTree(
                f=evaluate,
                ftol=self.t_tolerance,
                xbounds=x_bounds,
                xtols=x_tolerance)

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

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

        n_x = len(distances)

        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
83
84
85
86
            distances = grid.lateral_distances([r])
            coords = num.array((num.tile(s_depths, len(distances)),
                                num.repeat(distances, len(s_depths)))).T
            arr = self.timing.evaluate(interpolate, coords)
            arrivals[:, i_r] = arr
87

88
        return arrivals
89
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
]