Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/targets/insar.py: 93%
236 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-12 12:19 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Synthetic InSAR data generator.
8'''
10import numpy as num
11import logging
12from os import path as op
14from datetime import datetime
15from pyrocko import gf, util
16from pyrocko import orthodrome as od
17from pyrocko.guts import Float, Timestamp, Tuple, StringChoice, Bool, Object,\
18 String
19from pyrocko.util import num_full
21from .base import TargetGenerator, NoiseGenerator
22from ..base import get_gsshg
24DEFAULT_STORE_ID = 'ak135_static'
26km = 1e3
27d2r = num.pi/180.
29logger = logging.getLogger('pyrocko.scenario.targets.insar')
30guts_prefix = 'pf.scenario'
33class ScenePatch(Object):
34 center_lat = Float.T(
35 help='Center latitude anchor.')
36 center_lon = Float.T(
37 help='center longitude anchor.')
38 time_master = Timestamp.T(
39 help='Timestamp of the master.')
40 time_slave = Timestamp.T(
41 help='Timestamp of the slave.')
42 inclination = Float.T(
43 help='Orbital inclination towards the equatorial plane [deg].')
44 apogee = Float.T(
45 help='Apogee of the satellite in [m].')
46 swath_width = Float.T(
47 help='Swath width in [m].')
48 track_length = Float.T(
49 help='Track length in [m].')
50 incident_angle = Float.T(
51 help='Ground incident angle in [deg].')
52 resolution = Tuple.T(
53 help='Resolution of raster in east x north [px].')
54 orbital_node = StringChoice.T(
55 ['Ascending', 'Descending'],
56 help='Orbit heading.')
57 mask_water = Bool.T(
58 default=True,
59 help='Mask water bodies.')
61 class SatelliteGeneratorTarget(gf.SatelliteTarget):
63 def __init__(self, scene_patch, *args, **kwargs):
64 gf.SatelliteTarget.__init__(self, *args, **kwargs)
66 self.scene_patch = scene_patch
68 def post_process(self, *args, **kwargs):
69 resp = gf.SatelliteTarget.post_process(self, *args, **kwargs)
71 from kite import Scene
72 from kite.scene import SceneConfig, FrameConfig, Meta
74 patch = self.scene_patch
76 grid, _ = patch.get_grid()
78 displacement = num.empty_like(grid)
79 displacement.fill(num.nan)
80 displacement[patch.get_mask()] = resp.result['displacement.los']
82 theta, phi = patch.get_incident_angles()
84 llLat, llLon = patch.get_ll_anchor()
85 urLat, urLon = patch.get_ur_anchor()
86 dLon = num.abs(llLon - urLon) / patch.resolution[0]
87 dLat = num.abs(llLat - urLat) / patch.resolution[1]
89 scene_config = SceneConfig(
90 meta=Meta(
91 scene_title='Pyrocko Scenario Generator - {orbit} ({time})'
92 .format(orbit=self.scene_patch.orbital_node,
93 time=datetime.now()),
94 orbital_node=patch.orbital_node,
95 scene_id='pyrocko_scenario_%s'
96 % self.scene_patch.orbital_node,
97 satellite_name='Sentinel-1 (pyrocko-scenario)'),
98 frame=FrameConfig(
99 llLat=float(llLat),
100 llLon=float(llLon),
101 dN=float(dLat),
102 dE=float(dLon),
103 spacing='degree'))
105 scene = Scene(
106 displacement=displacement,
107 theta=theta,
108 phi=phi,
109 config=scene_config)
111 resp.scene = scene
113 return resp
115 def __init__(self, *args, **kwargs):
116 Object.__init__(self, *args, **kwargs)
117 self._mask_water = None
119 @property
120 def width(self):
121 track_shift = num.abs(
122 num.cos(self.inclination*d2r) * self.track_length)
124 return self.swath_width + track_shift
126 def get_ll_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_ur_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_ul_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_lr_anchor(self):
139 return od.ne_to_latlon(self.center_lat, self.center_lon,
140 -self.track_length/2, self.width/2)
142 def get_corner_coordinates(self):
143 inc = self.inclination
145 llLat, llLon = self.get_ll_anchor()
146 urLat, urLon = self.get_ur_anchor()
148 if self.orbital_node == 'Ascending':
150 ulLat, ulLon = od.ne_to_latlon(
151 self.center_lat, self.center_lon,
152 self.track_length/2,
153 -num.tan(inc*d2r) * self.width/2)
154 lrLat, lrLon = od.ne_to_latlon(
155 self.center_lat, self.center_lon,
156 -self.track_length/2,
157 num.tan(inc*d2r) * self.width/2)
159 elif self.orbital_node == 'Descending':
160 urLat, urLon = od.ne_to_latlon(
161 self.center_lat, self.center_lon,
162 self.track_length/2,
163 num.tan(inc*d2r) * self.width/2)
164 llLat, llLon = od.ne_to_latlon(
165 self.center_lat, self.center_lon,
166 -self.track_length/2,
167 -num.tan(inc*d2r) * self.width/2)
169 return ((llLat, llLon), (ulLat, ulLon),
170 (urLat, urLon), (lrLat, lrLon))
172 def get_grid(self):
173 '''
174 Return relative positions of scatterer.
176 :param track: Acquisition track, from `'asc'` or `'dsc'`.
177 :type track: str
178 '''
179 easts = num.linspace(0, self.width,
180 self.resolution[0])
181 norths = num.linspace(0, self.track_length,
182 self.resolution[1])
184 return num.meshgrid(easts, norths)
186 def get_mask_track(self):
187 east_shifts, north_shifts = self.get_grid()
188 norths = north_shifts[:, 0]
189 track = num.abs(num.cos(self.inclination*d2r)) * norths
191 track_mask = num.logical_and(
192 east_shifts > track[:, num.newaxis],
193 east_shifts < (track + self.swath_width)[:, num.newaxis])
195 if self.orbital_node == 'Ascending':
196 track_mask = num.fliplr(track_mask)
198 return track_mask
200 def get_mask_water(self):
201 if self._mask_water is None:
202 east_shifts, north_shifts = self.get_grid()
204 east_shifts -= east_shifts[0, -1]/2
205 north_shifts -= north_shifts[-1, -1]/2
207 latlon = od.ne_to_latlon(self.center_lat, self.center_lon,
208 north_shifts.ravel(), east_shifts.ravel())
209 points = num.array(latlon).T
210 self._mask_water = get_gsshg().get_land_mask(points)\
211 .reshape(*east_shifts.shape)
213 return self._mask_water
215 def get_mask(self):
216 mask_track = self.get_mask_track()
217 if self.mask_water:
218 mask_water = self.get_mask_water()
219 return num.logical_and(mask_track, mask_water)
220 return mask_track
222 def get_incident_angles(self):
223 # theta: elevation angle towards satellite from horizon in radians.
224 # phi: Horizontal angle towards satellite' :abbr:`line of sight (LOS)`
225 # in [rad] from East.
226 east_shifts, _ = self.get_grid()
228 phi = num.empty_like(east_shifts)
229 theta = num.empty_like(east_shifts)
231 east_shifts += num.tan(self.incident_angle*d2r) * self.apogee
232 theta = num.arctan(east_shifts/self.apogee)
234 if self.orbital_node == 'Ascending':
235 phi.fill(self.inclination*d2r + num.pi/2)
236 elif self.orbital_node == 'Descending':
237 phi.fill(2*num.pi-(self.inclination*d2r + 3/2*num.pi))
238 theta = num.fliplr(theta)
239 else:
240 raise AttributeError(
241 'Orbital node %s not defined!' % self.orbital_node)
243 theta[~self.get_mask()] = num.nan
244 phi[~self.get_mask()] = num.nan
246 return theta, phi
248 def get_target(self):
249 gE, gN = self.get_grid()
250 mask = self.get_mask()
252 east_shifts = gE[mask].ravel()
253 north_shifts = gN[mask].ravel()
254 llLat, llLon = self.get_ll_anchor()
256 ncoords = east_shifts.size
258 theta, phi = self.get_incident_angles()
260 theta = theta[mask].ravel()
261 phi = phi[mask].ravel()
263 if ncoords == 0:
264 logger.warning('InSAR taget has no valid points,'
265 " maybe it's all water?")
267 return self.SatelliteGeneratorTarget(
268 scene_patch=self,
269 lats=num_full(ncoords, fill_value=llLat),
270 lons=num_full(ncoords, fill_value=llLon),
271 east_shifts=east_shifts,
272 north_shifts=north_shifts,
273 theta=theta,
274 phi=phi)
277class AtmosphericNoiseGenerator(NoiseGenerator):
279 amplitude = Float.T(
280 default=1.,
281 help='Amplitude of the atmospheric noise.')
283 beta = [5./3, 8./3, 2./3]
284 regimes = [.15, .99, 1.]
286 def get_noise(self, scene):
287 nE = scene.frame.rows
288 nN = scene.frame.cols
290 if (nE+nN) % 2 != 0:
291 raise ArithmeticError('Dimensions of synthetic scene must '
292 'both be even!')
294 dE = scene.frame.dE
295 dN = scene.frame.dN
297 rfield = num.random.rand(nE, nN)
298 spec = num.fft.fft2(rfield)
300 kE = num.fft.fftfreq(nE, dE)
301 kN = num.fft.fftfreq(nN, dN)
302 k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)
304 regimes = num.array(self.regimes)
306 k0 = 0.
307 k1 = regimes[0] * k_rad.max()
308 k2 = regimes[1] * k_rad.max()
310 r0 = num.logical_and(k_rad > k0, k_rad < k1)
311 r1 = num.logical_and(k_rad >= k1, k_rad < k2)
312 r2 = k_rad >= k2
314 beta = num.array(self.beta)
315 # From Hanssen (2001)
316 # beta+1 is used as beta, since, the power exponent
317 # is defined for a 1D slice of the 2D spectrum:
318 # austin94: "Adler, 1981, shows that the surface profile
319 # created by the intersection of a plane and a
320 # 2-D fractal surface is itself fractal with
321 # a fractal dimension equal to that of the 2D
322 # surface decreased by one."
323 beta += 1.
324 # From Hanssen (2001)
325 # The power beta/2 is used because the power spectral
326 # density is proportional to the amplitude squared
327 # Here we work with the amplitude, instead of the power
328 # so we should take sqrt( k.^beta) = k.^(beta/2) RH
329 # beta /= 2.
331 amp = num.zeros_like(k_rad)
332 amp[r0] = k_rad[r0] ** -beta[0]
333 amp[r0] /= amp[r0].max()
335 amp[r1] = k_rad[r1] ** -beta[1]
336 amp[r1] /= amp[r1].max() / amp[r0].min()
338 amp[r2] = k_rad[r2] ** -beta[2]
339 amp[r2] /= amp[r2].max() / amp[r1].min()
341 amp[k_rad == 0.] = amp.max()
343 spec *= self.amplitude * num.sqrt(amp)
344 noise = num.abs(num.fft.ifft2(spec))
345 noise -= num.mean(noise)
347 return noise
350class InSARGenerator(TargetGenerator):
351 # https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath
352 store_id = String.T(
353 default=DEFAULT_STORE_ID,
354 help='Store ID for these stations.')
356 inclination = Float.T(
357 default=98.2,
358 help='Inclination of the satellite orbit towards equatorial plane'
359 ' in [deg]. Defaults to Sentinel-1 (98.1 deg).')
360 apogee = Float.T(
361 default=693.*km,
362 help='Apogee of the satellite in [m]. '
363 'Defaults to Sentinel-1 (693 km).')
364 swath_width = Float.T(
365 default=250*km,
366 help='Swath width in [m]. '
367 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).'
368 ' (IW; 250 km).')
369 track_length = Float.T(
370 default=250*km,
371 help='Track length in [m]. Defaults to 200 km.')
372 incident_angle = Float.T(
373 default=29.1,
374 help='Near range incident angle in [deg]. Defaults to 29.1 deg;'
375 ' Sentinel IW mode (29.1 - 46.0 deg).')
376 resolution = Tuple.T(
377 default=(250, 250),
378 help='Resolution of raster in east x north [px].')
379 mask_water = Bool.T(
380 default=True,
381 help='Mask out water bodies.')
382 noise_generator = NoiseGenerator.T(
383 default=AtmosphericNoiseGenerator.D(),
384 help='Add atmospheric noise model after Hansen, 2001.')
386 def get_scene_patches(self):
387 center_lat, center_lon = self.get_center_latlon()
389 scene_patches = []
390 for direction in ('Ascending', 'Descending'):
391 patch = ScenePatch(
392 center_lat=center_lat,
393 center_lon=center_lon,
394 time_master=0,
395 time_slave=0,
396 inclination=self.inclination,
397 apogee=self.apogee,
398 swath_width=self.swath_width,
399 track_length=self.track_length,
400 incident_angle=self.incident_angle,
401 resolution=self.resolution,
402 orbital_node=direction,
403 mask_water=self.mask_water)
404 scene_patches.append(patch)
406 return scene_patches
408 def get_targets(self):
409 targets = [s.get_target() for s in self.get_scene_patches()]
411 for t in targets:
412 t.store_id = self.store_id
414 return targets
416 def get_insar_scenes(self, engine, sources, tmin=None, tmax=None):
417 logger.debug('Forward modelling InSAR displacement...')
419 scenario_tmin, scenario_tmax = self.get_time_range(sources)
421 try:
422 resp = engine.process(
423 sources,
424 self.get_targets(),
425 nthreads=0)
426 except gf.meta.OutOfBounds:
427 logger.warning('Could not calculate InSAR displacements'
428 " - the GF store's extend is too small!")
429 return []
431 scenes = [res.scene for res in resp.static_results()]
433 tmin, tmax = self.get_time_range(sources)
434 for sc in scenes:
435 sc.meta.time_master = util.to_time_float(tmin)
436 sc.meta.time_slave = util.to_time_float(tmax)
438 scenes_asc = [sc for sc in scenes
439 if sc.config.meta.orbital_node == 'Ascending']
440 scenes_dsc = [sc for sc in scenes
441 if sc.config.meta.orbital_node == 'Descending']
443 def stack_scenes(scenes):
444 base = scenes[0]
445 for sc in scenes[1:]:
446 base += sc
447 return base
449 scene_asc = stack_scenes(scenes_asc)
450 scene_dsc = stack_scenes(scenes_dsc)
452 if self.noise_generator:
453 scene_asc.displacement += self.noise_generator.get_noise(scene_asc)
454 scene_dsc.displacement += self.noise_generator.get_noise(scene_dsc)
456 return scene_asc, scene_dsc
458 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
460 path_insar = op.join(path, 'insar')
461 util.ensuredir(path_insar)
463 tmin, tmax = self.get_time_range(sources)
464 tts = util.time_to_str
466 fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s'
467 % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d')))
469 def scene_fn(track):
470 return fn_tpl.format(orbital_node=track.lower())
472 for track in ('ascending', 'descending'):
473 fn = '%s.yml' % scene_fn(track)
475 if op.exists(fn):
476 logger.debug('Scene exists: %s' % fn)
477 continue
479 scenes = self.get_insar_scenes(engine, sources, tmin, tmax)
480 for sc in scenes:
481 fn = scene_fn(sc.config.meta.orbital_node)
482 logger.debug('Writing %s' % fn)
483 sc.save('%s.npz' % fn)
485 def add_map_artists(self, engine, sources, automap):
486 logger.warning('InSAR mapping is not implemented!')
487 return None