Repository: holoviz/holoviews

Issue:
ImageStack doesn't allow non-square data
#### ALL software version info
holoviews 1.18.0
bokeh 3.3.0

#### Description of expected behavior and the observed behavior
Looks like the x/y dimensions got swapped somewhere.

```python
hv.ImageStack(np.random.random((10,10,3))) # works
hv.ImageStack(np.random.random((10,15,3))) # breaks
```

gives the following traceback

<details>

```
---------------------------------------------------------------------------
DataError                                 Traceback (most recent call last)
Cell In[107], line 1
----> 1 hv.ImageStack(np.random.random((10,15,3)))

File ~/micromamba/envs/image_analysis/lib/python3.11/site-packages/holoviews/element/raster.py:547, in ImageStack.__init__(self, data, kdims, vdims, **params)
    545     elif isinstance(data, dict):
    546         vdims = [Dimension(key) for key in data.keys() if key not in self.kdims]
--> 547 super().__init__(data, kdims=kdims, vdims=vdims, **params)

File ~/micromamba/envs/image_analysis/lib/python3.11/site-packages/holoviews/element/raster.py:279, in Image.__init__(self, data, kdims, vdims, bounds, extents, xdensity, ydensity, rtol, **params)
    276 else:
    277     params['rtol'] = config.image_rtol
--> 279 Dataset.__init__(self, data, kdims=kdims, vdims=vdims, extents=extents, **params)
    280 if not self.interface.gridded:
    281     raise DataError("{} type expects gridded data, {} is columnar. "
    282                     "To display columnar data as gridded use the HeatMap "
    283                     "element or aggregate the data (e.g. using rasterize "
    284                     "or np.histogram2d).".format(type(self).__name__, self.interface.__name__))

File ~/micromamba/envs/image_analysis/lib/python3.11/site-packages/holoviews/core/data/__init__.py:329, in Dataset.__init__(self, data, kdims, vdims, **kwargs)
    326 kdims, vdims = kwargs.get('kdims'), kwargs.get('vdims')
    328 validate_vdims = kwargs.pop('_validate_vdims', True)
--> 329 initialized = Interface.initialize(type(self), data, kdims, vdims,
    330                                    datatype=kwargs.get('datatype'))
    331 (data, self.interface, dims, extra_kws) = initialized
    332 super().__init__(data, **dict(kwargs, **dict(dims, **extra_kws)))

File ~/micromamba/envs/image_analysis/lib/python3.11/site-packages/holoviews/core/data/interface.py:253, in Interface.initialize(cls, eltype, data, kdims, vdims, datatype)
    251     continue
    252 try:
--> 253     (data, dims, extra_kws) = interface.init(eltype, data, kdims, vdims)
    254     break
    255 except DataError:

File ~/micromamba/envs/image_analysis/lib/python3.11/site-packages/holoviews/core/data/grid.py:134, in GridInterface.init(cls, eltype, data, kdims, vdims)
    129         raise error('The shape of the %s value array does not '
    130                     'match the expected dimensionality indicated '
    131                     'by the key dimensions. Expected %d-D array, '
    132                     'found %d-D array.' % (vdim, len(expected), len(shape)))
    133     elif any((e not in (s, s + 1)) for s, e in zip(shape, valid_shape)):
--> 134         raise error(f'Key dimension values and value array {vdim} '
    135                     f'shapes do not match. Expected shape {valid_shape}, '
    136                     f'actual shape: {shape}', cls)
    137 return data, {'kdims':kdims, 'vdims':vdims}, {}

DataError: Key dimension values and value array level_0 shapes do not match. Expected shape (15, 10), actual shape: (10, 15)

GridInterface expects gridded data, for more information on supported datatypes see http://holoviews.org/user_guide/Gridded_Datasets.html
````

</details>


Discussion:
codecov-commenter: ## [Codecov](https://app.codecov.io/gh/holoviz/holoviews/pull/5961?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) Report
> Merging [#5961](https://app.codecov.io/gh/holoviz/holoviews/pull/5961?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) (0b60108) into [main](https://app.codecov.io/gh/holoviz/holoviews/commit/345dd001ce2186f773b11d5f9baf43ca3af79ab7?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) (345dd00) will **decrease** coverage by `65.44%`.
> The diff coverage is `0.00%`.

```diff
@@             Coverage Diff             @@
##             main    #5961       +/-   ##
===========================================
- Coverage   88.58%   23.14%   -65.44%     
===========================================
  Files         315      315               
  Lines       65523    65519        -4     
