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, U, N, and E.
37 :param source: :py:class:`pyrocko.gf.Location` object
38 :param target: :py:class:`pyrocko.gf.Location` object
39 :param component:
40 string ``'R'``, ``'T'``, ``'Z'``, ``'U'``, ``'N'`` or ``'E'``
41 '''
43 _, bazi = source.azibazi_to(target)
45 azi, dip = {
46 'T': (bazi + 270., 0.),
47 'R': (bazi + 180., 0.),
48 'N': (0., 0.),
49 'E': (90., 0.),
50 'Z': (0., -90.),
51 'U': (0., -90.)}[component]
53 return azi, dip
56class Target(meta.Receiver):
57 '''
58 A seismogram computation request for a single component, including
59 its post-processing parmeters.
60 '''
62 quantity = meta.QuantityType.T(
63 optional=True,
64 help='Measurement quantity type. If not given, it is guessed from the '
65 'channel code. For some common cases, derivatives of the stored '
66 'quantities are supported by using finite difference '
67 'approximations (e.g. displacement to velocity or acceleration). '
68 '4th order central FD schemes are used.')
70 codes = Tuple.T(
71 4, String.T(), default=('', 'STA', '', 'Z'),
72 help='network, station, location and channel codes to be set on '
73 'the response trace.')
75 elevation = Float.T(
76 default=0.0,
77 help='station surface elevation in [m]')
79 store_id = meta.StringID.T(
80 optional=True,
81 help='ID of Green\'s function store to use for the computation. '
82 'If not given, the processor may use a system default.')
84 sample_rate = Float.T(
85 optional=True,
86 help='sample rate to produce. '
87 'If not given the GF store\'s default sample rate is used. '
88 'GF store specific restrictions may apply.')
90 interpolation = meta.InterpolationMethod.T(
91 default='nearest_neighbor',
92 help='Interpolation method between Green\'s functions. Supported are'
93 ' ``nearest_neighbor`` and ``multilinear``')
95 optimization = OptimizationMethod.T(
96 default='enable',
97 optional=True,
98 help='disable/enable optimizations in weight-delay-and-sum operation')
100 tmin = Timestamp.T(
101 optional=True,
102 help='time of first sample to request in [s]. '
103 'If not given, it is determined from the Green\'s functions.')
105 tmax = Timestamp.T(
106 optional=True,
107 help='time of last sample to request in [s]. '
108 'If not given, it is determined from the Green\'s functions.')
110 azimuth = Float.T(
111 optional=True,
112 help='azimuth of sensor component in [deg], clockwise from north. '
113 'If not given, it is guessed from the channel code.')
115 dip = Float.T(
116 optional=True,
117 help='dip of sensor component in [deg], '
118 'measured downward from horizontal. '
119 'If not given, it is guessed from the channel code.')
121 filter = Filter.T(
122 optional=True,
123 help='frequency response filter.')
125 def base_key(self):
126 return (self.store_id, self.sample_rate, self.interpolation,
127 self.optimization,
128 self.tmin, self.tmax,
129 self.elevation, self.depth, self.north_shift, self.east_shift,
130 self.lat, self.lon)
132 def effective_quantity(self):
133 if self.quantity is not None:
134 return self.quantity
136 # guess from channel code
137 cha = self.codes[-1].upper()
138 if len(cha) == 3:
139 # use most common SEED conventions here, however units have to be
140 # guessed, because they are not uniquely defined by the conventions
141 if cha[-2] in 'HL': # high gain, low gain seismometer
142 return 'velocity'
143 if cha[-2] == 'N': # accelerometer
144 return 'acceleration'
145 if cha[-2] == 'D': # hydrophone, barometer, ...
146 return 'pressure'
147 if cha[-2] == 'A': # tiltmeter
148 return 'tilt'
149 elif len(cha) == 2:
150 if cha[-2] == 'U':
151 return 'displacement'
152 if cha[-2] == 'V':
153 return 'velocity'
154 elif len(cha) == 1:
155 if cha in 'NEZU':
156 return 'displacement'
157 if cha == 'P':
158 return 'pressure'
160 raise BadTarget('cannot guess measurement quantity type from channel '
161 'code "%s"' % cha)
163 def receiver(self, store):
164 rec = meta.Receiver(**dict(meta.Receiver.T.inamevals(self)))
165 return rec
167 def component_code(self):
168 cha = self.codes[-1].upper()
169 if cha:
170 return cha[-1]
171 else:
172 return ' '
174 def effective_azimuth(self):
175 if self.azimuth is not None:
176 return self.azimuth
177 elif self.component_code() in 'NEZ':
178 return {'N': 0., 'E': 90., 'Z': 0., 'U': 0.}[self.component_code()]
180 raise BadTarget('cannot determine sensor component azimuth for '
181 '%s.%s.%s.%s' % self.codes)
183 def effective_dip(self):
184 if self.dip is not None:
185 return self.dip
186 elif self.component_code() in 'NEZ':
187 return {
188 'N': 0., 'E': 0., 'Z': -90., 'U': -90.}[self.component_code()]
190 raise BadTarget('cannot determine sensor component dip')
192 def get_sin_cos_factors(self):
193 azi = self.effective_azimuth()
194 dip = self.effective_dip()
195 sa = math.sin(azi*d2r)
196 ca = math.cos(azi*d2r)
197 sd = math.sin(dip*d2r)
198 cd = math.cos(dip*d2r)
199 return sa, ca, sd, cd
201 def get_factor(self):
202 return 1.0
204 def post_process(self, engine, source, tr):
205 return meta.Result(trace=tr)
208class StaticTarget(meta.MultiLocation):
209 '''
210 A computation request for a spatial multi-location target of
211 static/geodetic quantities.
212 '''
213 quantity = meta.QuantityType.T(
214 optional=True,
215 default='displacement',
216 help='Measurement quantity type, for now only `displacement` is'
217 'supported.')
219 interpolation = meta.InterpolationMethod.T(
220 default='nearest_neighbor',
221 help='Interpolation method between Green\'s functions. Supported are'
222 ' ``nearest_neighbor`` and ``multilinear``')
224 tsnapshot = Timestamp.T(
225 optional=True,
226 help='time of the desired snapshot in [s], '
227 'If not given, the first sample is taken. If the desired sample'
228 ' exceeds the length of the Green\'s function store,'
229 ' the last sample is taken.')
231 store_id = meta.StringID.T(
232 optional=True,
233 help='ID of Green\'s function store to use for the computation. '
234 'If not given, the processor may use a system default.')
236 def base_key(self):
237 return (self.store_id,
238 self.coords5.shape,
239 self.quantity,
240 self.tsnapshot,
241 self.interpolation)
243 @property
244 def ntargets(self):
245 '''
246 Number of targets held by instance.
247 '''
248 return self.ncoords
250 def get_targets(self):
251 '''
252 Discretizes the multilocation target into a list of
253 :class:`Target:`
255 :returns: :class:`Target`
256 :rtype: list
257 '''
258 targets = []
259 for i in range(self.ntargets):
260 targets.append(
261 Target(
262 lat=float(self.coords5[i, 0]),
263 lon=float(self.coords5[i, 1]),
264 north_shift=float(self.coords5[i, 2]),
265 east_shift=float(self.coords5[i, 3]),
266 elevation=float(self.coords5[i, 4])))
267 return targets
269 def distance_to(self, source):
270 src_lats = num_full_like(self.lats, fill_value=source.lat)
271 src_lons = num_full_like(self.lons, fill_value=source.lon)
273 target_coords = self.get_latlon()
274 target_lats = target_coords[:, 0]
275 target_lons = target_coords[:, 1]
276 return distance_accurate50m_numpy(
277 src_lats, src_lons, target_lats, target_lons)
279 def post_process(self, engine, source, statics):
280 return meta.StaticResult(result=statics)
283class SatelliteTarget(StaticTarget):
284 '''
285 A computation request for a spatial multi-location target of
286 static/geodetic quantities measured from a satellite instrument.
287 The line of sight angles are provided and projecting
288 post-processing is applied.
289 '''
290 theta = Array.T(
291 shape=(None,),
292 dtype=float,
293 serialize_as='base64-compat',
294 help='Horizontal angle towards satellite\'s line of sight in radians.'
295 '\n\n .. important::\n\n'
296 ' :math:`0` is **east** and'
297 ' :math:`\\frac{\\pi}{2}` is **north**.\n\n')
299 phi = Array.T(
300 shape=(None,),
301 dtype=float,
302 serialize_as='base64-compat',
303 help='Theta is look vector elevation angle towards satellite from'
304 ' horizon in radians. Matrix of theta towards satellite\'s'
305 ' line of sight.'
306 '\n\n .. important::\n\n'
307 ' :math:`-\\frac{\\pi}{2}` is **down** and'
308 ' :math:`\\frac{\\pi}{2}` is **up**.\n\n')
310 def __init__(self, *args, **kwargs):
311 super(SatelliteTarget, self).__init__(*args, **kwargs)
312 self._los_factors = None
314 def get_los_factors(self):
315 if (self.theta.size != self.phi.size != self.lats.size):
316 raise AttributeError('LOS angles inconsistent with provided'
317 ' coordinate shape.')
318 if self._los_factors is None:
319 self._los_factors = num.empty((self.theta.shape[0], 3))
320 self._los_factors[:, 0] = num.sin(self.theta)
321 self._los_factors[:, 1] = num.cos(self.theta) * num.cos(self.phi)
322 self._los_factors[:, 2] = num.cos(self.theta) * num.sin(self.phi)
323 return self._los_factors
325 def post_process(self, engine, source, statics):
326 return meta.SatelliteResult(
327 result=statics,
328 theta=self.theta, phi=self.phi)
331class KiteSceneTarget(SatelliteTarget):
333 shape = Tuple.T(
334 2, Int.T(),
335 optional=False,
336 help='Shape of the displacement vectors.')
338 def __init__(self, scene, **kwargs):
339 size = scene.displacement.size
341 if scene.frame.spacing == 'meter':
342 lats = num_full(size, scene.frame.llLat)
343 lons = num_full(size, scene.frame.llLon)
344 north_shifts = scene.frame.gridN.data.flatten()
345 east_shifts = scene.frame.gridE.data.flatten()
347 elif scene.frame.spacing == 'degree':
348 lats = scene.frame.gridN.data.flatten() + scene.frame.llLat
349 lons = scene.frame.gridE.data.flatten() + scene.frame.llLon
350 north_shifts = num.zeros(size)
351 east_shifts = num.zeros(size)
353 self.scene = scene
355 super(KiteSceneTarget, self).__init__(
356 lats=lats, lons=lons,
357 north_shifts=north_shifts, east_shifts=east_shifts,
358 theta=scene.theta.flatten(),
359 phi=scene.phi.flatten(),
360 shape=scene.shape, **kwargs)
362 def post_process(self, engine, source, statics):
363 res = meta.KiteSceneResult(
364 result=statics,
365 theta=self.theta, phi=self.phi,
366 shape=self.scene.shape)
367 res.config = self.scene.config
368 return res
371class GNSSCampaignTarget(StaticTarget):
373 def post_process(self, engine, source, statics):
374 campaign = gnss.GNSSCampaign()
376 for ista in range(self.ntargets):
377 north = gnss.GNSSComponent(
378 shift=float(statics['displacement.n'][ista]))
379 east = gnss.GNSSComponent(
380 shift=float(statics['displacement.e'][ista]))
381 up = gnss.GNSSComponent(
382 shift=-float(statics['displacement.d'][ista]))
384 coords = self.coords5
385 station = gnss.GNSSStation(
386 lat=float(coords[ista, 0]),
387 lon=float(coords[ista, 1]),
388 east_shift=float(coords[ista, 2]),
389 north_shift=float(coords[ista, 3]),
390 elevation=float(coords[ista, 4]),
391 north=north,
392 east=east,
393 up=up)
395 campaign.add_station(station)
397 return meta.GNSSCampaignResult(result=statics, campaign=campaign)