diff --git a/NAMESPACE b/NAMESPACE
index cf5f75db..2820a59d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -194,7 +194,6 @@ importFrom(methods,"slot<-")
 importFrom(methods,as)
 importFrom(methods,callNextMethod)
 importFrom(methods,is)
-importFrom(methods,new)
 importFrom(methods,setClass)
 importFrom(methods,setClassUnion)
 importFrom(methods,setMethod)
diff --git a/R/CTgraph.R b/R/CTgraph.R
index 8a2b6104..990cb168 100644
--- a/R/CTgraph.R
+++ b/R/CTgraph.R
@@ -76,13 +76,11 @@ setMethod("CTgraph", "ANY", \(x) stop("'x' should be a",
 .make_g <- \(md) {
     g <- .init_g()
     for (l in names(md)) for (e in names(md[[l]])) {
-        .md <- md[[l]][[e]]
-        ms <- multiscales(.md)
-        if (!is.null(ms)) .md <- ms[[1]]
         .e <- paste0("_", e)
         g <- addNode(.e, g)
         nodeData(g, .e, "type") <- "element"
-        ct <- .md$coordinateTransformations
+        ms <- .get_ms(md[[l]][[e]])
+        ct <- ms$coordinateTransformations
         for (i in seq_along(ct)) {
             n <- ct[[i]]$output$name
             if (!n %in% nodes(g)) {
diff --git a/R/CTutils.R b/R/CTutils.R
index 368d17ce..0dbe1155 100644
--- a/R/CTutils.R
+++ b/R/CTutils.R
@@ -55,9 +55,8 @@ NULL
 #' @rdname CTutils
 #' @export
 setMethod("axes", "SpatialDataAttrs", \(x, y=NULL, ...) {
-    ms <- multiscales(x)
-    if (!is.null(ms)) x <- ms[[1]]
-    if (is.null(x <- x$axes)) stop("couldn't find 'axes'") 
+    x <- .get_ms(x)$axes
+    if (is.null(x)) stop("couldn't find 'axes'") 
     if (is.null(y)) return(x)
     y <- match.arg(y, c("name", "type", "unit"))
     vapply(x, `[[`, character(1), y)
@@ -68,24 +67,14 @@ setMethod("axes", "SpatialDataAttrs", \(x, y=NULL, ...) {
 #' @rdname CTutils
 #' @export
 setMethod("CTlist", "SpatialDataAttrs", \(x, ...) {
-    ms <- multiscales(x)
     ct <- "coordinateTransformations"
-    if (is.null(ms)) return(x[[ct]])
-    ms[[1]][[ct]]
+    .get_ms(x)[[ct]]
 })
 
 #' @rdname CTutils
 #' @export
 setMethod("CTdata", "SpatialDataAttrs", \(x, i=1, ...) {
-    stopifnot(length(i) == 1)
-    if (is.character(i)) {
-        match.arg(i, CTname(x))
-        i <- match(i, CTname(x))
-    } else if (is.numeric(i)) {
-        stopifnot(
-            i == round(i), 
-            i %in% seq_along(CTlist(x)))
-    } else stop("Invalid 'i'; should be a scalar character or integer")
+    i <- .resolve_id(i, CTname(x))
     t <- CTtype(x)[i]
     if (t != "sequence") 
         return(CTlist(x)[[i]][[t]])
diff --git a/R/centroids.R b/R/centroids.R
index b9fb3d28..a1dff0e9 100644
--- a/R/centroids.R
+++ b/R/centroids.R
@@ -34,13 +34,12 @@ setMethod("centroids", "ANY", \(x, ...) stop("'centroids' ",
 #' @importFrom Matrix summary
 setMethod("centroids", "SpatialDataLabel", \(x,
     as=c("data.frame", "matrix")) {
-    as <- match.arg(as)
     y <- data(x)
-    if (length(dim(y)) > 2) {
-        # max-projection
-        ax <- match(c("y", "x"), axes(x, "name"))
-        y <- apply(y, ax, max)
-    }
+    as <- match.arg(as)
+    ax <- .get_xy_axes(x)
+    # max-projection
+    if (length(dim(y)) > 2) 
+        y <- apply(y, c(ax$y, ax$x), max)
     y <- as(y, "dgCMatrix")
     i <- summary(y)
     # flip dimensions so that columns=x, rows=y
@@ -50,9 +49,9 @@ setMethod("centroids", "SpatialDataLabel", \(x,
     xy <- cbind(xy, as.integer(rownames(xy)))
     dimnames(xy) <- list(NULL, c("x", "y", "i"))
     # multi-scale adjustment
-    sf <- .get_multiscale_scale(x)
-    xy[,1] <- xy[,1]*tail(sf, 1)
-    xy[,2] <- xy[,2]*tail(sf, 2)[1]
+    sf <- .get_ms_scale(x)
+    xy[,1] <- xy[,1]*sf[ax$x]
+    xy[,2] <- xy[,2]*sf[ax$y]
     # offset
     wh <- metadata(x)$wh
     if (!is.null(wh)) {
diff --git a/R/crop.R b/R/crop.R
index 0eeaabce..be61d8e6 100644
--- a/R/crop.R
+++ b/R/crop.R
@@ -106,10 +106,10 @@ NULL
         x=c(y$xmin, y$xmax, y$xmax, y$xmin, y$xmin),
         y=c(y$ymin, y$ymin, y$ymax, y$ymax, y$ymin),
         id=seq_len(5))
-    # get transformation for space j
-    if (is.numeric(j)) j <- CTname(x)[j]
-    ct <- CTlist(x)[[match(j, CTname(x))]]
-    # identify spatial axes
+    # get transformation for space 'j'
+    j <- .resolve_id(j, CTname(x))
+    ct <- CTlist(x)[[j]]
+    # helper to adapt transformation data to spatial (XY) dims
     axs <- axes(x)
     nms <- vapply(axs, \(.) .$name, character(1))
     ix <- match("x", nms)
@@ -119,11 +119,11 @@ NULL
         n <- length(nms)
         ix <- n; iy <- n-1
     }
-    # helper to adapt transformation data to spatial (XY) dims
+    ax <- .get_xy_axes(x)
     .adapt <- \(t, type) {
         if (is.null(t)) return(NULL)
         if (type %in% c("scale", "translation"))
-            return(c(t[ix], t[iy]))
+            return(c(t[ax$x], t[ax$y]))
         if (type == "rotate") 
             return(t[1])
         return(t)
@@ -195,7 +195,7 @@ setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
     }
     metadata(x)$wh <- wh
     # multi-scale adjustment
-    t <- .get_multiscale_scale(x)
+    t <- .get_ms_scale(x)
     tx <- tail(t, 1)
     ty <- tail(t, 2)[1]
     z$xmin <- floor(z$xmin/tx)
@@ -209,13 +209,6 @@ setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
     if (ii) x[, i, j] else x[i, j] 
 })
 
-.get_multiscale_scale <- \(x) {
-    ms <- multiscales(meta(x))[[1]]
-    ds <- ms$datasets[[1]]
-    ct <- ds$coordinateTransformations[[1]]
-    return(unlist(ct$scale))
-}
-
 #' @export
 #' @rdname crop
 #' @importFrom dplyr pull .data
diff --git a/R/mask.R b/R/mask.R
index 10dde115..4bed9970 100644
--- a/R/mask.R
+++ b/R/mask.R
@@ -56,16 +56,7 @@ setMethod("mask", c("SpatialData", "ANY", "ANY"), \(x, i, j, k,
     if (!length(ct)) stop(
         "can't mask; found no common ",
         "coordinates between 'i' and 'j'")
-    if (missing(k)) {
-        k <- 1
-    } else {
-        if (is.character(k)) {
-            k <- match.arg(k, ct)
-            k <- match(k, ct)
-        } else if (is.numeric(k)) {
-            stopifnot(k > 0, k <= length(ct))
-        }
-    }
+    k <- if (missing(k)) 1 else .resolve_id(k, ct)
     .i <- transform(.i, ct[k])
     .j <- transform(.j, ct[k])
     t <- tryCatch(error=\(.) NULL, getTable(x, i))
@@ -98,18 +89,7 @@ setMethod("mask_i_by_j",
         how <- "mean"
     }
     # default to 1st matching scale
-    di <- lapply(data(i, NULL), dim)
-    dj <- lapply(data(j, NULL), dim)
-    ai <- axes(i, "type") == "space"
-    aj <- axes(j, "type") == "space"
-    ks <- outer(
-        seq_along(di),
-        seq_along(dj),
-        Vectorize(\(i, j) identical(di[[i]][ai], dj[[j]][aj])))
-    ks <- which(ks, arr.ind=TRUE)
-    if (nrow(ks) == 0)
-        stop("couldn't find shared multiscales level between label/image;",
-            " need at least one data() pair with identical dimensions")
+    ks <- .get_ms_match(i, j)
     di <- data(i, ks[1, 1])
     dj <- data(j, ks[1, 2])
     # utility to aggregate 'i' channels by instance in 'j'
diff --git a/R/read.R b/R/read.R
index 5e9fea41..78874855 100644
--- a/R/read.R
+++ b/R/read.R
@@ -136,15 +136,13 @@ readSpatialData <- function(x,
     
     # helper for layer reading
     .readLayer <- \(l) {
+        # 'j' are the paths on disk, 'nms' are their basenames
         j <- list.dirs(file.path(x, l), recursive=FALSE, full.names=TRUE)
-        names(j) <- basename(j)
+        nms <- names(j) <- basename(j)
         opt <- args[[l]]
         if (!isTRUE(opt)) {
-            if (is.numeric(opt) && opt > (. <- length(j)))
-                stop("'", l, "=", opt, "', but only ", ., " elements found")
-            if (is.character(opt) && length(. <- setdiff(opt, basename(j))))
-                stop("couldn't find ", l, " of name", .)
-            j <- j[opt]
+            # validate each requested element
+            j <- j[vapply(opt, .resolve_id, integer(1), ok=nms, nm=l)]
         }
         f <- paste0("read", toupper(substr(l, 1, 1)), substr(l, 2, nchar(l)-1))
         lapply(j, f)
diff --git a/R/sdArray.R b/R/sdArray.R
index 82c73701..937546ae 100644
--- a/R/sdArray.R
+++ b/R/sdArray.R
@@ -63,21 +63,10 @@ NULL
 
 # new ----
 
-.new_sda <- \(type, data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
-    if (is.array(data)) data <- list(data)
-    x <- new(type, data=data, meta=meta, ...)
-    metadata(x) <- metadata
-    return(x)
-}
-
-SpatialDataImage <- \(...) .new_sda("SpatialDataImage", ...)
-SpatialDataLabel <- \(...) .new_sda("SpatialDataLabel", ...)
-
 #' @export
 #' @rdname SpatialDataArray
-#' @importFrom methods new
 #' @importFrom S4Vectors metadata<-
-SpatialDataImage <- \(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
+SpatialDataImage <- \(data=list(), meta=SpatialDataAttrs(type="image"), metadata=list(), ...) {
     if (is.array(data)) data <- list(data)
     x <- .SpatialDataImage(data=data, meta=meta, ...)
     metadata(x) <- metadata
@@ -86,9 +75,8 @@ SpatialDataImage <- \(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...
 
 #' @export
 #' @rdname SpatialDataArray
-#' @importFrom methods new
 #' @importFrom S4Vectors metadata<-
-SpatialDataLabel <- \(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
+SpatialDataLabel <- \(data=list(), meta=SpatialDataAttrs(type="label"), metadata=list(), ...) {
     if (is.array(data)) data <- list(data)
     x <- .SpatialDataLabel(data=data, meta=meta, ...)
     metadata(x) <- metadata
@@ -178,12 +166,14 @@ setMethod("channels", "SpatialDataElement", \(x, ...) stop("only 'images' have c
 
 #' @importFrom utils head tail
 .sub_sda <- \(x, yx, z=list()) {
-    #x <- label(sd); yx <- list(1:10, 1:10); z <- list()
-    # yx: spatial; z: channels
+    # yx: user-provided spatial slices (list of 2: y, x)
+    # z: user-provided channel slices (list of 1)
+    ax <- .get_xy_axes(x)
     ls <- seq_along(data(x, NULL))
     data(x) <- lapply(ls, \(l) {
-        sf <- 2^(l-1)   
-        rc <- tail(dim(data(x, l)), 2) 
+        sf <- 2^(l-1) # scale factor for current level
+        ds <- dim(data(x, l))
+        rc <- ds[c(ax$y, ax$x)] # numbers rows/cols (YX)
         # get spatial indices
         .yx <- lapply(seq_along(yx), \(a) {
             ix <- yx[[a]]
diff --git a/R/sdAttrs.R b/R/sdAttrs.R
index 143e1f75..e70c4caf 100644
--- a/R/sdAttrs.R
+++ b/R/sdAttrs.R
@@ -15,7 +15,10 @@
 #' @param trans list of coordinate transformations; defaults to identity only.
 #' @param value character string (for one \code{region} and \code{_key}s), 
 #'   or vector (for many \code{region}s, \code{instances} and \code{regions}).
-#' @param ver character string; specified the .zarr version to comply with.
+#' @param ver character string; specifies the OME version to comply with.
+#' @param dim scalar integer in 2-4;
+#'   number of dimensions: 2 = XY, 3 adds Z, 4 adds T (time); 
+#'   when \code{type="image"}, C (channel) will be added (for any \code{dim}).
 #' @param nch scalar integer; how many channels should there be?
 #'   (ignored unless \code{type="frame"} and \code{label=FALSE}). 
 #' @param ... additional attributes (e.g., version, feature_key).
@@ -53,31 +56,24 @@
 #' 
 #' # constructor
 #' SpatialDataAttrs(type="frame")
-#' SpatialDataAttrs(type="array")
-#' SpatialDataAttrs(type="array", nch=7)
-#' SpatialDataAttrs(type="array", label=TRUE)
+#' SpatialDataAttrs(type="image", nch=7)
+#' SpatialDataAttrs(type="label", dim=3)
 #' 
 #' @export
-SpatialDataAttrs <- \(x, type=c("array", "frame"), 
-    label=FALSE, trans=NULL, ver="0.4", nch=3, ...) 
+SpatialDataAttrs <- \(x, type=c("image", "label", "frame"), 
+    trans=NULL, ver="0.4", dim=2, nch=3, ...) 
 {
+    stopifnot(
+        length(dim) == 1, is.numeric(dim), dim %in% seq(2, 4),
+        length(nch) == 1, is.numeric(nch), round(nch) == nch, nch > 0)
     if (!missing(x)) return(.SpatialDataAttrs(x))
     type <- match.arg(type)
-    # axes:
-    # xy for points/shapes
-    ax <-  list(
-        list(name="x", type="space"), 
-        list(name="y", type="space"))
-    if (type == "array") {
-        # yx for labels
-        ax <- rev(ax)
-        # cyx for images
-        if (!label) ax <- c(list(list(name="c", type="channel")), ax)
-    }
+    ver <- .val_ome_ver(ver)
+    ax <- .default_ax(type, dim)
     # transformations:
     ct <- trans %||% .default_ct(ax)
     # .zattrs list:
-    if (type == "array") {
+    if (type != "frame") {
         # default structure
         res <- list(
             omero=list(channels=list(label=letters[seq_len(nch)])),
@@ -96,17 +92,32 @@ SpatialDataAttrs <- \(x, type=c("array", "frame"),
 }
 
 # Internal helper to generate OME-NGFF axes
-.default_ax <- \(type=c("array", "frame")) {
-    switch(match.arg(type),
-        # cyx for images/labels
-        array=list(
-            list(name="c", type="channel"),
-            list(name="y", type="space"),
-            list(name="x", type="space")),
-        # xy for points/shapes
-        list(
-            list(name="x", type="space"),
-            list(name="y", type="space")))
+.default_ax <- \(type=c("image", "label", "frame"), dim=2) {
+    c <- list(name="c", type="channel")
+    t <- list(name="t", type="time")
+    z <- list(name="z", type="space")
+    y <- list(name="y", type="space")
+    x <- list(name="x", type="space")
+    switch(match.arg(type), 
+        # xyzt for points/shapes
+        frame={
+            ax <- list(x, y)
+            if (dim > 2) {
+                ax <- c(ax, list(z))
+                if (dim > 3) ax <- c(ax, list(t))
+            }
+        },
+        # tczyx for images/labels
+        {
+            ax <- list(y, x)
+            if (dim > 2) {
+                ax <- c(list(z), ax)
+                if (dim > 3) ax <- c(list(t), ax)
+            }
+            if (type == "image") ax <- c(list(c), ax)
+        }
+    )
+    return(ax)
 }
 
 # Internal helper to generate coordinate transformations
@@ -132,8 +143,7 @@ setMethod("$", "SpatialDataAttrs", \(x, name) x[[name]])
         x$omero$version %||% 
         x$ome$version
     if (!length(v)) stop("couldn't find 'version' in 'spatialdata_attrs'")
-    ok <- length(v) == 1 && is.character(v) && (v <- gsub("-.*", "", v)) %in% sprintf("0.%d", seq_len(6))
-    if (!ok) stop("invalid OME 'version'; expected '0.x' where x is 1-6")
+    v <- .val_ome_ver(v)
     return(v)
 }
 
@@ -148,9 +158,8 @@ setMethod("multiscales", "list", \(x) {
 # internal use only!
 #' @noRd 
 setMethod("datasets", "list", \(x, ...) {
-  vapply(multiscales(x)[[1]]$datasets, \(.){
-    .$path
-  }, character(1))
+    ds <- .get_ms(x)$datasets
+    vapply(ds, \(.) .$path, character(1))
 })
 
 # features ----
diff --git a/R/trans.R b/R/trans.R
index e4109511..4dd504a3 100644
--- a/R/trans.R
+++ b/R/trans.R
@@ -52,13 +52,7 @@ NULL
 #' @rdname trans
 #' @importFrom BiocGenerics transform
 setMethod("transform", "SpatialDataElement", \(x, i=1, ...) {
-    stopifnot(
-        length(i) == 1, is.character(i) | 
-        (is.numeric(i) && i == round(i)))
-    if (is.character(i)) {
-        i <- match.arg(i, CTname(x))
-        i <- match(i, CTname(x))
-    }
+    i <- .resolve_id(i, CTname(x))
     f <- CTtype(x)[i]
     t <- CTdata(x, i)
     if (f == "sequence") {
@@ -129,7 +123,8 @@ setMethod("rotate", "SpatialDataArray", \(x, t, k=1, ..., rev=FALSE) {
 
 .trans_a <- \(x, t, f=c("scale", "translation"), k=1, rev=FALSE) {
     f <- match.arg(f)
-    n <- length(d <- dim(data(x, k)))
+    d <- dim(data(x, k))
+    n <- length(d)
     
     # setup: identity, operator
     map <- list(
@@ -142,13 +137,18 @@ setMethod("rotate", "SpatialDataArray", \(x, t, k=1, ..., rev=FALSE) {
     if (rev) t <- if (f == "scale") 1/t else -t
     
     # project to spatial (XY) dims
-    if (n == 3) { t <- t[-1]; d <- d[-1] }
-    t <- rev(t); d <- rev(d)
+    ax <- .get_xy_axes(x)
+    xy <- c(ax$x, ax$y)
+    .t <- t[xy]
+    .d <- d[xy]
+    # t <- t[c(ax$x, ax$y)]
+    # t <- rev(t); d <- rev(d)
+    # d <- rev(dim(data(x, k)))
     
     # update 'wh' metadata
-    wh <- metadata(x)$wh %||% list(c(0, d[1]), c(0, d[2]))
+    wh <- metadata(x)$wh %||% list(c(0, .d[1]), c(0, .d[2]))
     op <- get(map$ops[f])
-    metadata(x)$wh <- mapply(op, t, wh, SIMPLIFY=FALSE)
+    metadata(x)$wh <- mapply(op, .t, wh, SIMPLIFY=FALSE)
     return(x)
 }
 
diff --git a/R/utils.R b/R/utils.R
index 409f106a..db86ae46 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -10,6 +10,22 @@
     .GlobalEnv[[nm]]
 }
 
+# internal helper to resolve name/index to integer index
+.resolve_id <- \(i, ok, nm=deparse1(substitute(i))) {
+    nm <- sprintf("'%s'", nm)
+    if (is.character(i)) {
+        i <- match.arg(i, ok)
+        return(match(i, ok))
+    }
+    if (is.numeric(i) && i == round(i) && length(i) == 1) {
+        if (i < 1 || i > length(ok)) {
+            stop(sprintf("invalid %s index: %d (max: %d)", nm, i, length(ok)))
+        }
+        return(as.integer(i))
+    }
+    stop(sprintf("invalid %s; expected character or integer index", nm))
+}
+
 # internal helper for null-coalescing
 `%||%` <- \(a, b) if (is.null(a)) b else a
 
@@ -103,3 +119,65 @@
     slot(x, "tables") <- ts
     return(x)
 }
+
+# internal helper to resolve spatial (XY) axis indices
+.get_xy_axes <- \(x) {
+    nm <- axes(x, "name")
+    ix <- match("x", nm)
+    iy <- match("y", nm)
+    # fallback: OME-NGFF usually places spatial dimensions at the end (YX)
+    if (is.na(ix) || is.na(iy)) {
+        n <- length(nm)
+        ix <- n
+        iy <- n-1
+    }
+    return(list(x=ix, y=iy))
+}
+
+# validation ----
+
+# validate OME version
+.val_ome_ver <- \(v) {
+    ok <- length(v) == 1 && is.character(v) && (v <- gsub("-.*", "", v)) %in% sprintf("0.%d", seq_len(6))
+    if (!ok) stop("invalid OME 'version'; expected '0.x' where x is 1-6")
+    return(v)
+}
+
+# multiscales ----
+
+# internal helper to get the 'active' metadata level
+# (drills into 'multiscales' if present, else returns the list itself)
+.get_ms <- \(x) {
+    # if 'x' is an element, get its attributes first
+    if (is(x, "SpatialDataElement")) x <- meta(x)
+    # check for 'multiscales' (handles OME version via 'multiscales()'
+    ms <- multiscales(x)
+    if (is.null(ms)) return(x)
+    # return the first (usually only) multiscale level's metadata;
+    # this contains the 'axes' & 'coordinateTransformations' we need
+    return(ms[[1]])
+}
+
+# get scale factors between 'multiscales' levels
+# (returns numeric vector, one value per dimension)
+.get_ms_scale <- \(x) {
+    ms <- .get_ms(x)
+    ds <- ms$datasets[[1]]
+    ct <- ds$coordinateTransformations[[1]]
+    return(unlist(ct$scale))
+}
+
+# find indices with equal spatial extents
+# (returns array: rows = matches, cols = x/y indices)
+.get_ms_match <- \(x, y) {
+    ax <- axes(x, "type") == "space"
+    ay <- axes(y, "type") == "space"
+    dx <- lapply(data(x, NULL), \(d) dim(d)[ax])
+    dy <- lapply(data(y, NULL), \(d) dim(d)[ay])
+    ks <- outer(dx, dy, Vectorize(identical))
+    ks <- which(ks, arr.ind=TRUE)
+    if (nrow(ks) == 0)
+        stop("couldn't find shared multiscales level; need at",
+            " least one data() pair with identical dimensions")
+    return(ks)
+}
diff --git a/man/SpatialDataArray.Rd b/man/SpatialDataArray.Rd
index 4cbe62a8..8b26c4ce 100644
--- a/man/SpatialDataArray.Rd
+++ b/man/SpatialDataArray.Rd
@@ -19,14 +19,14 @@
 \usage{
 SpatialDataImage(
   data = list(),
-  meta = SpatialDataAttrs(),
+  meta = SpatialDataAttrs(type = "image"),
   metadata = list(),
   ...
 )
 
 SpatialDataLabel(
   data = list(),
-  meta = SpatialDataAttrs(),
+  meta = SpatialDataAttrs(type = "label"),
   metadata = list(),
   ...
 )
diff --git a/man/SpatialDataAttrs.Rd b/man/SpatialDataAttrs.Rd
index 25f77d08..839ba4e2 100644
--- a/man/SpatialDataAttrs.Rd
+++ b/man/SpatialDataAttrs.Rd
@@ -38,10 +38,10 @@
 \usage{
 SpatialDataAttrs(
   x,
-  type = c("array", "frame"),
-  label = FALSE,
+  type = c("image", "label", "frame"),
   trans = NULL,
   ver = "0.4",
+  dim = 2,
   nch = 3,
   ...
 )
@@ -91,11 +91,13 @@ SpatialDataAttrs(
 
 \item{type}{character string; either "array" (image/label) or "frame" (point/shape).}
 
-\item{label}{flag; when \code{type="frame"}, should attributes be for a label?}
-
 \item{trans}{list of coordinate transformations; defaults to identity only.}
 
-\item{ver}{character string; specified the .zarr version to comply with.}
+\item{ver}{character string; specifies the OME version to comply with.}
+
+\item{dim}{scalar integer in 2-4;
+number of dimensions: 2 = XY, 3 adds Z, 4 adds T (time); 
+when \code{type="image"}, C (channel) will be added (for any \code{dim}).}
 
 \item{nch}{scalar integer; how many channels should there be?
 (ignored unless \code{type="frame"} and \code{label=FALSE}).}
@@ -106,6 +108,8 @@ SpatialDataAttrs(
 
 \item{value}{character string (for one \code{region} and \code{_key}s), 
 or vector (for many \code{region}s, \code{instances} and \code{regions}).}
+
+\item{label}{flag; when \code{type="frame"}, should attributes be for a label?}
 }
 \value{
 character string
@@ -144,8 +148,7 @@ CTdata(z, "scale")
 
 # constructor
 SpatialDataAttrs(type="frame")
-SpatialDataAttrs(type="array")
-SpatialDataAttrs(type="array", nch=7)
-SpatialDataAttrs(type="array", label=TRUE)
+SpatialDataAttrs(type="image", nch=7)
+SpatialDataAttrs(type="label", dim=3)
 
 }
diff --git a/tests/testthat/test-ctutils.R b/tests/testthat/test-ctutils.R
index b03ec1be..5075e441 100644
--- a/tests/testthat/test-ctutils.R
+++ b/tests/testthat/test-ctutils.R
@@ -2,6 +2,27 @@ x <- file.path("extdata", "blobs.zarr")
 x <- system.file(x, package="spatialdataR")
 x <- readSpatialData(x)
 
+test_that("axes", {
+    es <- list(image(x), label(x))
+    for (e in es) {
+        z <- axes(e)
+        d <- length(dim(e))
+        expect_is(z, "list")
+        expect_length(z, d)
+        expect_error(axes(e, "bad"))
+        # name
+        expect_silent(z <- axes(e, "name"))
+        expect_is(z, "character")
+        expect_length(z, d)
+        expect_in(z, c("t","c","z","y","x"))
+        # type
+        expect_silent(z <- axes(e, "type"))
+        expect_is(z, "character")
+        expect_length(z, d)
+        expect_in(z, c("time","channel","space"))
+    }
+})
+
 .CTtype <- c(
     "identity", "scale", "rotate", 
     "translation", "affine", "sequence")
diff --git a/tests/testthat/test-misc.R b/tests/testthat/test-misc.R
new file mode 100644
index 00000000..4159d027
--- /dev/null
+++ b/tests/testthat/test-misc.R
@@ -0,0 +1,98 @@
+zs <- file.path("extdata", "blobs.zarr")
+zs <- system.file(zs, package="spatialdataR")
+sd <- readSpatialData(zs)
+
+fn <- \(x, y) {
+    z <- paste(capture.output(show(x)), collapse="\n")
+    for (. in y) expect_match(z, .)
+}
+
+test_that("show(SpatialData)", {
+    # element counts
+    ni <- length(imageNames(sd))
+    nl <- length(labelNames(sd))
+    np <- length(pointNames(sd))
+    ns <- length(shapeNames(sd))
+    nt <- length(tableNames(sd))
+    
+    # coordinate systems
+    cg <- CTgraph(sd)
+    ts <- graph::nodeData(cg, graph::nodes(cg), "type")
+    cs <- graph::nodes(cg)[ts == "space"]
+    nc <- length(cs)
+    
+    # expected patterns 
+    ok <- c(
+        "class: SpatialData",
+        sprintf("- images\\(%d\\):", ni),
+        sprintf("- labels\\(%d\\):", nl),
+        sprintf("- points\\(%d\\):", np),
+        sprintf("- shapes\\(%d\\):", ns),
+        sprintf("- tables\\(%d\\):", nt),
+        sprintf("coordinate systems\\(%d\\):", nc))
+    
+    # include element names
+    el <- unname(unlist(colnames(sd[-5])))
+    ok <- c(ok, el, tableNames(sd))
+    
+    # add coordinate systems with element counts
+    for (c in cs) {
+        # check connectivity b/w spatial elements & coordinate systems
+        pa <- suppressWarnings(RBGL::sp.between(cg, paste0("_", el), c))
+        n <- sum(vapply(pa, \(.) !is.na(.$length), logical(1)))
+        ok <- c(ok, sprintf("- %s\\(%d\\):", c, n))
+    }
+    fn(sd, ok)
+})
+
+test_that("show(SpatialDataElement)", {
+    # image
+    x <- image(sd, 1)
+    ok <- c(
+        "class:  SpatialDataImage",
+        sprintf("Scales \\(%d\\):", length(data(x, NULL))),
+        sprintf("(%s)", paste(dim(x), collapse=",")))
+    fn(x, ok)
+    
+    # label
+    x <- label(sd, 1)
+    ok <- c(
+        "class:  SpatialDataLabel",
+        sprintf("Scales \\(%d\\):", length(data(x, NULL))),
+        sprintf("(%s)", paste(dim(x), collapse=",")))
+    fn(x, ok)
+    
+    # point
+    x <- point(sd, 1)
+    ok <- c(
+        "class: SpatialDataPoint",
+        sprintf("count: %d", length(x)),
+        sprintf("data\\(%d\\):", length(names(x))),
+        paste(names(x), collapse=" "))
+    fn(x, ok)
+    
+    # shape
+    x <- shape(sd, 1)
+    ok <- c(
+        "class: SpatialDataShape",
+        sprintf("count: %d", length(x)),
+        sprintf("data\\(%d\\):", length(names(x))),
+        paste(names(x), collapse=" "))
+    fn(x, ok)
+})
+
+test_that("show(SpatialDataAttrs)", {
+    l <- list(image(sd), label(sd))
+    for (x in l) {
+        x <- meta(x)
+        ok <- c(
+            "class: SpatialDataAttrs",
+            sprintf("axes\\(%d\\):", length(axes(x))),
+            sprintf("coordTrans\\(%d\\):", length(CTlist(x))))
+        if (!is.null(chs <- channels(x))) 
+            ok <- c(ok, sprintf("channels\\(%d\\):", length(chs)))
+        if (!is.null(ms <- multiscales(x)[[1]]))
+            ok <- c(ok, sprintf("datasets\\(%d\\):", length(ms$datasets)))
+        fn(x, ok)
+    }
+})
diff --git a/tests/testthat/test-sdarray.R b/tests/testthat/test-sdarray.R
index ce255c2a..35cf8629 100644
--- a/tests/testthat/test-sdarray.R
+++ b/tests/testthat/test-sdarray.R
@@ -31,7 +31,7 @@ test_that("SpatialDataImage()", {
     expect_silent(SpatialDataImage(list()))
     expect_silent(SpatialDataImage(mat))
     expect_silent(SpatialDataImage(list(mat)))
-    expect_silent(SpatialDataImage(list(mat), SpatialDataAttrs()))
+    expect_silent(SpatialDataImage(list(mat), SpatialDataAttrs(type="image")))
     # multiscale
     dim <- lapply(c(20, 10, 5), \(.) c(3, rep(., 2)))
     lys <- lapply(dim, \(.) array(rgb(prod(.)), dim=.))
@@ -54,7 +54,7 @@ test_that("data(),SpatialDataImage", {
 
 test_that("SpatialDataLabel()", {
     val <- sample(seq_len(12), 20*20, replace=TRUE)
-    mat <- array(val, dim=c(20, 20, 20))
+    mat <- array(val, dim=c(20, 20))
     # invalid
     expect_error(SpatialDataLabel(mat, 1))
     expect_error(SpatialDataLabel(mat, list()))
@@ -62,15 +62,15 @@ test_that("SpatialDataLabel()", {
     expect_silent(SpatialDataLabel(list()))
     expect_silent(SpatialDataLabel(mat))
     expect_silent(SpatialDataLabel(list(mat)))
-    expect_silent(SpatialDataLabel(list(mat), SpatialDataAttrs()))
+    expect_silent(SpatialDataLabel(list(mat), SpatialDataAttrs(type="label")))
     # multiscale
-    dim <- lapply(c(20, 10, 5), \(.) rep(., 3))
+    dim <- lapply(c(20, 10, 5), \(.) rep(., 2))
     lys <- lapply(dim, \(.) array(sample(seq_len(12), prod(.), replace=TRUE), dim=.))
     expect_silent(SpatialDataLabel(lys))
 })
 
 test_that("data(),SpatialDataLabel", {
-    dim <- lapply(c(8, 4, 2), \(.) rep(., 3))
+    dim <- lapply(c(8, 4, 2), \(.) rep(., 2))
     lys <- lapply(dim, \(.) array(0L, dim=.))
     lab <- SpatialDataLabel(lys)
     for (. in seq_along(lys))
diff --git a/tests/testthat/test-sdattrs.R b/tests/testthat/test-sdattrs.R
index 63cc1fca..8cdee9db 100644
--- a/tests/testthat/test-sdattrs.R
+++ b/tests/testthat/test-sdattrs.R
@@ -45,5 +45,64 @@ for (v in names(z)) {
         expect_silent(z <- channels(y <- image(x)))
         expect_length(z, dim(y)[1])
     })
-
 }
+
+test_that(".val_ome_ver()", {
+    # invalid
+    expect_error(.val_ome_ver(1))
+    expect_error(.val_ome_ver(TRUE))
+    expect_error(.val_ome_ver("0.0"))
+    expect_error(.val_ome_ver("0.30"))
+    expect_error(.val_ome_ver(c("0.3", "0.4")))
+    # valid
+    expect_silent(.val_ome_ver(v <- "0.3-x"))
+    expect_silent(x <- .val_ome_ver(v <- "0.3"))
+    expect_is(x, "character")
+    expect_length(x, 1)
+    expect_identical(x, v)
+})
+test_that("SpatialDataAttrs()", {
+    # invalid
+    expect_error(SpatialDataAttrs(nch=0))
+    expect_error(SpatialDataAttrs(dim=7))
+    expect_error(SpatialDataAttrs(ver="0.0"))
+    expect_error(SpatialDataAttrs(type="bad"))
+    # 2-4D image
+    nms <- c("c", "t", "z", "y", "x")
+    for (d in seq(2, 4)) {
+        x <- SpatialDataAttrs(type="image", dim=d, nch=7)
+        ok <- if (d == 2) nms[-c(2,3)] else if (d == 3) nms[-2] else nms
+        # axes name
+        y <- axes(x, "name")
+        expect_length(y, 1+d)
+        expect_is(y, "character")
+        expect_identical(y, ok)
+        # axes type
+        y <- axes(x, "type")
+        expect_equal(sum(y == "time"), ifelse(d == 4, 1, 0))
+        expect_equal(sum(y == "space"), ifelse(d == 2, 2, 3))
+        expect_equal(sum(y == "channel"), 1)
+        # channels
+        y <- channels(x)
+        expect_length(y, 7)
+        expect_is(y, "character")
+        expect_all_true(!duplicated(y))
+    }
+    # 2-4D label
+    for (d in seq(2, 4)) {
+        x <- SpatialDataAttrs(type="label", dim=d)
+        y <- axes(x, "type")
+        expect_length(y, d)
+        expect_equal(sum(y == "time"), ifelse(d == 4, 1, 0))
+        expect_equal(sum(y == "space"), ifelse(d == 2, 2, 3))
+    }
+    # 3-4D shape/point
+    for (d in seq(2, 4)) {
+        x <- SpatialDataAttrs(type="frame", dim=d)
+        y <- axes(x, "type")
+        expect_length(y, d)
+        expect_null(channels(x))
+        expect_equal(sum(y == "time"), ifelse(d == 4, 1, 0))
+        expect_equal(sum(y == "space"), ifelse(d == 2, 2, 3))
+    }
+})
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index 6457f9e5..e3d1bc22 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -6,7 +6,12 @@ x <- readSpatialData(x, tables=FALSE)
 
 # centroids ----
 
+test_that("centroids,invalid", {
+    expect_error(centroids(x), "supported")
+    expect_error(centroids(image(x)), "supported")
+})
 test_that("centroids,sdLabel", {
+    # 2D
     y <- label(x)
     z <- centroids(y, "data.frame")
     expect_is(z, "data.frame")
@@ -17,6 +22,15 @@ test_that("centroids,sdLabel", {
     expect_is(.z, "matrix")
     z$i <- as.integer(as.character(z$i))
     expect_identical(.z, as.matrix(z))
+    # 3D
+    n <- 7; m <- 8
+    u <- array(runif(2*n*m), c(1, 5, n, m))
+    v <- array(1L, c(5, n, m))
+    i <- SpatialDataImage(u, SpatialDataAttrs(dim=3, nch=1))
+    l <- SpatialDataLabel(v, SpatialDataAttrs(type="label", dim=3))
+    sd <- SpatialData(images=list(a=i), labels=list(b=l))
+    sd <- expect_silent(mask(sd, "a", "b", how="mean"))
+    expect_identical(as.numeric(assay(table(sd))), mean(u))
 })
 test_that("centroids,sdPoint", {
     i <- feature_key(y <- point(x))