diff --git a/DESCRIPTION b/DESCRIPTION
index a541151a..d91aad48 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: spatialdataR
 Title: Representation of Python's spatialdata in R
 Depends: R (>= 4.6)
-Version: 0.99.41
+Version: 0.99.42
 Description: R interface to Python/scverse's 'spatialdata' framework for 
     unified spatial omics data handling. Adheres to OME-NGFF standards, 
     providing lazy, on-disk representations for multiscale images and 
diff --git a/R/crop.R b/R/crop.R
index 79e7d58c..dc3f4b41 100644
--- a/R/crop.R
+++ b/R/crop.R
@@ -155,9 +155,11 @@ NULL
 
 #' @export
 #' @rdname crop
+#' @importFrom utils tail
 #' @importFrom methods is
 #' @importFrom sf st_bbox
 setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
+    #x <- label(sd); y <- bb; j <- 1
     if (is.matrix(y)) {
         y <- .check_pol(y)
         y <- st_bbox(st_polygon(list(y)))
@@ -193,12 +195,28 @@ setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
         wh[[2]] <- wh[[2]][1] + c(z$ymin, z$ymax)
     }
     metadata(x)$wh <- wh
+    # multi-scale adjustment
+    t <- .get_multiscale_scale(x)
+    tx <- tail(t, 1)
+    ty <- tail(t, 2)[1]
+    z$xmin <- floor(z$xmin/tx)
+    z$ymin <- floor(z$ymin/ty)
+    z$xmax <- ceiling(z$xmax/tx)
+    z$ymax <- ceiling(z$ymax/ty)
     # subset array
     i <- seq(z$ymin+1, z$ymax)
-    j <- seq(z$xmin+1, z$xmax) 
-    if (n == 3) x[, i, j] else x[i, j]
+    j <- seq(z$xmin+1, z$xmax)
+    ii <- is(x, "SpatialDataImage")
+    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/sdArray.R b/R/sdArray.R
index a916c87b..a299f8b1 100644
--- a/R/sdArray.R
+++ b/R/sdArray.R
@@ -63,6 +63,16 @@ 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
@@ -160,46 +170,53 @@ setMethod("channels", "SpatialDataElement", \(x, ...) stop("only 'images' have c
     )
 }
 
+# https://github.com/Huber-group-EMBL/Rarr/blob/1795c676e2ac81a9ba2a592c7210cc59036544b6/R/utils.R#L74-L79
+.sub <- \(x, ix) rlang::inject(x[!!!ix, drop=FALSE])
+
 #' @exportMethod [
 #' @rdname SpatialDataArray
 #' @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
+    ls <- seq_along(data(x, NULL))
+    data(x) <- lapply(ls, \(l) {
+        sf <- 2^(l-1)   
+        rc <- tail(dim(data(x, l)), 2) 
+        # get spatial indices
+        .yx <- lapply(seq_along(yx), \(a) {
+            ix <- yx[[a]]
+            if (isTRUE(ix)) return(seq_len(rc[a]))
+            if (is.numeric(ix)) {
+                return(seq(
+                    floor(head(ix, 1)/sf),
+                    min(ceiling(tail(ix, 1)/sf), rc[a])))
+            }
+            ix 
+        })
+        # combine leading & spatial indices
+        ix <- c(z, .yx)
+        # (optional) prepend additional indices
+        nd <- length(dim(data(x)))
+        na <- nd-length(ix)
+        if (na > 0) {
+            na <- !logical(na)
+            ix <- c(as.list(na), ix)
+        }
+        .sub(data(x, l), ix)
+    })
+    x
+}
+
 setMethod("[", "SpatialDataImage", \(x, i, j, k, ..., drop=FALSE) {
     if (missing(i)) i <- TRUE
     if (missing(j)) j <- TRUE else if (isFALSE(j)) j <- 0 else .check_jk(j, "j")
     if (missing(k)) k <- TRUE else if (isFALSE(k)) k <- 0 else .check_jk(k, "k")
-    ijk <- list(i, j, k)
-    n <- length(data(x, NULL))
-    d <- dim(data(x))
-    data(x) <- lapply(seq_len(n), \(.) {
-        j <- if (isTRUE(j)) seq_len(d[2]) else j
-        k <- if (isTRUE(k)) seq_len(d[3]) else k
-        jk <- lapply(list(j, k), \(jk) {
-            fac <- 2^(.-1)
-            seq(floor(head(jk, 1)/fac), 
-                ceiling(tail(jk, 1)/fac))
-        })
-        data(x, .)[i, jk[[1]], jk[[2]], drop=FALSE]
-    })
-    x
+    .sub_sda(x, yx=list(j, k), z=list(i))
 })
 
