diff --git a/NAMESPACE b/NAMESPACE
index 459ac95..5536319 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -48,5 +48,6 @@ importFrom(sf,st_coordinates)
 importFrom(sf,st_geometry_type)
 importFrom(spatialdataR,channels)
 importFrom(spatialdataR,data_type)
+importFrom(spatialdataR,meta)
 importFrom(spatialdataR,transform)
 importFrom(utils,tail)
diff --git a/R/plotImage.R b/R/plotImage.R
index ac115bb..86cbfa6 100644
--- a/R/plotImage.R
+++ b/R/plotImage.R
@@ -161,21 +161,29 @@ NULL
 #' @importFrom spatialdataR data_type
 .df_i <- \(x, k=NULL, ch=NULL, c=NULL, cl=NULL) {
     a <- .get_multiscale_data(x, k)
+    # max-projection over z-stacks
+    d <- length(dim(x))
+    if (d == 4) a <- apply(a, c(1, 3, 4), max)
+    # subset channels of interest
     a <- a[.ch_idx(x, ch),,,drop=FALSE]
     a <- .norm_ia(a, data_type(x))
+    # color merging & contrasts
     a <- .prep_ia(a, c, cl)
 }
 
-#' @importFrom spatialdataR transform
+#' @importFrom utils tail
 .get_wh <- \(x) {
     wh <- metadata(x)$wh
     if (!is.null(wh)) {
         df <- data.frame(x=wh[[1]], y=wh[[2]])
     } else {
         ds <- dim(data(x, 1))
-        df <- data.frame(x=c(0, ds[3]), y=c(0, ds[2]))
+        df <- data.frame(
+            x=c(0, tail(ds, 1)), 
+            y=c(0, tail(ds, 2)[1]))
     }
-    list(w=df[, 1], h=df[, 2])
+    wh <- list(w=df$x, h=df$y)
+    return(wh)
 }
 
 #' @importFrom ggplot2 guides geom_point geom_blank annotation_raster 
@@ -200,7 +208,6 @@ setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=NULL, cl
     if (is.numeric(j))
         j <- CTname(y)[j]
     y <- transform(y, j)
-    wh <- .get_wh(y)
     if (.is_rgb(y)) {
         # RGB: we plot everything by default and we don't normalize
         ch <- ch %||% channels(y)
@@ -212,6 +219,16 @@ setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=NULL, cl
         nms <- unlist(channels(y))[idx <- .ch_idx(y, ch)]
         pal <- pal[seq_along(idx)]; names(pal) <- nms
     }
+    # multi-scale adjustment
+    wh <- .get_wh(y)
+    if (wh$w[2] == tail(dim(y), 1) ||
+        wh$h[2] == tail(dim(y), 2)[1]) {
+        ts <- .get_multiscale_scale(y)
+        tx <- tail(ts, 1)
+        ty <- tail(ts, 2)[1]
+    } else tx <- ty <- 1
+    wh$w[2] <- wh$w[2]*tx
+    wh$h[2] <- wh$h[2]*ty
     .gg_i(df, wh$w, wh$h, pal)
 })
 
diff --git a/R/plotLabel.R b/R/plotLabel.R
index 534c7dc..f337def 100644
--- a/R/plotLabel.R
+++ b/R/plotLabel.R
@@ -16,6 +16,9 @@
 #' @param pal character vector; color for discrete/continuous values
 #'   (interpolated automatically when insufficient values are provided).
 #' @param nan character string; color for missing values (hidden by default).
+#' @param z scalar integer; 
+#'   specifies which z-slice to plot when \code{label(x, i)} is 3D; 
+#'   by default (NULL), will apply a max-projection across all z-slices.
 #' 
 #' @examples
 #' x <- system.file("extdata", "blobs.zarr", package="spatialdataR")
