Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/satellite/plot.py: 13%

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import logging 

5import numpy as num 

6from matplotlib import cm, gridspec 

7 

8from grond.plot.config import PlotConfig 

9from grond.plot.collection import PlotItem 

10 

11from matplotlib import pyplot as plt 

12from matplotlib.ticker import MaxNLocator, FuncFormatter 

13from matplotlib import patches 

14from pyrocko.guts import Tuple, Float, String, Int, Bool, StringChoice 

15 

16logger = logging.getLogger('grond.targets.satellite.plot') 

17 

18km = 1e3 

19d2r = num.pi/180. 

20guts_prefix = 'grond' 

21 

22 

23def drape_displacements( 

24 displacement, shad_data, mappable, 

25 shad_lim=(.4, .98), contrast=1., mask=None): 

26 '''Map color data (displacement) on shaded relief.''' 

27 

28 from scipy.ndimage import convolve as im_conv 

29 # Light source from somewhere above - psychologically the best choice 

30 # from upper left 

31 ramp = num.array([[1, 0], [0, -1.]]) * contrast 

32 

33 # convolution of two 2D arrays 

34 shad = im_conv(shad_data*km, ramp.T) 

35 shad *= -1. 

36 

37 # if there are strong artifical edges in the data, shades get 

38 # dominated by them. Cutting off the largest and smallest 2% of 

39 # shades helps 

40 percentile2 = num.quantile(shad, 0.02) 

41 percentile98 = num.quantile(shad, 0.98) 

42 shad[shad > percentile98] = percentile98 

43 shad[shad < percentile2] = percentile2 

44 

45 # normalize shading 

46 shad -= num.nanmin(shad) 

47 shad /= num.nanmax(shad) 

48 

49 if mask is not None: 

50 shad[mask] = num.nan 

51 

52 # reduce range to balance gray color 

53 shad *= shad_lim[1] - shad_lim[0] 

54 shad += shad_lim[0] 

55 

56 rgb_map = mappable.to_rgba(displacement) 

57 rgb_map[num.isnan(displacement)] = 1. 

58 rgb_map[:, :, :3] *= shad[:, :, num.newaxis] 

59 

60 return rgb_map 

61 

62 

63def displ2rad(displ, wavelength): 

64 return (displ % wavelength) / wavelength * num.pi 

65 

66 

67def scale_axes(axis, scale, offset=0., suffix=''): 

68 from matplotlib.ticker import ScalarFormatter 

69 

70 class FormatScaled(ScalarFormatter): 

71 

72 @staticmethod 

73 def __call__(value, pos): 

74 return '{:,.1f}{:}'.format((offset + value) * scale, suffix)\ 

75 .replace(',', ' ') 

76 

77 axis.set_major_formatter(FormatScaled()) 

78 

79 

80class SatelliteTargetDisplacement(PlotConfig): 

81 ''' Maps showing surface displacements from satellite and modelled data ''' 

82 

83 name = 'satellite' 

84 dpi = Int.T( 

85 default=250) 

86 size_cm = Tuple.T( 

87 2, Float.T(), 

88 default=(22., 12.)) 

89 colormap = String.T( 

90 default='RdBu', 

91 help='Colormap for the surface displacements') 

92 relative_coordinates = Bool.T( 

93 default=False, 

94 help='Show relative coordinates, initial location centered at 0N, 0E') 

95 fit = StringChoice.T( 

96 default='best', choices=['best', 'mean'], 

97 help='Show the \'best\' or \'mean\' fits and source model from the' 

98 ' ensamble.') 

99 

100 show_topo = Bool.T( 

101 default=True, 

102 help='Drape displacements over the topography.') 

103 displacement_unit = StringChoice.T( 

104 default='m', 

105 choices=['m', 'mm', 'cm', 'rad'], 

106 help="Show results in 'm', 'cm', 'mm' or 'rad' for radians.") 