===========================================
- Hits        58043    15165    -42878     
- Misses       7480    50354    +42874     
```

| [Flag](https://app.codecov.io/gh/holoviz/holoviews/pull/5961/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) | Coverage Δ | |
|---|---|---|
| [ui-tests](https://app.codecov.io/gh/holoviz/holoviews/pull/5961/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) | `23.14% <0.00%> (ø)` | |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz#carryforward-flags-in-the-pull-request-comment) to find out more.

| [Files](https://app.codecov.io/gh/holoviz/holoviews/pull/5961?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) | Coverage Δ | |
|---|---|---|
| [holoviews/element/raster.py](https://app.codecov.io/gh/holoviz/holoviews/pull/5961?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz#diff-aG9sb3ZpZXdzL2VsZW1lbnQvcmFzdGVyLnB5) | `19.77% <0.00%> (-60.23%)` | :arrow_down: |
| [holoviews/tests/plotting/bokeh/test\_rasterplot.py](https://app.codecov.io/gh/holoviz/holoviews/pull/5961?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz#diff-aG9sb3ZpZXdzL3Rlc3RzL3Bsb3R0aW5nL2Jva2VoL3Rlc3RfcmFzdGVycGxvdC5weQ==) | `9.09% <0.00%> (-89.62%)` | :arrow_down: |

... and [298 files with indirect coverage changes](https://app.codecov.io/gh/holoviz/holoviews/pull/5961/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz)

:mega: We’re building smart automated test selection to slash your CI/CD build times. [Learn more](https://about.codecov.io/iterative-testing/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz)

Hoxbro: Added conversion of different inputs:
``` python
import holoviews as hv
import numpy as np

hv.extension("bokeh")

np.random.seed(42)
data = np.random.random((10, 15, 3))

p0 = hv.ImageStack(data)
p1 = hv.ImageStack([data[:, :, i] for i in range(3)])
p2 = hv.ImageStack({str(i): data[:, :, i] for i in range(3)})
p0 + p1 + p2
```

![image](https://github.com/holoviz/holoviews/assets/19758978/4ab4d7cc-e168-404f-9c4a-0d1628bc85f0)


Hoxbro: I'm not sure about the transposing, needs to think more about it.

With the following code:
``` python
import xarray as xr
import holoviews as hv
import numpy as np

hv.extension("bokeh")
x = np.arange(0, 3)
y = np.arange(5, 8)
a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])

ds = (x, y, a, b, c)
im_tuple = hv.ImageStack(ds).opts(title="Tuple")

ds = {"x": x, "y": y, "a": a, "b": b, "c": c}
im_dict = hv.ImageStack(ds).opts(title="Dict")

ds = xr.Dataset(
    {"a": (["x", "y"], a), "b": (["x", "y"], b), "c": (["x", "y"], c)},
    coords={"x": x, "y": y},
)
im_xr = hv.ImageStack(ds).opts(title="xarray")

im_tuple + im_dict + im_xr
```

I get this on the current branch:
![image](https://github.com/holoviz/holoviews/assets/19758978/a21ca27e-dedf-4fed-98d0-aea58fbc8c73)

And this on main:
![image](https://github.com/holoviz/holoviews/assets/19758978/152428d7-e3ac-4d0d-adac-72505c09f29f)

Hoxbro: This work on the current branch but not on main:
``` python
import xarray as xr
import holoviews as hv
import numpy as np

hv.extension("bokeh")
x = np.arange(0, 3)
y = np.arange(5, 9)
a = np.array([[np.nan, np.nan, 1, np.nan], [np.nan] * 4, [np.nan] * 4])
b = np.array([[np.nan] * 4, [1, 1, np.nan, np.nan], [np.nan] * 4])
c = np.array([[np.nan] * 4, [np.nan] * 4, [1, 1, 1, np.nan]])

ds = {"x": x, "y": y, "a": a, "b": b, "c": c}
im_dict = hv.ImageStack(ds).opts(title="Dict")

ds = xr.Dataset(
    {"a": (["x", "y"], a), "b": (["x", "y"], b), "c": (["x", "y"], c)},
    coords={"x": x, "y": y},
)
im_xr = hv.ImageStack(ds).opts(title="xarray")

