1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import numpy as num
8import logging
9from os import path as op
11from datetime import datetime
12from pyrocko import gf, util
13from pyrocko import orthodrome as od
14from pyrocko.guts import Float, Timestamp, Tuple, StringChoice, Bool, Object,\
15 String
16from pyrocko.util import num_full
18from .base import TargetGenerator, NoiseGenerator
19from ..base import get_gsshg
21DEFAULT_STORE_ID = 'ak135_static'
23km = 1e3
24d2r = num.pi/180.
26logger = logging.getLogger('pyrocko.scenario.targets.insar')
27guts_prefix = 'pf.scenario'
30class ScenePatch(Object):
31 lat_center = Float.T(
32 help='Center latitude anchor.')
33 lon_center = Float.T(
34 help='center longitude anchor.')
35 time_master = Timestamp.T(
36 help='Timestamp of the master.')
37 time_slave = Timestamp.T(
38 help='Timestamp of the slave.')
39 inclination = Float.T(
40 help='Orbital inclination towards the equatorial plane [deg].')
41 apogee = Float.T(
42 help='Apogee of the satellite in [m].')
43 swath_width = Float.T(
44 help='Swath width in [m].')
45 track_length = Float.T(
46 help='Track length in [m].')
47 incident_angle = Float.T(
48 help='Ground incident angle in [deg].')
49 resolution = Tuple.T(
50 help='Resolution of raster in east x north [px].')
51 orbital_node = StringChoice.T(
52 ['Ascending', 'Descending'],
53 help='Orbit heading.')
54 mask_water = Bool.T(
55 default=True,
56 help='Mask water bodies.')
58 class SatelliteGeneratorTarget(gf.SatelliteTarget):
60 def __init__(self, scene_patch, *args, **kwargs):
61 gf.SatelliteTarget.__init__(self, *args, **kwargs)
63 self.scene_patch = scene_patch
65 def post_process(self, *args, **kwargs):
66 resp = gf.SatelliteTarget.post_process(self, *args, **kwargs)
68 from kite import Scene
69 from kite.scene import SceneConfig, FrameConfig, Meta
71 patch = self.scene_patch
73 grid, _ = patch.get_grid()
75 displacement = num.empty_like(grid)
76 displacement.fill(num.nan)
77 displacement[patch.get_mask()] = resp.result['displacement.los']
79 theta, phi = patch.get_incident_angles()
81 llLat, llLon = patch.get_ll_anchor()
82 urLat, urLon = patch.get_ur_anchor()
83 dLon = num.abs(llLon - urLon) / patch.resolution[0]
84 dLat = num.abs(llLat - urLat) / patch.resolution[1]
86 scene_config = SceneConfig(
87 meta=Meta(
88 scene_title='Pyrocko Scenario Generator - {orbit} ({time})'
89 .format(orbit=self.scene_patch.orbital_node,
90 time=datetime.now()),
91 orbital_node=patch.orbital_node,
92 scene_id='pyrocko_scenario_%s'
93 % self.scene_patch.orbital_node,
94 satellite_name='Sentinel-1 (pyrocko-scenario)'),
95 frame=FrameConfig(
96 llLat=float(llLat),
97 llLon=float(llLon),
98 dN=float(dLat),
99 dE=float(dLon),
100 spacing='degree'))
102 scene = Scene(
103 displacement=displacement,
104 theta=theta,
105 phi=phi,
106 config=scene_config)
108 resp.scene = scene
110 return resp
112 def __init__(self, *args, **kwargs):
113 Object.__init__(self, *args, **kwargs)
114 self._mask_water = None
116 @property
117 def width(self):
118 track_shift = num.abs(
119 num.cos(self.inclination*d2r) * self.track_length)
121 return self.swath_width + track_shift
123 def get_ll_anchor(self):
124 return od.ne_to_latlon(self.lat_center, self.lon_center,
125 -self.track_length/2, -self.width/2)
127 def get_ur_anchor(self):
128 return od.ne_to_latlon(self.lat_center, self.lon_center,
129 self.track_length/2, self.width/2)
131 def get_ul_anchor(self):
132 return od.ne_to_latlon(self.lat_center, self.lon_center,
133 self.track_length/2, -self.width/2)
135 def get_lr_anchor(self):
136 return od.ne_to_latlon(self.lat_center, self.lon_center,
137 -self.track_length/2, self.width/2)
139 def get_corner_coordinates(self):
140 inc = self.inclination
142 llLat, llLon = self.get_ll_anchor()
143 urLat, urLon = self.get_ur_anchor()
145 if self.orbital_node == 'Ascending':
147 ulLat, ulLon = od.ne_to_latlon(
148 self.lat_center, self.lon_center,
149 self.track_length/2,
150 -num.tan(inc*d2r) * self.width/2)
151 lrLat, lrLon = od.ne_to_latlon(
152 self.lat_center, self.lon_center,
153 -self.track_length/2,
154 num.tan(inc*d2r) * self.width/2)
156 elif self.orbital_node == 'Descending':
157 urLat, urLon = od.ne_to_latlon(
158 self.lat_center, self.lon_center,
159 self.track_length/2,
160 num.tan(inc*d2r) * self.width/2)
161 llLat, llLon = od.ne_to_latlon(
162 self.lat_center, self.lon_center,
163 -self.track_length/2,
164 -num.tan(inc*d2r) * self.width/2)
166 return ((llLat, llLon), (ulLat, ulLon),
167 (urLat, urLon), (lrLat, lrLon))
169 def get_grid(self):
170 '''
171 Return relative positions of scatterer.
173 :param track: Acquisition track, from `'asc'` or `'dsc'`.
174 :type track: string
175 '''
176 easts = num.linspace(0, self.width,
177 self.resolution[0])
178 norths = num.linspace(0, self.track_length,
179 self.resolution[1])
181 return num.meshgrid(easts, norths)
183 def get_mask_track(self):
184 east_shifts, north_shifts = self.get_grid()
185 norths = north_shifts[:, 0]
186 track = num.abs(num.cos(self.inclination*d2r)) * norths
188 track_mask = num.logical_and(
189 east_shifts > track[:, num.newaxis],
190 east_shifts < (track + self.swath_width)[:, num.newaxis])
192 if self.orbital_node == 'Ascending':
193 track_mask = num.fliplr(track_mask)
195 return track_mask
197 def get_mask_water(self):
198 if self._mask_water is None:
199 east_shifts, north_shifts = self.get_grid()
201 east_shifts -= east_shifts[0, -1]/2
202 north_shifts -= north_shifts[-1, -1]/2
204 latlon = od.ne_to_latlon(self.lat_center, self.lon_center,
205 north_shifts.ravel(), east_shifts.ravel())
206 points = num.array(latlon).T
207 self._mask_water = get_gsshg().get_land_mask(points)\
208 .reshape(*east_shifts.shape)
210 return self._mask_water
212 def get_mask(self):
213 mask_track = self.get_mask_track()
214 if self.mask_water:
215 mask_water = self.get_mask_water()
216 return num.logical_and(mask_track, mask_water)
217 return mask_track
219 def get_incident_angles(self):
220 # theta: elevation angle towards satellite from horizon in radians.
221 # phi: Horizontal angle towards satellite' :abbr:`line of sight (LOS)`
222 # in [rad] from East.
223 east_shifts, _ = self.get_grid()
225 phi = num.empty_like(east_shifts)
226 theta = num.empty_like(east_shifts)
228 east_shifts += num.tan(self.incident_angle*d2r) * self.apogee
229 theta = num.arctan(east_shifts/self.apogee)
231 if self.orbital_node == 'Ascending':
232 phi.fill(self.inclination*d2r + num.pi/2)
233 elif self.orbital_node == 'Descending':
234 phi.fill(2*num.pi-(self.inclination*d2r + 3/2*num.pi))
235 theta = num.fliplr(theta)
236 else:
237 raise AttributeError(
238 'Orbital node %s not defined!' % self.orbital_node)
240 theta[~self.get_mask()] = num.nan
241 phi[~self.get_mask()] = num.nan
243 return theta, phi
245 def get_target(self):
246 gE, gN = self.get_grid()
247 mask = self.get_mask()
249 east_shifts = gE[mask].ravel()
250 north_shifts = gN[mask].ravel()
251 llLat, llLon = self.get_ll_anchor()
253 ncoords = east_shifts.size
255 theta, phi = self.get_incident_angles()
257 theta = theta[mask].ravel()
258 phi = phi[mask].ravel()
260 if ncoords == 0:
261 logger.warning('InSAR taget has no valid points,'
262 ' maybe it\'s all water?')
264 return self.SatelliteGeneratorTarget(
265 scene_patch=self,
266 lats=num_full(ncoords, fill_value=llLat),
267 lons=num_full(ncoords, fill_value=llLon),
268 east_shifts=east_shifts,
269 north_shifts=north_shifts,
270 theta=theta,
271 phi=phi)
274class AtmosphericNoiseGenerator(NoiseGenerator):
276 amplitude = Float.T(
277 default=1.,
278 help='Amplitude of the atmospheric noise.')
280 beta = [5./3, 8./3, 2./3]
281 regimes = [.15, .99, 1.]
283 def get_noise(self, scene):
284 nE = scene.frame.rows
285 nN = scene.frame.cols
287 if (nE+nN) % 2 != 0:
288 raise ArithmeticError('Dimensions of synthetic scene must '
289 'both be even!')
291 dE = scene.frame.dE
292 dN = scene.frame.dN
294 rfield = num.random.rand(nE, nN)
295 spec = num.fft.fft2(rfield)
297 kE = num.fft.fftfreq(nE, dE)
298 kN = num.fft.fftfreq(nN, dN)
299 k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)
301 regimes = num.array(self.regimes)
303 k0 = 0.
304 k1 = regimes[0] * k_rad.max()
305 k2 = regimes[1] * k_rad.max()
307 r0 = num.logical_and(k_rad > k0, k_rad < k1)
308 r1 = num.logical_and(k_rad >= k1, k_rad < k2)
309 r2 = k_rad >= k2
311 beta = num.array(self.beta)
312 # From Hanssen (2001)
313 # beta+1 is used as beta, since, the power exponent
314 # is defined for a 1D slice of the 2D spectrum:
315 # austin94: "Adler, 1981, shows that the surface profile
316 # created by the intersection of a plane and a
317 # 2-D fractal surface is itself fractal with
318 # a fractal dimension equal to that of the 2D
319 # surface decreased by one."
320 beta += 1.
321 # From Hanssen (2001)
322 # The power beta/2 is used because the power spectral
323 # density is proportional to the amplitude squared
324 # Here we work with the amplitude, instead of the power
325 # so we should take sqrt( k.^beta) = k.^(beta/2) RH
326 # beta /= 2.
328 amp = num.zeros_like(k_rad)
329 amp[r0] = k_rad[r0] ** -beta[0]
330 amp[r0] /= amp[r0].max()
332 amp[r1] = k_rad[r1] ** -beta[1]
333 amp[r1] /= amp[r1].max() / amp[r0].min()
335 amp[r2] = k_rad[r2] ** -beta[2]
336 amp[r2] /= amp[r2].max() / amp[r1].min()
338 amp[k_rad == 0.] = amp.max()
340 spec *= self.amplitude * num.sqrt(amp)
341 noise = num.abs(num.fft.ifft2(spec))
342 noise -= num.mean(noise)
344 return noise
347class InSARGenerator(TargetGenerator):
348 # https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath
349 store_id = String.T(
350 default=DEFAULT_STORE_ID,
351 help='Store ID for these stations.')
353 inclination = Float.T(
354 default=98.2,
355 help='Inclination of the satellite orbit towards equatorial plane'
356 ' in [deg]. Defaults to Sentinel-1 (98.1 deg).')
357 apogee = Float.T(
358 default=693.*km,
359 help='Apogee of the satellite in [m]. '
360 'Defaults to Sentinel-1 (693 km).')
361 swath_width = Float.T(
362 default=250*km,
363 help='Swath width in [m]. '
364 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).'
365 ' (IW; 250 km).')
366 track_length = Float.T(
367 default=250*km,
368 help='Track length in [m]. Defaults to 200 km.')
369 incident_angle = Float.T(
370 default=29.1,
371 help='Near range incident angle in [deg]. Defaults to 29.1 deg;'
372 ' Sentinel IW mode (29.1 - 46.0 deg).')
373 resolution = Tuple.T(
374 default=(250, 250),
375 help='Resolution of raster in east x north [px].')
376 mask_water = Bool.T(
377 default=True,
378 help='Mask out water bodies.')
379 noise_generator = NoiseGenerator.T(
380 default=AtmosphericNoiseGenerator.D(),
381 help='Add atmospheric noise model after Hansen, 2001.')
383 def get_scene_patches(self):
384 lat_center, lon_center = self.get_center_latlon()
386 scene_patches = []
387 for direction in ('Ascending', 'Descending'):
388 patch = ScenePatch(
389 lat_center=lat_center,
390 lon_center=lon_center,
391 time_master=0,
392 time_slave=0,
393 inclination=self.inclination,
394 apogee=self.apogee,
395 swath_width=self.swath_width,
396 track_length=self.track_length,
397 incident_angle=self.incident_angle,
398 resolution=self.resolution,
399 orbital_node=direction,
400 mask_water=self.mask_water)
401 scene_patches.append(patch)
403 return scene_patches
405 def get_targets(self):
406 targets = [s.get_target() for s in self.get_scene_patches()]
408 for t in targets:
409 t.store_id = self.store_id
411 return targets
413 def get_insar_scenes(self, engine, sources, tmin=None, tmax=None):
414 logger.debug('Forward modelling InSAR displacement...')
416 scenario_tmin, scenario_tmax = self.get_time_range(sources)
418 try:
419 resp = engine.process(
420 sources,
421 self.get_targets(),
422 nthreads=0)
423 except gf.meta.OutOfBounds:
424 logger.warning('Could not calculate InSAR displacements'
425 ' - the GF store\'s extend is too small!')
426 return []
428 scenes = [res.scene for res in resp.static_results()]
430 tmin, tmax = self.get_time_range(sources)
431 for sc in scenes:
432 sc.meta.time_master = util.to_time_float(tmin)
433 sc.meta.time_slave = util.to_time_float(tmax)
435 scenes_asc = [sc for sc in scenes
436 if sc.config.meta.orbital_node == 'Ascending']
437 scenes_dsc = [sc for sc in scenes
438 if sc.config.meta.orbital_node == 'Descending']
440 def stack_scenes(scenes):
441 base = scenes[0]
442 for sc in scenes[1:]:
443 base += sc
444 return base
446 scene_asc = stack_scenes(scenes_asc)
447 scene_dsc = stack_scenes(scenes_dsc)
449 if self.noise_generator:
450 scene_asc.displacement += self.noise_generator.get_noise(scene_asc)
451 scene_dsc.displacement += self.noise_generator.get_noise(scene_dsc)
453 return scene_asc, scene_dsc
455 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
457 path_insar = op.join(path, 'insar')
458 util.ensuredir(path_insar)
460 tmin, tmax = self.get_time_range(sources)
461 tts = util.time_to_str
463 fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s'
464 % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d')))
466 def scene_fn(track):
467 return fn_tpl.format(orbital_node=track.lower())
469 for track in ('ascending', 'descending'):
470 fn = '%s.yml' % scene_fn(track)
472 if op.exists(fn):
473 logger.debug('Scene exists: %s' % fn)
474 continue
476 scenes = self.get_insar_scenes(engine, sources, tmin, tmax)
477 for sc in scenes:
478 fn = scene_fn(sc.config.meta.orbital_node)
479 logger.debug('Writing %s' % fn)
480 sc.save('%s.npz' % fn)
482 def add_map_artists(self, engine, sources, automap):
483 logger.warning('InSAR mapping is not implemented!')
484 return None