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.')
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()
129 environ.setup_modelling()
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).
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.
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''')
152 def draw_static_fits(self, ds, history, optimiser, closeup=False):
153 from pyrocko.orthodrome import latlon_to_ne_numpy
154 problem = history.problem
156 sat_targets = problem.satellite_targets
157 for target in sat_targets:
158 target.set_dataset(ds)
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()
167 results = problem.evaluate(model, targets=sat_targets)
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)
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
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')
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
198 ax.set_aspect(1./num.cos(source.effective_lat*d2r))
200 if self.relative_coordinates:
201 scale_x['offset'] = 0.
202 scale_y['offset'] = 0.
204 nticks_x = 4 if abs(scene.frame.llLon) >= 100 else 5
206 ax.xaxis.set_major_locator(MaxNLocator(self.nticks_x or nticks_x))
207 ax.yaxis.set_major_locator(MaxNLocator(5))
209 ax.scale_x = scale_x
210 ax.scale_y = scale_y
212 scale_axes(ax.get_xaxis(), **scale_x)
213 scale_axes(ax.get_yaxis(), **scale_y)
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
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)
233 def get_displacement_rgba(displacements, scene, mappable):
234 arr = num.full_like(scene.displacement, fill_value=num.nan)
235 qt = scene.quadtree
237 for syn_v, leaf in zip(displacements, qt.leaves):
238 arr[leaf._slice_rows, leaf._slice_cols] = syn_v
240 arr[scene.displacement_mask] = num.nan
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)
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.
258 return rgb_arr
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)
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)
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
280 az_dx = num.cos(phi - num.pi/2) * .125
281 az_dy = num.sin(phi - num.pi/2) * .125
283 anchor_x = .9 if los_dx < 0 else .1
284 anchor_y = .85 if los_dx < 0 else .975
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)
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)
304 ax.add_artist(az_arrow)
305 ax.add_artist(los_arrow)
307 urE, urN, llE, llN = (0., 0., 0., 0.)
308 for target in sat_targets:
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
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])
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)))
330 def generate_plot(sat_target, result, ifig):
332 scene = sat_target.scene
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])
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''')
353 stat_obs = result.statics_obs['displacement.los']
354 stat_syn = result.statics_syn['displacement.los']
355 res = stat_obs - stat_syn
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
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)
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')
377 stat_obs = displ2rad(stat_obs, wavelength)
378 stat_syn = displ2rad(stat_syn, wavelength)
379 res = displ2rad(res, wavelength)
381 self.colormap = 'hsv'
382 data_range = (0., num.pi)
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)
390 cmw = cm.ScalarMappable(cmap=self.colormap)
391 cmw.set_clim(*data_range)
392 cmw.set_array(stat_obs)
394 axes = [fig.add_subplot(gs[0, 0]),
395 fig.add_subplot(gs[0, 1]),
396 fig.add_subplot(gs[0, 2])]
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')
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)
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)
423 ax = axes[2]
424 ax.imshow(
425 get_displacement_rgba(res, scene, cmw),
426 extent=im_extent, origin='lower')
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)
434 for ax in axes:
435 ax.set_xlim(*im_extent[:2])
436 ax.set_ylim(*im_extent[2:])
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
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]
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
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)
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'
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'])
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)
493 cbar_args = dict(
494 orientation='horizontal',
495 format=FuncFormatter(cfmt),
496 use_gridspec=True)
497 cbar_label = 'LOS Displacement [%s]' % self.displacement_unit
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)
509 if not self.displacement_unit == 'rad':
510 abs_displ = num.abs(data).max()
511 cmw.set_clim(-abs_displ, abs_displ)
513 cbar = fig.colorbar(cmw, cax=cax, **cbar_args)
514 cbar.set_label(cbar_label)
516 return (item, fig)
518 for ifig, (sat_target, result) in enumerate(zip(sat_targets, results)):
519 yield generate_plot(sat_target, result, ifig)
522class SatelliteTargetDisplacementCloseup(SatelliteTargetDisplacement):
523 ''' Close-up of satellite surface displacements and modelled data. '''
524 name = 'satellite_closeup'
526 map_scale = Float.T(
527 default=2.,
528 help='Scale the map surroundings, larger value zooms out.')
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()
536 environ.setup_modelling()
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).
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.
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''')
560def get_plot_classes():
561 return [SatelliteTargetDisplacement, SatelliteTargetDisplacementCloseup]