Repository: holoviz/holoviews
Pull request #5840 (2023-08-02 14:56:12)

Implement support for subcoordinate systems in the y-axis
Implement support for subcoordinate systems in the y-axis

<img width="1253" alt="Screen Shot 2023-08-02 at 16 50 54" src="https://github.com/holoviz/holoviews/assets/1550771/20b312ef-b9a9-4523-a072-0fec10e102e0">



droumis: This is great!

From a first pass, I've identified a few related issues:

1. **Different scaling across subplots:** From my understanding, the subcoordinate_y parameter allows you to define a specific y-axis range for a subplot. The data that is plotted within each subplot will be scaled to fit within the y-axis range defined by subcoordinate_y, correct? This means that if the data in each subplot has a different min/max data amplitude (always the case), then each subplot will have its own y-axis scale compared to the other subplots, which is not desired. We want there to be a single scale bar applicable for a set of subplots.

2. **Floating y origin:** The first value of each trace is floating within the subplot range, but the desired behavior would be for the first value to align with its ytick.

3. Because each curve element is now on it's own subplot, and due to the issue of the [RangeTool not being able to use an overlay for the target](https://github.com/holoviz/holoviews/issues/4472#issuecomment-1592685063), it seems like the workaround of referencing a single element as the target no longer works to have the RangeTool control the y-axis linking. See video below; You can see the plot shows all the channels rather than a few, and that RangeTool works for the x axis but not the y.

https://github.com/holoviz/holoviews/assets/6613202/a42dbcae-df1a-4ba8-8a4b-91dc22517c96

4. **y_range.bounds:** The y_range.bounds applied to the overlay via backend_opts are no longer not respected; I can pan and zoom beyond them.





droumis: **UPDATE:** I am not so sure what's happening anymore.. the data for a particular curve seems to be allowed to go beyond the assigned subcoordinate_y range so I am just confused now.

