Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/targets/insar.py: 29%
235 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-03-07 11:54 +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
20from .base import TargetGenerator, NoiseGenerator
21from ..base import get_gsshg
23DEFAULT_STORE_ID = 'ak135_static'
25km = 1e3
26d2r = num.pi/180.
28logger = logging.getLogger('pyrocko.scenario.targets.insar')
29guts_prefix = 'pf.scenario'
32class ScenePatch(Object):
33 center_lat = Float.T(
34 help='Center latitude anchor.')
35 center_lon = Float.T(
36 help='center longitude anchor.')
37 time_master = Timestamp.T(
38 help='Timestamp of the master.')
39 time_slave = Timestamp.T(
40 help='Timestamp of the slave.')
41 inclination = Float.T(
42 help='Orbital inclination towards the equatorial plane [deg].')
43 apogee = Float.T(
44 help='Apogee of the satellite in [m].')
45 swath_width = Float.T(
46 help='Swath width in [m].')
47 track_length = Float.T(
48 help='Track length in [m].')
49 incident_angle = Float.T(
50 help='Ground incident angle in [deg].')
51 resolution = Tuple.T(
52 help='Resolution of raster in east x north [px].')
53 orbital_node = StringChoice.T(
54 ['Ascending', 'Descending'],
55 help='Orbit heading.')
56 mask_water = Bool.T(
57 default=True,
58 help='Mask water bodies.')
60 class SatelliteGeneratorTarget(gf.SatelliteTarget):
62 def __init__(self, scene_patch, *args, **kwargs):
63 gf.SatelliteTarget.__init__(self, *args, **kwargs)
65 self.scene_patch = scene_patch
67 def post_process(self, *args, **kwargs):
68 resp = gf.SatelliteTarget.post_process(self, *args, **kwargs)
70 from kite import Scene
71 from kite.scene import SceneConfig, FrameConfig, Meta
73 patch = self.scene_patch
75 grid, _ = patch.get_grid()
77 displacement = num.empty_like(grid)
78 displacement.fill(num.nan)
79 displacement[patch.get_mask()] = resp.result['displacement.los']
81 theta, phi = patch.get_incident_angles()
83 llLat, llLon = patch.get_ll_anchor()
84 urLat, urLon = patch.get_ur_anchor()
85 dLon = num.abs(llLon - urLon) / patch.resolution[0]
86 dLat = num.abs(llLat - urLat) / patch.resolution[1]
88 scene_config = SceneConfig(
89 meta=Meta(
90 scene_title='Pyrocko Scenario Generator - {orbit} ({time})'
91 .format(orbit=self.scene_patch.orbital_node,
92 time=datetime.now()),
93 orbital_node=patch.orbital_node,
94 scene_id='pyrocko_scenario_%s'
95 % self.scene_patch.orbital_node,
96 satellite_name='Sentinel-1 (pyrocko-scenario)'),
97 frame=FrameConfig(
98 llLat=float(llLat),
99 llLon=float(llLon),
100 dN=float(dLat),
101 dE=float(dLon),
102 spacing='degree'))
104 scene = Scene(
105 displacement=displacement,
106 theta=theta,
107 phi=phi,
108 config=scene_config)
110 resp.scene = scene
112 return resp
114 def __init__(self, *args, **kwargs):
115 Object.__init__(self, *args, **kwargs)
116 self._mask_water = None
118 @property
119 def width(self):
120 track_shift = num.abs(
121 num.cos(self.inclination*d2r) * self.track_length)
123 return self.swath_width + track_shift
125 def get_ll_anchor(self):
126 return od.ne_to_latlon(self.center_lat, self.center_lon,
127 -self.track_length/2, -self.width/2)
129 def get_ur_anchor(self):
130 return od.ne_to_latlon(self.center_lat, self.center_lon,
131 self.track_length/2, self.width/2)
133 def get_ul_anchor(self):
134 return od.ne_to_latlon(self.center_lat, self.center_lon,
135 self.track_length/2, -self.width/2)
137 def get_lr_anchor(self):
138 return od.ne_to_latlon(self.center_lat, self.center_lon,
139 -self.track_length/2, self.width/2)
141 def get_corner_coordinates(self):
142 inc = self.inclination
144 llLat, llLon = self.get_ll_anchor()
145 urLat, urLon = self.get_ur_anchor()
147 if self.orbital_node == 'Ascending':
149 ulLat, ulLon = od.ne_to_latlon(
150 self.center_lat, self.center_lon,
151 self.track_length/2,
152 -num.tan(inc*d2r) * self.width/2)
153 lrLat, lrLon = od.ne_to_latlon(
154 self.center_lat, self.center_lon,
155 -self.track_length/2,
156 num.tan(inc*d2r) * self.width/2)
158 elif self.orbital_node == 'Descending':
159 urLat, urLon = od.ne_to_latlon(
160 self.center_lat, self.center_lon,
161 self.track_length/2,
162 num.tan(inc*d2r) * self.width/2)
163 llLat, llLon = od.ne_to_latlon(
164 self.center_lat, self.center_lon,
165 -self.track_length/2,
166 -num.tan(inc*d2r) * self.width/2)
168 return ((llLat, llLon), (ulLat, ulLon),
169 (urLat, urLon), (lrLat, lrLon))
171 def get_grid(self):
172 '''
173 Return relative positions of scatterer.
175 :param track: Acquisition track, from `'asc'` or `'dsc'`.
176 :type track: str
177 '''
178 easts = num.linspace(0, self.width,
179 self.resolution[0])
180 norths = num.linspace(0, self.track_length,
181 self.resolution[1])
183 return num.meshgrid(easts, norths)
185 def get_mask_track(self):
186 east_shifts, north_shifts = self.get_grid()
187 norths = north_shifts[:, 0]
188 track = num.abs(num.cos(self.inclination*d2r)) * norths
190 track_mask = num.logical_and(
191 east_shifts > track[:, num.newaxis],
192 east_shifts < (track + self.swath_width)[:, num.newaxis])
194 if self.orbital_node == 'Ascending':
195 track_mask = num.fliplr(track_mask)
197 return track_mask
199 def get_mask_water(self):
200 if self._mask_water is None:
201 east_shifts, north_shifts = self.get_grid()
203 east_shifts -= east_shifts[0, -1]/2
204 north_shifts -= north_shifts[-1, -1]/2
206 latlon = od.ne_to_latlon(self.center_lat, self.center_lon,
207 north_shifts.ravel(), east_shifts.ravel())
208 points = num.array(latlon).T
209 self._mask_water = get_gsshg().get_land_mask(points)\
210 .reshape(*east_shifts.shape)
212 return self._mask_water
214 def get_mask(self):
215 mask_track = self.get_mask_track()
216 if self.mask_water:
217 mask_water = self.get_mask_water()
218 return num.logical_and(mask_track, mask_water)
219 return mask_track
221 def get_incident_angles(self):
222 # theta: elevation angle towards satellite from horizon in radians.
223 # phi: Horizontal angle towards satellite' :abbr:`line of sight (LOS)`
224 # in [rad] from East.
225 east_shifts, _ = self.get_grid()
227 phi = num.empty_like(east_shifts)
228 theta = num.empty_like(east_shifts)
230 east_shifts += num.tan(self.incident_angle*d2r) * self.apogee
231 theta = num.arctan(east_shifts/self.apogee)
233 if self.orbital_node == 'Ascending':
234 phi.fill(self.inclination*d2r + num.pi/2)
235 elif self.orbital_node == 'Descending':
236 phi.fill(2*num.pi-(self.inclination*d2r + 3/2*num.pi))
237 theta = num.fliplr(theta)
238 else:
239 raise AttributeError(
240 'Orbital node %s not defined!' % self.orbital_node)
242 theta[~self.get_mask()] = num.nan
243 phi[~self.get_mask()] = num.nan
245 return theta, phi
247 def get_target(self):
248 gE, gN = self.get_grid()
249 mask = self.get_mask()
251 east_shifts = gE[mask].ravel()
252 north_shifts = gN[mask].ravel()
253 llLat, llLon = self.get_ll_anchor()
255 ncoords = east_shifts.size
257 theta, phi = self.get_incident_angles()
259 theta = theta[mask].ravel()
260 phi = phi[mask].ravel()
262 if ncoords == 0:
263 logger.warning('InSAR taget has no valid points,'
264 " maybe it's all water?")
266 return self.SatelliteGeneratorTarget(
267 scene_patch=self,
268 lats=num.full(ncoords, fill_value=llLat),
269 lons=num.full(ncoords, fill_value=llLon),
270 east_shifts=east_shifts,
271 north_shifts=north_shifts,
272 theta=theta,
273 phi=phi)
276class AtmosphericNoiseGenerator(NoiseGenerator):
278 amplitude = Float.T(
279 default=1.,
280 help='Amplitude of the atmospheric noise.')
282 beta = [5./3, 8./3, 2./3]
283 regimes = [.15, .99, 1.]
285 def get_noise(self, scene):
286 nE = scene.frame.rows
287 nN = scene.frame.cols
289 if (nE+nN) % 2 != 0:
290 raise ArithmeticError('Dimensions of synthetic scene must '
291 'both be even!')
293 dE = scene.frame.dE
294 dN = scene.frame.dN
296 rfield = num.random.rand(nE, nN)
297 spec = num.fft.fft2(rfield)
299 kE = num.fft.fftfreq(nE, dE)
300 kN = num.fft.fftfreq(nN, dN)
301 k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)
303 regimes = num.array(self.regimes)
305 k0 = 0.
306 k1 = regimes[0] * k_rad.max()
307 k2 = regimes[1] * k_rad.max()
309 r0 = num.logical_and(k_rad > k0, k_rad < k1)
310 r1 = num.logical_and(k_rad >= k1, k_rad < k2)
311 r2 = k_rad >= k2
313 beta = num.array(self.beta)
314 # From Hanssen (2001)
315 # beta+1 is used as beta, since, the power exponent
316 # is defined for a 1D slice of the 2D spectrum:
317 # austin94: "Adler, 1981, shows that the surface profile
318 # created by the intersection of a plane and a
319 # 2-D fractal surface is itself fractal with
320 # a fractal dimension equal to that of the 2D
321 # surface decreased by one."
322 beta += 1.
323 # From Hanssen (2001)
324 # The power beta/2 is used because the power spectral
325 # density is proportional to the amplitude squared
326 # Here we work with the amplitude, instead of the power
327 # so we should take sqrt( k.^beta) = k.^(beta/2) RH
328 # beta /= 2.
330 amp = num.zeros_like(k_rad)
331 amp[r0] = k_rad[r0] ** -beta[0]
332 amp[r0] /= amp[r0].max()
334 amp[r1] = k_rad[r1] ** -beta[1]
335 amp[r1] /= amp[r1].max() / amp[r0].min()
337 amp[r2] = k_rad[r2] ** -beta[2]
338 amp[r2] /= amp[r2].max() / amp[r1].min()
340 amp[k_rad == 0.] = amp.max()
342 spec *= self.amplitude * num.sqrt(amp)
343 noise = num.abs(num.fft.ifft2(spec))
344 noise -= num.mean(noise)
346 return noise
349class InSARGenerator(TargetGenerator):
350 # https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath
351 store_id = String.T(
352 default=DEFAULT_STORE_ID,
353 help='Store ID for these stations.')
355 inclination = Float.T(
356 default=98.2,
357 help='Inclination of the satellite orbit towards equatorial plane'
358 ' in [deg]. Defaults to Sentinel-1 (98.1 deg).')
359 apogee = Float.T(
360 default=693.*km,
361 help='Apogee of the satellite in [m]. '
362 'Defaults to Sentinel-1 (693 km).')
363 swath_width = Float.T(
364 default=250*km,
365 help='Swath width in [m]. '
366 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).'
367 ' (IW; 250 km).')
368 track_length = Float.T(
369 default=250*km,
370 help='Track length in [m]. Defaults to 200 km.')
371 incident_angle = Float.T(
372 default=29.1,
373 help='Near range incident angle in [deg]. Defaults to 29.1 deg;'
374 ' Sentinel IW mode (29.1 - 46.0 deg).')
375 resolution = Tuple.T(
376 default=(250, 250),
377 help='Resolution of raster in east x north [px].')
378 mask_water = Bool.T(
379 default=True,
380 help='Mask out water bodies.')
381 noise_generator = NoiseGenerator.T(
382 default=AtmosphericNoiseGenerator.D(),
383 help='Add atmospheric noise model after Hansen, 2001.')
385 def get_scene_patches(self):
386 center_lat, center_lon = self.get_center_latlon()
388 scene_patches = []
389 for direction in ('Ascending', 'Descending'):
390 patch = ScenePatch(
391 center_lat=center_lat,
392 center_lon=center_lon,
393 time_master=0,
394 time_slave=0,
395 inclination=self.inclination,
396 apogee=self.apogee,
397 swath_width=self.swath_width,
398 track_length=self.track_length,
399 incident_angle=self.incident_angle,
400 resolution=self.resolution,
401 orbital_node=direction,
402 mask_water=self.mask_water)
403 scene_patches.append(patch)
405 return scene_patches
407 def get_targets(self):
408 targets = [s.get_target() for s in self.get_scene_patches()]
410 for t in targets:
411 t.store_id = self.store_id
413 return targets
415 def get_insar_scenes(self, engine, sources, tmin=None, tmax=None):
416 logger.debug('Forward modelling InSAR displacement...')
418 scenario_tmin, scenario_tmax = self.get_time_range(sources)
420 try:
421 resp = engine.process(
422 sources,
423 self.get_targets(),
424 nthreads=0)
425 except gf.meta.OutOfBounds:
426 logger.warning('Could not calculate InSAR displacements'
427 " - the GF store's extend is too small!")
428 return []
430 scenes = [res.scene for res in resp.static_results()]
432 tmin, tmax = self.get_time_range(sources)
433 for sc in scenes:
434 sc.meta.time_master = util.to_time_float(tmin)
435 sc.meta.time_slave = util.to_time_float(tmax)
437 scenes_asc = [sc for sc in scenes
438 if sc.config.meta.orbital_node == 'Ascending']
439 scenes_dsc = [sc for sc in scenes
440 if sc.config.meta.orbital_node == 'Descending']
442 def stack_scenes(scenes):
443 base = scenes[0]
444 for sc in scenes[1:]:
445 base += sc
446 return base
448 scene_asc = stack_scenes(scenes_asc)
449 scene_dsc = stack_scenes(scenes_dsc)
451 if self.noise_generator:
452 scene_asc.displacement += self.noise_generator.get_noise(scene_asc)
453 scene_dsc.displacement += self.noise_generator.get_noise(scene_dsc)
455 return scene_asc, scene_dsc
457 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
459 path_insar = op.join(path, 'insar')
460 util.ensuredir(path_insar)
462 tmin, tmax = self.get_time_range(sources)
463 tts = util.time_to_str
465 fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s'
466 % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d')))
468 def scene_fn(track):
469 return fn_tpl.format(orbital_node=track.lower())
471 for track in ('ascending', 'descending'):
472 fn = '%s.yml' % scene_fn(track)
474 if op.exists(fn):
475 logger.debug('Scene exists: %s' % fn)
476 continue
478 scenes = self.get_insar_scenes(engine, sources, tmin, tmax)
479 for sc in scenes:
480 fn = scene_fn(sc.config.meta.orbital_node)
481 logger.debug('Writing %s' % fn)
482 sc.save('%s.npz' % fn)
484 def add_map_artists(self, engine, sources, automap):
485 logger.warning('InSAR mapping is not implemented!')
486 return None