shifter.py 5.72 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
90
91
92
93
94
95
        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):
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
96
97
        f = BytesIO()
        earthmodel.profile('z').dump(f)
Sebastian Heimann's avatar
Sebastian Heimann committed
98
99
100
        earthmodel.profile('vp').dump(f)
        earthmodel.profile('vs').dump(f)
        earthmodel.profile('rho').dump(f)
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
101
102

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

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

124
125
        interpolated_tts = {}

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

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

Sebastian Heimann's avatar
Sebastian Heimann committed
168
        arrivals = num.zeros(distances.shape)
169
170
171
172

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

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

182
        return arrivals * self.factor + self.offset
183
184


Sebastian Heimann's avatar
Sebastian Heimann committed
185
186
187
__all__ = [
    'Shifter',
    'VelocityShifter',
188
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
189
]