107 show_leaf_centres = Bool.T( 

108 default=True, 

109 help='show the center points of Quadtree leaves') 

110 source_outline_color = String.T( 

111 default='grey', 

112 help='Choose color of source outline from named matplotlib Colors') 

113 common_color_scale = Bool.T( 

114 default=True, 

115 help='Results shown with common color scale for all satellite ' 

116 'data sets (based on the data)') 

117 map_limits = Tuple.T( 

118 4, Float.T(), 

119 optional=True, 

120 help='Overwrite map limits in native coordinates. ' 

121 'Use (xmin, xmax, ymin, ymax)') 

122 nticks_x = Int.T( 

123 optional=True, 

124 help='Number of ticks on the x-axis.') 

125 

126 def make(self, environ): 

127 cm = environ.get_plot_collection_manager() 

128 history = environ.get_history(subset='harvest') 

129 optimiser = environ.get_optimiser() 

130 ds = environ.get_dataset() 

131 

132 environ.setup_modelling() 

133 

134 cm.create_group_mpl( 

135 self, 

136 self.draw_static_fits(ds, history, optimiser), 

137 title=u'InSAR Displacements', 

138 section='fits', 

139 feather_icon='navigation', 

140 description=u''' 

141Maps showing subsampled surface displacements as observed, modelled and the 

142residual (observed minus modelled). 

143 

144The displacement values predicted by the orbit-ambiguity ramps are added to the 

145modelled displacements (middle panels). The color shows the LOS displacement 

146values associated with, and the extent of, every quadtree box. The light grey 

147dots show the focal point of pixels combined in the quadtree box. This point 

148corresponds to the position of the modelled data point. 

149 

150The large dark grey dot shows the reference source position. The grey filled 

151box shows the surface projection of the modelled source, with the thick-lined 

152edge marking the upper fault edge. Complete data extent is shown. 

153''') 

154 

155 def draw_static_fits(self, ds, history, optimiser, closeup=False): 

156 from pyrocko.orthodrome import latlon_to_ne_numpy 

157 problem = history.problem 

158 

159 sat_targets = problem.satellite_targets 

160 for target in sat_targets: 

161 target.set_dataset(ds) 

162 

163 if self.fit == 'best': 

164 source = history.get_best_source() 

165 model = history.get_best_model() 

166 elif self.fit == 'mean': 

167 source = history.get_mean_source() 

168 model = history.get_mean_model() 

169 

170 results = problem.evaluate(model, targets=sat_targets) 

171 

172 def init_axes(ax, scene, title, last_axes=False): 

173 ax.set_title(title, fontsize=self.font_size) 

174 ax.tick_params(length=2) 

175 

176 if scene.frame.isMeter(): 

177 import utm 

178 ax.set_xlabel('Easting [km]', fontsize=self.font_size) 

179 scale_x = dict(scale=1./km) 

180 scale_y = dict(scale=1./km) 

181 utm_E, utm_N, utm_zone, utm_zone_letter =\ 

182 utm.from_latlon(source.effective_lat, 

183 source.effective_lon) 

184 scale_x['offset'] = utm_E 

185 scale_y['offset'] = utm_N 

186 

187 if last_axes: 

188 ax.text(0.975, 0.025, 

189 'UTM Zone %d%s' % (utm_zone, utm_zone_letter), 

190 va='bottom', ha='right', 

191 fontsize=8, alpha=.7, 

192 transform=ax.transAxes) 

193 ax.set_aspect('equal') 

194 

195 elif scene.frame.isDegree(): 

196 scale_x = dict(scale=1., suffix='°') 

197 scale_y = dict(scale=1., suffix='°') 

198 scale_x['offset'] = source.effective_lon 

199 scale_y['offset'] = source.effective_lat 

200 

201 ax.set_aspect(1./num.cos(source.effective_lat*d2r)) 

202 

203 if self.relative_coordinates: 

