shifter.py 5.16 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
61
62
63
64
65
66
67
68
69
70
71
72
    def setup(self, config):
        Shifter.setup(self, config)
        self._earthmodels = config.earthmodels
        self._tabulated_phases = config.tabulated_phases
        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
73
    def get_vmin(self):
Sebastian Heimann's avatar
Sebastian Heimann committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
        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
93

94
95
    def get_table(self, grid, receivers):
        distances = grid.lateral_distances(receivers)
Sebastian Heimann's avatar
Sebastian Heimann committed
96
97
98
99
100
101
102
103
104
105
        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()
106

107
108
        interpolated_tts = {}

Sebastian Heimann's avatar
Sebastian Heimann committed
109
110
111
112
113
114
115
116
117
118
119
120
121
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
        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)
148

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

Sebastian Heimann's avatar
Sebastian Heimann committed
151
        arrivals = num.zeros(distances.shape)
152
153
154
155

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

156
        for i_r, r in enumerate(receivers):
Sebastian Heimann's avatar
Sebastian Heimann committed
157
158
159
160
161
            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]
162
            arr = self.timing.evaluate(interpolate, coords)
Sebastian Heimann's avatar
Sebastian Heimann committed
163
            arrivals[:, i_r] = arr
164

165
        return arrivals
166
167


Sebastian Heimann's avatar
Sebastian Heimann committed
168
169
170
__all__ = [
    'Shifter',
    'VelocityShifter',
171
    'CakePhaseShifter',
Sebastian Heimann's avatar
Sebastian Heimann committed
172
]