![image](https://github.com/holoviz/holoviews/assets/6613202/77b28c70-4aef-4267-b062-5a79f2397123)


Here is the code to reproduce:

```
import numpy as np
import pandas as pd
import holoviews as hv; hv.extension('bokeh')
from holoviews import opts
from holoviews import Dataset
from holoviews.plotting.links import RangeToolLink
from bokeh.models import HoverTool
import panel as pn; pn.extension()
from scipy.stats import zscore

n_channels = 10
n_seconds = 5
fs = 512
max_ch_disp = 5  # max channels to initially display
max_t_disp = 3 # max time in seconds to initially display

total_samples = n_seconds * fs
time = np.linspace(0, n_seconds, total_samples)
data = np.random.randn(n_channels, total_samples).cumsum(axis=1)
channels = ['EEG {}'.format(i) for i in range(n_channels)]

channel_curves = []
ch_subcoord = 1./n_channels

hover = HoverTool(tooltips=[
    ("Channel", "@channel"),
    ("Time", "$x s"),
    ("Amplitude", "$y µV")])

for i, channel_data in enumerate(data):
    ds = Dataset((time, channel_data, channels[i]), ["Time", "Amplitude", "channel"])
    channel_curves.append(hv.Curve(ds, "Time", ["Amplitude", "channel"], label=channels[i]).opts(
            subcoordinate_y=(i*ch_subcoord, (i+1)*ch_subcoord), color="black", line_width=1, 
            tools=[hover, 'xwheel_zoom'], shared_axes=False))

eeg_viewer = hv.Overlay(channel_curves, kdims="Channel").opts(
    padding=0, xlabel="Time (s)", ylabel="Channel",
    show_legend=False, aspect=1.5, responsive=True, # yticks=yticks,
    shared_axes=False, backend_opts={
        "x_range.bounds": (time.min(), time.max()),
        "y_range.bounds": (data.min(), data.max())})

yticks = [( (i*ch_subcoord + (i+1)*ch_subcoord) / 2, ich) for i, ich in enumerate(channels)]
y_positions, _ = zip(*yticks)

z_data = zscore(data, axis=1)

minimap = hv.Image((time, y_positions, z_data), ["Time (s)", "Channel"], "Amplitude (uV)")
minimap = minimap.opts(
    cmap="RdBu_r", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],
    height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))


if len(channels) < max_ch_disp:
    max_ch_disp = len(channels)
max_y_disp = (max_ch_disp+2)*ch_subcoord

time_s = len(time)/fs
if time_s < max_t_disp:
    max_t_disp = time_s
    
RangeToolLink(minimap, list(eeg_viewer.values())[0], axes=["x", "y"],
              boundsx=(None, max_t_disp),
              boundsy=(None, max_y_disp))

eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500), minimap).servable(target='main')
eeg_app
```
droumis: pinging @mattpap so we can stay aligned on this effort
droumis: 
5. **VSpan:** It seems like I cannot use hv.VSpan on an element or overlay of elements that use `subcoordinate_y`

```python
import numpy as np
import holoviews as hv
hv.extension('bokeh')

x = np.linspace(0, 10, 100)
y = np.sin(x)
curve = hv.Curve((x, y)).opts(subcoordinate_y=(0,1))
vspan = hv.VSpan(2, 4)
curve *= vspan
curve

```
Traceback:
```
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/IPython/core/formatters.py:974, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    971     method = get_real_method(obj, self.print_method)
    973     if method is not None:
--> 974         return method(include=include, exclude=exclude)
    975     return None
    976 else:

File ~/src/holoviews/holoviews/core/dimension.py:1287, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1280 def _repr_mimebundle_(self, include=None, exclude=None):
   1281     """
   1282     Resolves the class hierarchy for the class rendering the
   1283     object using any display hooks registered on Store.display
   1284     hooks.  The output of all registered display_hooks is then
   1285     combined and returned.
   1286     """
-> 1287     return Store.render(self)

File ~/src/holoviews/holoviews/core/options.py:1423, in Store.render(cls, obj)
   1421 data, metadata = {}, {}
   1422 for hook in hooks:
-> 1423     ret = hook(obj)
   1424     if ret is None:
   1425         continue

File ~/src/holoviews/holoviews/ipython/display_hooks.py:280, in pprint_display(obj)
    278 if not ip.display_formatter.formatters['text/plain'].pprint:
    279     return None
--> 280 return display(obj, raw_output=True)

File ~/src/holoviews/holoviews/ipython/display_hooks.py:248, in display(obj, raw_output, **kwargs)
    246 elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    247     with option_state(obj):
--> 248         output = element_display(obj)
    249 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    250     with option_state(obj):

File ~/src/holoviews/holoviews/ipython/display_hooks.py:142, in display_hook.<locals>.wrapped(element)
    140 try:
    141     max_frames = OutputSettings.options['max_frames']
--> 142     mimebundle = fn(element, max_frames=max_frames)
    143     if mimebundle is None:
    144         return {}, {}

File ~/src/holoviews/holoviews/ipython/display_hooks.py:188, in element_display(element, max_frames)
    185 if type(element) not in Store.registry[backend]:
    186     return None
--> 188 return render(element)

File ~/src/holoviews/holoviews/ipython/display_hooks.py:69, in render(obj, **kwargs)
     66 if renderer.fig == 'pdf':
     67     renderer = renderer.instance(fig='png')
---> 69 return renderer.components(obj, **kwargs)

File ~/src/holoviews/holoviews/plotting/renderer.py:397, in Renderer.components(self, obj, fmt, comm, **kwargs)
    394 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    396 if embed or config.comms == 'default':
--> 397     return self._render_panel(plot, embed, comm)
    398 return self._render_ipywidget(plot)

File ~/src/holoviews/holoviews/plotting/renderer.py:404, in Renderer._render_panel(self, plot, embed, comm)
    402 doc = Document()
    403 with config.set(embed=embed):
--> 404     model = plot.layout._render_model(doc, comm)
    405 if embed:
    406     return render_model(model, comm)

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/viewable.py:736, in Viewable._render_model(self, doc, comm)
    734 if comm is None:
    735     comm = state._comm_manager.get_server_comm()
--> 736 model = self.get_root(doc, comm)
    738 if self._design and self._design.theme.bokeh_theme:
    739     doc.theme = self._design.theme.bokeh_theme

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/layout/base.py:305, in Panel.get_root(self, doc, comm, preprocess)
    301 def get_root(
    302     self, doc: Optional[Document] = None, comm: Optional[Comm] = None,
    303     preprocess: bool = True
    304 ) -> Model:
--> 305     root = super().get_root(doc, comm, preprocess)
    306     # ALERT: Find a better way to handle this
    307     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/viewable.py:658, in Renderable.get_root(self, doc, comm, preprocess)
    656 wrapper = self._design._wrapper(self)
    657 if wrapper is self:
--> 658     root = self._get_model(doc, comm=comm)
    659     if preprocess:
    660         self._preprocess(root)

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/layout/base.py:173, in Panel._get_model(self, doc, root, parent, comm)
    171 root = root or model
    172 self._models[root.ref['id']] = (model, parent)
--> 173 objects, _ = self._get_objects(model, [], doc, root, comm)
    174 props = self._get_properties(doc)
    175 props[self._property_mapping['objects']] = objects

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/layout/base.py:155, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    153 else:
    154     try:
--> 155         child = pane._get_model(doc, root, model, comm)
    156     except RerenderError as e:
    157         if e.layout is not None and e.layout is not self:

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/pane/holoviews.py:411, in HoloViews._get_model(self, doc, root, parent, comm)
    409     plot = self.object
    410 else:
--> 411     plot = self._render(doc, comm, root)
    413 plot.pane = self
    414 backend = plot.renderer.backend

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/panel/pane/holoviews.py:506, in HoloViews._render(self, doc, comm, root)
    503     if comm:
    504         kwargs['comm'] = comm
--> 506 return renderer.get_plot(self.object, **kwargs)

File ~/src/holoviews/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/src/holoviews/holoviews/plotting/renderer.py:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

File ~/src/holoviews/holoviews/plotting/plot.py:943, in DimensionedPlot.update(self, key)
    941 def update(self, key):
    942     if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
--> 943         return self.initialize_plot()
    944     item = self.__getitem__(key)
    945     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/src/holoviews/holoviews/plotting/bokeh/element.py:2782, in OverlayPlot.initialize_plot(self, ranges, plot, plots)
   2779 self.handles['plot'] = plot
   2781 if plot and not self.overlaid:
-> 2782     self._update_plot(key, plot, element)
   2783     self._update_ranges(element, ranges)
   2785 panels = []

File ~/src/holoviews/holoviews/plotting/bokeh/element.py:877, in ElementPlot._update_plot(self, key, plot, element)
    875 plot.update(**self._plot_properties(key, element))
    876 if not self.multi_y:
--> 877     self._update_labels(key, plot, element)
    878 self._update_title(key, plot, element)
    879 self._update_grid(plot)

File ~/src/holoviews/holoviews/plotting/bokeh/element.py:885, in ElementPlot._update_labels(self, key, plot, element)
    883 el = el[0] if el else element
    884 dimensions = self._get_axis_dims(el)
--> 885 props = {axis: self._axis_properties(axis, key, plot, dim)
    886          for axis, dim in zip(['x', 'y'], dimensions)}
    887 xlabel, ylabel, zlabel = self._get_axis_labels(dimensions)
    888 if self.invert_axes:

File ~/src/holoviews/holoviews/plotting/bokeh/element.py:885, in <dictcomp>(.0)
    883 el = el[0] if el else element
    884 dimensions = self._get_axis_dims(el)
--> 885 props = {axis: self._axis_properties(axis, key, plot, dim)
    886          for axis, dim in zip(['x', 'y'], dimensions)}
    887 xlabel, ylabel, zlabel = self._get_axis_labels(dimensions)
    888 if self.invert_axes:

File ~/src/holoviews/holoviews/plotting/bokeh/element.py:815, in ElementPlot._axis_properties(self, axis, key, plot, dimension, ax_mapping)
    813 for i, (frame, coords) in enumerate(reversed(list(zip(self.current_frame, subcoords)))):
    814     labels.append(frame.label or f'Trace {i}')
--> 815     ticks.append(np.mean(coords))
    816 axis_props['ticker'] = FixedTicker(ticks=ticks)
    817 if labels is not None:

File <__array_function__ internals>:200, in mean(*args, **kwargs)

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3464, in mean(a, axis, dtype, out, keepdims, where)
   3461     else:
   3462         return mean(axis=axis, dtype=dtype, out=out, **kwargs)
-> 3464 return _methods._mean(a, axis=axis, dtype=dtype,
   3465                       out=out, **kwargs)

File ~/opt/miniconda3/envs/neuro-eeg-viewer/lib/python3.9/site-packages/numpy/core/_methods.py:194, in _mean(a, axis, dtype, out, keepdims, where)
    192         ret = ret.dtype.type(ret / rcount)
    193 else:
--> 194     ret = ret / rcount
    196 return ret

TypeError: unsupported operand type(s) for /: 'NoneType' and 'int'
:Overlay
   .Curve.I :Curve   [x]   (y)
   .VSpan.I :VSpan   [x,y]
```
droumis: @philippjfr's awesome latest pushes resolve a number of issues:

1. **Different scaling across subplots:** The subcoordinate ranges are now computed at the OverlayPlot level
2. The rangetoollink can now use an Overlay as a target. (fixes  #4472)
3. The workaround prior to # 2 was to link to a nested element.. but that was the cause of the first element being scaled in the subplot incorrectly.. so # 2 also fixed [my observation](https://github.com/holoviz/holoviews/pull/5840#issuecomment-1663507724) above about the first curve going beyond the subplot y_range.
4. **VSpan:** now works on elements with subcoordinate_y defined

Remaining issues:
1. **Floating y origin:** The first value of each trace is floating within the subplot range, but the desired behavior would be for the first value to align with its ytick.
2. **y_range.bounds:** The y_range.bounds applied to the overlay via backend_opts are no longer not respected; I can pan and zoom beyond them.


I think resolving **Floating y origin:** is the main thing left to address for this PR


<details>
<summary> Code: </summary>
```
import numpy as np
import pandas as pd
import holoviews as hv; hv.extension('bokeh')
from holoviews import opts
from holoviews import Dataset
from holoviews.plotting.links import RangeToolLink
from bokeh.models import HoverTool
import panel as pn; pn.extension()
from scipy.stats import zscore

n_channels = 10
n_seconds = 5
fs = 512
max_ch_disp = 5  # max channels to initially display
max_t_disp = 3 # max time in seconds to initially display

total_samples = n_seconds * fs
time = np.linspace(0, n_seconds, total_samples)
data = np.random.randn(n_channels, total_samples).cumsum(axis=1)
channels = ['EEG {}'.format(i) for i in range(n_channels)]

channel_curves = []
ch_subcoord = 1./n_channels

hover = HoverTool(tooltips=[
    ("Channel", "@channel"),
    ("Time", "$x s"),
    ("Amplitude", "$y µV")])

for i, channel_data in enumerate(data):
    ds = Dataset((time, channel_data, channels[i]), ["Time", "Amplitude", "channel"])
    channel_curves.append(hv.Curve(ds, "Time", ["Amplitude", "channel"], label=channels[i]).opts(
            subcoordinate_y=(i*ch_subcoord, (i+1)*ch_subcoord), color="black", line_width=1, 
            tools=[hover, 'xwheel_zoom'], shared_axes=False))

annotation = hv.VSpan(0.3, 0.5)
eeg_viewer = (hv.Overlay(channel_curves, kdims="Channel") * annotation).opts(
    padding=0, xlabel="Time (s)", ylabel="Channel",
    show_legend=False, aspect=1.5, min_height=500, responsive=True,
    shared_axes=False, backend_opts={
        "x_range.bounds": (time.min(), time.max()),
        "y_range.bounds": (data.min(), data.max())}) 

yticks = [( (i*ch_subcoord + (i+1)*ch_subcoord) / 2, ich) for i, ich in enumerate(channels)]
y_positions, _ = zip(*yticks)

z_data = zscore(data, axis=1)

minimap = hv.Image((time, y_positions, z_data), ["Time (s)", "Channel"], "Amplitude (uV)")
minimap = minimap.opts(
    cmap="RdBu_r", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],
    height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))

if len(channels) < max_ch_disp:
    max_ch_disp = len(channels)
max_y_disp = (max_ch_disp+2)*ch_subcoord

time_s = len(time)/fs
if time_s < max_t_disp:
    max_t_disp = time_s
    
RangeToolLink(minimap, eeg_viewer, axes=["x", "y"],
              boundsx=(None, max_t_disp),
              boundsy=(None, max_y_disp))

eeg_app = pn.Column((eeg_viewer + minimap * annotation).cols(1), min_height=650).servable(target='main')
eeg_app
```
</details>

![image](https://github.com/holoviz/holoviews/assets/6613202/6bc18681-b4dc-4efb-a7d4-c34ae3cd08d8)
droumis: > Very cool! I'm confused about the floating y origin; that already seems to be implemented in the plot above? In any case I don't think that is a must-have, but certainly a nice to have.

I think it just appears to be aligned with the y axis when the individual traces are scaled down a lot. Here is the same code with some real data.. it's scaled down so much that it isn't that useful of a view:

(Relevant line per `i` channel: `subcoordinate_y=(i*(1./n_channels), (i+1)*(1./n_channels))`)
![image](https://github.com/holoviz/holoviews/assets/6613202/c4a4bb4a-4566-4c1a-9bf6-568f3bb7a92f)


However, when I 'scale up' the traces by adjusting the range that each subplot is allowed to be within, the data no longer aligns with the y-tick, and it's especially bad when panning across the channels.

(Relevant line per `i` channel: `subcoordinate_y=(i*(1./n_channels), (i+10)*(1./n_channels))`)

https://github.com/holoviz/holoviews/assets/6613202/97be1c95-ab81-487b-94aa-33a9ca3171dc

I think instead of defining the subplot by the upper and lower y_range that it's allowed to be in, we will need to define it by the offset for the y-origin, and then have a separate parameter for the scaling to use for each group of curves. @philippjfr, do you think this would be possible without major changes to Bokeh?

mattpap: > However, when I 'scale up' the traces by adjusting the range that each subplot is allowed to be within, the data no longer aligns with the y-tick, and it's especially bad when panning across the channels.

If I understand the issue correctly, then the solution is to compute y source range for each sub-coordinate based on the data indices in the viewport (i.e. indices contained within the cartesian frame) and not all the available indices as it is done now. That would mean using data ranges for y source ranges and changes to bokeh, to allow data ranges optionally to consider the viewport (so called masked indices). I doubt this can be hacked around in Holoviews, because masked indices are not exposed to Python APIs in bokeh.
droumis: @philippjfr, @mattpap. To simplify things, I've decided to now drop the requested alignment of the y-tick with the left-most data point. The type of data we are working with for CZI should have very low frequencies ('slow drift') removed prior to plotting, so I think it's good enough as long as the y-tick is reasonably close to the trace. However, reasonable alignment of the ytick when panning with overlapping ranges is still an issue with this PR, as was shown in this video:

https://github.com/holoviz/holoviews/assets/6613202/97be1c95-ab81-487b-94aa-33a9ca3171dc

Any idea how to resolve this, @philippjfr? As @mattpap mentioned in the meeting yesterday, he does not see such behavior in a pure Bokeh implementation. 

If we can get this resolved (and also have [zoom tools to scale subcoordinate ranges](https://github.com/bokeh/bokeh/issues/13339)), then I think that's all we need for now.
mattpap: Reproduction in pure bokeh for both categorical and linear axis approaches (change `categorical` global variable):
```python
import numpy as np

from bokeh.core.properties import field
from bokeh.io import show
from bokeh.models import (ColumnDataSource, DataRange1d, FactorRange, HoverTool,
                          Range1d, WheelZoomTool, ZoomInTool, ZoomOutTool)
from bokeh.palettes import Category10
from bokeh.plotting import figure

np.random.seed(0)

categorical = False

n_channels = 10
n_seconds = 15
fs = 512
max_ch_disp = 5  # max channels to initially display
max_t_disp = 3 # max time in seconds to initially display

total_samples = n_seconds * fs
time = np.linspace(0, n_seconds, total_samples)
data = np.random.randn(n_channels, total_samples).cumsum(axis=1)
channels = [f'EEG {i}' for i in range(n_channels)]

hover = HoverTool(tooltips=[
    ("Channel", "$name"),
    ("Time", "$x s"),
    ("Amplitude", "$y µV"),
])

x_range = Range1d(start=time.min(), end=time.max())
#x_range = DataRange1d()
#y_range = DataRange1d() # bug
if categorical:
    y_range = FactorRange(factors=channels)
else:
    y_range = Range1d(start=-0.5, end=len(channels) - 1 + 0.5)

p = figure(x_range=x_range, y_range=y_range, lod_threshold=None)

source = ColumnDataSource(data=dict(time=time))
renderers = []

for i, channel in enumerate(channels):
    if categorical:
        y_target=Range1d(start=i, end=i + 1)
    else:
        y_target=Range1d(start=i - 0.5, end=i + 0.5)

    xy = p.subplot(
        x_source=p.x_range,
        y_source=Range1d(start=data[i].min(), end=data[i].max()),
        #y_source=DataRange1d(),
        x_target=p.x_range,
        y_target=y_target,
    )

    source.data[channel] = data[i]
    line = xy.line(field("time"), field(channel), color=Category10[10][i], source=source, name=channel)
    renderers.append(line)

if not categorical:
    from bokeh.models import FixedTicker
    ticks = list(range(len(channels)))
    p.yaxis.ticker = FixedTicker(ticks=ticks)
    p.yaxis.major_label_overrides = {i: f"EEG {i}" for i in ticks}

wheel_zoom = WheelZoomTool()#renderers=renderers)
zoom_in = ZoomInTool(renderers=renderers)
zoom_out = ZoomOutTool(renderers=renderers)

p.add_tools(wheel_zoom, zoom_in, zoom_out, hover)
p.toolbar.active_scroll = wheel_zoom

show(p)

```
philippjfr: Interesting, so in your pure-bokeh version I can't reproduce the panning issue.
maximlt: Just for reference some pure Bokeh code to reproduce what we want i.e. a plot with traces on subcoordinates with its minimap.

![Kapture 2023-09-06 at 17 07 12](https://github.com/holoviz/holoviews/assets/35924738/e1b796a6-5e3c-4159-bd62-c14b1af84621)


<details>
<summary>Code</summary>

```python
import numpy as np

from bokeh.layouts import column
from bokeh.core.properties import field
from bokeh.io import show
from bokeh.models import (ColumnDataSource, DataRange1d, FactorRange, HoverTool,
                          Range1d, WheelZoomTool, ZoomInTool, ZoomOutTool, RangeTool)
from bokeh.palettes import Category10
from bokeh.plotting import figure

from scipy.stats import zscore

np.random.seed(0)

categorical = False

n_channels = 10
n_seconds = 15
fs = 512
max_ch_disp = 5  # max channels to initially display
max_t_disp = 3 # max time in seconds to initially display

total_samples = n_seconds * fs
time = np.linspace(0, n_seconds, total_samples)
data = np.random.randn(n_channels, total_samples).cumsum(axis=1)
channels = [f'EEG {i}' for i in range(n_channels)]

hover = HoverTool(tooltips=[
    ("Channel", "$name"),
    ("Time", "$x s"),
    ("Amplitude", "$y µV"),
])

x_range = Range1d(start=time.min(), end=time.max())
#x_range = DataRange1d()
#y_range = DataRange1d() # bug
if categorical:
    y_range = FactorRange(factors=channels)
else:
    y_range = Range1d(start=-0.5, end=len(channels) - 1 + 0.5)

p = figure(x_range=x_range, y_range=y_range, lod_threshold=None)

source = ColumnDataSource(data=dict(time=time))
renderers = []

for i, channel in enumerate(channels):
    if categorical:
        y_target=Range1d(start=i, end=i + 1)
    else:
        y_target=Range1d(start=i - 0.5, end=i + 0.5)

    xy = p.subplot(
        x_source=p.x_range,
        y_source=Range1d(start=data[i].min(), end=data[i].max()),
        #y_source=DataRange1d(),
        x_target=p.x_range,
        y_target=y_target,
    )

    source.data[channel] = data[i]
    line = xy.line(field("time"), field(channel), color=Category10[10][i], source=source, name=channel)
    renderers.append(line)

if not categorical:
    from bokeh.models import FixedTicker
    ticks = list(range(len(channels)))
    p.yaxis.ticker = FixedTicker(ticks=ticks)
    p.yaxis.major_label_overrides = {i: f"EEG {i}" for i in ticks}

wheel_zoom = WheelZoomTool()#renderers=renderers)
zoom_in = ZoomInTool(renderers=renderers)
zoom_out = ZoomOutTool(renderers=renderers)

p.add_tools(wheel_zoom, zoom_in, zoom_out, hover)
p.toolbar.active_scroll = wheel_zoom

z_data = zscore(data, axis=1)

range_tool = RangeTool(x_range=p.x_range, y_range=p.y_range)
range_tool.x_range.update(start=4, end=6)
range_tool.y_range.update(start=2, end=8)

select = figure(height=130, #y_range=p.y_range,
                tools="", toolbar_location=None)
select.image(image=[z_data], x=0, y=0, dw=15, dh=11, palette="Sunset11")
select.add_tools(range_tool)

show(column(p, select))
```
</details>

maximlt: The original implementation had some code to handle when no `label` is set for an element, supposedly displaying each element with a tick value of `Trace <i>`. However, this code didn't seem to work, I restored support for it in https://github.com/holoviz/holoviews/pull/5840/commits/3b19d424b31f783322f0a9889ab5d7ab4400be54. The approach I've used in the `else` block seems a little fragile to me, I couldn't find any other way though. If that's deemed not correct or too fragile, I think forcing users to provide a `label` value per element isn't too bad, after all `Trace <i>` isn't that helpful as a tick value!

I've added various tests, some inspired by the multi-axis tests and some not.  The multi-axis tests were exerting features like log y-axis and inverted y-axis, which I guess are out of scope for this PR (if it anyway makes any sense to support it, and if that's even possible). I'm interested in feedback on the test, in particular to test features that are expected to be delivered for CZI if they aren't currently tested.

There is more "padding" over the last trace than below the first one, I'm not sure that was intentional:
```python
import holoviews as hv
from holoviews import Curve, Overlay
import numpy as np

hv.extension('bokeh')

N = 5
data = [np.random.randn(10000).cumsum() for _ in range(N)]
Overlay([Curve(data[i], label=f'Data {i}').opts(subcoordinate_y=True) for i in range(N)]).opts(show_legend=False)
```
<img width="316" alt="image" src="https://github.com/holoviz/holoviews/assets/35924738/1f094612-44b2-456f-b565-f46cd7c325be">

I've re-added some basic support for setting `subcoordinate_y` with a range, ignoring the scale. The range values are expected to be between 0 and 1. This is different compared to `subcoordinate_y=True`, where the y-range goes from _-scale/2_ to _n_curves + (scale/2)_, with each element center position on the y-axis being located at _i_ (0, 1, 2, etc.).

When an element is set with `subcoordinate_y=True`, it's no longer possible to display it on its own, unless you do `element.opts(subcoordinate=False)` of course. Not sure how big of an issue this is, I thought it was worth mentioning.

I have updated the EEG example to work with this PR, it also needs https://github.com/holoviz/holoviews/pull/5881. I think it's a nice example, which I suggest we add to the gallery? One trick to get it to work as expected was to make sure the y-range of both plots is the same.

<details open>
<summary>Code</summary>

```python
import numpy as np
import pandas as pd
import holoviews as hv
from holoviews import opts
from holoviews import Dataset
from holoviews.plotting.links import RangeToolLink
from bokeh.models import HoverTool
from scipy.stats import zscore

hv.extension('bokeh')

n_channels = 10
n_seconds = 5
fs = 512
max_ch_disp = 5  # max channels to initially display
max_t_disp = 3 # max time in seconds to initially display

total_samples = n_seconds * fs
time = np.linspace(0, n_seconds, total_samples)
data = np.random.randn(n_channels, total_samples).cumsum(axis=1)
channels = ['EEG {}'.format(i) for i in range(n_channels)]

channel_curves = []

hover = HoverTool(tooltips=[
    ("Channel", "@channel"),
    ("Time", "$x s"),
    ("Amplitude", "$y µV")])

for i, channel_data in enumerate(data):
    ds = Dataset((time, channel_data, channels[i]), ["Time", "Amplitude", "channel"])
    channel_curves.append(hv.Curve(ds, "Time", ["Amplitude", "channel"], label=channels[i]).opts(
            subcoordinate_y=True, color="black", line_width=1, 
            tools=[hover, 'xwheel_zoom'], shared_axes=False))

eeg_viewer = hv.Overlay(channel_curves, kdims="Channel").opts(
    padding=0, xlabel="Time (s)", ylabel="Channel",
    show_legend=False, aspect=4, responsive=True,# yticks=yticks,
    shared_axes=False,
    backend_opts={
        "x_range.bounds": (time.min(), time.max()),
        "y_range.bounds": (data.min(), data.max())
    }
)

yticks = [(i , ich) for i, ich in enumerate(channels)]
y_positions, _ = zip(*yticks)

z_data = zscore(data, axis=1)

minimap = hv.Image((time, y_positions, z_data), ["Time (s)", "Channel"], "Amplitude (uV)")
minimap = minimap.opts(
    cmap="RdBu_r", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],
    height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))

# The subcoordinate scale should be accounted for here if set.
if len(channels) < max_ch_disp:
    max_ch_disp = len(channels)
max_y_disp = max_ch_disp - 0.1

time_s = len(time)/fs
if time_s < max_t_disp:
    max_t_disp = time_s
    
RangeToolLink(minimap, eeg_viewer, axes=["x", "y"],
              boundsx=(None, max_t_disp),
              boundsy=(None, max_y_disp)
)


app = (eeg_viewer + minimap).opts(merge_tools=False).cols(1)
app
```

</details>

<img width="1178" alt="image" src="https://github.com/holoviz/holoviews/assets/35924738/6e9baa5d-217f-45cd-bd35-94ff15465aa7">


Here are a few places where this feature could be documented:

- Where the `Overlay` is introduced in the _Composing Elements_ guide: https://holoviews.org/user_guide/Composing_Elements.html#overlay
- Next to `multi_y` in the _Customizing Plots_ guide: https://holoviews.org/user_guide/Customizing_Plots.html#twin-axes
- In the gallery with the EEG example as suggested above

philippjfr: >The original implementation had some code to handle when no label is set for an element, supposedly displaying each element with a tick value of Trace <i>. However, this code didn't seem to work, I restored support for it in 3b19d42. The approach I've used in the else block seems a little fragile to me, I couldn't find any other way though. If that's deemed not correct or too fragile, I think forcing users to provide a label value per element isn't too bad, after all Trace <i> isn't that helpful as a tick value!

Could we maybe just use the vdim name instead? Even if it's duplicated that's at least something.

>I'm interested in feedback on the test, in particular to test features that are expected to be delivered for CZI if they aren't currently tested.

I agree that log axes don't make sense and I'm also okay with skipping inverted axes.

>There is more "padding" over the last trace than below the first one, I'm not sure that was intentional:

That looks/sounds like an off-by-one error.

>I've re-added some basic support for setting subcoordinate_y with a range, ignoring the scale. The range values are expected to be between 0 and 1. This is different compared to subcoordinate_y=True, where the y-range goes from -scale/2 to n_curves + (scale/2), with each element center position on the y-axis being located at i (0, 1, 2, etc.). Maybe the two approaches should better aligned?

Definitely, ideally I was hoping that there would be an explicit way to give a scale to the target range so you could actually make it meaningful, e.g. by labeling probe depth or similar along the axis. Did not figure out the best approach for doing this though.

>When an element is set with subcoordinate_y=True, it's no longer possible to display it on its own, unless you do element.opts(subcoordinate=False) of course. Not sure how big of an issue this is, I thought it was worth mentioning.

I think in the non-overlaid case we need to simply need to ignore subcoordinate_y.

>I have updated the EEG example to work with this PR, it also needs #5881. I think it's a nice example, which I suggest we add to the gallery? One trick to get it to work as expected was to make sure the y-range of both plots is the same.

Yes please add it.

>Next to multi_y in the Customizing Plots guide: https://holoviews.org/user_guide/Customizing_Plots.html#twin-axes

This is where I'd put it.
droumis: Hey @maximlt , great work! I hope it's ok that I went ahead and committed my revisions to the documentation. 

For the user guide, I added that the label is required, fixed typos, rephrased, and added some clarification that the scale of 2 will overlap 50% of the adjacent if the adjacent curves are not scaled. I also added a bit of scaling to the code example. 

For the eeg demo, I changed the data from a random walk to incrementing sine waves. Although it's still very fake, it's a bit more useful of a view for EEG researchers, it takes up more of the space, and it makes for a pretty minimap image. I also rephrased the text and simplified the minimap image clims. 

![image](https://github.com/holoviz/holoviews/assets/6613202/e133738a-0f81-4821-8cdc-d1f103efac5f)

jbednar: >the label is required

Aren't there scenarios that are equivalent to an NdOverlay where the label can be computed from the value of a dimension? I'd think this would have the same interface as an NdOverlay, but maybe I'm missing something.
droumis: > > the label is required
> 
> Aren't there scenarios that are equivalent to an NdOverlay where the label can be computed from the value of a dimension? I'd think this would have the same interface as an NdOverlay, but maybe I'm missing something.

Probably, but we concluded that it was easy enough for users to set a label for now and that it covers more use cases to just be explicit. Some sort of auto-labeling can be addressed in a future PR.
maximlt: @droumis yes sure no problem! I looked at your changes and they're all fine by me.
philippjfr: For some of those edge cases, I don't think they are edge cases. Specifically the "Same subcordinate_y case" is perhaps silly but we are not in the business of preventing users doing silly things. The "Only some curves have enabled subcoordinate_y" is also perfectly valid, in fact this is an important feature because you want to be able to render into the main axis. There is some work to do, to allow the target axis to represent a meaningful coordinate space but even in the current form you might want to render annotations into the target range instead of a subcoordinate range.
Hoxbro: Edge case was maybe not the right word. I wrote it when I only had the first examples. 

> Specifically the "Same subcordinate_y case" is perhaps silly but we are not in the business of preventing users doing silly things

Fair point. But in this case, the labels are different in the code but the same in the figure.   

> The "Only some curves have enabled subcoordinate_y" is also perfectly valid, in fact this is an important feature because you want to be able to render into the main axis. There is some work to do, to allow the target axis to represent a meaningful coordinate space but even in the current form you might want to render annotations into the target range instead of a subcoordinate range.

Agree, but an annotation will still only have a valid meaning in the x-dimension, correct? My main worry was that If we enable something that is not well defined to start with, it will become harder in the future to limit the scope of it.
philippjfr: >Fair point. But in this case, the labels are different in the code but the same in the figure.

I thought they were merely overlapping but I think this must be a consequence of them being stored in a dictionary.

>Agree, but an annotation will still only have a valid meaning in the x-dimension, correct? My main worry was that If we enable something that is not well defined to start with, it will become harder in the future to limit the scope of it.

The aim will be to broaden the scope in future, so I'm not worried about this.
maximlt: I have made changes based on your comments @Hoxbro , thanks for the thorough review!

On the edge cases you shared, the error displayed in the first one is indeed a little terse 🙃 although the error is pretty obvious in your example. After your discussion with Philipp, let me know if you want me to investigate how to raise a more informative error (I guess I'll have to traverse all the elements to check they have unique labels, maybe there's an API for that?). The "same subcordinate_y case" is a feature request to me, please somehow explode the ticks when they overlap. The last case with "only some curves have enabled subcoordinate_y" will hopefully be supported soon. I'm thinking that users could provide a range for each element that doesn't have to be between 0 and 1, we'd traverse all the elements and collect these ranges to compute v0 and v1.

On your comment on whether multi axes are subcoordinates are supported together, I believe that is not the case (haven't tried though!) and that it's very much on the edge.
maximlt: @Hoxbro @philippjfr is there any remaining changes required for this PR to be merged?
Hoxbro: Two things for me:

Is it too tedious to add a way to detect if the labels are the same? If it is please open a small issue so we can keep track of it. 


The other thing is the twin axis definitely looks wrong:
``` python
import holoviews as hv
import numpy as np

hv.extension("bokeh")

x = np.linspace(0, 10*np.pi)

curves = [
    hv.Curve((x + i*np.pi/2, np.sin(x)), label=f'Label {i}').opts(subcoordinate_y=True, subcoordinate_scale=1.2)
    for i in range(3)
]

curves2 = [
    hv.Curve((x + i*np.pi/2, np.sin(x)), label=f'Label {i}', vdims=["twin"]).opts(subcoordinate_y=True, subcoordinate_scale=1.2)
    for i in range(3, 6)
]

hv.Overlay(curves + curves2).opts(show_legend=False, multi_y=True)
```
Outputs (no labels, and only show the first curve)
![image](https://github.com/holoviz/holoviews/assets/19758978/a6307be6-ac31-42ce-bf88-db55fc95b771)


Zooming in on the y-axis, and then zoom on the plot you get:
![image](https://github.com/holoviz/holoviews/assets/19758978/79baf173-cac2-4274-b8c4-52c7f3cb199d)

Can we make a check in `_subcoord_overlaid` that `multi_y` and raise an error if both are enabled?


maximlt: > Is it too tedious to add a way to detect if the labels are the same?

Not sure it's tedious, I just haven't been able to find the internal API to do it, I'll look for it but welcome any pointer! I hope having to traverse the objects won't be too costly, since an aspect of this review was performance optimizations.

> Can we make a check in _subcoord_overlaid that multi_y and raise an error if both are enabled?

I'll add a check. Not sure the property is the right place to make that check though, I guess those kind of checks should all be made in the same place, at some point in the initialization phase, I just have to find where this is.

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 524281a25..0131ab47a 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -17,7 +17,7 @@ repos:
         exclude: \.min\.js$
       - id: trailing-whitespace
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.0.290
+    rev: v0.0.291
     hooks:
       - id: ruff
         args: [holoviews]
diff --git a/examples/gallery/demos/bokeh/eeg_viewer.ipynb b/examples/gallery/demos/bokeh/eeg_viewer.ipynb
new file mode 100644
index 000000000..01f3cb98c
--- /dev/null
+++ b/examples/gallery/demos/bokeh/eeg_viewer.ipynb
@@ -0,0 +1,172 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "549b47a4",
+   "metadata": {},
+   "source": [
+    "This example demonstrates advanced visualization techniques using HoloViews with the Bokeh plotting backend. You'll learn how to:\n",
+    "\n",
+    "1. Display multiple timeseries in a single plot using `subcoordinate_y`.\n",
+    "2. Create and link a minimap to the main plot with `RangeToolLink`.\n",
+    "\n",
+    "Specifically, we'll simulate [Electroencephalography](https://en.wikipedia.org/wiki/Electroencephalography) (EEG) data, plot it, and then create a minimap based on the [z-score](https://en.wikipedia.org/wiki/Standard_score) of the data for easier navigation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8109537b-5fba-4f07-aba4-91a56f7e95c7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import holoviews as hv\n",
+    "from bokeh.models import HoverTool\n",
+    "from holoviews.plotting.links import RangeToolLink\n",
+    "from scipy.stats import zscore\n",
+    "\n",
+    "hv.extension('bokeh')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1c95f241-2314-42b0-b6cb-2c0baf332686",
+   "metadata": {},
+   "source": [
+    "## Generating EEG data\n",
+    "\n",
+    "Let's start by simulating some EEG data. We'll create a timeseries for each channel using sine waves with varying frequencies."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5f4a9dbe",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "N_CHANNELS = 10\n",
+    "N_SECONDS = 5\n",
+    "SAMPLING_RATE = 200\n",
+    "INIT_FREQ = 2  # Initial frequency in Hz\n",
+    "FREQ_INC = 5  # Frequency increment\n",
+    "AMPLITUDE = 1\n",
+    "\n",
+    "# Generate time and channel labels\n",
+    "total_samples = N_SECONDS * SAMPLING_RATE\n",
+    "time = np.linspace(0, N_SECONDS, total_samples)\n",
+    "channels = [f'EEG {i}' for i in range(N_CHANNELS)]\n",
+    "\n",
+    "# Generate sine wave data\n",
+    "data = np.array([AMPLITUDE * np.sin(2 * np.pi * (INIT_FREQ + i * FREQ_INC) * time)\n",
+    "                     for i in range(N_CHANNELS)])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ec9e71b8-a995-4c0f-bdbb-5d148d8fa138",
+   "metadata": {},
+   "source": [
+    "## Visualizing EEG Data\n",
+    "\n",
+    "Next, let's dive into visualizing the EEG data. We construct each timeseries using a `Curve` element, assigning it a `label` and setting `subcoordinate_y=True`. All these curves are then aggregated into a list, which serves as the input for an `Overlay` element. Rendering this `Overlay` produces a plot where the timeseries are stacked vertically.\n",
+    "\n",
+    "Additionally, we'll enhance user interaction by implementing a custom hover tool. This will display key information—channel, time, and amplitude—when you hover over any of the curves."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9476769f-3935-4236-b010-1511d1a1e77f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hover = HoverTool(tooltips=[\n",
+    "    (\"Channel\", \"@channel\"),\n",
+    "    (\"Time\", \"$x s\"),\n",
+    "    (\"Amplitude\", \"$y µV\")\n",
+    "])\n",
+    "\n",
+    "channel_curves = []\n",
+    "for channel, channel_data in zip(channels, data):\n",
+    "    ds = hv.Dataset((time, channel_data, channel), [\"Time\", \"Amplitude\", \"channel\"])\n",
+    "    curve = hv.Curve(ds, \"Time\", [\"Amplitude\", \"channel\"], label=channel)\n",
+    "    curve.opts(\n",
+    "        subcoordinate_y=True, color=\"black\", line_width=1, tools=[hover],\n",
+    "    )\n",
+    "    channel_curves.append(curve)\n",
+    "\n",
+    "eeg = hv.Overlay(channel_curves, kdims=\"Channel\").opts(\n",
+    "    xlabel=\"Time (s)\", ylabel=\"Channel\", show_legend=False, aspect=3, responsive=True,\n",
+    ")\n",
+    "eeg"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b4f603e2-039d-421a-ba9a-ed9e77efab99",
+   "metadata": {},
+   "source": [
+    "## Creating the Minimap\n",
+    "\n",
+    "A minimap can provide a quick overview of the data and help you navigate through it. We'll compute the z-score for each channel and represent it as an image; the z-score will normalize the data and bring out the patterns more clearly. To enable linking in the next step between the EEG `Overlay` and the minimap `Image`, we ensure they share the same y-axis range."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "40fa2198-c3b5-41e1-944f-f8b812612168",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y_positions = range(N_CHANNELS)\n",
+    "yticks = [(i , ich) for i, ich in enumerate(channels)]\n",
+    "\n",
+    "z_data = zscore(data, axis=1)\n",
+    "\n",
+    "minimap = hv.Image((time, y_positions , z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n",
+    "minimap = minimap.opts(\n",
+    "    cmap=\"RdBu_r\", xlabel='Time (s)', alpha=.5, yticks=[yticks[0], yticks[-1]],\n",
+    "    height=150, responsive=True, default_tools=[], clim=(-z_data.std(), z_data.std())\n",
+    ")\n",
+    "minimap"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a5b77970-342f-4428-bd1c-4dbef1e6a2b5",
+   "metadata": {},
+   "source": [
+    "## Building the dashboard\n",
+    "\n",
+    "Finally, we use [`RangeToolLink`](../../../user_guide/Linking_Plots.ipynb) to connect the minimap `Image` and the EEG `Overlay`, setting bounds for the initial viewable area. Once the plots are linked and assembled into a unified dashboard, you can interact with it. Experiment by dragging the selection box on the minimap or resizing it by clicking and dragging its edges."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "260489eb-2dbf-4c88-ba83-dd1cba0e547b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "RangeToolLink(\n",
+    "    minimap, eeg, axes=[\"x\", \"y\"],\n",
+    "    boundsx=(None, 2), boundsy=(None, 6.5)\n",
+    ")\n",
+    "\n",
+    "dashboard = (eeg + minimap).opts(merge_tools=False).cols(1)\n",
+    "dashboard"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/user_guide/Customizing_Plots.ipynb b/examples/user_guide/Customizing_Plots.ipynb
index 24c5b8d91..5fea84edb 100644
--- a/examples/user_guide/Customizing_Plots.ipynb
+++ b/examples/user_guide/Customizing_Plots.ipynb
@@ -745,6 +745,34 @@
     "\n",
     "Note that as of HoloViews 1.17.0, `multi_y` does not have streaming plot support, extra axis labels are not dynamic and only the `RangeXY` linked stream is aware of additional y-axes."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Subcoordinate y-axis\n",
+    "*(Available in HoloViews >= 1.18)*\n",
+    "\n",
+    "HoloViews enables you to create overlays where each element has its own distinct y-axis subcoordinate system. To activate this feature, set the `subcoordinate_y` keyword to True for **each** overlay element; the default is False. When using `subcoordinate_y=True`, setting a `label` for each element is required for proper rendering and identification.This will automatically distribute overlay elements along the y-axis.\n",
+    "\n",
+    "For more fine-grained control over y-axis positioning, you can specify a numerical 2-tuple for subcoordinate_y with values ranging from 0 to 1. Additionally, the `subcoordinate_scale` keyword, which defaults to 1, allows you to adjust the vertical scale of each element. This option is only applicable when `subcoordinate_y=True`. For example, setting a single Curve's `subcoordinate_scale` to 2 will result in it overlapping 50% with its adjacent elements."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = np.linspace(0, 10*np.pi)\n",
+    "\n",
+    "curves = [\n",
+    "    hv.Curve((x + i*np.pi/2, np.sin(x)), label=f'Line {i}').opts(subcoordinate_y=True, subcoordinate_scale=1.2)\n",
+    "    for i in range(3)\n",
+    "]\n",
+    "\n",
+    "hv.Overlay(curves).opts(show_legend=False)"
+   ]
   }
  ],
  "metadata": {
diff --git a/holoviews/plotting/bokeh/element.py b/holoviews/plotting/bokeh/element.py
index 079311b93..c1d20faa2 100644
--- a/holoviews/plotting/bokeh/element.py
+++ b/holoviews/plotting/bokeh/element.py
@@ -154,6 +154,14 @@ class ElementPlot(BokehPlot, GenericElementPlot):
        elements and the overlay container, allowing customization on a
        per-axis basis.""")
 
+    subcoordinate_y = param.ClassSelector(default=False, class_=(bool, tuple), doc="""
+       Enables sub-coordinate systems for this plot. Accepts also a numerical
+       two-tuple that must be a range between 0 and 1, the plot will be
+       rendered on this vertical range of the axis.""")
+
+    subcoordinate_scale = param.Number(default=1, bounds=(0, None), inclusive_bounds=(False, True), doc="""
+       Scale factor for subcoordinate ranges to control the level of overlap.""")
+
     responsive = param.ObjectSelector(default=False, objects=[False, True, 'width', 'height'])
 
     fontsize = param.Parameter(default={'title': '12pt'}, allow_None=True,  doc="""
@@ -229,6 +237,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
     _stream_data = True
 
     def __init__(self, element, plot=None, **params):
+        self._subcoord_standalone_ = None
         self.current_ranges = None
         super().__init__(element, **params)
         self.handles = {} if plot is None else self.handles['plot']
@@ -401,6 +410,21 @@ class ElementPlot(BokehPlot, GenericElementPlot):
 
         return dim_range
 
+    @property
+    def _subcoord_overlaid(self):
+        """
+        Indicates when the context is a subcoordinate plot, either from within
+        the overlay rendering or one of its subplots. Used to skip code paths
+        when rendering an element outside of an overlay.
+        """
+        if self._subcoord_standalone_ is not None:
+            return self._subcoord_standalone_
+        self._subcoord_standalone_ = (
+            (isinstance(self, OverlayPlot) and self.subcoordinate_y) or
+            (not isinstance(self, OverlayPlot) and self.overlaid and self.subcoordinate_y)
+        )
+        return self._subcoord_standalone_
+
     def _axis_props(self, plots, subplots, element, ranges, pos, *, dim=None,
                     range_tags_extras=[], extra_range_name=None):
 
@@ -416,7 +440,10 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         # a synthetic dimension and exclusively supported for y-axes.
         if pos == 1 and dim:
             dims = [dim]
-            v0, v1 = util.max_range([elrange.get(dim, {'combined': (None, None)})['combined'] for elrange in ranges.values()])
+            v0, v1 = util.max_range([
+                elrange.get(dim.name, {'combined': (None, None)})['combined']
+                for elrange in ranges.values()
+            ])
             axis_label = str(dim)
             specs = ((dim.name, dim.label, dim.unit),)
         else:
@@ -428,7 +455,17 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                 l, b, r, t = self.get_extents(range_el, ranges)
             if self.invert_axes:
                 l, b, r, t = b, l, t, r
-            v0, v1 = (l, r) if pos == 0 else (b, t)
+            if pos == 1 and self._subcoord_overlaid:
+                if isinstance(self.subcoordinate_y, bool):
+                    offset = self.subcoordinate_scale / 2.
+                    # This sum() is equal to n+1, n being the number of elements contained
+                    # in the overlay with subcoordinate_y=True, as the traversal goes through
+                    # the root overlay that has subcoordinate_y=True too since it's propagated.
+                    v0, v1 = 0-offset, sum(self.traverse(lambda p: p.subcoordinate_y))-2+offset
+                else:
+                    v0, v1 = 0, 1
+            else:
+                v0, v1 = (l, r) if pos == 0 else (b, t)
 
             axis_dims = list(self._get_axis_dims(el))
             if self.invert_axes:
@@ -452,7 +489,9 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                 dims = dims[:2][::-1]
 
         categorical = any(self.traverse(lambda plot: plot._categorical))
-        if dims is not None and any(dim.name in ranges and 'factors' in ranges[dim.name] for dim in dims):
+        if self.subcoordinate_y:
+            categorical = False
+        elif dims is not None and any(dim.name in ranges and 'factors' in ranges[dim.name] for dim in dims):
             categorical = True
         else:
             categorical = any(isinstance(v, (str, bytes)) for v in (v0, v1))
@@ -474,12 +513,12 @@ class ElementPlot(BokehPlot, GenericElementPlot):
 
         norm_opts = self.lookup_options(el, 'norm').options
         shared_name = extra_range_name or ('x-main-range' if pos == 0 else 'y-main-range')
-        if plots and self.shared_axes and not norm_opts.get('axiswise', False):
+        if plots and self.shared_axes and not norm_opts.get('axiswise', False) and not dim:
             dim_range = self._shared_axis_range(plots, specs, range_type, axis_type, pos)
             if dim_range:
                 self._shared[shared_name] = True
 
-        if self._shared.get(shared_name):
+        if self._shared.get(shared_name) and not dim:
             pass
         elif categorical:
             axis_type = 'auto'
@@ -507,16 +546,24 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             axpos0, axpos1 = 'left', 'right'
 
         ax_specs, yaxes, dimensions = {}, {}, {}
+        subcoordinate_axes = 0
         for el, sp in zip(element, self.subplots.values()):
             ax_dims = sp._get_axis_dims(el)[:2]
             if sp.invert_axes:
                 ax_dims[::-1]
             yd = ax_dims[1]
-            dimensions[yd.name] = yd
             opts = el.opts.get('plot', backend='bokeh').kwargs
             if not isinstance(yd, Dimension) or yd.name in yaxes:
                 continue
-            yaxes[yd.name] = {
+            if self._subcoord_overlaid:
+                if opts.get('subcoordinate_y') is None:
+                    continue
+                ax_name = el.label
+                subcoordinate_axes += 1
+            else:
+                ax_name = yd.name
+            dimensions[ax_name] = yd
+            yaxes[ax_name] = {
                 'position': opts.get('yaxis', axpos1 if len(yaxes) else axpos0),
                 'autorange': opts.get('autorange', None),
                 'logx': opts.get('logx', False),
@@ -528,11 +575,14 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                 'fontsize': {
                     'axis_label_text_font_size': sp._fontsize('ylabel').get('fontsize'),
                     'major_label_text_font_size': sp._fontsize('yticks').get('fontsize')
-                }
+                },
+                'subcoordinate_y': subcoordinate_axes-1 if self._subcoord_overlaid else None
             }
 
         for ydim, info in yaxes.items():
             range_tags_extras = {'invert_yaxis': info['invert_yaxis']}
+            if info['subcoordinate_y'] is not None:
+                range_tags_extras['subcoordinate_y'] = info['subcoordinate_y']
             if info['autorange'] == 'y':
                 range_tags_extras['autorange'] = True
                 lowerlim, upperlim = info['ylim'][0], info['ylim'][1]
@@ -542,7 +592,6 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                     range_tags_extras['y-upperlim'] = upperlim
             else:
                 range_tags_extras['autorange'] = False
-
             ax_props = self._axis_props(
                 plots, subplots, element, ranges, pos=1, dim=dimensions[ydim],
                 range_tags_extras=range_tags_extras,
@@ -582,18 +631,23 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             else:
                 range_tags_extras['autorange'] = False
             axis_specs['y']['y'] = self._axis_props(
-                plots, subplots, element, ranges, pos=1, range_tags_extras = range_tags_extras
+                plots, subplots, element, ranges, pos=1, range_tags_extras=range_tags_extras
             ) + (self.yaxis, {})
 
+        if self._subcoord_overlaid:
+            _, extra_axis_specs = self._create_extra_axes(plots, subplots, element, ranges)
+            axis_specs['y'].update(extra_axis_specs)
+
         properties, axis_props = {}, {'x': {}, 'y': {}}
         for axis, axis_spec in axis_specs.items():
             for (axis_dim, (axis_type, axis_label, axis_range, axis_position, fontsize)) in axis_spec.items():
                 scale = get_scale(axis_range, axis_type)
                 if f'{axis}_range' in properties:
                     properties[f'extra_{axis}_ranges'] = extra_ranges = properties.get(f'extra_{axis}_ranges', {})
-                    properties[f'extra_{axis}_scales'] = extra_scales = properties.get(f'extra_{axis}_scales', {})
                     extra_ranges[axis_dim] = axis_range
-                    extra_scales[axis_dim] = scale
+                    if not self.subcoordinate_y:
+                        properties[f'extra_{axis}_scales'] = extra_scales = properties.get(f'extra_{axis}_scales', {})
+                        extra_scales[axis_dim] = scale
                 else:
                     properties[f'{axis}_range'] = axis_range
                     properties[f'{axis}_scale'] = scale
@@ -642,6 +696,10 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         fig.xaxis[0].update(**axis_props['x'])
         fig.yaxis[0].update(**axis_props['y'])
 
+        # Do not add the extra axes to the layout if subcoordinates are used
+        if self._subcoord_overlaid:
+            return fig
+
         multi_ax = 'x' if self.invert_axes else 'y'
         for axis_dim, range_obj in properties.get(f'extra_{multi_ax}_ranges', {}).items():
             axis_type, axis_label, _, axis_position, fontsize = axis_specs[multi_ax][axis_dim]
@@ -805,6 +863,17 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                 axis_props['ticker'] = FixedTicker(ticks=ticks)
                 if labels is not None:
                     axis_props['major_label_overrides'] = dict(zip(ticks, labels))
+            elif self._subcoord_overlaid and axis == 'y':
+                ticks, labels = [], []
+                for i, (el, sp) in enumerate(zip(self.current_frame, self.subplots.values())):
+                    if not sp.subcoordinate_y:
+                        continue
+                    ycenter = i if isinstance(sp.subcoordinate_y, bool) else 0.5 * sum(sp.subcoordinate_y)
+                    ticks.append(ycenter)
+                    labels.append(el.label)
+                axis_props['ticker'] = FixedTicker(ticks=ticks)
+                if labels is not None:
+                    axis_props['major_label_overrides'] = dict(zip(ticks, labels))
         formatter = self.xformatter if axis == 'x' else self.yformatter
         if formatter:
             formatter = wrap_formatter(formatter, axis)
@@ -956,6 +1025,9 @@ class ElementPlot(BokehPlot, GenericElementPlot):
 
         self._update_main_ranges(element, x_range, y_range, ranges)
 
+        if self._subcoord_overlaid:
+            return
+
         # ALERT: stream handling not handled
         streaming = False
         multi_dim = 'x' if self.invert_axes else 'y'
@@ -1002,8 +1074,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                     and (framewise or streaming))
                    or xfactors is not None)
         yupdate = ((not (self.model_changed(x_range) or self.model_changed(plot))
-                    and (framewise or streaming))
-                   or yfactors is not None)
+                    and (framewise or streaming) or yfactors is not None) and not self.subcoordinate_y)
 
         options = self._traverse_options(element, 'plot', ['width', 'height'], defaults=False)
         fixed_width = (self.frame_width or options.get('width'))
@@ -1115,7 +1186,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         if not self.drawn or xupdate:
             self._update_range(x_range, l, r, xfactors, self.invert_xaxis,
                                self._shared['x-main-range'], self.logx, streaming)
-        if not self.drawn or yupdate:
+        if not (self.drawn or self.subcoordinate_y) or yupdate:
             self._update_range(
                 y_range, b, t, yfactors, self._get_tag(y_range, 'invert_yaxis'),
                 self._shared['y-main-range'], self.logy, streaming
@@ -1394,13 +1465,28 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         if 'legend_field' in properties and 'legend_label' in properties:
             del properties['legend_label']
 
-        if self.handles['x_range'].name in plot.extra_x_ranges:
-            properties['x_range_name'] = self.handles['y_range'].name
-        if self.handles['y_range'].name in plot.extra_y_ranges:
+        if self.handles['x_range'].name in plot.extra_x_ranges and not self.subcoordinate_y:
+            properties['x_range_name'] = self.handles['x_range'].name
+        if self.handles['y_range'].name in plot.extra_y_ranges and not self.subcoordinate_y:
             properties['y_range_name'] = self.handles['y_range'].name
 
         if "name" not in properties:
             properties["name"] = properties.get("legend_label") or properties.get("legend_field")
+
+        if self._subcoord_overlaid:
+            y_source_range = self.handles['y_range']
+            if isinstance(self.subcoordinate_y, bool):
+                center = y_source_range.tags[1]['subcoordinate_y']
+                offset = self.subcoordinate_scale/2.
+                ytarget_range = dict(start=center-offset, end=center+offset)
+            else:
+                ytarget_range = dict(start=self.subcoordinate_y[0], end=self.subcoordinate_y[1])
+            plot = plot.subplot(
+                x_source=plot.x_range,
+                x_target=plot.x_range,
+                y_source=y_source_range,
+                y_target=Range1d(**ytarget_range),
+            )
         renderer = getattr(plot, plot_method)(**dict(properties, **mapping))
         return renderer, renderer.glyph
 
@@ -1773,6 +1859,9 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             axes, plot_ranges = self._find_axes(plot, element)
             self.handles['xaxis'], self.handles['yaxis'] = axes
             self.handles['x_range'], self.handles['y_range'] = plot_ranges
+            if self._subcoord_overlaid:
+                if style_element.label in plot.extra_y_ranges:
+                    self.handles['y_range'] = plot.extra_y_ranges.pop(style_element.label)
         self.handles['plot'] = plot
 
         if self.autorange:
@@ -2526,7 +2615,8 @@ class OverlayPlot(GenericOverlayPlot, LegendPlot):
                           'xformatter', 'yformatter', 'active_tools',
                           'min_height', 'max_height', 'min_width', 'min_height',
                           'margin', 'aspect', 'data_aspect', 'frame_width',
-                          'frame_height', 'responsive', 'fontscale']
+                          'frame_height', 'responsive', 'fontscale', 'subcoordinate_y',
+                          'subcoordinate_scale']
 
     def __init__(self, overlay, **kwargs):
         self._multi_y_propagation = self.lookup_options(overlay, 'plot').options.get('multi_y', False)
@@ -2750,6 +2840,22 @@ class OverlayPlot(GenericOverlayPlot, LegendPlot):
         return super()._get_axis_dims(element)
 
     def initialize_plot(self, ranges=None, plot=None, plots=None):
+        if self.multi_y and self.subcoordinate_y:
+            raise ValueError('multi_y and subcoordinate_y are not supported together.')
+        if self.subcoordinate_y:
+            labels = self.hmap.last.traverse(lambda x: x.label, [
+                lambda el: isinstance(el, Element) and el.opts.get('plot').kwargs.get('subcoordinate_y', False)
+            ])
+            if any(not label for label in labels):
+                raise ValueError(
+                    'Every element wrapped in a subcoordinate_y overlay must have '
+                    'a label.'
+                )
+            if len(set(labels)) == 1:
+                raise ValueError(
+                    'Elements wrapped in a subcoordinate_y overlay must all have '
+                    'a unique label.'
+                )
         key = util.wrap_tuple(self.hmap.last_key)
         nonempty = [(k, el) for k, el in self.hmap.data.items() if el]
         if not nonempty:


Test patch:
diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index bdfbdadd6..553bc7264 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -33,7 +33,7 @@ jobs:
     name: Run pre-commit
     runs-on: 'ubuntu-latest'
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pre-commit@v0.1a16
+      - uses: holoviz-dev/holoviz_tasks/pre-commit@v0.1a17
 
   unit_test_suite:
     name: Unit tests on Python ${{ matrix.python-version }}, ${{ matrix.os }}
@@ -57,7 +57,7 @@ jobs:
       DESC: "Python ${{ matrix.python-version }}, ${{ matrix.os }} unit tests"
       PYTHON_VERSION: ${{ matrix.python-version }}
     steps:
-      - uses: holoviz-dev/holoviz_tasks/install@v0.1a16
+      - uses: holoviz-dev/holoviz_tasks/install@v0.1a17
         with:
           name: unit_test_suite
           python-version: ${{ matrix.python-version }}
@@ -104,7 +104,7 @@ jobs:
       # it as one of the sources.
       PYCTDEV_SELF_CHANNEL: "pyviz/label/dev"
     steps:
-      - uses: holoviz-dev/holoviz_tasks/install@v0.1a16
+      - uses: holoviz-dev/holoviz_tasks/install@v0.1a17
         with:
           name: ui_test_suite
           python-version: ${{ matrix.python-version }}
@@ -144,7 +144,7 @@ jobs:
       DESC: "Python ${{ matrix.python-version }}, ${{ matrix.os }} core tests"
       PYTHON_VERSION: ${{ matrix.python-version }}
     steps:
-      - uses: holoviz-dev/holoviz_tasks/install@v0.1a16
+      - uses: holoviz-dev/holoviz_tasks/install@v0.1a17
         with:
           name: core_test_suite
           python-version: ${{ matrix.python-version }}
diff --git a/holoviews/tests/plotting/bokeh/test_subcoordy.py b/holoviews/tests/plotting/bokeh/test_subcoordy.py
new file mode 100644
index 000000000..d36ae6447
--- /dev/null
+++ b/holoviews/tests/plotting/bokeh/test_subcoordy.py
@@ -0,0 +1,190 @@
+import numpy as np
+import pytest
+
+from holoviews.core import Overlay
+from holoviews.element import Curve
+from holoviews.element.annotation import VSpan
+
+from .test_plot import TestBokehPlot, bokeh_renderer
+
+
+class TestSubcoordinateY(TestBokehPlot):
+
+    # With subcoordinate_y set to True
+
+    def test_bool_base(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        plot = bokeh_renderer.get_plot(overlay)
+        # subcoordinate_y is propagated to the overlay
+        assert plot.subcoordinate_y is True
+        # the figure has only one yaxis
+        assert len(plot.state.yaxis) == 1
+        # the overlay has two subplots
+        assert len(plot.subplots) == 2
+        assert ('Curve', 'Data_0') in plot.subplots
+        assert ('Curve', 'Data_1') in plot.subplots
+        # the range per subplots are correctly computed
+        sp1 = plot.subplots[('Curve', 'Data_0')]
+        assert sp1.handles['glyph_renderer'].coordinates.y_target.start == -0.5
+        assert sp1.handles['glyph_renderer'].coordinates.y_target.end == 0.5
+        sp2 = plot.subplots[('Curve', 'Data_1')]
+        assert sp2.handles['glyph_renderer'].coordinates.y_target.start == 0.5
+        assert sp2.handles['glyph_renderer'].coordinates.y_target.end == 1.5
+        # y_range is correctly computed
+        assert plot.handles['y_range'].start == -0.5
+        assert plot.handles['y_range'].end == 1.5
+        # extra_y_range is empty
+        assert plot.handles['extra_y_ranges'] == {}
+        # the ticks show the labels
+        assert plot.state.yaxis.ticker.ticks == [0, 1]
+        assert plot.state.yaxis.major_label_overrides == {0: 'Data 0', 1: 'Data 1'}
+
+    def test_bool_scale(self):
+        test_data = [
+            (0.5, (-0.25, 0.25), (0.75, 1.25), (-0.25, 1.25)),
+            (1, (-0.5, 0.5), (0.5, 1.5), (-0.5, 1.5)),
+            (2, (-1, 1), (0, 2), (-1, 2)),
+            (5, (-2.5, 2.5), (-1.5, 3.5), (-2.5, 3.5)),
+        ]
+        for scale, ytarget1, ytarget2, ytarget in test_data:
+            overlay = Overlay([
+                Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True, subcoordinate_scale=scale)
+                for i in range(2)
+            ])
+            plot = bokeh_renderer.get_plot(overlay)
+            # the range per subplots are correctly computed
+            sp1 = plot.subplots[('Curve', 'Data_0')]
+            assert sp1.handles['glyph_renderer'].coordinates.y_target.start == ytarget1[0]
+            assert sp1.handles['glyph_renderer'].coordinates.y_target.end == ytarget1[1]
+            sp2 = plot.subplots[('Curve', 'Data_1')]
+            assert sp2.handles['glyph_renderer'].coordinates.y_target.start == ytarget2[0]
+            assert sp2.handles['glyph_renderer'].coordinates.y_target.end == ytarget2[1]
+            # y_range is correctly computed
+            assert plot.handles['y_range'].start == ytarget[0]
+            assert plot.handles['y_range'].end == ytarget[1]
+
+    def test_no_label(self):
+        overlay = Overlay([Curve(range(10)).opts(subcoordinate_y=True) for i in range(2)])
+        with pytest.raises(
+            ValueError,
+            match='Every element wrapped in a subcoordinate_y overlay must have a label'
+        ):
+            bokeh_renderer.get_plot(overlay)
+
+    def test_overlaid_without_label_no_error(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        with_span = overlay * VSpan(1, 2)
+        bokeh_renderer.get_plot(with_span)
+
+    def test_custom_ylabel(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        overlay.opts(ylabel='Y label')
+        plot = bokeh_renderer.get_plot(overlay)
+        # the figure axis has the label set
+        assert plot.state.yaxis.axis_label == 'Y label'
+        # the ticks show the labels
+        assert plot.state.yaxis.ticker.ticks == [0, 1]
+        assert plot.state.yaxis.major_label_overrides == {0: 'Data 0', 1: 'Data 1'}
+
+    def test_legend_label(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        plot = bokeh_renderer.get_plot(overlay)
+        legend_labels = [l.label['value'] for l in plot.state.legend[0].items]
+        # the legend displays the labels
+        assert legend_labels == ['Data 0', 'Data 1']
+
+    def test_shared_multi_axes(self):
+        overlay1 = Overlay([Curve(np.arange(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        overlay2 = Overlay([Curve(np.arange(10) + 5, label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+
+        plot = bokeh_renderer.get_plot(overlay1 + overlay2)
+
+        oplot1 = plot.subplots[(0, 0)].subplots['main']
+        oplot2 = plot.subplots[(0, 1)].subplots['main']
+        assert (oplot1.handles['y_range'].start, oplot1.handles['y_range'].end) == (-0.5, 1.5)
+        assert oplot1.handles['extra_y_ranges'] == {}
+        assert (oplot2.handles['y_range'].start, oplot2.handles['y_range'].end) == (-0.5, 1.5)
+        assert oplot2.handles['extra_y_ranges'] == {}
+
+    def test_invisible_yaxis(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        overlay.opts(yaxis=None)
+        plot = bokeh_renderer.get_plot(overlay)
+        assert not plot.state.yaxis.visible
+
+    def test_axis_labels(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        plot = bokeh_renderer.get_plot(overlay)
+
+        assert plot.state.xaxis.axis_label == 'x'
+        assert plot.state.yaxis.axis_label == 'y'
+
+    def test_only_x_axis_labels(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        overlay.opts(labelled=['x'])
+        plot = bokeh_renderer.get_plot(overlay)
+
+        assert plot.state.xaxis.axis_label == 'x'
+        assert plot.state.yaxis.axis_label == ''
+
+    def test_none_x_axis_labels(self):
+        overlay = Overlay([Curve(range(10), vdims=['A'], label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        plot = bokeh_renderer.get_plot(overlay)
+
+        assert plot.state.xaxis.axis_label == 'x'
+        assert plot.state.yaxis.axis_label == 'A'
+
+    # With subcoordinate_y set to a range
+
+    def test_range_base(self):
+        overlay = Overlay([
+            Curve(range(10), label='Data 0').opts(subcoordinate_y=(0, 0.5)),
+            Curve(range(10), label='Data 1').opts(subcoordinate_y=(0.5, 1)),
+        ])
+        plot = bokeh_renderer.get_plot(overlay)
+        # subcoordinate_y is propagated to the overlay, just the first one though :(
+        assert plot.subcoordinate_y == (0, 0.5)
+        # the figure has only one yaxis
+        assert len(plot.state.yaxis) == 1
+        # the overlay has two subplots
+        assert len(plot.subplots) == 2
+        assert ('Curve', 'Data_0') in plot.subplots
+        assert ('Curve', 'Data_1') in plot.subplots
+        # the range per subplots are correctly computed
+        sp1 = plot.subplots[('Curve', 'Data_0')]
+        assert sp1.handles['glyph_renderer'].coordinates.y_target.start == 0
+        assert sp1.handles['glyph_renderer'].coordinates.y_target.end == 0.5
+        sp2 = plot.subplots[('Curve', 'Data_1')]
+        assert sp2.handles['glyph_renderer'].coordinates.y_target.start == 0.5
+        assert sp2.handles['glyph_renderer'].coordinates.y_target.end == 1
+        # y_range is correctly computed
+        assert plot.handles['y_range'].start == 0
+        assert plot.handles['y_range'].end == 1
+        # extra_y_range is empty
+        assert plot.handles['extra_y_ranges'] == {}
+        # the ticks show the labels
+        assert plot.state.yaxis.ticker.ticks == [0.25, 0.75]
+        assert plot.state.yaxis.major_label_overrides == {0.25: 'Data 0', 0.75: 'Data 1'}
+
+    def test_plot_standalone(self):
+        standalone = Curve(range(10), label='Data 0').opts(subcoordinate_y=True)
+        plot = bokeh_renderer.get_plot(standalone)
+        assert (plot.state.x_range.start, plot.state.x_range.end) == (0, 9)
+        assert (plot.state.y_range.start, plot.state.y_range.end) == (0, 9)
+
+    def test_multi_y_error(self):
+        overlay = Overlay([Curve(range(10), label=f'Data {i}').opts(subcoordinate_y=True) for i in range(2)])
+        overlay.opts(multi_y=True)
+        with pytest.raises(
+            ValueError,
+            match='multi_y and subcoordinate_y are not supported together'
+        ):
+            bokeh_renderer.get_plot(overlay)
+
+    def test_same_label_error(self):
+        overlay = Overlay([Curve(range(10), label='Same').opts(subcoordinate_y=True) for _ in range(2)])
+        with pytest.raises(
+            ValueError,
+            match='Elements wrapped in a subcoordinate_y overlay must all have a unique label',
+        ):
+            bokeh_renderer.get_plot(overlay)