shifter.py 5.82 KB
Newer Older
1
import logging
Sebastian Heimann's avatar
Sebastian Heimann committed
2
import hashlib
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
3
from io import BytesIO
Sebastian Heimann's avatar
Sebastian Heimann committed
4
import os.path as op
5
import numpy as num
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
6
from pyrocko.guts import Object, Float, String   # noqa
Sebastian Heimann's avatar
Sebastian Heimann committed
7
8
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
58
    offset = Float.T(default=0.0)
    factor = Float.T(default=1.0)
59

Sebastian Heimann's avatar
Sebastian Heimann committed
60
61
62
    def setup(self, config):
        Shifter.setup(self, config)
        self._earthmodels = config.earthmodels
63
64
65
66
67
68
69
70
71
        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
72
        self._tabulated_phases = config.tabulated_phases
73
74
75
76

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

Sebastian Heimann's avatar
Sebastian Heimann committed
77
        self._cache_path = config.expand_path(config.cache_path)
Sebastian Heimann's avatar
Sebastian Heimann committed
78
79
80
81
82
83
84
85
86
87

    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
88
    def get_vmin(self):
Sebastian Heimann's avatar
Sebastian Heimann committed
89
        vs = self.get_earthmodel().profile('vs')
90
91
92
93
        vp = self.get_earthmodel().profile('vp')
        v = num.concatenate((vs, vp))
        vmin = num.min(v[v != 0.0])
        return vmin
Sebastian Heimann's avatar
Sebastian Heimann committed
94
95
96
97
98

    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):
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
99
100
        f = BytesIO()
        earthmodel.profile('z').dump(f)
Sebastian Heimann's avatar
Sebastian Heimann committed
101
102
103
        earthmodel.profile('vp').dump(f)
        earthmodel.profile('vs').dump(f)
        earthmodel.profile('rho').dump(f)
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
104
105

        f.write(b','.join(phase.definition().encode() for phase in phases))
Sebastian Heimann's avatar
Sebastian Heimann committed
106
107
        x_bounds.dump(f)
        x_tolerance.dump(f)
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
108
        f.write(str(t_tolerance).encode())
Sebastian Heimann's avatar
Sebastian Heimann committed
109
110
111
112
        s = f.getvalue()
        h = hashlib.sha1(s).hexdigest()
        f.close()
        return h
Marius Kriegerowski's avatar
Marius Kriegerowski committed
113

114
115
    def get_table(self, grid, receivers):
        distances = grid.lateral_distances(receivers)
Sebastian Heimann's avatar
Sebastian Heimann committed
116
117
118
119
120
121
122
123
124
125
        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()
126

127
128
        interpolated_tts = {}

Sebastian Heimann's avatar
Sebastian Heimann committed
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
161
162
163
164
165
166
167
        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)
168

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

Sebastian Heimann's avatar
Sebastian Heimann committed
171
        arrivals = num.zeros(distances.shape)
172
173
174
175

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

176
        for i_r, r in enumerate(receivers):
Sebastian Heimann's avatar
Sebastian Heimann committed
177
178
179
180
181
            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]
182
            arr = self.timing.evaluate(interpolate, coords)
Sebastian Heimann's avatar
Sebastian Heimann committed
183
            arrivals[:, i_r] = arr
184

185
        return arrivals * self.factor + self.offset
186
187


Sebastian Heimann's avatar
Sebastian Heimann committed
188
189
190
__all__ = [
    'Shifter',
    'VelocityShifter',
191
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
192
]