shifter.py 3.13 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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
        x_bounds = num.array(((grid.zmin, grid.zmax),
                              (min(num.ravel(distances)),
                               max(num.ravel(distances)))))
        x_tolerance = num.array((grid.dx, grid.dz))
        interpolated_tts = {}
        for phase_def in self.tabulated_phases:
            v_horizontal = phase_def.horizontal_velocities

            def evaluate(args):
                '''Calculate arrival using source and receiver location defined by *args*.
                To be evaluated by the SPTree instance.'''
                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

            sptree = spit.SPTree(
                f=evaluate,
                ftol=self.t_tolerance,
                xbounds=x_bounds,
                xtols=x_tolerance)

            interpolated_tts[phase_def.id] = sptree

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

        n_x = len(distances)

        arrivals = num.zeros((len(distances)*len(s_depths), len(receivers)))
        for i_r, r in enumerate(receivers):
78
            logger.info('interpolate arrivals at receiver %s / %s' % (i_r+1, len(receivers)))
79
80
81
            for i_z, z in enumerate(s_depths):
                distances = grid.lateral_distances([r])

82
83
84
85
86
87
88
                #coords = num.array((distances*cake.m2d,
                                    #num.repeat(z, len(distances))))
                #arr = interpolated_tts[phase_def.id].interpolate_many(coords)
                arr = num.zeros(len(distances))
                for i_d, d in enumerate(distances):
                    arr[i_d] = self.timing.evaluate(lambda x:
                                                    interpolated_tts[x], (z, d))
89
                assert len(arr) == len(distances)
90
                arrivals[i_z * n_x: (i_z+1) * n_x, i_r] = arr
91

92
        return arrivals
93
94


Sebastian Heimann's avatar
Sebastian Heimann committed
95
96
97
__all__ = [
    'Shifter',
    'VelocityShifter',
98
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
99
]