204 scale_x['offset'] = 0. 

205 scale_y['offset'] = 0. 

206 

207 nticks_x = 4 if abs(scene.frame.llLon) >= 100 else 5 

208 

209 ax.xaxis.set_major_locator(MaxNLocator(self.nticks_x or nticks_x)) 

210 ax.yaxis.set_major_locator(MaxNLocator(5)) 

211 

212 ax.scale_x = scale_x 

213 ax.scale_y = scale_y 

214 

215 scale_axes(ax.get_xaxis(), **scale_x) 

216 scale_axes(ax.get_yaxis(), **scale_y) 

217 

218 def draw_source(ax, scene): 

219 if scene.frame.isMeter(): 

220 fn, fe = source.outline(cs='xy').T 

221 fn -= fn.mean() 

222 fe -= fe.mean() 

223 elif scene.frame.isDegree(): 

224 fn, fe = source.outline(cs='latlon').T 

225 fn -= source.effective_lat 

226 fe -= source.effective_lon 

227 

228 # source is centered 

229 ax.scatter(0., 0., color='black', s=3, alpha=.5, marker='o') 

230 ax.fill(fe, fn, 

231 edgecolor=(0., 0., 0.), 

232 facecolor=self.source_outline_color, 

233 alpha=0.7) 

234 ax.plot(fe[0:2], fn[0:2], 'k', linewidth=1.3) 

235 

236 def get_displacement_rgba(displacements, scene, mappable): 

237 arr = num.full_like(scene.displacement, fill_value=num.nan) 

238 qt = scene.quadtree 

239 

240 for syn_v, leaf in zip(displacements, qt.leaves): 

241 arr[leaf._slice_rows, leaf._slice_cols] = syn_v 

242 

243 arr[scene.displacement_mask] = num.nan 

244 

245 if not self.common_color_scale \ 

246 and not self.displacement_unit == 'rad': 

247 abs_displ = num.abs(displacements).max() 

248 mappable.set_clim(-abs_displ, abs_displ) 

249 

250 if self.show_topo: 

251 try: 

252 elevation = scene.get_elevation() 

253 return drape_displacements(arr, elevation, mappable) 

254 except Exception as e: 

255 logger.warning('could not plot hillshaded topo') 

256 logger.exception(e) 

257 rgb_arr = mappable.to_rgba(arr) 

258 rgb_arr[num.isnan(arr)] = 1. 

259 rgb_arr[scene.displacement_mask] = 1. 

260 

261 return rgb_arr 

262 

263 def draw_leaves(ax, scene, offset_e=0., offset_n=0.): 

264 rects = scene.quadtree.getMPLRectangles() 

265 for r in rects: 

266 r.set_edgecolor((.4, .4, .4)) 

267 r.set_linewidth(.5) 

268 r.set_facecolor('none') 

269 r.set_x(r.get_x() - offset_e) 

270 r.set_y(r.get_y() - offset_n) 

271 map(ax.add_artist, rects) 

272 

273 if self.show_leaf_centres: 

274 ax.scatter(scene.quadtree.leaf_coordinates[:, 0] - offset_e, 

275 scene.quadtree.leaf_coordinates[:, 1] - offset_n, 

276 s=.25, c='black', alpha=.1) 

277 

278 def add_arrow(ax, scene): 

279 phi = num.nanmean(scene.phi) 

280 los_dx = num.cos(phi + num.pi) * .0625 

281 los_dy = num.sin(phi + num.pi) * .0625 

282 

283 az_dx = num.cos(phi - num.pi/2) * .125 

284 az_dy = num.sin(phi - num.pi/2) * .125 

285 

286 anchor_x = .9 if los_dx < 0 else .1 

287 anchor_y = .85 if los_dx < 0 else .975 

288 

