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

1import logging

2import numpy as num

3from matplotlib import cm, gridspec

5from grond.plot.config import PlotConfig

6from grond.plot.collection import PlotItem

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

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

15km = 1e3

16d2r = num.pi/180.

17guts_prefix = 'grond'

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

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

30 # convolution of two 2D arrays

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

32 shad *= -1.

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

42 # normalize shading

43 shad -= num.nanmin(shad)

44 shad /= num.nanmax(shad)

46 if mask is not None:

47 shad[mask] = num.nan

49 # reduce range to balance gray color

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

51 shad += shad_lim[0]

53 rgb_map = mappable.to_rgba(displacement)

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

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

57 return rgb_map

60def displ2rad(displ, wavelength):

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

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

65 from matplotlib.ticker import ScalarFormatter

67 class FormatScaled(ScalarFormatter):

69 @staticmethod

70 def __call__(value, pos):

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

72 .replace(',', ' ')

74 axis.set_major_formatter(FormatScaled())

77class SatelliteTargetDisplacement(PlotConfig):

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

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

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]