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

298 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-27 13:27 +0000

1import logging 

2import numpy as num 

3from matplotlib import cm, gridspec 

4 

5from grond.plot.config import PlotConfig 

6from grond.plot.collection import PlotItem 

7 

8from matplotlib import pyplot as plt 

9from matplotlib.ticker import MaxNLocator, FuncFormatter 

10from matplotlib import patches 

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

12 

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

14 

15km = 1e3 

16d2r = num.pi/180. 

17guts_prefix = 'grond' 

18 

19 

20def drape_displacements( 

21 displacement, shad_data, mappable, 

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

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

24 

25 from scipy.ndimage import convolve as im_conv 

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

27 # from upper left 

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

29 

30 # convolution of two 2D arrays 

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

32 shad *= -1. 

33 

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

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

36 # shades helps 

37 percentile2 = num.quantile(shad, 0.02) 

38 percentile98 = num.quantile(shad, 0.98) 

39 shad[shad > percentile98] = percentile98 

40 shad[shad < percentile2] = percentile2 

41 

42 # normalize shading 

43 shad -= num.nanmin(shad) 

44 shad /= num.nanmax(shad) 

45 

46 if mask is not None: 

47 shad[mask] = num.nan 

48 

49 # reduce range to balance gray color 

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

51 shad += shad_lim[0] 

52 

53 rgb_map = mappable.to_rgba(displacement) 

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

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

56 

57 return rgb_map 

58 

59 

60def displ2rad(displ, wavelength): 

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

62 

63 

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

65 from matplotlib.ticker import ScalarFormatter 

66 

67 class FormatScaled(ScalarFormatter): 

68 

69 @staticmethod 

70 def __call__(value, pos): 

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

72 .replace(',', ' ') 

73 

74 axis.set_major_formatter(FormatScaled()) 

75 

76 

77class SatelliteTargetDisplacement(PlotConfig): 

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

79 

80 name = 'satellite' 

81 dpi = Int.T( 

82 default=250) 

83 size_cm = Tuple.T( 

84 2, Float.T(), 

85 default=(22., 12.)) 

86 colormap = String.T( 

87 default='RdBu', 

88 help='Colormap for the surface displacements') 

89 relative_coordinates = Bool.T( 

90 default=False, 

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

92 fit = StringChoice.T( 

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

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

95 ' ensamble.') 

96 

97 show_topo = Bool.T( 

98 default=True, 

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

100 displacement_unit = StringChoice.T( 

101 default='m', 

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

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

104 show_leaf_centres = Bool.T( 

105 default=True, 

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

107 source_outline_color = String.T( 

108 default='grey', 

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

110 common_color_scale = Bool.T( 

111 default=True, 

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

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

114 map_limits = Tuple.T( 

115 4, Float.T(), 

116 optional=True, 

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

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

119 nticks_x = Int.T( 

120 optional=True, 

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

122 

123 def make(self, environ): 

124 cm = environ.get_plot_collection_manager() 

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

126 optimiser = environ.get_optimiser() 

127 ds = environ.get_dataset() 

128 

129 environ.setup_modelling() 

130 

131 cm.create_group_mpl( 

132 self, 

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

134 title=u'InSAR Displacements', 

135 section='fits', 

136 feather_icon='navigation', 

137 description=u''' 

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

139residual (observed minus modelled). 

140 

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

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

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

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

145corresponds to the position of the modelled data point. 

146 

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

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

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

150''') 

151 

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

153 from pyrocko.orthodrome import latlon_to_ne_numpy 

154 problem = history.problem 

155 

156 sat_targets = problem.satellite_targets 

157 for target in sat_targets: 

158 target.set_dataset(ds) 

159 

160 if self.fit == 'best': 

161 source = history.get_best_source() 

162 model = history.get_best_model() 

163 elif self.fit == 'mean': 

164 source = history.get_mean_source() 

165 model = history.get_mean_model() 

166 

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

168 

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

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

171 ax.tick_params(length=2) 

172 

173 if scene.frame.isMeter(): 

174 import utm 

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

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

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

178 utm_E, utm_N, utm_zone, utm_zone_letter =\ 

179 utm.from_latlon(source.effective_lat, 

180 source.effective_lon) 

181 scale_x['offset'] = utm_E 

182 scale_y['offset'] = utm_N 

183 

184 if last_axes: 

185 ax.text(0.975, 0.025, 

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

187 va='bottom', ha='right', 

188 fontsize=8, alpha=.7, 

189 transform=ax.transAxes) 

190 ax.set_aspect('equal') 

191 

192 elif scene.frame.isDegree(): 

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

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

195 scale_x['offset'] = source.effective_lon 

196 scale_y['offset'] = source.effective_lat 

197 

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

199 

200 if self.relative_coordinates: 

201 scale_x['offset'] = 0. 

202 scale_y['offset'] = 0. 

203 

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

205 

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

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

208 

209 ax.scale_x = scale_x 

210 ax.scale_y = scale_y 

211 

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

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

214 

215 def draw_source(ax, scene): 

216 if scene.frame.isMeter(): 

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

218 fn -= fn.mean() 

219 fe -= fe.mean() 

220 elif scene.frame.isDegree(): 

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

222 fn -= source.effective_lat 

223 fe -= source.effective_lon 

224 

225 # source is centered 

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

227 ax.fill(fe, fn, 

228 edgecolor=(0., 0., 0.), 

229 facecolor=self.source_outline_color, 

230 alpha=0.7) 

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

232 

233 def get_displacement_rgba(displacements, scene, mappable): 

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

235 qt = scene.quadtree 

236 

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

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

239 

240 arr[scene.displacement_mask] = num.nan 

241 

242 if not self.common_color_scale \ 

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

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

245 mappable.set_clim(-abs_displ, abs_displ) 

246 

247 if self.show_topo: 

248 try: 

249 elevation = scene.get_elevation() 

250 return drape_displacements(arr, elevation, mappable) 

251 except Exception as e: 

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

253 logger.exception(e) 

254 rgb_arr = mappable.to_rgba(arr) 

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

256 rgb_arr[scene.displacement_mask] = 1. 

257 

258 return rgb_arr 

259 

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

261 rects = scene.quadtree.getMPLRectangles() 

262 for r in rects: 

263 r.set_edgecolor((.4, .4, .4)) 

264 r.set_linewidth(.5) 

265 r.set_facecolor('none') 

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

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

268 map(ax.add_artist, rects) 

269 

270 if self.show_leaf_centres: 

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

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

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

274 

275 def add_arrow(ax, scene): 

276 phi = num.nanmean(scene.phi) 

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

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

279 

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

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

282 

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

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

285 

286 az_arrow = patches.FancyArrow( 

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

288 dx=az_dx, dy=az_dy, 

289 head_width=.025, 

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

291 head_starts_at_zero=False, 

292 length_includes_head=True, 

293 transform=ax.transAxes) 

294 

295 los_arrow = patches.FancyArrow( 

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

297 dx=los_dx, dy=los_dy, 

298 head_width=.02, 

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

300 head_starts_at_zero=False, 

301 length_includes_head=True, 

302 transform=ax.transAxes) 

303 

304 ax.add_artist(az_arrow) 

305 ax.add_artist(los_arrow) 

306 

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

308 for target in sat_targets: 

309 

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

311 off_n, off_e = map(float, latlon_to_ne_numpy( 

312 target.scene.frame.llLat, target.scene.frame.llLon, 

313 source.effective_lat, source.effective_lon)) 

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

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

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

317 

318 turE, turN, tllE, tllN = zip( 

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

320 leaf.gridN.max()-off_n, 

321 leaf.gridE.min()-off_e, 

322 leaf.gridN.min()-off_n) 

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

324 

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

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

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

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

329 

330 def generate_plot(sat_target, result, ifig): 

331 

332 scene = sat_target.scene 

333 

334 fig = plt.figure() 

335 fig.set_size_inches(*self.size_inch) 

336 gs = gridspec.GridSpec( 

337 2, 3, 

338 wspace=.15, hspace=.2, 

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

340 height_ratios=[12, 1]) 

341 

342 item = PlotItem( 

343 name='fig_%i' % ifig, 

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

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

346 % scene.meta.scene_title, 

347 description=u''' 

348Surface displacements derived from satellite data. 

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

350data and (right) the model residual. 

351''') 

352 

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

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

355 res = stat_obs - stat_syn 

356 

357 if scene.frame.isMeter(): 

358 offset_n, offset_e = map(float, latlon_to_ne_numpy( 

359 scene.frame.llLat, scene.frame.llLon, 

360 source.effective_lat, source.effective_lon)) 

361 elif scene.frame.isDegree(): 

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

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

364 

365 im_extent = ( 

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

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

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

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

370 

371 if self.displacement_unit == 'rad': 

372 wavelength = scene.meta.wavelength 

373 if wavelength is None: 

374 raise AttributeError( 

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

376 

377 stat_obs = displ2rad(stat_obs, wavelength) 

378 stat_syn = displ2rad(stat_syn, wavelength) 

379 res = displ2rad(res, wavelength) 

380 

381 self.colormap = 'hsv' 

382 data_range = (0., num.pi) 

383 

384 else: 

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

386 stat_syn.min(), stat_syn.max(), 

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

388 data_range = (-abs_displ, abs_displ) 

389 

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

391 cmw.set_clim(*data_range) 

392 cmw.set_array(stat_obs) 

393 

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

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

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

397 

398 ax = axes[0] 

399 ax.imshow( 

400 get_displacement_rgba(stat_obs, scene, cmw), 

401 extent=im_extent, origin='lower') 

402 draw_leaves(ax, scene, offset_e, offset_n) 

403 draw_source(ax, scene) 

404 add_arrow(ax, scene) 

405 init_axes(ax, scene, 'Observed') 

406 

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

408 fontsize=8, alpha=.7, 

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

410 if scene.frame.isMeter(): 

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

412 

413 ax = axes[1] 

414 ax.imshow( 

415 get_displacement_rgba(stat_syn, scene, cmw), 

416 extent=im_extent, origin='lower') 

417 draw_leaves(ax, scene, offset_e, offset_n) 

418 draw_source(ax, scene) 

419 add_arrow(ax, scene) 

420 init_axes(ax, scene, 'Model') 

421 ax.get_yaxis().set_visible(False) 

422 

423 ax = axes[2] 

424 ax.imshow( 

425 get_displacement_rgba(res, scene, cmw), 

426 extent=im_extent, origin='lower') 

427 

428 draw_leaves(ax, scene, offset_e, offset_n) 

429 draw_source(ax, scene) 

430 add_arrow(ax, scene) 

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

432 ax.get_yaxis().set_visible(False) 

433 

434 for ax in axes: 

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

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

437 

438 if closeup: 

439 if scene.frame.isMeter(): 

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

441 elif scene.frame.isDegree(): 

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

443 fn -= source.effective_lat 

444 fe -= source.effective_lon 

445 

446 if fn.size > 1: 

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

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

449 else: 

450 off_n = fn[0] 

451 off_e = fe[0] 

452 

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

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

455 fault_size *= self.map_scale 

456 if fault_size == 0.0: 

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

458 fault_size = extent * .25 

459 

460 for ax in axes: 

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

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

463 

464 if self.map_limits is not None: 

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

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

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

468 

469 for ax in axes: 

470 ax.set_xlim( 

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

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

473 ax.set_ylim( 

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

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

476 

477 if self.displacement_unit == 'm': 

478 def cfmt(x, p): 

479 return '%.2f' % x 

480 elif self.displacement_unit == 'cm': 

481 def cfmt(x, p): 

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

483 elif self.displacement_unit == 'mm': 

484 def cfmt(x, p): 

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

486 elif self.displacement_unit == 'rad': 

487 def cfmt(x, p): 

488 return '%.2f' % x 

489 else: 

490 raise AttributeError( 

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

492 

493 cbar_args = dict( 

494 orientation='horizontal', 

495 format=FuncFormatter(cfmt), 

496 use_gridspec=True) 

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

498 

499 if self.common_color_scale: 

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

501 cax.set_aspect(.05) 

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

503 cbar.set_label(cbar_label) 

504 else: 

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

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

507 cax.set_aspect(.05) 

508 

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

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

511 cmw.set_clim(-abs_displ, abs_displ) 

512 

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

514 cbar.set_label(cbar_label) 

515 

516 return (item, fig) 

517 

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

519 yield generate_plot(sat_target, result, ifig) 

520 

521 

522class SatelliteTargetDisplacementCloseup(SatelliteTargetDisplacement): 

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

524 name = 'satellite_closeup' 

525 

526 map_scale = Float.T( 

527 default=2., 

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

529 

530 def make(self, environ): 

531 cm = environ.get_plot_collection_manager() 

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

533 optimiser = environ.get_optimiser() 

534 ds = environ.get_dataset() 

535 

536 environ.setup_modelling() 

537 

538 cm.create_group_mpl( 

539 self, 

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

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

542 section='fits', 

543 feather_icon='zoom-in', 

544 description=u''' 

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

546residual (observed minus modelled). 

547 

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

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

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

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

552corresponds to the position of the modelled data point. 

553 

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

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

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

557''') 

558 

559 

560def get_plot_classes(): 

561 return [SatelliteTargetDisplacement, SatelliteTargetDisplacementCloseup]