@@ -54,21 +57,52 @@ NULL
 #' @importFrom SingleCellExperiment colData
 #' @export
 setMethod("plotLabel", "SpatialData", \(x, i=1, j=1, k=NULL, c=NULL, 
-    a=0.5, pal=c("red", "green"), nan=NA, assay=1) {
+    a=0.5, pal=c("red", "green"), nan=NA, assay=1, z=NULL) {
     if (is.numeric(i)) i <- labelNames(x)[i]
     i <- match.arg(i, labelNames(x))
     y <- label(x, i)
+    
     # transformation
     if (is.numeric(j))
       j <- CTname(y)[j]
     y <- transform(y, j)
+
+    # get array data
     ym <- .get_multiscale_data(y, k)
-    wh <- .get_wh(y)
+    if (length(dim(ym)) > 2) {
+        if (is.null(z)) {
+            # max-projection across z-slices
+            nm <- vapply(axes(y), \(.) .$name, character(1))
+            yx <- match(c("y", "x"), nm)
+            ym <- apply(ym, yx, max)
+        } else {
+            # subset target z-slice
+            ym <- ym[z,,]
+        }
+    }
   
-    # Keep only indices != 0 since labels might be sparse and thus save memory by not plotting all pixels
+    # keep only indices != 0 since labels might be sparse 
+    # and thus save memory by not plotting all pixels
     idx <- BiocGenerics::which(ym != 0L, arr.ind=TRUE)
-    # All other SD elements are flipped when plotted. Let's keep the same convention here.
-    df <- data.frame(x=idx[,2L]+wh$w[1], y=idx[,1L]+wh$h[1], z=ym[idx])
+    
+    # offset & multi-scale adjustment
+    ds <- dim(ym)
+    wh <- .get_wh(y)
+    if (wh$w[2] == tail(dim(y), 1) ||
+        wh$h[2] == tail(dim(y), 2)[1]) {
+        ts <- .get_multiscale_scale(y)
+        tx <- tail(ts, 1)
+        ty <- tail(ts, 2)[1]
+    } else tx <- ty <- 1
+    nx <- tail(ds, 1)
+    ny <- tail(ds, 2)[1]
+    sx <- (diff(wh$w)/nx)*tx
+    sy <- (diff(wh$h)/ny)*ty
+    df <- data.frame(
+        x=wh$w[1]+idx[,2L]*sx, 
+        y=wh$h[1]+idx[,1L]*sy, 
+        z=ym[idx])
+    
     aes <- aes(.data[["x"]], .data[["y"]])
     if (!is.null(c)) {
         stopifnot(length(c) == 1, is.character(c))
@@ -77,7 +111,7 @@ setMethod("plotLabel", "SpatialData", \(x, i=1, j=1, k=NULL, c=NULL,
         # TODO: search ik in both internal and regular colData for now
         # thus perhaps update, spatialdataR::valTable instead
         # idx <- match(df$z, int_colData(t)[[ik]])
-        if(ik %in% names(int_colData(t))){
+        if (ik %in% names(int_colData(t))){
           coldata <- int_colData(t)[[ik]]
         } else {
           coldata <- colData(t)[[ik]]
@@ -105,4 +139,4 @@ setMethod("plotLabel", "SpatialData", \(x, i=1, j=1, k=NULL, c=NULL,
             scale_fill_manual(NULL, values=pal))
     }
     list(thm, do.call(geom_tile, list(data=df, mapping=aes, alpha=a)))
-})
\ No newline at end of file
+})
diff --git a/R/utils.R b/R/utils.R
index e4b7f0c..216d918 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -47,9 +47,8 @@
 
 # guess scale of image or label
 .guess_scale <- \(x, w, h) {
-  n <- length(dim(x))
-  i <- ifelse(n == 3, -1, TRUE)
-  d <- vapply(x@data, dim, numeric(n))
+  i <- match(c("y", "x"), vapply(axes(x), \(.) .$name, character(1)))
+  d <- vapply(x@data, dim, numeric(length(dim(x))))
   d <- apply(d, 2, \(.) sum(abs(.[i]-c(h, w))))
   which.min(d)
 }
@@ -58,4 +57,12 @@
 .get_multiscale_data <- \(x, k=NULL, w=800, h=800) {
   if (!is.null(k)) return(data(x, k))
   data(x, .guess_scale(x, w, h))
-}
\ No newline at end of file
+}
+
+#' @importFrom spatialdataR meta
+.get_multiscale_scale <- \(x) {
+    ms <- spatialdataR:::multiscales(meta(x))[[1]]
+    ds <- ms$datasets[[1]]
+    ct <- ds$coordinateTransformations[[1]]
+    return(unlist(ct$scale))
+}
diff --git a/man/plotLabel.Rd b/man/plotLabel.Rd
index cef3af7..eb54847 100644
--- a/man/plotLabel.Rd
+++ b/man/plotLabel.Rd
@@ -14,7 +14,8 @@
   a = 0.5,
   pal = c("red", "green"),
   nan = NA,
-  assay = 1
+  assay = 1,
+  z = NULL
 )
 }
 \arguments{
@@ -41,6 +42,10 @@ a \code{colData} column or row name in a \code{table} annotating \code{i}.}
 
 \item{assay}{character string; in case of \code{c} denoting a row name,
 specifies which \code{assay} data to use (see \code{\link{valTable}}).}
+
+\item{z}{scalar integer; 
+specifies which z-slice to plot when \code{label(x, i)} is 3D; 
+by default (NULL), will apply a max-projection across all z-slices.}
 }
 \description{
 \code{SpatialData} label viz.