im_dict + im_xr
```

![image](https://github.com/holoviz/holoviews/assets/19758978/f1811856-3536-4e94-b97e-f0f81d53adaf)


Patch:
diff --git a/examples/reference/elements/bokeh/ImageStack.ipynb b/examples/reference/elements/bokeh/ImageStack.ipynb
index 903364b3c..7db2084ca 100644
--- a/examples/reference/elements/bokeh/ImageStack.ipynb
+++ b/examples/reference/elements/bokeh/ImageStack.ipynb
@@ -9,7 +9,7 @@
     "  <dt>Title</dt> <dd> ImageStack Element</dd>\n",
     "  <dt>Dependencies</dt> <dd>Bokeh</dd>\n",
     "  <dt>Backends</dt>\n",
-    "    <dd><a href='./ImageStack.ipynb'>Bokeh</a></dd>\n",
+    "    <dd><a href='../bokeh/ImageStack.ipynb'>Bokeh</a></dd>\n",
     "    <dd><a href='../matplotlib/ImageStack.ipynb'>Matplotlib</a></dd>\n",
     "    <dd><a href='../plotly/ImageStack.ipynb'>Plotly</a></dd>\n",
     "</dl>\n",
diff --git a/examples/reference/elements/matplotlib/ImageStack.ipynb b/examples/reference/elements/matplotlib/ImageStack.ipynb
index ab5e3cf97..02004d578 100644
--- a/examples/reference/elements/matplotlib/ImageStack.ipynb
+++ b/examples/reference/elements/matplotlib/ImageStack.ipynb
@@ -9,7 +9,7 @@
     "  <dt>Title</dt> <dd> ImageStack Element</dd>\n",
     "  <dt>Dependencies</dt> <dd>Bokeh</dd>\n",
     "  <dt>Backends</dt>\n",
-    "    <dd><a href='./ImageStack.ipynb'>Bokeh</a></dd>\n",
+    "    <dd><a href='../bokeh/ImageStack.ipynb'>Bokeh</a></dd>\n",
     "    <dd><a href='../matplotlib/ImageStack.ipynb'>Matplotlib</a></dd>\n",
     "    <dd><a href='../plotly/ImageStack.ipynb'>Plotly</a></dd>\n",
     "</dl>\n",
diff --git a/examples/reference/elements/plotly/ImageStack.ipynb b/examples/reference/elements/plotly/ImageStack.ipynb
index 8880c51e5..162ad4cef 100644
--- a/examples/reference/elements/plotly/ImageStack.ipynb
+++ b/examples/reference/elements/plotly/ImageStack.ipynb
@@ -9,7 +9,7 @@
     "  <dt>Title</dt> <dd> ImageStack Element</dd>\n",
     "  <dt>Dependencies</dt> <dd>Bokeh</dd>\n",
     "  <dt>Backends</dt>\n",
-    "    <dd><a href='./ImageStack.ipynb'>Bokeh</a></dd>\n",
+    "    <dd><a href='../bokeh/ImageStack.ipynb'>Bokeh</a></dd>\n",
     "    <dd><a href='../matplotlib/ImageStack.ipynb'>Matplotlib</a></dd>\n",
     "    <dd><a href='../plotly/ImageStack.ipynb'>Plotly</a></dd>\n",
     "</dl>\n",
diff --git a/holoviews/core/util.py b/holoviews/core/util.py
index 5b6ecd4bd..f7c651084 100644
--- a/holoviews/core/util.py
+++ b/holoviews/core/util.py
@@ -1218,8 +1218,10 @@ def unique_array(arr):
         if (isinstance(v, datetime_types) and
             not isinstance(v, cftime_types)):
             v = pd.Timestamp(v).to_datetime64()
+        elif isinstance(getattr(v, "dtype", None), pd.CategoricalDtype):
+            v = v.dtype.categories
         values.append(v)
-    return pd.unique(np.asarray(values))
+    return pd.unique(np.asarray(values).ravel())
 
 
 def match_spec(element, specification):
diff --git a/holoviews/element/raster.py b/holoviews/element/raster.py
index ae7720c10..d8ac23e66 100644
--- a/holoviews/element/raster.py
+++ b/holoviews/element/raster.py
@@ -508,16 +508,26 @@ class Image(Selection2DExpr, Dataset, Raster, SheetCoordinateSystem):
 
 class ImageStack(Image):
     """
-    Supports the same constructor RGB and HSV elements, but without the
-    limit of 3/4 channels (one of more channels).
-
-    If only one channel it should behave the same as an Image.
-
-    Type of data inputs:
-    - 3D ndarray (x,y,level)
-    - A list of 2D ndarrays
-    - A dict of 2D ndarrays (key=level: value=2D ndarray)
-    - xarray with all the whistles
+    ImageStack expands the capabilities of Image to by supporting
+    multiple layers of images.
+
+    As there is many ways to represent multiple layers of images,
+    the following options are supported:
+
+        1) A 3D Numpy array with the shape (y, x, level)
+        2) A list of 2D Numpy arrays with identical shape (y, x)
+        3) A dictionary where the keys will be set as the vdims and the
+            values are 2D Numpy arrays with identical shapes (y, x).
+            If the dictionary's keys matches the kdims of the element,
+            they need to be 1D arrays.
+        4) A tuple containing (x, y, level_0, level_1, ...),
+            where the level is a 2D Numpy array in the shape of (y, x).
+        5) An xarray DataArray or Dataset where its `coords` contain the kdims.
+
+    If no kdims are supplied, x and y are used.
+
+    If no vdims are supplied, and the naming can be inferred like with a dictionary
+    the levels will be named level_0, level_1, etc.
     """
 
     vdims = param.List(doc="""
@@ -530,20 +540,35 @@ class ImageStack(Image):
     _vdim_reductions = {1: Image}
 
     def __init__(self, data, kdims=None, vdims=None, **params):
-        if isinstance(data, np.ndarray) and data.ndim == 3:
-            x = np.arange(data.shape[0])
-            y = np.arange(data.shape[1])
-            data = (x, y, *(data[:, :, n] for n in range(data.shape[2])))
+        _kdims = kdims or self.kdims
+        if isinstance(data, list) and len(data):
+            x = np.arange(data[0].shape[1])
+            y = np.arange(data[0].shape[0])
+            data = (x, y, *data)
+        elif isinstance(data, dict):
+            first = next(v for k, v in data.items() if k not in _kdims)
+            xdim, ydim = map(str, _kdims)
+            if xdim not in data:
+                data[xdim] = np.arange(first.shape[1])
+            if ydim not in data:
+                data[ydim] = np.arange(first.shape[0])
+        elif isinstance(data, np.ndarray) and data.ndim == 3:
+            x = np.arange(data.shape[1])
+            y = np.arange(data.shape[0])
+            arr = (data[:, :, n] for n in range(data.shape[2]))
+            data = (x, y, *arr)
         elif (
             isinstance(data, tuple) and len(data) == 3
             and isinstance(data[2], np.ndarray) and data[2].ndim == 3
         ):
-            data = (data[0], data[1], *(data[2][:, :,n] for n in range(data[2].shape[2])))
+            arr = (data[2][:, :, n] for n in range(data[2].shape[2]))
+            data = (data[0], data[1], *arr)
+
         if vdims is None:
             if isinstance(data, tuple):
                 vdims = [Dimension(f"level_{i}") for i in range(len(data[2:]))]
             elif isinstance(data, dict):
-                vdims = [Dimension(key) for key in data.keys() if key not in self.kdims]
+                vdims = [Dimension(key) for key in data.keys() if key not in _kdims]
         super().__init__(data, kdims=kdims, vdims=vdims, **params)
 
 
diff --git a/pyproject.toml b/pyproject.toml
index 0cb9017ac..81ea95642 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -102,6 +102,6 @@ known-first-party = ["holoviews"]
 combine-as-imports = true
 
 [tool.codespell]
-ignore-words-list = "lod,nd,ndoes,reenabled,spreaded,whn"
+ignore-words-list = "lod,nd,ndoes,reenabled,spreaded,whn,ser"
 skip = "doc/generate_modules.py"
 write-changes = true


Test patch:
diff --git a/holoviews/tests/core/test_utils.py b/holoviews/tests/core/test_utils.py
index 1aa628f27..560f11245 100644
--- a/holoviews/tests/core/test_utils.py
+++ b/holoviews/tests/core/test_utils.py
@@ -28,6 +28,7 @@ from holoviews.core.util import (
     sanitize_identifier_fn,
     search_indices,
     tree_attribute,
+    unique_array,
     wrap_tuple_streams,
 )
 from holoviews.element.comparison import ComparisonTestCase
@@ -799,3 +800,9 @@ def test_seach_indices_dtype_object():
     values = np.array(["c0", "c0", np.nan], dtype=object)
     source = np.array(["c0", np.nan], dtype=object)
     search_indices(values, source)
+
+
+def test_unique_array_categorial():
+    ser = pd.Series(np.random.choice(["a", "b", "c"], 100)).astype("category")
+    res = unique_array([ser])
+    assert sorted(res) == ["a", "b", "c"]
diff --git a/holoviews/tests/plotting/bokeh/test_rasterplot.py b/holoviews/tests/plotting/bokeh/test_rasterplot.py
index 7e9c2c8be..c12f96912 100644
--- a/holoviews/tests/plotting/bokeh/test_rasterplot.py
+++ b/holoviews/tests/plotting/bokeh/test_rasterplot.py
@@ -129,83 +129,110 @@ class TestRasterPlot(TestBokehPlot):
         assert cdata["dw"] == [1.0]
         assert cdata["y"] == [-0.5]
 
-    def test_image_stack_tuple(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
 
+class _ImageStackBase(TestRasterPlot):
+    __test__ = False
+
+    def test_image_stack_tuple(self):
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
         img_stack = ImageStack((x, y, a, b, c), kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_tuple_unspecified_dims(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         img_stack = ImageStack((x, y, a, b, c), kdims=["x", "y"])
         assert img_stack.vdims == ["level_0", "level_1", "level_2"]
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_dict(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         ds = {"x": x, "y": y, "a": a, "b": b, "c": c}
         img_stack = ImageStack(ds, kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_dict_unspecified_dims(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         ds = {"x": x, "y": y, "a": a, "b": b, "c": c}
         img_stack = ImageStack(ds)
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
+        assert isinstance(plot, ImageStackPlot)
+
+    def test_image_stack_dict_no_kdims(self):
+        a, b, c = self.a, self.b, self.c
+
+        ds = {"a": a, "b": b, "c": c}
+        img_stack = ImageStack(ds)
+        plot = bokeh_renderer.get_plot(img_stack)
+        source = plot.handles["source"]
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
+        assert source.data["x"][0] == -0.5
+        assert source.data["y"][0] == -0.5
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
+        assert isinstance(plot, ImageStackPlot)
+
+    def test_image_stack_list(self):
+        a, b, c = self.a, self.b, self.c
+
+        ds = [a, b, c]
+        img_stack = ImageStack(ds)
+        plot = bokeh_renderer.get_plot(img_stack)
+        source = plot.handles["source"]
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
+        assert source.data["x"][0] == -0.5
+        assert source.data["y"][0] == -0.5
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_xarray_dataset(self):
@@ -214,27 +241,24 @@ class TestRasterPlot(TestBokehPlot):
         except ImportError:
             raise SkipTest("xarray not available for core tests")
 
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         ds = xr.Dataset(
-            {"a": (["x", "y"], a), "b": (["x", "y"], b), "c": (["x", "y"], c)},
+            {"a": (["y", "x"], a), "b": (["y", "x"], b), "c": (["y", "x"], c)},
             coords={"x": x, "y": y},
         )
-        at, bt, ct = np.dstack([a, b, c]).reshape(3, 9).T.reshape(3, 3, 3)
         img_stack = ImageStack(ds, kdims=["x", "y"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], at)
-        np.testing.assert_equal(source.data["image"][0][1], bt)
-        np.testing.assert_equal(source.data["image"][0][2], ct)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_xarray_dataarray(self):
@@ -243,121 +267,145 @@ class TestRasterPlot(TestBokehPlot):
         except ImportError:
             raise SkipTest("xarray not available for core tests")
 
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
-        ds = xr.Dataset(
-            {"a": (["x", "y"], a), "b": (["x", "y"], b), "c": (["x", "y"], c)},
+        da = xr.Dataset(
+            {"a": (["y", "x"], a), "b": (["y", "x"], b), "c": (["y", "x"], c)},
             coords={"x": x, "y": y},
         ).to_array("level")
-        at, bt, ct = np.dstack([a, b, c]).reshape(3, 9).T.reshape(3, 3, 3)
-        img_stack = ImageStack(ds, vdims=["level"])
+        img_stack = ImageStack(da, vdims=["level"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], at)
-        np.testing.assert_equal(source.data["image"][0][1], bt)
-        np.testing.assert_equal(source.data["image"][0][2], ct)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_invert_xaxis(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         img_stack = ImageStack((x, y, a, b, c), kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack.opts(invert_xaxis=True))
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
 
     def test_image_stack_invert_yaxis(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         img_stack = ImageStack((x, y, a, b, c), kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack.opts(invert_yaxis=True))
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
 
     def test_image_stack_invert_axes(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         img_stack = ImageStack((x, y, a, b, c), kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack.opts(invert_axes=True))
         source = plot.handles["source"]
 
-        at, bt, ct = np.dstack([a, b, c]).reshape(3, 9).T.reshape(3, 3, 3)
-        np.testing.assert_equal(source.data["image"][0][0], at)
-        np.testing.assert_equal(source.data["image"][0][1], bt)
-        np.testing.assert_equal(source.data["image"][0][2], ct)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0].T, a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1].T, b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2].T, c)
         assert source.data["x"][0] == 4.5
         assert source.data["y"][0] == -0.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.ysize
+        assert source.data["dh"][0] == self.xsize
 
     def test_image_stack_array(self):
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        a, b, c = self.a, self.b, self.c
 
         data = np.dstack((a, b, c))
 
         img_stack = ImageStack(data, kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == -0.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
 
     def test_image_stack_tuple_single_3darray(self):
-        x = np.arange(0, 3)
-        y = np.arange(5, 8)
-        a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
-        b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
-        c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        x = np.arange(self.xsize)
+        y = np.arange(self.ysize) + 5
+        a, b, c = self.a, self.b, self.c
 
         data = (x, y, np.dstack((a, b, c)))
 
         img_stack = ImageStack(data, kdims=["x", "y"], vdims=["a", "b", "c"])
         plot = bokeh_renderer.get_plot(img_stack)
         source = plot.handles["source"]
-        np.testing.assert_equal(source.data["image"][0][0], a.T)
-        np.testing.assert_equal(source.data["image"][0][1], b.T)
-        np.testing.assert_equal(source.data["image"][0][2], c.T)
+        np.testing.assert_equal(source.data["image"][0][:, :, 0], a)
+        np.testing.assert_equal(source.data["image"][0][:, :, 1], b)
+        np.testing.assert_equal(source.data["image"][0][:, :, 2], c)
         assert source.data["x"][0] == -0.5
         assert source.data["y"][0] == 4.5
-        assert source.data["dw"][0] == 3
-        assert source.data["dh"][0] == 3
+        assert source.data["dw"][0] == self.xsize
+        assert source.data["dh"][0] == self.ysize
         assert isinstance(plot, ImageStackPlot)
+
+
+class TestImageStackEven(_ImageStackBase):
+    __test__ = True
+
+    def setUp(self):
+        self.a = np.array([[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3])
+        self.b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3])
+        self.c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1]])
+        self.ysize = 3
+        self.xsize = 3
+        super().setUp()
+
+
+class TestImageStackUneven1(_ImageStackBase):
+    __test__ = True
+
+    def setUp(self):
+        self.a = np.array(
+            [[np.nan, np.nan, 1], [np.nan] * 3, [np.nan] * 3, [np.nan] * 3]
+        )
+        self.b = np.array([[np.nan] * 3, [1, 1, np.nan], [np.nan] * 3, [np.nan] * 3])
+        self.c = np.array([[np.nan] * 3, [np.nan] * 3, [1, 1, 1], [np.nan] * 3])
+        self.ysize = 4
+        self.xsize = 3
+        super().setUp()
+
+
+class TestImageStackUneven2(_ImageStackBase):
+    __test__ = True
+
+    def setUp(self):
+        self.a = np.array([[np.nan, np.nan, 1, np.nan], [np.nan] * 4, [np.nan] * 4])
+        self.b = np.array([[np.nan] * 4, [1, 1, np.nan, np.nan], [np.nan] * 4])
+        self.c = np.array([[np.nan] * 4, [np.nan] * 4, [1, 1, 1, 1]])
+        self.ysize = 3
+        self.xsize = 4
+        super().setUp()