1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import numpy as num
8import math
10from . import meta
11from pyrocko.guts import Timestamp, Tuple, String, Float, Object,\
12 StringChoice, Int
13from pyrocko.guts_array import Array
14from pyrocko.model import gnss
15from pyrocko.orthodrome import distance_accurate50m_numpy
16from pyrocko.util import num_full_like, num_full
18d2r = num.pi / 180.
21class BadTarget(Exception):
22 pass
25class Filter(Object):
26 pass
29class OptimizationMethod(StringChoice):
30 choices = ['enable', 'disable']
33def component_orientation(source, target, component):
34 '''
35 Get component and azimuth for standard components R, T, Z, N, and E.
37 :param source: :py:class:`pyrocko.gf.Location` object
38 :param target: :py:class:`pyrocko.gf.Location` object
39 :param component: string ``'R'``, ``'T'``, ``'Z'``, ``'N'`` or ``'E'``
40 '''
42 _, bazi = source.azibazi_to(target)
44 azi, dip = {
45 'T': (bazi + 270., 0.),
46 'R': (bazi + 180., 0.),
47 'N': (0., 0.),
48 'E': (90., 0.),
49 'Z': (0., -90.)}[component]
51 return azi, dip
54class Target(meta.Receiver):
55 '''
56 A seismogram computation request for a single component, including
57 its post-processing parmeters.
58 '''
60 quantity = meta.QuantityType.T(
61 optional=True,
62 help='Measurement quantity type. If not given, it is guessed from the '
63 'channel code. For some common cases, derivatives of the stored '
64 'quantities are supported by using finite difference '
65 'approximations (e.g. displacement to velocity or acceleration). '
66 '4th order central FD schemes are used.')
68 codes = Tuple.T(
69 4, String.T(), default=('', 'STA', '', 'Z'),
70 help='network, station, location and channel codes to be set on '
71 'the response trace.')
73 elevation = Float.T(
74 default=0.0,
75 help='station surface elevation in [m]')
77 store_id = meta.StringID.T(
78 optional=True,
79 help='ID of Green\'s function store to use for the computation. '
80 'If not given, the processor may use a system default.')
82 sample_rate = Float.T(
83 optional=True,
84 help='sample rate to produce. '
85 'If not given the GF store\'s default sample rate is used. '
86 'GF store specific restrictions may apply.')
88 interpolation = meta.InterpolationMethod.T(
89 default='nearest_neighbor',
90 help='Interpolation method between Green\'s functions. Supported are'
91 ' ``nearest_neighbor`` and ``multilinear``')
93 optimization = OptimizationMethod.T(
94 default='enable',
95 optional=True,
96 help='disable/enable optimizations in weight-delay-and-sum operation')
98 tmin = Timestamp.T(
99 optional=True,
100 help='time of first sample to request in [s]. '
101 'If not given, it is determined from the Green\'s functions.')
103 tmax = Timestamp.T(
104 optional=True,
105 help='time of last sample to request in [s]. '
106 'If not given, it is determined from the Green\'s functions.')
108 azimuth = Float.T(
109 optional=True,
110 help='azimuth of sensor component in [deg], clockwise from north. '
111 'If not given, it is guessed from the channel code.')
113 dip = Float.T(
114 optional=True,
115 help='dip of sensor component in [deg], '
116 'measured downward from horizontal. '
117 'If not given, it is guessed from the channel code.')
119 filter = Filter.T(
120 optional=True,
121 help='frequency response filter.')
123 def base_key(self):
124 return (self.store_id, self.sample_rate, self.interpolation,
125 self.optimization,
126 self.tmin, self.tmax,
127 self.elevation, self.depth, self.north_shift, self.east_shift,
128 self.lat, self.lon)
130 def effective_quantity(self):
131 if self.quantity is not None:
132 return self.quantity
134 # guess from channel code
135 cha = self.codes[-1].upper()
136 if len(cha) == 3:
137 # use most common SEED conventions here, however units have to be
138 # guessed, because they are not uniquely defined by the conventions
139 if cha[-2] in 'HL': # high gain, low gain seismometer
140 return 'velocity'
141 if cha[-2] == 'N': # accelerometer
142 return 'acceleration'
143 if cha[-2] == 'D': # hydrophone, barometer, ...
144 return 'pressure'
145 if cha[-2] == 'A': # tiltmeter
146 return 'tilt'
147 elif len(cha) == 2:
148 if cha[-2] == 'U':
149 return 'displacement'
150 if cha[-2] == 'V':
151 return 'velocity'
152 elif len(cha) == 1:
153 if cha in 'NEZ':
154 return 'displacement'
155 if cha == 'P':
156 return 'pressure'
158 raise BadTarget('cannot guess measurement quantity type from channel '
159 'code "%s"' % cha)
161 def receiver(self, store):
162 rec = meta.Receiver(**dict(meta.Receiver.T.inamevals(self)))
163 return rec
165 def component_code(self):
166 cha = self.codes[-1].upper()
167 if cha:
168 return cha[-1]
169 else:
170 return ' '
172 def effective_azimuth(self):
173 if self.azimuth is not None:
174 return self.azimuth
175 elif self.component_code() in 'NEZ':
176 return {'N': 0., 'E': 90., 'Z': 0.}[self.component_code()]
178 raise BadTarget('cannot determine sensor component azimuth for '
179 '%s.%s.%s.%s' % self.codes)
181 def effective_dip(self):
182 if self.dip is not None:
183 return self.dip
184 elif self.component_code() in 'NEZ':
185 return {'N': 0., 'E': 0., 'Z': -90.}[self.component_code()]
187 raise BadTarget('cannot determine sensor component dip')
189 def get_sin_cos_factors(self):
190 azi = self.effective_azimuth()
191 dip = self.effective_dip()
192 sa = math.sin(azi*d2r)
193 ca = math.cos(azi*d2r)
194 sd = math.sin(dip*d2r)
195 cd = math.cos(dip*d2r)
196 return sa, ca, sd, cd
198 def get_factor(self):
199 return 1.0
201 def post_process(self, engine, source, tr):
202 return meta.Result(trace=tr)
205class StaticTarget(meta.MultiLocation):
206 '''
207 A computation request for a spatial multi-location target of
208 static/geodetic quantities.
209 '''
210 quantity = meta.QuantityType.T(
211 optional=True,
212 default='displacement',
213 help='Measurement quantity type, for now only `displacement` is'
214 'supported.')
216 interpolation = meta.InterpolationMethod.T(
217 default='nearest_neighbor',
218 help='Interpolation method between Green\'s functions. Supported are'
219 ' ``nearest_neighbor`` and ``multilinear``')
221 tsnapshot = Timestamp.T(
222 optional=True,
223 help='time of the desired snapshot in [s], '
224 'If not given, the first sample is taken. If the desired sample'
225 ' exceeds the length of the Green\'s function store,'
226 ' the last sample is taken.')
228 store_id = meta.StringID.T(
229 optional=True,
230 help='ID of Green\'s function store to use for the computation. '
231 'If not given, the processor may use a system default.')
233 def base_key(self):
234 return (self.store_id,
235 self.coords5.shape,
236 self.quantity,
237 self.tsnapshot,
238 self.interpolation)
240 @property
241 def ntargets(self):
242 '''
243 Number of targets held by instance.
244 '''
245 return self.ncoords
247 def get_targets(self):
248 '''
249 Discretizes the multilocation target into a list of
250 :class:`Target:`
252 :returns: :class:`Target`
253 :rtype: list
254 '''
255 targets = []
256 for i in range(self.ntargets):
257 targets.append(
258 Target(
259 lat=float(self.coords5[i, 0]),
260 lon=float(self.coords5[i, 1]),
261 north_shift=float(self.coords5[i, 2]),
262 east_shift=float(self.coords5[i, 3]),
263 elevation=float(self.coords5[i, 4])))
264 return targets
266 def distance_to(self, source):
267 src_lats = num_full_like(self.lats, fill_value=source.lat)
268 src_lons = num_full_like(self.lons, fill_value=source.lon)
270 target_coords = self.get_latlon()
271 target_lats = target_coords[:, 0]
272 target_lons = target_coords[:, 1]
273 return distance_accurate50m_numpy(
274 src_lats, src_lons, target_lats, target_lons)
276 def post_process(self, engine, source, statics):
277 return meta.StaticResult(result=statics)
280class SatelliteTarget(StaticTarget):
281 '''
282 A computation request for a spatial multi-location target of
283 static/geodetic quantities measured from a satellite instrument.
284 The line of sight angles are provided and projecting
285 post-processing is applied.
286 '''
287 theta = Array.T(
288 shape=(None,),
289 dtype=float,
290 serialize_as='base64-compat',
291 help='Horizontal angle towards satellite\'s line of sight in radians.'
292 '\n\n .. important::\n\n'
293 ' :math:`0` is **east** and'
294 ' :math:`\\frac{\\pi}{2}` is **north**.\n\n')
296 phi = Array.T(
297 shape=(None,),
298 dtype=float,
299 serialize_as='base64-compat',
300 help='Theta is look vector elevation angle towards satellite from'
301 ' horizon in radians. Matrix of theta towards satellite\'s'
302 ' line of sight.'
303 '\n\n .. important::\n\n'
304 ' :math:`-\\frac{\\pi}{2}` is **down** and'
305 ' :math:`\\frac{\\pi}{2}` is **up**.\n\n')
307 def __init__(self, *args, **kwargs):
308 super(SatelliteTarget, self).__init__(*args, **kwargs)
309 self._los_factors = None
311 def get_los_factors(self):
312 if (self.theta.size != self.phi.size != self.lats.size):
313 raise AttributeError('LOS angles inconsistent with provided'
314 ' coordinate shape.')
315 if self._los_factors is None:
316 self._los_factors = num.empty((self.theta.shape[0], 3))
317 self._los_factors[:, 0] = num.sin(self.theta)
318 self._los_factors[:, 1] = num.cos(self.theta) * num.cos(self.phi)
319 self._los_factors[:, 2] = num.cos(self.theta) * num.sin(self.phi)
320 return self._los_factors
322 def post_process(self, engine, source, statics):
323 return meta.SatelliteResult(
324 result=statics,
325 theta=self.theta, phi=self.phi)
328class KiteSceneTarget(SatelliteTarget):
330 shape = Tuple.T(
331 2, Int.T(),
332 optional=False,
333 help='Shape of the displacement vectors.')
335 def __init__(self, scene, **kwargs):
336 size = scene.displacement.size
338 if scene.frame.spacing == 'meter':
339 lats = num_full(size, scene.frame.llLat)
340 lons = num_full(size, scene.frame.llLon)
341 north_shifts = scene.frame.gridN.data.flatten()
342 east_shifts = scene.frame.gridE.data.flatten()
344 elif scene.frame.spacing == 'degree':
345 lats = scene.frame.gridN.data.flatten() + scene.frame.llLat
346 lons = scene.frame.gridE.data.flatten() + scene.frame.llLon
347 north_shifts = num.zeros(size)
348 east_shifts = num.zeros(size)
350 self.scene = scene
352 super(KiteSceneTarget, self).__init__(
353 lats=lats, lons=lons,
354 north_shifts=north_shifts, east_shifts=east_shifts,
355 theta=scene.theta.flatten(),
356 phi=scene.phi.flatten(),
357 shape=scene.shape, **kwargs)
359 def post_process(self, engine, source, statics):
360 res = meta.KiteSceneResult(
361 result=statics,
362 theta=self.theta, phi=self.phi,
363 shape=self.scene.shape)
364 res.config = self.scene.config
365 return res
368class GNSSCampaignTarget(StaticTarget):
370 def post_process(self, engine, source, statics):
371 campaign = gnss.GNSSCampaign()
373 for ista in range(self.ntargets):
374 north = gnss.GNSSComponent(
375 shift=float(statics['displacement.n'][ista]))
376 east = gnss.GNSSComponent(
377 shift=float(statics['displacement.e'][ista]))
378 up = gnss.GNSSComponent(
379 shift=-float(statics['displacement.d'][ista]))
381 coords = self.coords5
382 station = gnss.GNSSStation(
383 lat=float(coords[ista, 0]),
384 lon=float(coords[ista, 1]),
385 east_shift=float(coords[ista, 2]),
386 north_shift=float(coords[ista, 3]),
387 elevation=float(coords[ista, 4]),
388 north=north,
389 east=east,
390 up=up)
392 campaign.add_station(station)
394 return meta.GNSSCampaignResult(result=statics, campaign=campaign)