-#' @exportMethod [
-#' @rdname SpatialDataArray
-#' @importFrom utils head tail
 setMethod("[", "SpatialDataLabel", \(x, i, j, ..., drop=FALSE) {
     if (missing(i)) i <- TRUE else if (isFALSE(i)) i <- 0 else .check_jk(i, "i")
     if (missing(j)) j <- TRUE else if (isFALSE(j)) j <- 0 else .check_jk(j, "j")
-    n <- length(data(x, NULL))
-    d <- dim(data(x, 1))
-    data(x) <- lapply(seq_len(n), \(.) {
-        i <- if (isTRUE(i)) seq_len(d[1]) else i
-        j <- if (isTRUE(j)) seq_len(d[2]) else j
-        ij <- lapply(list(i, j), \(ij) {
-            fac <- 2^(.-1)
-            seq(floor(head(ij, 1)/fac), 
-                ceiling(tail(ij, 1)/fac))
-        })
-        data(x, .)[ij[[1]], ij[[2]], drop=FALSE]
-    })
-    x
+    .sub_sda(x, yx=list(i, j), z=list())
 })
diff --git a/README.md b/README.md
index f6e4c8a5..2287affe 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # SpatialData
 
-[![Bioc Check](https://github.com/HelenaLC/spatialdataR/actions/workflows/check-bioc.yml/badge.svg?branch=main&event=push)](https://github.com/HelenaLC/SpatialData/actions/workflows/check-bioc.yml)
+[![R-universe](https://github.com/HelenaLC/spatialdataR/actions/workflows/r-universe.yaml/badge.svg?branch=main&event=push)](https://github.com/HelenaLC/spatialdataR/actions/workflows/r-universe.yaml)
 
 `spatialdataR` provides an R interface to Python's [spatialdata](https://spatialdata.scverse.org) framework. 
 It enables the representation, handling, and integration of diverse spatial omics datasets 
diff --git a/inst/NEWS b/inst/NEWS
index 5f9fa546..171df13b 100644
--- a/inst/NEWS
+++ b/inst/NEWS
@@ -1,3 +1,8 @@
+changes in version 0.99.42
+
+- revised 'crop()' to adjust for array multiscales with scale factor != 1
+- revised subsetting to support >2/3D arrays (t & z dims. are protected)
+
 changes in version 0.99.41
 
 - fix vignette to call 'spatialdataR::transform()' explicitly
diff --git a/man/SpatialDataArray.Rd b/man/SpatialDataArray.Rd
index 4cbe62a8..7549c6e2 100644
--- a/man/SpatialDataArray.Rd
+++ b/man/SpatialDataArray.Rd
@@ -13,8 +13,7 @@
 \alias{channels,SpatialDataAttrs-method}
 \alias{channels,SpatialDataImage-method}
 \alias{channels,SpatialDataElement-method}
-\alias{[,SpatialDataImage,ANY,ANY,ANY-method}
-\alias{[,SpatialDataLabel,ANY,ANY,ANY-method}
+\alias{.sub_sda}
 \title{\code{SpatialDataArray}}
 \usage{
 SpatialDataImage(
@@ -45,9 +44,7 @@ SpatialDataLabel(
 
 \S4method{channels}{SpatialDataElement}(x, ...)
 
-\S4method{[}{SpatialDataImage,ANY,ANY,ANY}(x, i, j, k, ..., drop = FALSE)
-
-\S4method{[}{SpatialDataLabel,ANY,ANY,ANY}(x, i, j, ..., drop = FALSE)
+.sub_sda(x, yx, z = list())
 }
 \arguments{
 \item{data}{list of \code{ZarrArray}s}