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---------- 

5 

6''' 

7Synthetic InSAR data generator. 

8''' 

9 

10import numpy as num 

11import logging 

12from os import path as op 

13 

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 

19 

20from .base import TargetGenerator, NoiseGenerator 

21from ..base import get_gsshg 

22 

23DEFAULT_STORE_ID = 'ak135_static' 

24 

25km = 1e3 

26d2r = num.pi/180. 

27 

28logger = logging.getLogger('pyrocko.scenario.targets.insar') 

29guts_prefix = 'pf.scenario' 

30 

31 

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.') 

59 

60 class SatelliteGeneratorTarget(gf.SatelliteTarget): 

61 

62 def __init__(self, scene_patch, *args, **kwargs): 

63 gf.SatelliteTarget.__init__(self, *args, **kwargs) 

64 

65 self.scene_patch = scene_patch 

66 

67 def post_process(self, *args, **kwargs): 

68 resp = gf.SatelliteTarget.post_process(self, *args, **kwargs) 

69 

70 from kite import Scene 

71 from kite.scene import SceneConfig, FrameConfig, Meta 

72 

73 patch = self.scene_patch 

74 

75 grid, _ = patch.get_grid() 

76 

77 displacement = num.empty_like(grid) 

78 displacement.fill(num.nan) 

79 displacement[patch.get_mask()] = resp.result['displacement.los'] 

80 

81 theta, phi = patch.get_incident_angles() 

82 

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] 

87 

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')) 

103 

104 scene = Scene( 

105 displacement=displacement, 

106 theta=theta, 

107 phi=phi, 

108 config=scene_config) 

109 

110 resp.scene = scene 

111 

112 return resp 

113 

114 def __init__(self, *args, **kwargs): 

115 Object.__init__(self, *args, **kwargs) 

116 self._mask_water = None 

117 

118 @property 

119 def width(self): 

120 track_shift = num.abs( 

121 num.cos(self.inclination*d2r) * self.track_length) 

122 

123 return self.swath_width + track_shift 

124 

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) 

128 

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) 

132 

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) 

136 

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) 

140 

141 def get_corner_coordinates(self): 

142 inc = self.inclination 

143 

144 llLat, llLon = self.get_ll_anchor() 

145 urLat, urLon = self.get_ur_anchor() 

146 

147 if self.orbital_node == 'Ascending': 

148 

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) 

157 

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) 

167 

168 return ((llLat, llLon), (ulLat, ulLon), 

169 (urLat, urLon), (lrLat, lrLon)) 

170 

171 def get_grid(self): 

