This is the first question I have ever asked on Stack Overflow. Please forgive me if I make any mistakes when asking this question.

I am attempting to view several large images using the Python library, Datashader with Holoviews, using Bokeh as the backend. I would like the color map to be automatically re-scaled as I zoom into different portions of the image. I plan to view several versions of a similar image and also need to make it so that the panning/scrolling is simultaneous in all subplots. I understand that is the default when using Holoviews with Bokeh.

The main issue I have noticed is that while the resampling works with my images as I pan and zoom around, only the color map range in the last subplot is updated. All the others retain their original color map ranges and do not update. See the example images below.

Image shows the default view of two similar images when the browser loads.

The above image shows the default view of the plots once the browser loads. If I pan closer to the left side of the images, only the right image is re-stretched so that I can see the dimmer peak more easily.

Image shows only the color map range for the right image being updated.

Please see the example code below. Note that if axiswise is set to True, the auto color map re-scaling works in both images, instead of the just the right, but I loose the ability to pan and zoom simultaneously. Also, note that I am running the script by using bokeh serve [file_name].

I have been trying to get this to work for some time now, and I have noticed that there are no other posts on this. I would appreciate any help with this. Thank you!

import numpy as np
import holoviews as hv
import panel as pn
from holoviews.operation.datashader import rasterize
hv.extension('bokeh')


def gaus2d(x=0, y=0, mx=0, my=0, sx=1, sy=1):
    """
    2D Gaussian distribution
    """
    return 1. / (2. * np.pi * sx * sy) * np.exp(-((x - mx)**2. / (2. * sx**2.) + (y - my)**2. / (2. * sy**2.)))


# Get array of x and y values
N = 1_000
x = np.linspace(-1, 1, N)
y = np.linspace(-1, 1, N)
x, y = np.meshgrid(x, y)

# Calculate 2D Gaussian distribution for peaks 1 and 2
peak1 = 0.7*gaus2d(x=x, y=y, mx=-0.5, my=0, sx=0.1, sy=0.1)
peak2 = gaus2d(x=x, y=y, mx=0.5, my=0, sx=0.05, sy=0.05)
peaks = peak1 + peak2

# Plot parms
figsize = 300
axiswise = False  # If True, makes independent axes for each subplot

# Construct plot
img1 = rasterize(hv.Image(peaks)).opts(cmap="inferno", axiswise=axiswise, 
                                       height=figsize, width=figsize)
img2 = rasterize(hv.Image(peaks)).opts(cmap="inferno", axiswise=axiswise, 
                                       height=figsize, width=figsize)

fig_comb = img1 + img2

# Use this line if running "bokeh serve [file_name]" from terminal window
pn.panel(hv.Layout(fig_comb)).servable(title='test')