289 az_arrow = patches.FancyArrow( 

290 x=anchor_x-az_dx, y=anchor_y-az_dy, 

291 dx=az_dx, dy=az_dy, 

292 head_width=.025, 

293 alpha=.5, fc='k', 

294 head_starts_at_zero=False, 

295 length_includes_head=True, 

296 transform=ax.transAxes) 

297 

298 los_arrow = patches.FancyArrow( 

299 x=anchor_x-az_dx/2, y=anchor_y-az_dy/2, 

300 dx=los_dx, dy=los_dy, 

301 head_width=.02, 

302 alpha=.5, fc='k', 

303 head_starts_at_zero=False, 

304 length_includes_head=True, 

305 transform=ax.transAxes) 

306 

307 ax.add_artist(az_arrow) 

308 ax.add_artist(los_arrow) 

309 

310 urE, urN, llE, llN = (0., 0., 0., 0.) 

311 for target in sat_targets: 

312 

313 if target.scene.frame.isMeter(): 

314 off_n, off_e = map(float, latlon_to_ne_numpy( 

315 target.scene.frame.llLat, target.scene.frame.llLon, 

316 source.effective_lat, source.effective_lon)) 

317 if target.scene.frame.isDegree(): 

318 off_n = source.effective_lat - target.scene.frame.llLat 

319 off_e = source.effective_lon - target.scene.frame.llLon 

320 

321 turE, turN, tllE, tllN = zip( 

322 *[(leaf.gridE.max()-off_e, 

323 leaf.gridN.max()-off_n, 

324 leaf.gridE.min()-off_e, 

325 leaf.gridN.min()-off_n) 

326 for leaf in target.scene.quadtree.leaves]) 

327 

328 turE, turN = map(max, (turE, turN)) 

329 tllE, tllN = map(min, (tllE, tllN)) 

330 urE, urN = map(max, ((turE, urE), (urN, turN))) 

331 llE, llN = map(min, ((tllE, llE), (llN, tllN))) 

332 

333 def generate_plot(sat_target, result, ifig): 

334 

335 scene = sat_target.scene 

336 

337 fig = plt.figure() 

338 fig.set_size_inches(*self.size_inch) 

339 gs = gridspec.GridSpec( 

340 2, 3, 

341 wspace=.15, hspace=.2, 

342 left=.1, right=.975, top=.95, 

343 height_ratios=[12, 1]) 

344 

345 item = PlotItem( 

346 name='fig_%i' % ifig, 

347 attributes={'targets': [sat_target.path]}, 

348 title=u'Satellite Surface Displacements - %s' 

349 % scene.meta.scene_title, 

350 description=u''' 

351Surface displacements derived from satellite data. 

352(Left) the input data, (center) the modelled 

353data and (right) the model residual. 

354''') 

355 

356 stat_obs = result.statics_obs['displacement.los'] 

357 stat_syn = result.statics_syn['displacement.los'] 

358 res = stat_obs - stat_syn 

359 

360 if scene.frame.isMeter(): 

361 offset_n, offset_e = map(float, latlon_to_ne_numpy( 

362 scene.frame.llLat, scene.frame.llLon, 

363 source.effective_lat, source.effective_lon)) 

364 elif scene.frame.isDegree(): 

365 offset_n = source.effective_lat - scene.frame.llLat 

366 offset_e = source.effective_lon - scene.frame.llLon 

367 

368 im_extent = ( 

369 scene.frame.E.min() - offset_e, 

370 scene.frame.E.max() - offset_e, 

371 scene.frame.N.min() - offset_n, 

372 scene.frame.N.max() - offset_n) 

373 

374 if self.displacement_unit == 'rad': 

375 wavelength = scene.meta.wavelength 

376 if wavelength is None: 

377 raise AttributeError( 

378 'The satellite\'s wavelength is not set') 

379 

380 stat_obs = displ2rad(stat_obs, wavelength) 

381 stat_syn = displ2rad(stat_syn, wavelength) 

382 res = displ2rad(res, wavelength) 

