1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
7import logging
8from os import path as op
10from datetime import datetime
11from pyrocko import gf, util
12from pyrocko import orthodrome as od
13from pyrocko.guts import Float, Timestamp, Tuple, StringChoice, Bool, Object,\
14 String
15from pyrocko.util import num_full
17from .base import TargetGenerator, NoiseGenerator
18from ..base import get_gsshg
20DEFAULT_STORE_ID = 'ak135_static'
22km = 1e3
23d2r = num.pi/180.
25logger = logging.getLogger('pyrocko.scenario.targets.insar')
26guts_prefix = 'pf.scenario'
29class ScenePatch(Object):
30 center_lat = Float.T(
31 help='Center latitude anchor.')
32 center_lon = Float.T(
33 help='center longitude anchor.')
34 time_master = Timestamp.T(
35 help='Timestamp of the master.')
36 time_slave = Timestamp.T(
37 help='Timestamp of the slave.')
38 inclination = Float.T(
39 help='Orbital inclination towards the equatorial plane [deg].')
40 apogee = Float.T(
41 help='Apogee of the satellite in [m].')
42 swath_width = Float.T(
43 help='Swath width in [m].')
44 track_length = Float.T(
45 help='Track length in [m].')
46 incident_angle = Float.T(
47 help='Ground incident angle in [deg].')
48 resolution = Tuple.T(
49 help='Resolution of raster in east x north [px].')
50 orbital_node = StringChoice.T(
51 ['Ascending', 'Descending'],
52 help='Orbit heading.')
53 mask_water = Bool.T(
54 default=True,
55 help='Mask water bodies.')
57 class SatelliteGeneratorTarget(gf.SatelliteTarget):
59 def __init__(self, scene_patch, *args, **kwargs):
60 gf.SatelliteTarget.__init__(self, *args, **kwargs)
62 self.scene_patch = scene_patch
64 def post_process(self, *args, **kwargs):
65 resp = gf.SatelliteTarget.post_process(self, *args, **kwargs)
67 from kite import Scene
68 from kite.scene import SceneConfig, FrameConfig, Meta
70 patch = self.scene_patch
72 grid, _ = patch.get_grid()
74 displacement = num.empty_like(grid)
75 displacement.fill(num.nan)
76 displacement[patch.get_mask()] = resp.result['displacement.los']
78 theta, phi = patch.get_incident_angles()
80 llLat, llLon = patch.get_ll_anchor()
81 urLat, urLon = patch.get_ur_anchor()
82 dLon = num.abs(llLon - urLon) / patch.resolution[0]
83 dLat = num.abs(llLat - urLat) / patch.resolution[1]
85 scene_config = SceneConfig(
86 meta=Meta(
87 scene_title='Pyrocko Scenario Generator - {orbit} ({time})'
88 .format(orbit=self.scene_patch.orbital_node,
89 time=datetime.now()),
90 orbital_node=patch.orbital_node,
91 scene_id='pyrocko_scenario_%s'
92 % self.scene_patch.orbital_node,
93 satellite_name='Sentinel-1 (pyrocko-scenario)'),
94 frame=FrameConfig(
95 llLat=float(llLat),
96 llLon=float(llLon),
97 dN=float(dLat),
98 dE=float(dLon),
99 spacing='degree'))
101 scene = Scene(
102 displacement=displacement,
103 theta=theta,
104 phi=phi,
105 config=scene_config)
107 resp.scene = scene
109 return resp
111 def __init__(self, *args, **kwargs):
112 Object.__init__(self, *args, **kwargs)
113 self._mask_water = None
115 @property
116 def width(self):
117 track_shift = num.abs(
118 num.cos(self.inclination*d2r) * self.track_length)
120 return self.swath_width + track_shift
122 def get_ll_anchor(self):
123 return od.ne_to_latlon(self.center_lat, self.center_lon,
124 -self.track_length/2, -self.width/2)
126 def get_ur_anchor(self):
127 return od.ne_to_latlon(self.center_lat, self.center_lon,
128 self.track_length/2, self.width/2)
130 def get_ul_anchor(self):
131 return od.ne_to_latlon(self.center_lat, self.center_lon,
132 self.track_length/2, -self.width/2)
134 def get_lr_anchor(self):
135 return od.ne_to_latlon(self.center_lat, self.center_lon,
136 -self.track_length/2, self.width/2)
138 def get_corner_coordinates(self):
139 inc = self.inclination
141 llLat, llLon = self.get_ll_anchor()
142 urLat, urLon = self.get_ur_anchor()
144 if self.orbital_node == 'Ascending':
146 ulLat, ulLon = od.ne_to_latlon(
147 self.center_lat, self.center_lon,
148 self.track_length/2,
149 -num.tan(inc*d2r) * self.width/2)
150 lrLat, lrLon = od.ne_to_latlon(
151 self.center_lat, self.center_lon,
152 -self.track_length/2,
153 num.tan(inc*d2r) * self.width/2)
155 elif self.orbital_node == 'Descending':
156 urLat, urLon = od.ne_to_latlon(
157 self.center_lat, self.center_lon,
158 self.track_length/2,
159 num.tan(inc*d2r) * self.width/2)
160 llLat, llLon = od.ne_to_latlon(
161 self.center_lat, self.center_lon,
162 -self.track_length/2,
163 -num.tan(inc*d2r) * self.width/2)
165 return ((llLat, llLon), (ulLat, ulLon),
166 (urLat, urLon), (lrLat, lrLon))
168 def get_grid(self):
169 '''
170 Return relative positions of scatterer.
172 :param track: Acquisition track, from `'asc'` or `'dsc'`.
173 :type track: string
174 '''
175 easts = num.linspace(0, self.width,
176 self.resolution[0])
177 norths = num.linspace(0, self.track_length,
178 self.resolution[1])
180 return num.meshgrid(easts, norths)
182 def get_mask_track(self):
183 east_shifts, north_shifts = self.get_grid()
184 norths = north_shifts[:, 0]
185 track = num.abs(num.cos(self.inclination*d2r)) * norths
187 track_mask = num.logical_and(
188 east_shifts > track[:, num.newaxis],
189 east_shifts < (track + self.swath_width)[:, num.newaxis])
191 if self.orbital_node == 'Ascending':
192 track_mask = num.fliplr(track_mask)
194 return track_mask
196 def get_mask_water(self):
197 if self._mask_water is None:
198 east_shifts, north_shifts = self.get_grid()
200 east_shifts -= east_shifts[0, -1]/2
201 north_shifts -= north_shifts[-1, -1]/2
203 latlon = od.ne_to_latlon(self.center_lat, self.center_lon,
204 north_shifts.ravel(), east_shifts.ravel())
205 points = num.array(latlon).T
206 self._mask_water = get_gsshg().get_land_mask(points)\
207 .reshape(*east_shifts.shape)
209 return self._mask_water
211 def get_mask(self):
212 mask_track = self.get_mask_track()
213 if self.mask_water:
214 mask_water = self.get_mask_water()
215 return num.logical_and(mask_track, mask_water)
216 return mask_track
218 def get_incident_angles(self):
219 # theta: elevation angle towards satellite from horizon in radians.
220 # phi: Horizontal angle towards satellite' :abbr:`line of sight (LOS)`
221 # in [rad] from East.
222 east_shifts, _ = self.get_grid()
224 phi = num.empty_like(east_shifts)
225 theta = num.empty_like(east_shifts)
227 east_shifts += num.tan(self.incident_angle*d2r) * self.apogee
228 theta = num.arctan(east_shifts/self.apogee)
230 if self.orbital_node == 'Ascending':
231 phi.fill(self.inclination*d2r + num.pi/2)
232 elif self.orbital_node == 'Descending':
233 phi.fill(2*num.pi-(self.inclination*d2r + 3/2*num.pi))
234 theta = num.fliplr(theta)
235 else:
236 raise AttributeError(
237 'Orbital node %s not defined!' % self.orbital_node)
239 theta[~self.get_mask()] = num.nan
240 phi[~self.get_mask()] = num.nan
242 return theta, phi
244 def get_target(self):
245 gE, gN = self.get_grid()
246 mask = self.get_mask()
248 east_shifts = gE[mask].ravel()
249 north_shifts = gN[mask].ravel()
250 llLat, llLon = self.get_ll_anchor()
252 ncoords = east_shifts.size
254 theta, phi = self.get_incident_angles()
256 theta = theta[mask].ravel()
257 phi = phi[mask].ravel()
259 if ncoords == 0:
260 logger.warning('InSAR taget has no valid points,'
261 ' maybe it\'s all water?')
263 return self.SatelliteGeneratorTarget(
264 scene_patch=self,
265 lats=num_full(ncoords, fill_value=llLat),
266 lons=num_full(ncoords, fill_value=llLon),
267 east_shifts=east_shifts,
268 north_shifts=north_shifts,
269 theta=theta,
270 phi=phi)
273class AtmosphericNoiseGenerator(NoiseGenerator):
275 amplitude = Float.T(
276 default=1.,
277 help='Amplitude of the atmospheric noise.')
279 beta = [5./3, 8./3, 2./3]
280 regimes = [.15, .99, 1.]
282 def get_noise(self, scene):
283 nE = scene.frame.rows
284 nN = scene.frame.cols
286 if (nE+nN) % 2 != 0:
287 raise ArithmeticError('Dimensions of synthetic scene must '
288 'both be even!')
290 dE = scene.frame.dE
291 dN = scene.frame.dN
293 rfield = num.random.rand(nE, nN)
294 spec = num.fft.fft2(rfield)
296 kE = num.fft.fftfreq(nE, dE)
297 kN = num.fft.fftfreq(nN, dN)
298 k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)
300 regimes = num.array(self.regimes)
302 k0 = 0.
303 k1 = regimes[0] * k_rad.max()
304 k2 = regimes[1] * k_rad.max()
306 r0 = num.logical_and(k_rad > k0, k_rad < k1)
307 r1 = num.logical_and(k_rad >= k1, k_rad < k2)
308 r2 = k_rad >= k2
310 beta = num.array(self.beta)
311 # From Hanssen (2001)
312 # beta+1 is used as beta, since, the power exponent
313 # is defined for a 1D slice of the 2D spectrum:
314 # austin94: "Adler, 1981, shows that the surface profile
315 # created by the intersection of a plane and a
316 # 2-D fractal surface is itself fractal with
317 # a fractal dimension equal to that of the 2D
318 # surface decreased by one."
319 beta += 1.
320 # From Hanssen (2001)
321 # The power beta/2 is used because the power spectral
322 # density is proportional to the amplitude squared
323 # Here we work with the amplitude, instead of the power
324 # so we should take sqrt( k.^beta) = k.^(beta/2) RH
325 # beta /= 2.
327 amp = num.zeros_like(k_rad)
328 amp[r0] = k_rad[r0] ** -beta[0]
329 amp[r0] /= amp[r0].max()
331 amp[r1] = k_rad[r1] ** -beta[1]
332 amp[r1] /= amp[r1].max() / amp[r0].min()
334 amp[r2] = k_rad[r2] ** -beta[2]
335 amp[r2] /= amp[r2].max() / amp[r1].min()
337 amp[k_rad == 0.] = amp.max()
339 spec *= self.amplitude * num.sqrt(amp)
340 noise = num.abs(num.fft.ifft2(spec))
341 noise -= num.mean(noise)
343 return noise
346class InSARGenerator(TargetGenerator):
347 # https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath
348 store_id = String.T(
349 default=DEFAULT_STORE_ID,
350 help='Store ID for these stations.')
352 inclination = Float.T(
353 default=98.2,
354 help='Inclination of the satellite orbit towards equatorial plane'
355 ' in [deg]. Defaults to Sentinel-1 (98.1 deg).')
356 apogee = Float.T(
357 default=693.*km,
358 help='Apogee of the satellite in [m]. '
359 'Defaults to Sentinel-1 (693 km).')
360 swath_width = Float.T(
361 default=250*km,
362 help='Swath width in [m]. '
363 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).'
364 ' (IW; 250 km).')
365 track_length = Float.T(
366 default=250*km,
367 help='Track length in [m]. Defaults to 200 km.')
368 incident_angle = Float.T(
369 default=29.1,
370 help='Near range incident angle in [deg]. Defaults to 29.1 deg;'
371 ' Sentinel IW mode (29.1 - 46.0 deg).')
372 resolution = Tuple.T(
373 default=(250, 250),
374 help='Resolution of raster in east x north [px].')
375 mask_water = Bool.T(
376 default=True,
377 help='Mask out water bodies.')
378 noise_generator = NoiseGenerator.T(
379 default=AtmosphericNoiseGenerator.D(),
380 help='Add atmospheric noise model after Hansen, 2001.')
382 def get_scene_patches(self):
383 center_lat, center_lon = self.get_center_latlon()
385 scene_patches = []
386 for direction in ('Ascending', 'Descending'):
387 patch = ScenePatch(
388 center_lat=center_lat,
389 center_lon=center_lon,
390 time_master=0,
391 time_slave=0,
392 inclination=self.inclination,
393 apogee=self.apogee,
394 swath_width=self.swath_width,
395 track_length=self.track_length,
396 incident_angle=self.incident_angle,
397 resolution=self.resolution,
398 orbital_node=direction,
399 mask_water=self.mask_water)
400 scene_patches.append(patch)
402 return scene_patches
404 def get_targets(self):
405 targets = [s.get_target() for s in self.get_scene_patches()]
407 for t in targets:
408 t.store_id = self.store_id
410 return targets
412 def get_insar_scenes(self, engine, sources, tmin=None, tmax=None):
413 logger.debug('Forward modelling InSAR displacement...')
415 scenario_tmin, scenario_tmax = self.get_time_range(sources)
417 try:
418 resp = engine.process(
419 sources,
420 self.get_targets(),
421 nthreads=0)
422 except gf.meta.OutOfBounds:
423 logger.warning('Could not calculate InSAR displacements'
424 ' - the GF store\'s extend is too small!')
425 return []
427 scenes = [res.scene for res in resp.static_results()]
429 tmin, tmax = self.get_time_range(sources)
430 for sc in scenes:
431 sc.meta.time_master = util.to_time_float(tmin)
432 sc.meta.time_slave = util.to_time_float(tmax)
434 scenes_asc = [sc for sc in scenes
435 if sc.config.meta.orbital_node == 'Ascending']
436 scenes_dsc = [sc for sc in scenes
437 if sc.config.meta.orbital_node == 'Descending']
439 def stack_scenes(scenes):
440 base = scenes[0]
441 for sc in scenes[1:]:
442 base += sc
443 return base
445 scene_asc = stack_scenes(scenes_asc)
446 scene_dsc = stack_scenes(scenes_dsc)
448 if self.noise_generator:
449 scene_asc.displacement += self.noise_generator.get_noise(scene_asc)
450 scene_dsc.displacement += self.noise_generator.get_noise(scene_dsc)
452 return scene_asc, scene_dsc
454 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
456 path_insar = op.join(path, 'insar')
457 util.ensuredir(path_insar)
459 tmin, tmax = self.get_time_range(sources)
460 tts = util.time_to_str
462 fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s'
463 % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d')))
465 def scene_fn(track):
466 return fn_tpl.format(orbital_node=track.lower())
468 for track in ('ascending', 'descending'):
469 fn = '%s.yml' % scene_fn(track)
471 if op.exists(fn):
472 logger.debug('Scene exists: %s' % fn)
473 continue
475 scenes = self.get_insar_scenes(engine, sources, tmin, tmax)
476 for sc in scenes:
477 fn = scene_fn(sc.config.meta.orbital_node)
478 logger.debug('Writing %s' % fn)
479 sc.save('%s.npz' % fn)
481 def add_map_artists(self, engine, sources, automap):
482 logger.warning('InSAR mapping is not implemented!')
483 return None