Hue and intensity visualization from two bands in python

Question

I would like to visualize this raster image (SAR interferogram) in a way that the phase band gives the color to the pixel and the coherence gives the intensity (value). This is the image: https://drive.google.com/file/d/1aeG4uE9UxZRihh88AqsCsPE_Dh0NwIWO/view?usp=share_link. Don't know where to start and I couldn't find anything (in QGis I wasn't able to find a way neither). This is a plotting of the two bands:

import os
import rasterio
from rasterio.plot import show
from rasterio.windows import from_bounds
from rasterio.enums import Resampling

from matplotlib import cm
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, BoundaryNorm

import numpy as np
import colorsys
import seaborn as sns



img = r"path ... _1SDV_20191103T051928_20191103T051955_029742_0363BA_83F0_split_oa_Stack_esd_ifg_Deb_DInSAR_TC_clip.tif"

#-----------------------------------------------------------------------
with rasterio.open(img) as ds:
    win_limits = [251303.801, 4941032.520, 273000.191, 4957748.474]
    #get phase
    phase = ds.read(3, window=from_bounds(251303.801, 4941032.520, 273000.191, 4957748.474, ds.transform))
    #get coherence
    coh = ds.read(4, window=from_bounds(251303.801, 4941032.520, 273000.191, 4957748.474, ds.transform))

#define HUE color map
my_colors = ['#808080', "#A14D04", "#A46905", "#AF8A44", "#B3A053",
             '#AE9F5C','#A8A965', '#ADB55F', '#AAB86D', '#A5BC76',
             '#94B660','#95B645', '#7DAE38', '#7BA654', '#759E72',
            '#698C86', '#6D9298', '#568F96','#4B80A0', '#327CBB',
            '#316DC5', '#2158BC', '#000000', '#808080']
my_cmap = ListedColormap(my_colors)

# bounds of colors
bounds = [-2.64117, -2.640942, -2.56253, -2.475379,-2.373972,
          -2.270036, -2.15046,-2.023531,-1.895680,-1.758170,
          -1.63146,-1.492810,-1.327017,-1.089251,-0.83538,
          -0.39871,0.985112924,1.539288,1.868574,2.044485,
          2.15325057,2.268684654,2.6891]

my_norm = BoundaryNorm(bounds, ncolors=len(my_colors))


##plot
fig, ((axp, axc), (cbar_axp, cbar_axc)) = plt.subplots(2,2, gridspec_kw={'width_ratios':[1,1], 'wspace':0.1, 'height_ratios':[1.1,0.1],  'hspace': 0.1})
i1 = sns.heatmap(phase, 
                ax=axp,
                cmap=my_cmap,
                norm=my_norm,
                cbar_ax=cbar_axp,
                cbar_kws={"orientation": "horizontal"}).set(yticklabels=[],xticklabels=[], title="Phase")
colorbar = axp.collections[0].colorbar
i2 = sns.heatmap(coh,
                 cmap="Greys",
                 ax=axc,
                 cbar_ax = cbar_axc,
                 cbar_kws={"orientation": "horizontal"}).set(yticklabels=[],xticklabels=[], title="Coherence")
plt.show()

In case it is not possible I could also be happy with overlaying the coherence band (with some transparency) on top of the phase band.

Any ideas?