diff --git a/src/mintpy/utils/arg_utils.py b/src/mintpy/utils/arg_utils.py index 64a1f9979..70655e6e2 100644 --- a/src/mintpy/utils/arg_utils.py +++ b/src/mintpy/utils/arg_utils.py @@ -102,6 +102,8 @@ def add_dem_argument(parser): help='Mask out DEM pixels not coincident with valid data pixels') dem.add_argument('--dem-noshade', dest='disp_dem_shade', action='store_false', help='do not show DEM shaded relief') + 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('--dem-nocontour', dest='disp_dem_contour', action='store_false', help='do not show DEM contour lines') @@ -127,6 +129,10 @@ def add_dem_argument(parser): 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).') + dem.add_argument('--shade-frac', dest='shade_frac', type=float, default=0.5, + help='Only for --dem-blend. Increases/decreases the contrast of the hillshade (default: %(default)s).') + dem.add_argument('--base-color', dest='base_color', type=float, default=0.9, + help='Only for --dem-blend. DEM basemap greyish color ranges in [0,1] (default: %(default)s).') return parser diff --git a/src/mintpy/utils/plot.py b/src/mintpy/utils/plot.py index 59cdbfd00..cb7e117b5 100644 --- a/src/mintpy/utils/plot.py +++ b/src/mintpy/utils/plot.py @@ -14,8 +14,10 @@ import h5py import matplotlib as mpl +import matplotlib.colors as colors import numpy as np from matplotlib import dates as mdates, pyplot as plt, ticker +from matplotlib.colors import LightSource from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy import stats @@ -1877,17 +1879,20 @@ 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 prepare_dem_background(dem, inps, data=None, print_msg=True): """Prepare to plot DEM on background Parameters: dem : 2D np.int16 matrix, dem data - inps : Namespace with the following 4 items: + inps : Namespace with the following 5 items: 'disp_dem_shade' : bool, True/False 'disp_dem_contour' : bool, True/False 'dem_contour_step' : float, 200.0 'dem_contour_smooth': float, 3.0 + 'disp_dem_blend' : bool, True/False + data : 2D matrix, pass to shaded_image() if disp_dem_blend is True 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 + blend_img : 3D np.array of DEM-blended image rgbt Examples: from mintpy.cli import view from mintpy.utils import plot as pp @@ -1910,7 +1915,6 @@ def prepare_dem_background(dem, inps, print_msg=True): # prepare shade relief if inps.disp_dem_shade: - from matplotlib.colors import LightSource ls = LightSource(azdeg=inps.shade_azdeg, altdeg=inps.shade_altdeg) with warnings.catch_warnings(): warnings.simplefilter("ignore", category=RuntimeWarning) @@ -1922,8 +1926,24 @@ def prepare_dem_background(dem, inps, print_msg=True): vmax=inps.shade_max, ) dem_shade[np.isnan(dem_shade[:, :, 0])] = np.nan + + # fancier blended image: GMT-like data illuminated/blended with the shaded relief + blend_img = None + if inps.disp_dem_blend and not isinstance(data, type(None)): + vlim = inps.vlim + vlim = vlim if vlim is not None else [np.nanmin(data), np.nanmax(data)] + blend_img = shaded_image( + data, dem, ls, + vmin=vlim[0], vmax=vlim[1], + cmap=inps.colormap, base_color=inps.base_color, + shade_frac=inps.shade_frac, vert_exag=inps.shade_exag, + mask_dem_dataNan=inps.mask_dem) + if print_msg: - msg = f'show shaded relief DEM (min/max={inps.shade_min:.0f}/{inps.shade_max:.0f} m; ' + if inps.disp_dem_blend: + msg = f'show shaded relief DEM blended with data array (contrast={inps.shade_frac:.1f} base_color={inps.base_color:.1f}; ' + else: + msg = f'show shaded relief DEM (min/max={inps.shade_min:.0f}/{inps.shade_max:.0f} m; ' msg += f'exag={inps.shade_exag}; az/alt={inps.shade_azdeg}/{inps.shade_altdeg} deg)' print(msg) @@ -1957,11 +1977,11 @@ 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_sequence, blend_img def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_contour_seq=None, - dem=None, inps=None, print_msg=True): + dem=None, data=None, blend_img=None, inps=None, print_msg=True): """Plot DEM as background. Parameters: ax : matplotlib.pyplot.Axes or BasemapExt object geo_box : tuple of 4 float in order of (E, N, W, S), geo bounding box @@ -1969,11 +1989,14 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ dem_contour : 2D np.array in size of (length, width) dem_contour_sequence : 1D np.array dem : 2D np.array of DEM data - inps : Namespace with the following 4 items: + data : 2D np.array of data, for prepare_dem_background() if needed + blend_img : 3D np.array of DEM-blended image rgbt + inps : Namespace with the following 6 items: 'disp_dem_shade' : bool, True/False 'disp_dem_contour' : bool, True/False 'dem_contour_step' : float, 200.0 'dem_contour_smooth': float, 3.0 + 'disp_dem_blend' : bool, True/False 'pix_box' : 4-tuple of int, (x0, y0, x1, y1) Returns: ax : matplotlib.pyplot.Axes or BasemapExt object Examples: ax = pp.plot_dem_background(ax, geo_box=inps.geo_box, dem=dem, inps=inps) @@ -1984,10 +2007,11 @@ 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( + if any(i is None for i in [dem_shade, dem_contour, dem_contour_seq, blend_img]) and dem is not None: + dem_shade, dem_contour, dem_contour_seq, blend_img = prepare_dem_background( dem, inps=inps, + data=data, print_msg=print_msg, ) @@ -2005,17 +2029,31 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ pix_box[3]-0.5, pix_box[1]-0.5) # plot shaded relief + im = None if dem_shade is not None: # config kwargs = dict(interpolation='spline16', zorder=0, origin='upper') - # geo coordinates - if geo_box is not None: - ax.imshow(dem_shade, extent=geo_extent, **kwargs) + # style 1 (default style) + if blend_img is None: + # geo coordinates + if geo_box is not None: + ax.imshow(dem_shade, extent=geo_extent, **kwargs) + # radar coordinates + elif isinstance(ax, plt.Axes): + ax.imshow(dem_shade, extent=rdr_extent, **kwargs) - # radar coordinates - elif isinstance(ax, plt.Axes): - ax.imshow(dem_shade, extent=rdr_extent, **kwargs) + # style 2 (fancy GMT syle) + else: + vlim = inps.vlim + vlim = vlim if vlim is not None else [np.nanmin(data), np.nanmax(data)] + # geo coordinates + if geo_box is not None: + ax.imshow(blend_img, extent=geo_extent, **kwargs) + # radar coordinates + elif isinstance(ax, plt.Axes): + ax.imshow(blend_img, extent=rdr_extent, **kwargs) + im = plt.cm.ScalarMappable(norm=colors.Normalize(vlim[0], vlim[1]), cmap=inps.colormap) # plot topo contour if dem_contour is not None and dem_contour_seq is not None: @@ -2038,4 +2076,192 @@ def plot_dem_background(ax, geo_box=None, dem_shade=None, dem_contour=None, dem_ elif isinstance(ax, plt.Axes): ax.contour(dem_contour, dem_contour_seq, extent=rdr_extent, **kwargs) - return ax + return ax, im + + +########################################## Plot slice with DEM ########################################### +def shaded_colorbar(inps, cax, fraction=0.75, blend_mode='soft', vert_exag=6000): + """Create a shaded illuminated colorbar + Parameters : inps inps : Namespace with the following 5 items: + 'vlim' : list [float, float] display data limit + 'shade_azdeg' : float, True/False + 'shade_alt' : float, 200.0 + 'colormap' : string or matplotlib.colors.colormap class + 'cbar_label' : string + cax colorbar axis + 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 + + Returns : cax updated colorbar axis + + Examples : cax = pp.shaded_colorbar(inps, cax) + """ + # 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 + + ls = LightSource(azdeg=inps.shade_azdeg, altdeg=inps.shade_altdeg) + + num_x = 1000 # array size along the illumination profile + num_y = 255 # array size along the cmap + + arr = np.linspace(vmax, vmin, num_y) # data value array + arr = np.tile(arr, (num_x,1)).T + + x = np.linspace(-np.pi/4, np.pi/4, num_x) # x position along the illum. profile + ele = np.ones_like(arr) * np.cos(0.6*x) # altitude as a cosine along teh illum. profile + + dnorm = (arr - vmin) / (vmax - vmin) + cMap = plt.cm.ScalarMappable(norm=colors.Normalize(0,1), cmap=inps.colormap) + image = cMap.to_rgba(dnorm)[:, :, :3] + rgb = ls.shade_rgb(image, ele, fraction=fraction, blend_mode=blend_mode, vert_exag=vert_exag) + + + # show custom colorbar + cax.tick_params(axis='x', + which='both', # both major and minor ticks are affected + bottom=False, # ticks along the bottom edge are off + top=False, # ticks along the top edge are off + labelbottom=True) + cax.set_xticklabels([]) + cax.imshow(rgb, extent=[0, num_x, vmin, vmax], aspect='auto') + cax.set_ylim((vmin, vmax)) + cax.set_xlim((0, num_x)) + # can modify below outside this function for better visual + cax.set_ylabel(inps.cbar_label, rotation=-90, labelpad=15) + cax.yaxis.set_label_position("right") + cax.yaxis.tick_right() + return cax + + +def shaded_image(data, dem, ls, vmin=None, vmax=None, cmap='viridis', + base_color=0.9, shade_frac=0.5, blend_mode='hsv', vert_exag=0.5, + mask_dem_demNan=True, mask_dem_dataNan=False, mask_value=0): + """Add illumination to the rgb data array based on a DEM file. + Use shaded relief from a light source to adjust the colors of the data rgb array to have a cool impression. + Parameters : data 2D (m, n) np.array + dem 2D (m, n) np.int16 matrix, dem data + ls matplotlib.colors.LightSource class + vmin float; lower display limit of the data + vmax float; upper display limit of the data + cmap string 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 + vert_exag float; The amount to exaggerate the elevation values by when calculating illumination + mask_dem_demNan bool, True/False; whether to mask dem on nan dem pixels + mask_dem_dataNan bool, True/False; whether to mask dem on nan data pixels + mask_value float; set the masked pixels as alpha=mask_value (transparent) + + Returns : illum_rgb A 3D (m, n, 4) array of floats ranging between 0-1. + 1st to 3rd layers are the rgb values; 4th layer is the transparency + + Examples : illum_rgb = pp.shaded_image(data, dem, ls, vmin, vmax) + """ + # use numpy.ma to mask missing or invalid entries + data = np.ma.masked_invalid(data) + dem = np.ma.masked_invalid(dem) + + # data normalization + if not vmin: vmin = np.nanmin(data) + if not vmax: vmax = np.nanmax(data) + data_norm = (data-vmin) / (vmax-vmin) + + # cmap norm and ScalarMappable + mappable = plt.cm.ScalarMappable(norm=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 + if dem.shape != img_rgb.shape[:-1]: + print('DEM and Data dimension mismatch. Expecting errors in LightSource.shade_rgb()') + print(f'DEM shape {dem.shape}; Data shape {data.shape}') + + # add shaded relief to illuminate rgb image + 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_dem_dataNan: + # mask demShade on nan data pixels (default no) + illum_rgb[data.mask, -1] = mask_value + if mask_dem_demNan: + # mask demShade on nan dem pixels (default yes) + illum_rgb[dem.mask, -1] = mask_value + + return illum_rgb + + +def plot_image4view(ax, data, dem=None, extent=None, geo_box=None, inps=None): + """Plot image with DEM if provided + Parameters : ax matplotlib.pyplot.Axes or BasemapExt object + data 2D np.array + dem 2D np.int16 matrix, dem data + extent scalars (left, right, bottom, top) for pyplot.imshow() + geo_box tuple of 4 float in order of (E, N, W, S), geo bounding box + inps inps : Namespace with the following 7 items: + 'dem_file' : string; file path to the intput DEM + 'colormap' : string or matplotlib.colors.colormap class + 'vlim' : list, tuple of [float, float] + 'interpolation' : string + 'transparency' : float + 'animation' : string + 'disp_dem_blend' : bool, True/False + and other inps needed in plot_dem_background() + Returns : ax matplotlib.pyplot.Axes or BasemapExt object + im matplotlob.pyplot.AxesImage + + Examples : ax, im = pp.plot_image4view(ax, data, dem, inps=inps) + """ + imshow_data = True + + # print messages + global vprint + vprint = print if inps.print_msg else lambda *args, **kwargs: None + + # get plotting extent from geo_box if not given + if not extent is None: + if geo_box: + extent = (geo_box[0], geo_box[2], geo_box[3], geo_box[1]) # (W, E, S, N) + else: + print('WARNING: No extent and geo_box are given, expect to have issues in plot.') + + # shown dem style + if inps.dem_file: + if inps.disp_dem_blend: + dem_style = 'GMT-like image blended with DEM' + imshow_data = False + else: + dem_style = 'DEM background' + + # Plot DEM + if inps.dem_file: + vprint(f'plotting {dem_style} ...') + ax, im = plot_dem_background( + ax=ax, + geo_box=geo_box, + dem=dem, + data=data, + inps=inps, + print_msg=inps.print_msg, + ) + + if imshow_data: + # Plot Data + vprint('plotting image ...') + 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) + + return ax, im diff --git a/src/mintpy/view.py b/src/mintpy/view.py index 7b21e9439..539b7cbe2 100644 --- a/src/mintpy/view.py +++ b/src/mintpy/view.py @@ -485,6 +485,7 @@ def plot_slice(ax, data, metadata, inps): ).colormap # read DEM + dem = None # initialize for pp.plot_image4view() if inps.dem_file: dem, dem_metadata, dem_pix_box = pp.read_dem( inps.dem_file, @@ -512,20 +513,8 @@ def plot_slice(ax, data, metadata, inps): vprint(f'draw coast line with resolution: {inps.coastline}') ax.coastlines(resolution=inps.coastline, linewidth=inps.coastline_linewidth) - # Plot DEM - if inps.dem_file: - vprint('plotting DEM background ...') - pp.plot_dem_background( - ax=ax, - geo_box=inps.geo_box, - dem=dem, - inps=inps, - print_msg=inps.print_msg, - ) - - # Plot Data + # Whether reference to a gps 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 +526,8 @@ 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 the image + ax, im = pp.plot_image4view(ax, data, dem, extent=extent, geo_box=inps.geo_box, inps=inps) # Draw faultline using GMT lonlat file if inps.faultline_file: @@ -636,24 +624,12 @@ def format_coord(x, y): inps.fig_coord = 'yx' vprint('plotting in Y/X coordinate ...') - # Plot DEM - if inps.dem_file: - vprint('plotting DEM background ...') - pp.plot_dem_background( - ax=ax, - geo_box=None, - dem=dem, - inps=inps, - 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) + + ## Plot the image + ax, im = pp.plot_image4view(ax, data, dem, extent=extent, inps=inps) ax.tick_params(labelsize=inps.font_size) # Plot Seed Point @@ -737,7 +713,11 @@ def format_coord(x, y): if inps.disp_cbar: divider = make_axes_locatable(ax) cax = divider.append_axes(inps.cbar_loc, inps.cbar_size, pad=inps.cbar_size, axes_class=plt.Axes) - inps, cbar = pp.plot_colorbar(inps, im, cax) + + if not inps.disp_dem_blend: + inps, cbar = pp.plot_colorbar(inps, im, cax) + else: + cax = pp.shaded_colorbar(inps, cax) # 3.2 Title if inps.disp_title: @@ -1191,26 +1171,29 @@ def plot_subplot4figure(i, inps, ax, data, metadata): 1) Plot DEM, data and reference pixel 2) axes setting: tick, ticklabel, title, axis etc. """ + imshow_data = ~inps.disp_dem_blend # if disp_dem_blend, no need to imshow the data as a layer again + # Plot DEM if inps.dem_file: - pp.plot_dem_background( - ax=ax, - geo_box=None, - dem_shade=inps.dem_shade, - dem_contour=inps.dem_contour, - dem_contour_seq=inps.dem_contour_seq, - inps=inps, - 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) - - im = ax.imshow(data, cmap=inps.colormap, vmin=vlim[0], vmax=vlim[1], - interpolation=inps.interpolation, alpha=inps.transparency, - extent=extent, zorder=1) + ax, im = pp.plot_dem_background(ax, + dem_shade=inps.dem_shade, + dem_contour=inps.dem_contour, + dem_contour_seq=inps.dem_contour_seq, + dem=inps.dem, + data=data, + blend_img=inps.blend_img, + inps=inps, + print_msg=False) + + if imshow_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) + # Plot Data + im = ax.imshow(data, cmap=inps.colormap, vmin=vlim[0], vmax=vlim[1], + extent=extent, origin='upper', interpolation=inps.interpolation, + alpha=inps.transparency, zorder=1) # Plot Seed Point if inps.disp_ref_pixel: @@ -1406,7 +1389,11 @@ def adjust_subplots_layout(fig, inps): vprint('show colorbar') fig.subplots_adjust(right=0.93) cax = fig.add_axes([0.94, (1.0-cbar_length)/2, 0.005, cbar_length]) - inps, cbar = pp.plot_colorbar(inps, im, cax) + + if not inps.disp_dem_blend: + inps, cbar = pp.plot_colorbar(inps, im, cax) + else: + cax = pp.shaded_colorbar(inps, cax) else: adjust_subplots_layout(fig, inps) @@ -1494,7 +1481,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, inps.blend_img = pp.prepare_dem_background( dem=dem, inps=inps, print_msg=inps.print_msg, @@ -1507,6 +1494,7 @@ def prepare4multi_subplots(inps, metadata): msg += 'This feature is only supported for single subplot, and not for multi-subplots.' msg += '\n --> Ignore it and continue.' print(msg) + inps.dem = dem # store it for later call pp.plot_dem_background() return inps