shifter.py 5.55 KB
Newer Older
1
import logging
Sebastian Heimann's avatar
Sebastian Heimann committed
2
3
4
import hashlib
from StringIO import StringIO
import os.path as op
5
import numpy as num
Sebastian Heimann's avatar
Sebastian Heimann committed
6
7
8
from pyrocko.guts import Object, Float, String
from pyrocko import cake, spit, util
from pyrocko.gf import meta
9
from lassie import geo
Sebastian Heimann's avatar
Sebastian Heimann committed
10
from lassie.common import LassieError, CakeEarthmodel
Sebastian Heimann's avatar
Sebastian Heimann committed
11
12

guts_prefix = 'lassie'
13
logger = logging.getLogger('lassie.config')
Sebastian Heimann's avatar
Sebastian Heimann committed
14
15
16


class Shifter(Object):
Sebastian Heimann's avatar
Sebastian Heimann committed
17
18
19

    def setup(self, config):
        pass
Sebastian Heimann's avatar
Sebastian Heimann committed
20
21
22
23


class VelocityShifter(Shifter):
    velocity = Float.T()
24
    offset = Float.T(default=0.0)
Sebastian Heimann's avatar
Sebastian Heimann committed
25
26
27

    def get_table(self, grid, receivers):
        distances = grid.distances(receivers)
28
        return self.offset + distances / self.velocity
Sebastian Heimann's avatar
Sebastian Heimann committed
29

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
    def get_offset_table(self, grid, receivers, origin):
        distances = grid.distances(receivers)

        x, y = geo.points_coords([origin], system=('ne', grid.lat, grid.lon))
        x = x[0]
        y = y[0]
        z = origin.z

        rx, ry = geo.points_coords(
            receivers, system=('ne', grid.lat, grid.lon))

        rz = num.array([r.z for r in receivers], dtype=num.float)

        distances_origin = num.sqrt((x - rx)**2 + (y - ry)**2 + (z - rz)**2)

        na = num.newaxis
        offset_distances = distances - distances_origin[na, :]

        return offset_distances / self.velocity

Sebastian Heimann's avatar
Sebastian Heimann committed
50
51
52
53
    def get_vmin(self):
        return self.velocity


54
class CakePhaseShifter(Shifter):
Sebastian Heimann's avatar
Sebastian Heimann committed
55
    earthmodel_id = String.T()
56
    timing = meta.Timing.T()
57

Sebastian Heimann's avatar
Sebastian Heimann committed
58
59
60
    def setup(self, config):
        Shifter.setup(self, config)
        self._earthmodels = config.earthmodels
61
62
63
64
65
66
67
68
69
        self._earthmodels.extend(
            [
                CakeEarthmodel(
                    id=fn,
                    earthmodel_1d=cake.load_model(
                        cake.builtin_model_filename(fn)))
                for fn in cake.builtin_models()
            ]
        )
Sebastian Heimann's avatar
Sebastian Heimann committed
70
        self._tabulated_phases = config.tabulated_phases
71
72
73
74

        if not self._tabulated_phases:
            raise LassieError('missing tabulated phases in config')

Sebastian Heimann's avatar
Sebastian Heimann committed
75
76
77
78
79
80
81
82
83
84
85
        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)

Marius Kriegerowski's avatar
Marius Kriegerowski committed
86
    def get_vmin(self):
Sebastian Heimann's avatar
Sebastian Heimann committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
        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
Marius Kriegerowski's avatar
Marius Kriegerowski committed
106

107
108
    def get_table(self, grid, receivers):
        distances = grid.lateral_distances(receivers)
Sebastian Heimann's avatar
Sebastian Heimann committed
109
110
111
112
113
114
115
116
117
118
        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()
119

120
121
        interpolated_tts = {}

Sebastian Heimann's avatar
Sebastian Heimann committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
        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)
161

Sebastian Heimann's avatar
Sebastian Heimann committed
162
            interpolated_tts["stored:"+str(phase_def.id)] = sptree
163

Sebastian Heimann's avatar
Sebastian Heimann committed
164
        arrivals = num.zeros(distances.shape)
165
166
167
168

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

169
        for i_r, r in enumerate(receivers):
Sebastian Heimann's avatar
Sebastian Heimann committed
170
171
172
173
174
            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]
175
            arr = self.timing.evaluate(interpolate, coords)
Sebastian Heimann's avatar
Sebastian Heimann committed
176
            arrivals[:, i_r] = arr
177

178
        return arrivals
179
180


Sebastian Heimann's avatar
Sebastian Heimann committed
181
182
183
__all__ = [
    'Shifter',
    'VelocityShifter',
184
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
185
]