From b29cca73597eafc3f380d867471431ed9ab070b7 Mon Sep 17 00:00:00 2001 From: Yuan-Kai Liu <55262505+yuankailiu@users.noreply.github.com> Date: Sun, 24 Sep 2023 03:57:13 -0700 Subject: [PATCH] GMT-like shaded relief illumination via `view.py --dem-blend` (#1089) Rather than plotting DEM in the background and overlaying the data with transparency, I attempt to view data with color adjusted by the shaded relief DEM when provided. Just like `gmt grdimage -I` option. Hereafter, I call this DEM-blended data. This can be done via `matplotlib.colors.LightSource.shade_rgb()` module (as opposed to what has been used in MintPy, i.e., the `shade()` module, that only produces the shaded relief of the DEM itself). This PR adds this support for the single-subplot scenario via changes in `view.plot_slice()`, the multi-subplot scenarios are not supported. Detailed changes in MintPy code are as follows: + `utils.arg_utils.add_dem_argument()`: add the following new options to allow users to activate and tune the DEM-blended data plotting: `--dem-blend`, `--blend-mode`, `--shade-frac` and `--base-color` + `utils.plot.py`: - add `prep_blend_image()` to prepare the illuminated RGB array from a given data and DEM - add `plot_blend_image()` to call `prep_blend_image()` and `imshow` for plotting - add `blend_colorbar()` to create a DEM-blended colorbar, to be called within `plot_colorbar()` - `plot_colorbar()`: add generic customization support for the DEM-blended colorbar from `blend_colorbar()` - rename: `prepare_dem_background()` --> `prep_dem_background()` + `cli.view.cmd_line_parse()`: add checkings for `--dem-blend` option - turn off the disp_dem_shade option - check the existence of `--dem` option + `view.plot_slice()`: call `plot.plot_blend_image()` for the DEM-blended style. TODO: we could add the resampling capability in plot_blend_image(), to handle the different resolution between data and DEM. --------- Co-authored-by: Zhang Yunjun --- src/mintpy/cli/view.py | 6 + src/mintpy/utils/arg_utils.py | 40 +++-- src/mintpy/utils/plot.py | 304 +++++++++++++++++++++++++++++----- src/mintpy/view.py | 65 ++++---- 4 files changed, 333 insertions(+), 82 deletions(-) diff --git a/src/mintpy/cli/view.py b/src/mintpy/cli/view.py index 539530109..1a3306afe 100755 --- a/src/mintpy/cli/view.py +++ b/src/mintpy/cli/view.py @@ -143,6 +143,12 @@ def cmd_line_parse(iargs=None): if inps.flip_lr or inps.flip_ud: inps.auto_flip = False + if inps.disp_dem_blend: + inps.disp_dem_shade = False + # --dem-blend option requires --dem option + if inps.dem_file is None: + parser.error("--dem-blend requires -d/-dem.") + # check: conflicted options (geo-only options if inpput file is in radar-coordinates) geo_opt_names = ['--coord', '--show-gps', '--coastline', '--lalo-label', '--lalo-step', '--scalebar', '--faultline'] geo_opt_names = list(set(geo_opt_names) & set(inps.argv)) diff --git a/src/mintpy/utils/arg_utils.py b/src/mintpy/utils/arg_utils.py index 64a1f9979..3ef524a07 100644 --- a/src/mintpy/utils/arg_utils.py +++ b/src/mintpy/utils/arg_utils.py @@ -104,29 +104,43 @@ def add_dem_argument(parser): help='do not show DEM shaded relief') dem.add_argument('--dem-nocontour', dest='disp_dem_contour', action='store_false', help='do not show DEM contour lines') + dem.add_argument('--dem-blend', dest='disp_dem_blend', action='store_true', + help='blend the DEM shade with input image to have a GMT-like impression.') - dem.add_argument('--contour-smooth', dest='dem_contour_smooth', type=float, default=3.0, - help='Background topography contour smooth factor - sigma of Gaussian filter. \n' + # DEM contours + dem.add_argument('--contour-smooth', dest='dem_contour_smooth', type=float, default=3.0, metavar='NUM', + help='[Contour] Topography contour smooth factor - sigma of Gaussian filter. \n' 'Set to 0.0 for no smoothing; (default: %(default)s).') dem.add_argument('--contour-step', dest='dem_contour_step', metavar='NUM', type=float, default=200.0, - help='Background topography contour step in meters (default: %(default)s).') + help='[Contour] Topography contour step in meters (default: %(default)s).') dem.add_argument('--contour-lw','--contour-linewidth', dest='dem_contour_linewidth', metavar='NUM', type=float, default=0.5, - help='Background topography contour linewidth (default: %(default)s).') + help='[Contour] Topography contour linewidth (default: %(default)s).') + # DEM shade dem.add_argument('--shade-az', dest='shade_azdeg', type=float, default=315., metavar='DEG', - help='The azimuth (0-360, degrees clockwise from North) of the light source ' + help='[Shade] Azimuth angle (0-360, degrees clockwise from North) of the light source ' '(default: %(default)s).') dem.add_argument('--shade-alt', dest='shade_altdeg', type=float, default=45., metavar='DEG', - help='The altitude (0-90, degrees up from horizontal) of the light source ' + help='[Shade] Altitude (0-90, degrees up from horizontal) of the light source ' '(default: %(default)s).') - dem.add_argument('--shade-min', dest='shade_min', type=float, default=-4000., metavar='MIN', - help='The min height in m of colormap of shaded relief topography (default: %(default)s).') + help='[Shade] Minimum height of shaded relief topography (default: %(default)s m).') dem.add_argument('--shade-max', dest='shade_max', type=float, default=999., metavar='MAX', - help='The max height of colormap of shaded relief topography (default: max(DEM)+2000).') - dem.add_argument('--shade-exag', dest='shade_exag', type=float, default=0.5, - help='Vertical exaggeration ratio (default: %(default)s).') + help='[Shade] Maximum height of shaded relief topography (default: max(DEM)+2000 m).') + dem.add_argument('--shade-exag', dest='shade_exag', type=float, default=0.5, metavar='NUM', + help='[Shade] Vertical exaggeration ratio (default: %(default)s).') + + # DEM-blended image + dem.add_argument('--shade-frac', dest='shade_frac', type=float, default=0.5, metavar='NUM', + help='[Blend] Increases/decreases the contrast of the hillshade (default: %(default)s).') + dem.add_argument('--base-color', dest='base_color', type=float, default=0.9, metavar='NUM', + help='[Blend] Topograhpy basemap greyish color ranges in [0,1] (default: %(default)s).') + dem.add_argument('--blend-mode', dest='blend_mode', type=str, default='overlay', + choices={'hsv','overlay','soft'}, metavar='STR', + help='[Blend] Type of blending used to combine the colormapped data with illumated ' + 'topography.\n(choices: %(choices)s; default: %(default)s).\n' + 'https://matplotlib.org/stable/gallery/specialty_plots/topographic_hillshading.html') return parser @@ -277,8 +291,8 @@ def add_map_argument(parser): mapg.add_argument('--coastline', dest='coastline', type=str, choices={'10m', '50m', '110m'}, help="Draw coastline with specified resolution (default: %(default)s).\n" "This will enable --lalo-label option.\n" - "Link: https://scitools.org.uk/cartopy/docs/latest/matplotlib/geoaxes.html" - "#cartopy.mpl.geoaxes.GeoAxes.coastlines") + "https://scitools.org.uk/cartopy/docs/latest/reference/generated/" + "cartopy.mpl.geoaxes.GeoAxes.html") mapg.add_argument('--coastline-lw', '--coastline-linewidth', dest='coastline_linewidth', metavar='NUM', type=float, default=1, help='Coastline linewidth (default: %(default)s).') diff --git a/src/mintpy/utils/plot.py b/src/mintpy/utils/plot.py index 59cdbfd00..99656e4cb 100644 --- a/src/mintpy/utils/plot.py +++ b/src/mintpy/utils/plot.py @@ -1380,59 +1380,98 @@ def plot_insar_vs_gps_scatter(vel_file, csv_file='gps_enu2los.csv', msk_file=Non def plot_colorbar(inps, im, cax): + # orientation + if inps.cbar_loc in ['left', 'right']: + orientation = 'vertical' + else: + orientation = 'horizontal' + # expand vlim by 0.01% to account for potential numerical precision leak # e.g. wrapped phase epsilon = (inps.vlim[1] - inps.vlim[0]) * 0.0001 vmin = inps.vlim[0] - epsilon vmax = inps.vlim[1] + epsilon - # extend + # extend type if not inps.cbar_ext: if vmin <= inps.dlim[0] and vmax >= inps.dlim[1]: inps.cbar_ext='neither' elif vmin > inps.dlim[0] and vmax >= inps.dlim[1]: inps.cbar_ext='min' elif vmin <= inps.dlim[0] and vmax < inps.dlim[1]: inps.cbar_ext='max' else: inps.cbar_ext='both' - # orientation - if inps.cbar_loc in ['left', 'right']: - orientation = 'vertical' + # ticks for special cases + if abs(vmin + np.pi) / np.pi < 0.001 and abs(vmax - np.pi) / np.pi < 0.001: + ticks = [-np.pi, 0, np.pi] # special case 1: -pi/pi + elif hasattr(inps, 'unique_values') and inps.unique_values is not None and len(inps.unique_values) <= 5: + ticks = list(inps.unique_values) # special case 2: show finite exact tick values else: - orientation = 'horizontal' + ticks = None # plot colorbar - unique_values = getattr(inps, 'unique_values', None) - if abs(vmin + np.pi) / np.pi < 0.001 and abs(vmax - np.pi) / np.pi < 0.001: - cbar = plt.colorbar(im, cax=cax, orientation=orientation, ticks=[-np.pi, 0, np.pi]) - if orientation == 'vertical': - cbar.ax.set_yticklabels([r'-$\pi$', '0', r'$\pi$']) - else: - cbar.ax.set_xticklabels([r'-$\pi$', '0', r'$\pi$']) - - elif unique_values is not None and len(inps.unique_values) <= 5: - # show the exact tick values - cbar = plt.colorbar(im, cax=cax, orientation=orientation, ticks=inps.unique_values) - + if not inps.disp_dem_blend: + # regular colorbar + cbar = plt.colorbar(im, cax=cax, orientation=orientation, extend=inps.cbar_ext, ticks=ticks) + cbar_type = 'mpl' else: - cbar = plt.colorbar(im, cax=cax, orientation=orientation, extend=inps.cbar_ext) + # illuminated colorbar for DEM-blended images + blend_colorbar(cax, inps, vlim=[vmin, vmax], orientation=orientation, ticks=ticks) + cbar_type = 'img' + cbar = None + # ticks for generic cases if inps.cbar_nbins: if inps.cbar_nbins <= 2: # manually set tick for better positions when the color step is not a common number # e.g. for numInvIfgram.h5 - cbar.set_ticks(inps.dlim) + if cbar_type == 'mpl': + cbar.set_ticks(inps.dlim) + elif orientation == 'vertical': + cax.set_yticks(inps.dlim) + elif orientation == 'horizontal': + cax.set_xticks(inps.dlim) + else: - cbar.locator = ticker.MaxNLocator(nbins=inps.cbar_nbins) - cbar.update_ticks() + if cbar_type == 'mpl': + cbar.locator = ticker.MaxNLocator(nbins=inps.cbar_nbins) + cbar.update_ticks() + elif orientation == 'vertical': + cax.yaxis.set_major_locator(ticker.MaxNLocator(nbins=inps.cbar_nbins)) + elif orientation == 'horizontal': + cax.xaxis.set_major_locator(ticker.MaxNLocator(nbins=inps.cbar_nbins)) + elif inps.cbar_ticks: - cbar.set_ticks(inps.cbar_ticks) + if cbar_type == 'mpl': + cbar.set_ticks(inps.cbar_ticks) + elif orientation == 'vertical': + cax.set_yticks(inps.cbar_ticks) + elif orientation == 'horizontal': + cax.set_xticks(inps.cbar_ticks) + + # update tick labels for special symbol: pi + if ticks and len(ticks) == 3 and ticks == [-np.pi, 0, np.pi]: + if orientation == 'vertical': + cax.set_yticklabels([r'-$\pi$', '0', r'$\pi$']) + else: + cax.set_xticklabels([r'-$\pi$', '0', r'$\pi$']) - cbar.ax.tick_params(which='both', direction='out', labelsize=inps.font_size, colors=inps.font_color) + cax.tick_params(which='both', direction='out', labelsize=inps.font_size, colors=inps.font_color) - # add label + # colorbar label if inps.cbar_label: - cbar.set_label(inps.cbar_label, fontsize=inps.font_size, color=inps.font_color) + cbar_label = inps.cbar_label elif inps.disp_unit != '1': - cbar.set_label(inps.disp_unit, fontsize=inps.font_size, color=inps.font_color) + cbar_label = inps.disp_unit + else: + cbar_label = None + + if cbar_label is not None: + kwargs = dict(fontsize=inps.font_size, color=inps.font_color) + if cbar_type == 'mpl': + cbar.set_label(cbar_label, **kwargs) + elif orientation == 'vertical': + cax.set_ylabel(cbar_label, **kwargs) + elif orientation == 'horizontal': + cax.set_xlabel(cbar_label, **kwargs) return inps, cbar @@ -1877,7 +1916,7 @@ def read_dem(dem_file, pix_box=None, geo_box=None, print_msg=True): return dem, dem_meta, dem_pix_box -def prepare_dem_background(dem, inps, print_msg=True): +def prep_dem_background(dem, inps, print_msg=True): """Prepare to plot DEM on background Parameters: dem : 2D np.int16 matrix, dem data inps : Namespace with the following 4 items: @@ -1887,14 +1926,14 @@ def prepare_dem_background(dem, inps, print_msg=True): 'dem_contour_smooth': float, 3.0 Returns: dem_shade : 3D np.array in size of (length, width, 4) dem_contour : 2D np.array in size of (length, width) - dem_contour_sequence : 1D np.array + dem_contour_seq : 1D np.array Examples: from mintpy.cli import view from mintpy.utils import plot as pp inps = view.cmd_line_parse() dem = readfile.read('inputs/geometryRadar.h5')[0] - dem_shade, dem_contour, dem_contour_seq = pp.prepare_dem_background( + dem_shade, dem_contour, dem_contour_seq = pp.prep_dem_background( dem=dem, inps=inps, ) @@ -1902,7 +1941,7 @@ def prepare_dem_background(dem, inps, print_msg=True): # default returns dem_shade = None dem_contour = None - dem_contour_sequence = None + dem_contour_seq = None # default inputs if inps.shade_max == 999.: @@ -1938,7 +1977,7 @@ def prepare_dem_background(dem, inps, print_msg=True): else: from scipy import ndimage dem_contour = ndimage.gaussian_filter(dem, sigma=inps.dem_contour_smooth, order=0) - dem_contour_sequence = np.arange(inps.dem_contour_step, 9000, step=inps.dem_contour_step) + dem_contour_seq = np.arange(inps.dem_contour_step, 9000, step=inps.dem_contour_step) if print_msg: msg = f'show contour in step of {inps.dem_contour_step} m ' msg += f'with a smoothing factor of {inps.dem_contour_smooth}' @@ -1957,7 +1996,7 @@ def prepare_dem_background(dem, inps, print_msg=True): else: print('WARNING: DEM has different size than mask, ignore --mask-dem and continue.') - return dem_shade, dem_contour, dem_contour_sequence + return dem_shade, dem_contour, dem_contour_seq def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_contour_seq=None, @@ -1967,7 +2006,7 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ geo_box : tuple of 4 float in order of (E, N, W, S), geo bounding box dem_shade : 3D np.array in size of (length, width, 4) dem_contour : 2D np.array in size of (length, width) - dem_contour_sequence : 1D np.array + dem_contour_seq : 1D np.array dem : 2D np.array of DEM data inps : Namespace with the following 4 items: 'disp_dem_shade' : bool, True/False @@ -1985,7 +2024,7 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ # prepare DEM shade/contour datasets if all(i is None for i in [dem_shade, dem_contour, dem_contour_seq]) and dem is not None: - dem_shade, dem_contour, dem_contour_seq = prepare_dem_background( + dem_shade, dem_contour, dem_contour_seq = prep_dem_background( dem, inps=inps, print_msg=print_msg, @@ -2006,15 +2045,12 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ # plot shaded relief if dem_shade is not None: - # config kwargs = dict(interpolation='spline16', zorder=0, origin='upper') - - # geo coordinates if geo_box is not None: + # geo coordinates ax.imshow(dem_shade, extent=geo_extent, **kwargs) - - # radar coordinates elif isinstance(ax, plt.Axes): + # radar coordinates ax.imshow(dem_shade, extent=rdr_extent, **kwargs) # plot topo contour @@ -2039,3 +2075,193 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ ax.contour(dem_contour, dem_contour_seq, extent=rdr_extent, **kwargs) return ax + + +########################################## DEM-blended data ########################################### +def blend_colorbar(cax, inps, vlim, orientation='vertical', ticks=None, + fraction=0.75, blend_mode='soft', vert_exag=6000): + """Create a shade-illuminated colorbar. + + Parameters: cax - colorbar axis + inps - inps : Namespace with the following items: + 'shade_azdeg' : float, True/False + 'shade_alt' : float, 200.0 + 'colormap' : string or matplotlib.colors.colormap class + 'cbar_label' : string + vlim - list of 2 float, colorbar range + orientation - str, vertical or horizontal + ticks - list of float or None, colorbar ticks + fraction - float, increases or decreases the contrast of the hillshade + blend_mode - {'hsv', 'overlay', 'soft'} or callable + vert_exag - float, the amount to exaggerate the elevation values by + when calculating illumination + Examples: blend_colorbar(cax, inps, vlim=[-3, 6], orientation='vertical') + """ + from matplotlib.colors import LightSource + + # create normalized array for colorbar + nx = 1000 # array size along the illumination profile + ny = 255 # array size along the cmap + arr = np.tile(np.linspace(1, 0, ny).reshape(-1, 1), (1, nx)) + + # create an artificial topography for illumination + x = np.linspace(-np.pi/4, np.pi/4, nx) # x position along the illum. profile + elev = np.ones_like(arr) * np.cos(0.6*x) # altitude as a cosine along teh illum. profile + + if orientation == 'vertical': + extent = [0, nx, vlim[0], vlim[1]] # (left, right, bottom, top) + else: + extent = [vlim[0], vlim[1], 0, nx] # (left, right, bottom, top) + arr = arr.T + elev = elev.T + + # create illuminated RGB array + ls = LightSource(azdeg=inps.shade_azdeg, altdeg=inps.shade_altdeg) + mappable = plt.cm.ScalarMappable(norm=mpl.colors.Normalize(0,1), cmap=inps.colormap) + img_rgb = mappable.to_rgba(arr)[:, :, :3] + illum_rgb = ls.shade_rgb(img_rgb, elev, fraction=fraction, blend_mode=blend_mode, vert_exag=vert_exag) + + # plot colorbar as an image + cax.imshow(illum_rgb, extent=extent, aspect='auto') + + # axis format + tick_kwargs = dict(labelsize=inps.font_size, colors=inps.font_color) + label_kwargs = dict(fontsize=inps.font_size, color=inps.font_color) + if orientation == 'vertical': + cax.tick_params(which='both', bottom=False, top=False, labelbottom=False, **tick_kwargs) + cax.set_ylabel(inps.cbar_label, rotation=90, **label_kwargs) + cax.yaxis.set_label_position("right") + cax.yaxis.tick_right() + cax.set_xticks([]) + else: + cax.tick_params(which='both', left=False, right=False, labeltop=False, **tick_kwargs) + cax.set_xlabel(inps.cbar_label, **label_kwargs) + cax.xaxis.set_label_position("bottom") + cax.xaxis.tick_bottom() + cax.set_yticks([]) + + # update ticks if given + if ticks is not None: + if orientation == 'vertical': + cax.set_yticks(ticks) + else: + cax.set_xticks(ticks) + + return + + +def prep_blend_image(data, dem, vmin=None, vmax=None, cmap='viridis', + base_color=0.9, shade_frac=0.5, blend_mode='overlay', + azdeg=315, altdeg=45, vert_exag=0.5, mask_nan_dem=True, + mask_nan_data=False, fill_value=0): + """Prepare the illuminated RGB array for the data, using shaded relief DEM from a light source, + like the `gmt grdimage -I` feature, i.e. hillshade + DEM-blended data. + + Parameters: data - 2D np.ndarray in size of (m, n), data to be blended + dem - 2D np.ndarray in size of (m, n), dem data + vmin/max - float, lower/upper display limit of the data + cmap - str or matplotlib.colors.colormap class + base_color - float or color hex codes + shade_frac - float, increases or decreases the contrast of the hillshade + blend_mode - {'hsv', 'overlay', 'soft'} or callable + az/altdeg - float, azimuth/altitude angle of the light source + vert_exag - float, the amount to exaggerate the elevation values by + when calculating illumination + mask_nan_dem - bool, whether to mask blended image based on nan dem pixels + mask_nan_data - bool, whether to mask blended image based on nan data pixels + fill_value - float, set the masked pixels as alpha = fill_value (transparent) + Returns: illum_rgb - 3D np.ndarray of float32 in size of (m, n, 4), ranging between 0-1. + 1st to 3rd layers are the RGB values; 4th layer is the transparency + Examples: illum_rgb = pp.prep_blend_image(data, dem, vmin, vmax) + """ + from matplotlib.colors import LightSource + + # use numpy.ma to mask missing or invalid entries + data = np.ma.masked_invalid(data) + dem = np.ma.masked_invalid(dem) + + # data normalization + vmin = vmin if vmin else np.nanmin(data) + vmax = vmax if vmax else np.nanmax(data) + data_norm = (data - vmin) / (vmax - vmin) + + ## create data RGB array + # cmap norm and ScalarMappable + mappable = plt.cm.ScalarMappable(norm=mpl.colors.Normalize(0,1), cmap=cmap) + + # convert data norm to image and remove alpha channel (the fourth dimension) + img_rgb = mappable.to_rgba(data_norm)[:, :, :3] + + # assign a greyish basemap color to the masked pixels + img_rgb[data.mask, :] = base_color + + # check if the dimensions are coherent + # Note: in the future, resample one into the grid of the other one + if dem.shape != img_rgb.shape[:-1]: + raise ValueError(f'Dimension mismatch between DEM ({dem.shape}) and data ({data.shape})!') + + ## add shaded relief to illuminate the RGB array + ls = LightSource(azdeg=azdeg, altdeg=altdeg) + illum_rgb = ls.shade_rgb(img_rgb, dem, fraction=shade_frac, blend_mode=blend_mode, vert_exag=vert_exag) + + # add tranparency layer to the array (defualt: all ones = opaque) + illum_rgb = np.dstack([illum_rgb, np.ones_like(illum_rgb[:, :, 0])]) + + ## masking the shaded-relief image: + # - can set rgb (first 3 columns) to [0: black ,1: white, np.nan: transparent] + # - or can set the fourth column transparency to 0 (default) + if mask_nan_dem: + illum_rgb[dem.mask, -1] = fill_value + if mask_nan_data: + illum_rgb[data.mask, -1] = fill_value + + return illum_rgb + + +def plot_blend_image(ax, data, dem, inps, print_msg=True): + """Plot DEM-blended image. + + Parameters: ax - matplotlib.pyplot.Axes or BasemapExt object + data - 2D np.ndarray, image to be blended + dem - 2D np.ndarray, topography used for blending + inps - Namespace object with the following items: + 'base_color' : float + 'blend_mode' : str + 'colormap' : str or matplotlib.colors.colormap class + 'extent' : tuple of 4 float + 'shade_altdeg' : float + 'shade_azdeg' : float + 'shade_exag' : float + 'shade_frac' : float + 'mask_dem' : bool + 'vlim' : list of 2 float + print_msg - bool, print verbose message or not + Returns: im - matplotlob.pyplot.AxesImage + """ + if print_msg: + msg = 'plotting data ' + msg += f'blended by DEM shaded relief (contrast={inps.shade_frac:.1f}, ' + msg += f'base_color={inps.base_color:.1f}, exag={inps.shade_exag}, ' + msg += f'az/alt={inps.shade_azdeg}/{inps.shade_altdeg} deg) ...' + print(msg) + + # prepare + blend_img = prep_blend_image( + data, dem, + vmin=inps.vlim[0], + vmax=inps.vlim[1], + cmap=inps.colormap, + base_color=inps.base_color, + shade_frac=inps.shade_frac, + blend_mode=inps.blend_mode, + azdeg=inps.shade_azdeg, + altdeg=inps.shade_altdeg, + vert_exag=inps.shade_exag, + mask_nan_data=inps.mask_dem, + ) + + # plot + ax.imshow(blend_img, extent=inps.extent, interpolation='spline16', zorder=1, origin='upper') + im = plt.cm.ScalarMappable(norm=mpl.colors.Normalize(inps.vlim[0], inps.vlim[1]), cmap=inps.colormap) + + return im diff --git a/src/mintpy/view.py b/src/mintpy/view.py index 7b21e9439..e3ef94a0e 100644 --- a/src/mintpy/view.py +++ b/src/mintpy/view.py @@ -504,7 +504,7 @@ def plot_slice(ax, data, metadata, inps): vprint('plot in geo-coordinate') # extent info for matplotlib.imshow and other functions - extent = (inps.geo_box[0], inps.geo_box[2], inps.geo_box[3], inps.geo_box[1]) # (W, E, S, N) + inps.extent = (inps.geo_box[0], inps.geo_box[2], inps.geo_box[3], inps.geo_box[1]) # (W, E, S, N) SNWE = (inps.geo_box[3], inps.geo_box[1], inps.geo_box[0], inps.geo_box[2]) # Draw coastline using cartopy resolution parameters @@ -523,9 +523,8 @@ def plot_slice(ax, data, metadata, inps): print_msg=inps.print_msg, ) - # Plot Data + # Reference (InSAR) data to a GNSS site coord = ut.coordinate(metadata) - vprint('plotting image ...') if inps.disp_gps and inps.gps_component and inps.ref_gps_site: ref_site_lalo = GPS(site=inps.ref_gps_site).get_stat_lat_lon(print_msg=False) y, x = coord.geo2radar(ref_site_lalo[0], ref_site_lalo[1])[0:2] @@ -537,9 +536,14 @@ def plot_slice(ax, data, metadata, inps): # do not show the original InSAR reference point inps.disp_ref_pixel = False - im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], - extent=extent, origin='upper', interpolation=inps.interpolation, - alpha=inps.transparency, animated=inps.animation, zorder=1) + # Plot data + if inps.disp_dem_blend: + im = pp.plot_blend_image(ax, data, dem, inps, print_msg=inps.print_msg) + else: + vprint('plotting data ...') + im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], + extent=inps.extent, origin='upper', interpolation=inps.interpolation, + alpha=inps.transparency, animated=inps.animation, zorder=1) # Draw faultline using GMT lonlat file if inps.faultline_file: @@ -599,8 +603,6 @@ def plot_slice(ax, data, metadata, inps): # Show UNR GPS stations if inps.disp_gps: - SNWE = (inps.geo_box[3], inps.geo_box[1], - inps.geo_box[0], inps.geo_box[2]) ax = pp.plot_gps(ax, SNWE, inps, metadata, print_msg=inps.print_msg) # Status bar @@ -647,13 +649,18 @@ def format_coord(x, y): print_msg=inps.print_msg, ) - # Plot Data - vprint('plotting Data ...') # extent = (left, right, bottom, top) in data coordinates - extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, - inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) - im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], - extent=extent, interpolation=inps.interpolation, alpha=inps.transparency, zorder=1) + inps.extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, + inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) + + # Plot Data + if inps.disp_dem_blend: + im = pp.plot_blend_image(ax, data, dem, inps, print_msg=inps.print_msg) + else: + vprint('plotting data ...') + im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], + extent=inps.extent, interpolation=inps.interpolation, + alpha=inps.transparency, zorder=1) ax.tick_params(labelsize=inps.font_size) # Plot Seed Point @@ -688,8 +695,8 @@ def format_coord(x, y): ]) ax.plot(pts_yx[:, 1], pts_yx[:, 0], '-', ms=inps.ref_marker_size, mec='black', mew=1.) - ax.set_xlim(extent[0:2]) - ax.set_ylim(extent[2:4]) + ax.set_xlim(inps.extent[0:2]) + ax.set_ylim(inps.extent[2:4]) # Status bar @@ -761,7 +768,7 @@ def format_coord(x, y): if inps.disp_tick: # move x-axis tick label to the top if colorbar is at the bottom if inps.cbar_loc == 'bottom': - ax.tick_params(labelbottom=False, labeltop=True) + ax.tick_params(bottom=False, top=True, labelbottom=False, labeltop=True) # manually turn ON to enable tick labels for UTM with cartopy # link: https://github.com/SciTools/cartopy/issues/491 @@ -1203,14 +1210,12 @@ def plot_subplot4figure(i, inps, ax, data, metadata): print_msg=inps.print_msg) # Plot Data - vlim = inps.vlim - vlim = vlim if vlim is not None else [np.nanmin(data), np.nanmax(data)] - extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, - inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) - + vlim = inps.vlim if inps.vlim is not None else [np.nanmin(data), np.nanmax(data)] + inps.extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, + inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) im = ax.imshow(data, cmap=inps.colormap, vmin=vlim[0], vmax=vlim[1], interpolation=inps.interpolation, alpha=inps.transparency, - extent=extent, zorder=1) + extent=inps.extent, zorder=1) # Plot Seed Point if inps.disp_ref_pixel: @@ -1223,11 +1228,11 @@ def plot_subplot4figure(i, inps, ax, data, metadata): if ref_y and ref_x: ax.plot(ref_x, ref_y, inps.ref_marker, ms=inps.ref_marker_size) - ax.set_xlim(extent[0:2]) - ax.set_ylim(extent[2:4]) + ax.set_xlim(inps.extent[0:2]) + ax.set_ylim(inps.extent[2:4]) - # Subplot Setting - ## Tick and Label + ## Subplot Setting + # Tick and Label if not inps.disp_tick or inps.fig_row_num * inps.fig_col_num > 10: ax.get_xaxis().set_ticks([]) ax.get_yaxis().set_ticks([]) @@ -1395,11 +1400,11 @@ def adjust_subplots_layout(fig, inps): vprint('Note: different color scale for EACH subplot!') vprint('Adjust figsize for the colorbar of each subplot.') fig.set_size_inches(inps.fig_size[0] * 1.1, inps.fig_size[1]) - adjust_subplots_layout(fig, inps) + else: adjust_subplots_layout(fig, inps) - + # plot common colorbar cbar_length = 0.4 if inps.fig_size[1] > 8.0: cbar_length /= 2 @@ -1494,7 +1499,7 @@ def prepare4multi_subplots(inps, metadata): print_msg=False, )[0] - inps.dem_shade, inps.dem_contour, inps.dem_contour_seq = pp.prepare_dem_background( + inps.dem_shade, inps.dem_contour, inps.dem_contour_seq = pp.prep_dem_background( dem=dem, inps=inps, print_msg=inps.print_msg,