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-04 09:52 +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 

19from pyrocko.util import num_full 

20 

21from .base import TargetGenerator, NoiseGenerator 

22from ..base import get_gsshg 

23 

24DEFAULT_STORE_ID = 'ak135_static' 

25 

26km = 1e3 

27d2r = num.pi/180. 

28 

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

30guts_prefix = 'pf.scenario' 

31 

32 

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

60 

61 class SatelliteGeneratorTarget(gf.SatelliteTarget): 

62 

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

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

65 

66 self.scene_patch = scene_patch 

67 

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

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

70 

71 from kite import Scene 

72 from kite.scene import SceneConfig, FrameConfig, Meta 

73 

74 patch = self.scene_patch 

75 

76 grid, _ = patch.get_grid() 

77 

78 displacement = num.empty_like(grid) 

79 displacement.fill(num.nan) 

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

81 

82 theta, phi = patch.get_incident_angles() 

83 

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] 

88 

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

104 

105 scene = Scene( 

106 displacement=displacement, 

107 theta=theta, 

108 phi=phi, 

109 config=scene_config) 

110 

111 resp.scene = scene 

112 

113 return resp 

114 

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

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

117 self._mask_water = None 

118 

119 @property 

120 def width(self): 

121 track_shift = num.abs( 

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

123 

124 return self.swath_width + track_shift 

125 

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) 

129 

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) 

133 

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) 

137 

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) 

141 

142 def get_corner_coordinates(self): 

143 inc = self.inclination 

144 

145 llLat, llLon = self.get_ll_anchor() 

146 urLat, urLon = self.get_ur_anchor() 

147 

148 if self.orbital_node == 'Ascending': 

149 

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) 

158 

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) 

168 

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

170 (urLat, urLon), (lrLat, lrLon)) 

171 

172 def get_grid(self): 

173 ''' 

174 Return relative positions of scatterer. 

175 

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

183 

184 return num.meshgrid(easts, norths) 

185 

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 

190 

191 track_mask = num.logical_and( 

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

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

194 

195 if self.orbital_node == 'Ascending': 

196 track_mask = num.fliplr(track_mask) 

197 

198 return track_mask 

199 

200 def get_mask_water(self): 

201 if self._mask_water is None: 

202 east_shifts, north_shifts = self.get_grid() 

203 

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

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

206 

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) 

212 

213 return self._mask_water 

214 

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 

221 

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

227 

228 phi = num.empty_like(east_shifts) 

229 theta = num.empty_like(east_shifts) 

230 

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

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

233 

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) 

242 

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

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

245 

246 return theta, phi 

247 

248 def get_target(self): 

249 gE, gN = self.get_grid() 

250 mask = self.get_mask() 

251 

252 east_shifts = gE[mask].ravel() 

253 north_shifts = gN[mask].ravel() 

254 llLat, llLon = self.get_ll_anchor() 

255 

256 ncoords = east_shifts.size 

257 

258 theta, phi = self.get_incident_angles() 

259 

260 theta = theta[mask].ravel() 

261 phi = phi[mask].ravel() 

262 

263 if ncoords == 0: 

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

265 " maybe it's all water?") 

266 

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) 

275 

276 

277class AtmosphericNoiseGenerator(NoiseGenerator): 

278 

279 amplitude = Float.T( 

280 default=1., 

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

282 

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

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

285 

286 def get_noise(self, scene): 

287 nE = scene.frame.rows 

288 nN = scene.frame.cols 

289 

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

291 raise ArithmeticError('Dimensions of synthetic scene must ' 

292 'both be even!') 

293 

294 dE = scene.frame.dE 

295 dN = scene.frame.dN 

296 

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

298 spec = num.fft.fft2(rfield) 

299 

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) 

303 

304 regimes = num.array(self.regimes) 

305 

306 k0 = 0. 

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

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

309 

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 

313 

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. 

330 

331 amp = num.zeros_like(k_rad) 

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

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

334 

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

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

337 

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

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

340 

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

342 

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

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

345 noise -= num.mean(noise) 

346 

347 return noise 

348 

349 

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

355 

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

385 

386 def get_scene_patches(self): 

387 center_lat, center_lon = self.get_center_latlon() 

388 

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) 

405 

406 return scene_patches 

407 

408 def get_targets(self): 

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

410 

411 for t in targets: 

412 t.store_id = self.store_id 

413 

414 return targets 

415 

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

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

418 

419 scenario_tmin, scenario_tmax = self.get_time_range(sources) 

420 

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

430 

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

432 

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) 

437 

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

442 

443 def stack_scenes(scenes): 

444 base = scenes[0] 

445 for sc in scenes[1:]: 

446 base += sc 

447 return base 

448 

449 scene_asc = stack_scenes(scenes_asc) 

450 scene_dsc = stack_scenes(scenes_dsc) 

451 

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) 

455 

456 return scene_asc, scene_dsc 

457 

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

459 

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

461 util.ensuredir(path_insar) 

462 

463 tmin, tmax = self.get_time_range(sources) 

464 tts = util.time_to_str 

465 

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

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

468 

469 def scene_fn(track): 

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

471 

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

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

474 

475 if op.exists(fn): 

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

477 continue 

478 

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) 

484 

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

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

487 return None