172 ''' 

173 Return relative positions of scatterer. 

174 

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]) 

182 

183 return num.meshgrid(easts, norths) 

184 

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 

189 

190 track_mask = num.logical_and( 

191 east_shifts > track[:, num.newaxis], 

192 east_shifts < (track + self.swath_width)[:, num.newaxis]) 

193 

194 if self.orbital_node == 'Ascending': 

195 track_mask = num.fliplr(track_mask) 

196 

197 return track_mask 

198 

199 def get_mask_water(self): 

200 if self._mask_water is None: 

201 east_shifts, north_shifts = self.get_grid() 

202 

203 east_shifts -= east_shifts[0, -1]/2 

204 north_shifts -= north_shifts[-1, -1]/2 

205 

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) 

211 

212 return self._mask_water 

213 

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 

220 

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() 

226 

227 phi = num.empty_like(east_shifts) 

228 theta = num.empty_like(east_shifts) 

229 

230 east_shifts += num.tan(self.incident_angle*d2r) * self.apogee 

231 theta = num.arctan(east_shifts/self.apogee) 

232 

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) 

241 

242 theta[~self.get_mask()] = num.nan 

243 phi[~self.get_mask()] = num.nan 

244 

245 return theta, phi 

246 

247 def get_target(self): 

248 gE, gN = self.get_grid() 

249 mask = self.get_mask() 

250 

251 east_shifts = gE[mask].ravel() 

252 north_shifts = gN[mask].ravel() 

253 llLat, llLon = self.get_ll_anchor() 

254 

255 ncoords = east_shifts.size 

256 

257 theta, phi = self.get_incident_angles() 

258 

259 theta = theta[mask].ravel() 

260 phi = phi[mask].ravel() 

261 

262 if ncoords == 0: 

263 logger.warning('InSAR taget has no valid points,' 

264 " maybe it's all water?") 

265 

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) 

274 

275 

276class AtmosphericNoiseGenerator(NoiseGenerator): 

277 

278 amplitude = Float.T( 

279 default=1., 

280 help='Amplitude of the atmospheric noise.') 

281 

282 beta = [5./3, 8./3, 2./3] 

283 regimes = [.15, .99, 1.] 

284 

285 def get_noise(self, scene): 

286 nE = scene.frame.rows 

287 nN = scene.frame.cols 

288 

289 if (nE+nN) % 2 != 0: 

290 raise ArithmeticError('Dimensions of synthetic scene must ' 

291 'both be even!') 

292 

293 dE = scene.frame.dE 

294 dN = scene.frame.dN 

295 

296 rfield = num.random.rand(nE, nN) 

297 spec = num.fft.fft2(rfield) 

298 

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) 

302 

303 regimes = num.array(self.regimes) 

304 

305 k0 = 0. 

306 k1 = regimes[0] * k_rad.max() 

307 k2 = regimes[1] * k_rad.max() 

308 

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 

312 

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. 

329 

330 amp = num.zeros_like(k_rad) 

331 amp[r0] = k_rad[r0] ** -beta[0] 

332 amp[r0] /= amp[r0].max() 

333 

334 amp[r1] = k_rad[r1] ** -beta[1] 

335 amp[r1] /= amp[r1].max() / amp[r0].min() 

336 

337 amp[r2] = k_rad[r2] ** -beta[2] 

338 amp[r2] /= amp[r2].max() / amp[r1].min() 

339 

340 amp[k_rad == 0.] = amp.max() 

341 

342 spec *= self.amplitude * num.sqrt(amp) 

343 noise = num.abs(num.fft.ifft2(spec)) 

344 noise -= num.mean(noise) 

345 

346 return noise 

347 

348 

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.') 

354 

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.') 

384 

385 def get_scene_patches(self): 

386 center_lat, center_lon = self.get_center_latlon() 

387 

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) 

404 

405 return scene_patches 

406 

407 def get_targets(self): 

408 targets = [s.get_target() for s in self.get_scene_patches()] 

409 

410 for t in targets: 

411 t.store_id = self.store_id 

412 

413 return targets 

414 

415 def get_insar_scenes(self, engine, sources, tmin=None, tmax=None): 

416 logger.debug('Forward modelling InSAR displacement...') 

417 

418 scenario_tmin, scenario_tmax = self.get_time_range(sources) 

419 

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 [] 

429 

430 scenes = [res.scene for res in resp.static_results()] 

431 

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) 

436 

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'] 

441 

442 def stack_scenes(scenes): 

443 base = scenes[0] 

444 for sc in scenes[1:]: 

445 base += sc 

446 return base 

447 

448 scene_asc = stack_scenes(scenes_asc) 

449 scene_dsc = stack_scenes(scenes_dsc) 

450 

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) 

454 

455 return scene_asc, scene_dsc 

456 

457 def ensure_data(self, engine, sources, path, tmin=None, tmax=None): 

458 

459 path_insar = op.join(path, 'insar') 

460 util.ensuredir(path_insar) 

461 

462 tmin, tmax = self.get_time_range(sources) 

463 tts = util.time_to_str 

464 

465 fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s' 

466 % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d'))) 

467 

468 def scene_fn(track): 

469 return fn_tpl.format(orbital_node=track.lower()) 

470 

471 for track in ('ascending', 'descending'): 

472 fn = '%s.yml' % scene_fn(track) 

473 

474 if op.exists(fn): 

475 logger.debug('Scene exists: %s' % fn) 

476 continue 

477 

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) 

483 

484 def add_map_artists(self, engine, sources, automap): 

485 logger.warning('InSAR mapping is not implemented!') 

486 return None