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, l in zip(displacements, qt.leaves):

238 arr[l._slice_rows, l._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 *[(leave.gridE.max()-off_e,

320 leave.gridN.max()-off_n,

321 leave.gridE.min()-off_e,

322 leave.gridN.min()-off_n)

323 for leave 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

467 for ax in axes:

468 ax.set_xlim(

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

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

471 ax.set_ylim(

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

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

474

475 if self.displacement_unit == 'm':

476 def cfmt(x, p):

477 return '%.2f' % x

478 elif self.displacement_unit == 'cm':

479 def cfmt(x, p):

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

481 elif self.displacement_unit == 'mm':

482 def cfmt(x, p):

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

484 elif self.displacement_unit == 'rad':

485 def cfmt(x, p):

486 return '%.2f' % x

487 else:

488 raise AttributeError(

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

490

491 cbar_args = dict(

492 orientation='horizontal',

493 format=FuncFormatter(cfmt),

494 use_gridspec=True)

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

496

497 if self.common_color_scale:

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

499 cax.set_aspect(.05)

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

501 cbar.set_label(cbar_label)

502 else:

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

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

505 cax.set_aspect(.05)

506

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

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

509 cmw.set_clim(-abs_displ, abs_displ)

510

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

512 cbar.set_label(cbar_label)

513

514 return (item, fig)

515

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

517 yield generate_plot(sat_target, result, ifig)

518

519

520class SatelliteTargetDisplacementCloseup(SatelliteTargetDisplacement):

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

522 name = 'satellite_closeup'

523

524 map_scale = Float.T(

525 default=2.,

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

527

528 def make(self, environ):

529 cm = environ.get_plot_collection_manager()

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

531 optimiser = environ.get_optimiser()

532 ds = environ.get_dataset()

533

534 environ.setup_modelling()

535

536 cm.create_group_mpl(

537 self,

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

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

540 section='fits',

541 feather_icon='zoom-in',

542 description=u'''

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

544residual (observed minus modelled).

545

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

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

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

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

550corresponds to the position of the modelled data point.

551

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

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

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

555''')

556

557

558def get_plot_classes():

559 return [SatelliteTargetDisplacement, SatelliteTargetDisplacementCloseup]