383 

384 self.colormap = 'hsv' 

385 data_range = (0., num.pi) 

386 

387 else: 

388 abs_displ = num.abs([stat_obs.min(), stat_obs.max(), 

389 stat_syn.min(), stat_syn.max(), 

390 res.min(), res.max()]).max() 

391 data_range = (-abs_displ, abs_displ) 

392 

393 cmw = cm.ScalarMappable(cmap=self.colormap) 

394 cmw.set_clim(*data_range) 

395 cmw.set_array(stat_obs) 

396 

397 axes = [fig.add_subplot(gs[0, 0]), 

398 fig.add_subplot(gs[0, 1]), 

399 fig.add_subplot(gs[0, 2])] 

400 

401 ax = axes[0] 

402 ax.imshow( 

403 get_displacement_rgba(stat_obs, scene, cmw), 

404 extent=im_extent, origin='lower') 

405 draw_leaves(ax, scene, offset_e, offset_n) 

406 draw_source(ax, scene) 

407 add_arrow(ax, scene) 

408 init_axes(ax, scene, 'Observed') 

409 

410 ax.text(.025, .025, 'Scene ID: %s' % scene.meta.scene_id, 

411 fontsize=8, alpha=.7, 

412 va='bottom', transform=ax.transAxes) 

413 if scene.frame.isMeter(): 

414 ax.set_ylabel('Northing [km]', fontsize=self.font_size) 

415 

416 ax = axes[1] 

417 ax.imshow( 

418 get_displacement_rgba(stat_syn, scene, cmw), 

419 extent=im_extent, origin='lower') 

420 draw_leaves(ax, scene, offset_e, offset_n) 

421 draw_source(ax, scene) 

422 add_arrow(ax, scene) 

423 init_axes(ax, scene, 'Model') 

424 ax.get_yaxis().set_visible(False) 

425 

426 ax = axes[2] 

427 ax.imshow( 

428 get_displacement_rgba(res, scene, cmw), 

429 extent=im_extent, origin='lower') 

430 

431 draw_leaves(ax, scene, offset_e, offset_n) 

432 draw_source(ax, scene) 

433 add_arrow(ax, scene) 

434 init_axes(ax, scene, 'Residual', last_axes=True) 

435 ax.get_yaxis().set_visible(False) 

436 

437 for ax in axes: 

438 ax.set_xlim(*im_extent[:2]) 

439 ax.set_ylim(*im_extent[2:]) 

440 

441 if closeup: 

442 if scene.frame.isMeter(): 

443 fn, fe = source.outline(cs='xy').T 

444 elif scene.frame.isDegree(): 

445 fn, fe = source.outline(cs='latlon').T 

446 fn -= source.effective_lat 

447 fe -= source.effective_lon 

448 

449 if fn.size > 1: 

450 off_n = (fn[0] + fn[1]) / 2 

451 off_e = (fe[0] + fe[1]) / 2 

452 else: 

453 off_n = fn[0] 

454 off_e = fe[0] 

455 

456 fault_size = 2*num.sqrt(max(abs(fn-off_n))**2 

457 + max(abs(fe-off_e))**2) 

458 fault_size *= self.map_scale 

459 if fault_size == 0.0: 

460 extent = (scene.frame.N[-1] + scene.frame.E[-1]) / 2 

461 fault_size = extent * .25 

462 

463 for ax in axes: 

464 ax.set_xlim(-fault_size/2 + off_e, fault_size/2 + off_e) 

465 ax.set_ylim(-fault_size/2 + off_n, fault_size/2 + off_n) 

466 

467 if self.map_limits is not None: 

468 xmin, xmax, ymin, ymax = self.map_limits 

469 assert xmin < xmax, 'bad map_limits xmin > xmax' 

470 assert ymin < ymax, 'bad map_limits ymin > ymax' 

471 

472 for ax in axes: 

473 ax.set_xlim( 

474 xmin/ax.scale_x['scale'] - ax.scale_x['offset'], 

475 xmax/ax.scale_x['scale'] - ax.scale_x['offset'],) 

476 ax.set_ylim( 

477 ymin/ax.scale_y['scale'] - ax.scale_y['offset'], 

478 ymax/ax.scale_y['scale'] - ax.scale_y['offset']) 

479 

480 if self.displacement_unit == 'm': 

481 def cfmt(x, p): 

482 return '%.2f' % x 

483 elif self.displacement_unit == 'cm': 

484 def cfmt(x, p): 

485 return '%.1f' % (x * 1e2) 

486 elif self.displacement_unit == 'mm': 

487 def cfmt(x, p): 

488 return '%.0f' % (x * 1e3) 

489 elif self.displacement_unit == 'rad': 

490 def cfmt(x, p): 

491 return '%.2f' % x 

492 else: 

493 raise AttributeError( 

494 'unknown displacement unit %s' % self.displacement_unit) 

495 

496 cbar_args = dict( 

497 orientation='horizontal', 

498 format=FuncFormatter(cfmt), 

499 use_gridspec=True) 

500 cbar_label = 'LOS Displacement [%s]' % self.displacement_unit 

501 

502 if self.common_color_scale: 

503 cax = fig.add_subplot(gs[1, 1]) 

504 cax.set_aspect(.05) 

505 cbar = fig.colorbar(cmw, cax=cax, **cbar_args) 

506 cbar.set_label(cbar_label) 

507 else: 

508 for idata, data in enumerate((stat_syn, stat_obs, res)): 

509 cax = fig.add_subplot(gs[1, idata]) 

510 cax.set_aspect(.05) 

511 

512 if not self.displacement_unit == 'rad': 

513 abs_displ = num.abs(data).max() 

514 cmw.set_clim(-abs_displ, abs_displ) 

515 

516 cbar = fig.colorbar(cmw, cax=cax, **cbar_args) 

517 cbar.set_label(cbar_label) 

518 

519 return (item, fig) 

520 

521 for ifig, (sat_target, result) in enumerate(zip(sat_targets, results)): 

522 yield generate_plot(sat_target, result, ifig) 

523 

524 

525class SatelliteTargetDisplacementCloseup(SatelliteTargetDisplacement): 

526 ''' Close-up of satellite surface displacements and modelled data. ''' 

527 name = 'satellite_closeup' 

528 

529 map_scale = Float.T( 

530 default=2., 

531 help='Scale the map surroundings, larger value zooms out.') 

532 

533 def make(self, environ): 

534 cm = environ.get_plot_collection_manager() 

535 history = environ.get_history(subset='harvest') 

536 optimiser = environ.get_optimiser() 

537 ds = environ.get_dataset() 

538 

539 environ.setup_modelling() 

540 

541 cm.create_group_mpl( 

542 self, 

543 self.draw_static_fits(ds, history, optimiser, closeup=True), 

544 title=u'InSAR Displacements (Closeup)', 

545 section='fits', 

546 feather_icon='zoom-in', 

547 description=u''' 

548Maps showing subsampled surface displacements as observed, modelled and the 

549residual (observed minus modelled). 

550 

551The displacement values predicted by the orbit-ambiguity ramps are added to the 

552modelled displacements (middle panels). The color shows the LOS displacement 

553values associated with, and the extent of, every quadtree box. The light grey 

554dots show the focal point of pixels combined in the quadtree box. This point 

555corresponds to the position of the modelled data point. 

556 

557The large dark grey dot shows the reference source position. The grey filled 

558box shows the surface projection of the modelled source, with the thick-lined 

559edge marking the upper fault edge. Map is focused around the fault's extent. 

560''') 

561 

562 

563def get_plot_classes(): 

564 return [SatelliteTargetDisplacement, SatelliteTargetDisplacementCloseup]