diff --git a/NAMESPACE b/NAMESPACE index 07522e54..50257202 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,8 @@ export(SpatialDataImage) export(SpatialDataLabel) export(SpatialDataPoint) export(SpatialDataShape) +export(create_zarr) +export(create_zarr_group) export(filter) export(mutate) export(pull) @@ -26,6 +28,12 @@ export(readShape) export(readSpatialData) export(readTable) export(select) +export(writeImage) +export(writeLabel) +export(writePoint) +export(writeShape) +export(writeSpatialData) +export(writeTable) exportClasses(SpatialData) exportMethods("$") exportMethods("[") @@ -118,12 +126,15 @@ importFrom(BiocGenerics,colnames) importFrom(BiocGenerics,combine) importFrom(BiocGenerics,rownames) importFrom(DelayedArray,DelayedArray) +importFrom(EBImage,resize) importFrom(EBImage,rotate) importFrom(Matrix,sparseMatrix) importFrom(Matrix,sparseVector) importFrom(Matrix,summary) importFrom(RBGL,sp.between) importFrom(Rarr,read_zarr_attributes) +importFrom(Rarr,write_zarr_array) +importFrom(Rarr,write_zarr_attributes) importFrom(Rarr,zarr_overview) importFrom(S4Vectors,"metadata<-") importFrom(S4Vectors,coolcat) @@ -146,8 +157,10 @@ importFrom(ZarrArray,ZarrArray) importFrom(ZarrArray,path) importFrom(ZarrArray,type) importFrom(anndataR,read_zarr) +importFrom(anndataR,write_zarr) importFrom(dplyr,all_of) importFrom(dplyr,anti_join) +importFrom(dplyr,bind_cols) importFrom(dplyr,coalesce) importFrom(dplyr,collect) importFrom(dplyr,count) @@ -161,10 +174,12 @@ importFrom(dplyr,select) importFrom(dplyr,slice) importFrom(dplyr,sql) importFrom(dplyr,tally) +importFrom(dplyr,tibble) importFrom(duckspatial,as_duckspatial_df) importFrom(duckspatial,ddbs_bbox) importFrom(duckspatial,ddbs_intersects) importFrom(duckspatial,ddbs_open_dataset) +importFrom(duckspatial,ddbs_write_dataset) importFrom(graph,"edgeData<-") importFrom(graph,"edgeDataDefaults<-") importFrom(graph,"nodeData<-") @@ -196,6 +211,7 @@ importFrom(sf,st_geometry_type) importFrom(sf,st_polygon) importFrom(sf,st_sf) importFrom(sf,st_sfc) +importFrom(stats,setNames) importFrom(utils,.DollarNames) importFrom(utils,head) importFrom(utils,tail) diff --git a/R/AllClasses.R b/R/AllClasses.R index bbaafdb3..d63806a9 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -2,6 +2,20 @@ Class="SpatialDataAttrs", contains="list") +.sdFormat <- setClass( + Class = "sdFormat", + slots = list( + version = "character", + zarr_version = "integer", + ome_version = "character", + image = "character", + label = "character", + point = "character", + shape = "character", + table = "character" + ) +) + #' @importFrom methods setClassUnion #' @importClassesFrom S4Arrays Array setClassUnion( diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 7936cba1..dc365855 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -107,9 +107,17 @@ setGeneric("centroids", \(x, ...) standardGeneric("centroids")) setGeneric("data_type", \(x, ...) standardGeneric("data_type")) setGeneric("geom_type", \(x, ...) standardGeneric("geom_type")) setGeneric("multiscales", \(x, ...) standardGeneric("multiscales")) +setGeneric("datasets", \(x, ...) standardGeneric("datasets")) # tbl ---- setGeneric("hasTable", \(x, i, ...) standardGeneric("hasTable")) setGeneric("getTable", \(x, i, ...) standardGeneric("getTable")) setGeneric("setTable", \(x, i, ...) standardGeneric("setTable")) + +# zarr ---- + +setGeneric("version", \(x, ...) standardGeneric("version")) +setGeneric("version<-", \(x, value) standardGeneric("version<-")) +setGeneric("zarr_version", \(x, ...) standardGeneric("zarr_version")) +setGeneric("ome_version", \(x, ...) standardGeneric("ome_version")) diff --git a/R/SpatialData.R b/R/SpatialData.R index 204a0ca0..b998cf28 100644 --- a/R/SpatialData.R +++ b/R/SpatialData.R @@ -28,6 +28,7 @@ #' @param x \code{SpatialData} #' @param i,j character string, scalar or vector of indices #' specifying the element to extract from a given layer. +#' #' @param drop ignored. #' @param name character string for extraction (see \code{?base::`$`}). #' @param value (list of) element(s) with layer-compliant object(s), @@ -51,11 +52,11 @@ #' #' @export SpatialData <- \(images, labels, points, shapes, tables) { - if (missing(images)) images <- list() - if (missing(labels)) labels <- list() - if (missing(points)) points <- list() - if (missing(shapes)) shapes <- list() - if (missing(tables)) tables <- list() + if (missing(images)) images <- setNames(list(), character(0)) + if (missing(labels)) labels <- setNames(list(), character(0)) + if (missing(points)) points <- setNames(list(), character(0)) + if (missing(shapes)) shapes <- setNames(list(), character(0)) + if (missing(tables)) tables <- setNames(list(), character(0)) .SpatialData( images=images, labels=labels, points=points, shapes=shapes, tables=tables) diff --git a/R/format.R b/R/format.R new file mode 100644 index 00000000..4cd0499e --- /dev/null +++ b/R/format.R @@ -0,0 +1,49 @@ +#' @name sdFormat +#' @title The `sdFormat` class +#' +#' @param version SpatialData version: 0.1 or 0.2. +#' +#' @details +#' +#' @return \code{sdFormat} +#' +#' @noRd +sdFormat <- function(version = "0.1") { + switch(as.character(version), + "0.2" = { + .sdFormat( + version = "0.2", + zarr_version = 3L, + ome_version = "0.5", + image = "0.3", + label = "0.3", + shape = "0.3", + point = "0.2", + table = "0.2" + ) + }, + "0.1" = { + .sdFormat( + version = "0.1", + zarr_version = 2L, + ome_version = "0.5", + image = "0.2", + label = "0.2", + shape = "0.2", + point = "0.1", + table = "0.1" + ) + }, + stop("Incorrect SpatialData version. Must be 0.1 or 0.2!") + ) +} + +setMethod("image", "sdFormat", \(x) x@image) +setMethod("label", "sdFormat", \(x) x@label) +setMethod("shape", "sdFormat", \(x) x@shape) +setMethod("point", "sdFormat", \(x) x@point) +setMethod("table", "sdFormat", \(x) x@table) +setMethod("zarr_version", "sdFormat", \(x) x@zarr_version) +setMethod("zarr_version", "character", \(x) zarr_version(sdFormat(x))) +setMethod("ome_version", "sdFormat", \(x) x@ome_version) +setMethod("version", "sdFormat", \(x) x@version) \ No newline at end of file diff --git a/R/read.R b/R/read.R index 355ae1fd..23203383 100644 --- a/R/read.R +++ b/R/read.R @@ -56,14 +56,16 @@ NULL #' @export readImage <- function(x, ...) { l <- .readArray(x, ...) - SpatialDataImage(data=l$array, meta=SpatialDataAttrs(l$md), ...) + SpatialDataImage(data=l$array, meta=SpatialDataAttrs(l$md), + version = version(l$md), ...) } #' @rdname readSpatialData #' @export readLabel <- function(x, ...) { l <- .readArray(x, ...) - SpatialDataLabel(data=l$array, meta=SpatialDataAttrs(l$md), ...) + SpatialDataLabel(data=l$array, meta=SpatialDataAttrs(l$md), + version = version(l$md), ...) } #' @rdname readSpatialData @@ -79,7 +81,7 @@ readPoint <- function(x, ...) { mutate(geometry=sql(sprintf("ST_Point(%s, %s)", ax[1], ax[2]))) |> as_duckspatial_df(crs=NA_character_) |> select(-all_of(ax)) - SpatialDataPoint(data=df, meta=SpatialDataAttrs(md)) + SpatialDataPoint(data=df, meta=SpatialDataAttrs(md), version = version(md)) } #' @rdname readSpatialData @@ -90,7 +92,8 @@ readPoint <- function(x, ...) { readShape <- function(x, ...) { md <- read_zarr_attributes(x) pq <- list.files(x, "\\.parquet$", full.names=TRUE) - SpatialDataShape(data=ddbs_open_dataset(pq), meta=SpatialDataAttrs(md)) + SpatialDataShape(data=ddbs_open_dataset(pq), meta=SpatialDataAttrs(md), + version = version(md)) } #' @export @@ -105,7 +108,7 @@ readTable <- function(x) { }) # move these to 'int_metadata' nm <- "spatialdata_attrs" - md <- metadata(sce)[[nm]] + md <- read_zarr_attributes(x) int_metadata(sce)[[nm]] <- md metadata(sce)[[nm]] <- NULL # move these to 'int_colData' diff --git a/R/sdArray.R b/R/sdArray.R index 63b72a4f..964a419e 100644 --- a/R/sdArray.R +++ b/R/sdArray.R @@ -67,20 +67,59 @@ NULL #' @rdname SpatialDataArray #' @importFrom methods new #' @importFrom S4Vectors metadata<- -SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) { - x <- .SpatialDataImage(data=data, meta=meta, ...) - metadata(x) <- metadata - return(x) +SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(), + version = image(sdFormat(0.1)), + metadata=list(), + scale_factors = NULL, ...) { + if(!is.list(data)) + data <- list(data) + if(!is.null(scale_factors)){ + data <- .generate_multiscale(data[[1]], + axes = vapply(axes(meta), + \(.) .$name, + character(1)), + scale_factors = scale_factors, + method = "image") + # TODO: this supposed to update the scale_factors not write a new meta + meta <- SpatialDataAttrs(scale_factors = scale_factors) + } + # construct S4 object + x <- .SpatialDataImage(data=data, meta=meta, ...) + metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version + return(x) } #' @export #' @rdname SpatialDataArray #' @importFrom methods new #' @importFrom S4Vectors metadata<- -SpatialDataLabel <- function(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) { - x <- .SpatialDataLabel(data=data, meta=meta, ...) - metadata(x) <- metadata - return(x) +SpatialDataLabel <- function(data=list(), + meta=SpatialDataAttrs(label = TRUE), + version = image(sdFormat(0.1)), + metadata=list(), + scale_factors = NULL, ...) { + if(!is.list(data)) + data <- list(data) + if(!is.null(scale_factors)){ + data <- .generate_multiscale(data[[1]], + axes = vapply(axes(meta), + \(.) .$name, + character(1)), + scale_factors = scale_factors, + method = "label") + meta <- SpatialDataAttrs(scale_factors = scale_factors, label = TRUE) + } + x <- .SpatialDataLabel(data=data, meta=meta, ...) + metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version + return(x) } # utils ---- @@ -112,7 +151,7 @@ setMethod("length", "SpatialDataArray", \(x) length(data(x, NULL))) #' @importFrom S4Vectors metadata setMethod("data_type", "SpatialDataArray", \(x) { if (is(y <- data(x), "DelayedArray")) - data_type(y) else metadata(x)$data_type + data_type(y) else metadata(x)$data_type }) #' @export @@ -125,6 +164,74 @@ setMethod("data_type", "DelayedArray", \(x) { return(df$data_type) }) +#' @importFrom S4Vectors isSequence +.get_multiscales_paths <- function(x) { + ps <- list.files(x) + ps <- suppressWarnings(as.numeric(sort(ps, decreasing=FALSE))) + ps <- ps[!is.na(ps)] + if (length(ps)) { + qs <- seq(min(ps), max(ps)) + if (!isTRUE(all.equal(ps, qs))) + stop("SpatialDataImage paths are ill-defined, should", + " be an integer sequence, e.g., 0,1,...,n") + } else { + stop("SpatialDataImage path is empty") + } + return(ps) +} + +#' .create_mip +#' +#' Generate a downsampled pyramid of images. +#' +#' @param image image +#' @param scale_factors +#' +#' @importFrom EBImage resize +#' @importFrom stats setNames +#' +#' @inheritParams write_image +#' +#' @noRd +.generate_multiscale <- function(image, + scale_factors = c(2,2,2,2), + axes, + method = "image"){ + + # check dim + ndim <- length(dim(image)) + if (ndim > 3) { + stop("Only images of 5D or less are supported") + } + + # get x y dimensions for EBImage + dim_image <- stats::setNames(dim(image), axes) + dim_image <- dim_image[c("x", "y")] + + # downscale image + image_list <- list(image) + cur_image <- aperm(image, + perm = rev(seq_along(axes))) + for (i in seq_along(scale_factors)) { + dim_image <- ceiling(dim_image / scale_factors[i]) + image_list[[i+1]] <- + aperm(EBImage::resize(cur_image, + w = dim_image[1], + h = dim_image[2], + filter = switch(method, + image = "bilinear", + label = "none")), + perm = rev(seq_along(axes))) + } + if (method == "label") { + image_list <- lapply(image_list, function(x) { + storage.mode(x) <- "integer" + x + }) + } + image_list +} + # chs ---- # internal use only! @@ -146,34 +253,18 @@ setMethod("channels", "SpatialDataImage", \(x, ...) channels(meta(x))) #' @rdname SpatialDataArray setMethod("channels", "SpatialDataElement", \(x, ...) stop("only 'images' have channels")) -#' @importFrom S4Vectors isSequence -.get_multiscales_paths <- function(x) { - ps <- list.files(x) - ps <- suppressWarnings(as.numeric(sort(ps, decreasing=FALSE))) - ps <- ps[!is.na(ps)] - if (length(ps)) { - qs <- seq(min(ps), max(ps)) - if (!isTRUE(all.equal(ps, qs))) - stop("SpatialDataImage paths are ill-defined, should", - " be an integer sequence, e.g., 0,1,...,n") - } else { - stop("SpatialDataImage path is empty") - } - return(ps) -} - # sub ---- .check_jk <- \(x, .) { - if (isTRUE(x)) return() - tryCatch( - stopifnot( - is.numeric(x), x == round(x), - diff(range(x)) == length(x)-1, - (y <- abs(x)) == seq(min(y), max(y)) - ), - error=\(e) stop(sprintf("invalid '%s'", .)) - ) + if (isTRUE(x)) return() + tryCatch( + stopifnot( + is.numeric(x), x == round(x), + diff(range(x)) == length(x)-1, + (y <- abs(x)) == seq(min(y), max(y)) + ), + error=\(e) stop(sprintf("invalid '%s'", .)) + ) } #' @exportMethod [ @@ -218,4 +309,4 @@ setMethod("[", "SpatialDataLabel", \(x, i, j, ..., drop=FALSE) { data(x, .)[ij[[1]], ij[[2]], drop=FALSE] }) x -}) +}) \ No newline at end of file diff --git a/R/sdAttrs.R b/R/sdAttrs.R index 8944f798..0ac37561 100644 --- a/R/sdAttrs.R +++ b/R/sdAttrs.R @@ -59,63 +59,104 @@ #' #' @export SpatialDataAttrs <- \(x, type=c("array", "frame"), - label=FALSE, trans=NULL, ver="0.4", nch=3, ...) + label=FALSE, trans=NULL, ver="0.4", nch=3, scale_factors = NULL, ...) { 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) - } + ax <- .default_ax(type, label) # transformations: ct <- trans %||% .default_ct(ax) + ds <- .default_ds(.ax_names(ax), scale_factors) # .zattrs list: if (type == "array") { - # default structure - res <- list( - omero=list(channels=list(label=letters[seq_len(nch)])), - multiscales=list(list( - axes=ax, - version="0.4", - coordinateTransformations=ct, - datasets=list(list(path="0", coordinateTransformations=list(list(type="scale", scale=list(1, 1)))))))) - if (ver == "0.3") res <- list(ome=res) + # default structure + res <- list() + if(!label) + res <- c(res, + list(omero=list(channels=lapply(letters[seq_len(nch)], + \(.) list(label = .))))) + res <- c(res, + list( + multiscales= + list( + list( + axes=ax, + version="0.4", + coordinateTransformations=ct, + datasets=ds + ) + ) + ) + ) + if (ver == "0.3") res <- list(ome=res) } else { - # points/shapes - res <- list(axes=ax, coordinateTransformations=ct) + # points/shapes + res <- list(axes=ax, coordinateTransformations=ct) } res$spatialdata_attrs <- list(version=ver) SpatialDataAttrs(res) } # 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("array", "frame"), label = FALSE) { + switch(match.arg(type), + # (c)yx for images/labels + array={ + ax <- list( + list(name="x", type="space"), + list(name="y", type="space")) + if (type == "array") { + # yx for labels + ax <- rev(ax) + # yx for images, cyx if requested + if (!label) ax <- c(list(list(name="c", type="channel")), ax) + } + ax + }, + # xy for points/shapes + list("x", "y")) +} + +.ax_names <- function(ax){ + if (is.character(ax[[1]])) { + unlist(ax) + } else { + vapply(ax, \(.) .$name, character(1)) + } } # Internal helper to generate coordinate transformations .default_ct <- \(axes, name="global", type="identity", data=NULL) { - ct <- list(input=axes, output=list(name=name), type=type) - if (!is.null(data)) ct[[type]] <- data - list(ct) + ct <- list( + input=list(axes = axes, + name = if(length(axes) == 3) "cyx" else "yx"), + output=list(axes = axes, + name = name), + type = type) + if (!is.null(data)) ct[[type]] <- data + list(ct) } +.default_ds <- function(axes, scale_factors = NULL){ + scale_factors <- cumprod(c(1,scale_factors)) + paths <- paste0(seq_along(scale_factors) - 1) + mapply(\(p,s) { + list( + coordinateTransformations = list( + list( + scale = lapply( + axes, + \(.) if(. == "c") 1 else s), + type = "scale" + ) + ), + path = p + ) + }, paths, scale_factors, USE.NAMES = FALSE, SIMPLIFY = FALSE) +} + + #' @export #' @importFrom utils .DollarNames .DollarNames.SpatialDataAttrs <- \(x, pattern="") names(x) @@ -126,21 +167,19 @@ setMethod("$", "SpatialDataAttrs", \(x, name) x[[name]]) # internal use only! #' @noRd -.zv <- \(x) { - v <- x$spatialdata_attrs$version - if (!length(v)) stop("couldn't find 'version' in 'spatialdata_attrs'") - ok <- length(v) == 1 && is.character(v) && v %in% sprintf("0.%d", seq_len(5)) - if (!ok) stop("invalid 'version' in 'spatialdata_attrs'; expected '0.x' where x is 1-5") - return(v) -} +.ms <- \(x) switch(version(x), "0.3"=x$ome$multiscales, x$multiscales) # internal use only! #' @noRd -.ms <- \(x) switch(.zv(x), "0.3"=x$ome$multiscales, x$multiscales) +setMethod("multiscales", "list", .ms) # internal use only! #' @noRd -setMethod("multiscales", "list", .ms) +setMethod("datasets", "list", \(x, ...) { + vapply(multiscales(x)[[1]]$datasets, \(.){ + .$path + }, character(1)) +}) # features ---- @@ -300,3 +339,62 @@ setReplaceMethod("instances", c("SingleCellExperiment", "ANY"), \(x, value) { int_colData(x)[[ik]] <- value return(x) }) + +# elements ---- + +setMethod("version", c("SpatialDataElement"), \(x) { + version(meta(x)) +}) + +setMethod("version", c("SingleCellExperiment"), \(x) { + meta(x)$version +}) + +setMethod("version", "SpatialDataAttrs", \(x) .zv(x)) + +setMethod("version", "list", \(x) .zv(x)) + +.zv <- \(x) { + v <- x$spatialdata_attrs$version + if (!length(v)) stop("couldn't find 'version' in 'spatialdata_attrs'") + ok <- length(v) == 1 && is.character(v) && v %in% sprintf("0.%d", seq_len(5)) + if (!ok) stop("invalid 'version' in 'spatialdata_attrs'; expected '0.x' where x is 1-5") + return(v) +} + +setReplaceMethod("version", c("SpatialDataFrame"), \(x, value) { + if(!value %in% c("0.1", "0.2", "0.3")) + stop("Unknown version for shape/point! Must be 0.2 or 0.3.") + meta(x)$spatialdata_attrs$version <- value + x +}) + +setReplaceMethod("version", c("SpatialDataArray"), \(x, value) { + mt <- meta(x) + if(value == "0.3"){ + if(is.null(mt$ome)){ + mt$ome = list(omero = mt$omero, + multiscales = mt$multiscales) + mt$omero <- NULL + mt$multiscales <- NULL + } + } else if(value %in% c("0.1" ,"0.2")){ + if(is.null(mt$multiscales)){ + mt$omero <- mt$ome$omero + mt$multiscales <- mt$ome$multiscales + mt[["ome"]] <- NULL + } + } else { + stop("Unknown version for image/label! Must be 0.1, 0.2, 0.3.") + } + mt$spatialdata_attrs$version <- value + meta(x) <- mt + x +}) + +setReplaceMethod("version", c("SingleCellExperiment"), \(x, value) { + if(!value %in% c("0.1", "0.2")) + stop("Unknown version for table! Must be 0.1 or 0.2") + int_metadata(x)$spatialdata_attrs$version <- value + return(x) +}) \ No newline at end of file diff --git a/R/sdFrame.R b/R/sdFrame.R index 5488392a..9274d358 100644 --- a/R/sdFrame.R +++ b/R/sdFrame.R @@ -104,7 +104,9 @@ NULL #' @importFrom sf st_geometry_type #' @importFrom S4Vectors metadata<- #' @importFrom duckspatial as_duckspatial_df -SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=list(), ik=NULL, fk=NULL, ...) { +SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), + version = point(sdFormat(0.1)), + metadata=list(), ik=NULL, fk=NULL, ...) { data <- .df_to_sf(data, "POINT") # validate geometry type (must be points) if (isTRUE(nrow(data) > 0L)) { @@ -131,6 +133,10 @@ SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=l # construct S4 object x <- .SpatialDataPoint(data=data, meta=SpatialDataAttrs(za), ...) metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version return(x) } @@ -139,7 +145,9 @@ SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=l #' @importFrom methods is #' @importFrom S4Vectors metadata<- #' @importFrom duckspatial as_duckspatial_df -SpatialDataShape <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=list(), ...) { +SpatialDataShape <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), + version = shape(sdFormat(0.1)), + metadata=list(), ...) { data <- .df_to_sf(data, "POLYGON") # always ensure internal data is 'duckspatial_df' if (isTRUE(nrow(data) > 0L) && @@ -147,6 +155,10 @@ SpatialDataShape <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=l data <- as_duckspatial_df(data, crs=NA) x <- .SpatialDataShape(data=data, meta=meta, ...) metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version return(x) } diff --git a/R/tables.R b/R/tables.R index ae593733..ee982daa 100644 --- a/R/tables.R +++ b/R/tables.R @@ -158,8 +158,7 @@ setMethod("setTable", c("SpatialData", "ANY"), \(x, i, ..., name=NULL, rk="rk", #' @importFrom SingleCellExperiment SingleCellExperiment int_colData int_colData<- int_metadata<- #' @export setMethod("setTable", c("SpatialData", "character"), \(x, i, y, - name=NULL, rk="region", ik="instance_id") { - + name=NULL, rk="region", ik="instance_id", version = "0.1") { # validity stopifnot( is(y, "SingleCellExperiment"), @@ -184,6 +183,8 @@ setMethod("setTable", c("SpatialData", "character"), \(x, i, y, if (is.null(region_key(y))) region_key(y) <- rk if (is.null(instance_key(y))) instance_key(y) <- ik + version(y) <- version + if (is.null(region(y))) { regions(y) <- i } else { diff --git a/R/write.R b/R/write.R new file mode 100644 index 00000000..6726fec2 --- /dev/null +++ b/R/write.R @@ -0,0 +1,245 @@ +#' @name writeSpatialData +#' @title Writing `SpatialData` +#' +#' @aliases +#' writeSpatialData +#' writeImage writeLabel +#' writePoint writeShape writeTable +#' +#' @param x +#' For \code{writeSpatialData}, +#' a \code{SpatialData} +#' For \code{writeImage/Label/Point/Shape/Table}, +#' a \code{ImageArray},\code{LabelArray}, +#' \code{PointFrame}, \code{ShapeFrame} +#' @param path path to zarr store. +#' @param replace if TRUE, existing elements with the same name will be +#' replaced with the given element +#' @param version SpatialData version, 0.1 (zarr v2) or 0.2 (zarr v3) +#' @param ... option arguments passed to and from other methods. +#' +#' @return +#' \itemize{ +#' \item{For \code{writeSpatialData}, a \code{SpatialData}.}, +#' \item{For element writers, a \code{ImageArray}, \code{LabelArray}, +#' \code{PointFrame}, \code{ShapeFrame}, or \code{SingleCellExperiment}.}} +#' +NULL + +#' @rdname writeSpatialData +#' @export +writeSpatialData <- function(x, path, replace = TRUE, version = "0.2", + ...) { + format <- sdFormat(version) + zarr.path <- .replace_zarr(path, + replace, + version = zarr_version(format)) + + # write root-level spatialdata_attrs for v3 (Python uses this to pick the read path) + if (version == "0.2") + Rarr::write_zarr_attributes( + zarr.path, + new.zattrs = list( + spatialdata_attrs = list(version = version), + spatialdata_software_version = + paste0("SpatialData v", packageVersion("SpatialData")) + ) + ) + + # write points + . <- lapply(pointNames(x), \(.){ + writePoint(point(x, .),., path = zarr.path, + replace = replace, format = format) + }) + + # write shapes + . <- lapply(shapeNames(x), \(.){ + writeShape(shape(x, .),., path = zarr.path, + replace = replace, format = format) + }) + + # write images + . <- lapply(imageNames(x), \(.){ + writeImage(image(x, .),., path = zarr.path, + replace = replace, format = format) + }) + + # write labels + . <- lapply(labelNames(x), \(.){ + writeLabel(label(x, .),., path = zarr.path, + replace = replace, format = format) + }) + + # write tables + . <- lapply(tableNames(x), \(.){ + writeTable(table(x, .),., path = zarr.path, + replace = replace, format = format) + }) +} + +#' @rdname writeSpatialData +#' @export +writePoint <- function(x, name, path, replace = TRUE, + format = sdFormat("0.1")) { + + # if no PointFrames were written before, update zarr store + zarr.group <- .make_zarr_group(x, name, + file.path(path, "points"), + replace, + version = zarr_version(format)) + + # write meta + Rarr::write_zarr_attributes(zarr.group, new.zattrs = meta(x)) + + # version + version(x) <- point(format) + + # write data + arrow::write_dataset(.point_to_xy(data(x)), + file.path(zarr.group, "points.parquet"), + basename_template = "part.{i}.parquet") +} + +#' @importFrom dplyr bind_cols tibble +.point_to_xy <- function(data) { + data %>% + st_as_sf() %>% + { + coords <- st_coordinates(.) + + bind_cols( + tibble( + x = coords[,1], + y = coords[,2] + ), + . + ) + } %>% + select(-geometry) +} + +#' @rdname writeSpatialData +#' @importFrom duckspatial ddbs_write_dataset +#' @importFrom Rarr write_zarr_attributes +#' @export +writeShape <- function(x, name, path, replace = TRUE, + format = sdFormat("0.1")) { + + # if no ShapeFrames were written before, update zarr store + zarr.group <- .make_zarr_group(x, name, + file.path(path, "shapes"), + replace, + version = zarr_version(format)) + + # write meta + Rarr::write_zarr_attributes(zarr.group, new.zattrs = meta(x)) + + # version + version(x) <- shape(format) + + # write data as a single parquet file (matches Python spatialdata convention) + duckspatial::ddbs_write_dataset( + data(x), + file.path(zarr.group, "shapes.parquet"), + overwrite = TRUE, + quiet = TRUE + )} + +#' @rdname writeSpatialData +#' @importFrom Rarr write_zarr_array write_zarr_attributes +#' @export +writeImage <- function(x, name, path, replace = TRUE, + format = sdFormat("0.1")) { + + # if no ImageArray were written before, update zarr store + zarr.group <- .make_zarr_group(x, name, + file.path(path, "images"), + replace, + version = zarr_version(format)) + + # write meta: + Rarr::write_zarr_attributes(zarr.group, new.zattrs = meta(x)) + + # version + version(x) <- image(format) + + # write data + dimension_names <- vapply(axes(meta(x)), \(.) .$name, character(1)) + lapply( + as.numeric(datasets(meta(x))), + \(.){ + arr <- realize(data(x, . + 1)) + # Rarr reads names(dimnames(x)) to write dimension_names in v3 zarr.json + if (!is.null(dimension_names)) + dimnames(arr) <- setNames(vector("list", length(dim(arr))), dimension_names) + Rarr::write_zarr_array(arr, + zarr_array_path = file.path(zarr.group, .), + chunk_dim = dim(arr), + order = "C", + dimension_separator = "/", + zarr_version = zarr_version(format)) + } + ) +} + +#' @rdname writeSpatialData +#' @importFrom Rarr write_zarr_array write_zarr_attributes +#' @export +writeLabel <- function(x, name, path, replace = TRUE, + format = sdFormat("0.1")) { + + # if no LabelArray were written before, update zarr store + zarr.group <- .make_zarr_group(x, name, + file.path(path, "labels"), + replace, + version = zarr_version(format)) + + # write meta: + Rarr::write_zarr_attributes(zarr.group, new.zattrs = meta(x)) + + # version + version(x) <- label(format) + + # write data + dimension_names <- vapply(axes(meta(x)), \(.) .$name, character(1)) + lapply( + as.numeric(datasets(meta(x))), + \(.){ + arr <- realize(data(x, . + 1)) + if (!is.null(dimension_names)) + dimnames(arr) <- setNames(vector("list", length(dim(arr))), dimension_names) + Rarr::write_zarr_array(arr, + zarr_array_path = file.path(zarr.group, .), + chunk_dim = dim(arr), + order = "C", + dimension_separator = "/", + zarr_version = zarr_version(format)) + } + ) +} + +#' @rdname writeSpatialData +#' @importFrom Rarr write_zarr_attributes +#' @importFrom anndataR write_zarr +#' @export +writeTable <- function(x, name, path, replace = TRUE, + format = sdFormat("0.1")) { + + # if no Table were written before, update zarr store + zarr.group <- .make_zarr_group(x, name, + file.path(path, "tables"), + replace, + version = zarr_version(format)) + + # write meta: + Rarr::write_zarr_attributes(zarr.group, new.zattrs = meta(x)) + + # version + version(x) <- table(format) + + # write data + if(zarr_version(format) == 3) + stop("Write support for anndata v3 zarr is not supported yet!") + anndataR::write_zarr(x, path = zarr.group, mode = "a") +} + diff --git a/R/zarr_utils.R b/R/zarr_utils.R new file mode 100644 index 00000000..31ee0381 --- /dev/null +++ b/R/zarr_utils.R @@ -0,0 +1,154 @@ +#' create_zarr_group +#' +#' create zarr groups +#' +#' @param store the location of (zarr) store +#' @param name name of the group +#' @param version zarr version +#' @export +create_zarr_group <- function(store, name, version = 2){ + split.name <- strsplit(name, split = "\\/")[[1]] + if(length(split.name) > 1){ + split.name <- vapply(seq_along(split.name), + function(x) paste(split.name[seq_len(x)], collapse = "/"), + FUN.VALUE = character(1)) + split.name <- rev(tail(split.name,2)) + if(!dir.exists(file.path(store,split.name[2]))) + create_zarr_group(store = store, name = split.name[2], version = version) + } + dir.create(file.path(store, split.name[1]), showWarnings = FALSE) + switch(as.character(version), + "2" = { + write("{\"zarr_format\":2}", file = file.path(store, split.name[1], ".zgroup"))}, + "3" = { + write( + "{\"zarr_format\":3,\"node_type\":\"group\",\"attributes\":{}}", + file = file.path(store, split.name[1], "zarr.json")) + }, + stop("version must be '2' or '3'") + ) + +} + +#' create_zarr +#' +#' Create Zarr store +#' +#' @param store The location of the Zarr store +#' @param version Zarr version +#' +#' @return `NULL` +#' +#' @examples +#' store <- tempfile(fileext = ".zarr") +#' create_zarr(store = store) +#' dir.exists(store) +#' +#' @export +create_zarr <- function(store, version = 2) { + prefix <- basename(store) + dir <- gsub(paste0(prefix, "$"), "", store) + create_zarr_group(store = dir, name = prefix, version = version) +} + +.replace_zarr <- function(zarr.path, replace, version = 2) +{ + if (dir.exists(zarr.path) && !replace) + stop("zarr store with name ", zarr.path ," doesnt exist") + if (!replace) + stop("Directory \"", zarr.path, "\" already exists. ", + "Use 'replace=TRUE' to replace it. ", + "Its content will be lost!") + if (unlink(zarr.path, recursive=TRUE) != 0L) + stop("failed to delete directory \"", dir, "\"") + create_zarr(zarr.path, version = version) + return(zarr.path) +} + +.make_zarr_group <- function(x, name, path, replace, version){ + + # create element parent dir + if(!dir.exists(path)) + dir.create(path) + + # check element dir + ng <- file.path(path, name) + if(replace){ + unlink(ng, recursive = TRUE) + } else { + nms <- list.dirs(file.path(path), full.names = FALSE) + if(name %in% nms) + stop("Directory \"", ng, "\" already exists. ", + "Use 'replace=TRUE' to replace it. ", + "Its content will be lost!") + } + + # create group + create_zarr_group(path, name, version) + + return(ng) +} + + +# For zarr v3, OME-NGFF content (multiscales, omero, image-label) must be +# nested under an "ome" key inside "attributes"; spatialdata_attrs stays at top. +# If the metadata was read from a v3 store it already has "ome", so skip wrapping. +.wrap_ome_for_v3 <- function(zattrs, version) { + if (version != "v3" || "ome" %in% names(zattrs)) return(as.list(zattrs)) + ome_keys <- setdiff(names(zattrs), "spatialdata_attrs") + ome_content <- as.list(zattrs)[ome_keys] + # Strip v2-only fields from each multiscales entry + if (!is.null(ome_content$multiscales)) { + ome_content$multiscales <- lapply(ome_content$multiscales, function(ms) { + ms[setdiff(names(ms), c("version", "metadata"))] + }) + } + list( + ome = c(list(version = "0.5-dev-spatialdata"), ome_content), + spatialdata_attrs = zattrs[["spatialdata_attrs"]] + ) +} + +# .get_multiscale_axes <- function(zattrs) { +# multiscales <- zattrs[["multiscales"]] +# if (is.null(multiscales) && !is.null(zattrs[["ome"]])) +# multiscales <- zattrs[["ome"]][["multiscales"]] +# if (is.null(multiscales) || length(multiscales) == 0L) return(NULL) +# axes <- multiscales[[1]][["axes"]] +# if (is.null(axes) || length(axes) == 0L) return(NULL) +# vapply(axes, `[[`, character(1), "name") +# } + +# # Post-processes Rarr-written v3 array zarr.json: +# # 1. Sorts codecs to required order [array-array → array-bytes → bytes-bytes]. +# # Rarr currently serialises them as [transpose, zstd, bytes] which Python rejects. +# # 2. Adds "attributes": {} and "storage_transformers": [] which Python zarr expects +# # but Rarr does not emit. +# # dimension_names are handled upstream by setting names(dimnames()) before write_zarr_array. +# .normalize_v3_array_metadata <- function(zarr_array_path) { +# metadata_path <- file.path(zarr_array_path, "zarr.json") +# if (!file.exists(metadata_path)) return(invisible(FALSE)) +# +# metadata <- jsonlite::read_json(metadata_path, simplifyVector = FALSE) +# codecs <- metadata[["codecs"]] +# if (!is.null(codecs) && length(codecs) > 1L) { +# codec_names <- vapply(codecs, `[[`, character(1), "name") +# codec_stage <- ifelse( +# codec_names %in% "transpose", 1L, +# ifelse(codec_names %in% c("bytes", "vlen-utf8", "vlen_utf8"), 2L, 3L) +# ) +# metadata[["codecs"]] <- codecs[order(codec_stage)] +# } +# +# if (is.null(metadata[["attributes"]])) metadata[["attributes"]] <- list() +# if (is.null(metadata[["storage_transformers"]])) metadata[["storage_transformers"]] <- list() +# +# jsonlite::write_json( +# metadata, +# path = metadata_path, +# auto_unbox = TRUE, +# pretty = 4, +# null = "null" +# ) +# invisible(TRUE) +# } diff --git a/inst/scripts/legacy/PointFrame.R b/inst/scripts/legacy/PointFrame.R new file mode 100644 index 00000000..88b015f7 --- /dev/null +++ b/inst/scripts/legacy/PointFrame.R @@ -0,0 +1,156 @@ +#' @name PointFrame +#' @title The `PointFrame` class +#' +#' @description +#' The \code{PointFrame} class stores \code{SpatialData} elements from its +#' \code{"points"} layers. These are represented as \code{\link[arrow]{Table}} +#' (\code{data} slot) associated with .zattrs stored as \code{\link{Zattrs}} +#' (\code{meta} slot); a list of \code{metadata} stores other arbitrary info. +#' +#' Currently defined methods (here, \code{x} is a \code{PointFrame}): +#' \itemize{ +#' \item \code{data/meta(x)} to access underlying \code{Table/Zattrs} +#' \item \code{names(x)} returns the underlying table's column names +#' \item \code{dim(x)} returns the dimensions of \code{data(x)} +#' \item \code{`$`,`[[`} directly access columns of \code{data(x)} +#' \item \code{filter,select} to subset rows/columns à la \code{dplyr} +#' \item \code{as.data.frame} to coerce \code{x} to a \code{data.frame} +#' } +#' +#' @param x \code{PointFrame} +#' @param data \code{arrow}-derived table for on-disk, +#' \code{data.frame} for in-memory representation. +#' @param meta \code{\link{Zattrs}} +#' @param metadata optional list of arbitrary +#' content describing the overall object. +#' @param name character string for extraction (see \code{?base::`$`}). +#' @param i,j indices for subsetting (see \code{?base::Extract}). +#' @param drop ignored. +#' @param ... optional arguments passed to and from other methods. +#' +#' @return \code{PointFrame} +#' +#' @examples +#' library(SpatialData.data) +#' zs <- get_demo_SDdata("merfish") +#' x <- file.path(zs, "points", "single_molecule") +#' (p <- readPoint(x)) +#' +#' head(as.data.frame(data(p))) +#' (q <- dplyr::filter(p, cell_type == "VISp_wm")) +#' +#' @importFrom S4Vectors metadata<- +#' @importFrom methods new +#' @export +PointFrame <- function(data=data.frame(), meta=Zattrs(), metadata=list(), ...) { + if(length(meta) < 1){ + meta <- .make_pointshape_meta(data, + encoding_type = "ngff:points", + version = 0.1) + } + x <- .PointFrame(data=data, meta=meta, ...) + metadata(x) <- metadata + return(x) +} + +#' @rdname PointFrame +#' @export +setMethod("names", "PointFrame", \(x) { + setdiff(names(data(x)), "__null_dask_index__") }) + +#' @rdname PointFrame +#' @export +setMethod("dim", "PointFrame", \(x) c(nrow(data(x)), length(names(x)))) + +#' @rdname PointFrame +#' @export +setMethod("length", "PointFrame", \(x) nrow(data(x))) + +#' @rdname PointFrame +#' @importFrom dplyr select all_of collect +#' @exportMethod [[ +setMethod("[[", "PointFrame", \(x, i, ...) { + # x <- select(data(x), !"__null_dask_index__") + x <- select(data(x), names(x)) + collect(select(x, all_of(i)))[[1]] }) + +#' @importFrom utils .DollarNames +#' @export +.DollarNames.PointFrame <- \(x, pattern="") { + setdiff(names(data(x)), "__null_dask_index__") } + +#' @rdname PointFrame +#' @importFrom dplyr select all_of collect +#' @exportMethod $ +setMethod("$", "PointFrame", \(x, name) do.call(`[[`, list(x, name))) + # x <- select(data(x), !"__null_dask_index__") + # collect(select(x, all_of(name)))[[1]] }) + +# sub ---- + +#' @rdname PointFrame +#' @export +setMethod("[", c("PointFrame", "missing", "ANY"), + \(x, i, j, ...) x[seq_len(nrow(x)), j]) + +#' @rdname PointFrame +#' @export +setMethod("[", c("PointFrame", "ANY", "missing"), + \(x, i, j, ...) x[i, seq_len(ncol(x))]) + +#' @rdname PointFrame +#' @export +setMethod("[", c("PointFrame", "missing", "missing"), + \(x, i, j, ...) x[seq_len(nrow(x)), seq_len(ncol(x))]) + +#' @rdname PointFrame +#' @export +setMethod("[", c("PointFrame", "ANY", "character"), \(x, i, j, ...) { + stopifnot(all(j %in% names(x))) + x[i, match(j, names(x))] +}) + +#' @rdname PointFrame +#' @importFrom dplyr mutate filter select +#' @export +setMethod("[", c("PointFrame", "numeric", "numeric"), \(x, i, j, ...) { + if("__null_dask_index__" %in% names(data(x))){ + .i <- `__null_dask_index__` <- NULL # R CMD check + i <- seq_len(nrow(x))[i] + x@data <- data(x) |> + mutate(.i=1+`__null_dask_index__`) |> + filter(.i %in% i) |> + select(-.i) + # make sure this is kept in any case + ndi <- "__null_dask_index__" + ndi <- match(ndi, names(x@data), nomatch=0) + x@data <- x@data[, c(j, ndi)] + } else { + # TODO: can we avoid checking for __null_dask_index__ + x@data <- x@data[i,j,drop = FALSE] + } + return(x) +}) + +#' @rdname PointFrame +#' @importFrom BiocGenerics as.data.frame +#' @export +setMethod("as.data.frame", "PointFrame", \(x) as.data.frame(data(x))[names(x)]) + +setAs( + from="PointFrame", to="data.frame", + function(from) as.data.frame(from)) + +#' @importFrom dplyr filter +#' @export +filter.PointFrame <- \(.data, ...) { + .data@data <- filter(data(.data), ...) + return(.data) +} + +#' @importFrom dplyr select +#' @export +select.PointFrame <- \(.data, ...) { + .data@data <- select(data(.data), ...) + return(.data) +} diff --git a/inst/scripts/legacy/SDattrs_old.R b/inst/scripts/legacy/SDattrs_old.R new file mode 100644 index 00000000..8d606da6 --- /dev/null +++ b/inst/scripts/legacy/SDattrs_old.R @@ -0,0 +1,363 @@ +#' @name SpatialDataAttrs +#' @title The `SpatialDataAttrs` class +#' +#' @param x element or list extracted from a OME-NGFF compliant .zattrs file. +#' @param name character string for extraction (see ?base::`$`). +#' @param type character string; either "array" (image/label) or "frame" (point/shape). +#' @param axes list of axes; if NULL, defaults to cyx (array) or xy (frame). +#' @param transformations list of transformations; if NULL, defaults to global identity. +#' @param ... additional attributes (e.g., version, feature_key). +#' +#' @details +#' When \code{x} is a spatial element, the following applies: +#' \code{SpatialDataFrame}: \code{feature/instance_key}, +#' \code{SingleCellExperiment}: \code{region}, \code{region/instance_key}. +#' +#' When missing \code{x}, \code{SpatialDataAttrs} will generate a valid object +#' with default axes (array: cyx, frame: xy) and transformations (identify) +#' according to the specified type. +#' +#' @return character string +#' +#' @examples +#' x <- file.path("extdata", "blobs.zarr") +#' x <- system.file(x, package="SpatialData") +#' x <- readSpatialData(x) +#' +#' # tables +#' region(table(x)) +#' region_key(table(x)) +#' +#' # points +#' instance_key(point(x)) +#' fk <- feature_key(point(x)) +#' base::table(point(x)[[fk]]) +#' +#' # transformations +#' (z <- meta(label(x))) +#' CTname(z) +#' CTtype(z) +#' CTdata(z, "scale") +#' +#' # constructor +#' SpatialDataAttrs(type="frame") +#' SpatialDataAttrs(type="array") +#' SpatialDataAttrs(type="array", n=7) +#' SpatialDataAttrs(type="array", label=TRUE) +#' +#' @export +SpatialDataAttrs <- \(x, type=c("array", "frame"), + label=FALSE, trans=NULL, ver="0.4", n=3, ...) +{ + 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) + } + # transformations: + ct <- trans %||% .default_ct(ax) + # .zattrs list: + if (type == "array") { + # default structure + res <- list( + omero=list(channels=list(label=letters[seq_len(n)])), + multiscales=list(list( + axes=ax, + version="0.4", + coordinateTransformations=ct, + datasets=list(list(path="0", coordinateTransformations=list(list(type="scale", scale=list(1, 1)))))))) + if (ver == "0.3") res <- list(ome=res) + } else { + # points/shapes + res <- list(axes=ax, coordinateTransformations=ct) + } + res$spatialdata_attrs <- list(version=ver) + SpatialDataAttrs(res) +} + +# 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"))) +} + +# Internal helper to generate coordinate transformations +.default_ct <- \(axes, name="global", type="identity", data=NULL) { + ct <- list(input=axes, output=list(name=name), type=type) + if (!is.null(data)) ct[[type]] <- data + list(ct) +} + +#' @export +#' @importFrom utils .DollarNames +.DollarNames.SpatialDataAttrs <- \(x, pattern="") names(x) + +#' @rdname SpatialDataAttrs +#' @exportMethod $ +setMethod("$", "SpatialDataAttrs", \(x, name) x[[name]]) + +# internal use only! +#' @noRd +.zv <- \(x) { + v <- x$spatialdata_attrs$version + if (!length(v)) stop("couldn't find 'version' in 'spatialdata_attrs'") + ok <- length(v) == 1 && is.character(v) && v %in% sprintf("0.%d", seq_len(5)) + if (!ok) stop("invalid 'version' in 'spatialdata_attrs'; expected '0.x' where x is 1-5") + return(v) +} + +# internal use only! +#' @noRd +.ms <- \(x) switch(.zv(x), "0.3"=x$ome$multiscales, x$multiscales) + +# internal use only! +#' @noRd +.ch <- \(x) { + if (.zv(x) == "0.3") x <- x$ome + unlist(x$omero$channels) +} + +# internal use only! +#' @noRd +setMethod("multiscales", "list", .ms) + +#' @export +setMethod("channels", "SpatialDataAttrs", \(x, ...) .ch(x)) + +# features ---- + +#' @export +#' @rdname SpatialDataAttrs +setMethod("feature_key", "SpatialDataPoint", \(x) feature_key(meta(x))) +#' @export +#' @rdname SpatialDataAttrs +setMethod("feature_key", "SpatialDataAttrs", \(x) x$spatialdata_attrs$feature_key) +#' @export +#' @rdname SpatialDataAttrs +setReplaceMethod("feature_key", c("SpatialDataAttrs", "character"), + \(x, value) { x$spatialdata_attrs$feature_key <- value; x }) + +# region(s) ---- + +#' @export +#' @rdname SpatialDataAttrs +setMethod("region_key", "SingleCellExperiment", \(x) meta(x)$region_key) + +# internal use only! +#' @noRd +#' @importFrom SingleCellExperiment int_metadata<- +setReplaceMethod("region_key", c("SingleCellExperiment", "character"), \(x, value) { + stopifnot(length(value) == 1, nchar(value) > 0) + int_metadata(x)$spatialdata_attrs$region_key <- value + return(x) +}) + +# internal use only! +#' @noRd +#' @importFrom SingleCellExperiment int_metadata<- +setReplaceMethod("region_key", c("SingleCellExperiment", "NULL"), \(x, value) { + int_metadata(x)$spatialdata_attrs$region_key <- value + return(x) +}) + +#' @export +#' @rdname SpatialDataAttrs +setMethod("region", "SingleCellExperiment", \(x) { + rk <- region_key(x) + if (is.null(rk)) return(NULL) + meta(x)[[rk]] +}) + +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom SingleCellExperiment int_colData +setMethod("regions", "SingleCellExperiment", \(x) { + rk <- region_key(x) + if (is.null(rk)) return(NULL) + int_colData(x)[[rk]] +}) + +# internal use only! +#' @noRd +#' @importFrom SingleCellExperiment int_metadata<- +setReplaceMethod("region", c("SingleCellExperiment", "character"), \(x, value) { + stopifnot(all(nchar(value) > 0, na.rm=TRUE)) + if (is.null(rk <- region_key(x))) + rk <- region_key(x) <- "region" + int_metadata(x)$spatialdata_attrs[[rk]] <- sort(unique(value)) + return(x) +}) + +# internal use only! +#' @noRd +#' @importFrom SingleCellExperiment int_metadata<- +setReplaceMethod("region", c("SingleCellExperiment", "NULL"), \(x, value) { + if (!is.null(rk <- region_key(x))) + int_metadata(x)$spatialdata_attrs[[rk]] <- value + return(x) +}) + +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom SingleCellExperiment int_colData<- +setReplaceMethod("regions", c("SingleCellExperiment", "character"), \(x, value) { + stopifnot(length(value) %in% c(1, ncol(x))) + stopifnot(all(nchar(value) > 0, na.rm=TRUE)) + if (is.null(rk <- region_key(x))) region_key(x) <- "region" + int_metadata(x)$spatialdata_attrs[[rk]] <- sort(unique(value)) + int_colData(x)[[rk]] <- value + return(x) +}) + +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom SingleCellExperiment int_colData<- +setReplaceMethod("regions", c("SingleCellExperiment", "NULL"), \(x, value) { + if (!is.null(rk <- region_key(x))) { + int_metadata(x)$spatialdata_attrs[[rk]] <- value + int_colData(x)[[rk]] <- value + } + region_key(x) <- value + return(x) +}) + +# instances ---- + +# NOTE: does not apply to images +#' @export +#' @rdname SpatialDataAttrs +setMethod("instance_key", "list", \(x) x$instance_key) +#' @export +#' @rdname SpatialDataAttrs +setMethod("instance_key", "SingleCellExperiment", \(x) instance_key(meta(x))) +#' @export +#' @rdname SpatialDataAttrs +setMethod("instance_key", "SpatialDataFrame", \(x) instance_key(meta(x)$spatialdata_attrs)) +#' @export +#' @rdname SpatialDataAttrs +setMethod("instance_key", "SpatialDataLabel", \(x) instance_key(meta(x)$spatialdata_attrs)) +#' @export +#' @rdname SpatialDataAttrs +setReplaceMethod("instance_key", c("SpatialDataAttrs", "character"), \(x, value) { + x$spatialdata_attrs$instance_key <- value + return(x) +}) +#' @export +#' @rdname SpatialDataAttrs +setReplaceMethod("instance_key", c("SingleCellExperiment", "character"), \(x, value) { + int_metadata(x)$spatialdata_attrs$instance_key <- value + return(x) +}) + +#' @export +#' @rdname SpatialDataAttrs +setMethod("instances", "SpatialDataLabel", \(x) { + # unique values in first scale, excluding 0 + z <- data(x, 1) + as.integer(setdiff(unique(as.vector(z)), 0)) +}) +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom dplyr pull +setMethod("instances", "SpatialDataPoint", \(x) pull(data(x), instance_key(x))) +#' @export +#' @rdname SpatialDataAttrs +setMethod("instances", "SpatialDataShape", \(x) { + ik <- tryCatch(instance_key(x), error=\(e) NULL) + if (is.null(ik)) return(seq_len(nrow(x))) + pull(data(x), ik) +}) +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom SingleCellExperiment int_colData +setMethod("instances", "SingleCellExperiment", \(x) { + if (is.null(ik <- instance_key(x))) + stop("no 'instance_key' found in 'x'") + int_colData(x)[[ik]] +}) + +#' @export +#' @rdname SpatialDataAttrs +#' @importFrom SingleCellExperiment int_colData<- +setReplaceMethod("instances", c("SingleCellExperiment", "ANY"), \(x, value) { + ik <- instance_key(x) + if (is.null(ik)) + ik <- "instance_id" + int_colData(x)[[ik]] <- value + return(x) +}) + +# elements ---- + +setMethod("version", c("SpatialDataElement"), \(x) { + version(meta(x)) +}) + +setMethod("version", c("SingleCellExperiment"), \(x) { + meta(x)$version +}) + +setMethod("version", "Zattrs", \(x) .zv(x)) + +setMethod("version", "list", \(x) .zv(x)) + +.zv <- \(x) { + v <- x$spatialdata_attrs$version + if (!length(v)) stop("couldn't find 'version' in 'spatialdata_attrs'") + ok <- length(v) == 1 && is.character(v) && v %in% sprintf("0.%d", seq_len(5)) + if (!ok) stop("invalid 'version' in 'spatialdata_attrs'; expected '0.x' where x is 1-5") + return(v) +} + +setReplaceMethod("version", c("sdFrame"), \(x, value) { + if(!value %in% c("0.1", "0.2", "0.3")) + stop("Unknown version for shape/point! Must be 0.2 or 0.3.") + meta(x)$spatialdata_attrs$version <- value + x +}) + +setReplaceMethod("version", c("sdArray"), \(x, value) { + mt <- meta(x) + if(value == "0.3"){ + if(is.null(mt$ome)){ + mt$ome = list(omero = mt$omero, + multiscales = mt$multiscales) + mt$omero <- NULL + mt$multiscales <- NULL + } + } else if(value %in% c("0.1" ,"0.2")){ + if(is.null(mt$multiscales)){ + mt$omero <- mt$ome$omero + mt$multiscales <- mt$ome$multiscales + mt[["ome"]] <- NULL + } + } else { + stop("Unknown version for image/label! Must be 0.1, 0.2, 0.3.") + } + mt$spatialdata_attrs$version <- value + meta(x) <- mt + x +}) + +setReplaceMethod("version", c("SingleCellExperiment"), \(x, value) { + if(!value %in% c("0.1", "0.2")) + stop("Unknown version for table! Must be 0.1 or 0.2") + int_metadata(x)$spatialdata_attrs$version <- value + return(x) +}) \ No newline at end of file diff --git a/inst/scripts/legacy/ShapeFrame.R b/inst/scripts/legacy/ShapeFrame.R new file mode 100644 index 00000000..f6795436 --- /dev/null +++ b/inst/scripts/legacy/ShapeFrame.R @@ -0,0 +1,103 @@ +#' @name ShapeFrame +#' @title The `ShapeFrame` class +#' @aliases geom_type +#' +#' @param x \code{ShapeFrame} +#' @param data \code{arrow}-derived table for on-disk, +#' \code{data.frame} for in-memory representation. +#' @param meta \code{\link{Zattrs}} +#' @param metadata optional list of arbitrary +#' content describing the overall object. +#' @param name character string for extraction (see \code{?base::`$`}). +#' @param i,j indices specifying elements to extract. +#' @param drop,pattern ignored. +#' @param ... optional arguments passed to and from other methods. +#' +#' @return \code{ShapeFrame} +#' +#' @examples +#' library(SpatialData.data) +#' zs <- get_demo_SDdata("merfish") +#' +#' y <- file.path(zs, "shapes", "cells") +#' (s <- readShape(y)) +#' plot(sf::st_as_sf(data(s)), cex=0.2) +#' +#' y <- file.path(zs, "shapes", "anatomical") +#' (s <- readShape(y)) +#' plot(sf::st_as_sf(data(s)), cex=0.2) +#' +#' @importFrom S4Vectors metadata<- +#' @importFrom methods new +#' @export +ShapeFrame <- function(data=data.frame(), meta=Zattrs(), metadata=list(), ...) { + if(length(meta) < 1){ + meta <- .make_pointshape_meta(data, + encoding_type = "ngff:points", + version = 0.1) + } + x <- .ShapeFrame(data=data, meta=meta, ...) + metadata(x) <- metadata + return(x) +} + +# TODO: it's really annoying that this doesn't just inherit +# data.frame() operations, cuz data are in an extra slot... +# but else not sure how to assure validity, stash .zattrs etc. + +#' @rdname ShapeFrame +#' @export +setMethod("dim", "ShapeFrame", \(x) dim(data(x))) + +#' @rdname ShapeFrame +#' @export +setMethod("length", "ShapeFrame", \(x) nrow(data(x))) + +#' @rdname ShapeFrame +#' @export +setMethod("names", "ShapeFrame", \(x) names(data(x))) + +#' @export +#' @rdname ShapeFrame +#' @importFrom utils .DollarNames +.DollarNames.ShapeFrame <- \(x, pattern="") + grep(pattern, names(x), value=TRUE) + +#' @rdname ShapeFrame +#' @exportMethod $ +setMethod("$", "ShapeFrame", \(x, name) data(x)[[name]]) + +#' @export +#' @rdname ShapeFrame +#' @importFrom sf st_as_sf st_geometry_type +setMethod("geom_type", "ShapeFrame", \(x) { + y <- st_as_sf(data(x[1, ])) + z <- st_geometry_type(y) + return(as.character(z)) +}) + +# sub ---- + +#' @rdname ShapeFrame +#' @export +setMethod("[", c("ShapeFrame", "missing", "ANY"), + \(x, i, j, ...) x[seq_len(nrow(x)), j]) + +#' @rdname ShapeFrame +#' @export +setMethod("[", c("ShapeFrame", "ANY", "missing"), + \(x, i, j, ...) x[i, seq_len(ncol(x))]) + +#' @rdname ShapeFrame +#' @export +setMethod("[", c("ShapeFrame", "missing", "missing"), + \(x, i, j, ...) x[seq_len(nrow(x)), seq_len(ncol(x))]) + +#' @rdname ShapeFrame +#' @export +setMethod("[", c("ShapeFrame", "numeric", "numeric"), \(x, i, j, ...) { + i <- seq_len(nrow(x))[i] + j <- seq_len(ncol(x))[j] + x@data <- x@data[i, j] + return(x) +}) diff --git a/inst/scripts/legacy/Zattrs.R b/inst/scripts/legacy/Zattrs.R new file mode 100644 index 00000000..60a6742c --- /dev/null +++ b/inst/scripts/legacy/Zattrs.R @@ -0,0 +1,160 @@ +#' @name Zattrs +#' @title The `Zattrs` class +#' +#' @param x list extracted from a OME-NGFF compliant .zattrs file. +#' @param name character string for extraction (see ?base::`$`). +#' @param type character string; either "array" (image/label) or "frame" (point/shape). +#' @param axes list of axes; if NULL, defaults to cyx (array) or xy (frame). +#' @param transformations list of transformations; if NULL, defaults to global identity. +#' @param ... additional attributes (e.g., version, feature_key). +#' +#' @details +#' When missing \code{x}, \code{Zattrs} will generate a valid object with +#' default axes (array: cyx, frame: xy) and transformations (identify) +#' according to the specified type. +#' +#' @return \code{Zattrs} +#' +#' @examples +#' x <- file.path("extdata", "blobs.zarr") +#' x <- system.file(x, package="SpatialData") +#' x <- readSpatialData(x, tables=FALSE) +#' +#' (z <- meta(label(x))) +#' +#' CTname(z) +#' CTtype(z) +#' CTdata(z, "scale") +#' +#' feature_key(point(x)) +#' +#' # constructor +#' Zattrs(type="frame") +#' Zattrs(type="array") +#' Zattrs(type="array", n=7) +#' Zattrs(type="array", label=TRUE) +#' +#' @export +Zattrs <- function(x, type=c("array", "frame"), label=FALSE, trans=NULL, + ver="0.3", n=3, scale_factors = NULL, ...) { + if (!missing(x)) return(.Zattrs(x)) + type <- match.arg(type) + # axes: + ax <- .default_ax(type, label) + # transformations: + ct <- trans %||% .default_ct(ax) + ds <- .default_ds(.ax_names(ax), scale_factors) + # .zattrs list: + if (type == "array") { + # default structure + res <- list() + if(!label) + res <- c(res, + list(omero=list(channels=lapply(letters[seq_len(n)], + \(.) list(label = .))))) + res <- c(res, + list( + multiscales= + list( + list( + axes=ax, + version="0.4", + coordinateTransformations=ct, + datasets=ds + ) + ) + ) + ) + if (ver == "0.3") res <- list(ome=res) + } else { + # points/shapes + res <- list(axes=ax, coordinateTransformations=ct) + } + res$spatialdata_attrs <- list(version=ver) + Zattrs(res) +} + +# Internal helper to generate OME-NGFF axes +.default_ax <- \(type=c("array", "frame"), label = FALSE) { + switch(match.arg(type), + # (c)yx for images/labels + array={ + ax <- list( + list(name="x", type="space"), + list(name="y", type="space")) + if (type == "array") { + # yx for labels + ax <- rev(ax) + # yx for images, cyx if requested + if (!label) ax <- c(list(list(name="c", type="channel")), ax) + } + ax + }, + # xy for points/shapes + list("x", "y")) +} + +.ax_names <- function(ax){ + if (is.character(ax[[1]])) { + unlist(ax) + } else { + vapply(ax, \(.) .$name, character(1)) + } +} + +# Internal helper to generate coordinate transformations +.default_ct <- \(axes, name="global", type="identity", data=NULL) { + ct <- list( + input=list(axes = axes, + name = if(length(axes) == 3) "cyx" else "yx"), + output=list(axes = axes, + name = name), + type = type) + if (!is.null(data)) ct[[type]] <- data + list(ct) +} + +.default_ds <- function(axes, scale_factors = NULL){ + scale_factors <- cumprod(c(1,scale_factors)) + paths <- paste0(seq_along(scale_factors) - 1) + mapply(\(p,s) { + list( + coordinateTransformations = list( + list( + scale = lapply( + axes, + \(.) if(. == "c") 1 else s), + type = "scale" + ) + ), + path = p + ) + }, paths, scale_factors, USE.NAMES = FALSE, SIMPLIFY = FALSE) +} + + +#' @export +#' @importFrom utils .DollarNames +.DollarNames.Zattrs <- \(x, pattern="") names(x) + +#' @rdname Zattrs +#' @exportMethod $ +setMethod("$", "Zattrs", \(x, name) x[[name]]) + +# internal use only! +#' @noRd +.ms <- \(x) switch(version(x), "0.3"=x$ome$multiscales, x$multiscales) + +# internal use only! +#' @noRd +.ch <- \(x) { + if (version(x) == "0.3") x <- x$ome + unlist(x$omero$channels) +} + +# internal use only! +#' @noRd +setMethod("multiscales", "list", .ms) + +#' @export +setMethod("channels", "Zattrs", \(x, ...) .ch(x)) diff --git a/inst/scripts/legacy/metadata.R b/inst/scripts/legacy/metadata.R new file mode 100644 index 00000000..dbde8ee6 --- /dev/null +++ b/inst/scripts/legacy/metadata.R @@ -0,0 +1,285 @@ +# helpers ---- + +.make_axes_meta <- function(x, unit = FALSE){ + lapply(x, \(.){ + meta <- list(names = ., + type = if(. == "c") "channel" else "space") + if(unit) + meta <- c(meta, list(unit = "unit")) + meta + }) +} + +.make_empty_ct <- function(x){ + space <- .make_axes_meta(x, unit = TRUE) + input <- list(axes = space, + name = paste(x, collapse = "")) + output <- list(axes = space, name = "global") + meta <- list( + list(input = input, + output = output, + type = "identity") + ) + meta +} + +# .make_datasets <- function(x, axes){ +# paths <- paste0(seq_len(length(x)) - 1) +# mapply(\(p) { +# list( +# coordinateTransformations = list( +# list( +# scale = vapply(axes, \(.){ +# if(. == "c") 1 else (2^as.numeric(p)) +# }, numeric(1)), +# type = "scale" +# ) +# ), +# path = p +# ) +# }, paths, USE.NAMES = FALSE, SIMPLIFY = FALSE) +# } + +# metadata constructors ---- +#' @title Make point/shape metadata +#' @description Make point/shape metadata +#' @param x A points or shapes object +#' @param axes A character vector of axes names +#' @param encoding_type A string specifying the encoding type +#' @param feature_key A string specifying the feature key +#' @param instance_key A string specifying the instance key +#' @param version A string specifying the version +#' @return A list of metadata for the point/shape object +#' @importFrom jsonlite fromJSON toJSON +#' @noRd +.make_pointshape_meta <- function(x, + axes = NULL, + encoding_type = "ngff:points", + feature_key = NULL, + instance_key = NULL, + version = 0.1){ + meta <- list() + ax <- "axes" + ct <- "coordinateTransformations" + sa <- "spatialdata_attrs" + + # axis + # NOTE: rev dimensions since points and shapes want x, y + # whereas images and labels want y, x, etc. + meta[[ax]] <- rev(.get_valid_axes(x, axes = axes, image = FALSE)) + + # encoding type + meta[["encoding-type"]] <- encoding_type + + # spatialdata_attrs + meta[[sa]] <- list(version = version) + if(!is.null(feature_key)) + meta[[sa]][["feature_key"]] <- feature_key + if(!is.null(instance_key)) + meta[[sa]][["instance_key"]] <- instance_key + + # coordinate transformations + meta[[ct]] <- .make_empty_ct(meta[[ax]]) + + # update json list + meta <- fromJSON(toJSON(meta, auto_unbox = TRUE), simplifyVector = FALSE) + Zattrs(meta) +} +# TODO: make it the functions take a global option e.g. sd_zarr_version +# as an argument for the default zarr version + +#' @title Make image metadata +#' @description Make image metadata +#' @param x An image object +#' @param axes A character vector of axes names +#' @param version A string specifying the version +#' @return A list of metadata for the image object +#' @importFrom jsonlite fromJSON toJSON +#' @noRd +.make_image_meta <- function(x, + axes = NULL, + version = 0.4){ + meta <- list() + ax <- "axes" + ct <- "coordinateTransformations" + ds <- "datasets" + mt <- "metadata" + v <- "version" + n <- "name" + + # axis + axes <- .get_valid_axes(x, axes, image = TRUE) + meta[[ax]] <- .make_axes_meta(axes, unit = FALSE) + + # coordinate transformations + # TODO: shall we do coordinate transformations only + # without datasets:coordinateTransformations + # see https://ngff.openmicroscopy.org/0.4/index.html#multiscale-md + meta[[ct]] <- .make_empty_ct(axes) + + # datasets + meta[[ds]] <- .make_datasets(x, axes) + + # metadata + meta[[mt]] <- list(omero = list( + channels = lapply(seq_len(length(axes))-1, \(.) + list(label = .)) + )) + + # name + meta[[n]] <- "" + + # version + meta[[v]] <- list(version = version) + + # multiscales + meta <- list(multiscales = list(meta), + omero = list( + channels = lapply(seq_len(length(axes))-1, \(.) + list(label = .)) + ), + spatialdata_attrs = list(version = "0.1")) + + # update json list + # meta <- fromJSON(toJSON(meta, auto_unbox = TRUE), simplifyVector = FALSE) + Zattrs(meta) +} + +#' @title Make label metadata +#' @description Make label metadata +#' @param x A label object +#' @param axes A character vector of axes names +#' @param version A string specifying the version +#' @return A list of metadata for the label object +#' @importFrom jsonlite fromJSON toJSON +#' @noRd +.make_label_meta <- function(x, + axes = NULL, + version = 0.4){ + meta <- list() + ax <- "axes" + ct <- "coordinateTransformations" + ds <- "datasets" + v <- "version" + n <- "name" + + # axis + axes <- .get_valid_axes(x, axes, image = FALSE) + if(is.null(axes)){ + axes <- c("y", "x") + axes <- if(length(dim(x)) == 3) c("z", axes) else axes + } + meta[[ax]] <- .make_axes_meta(axes, unit = FALSE) + + # coordinate transformations + # TODO: shall we do coordinate transformations only + # without datasets:coordinateTransformations + # see https://ngff.openmicroscopy.org/0.4/index.html#multiscale-md + meta[[ct]] <- .make_empty_ct(axes) + + # datasets + meta[[ds]] <- .make_datasets(x, axes) + + # name + meta[[n]] <- "" + + # version + meta[[v]] <- list(version = version) + + # multiscales + meta <- list(`image-label`= list(version = version), + multiscales = list(meta), + spatialdata_attrs = list(version = "0.1")) + + # update json list + meta <- fromJSON(toJSON(meta, auto_unbox = TRUE), simplifyVector = FALSE) + Zattrs(meta) +} + +#' @title Get valid axes +#' @description Get validated axes +#' +#' @inheritParams write_image +#' +#' @noRd +.get_valid_axes <- function( + x, + axes = NULL, + image = FALSE +) { + + # axes may be string e.g. "tczyx" + if (is.character(axes) && length(axes) == 1L) + axes <- strsplit(axes, "", fixed = TRUE)[[1]] + + # We can guess axes for images, labels, points/shapes + ndim <- length(.get_dim(x)) + if (is.null(axes)) { + if (ndim == 2) { + axes <- c("y", "x") + } else { + if(image){ + stop("axes must be provided. Can't be guessed beyond 2D image", + call. = FALSE) + } else { + if(ndim == 3) { + axes <- c("z", "y", "x") + } else { + stop("axes must be provided. Can't be guessed beyond 2D or 3D data", + call. = FALSE) + } + } + } + } else { + if (length(axes) != ndim) { + stop( + sprintf("axes length (%d) must match number of dimensions (%d)", + length(axes), ndim), + call. = FALSE + ) + } + } + + axes +} + +# TODO: what is the best way to get the inherint dimension of geometry +# objects + +.get_dim <- function(x){ + if(is.list(x) && length(x) > 0 && + !is.matrix(x) && !is.data.frame(x)) + x <- x[[1]] + if("arrow_OR_df" %in% is(x)){ + return(.get_arrow_dim(x)) + } else if(!is.null(nd <- dim(x))) { + return(nd) + } else { + # TODO: I don't like this message! + stop("no dimensions!") + } +} + +#' @importFrom sf st_as_sf st_geometry +.get_arrow_dim <- function(x){ + if("geometry" %in% colnames(x)){ + sfx <- st_as_sf(x) + sfx <- st_geometry(sfx) + axes <- class(st_geometry(st_as_sf(sfx))[[1]]) + if("XY" %in% axes){ + n_col <- 2 + } else if("XYZ" %in% axes){ + n_col <- 3 + } else{ + stop("No geometry object is detected!") + } + } else { + axes <- colnames(x) + axes <- axes[axes %in% c("x", "y", "z")] + n_col <- length(axes) + } + return(c(nrow(x), n_col)) +} + + + diff --git a/inst/scripts/legacy/sdImage.R b/inst/scripts/legacy/sdImage.R new file mode 100644 index 00000000..7d6eb184 --- /dev/null +++ b/inst/scripts/legacy/sdImage.R @@ -0,0 +1,131 @@ +#' @name SpatialDataImage +#' @title The `SpatialDataImage` class +#' @aliases channels +#' +#' @param x \code{SpatialDataImage} +#' @param data list of \code{\link{SpatialDataAttrs}}s +#' @param meta \code{\link{SpatialDataAttrs}} +#' @param metadata optional list of arbitrary +#' content describing the overall object. +#' @param multiscale if TRUE (and \code{data} is not a list), +#' multiscale image will be generated. +#' @param axes axes +#' @param i,j indices specifying elements to extract. +#' @param k scalar index specifying which scale to extract. +#' @param drop ignored. +#' @param ... option arguments passed to and from other methods. +#' +#' @return \code{SpatialDataImage} +#' +#' @examples +#' zs <- file.path("extdata", "blobs.zarr") +#' zs <- system.file(zs, package="SpatialData") +#' +#' pa <- list.dirs( +#' file.path(zs, "images"), +#' recursive=FALSE, full.names=TRUE) +#' +#' # simple +#' readImage(pa[1]) +#' +#' # multi-scale +#' (x <- readImage(pa[2])) +#' +#' dim(data(x, 1)) # highest res. +#' dim(data(x, Inf)) # lowest res. +#' +#' rgb <- apply( +#' data(x, 1), c(2, 3), +#' \(.) rgb(.[1], .[2], .[3])) +#' plot( +#' row(rgb), col(rgb), col=rgb, +#' pch=15, asp=1, ylim=c(ncol(rgb), 0)) +#' +#' @importFrom S4Vectors metadata<- +#' @importFrom methods new +#' @importFrom DelayedArray DelayedArray +#' @export +SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(), + version = image(sdFormat(0.1)), + metadata=list(), + scale_factors = NULL, ...) { + if(!is.list(data)) + data <- list(data) + if(!is.null(scale_factors)){ + data <- .generate_multiscale(data[[1]], + axes = vapply(axes(meta), + \(.) .$name, + character(1)), + scale_factors = scale_factors, + method = "image") + # TODO: this supposed to update the scale_factors not write a new meta + meta <- SpatialDataAttrs(scale_factors = scale_factors) + } + # construct S4 object + x <- .SpatialDataImage(data=data, meta=meta, ...) + metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version + return(x) +} + +#' @export +#' @rdname SpatialDataImage +setMethod("channels", "SpatialDataImage", \(x, ...) channels(meta(x))) + +#' @export +#' @rdname SpatialDataImage +setMethod("channels", "ANY", \(x, ...) stop("only 'images' have channels")) + +#' @importFrom S4Vectors isSequence +.get_multiscales_paths <- function(x) { + ps <- list.files(x) + ps <- suppressWarnings(as.numeric(sort(ps, decreasing=FALSE))) + ps <- ps[!is.na(ps)] + if (length(ps)) { + qs <- seq(min(ps), max(ps)) + if (!isTRUE(all.equal(ps, qs))) + stop("SpatialDataImage paths are ill-defined, should", + " be an integer sequence, e.g., 0,1,...,n") + } else { + stop("SpatialDataImage path is empty") + } + return(ps) +} + +.check_jk <- \(x, .) { + if (isTRUE(x)) return() + tryCatch( + stopifnot( + is.numeric(x), x == round(x), + diff(range(x)) == length(x)-1, + (y <- abs(x)) == seq(min(y), max(y)) + ), + error=\(e) stop(sprintf("invalid '%s'", .)) + ) +} + +#' @exportMethod [ +#' @rdname SpatialDataImage +#' @importFrom utils head tail +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 +}) \ No newline at end of file diff --git a/inst/scripts/legacy/sdLabel.R b/inst/scripts/legacy/sdLabel.R new file mode 100644 index 00000000..9375f5ee --- /dev/null +++ b/inst/scripts/legacy/sdLabel.R @@ -0,0 +1,86 @@ +#' @name SpatialDataLabel +#' @title The \code{SpatialDataLabel} class +#' +#' @description +#' The \code{SpatialDataLabel} class stores \code{SpatialData} elements from its +#' \code{"labels"} layers. These are represented as a \code{ZarrMatrix} +#' (\code{data} slot) associated with \code{\link{SpatialDataAttrs}} +#' (\code{meta} slot); a list of \code{metadata} stores other arbitrary info. +#' +#' Currently defined methods (here, \code{x} is a \code{SpatialDataLabel}): +#' \itemize{ +#' \item \code{data/meta(x)} to access underlying \code{ZarrMatrix/SpatialDataAttrs} +#' \item \code{dim(x)} returns the dimensions of \code{data(x)} +#' } +#' +#' @param x \code{SpatialDataLabel} +#' @param data list of \code{\link[ZarrArray]{ZarrArray}}s +#' @param meta \code{\link{SpatialDataAttrs}} +#' @param metadata optional list of arbitrary +#' content describing the overall object. +#' @param multiscale if TRUE (and \code{data} is not a list), +#' multiscale image will be generated. +#' @param axes axes +#' @param i,j indices specifying elements to extract. +#' @param drop ignored. +#' @param ... option arguments passed to and from other methods. +#' +#' @return \code{SpatialDataLabel} +#' +#' @examples +#' x <- file.path("extdata", "blobs.zarr") +#' x <- system.file(x, package="SpatialData") +#' x <- file.path(x, "labels", "blobs_labels") +#' +#' (y <- readLabel(x)) +#' y[1:10, 1:10] +#' meta(y) +#' +#' @importFrom S4Vectors metadata<- +#' @importFrom methods new +#' @export +SpatialDataLabel <- function(data=list(), + meta=SpatialDataAttrs(label = TRUE), + version = image(sdFormat(0.1)), + metadata=list(), + scale_factors = NULL, ...) { + if(!is.list(data)) + data <- list(data) + if(!is.null(scale_factors)){ + data <- .generate_multiscale(data[[1]], + axes = vapply(axes(meta), + \(.) .$name, + character(1)), + scale_factors = scale_factors, + method = "label") + meta <- SpatialDataAttrs(scale_factors = scale_factors, label = TRUE) + } + x <- .SpatialDataLabel(data=data, meta=meta, ...) + metadata(x) <- metadata + + # update version if provided + if(!is.null(version)) + version(x) <- version + return(x) +} + +#' @rdname SpatialDataLabel +#' @importFrom utils head tail +#' @exportMethod [ +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 +}) \ No newline at end of file diff --git a/inst/scripts/legacy/test-imagearray.R b/inst/scripts/legacy/test-imagearray.R new file mode 100644 index 00000000..bcf25c41 --- /dev/null +++ b/inst/scripts/legacy/test-imagearray.R @@ -0,0 +1,36 @@ +require(DelayedArray, quietly = TRUE) +rgb <- seq_len(255) + +test_that("ImageArray()", { + val <- sample(rgb, 3*20*20, replace=TRUE) + mat <- array(val, dim=c(3, 20, 20)) + # invalid + imgarray <- ImageArray(mat) + expect_equal(dim(imgarray), dim(mat)) + expect_error(ImageArray(mat, 1)) + expect_error(ImageArray(mat, list())) + # single scale + expect_silent(ImageArray(list())) + expect_silent(ImageArray(list(mat))) + expect_silent(ImageArray(list(mat), Zattrs())) + + # multiscale + # only for ImageArray with 2 dimensions, we can guess the dimensions + dim <- lapply(c(20, 10, 5), \(.) c(3, rep(., 2))) + lys <- lapply(dim, \(.) array(sample(rgb, prod(.), replace=TRUE), dim=.)) + expect_silent(ImageArray(lys)) +}) + +test_that("data(),ImageArray", { + dim <- lapply(c(8, 4, 2), \(.) c(3, rep(., 2))) + lys <- lapply(dim, \(.) array(0, dim=.)) + img <- ImageArray(lys) + for (. in seq_along(lys)) + expect_identical(data(img, .), lys[[.]]) + expect_identical(data(img, Inf), lys[[3]]) + expect_error(data(img, 0)) + expect_error(data(img, -1)) + expect_error(data(img, 99)) + expect_error(data(img, "")) + expect_error(data(img, c(1,2))) +}) \ No newline at end of file diff --git a/inst/scripts/legacy/test-labelarray.R b/inst/scripts/legacy/test-labelarray.R new file mode 100644 index 00000000..337a5ceb --- /dev/null +++ b/inst/scripts/legacy/test-labelarray.R @@ -0,0 +1,32 @@ +test_that("LabelArray()", { + val <- sample(seq_len(12), 20*20, replace=TRUE) + mat <- array(val, dim=c(20, 20)) + lblarray <- LabelArray(mat) + expect_equal(dim(lblarray), dim(mat)) + # invalid + expect_error(LabelArray(mat, 1)) + expect_error(LabelArray(mat, list())) + # single scale + expect_silent(LabelArray(list())) + expect_silent(LabelArray(list(mat))) + expect_silent(LabelArray(list(mat), Zattrs())) + # multiscale + dim <- lapply(c(20, 10, 5), \(.) rep(., 2)) + lys <- lapply(dim, \(.) array(sample(seq_len(12), prod(.), replace=TRUE), dim=.)) + expect_silent(LabelArray(lys)) +}) + + +test_that("data(),LabelArray", { + dim <- lapply(c(8, 4, 2), \(.) rep(., 2)) + lys <- lapply(dim, \(.) array(0L, dim=.)) + lab <- LabelArray(lys) + for (. in seq_along(lys)) + expect_identical(data(lab, .), lys[[.]]) + expect_identical(data(lab, Inf), lys[[3]]) + expect_error(data(lab, 0)) + expect_error(data(lab, -1)) + expect_error(data(lab, 99)) + expect_error(data(lab, "")) + expect_error(data(lab, c(1,2))) +}) \ No newline at end of file diff --git a/man/SpatialData.Rd b/man/SpatialData.Rd index a9069374..9b7c6ed1 100644 --- a/man/SpatialData.Rd +++ b/man/SpatialData.Rd @@ -60,8 +60,6 @@ \alias{points,SpatialData-method} \alias{shapes,SpatialData-method} \alias{tables,SpatialData-method} -\alias{[[<-,SpatialData,numeric,ANY,ANY-method} -\alias{[[<-,SpatialData,character,ANY,ANY-method} \title{The `SpatialData` class} \usage{ SpatialData(images, labels, points, shapes, tables) @@ -104,9 +102,9 @@ SpatialData(images, labels, points, shapes, tables) \S4method{tables}{SpatialData}(x) -\S4method{[[}{SpatialData,numeric,ANY,ANY}(x, i) <- value +\S4method{[[}{SpatialData,numeric,ANY}(x, i) <- value -\S4method{[[}{SpatialData,character,ANY,ANY}(x, i) <- value +\S4method{[[}{SpatialData,character,ANY}(x, i) <- value } \arguments{ \item{images}{list of \code{\link{SpatialDataImage}}s} diff --git a/man/SpatialDataArray.Rd b/man/SpatialDataArray.Rd index a50874af..40fb70d8 100644 --- a/man/SpatialDataArray.Rd +++ b/man/SpatialDataArray.Rd @@ -21,14 +21,18 @@ SpatialDataImage( data = list(), meta = SpatialDataAttrs(), + version = image(sdFormat(0.1)), metadata = list(), + scale_factors = NULL, ... ) SpatialDataLabel( data = list(), - meta = SpatialDataAttrs(), + meta = SpatialDataAttrs(label = TRUE), + version = image(sdFormat(0.1)), metadata = list(), + scale_factors = NULL, ... ) diff --git a/man/SpatialDataAttrs.Rd b/man/SpatialDataAttrs.Rd index e2d2029c..9ebe8593 100644 --- a/man/SpatialDataAttrs.Rd +++ b/man/SpatialDataAttrs.Rd @@ -43,6 +43,7 @@ SpatialDataAttrs( trans = NULL, ver = "0.4", nch = 3, + scale_factors = NULL, ... ) diff --git a/man/SpatialDataFrame.Rd b/man/SpatialDataFrame.Rd index c175fb8d..98393a03 100644 --- a/man/SpatialDataFrame.Rd +++ b/man/SpatialDataFrame.Rd @@ -24,6 +24,7 @@ SpatialDataPoint( data = NULL, meta = SpatialDataAttrs(type = "frame"), + version = point(sdFormat(0.1)), metadata = list(), ik = NULL, fk = NULL, @@ -33,6 +34,7 @@ SpatialDataPoint( SpatialDataShape( data = NULL, meta = SpatialDataAttrs(type = "frame"), + version = shape(sdFormat(0.1)), metadata = list(), ... ) diff --git a/man/create_zarr.Rd b/man/create_zarr.Rd new file mode 100644 index 00000000..755810e7 --- /dev/null +++ b/man/create_zarr.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/zarr_utils.R +\name{create_zarr} +\alias{create_zarr} +\title{create_zarr} +\usage{ +create_zarr(store, version = 2) +} +\arguments{ +\item{store}{The location of the Zarr store} + +\item{version}{Zarr version} +} +\value{ +`NULL` +} +\description{ +Create Zarr store +} +\examples{ +store <- tempfile(fileext = ".zarr") +create_zarr(store = store) +dir.exists(store) + +} diff --git a/man/create_zarr_group.Rd b/man/create_zarr_group.Rd new file mode 100644 index 00000000..da8f2279 --- /dev/null +++ b/man/create_zarr_group.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/zarr_utils.R +\name{create_zarr_group} +\alias{create_zarr_group} +\title{create_zarr_group} +\usage{ +create_zarr_group(store, name, version = 2) +} +\arguments{ +\item{store}{the location of (zarr) store} + +\item{name}{name of the group} + +\item{version}{zarr version} +} +\description{ +create zarr groups +} diff --git a/man/table-utils.Rd b/man/table-utils.Rd index 5f00263b..66f3a8b0 100644 --- a/man/table-utils.Rd +++ b/man/table-utils.Rd @@ -26,7 +26,15 @@ \S4method{setTable}{SpatialData,ANY}(x, i, ..., name = NULL, rk = "rk", ik = "ik") -\S4method{setTable}{SpatialData,character}(x, i, y, name = NULL, rk = "region", ik = "instance_id") +\S4method{setTable}{SpatialData,character}( + x, + i, + y, + name = NULL, + rk = "region", + ik = "instance_id", + version = "0.1" +) } \arguments{ \item{x}{\code{\link{SpatialData}} object.} diff --git a/man/writeSpatialData.Rd b/man/writeSpatialData.Rd new file mode 100644 index 00000000..56e849c2 --- /dev/null +++ b/man/writeSpatialData.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/write.R +\name{writeSpatialData} +\alias{writeSpatialData} +\alias{writeImage} +\alias{writeLabel} +\alias{writePoint} +\alias{writeShape} +\alias{writeTable} +\title{Writing `SpatialData`} +\usage{ +writeSpatialData(x, path, replace = TRUE, version = "0.2", ...) + +writePoint(x, name, path, replace = TRUE, format = sdFormat("0.1")) + +writeShape(x, name, path, replace = TRUE, format = sdFormat("0.1")) + +writeImage(x, name, path, replace = TRUE, format = sdFormat("0.1")) + +writeLabel(x, name, path, replace = TRUE, format = sdFormat("0.1")) + +writeTable(x, name, path, replace = TRUE, format = sdFormat("0.1")) +} +\arguments{ +\item{x}{For \code{writeSpatialData}, +a \code{SpatialData} +For \code{writeImage/Label/Point/Shape/Table}, +a \code{ImageArray},\code{LabelArray}, +\code{PointFrame}, \code{ShapeFrame}} + +\item{path}{path to zarr store.} + +\item{replace}{if TRUE, existing elements with the same name will be +replaced with the given element} + +\item{version}{SpatialData version, 0.1 (zarr v2) or 0.2 (zarr v3)} + +\item{...}{option arguments passed to and from other methods.} +} +\value{ +\itemize{ +\item{For \code{writeSpatialData}, a \code{SpatialData}.}, +\item{For element writers, a \code{ImageArray}, \code{LabelArray}, +\code{PointFrame}, \code{ShapeFrame}, or \code{SingleCellExperiment}.}} +} +\description{ +Writing `SpatialData` +} diff --git a/tests/testthat/helper-examples.R b/tests/testthat/helper-examples.R new file mode 100644 index 00000000..2bf83990 --- /dev/null +++ b/tests/testthat/helper-examples.R @@ -0,0 +1,52 @@ +require(duckspatial, quietly=TRUE) +require(arrow, quietly=TRUE) + +# seed +set.seed(1) + +example_points <- function(){ + data.frame(x = runif(100), y = runif(100)) +} + +example_circles <- function(){ + duckspatial::as_duckspatial_df( + st_as_sf( + arrow::arrow_table( + geometry = geoarrow::as_geoarrow_vctr( + c( + "POINT (36.382774 24.6331748)", + "POINT (32.378292 46.4148383)", + "POINT (24.3715883 25.5517166)", + "POINT (18.7407733 23.5779362)" + ) + ), + radius = c(4,4,4,4) + ) + ), + conn = duckspatial::ddbs_create_conn(dbdir = "memory"), + wkt = "wkt", + geom_col = "geometry", + remove = TRUE + ) +} + +example_polygons <- function(){ + duckspatial::as_duckspatial_df( + st_as_sf( + arrow::arrow_table( + geometry = geoarrow::as_geoarrow_vctr( + c( + "POLYGON ((4.53 2.11, 5.55 1.43, 5.78 1.33, 6.89 9.10, 4.30 4.15, 3.06 4.29, 4.53 2.11))", + "POLYGON ((4.71 3.73, 7.62 2.48, 9.43 1.09, 9.33 4.99, 6.04 9.35, 4.60 4.85, 4.71 3.73))", + "POLYGON ((1.65 1.09, 5.24 0.64, 7.02 0.62, 7.88 1.70, 3.17 7.55, 2.78 6.20, 1.65 1.09))", + "POLYGON ((1.81 3.73, 2.99 0.28, 3.82 4.77, 2.57 8.80, 1.69 7.71, 1.92 5.27, 1.81 3.73))" + ) + ) + ) + ), + conn = duckspatial::ddbs_create_conn(dbdir = "memory"), + wkt = "wkt", + geom_col = "geometry", + remove = TRUE + ) +} \ No newline at end of file diff --git a/tests/testthat/test-sdarray.R b/tests/testthat/test-sdarray.R index 29923aeb..753afb8e 100644 --- a/tests/testthat/test-sdarray.R +++ b/tests/testthat/test-sdarray.R @@ -24,8 +24,8 @@ test_that("data_type()", { test_that("SpatialDataImage()", { rgb <- \(n) sample(seq_len(255), n, replace=TRUE) mat <- array(rgb(3*20*20), dim=c(3,20,20)) + SpatialDataImage(mat) # invalid - expect_error(SpatialDataImage(mat)) expect_error(SpatialDataImage(mat, 1)) expect_error(SpatialDataImage(mat, list())) # single scale @@ -52,11 +52,135 @@ test_that("data(),SpatialDataImage", { expect_error(data(img, c(1,2))) }) + +test_that("create, SpatialDataImage", { + + # create image + set.seed(1) + img <- array(sample(1:255, size = 100*100*3, replace = TRUE), + dim = c(3,100,100)) + + # make image array + imgarray <- SpatialDataImage(img) + expect_identical(data(imgarray), img) + expect_identical(dim(imgarray),dim(img)) + + # coordinate systems + expect_identical(CTname(imgarray), "global") + expect_identical(CTtype(imgarray), "identity") + imgarray_new <- addCT(imgarray, "test", "scale", c(1,2,2)) + expect_identical(CTname(imgarray_new), c("global", "test")) + expect_identical(CTtype(imgarray_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(images = list(test_image = imgarray)) + expect_identical(data(image(sd)), data(imgarray)) + expect_identical(image(sd), imgarray) + expect_identical(image(sd, 1), imgarray) +}) + +test_that("create multiscale, SpatialDataImage", { + + # create image + set.seed(1) + img <- array(sample(1:255, size = 100*100*3, replace = TRUE), + dim = c(3,100,100)) + + # make image array + imgarray <- SpatialDataImage(img, scale_factors = c(2,2,2)) + expect_identical(data(imgarray), img) + expect_identical(dim(imgarray),dim(img)) + + # coordinate systems + expect_identical(CTname(imgarray), "global") + expect_identical(CTtype(imgarray), "identity") + imgarray_new <- addCT(imgarray, "test", "scale", c(1,2,2)) + expect_identical(CTname(imgarray_new), c("global", "test")) + expect_identical(CTtype(imgarray_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(images = list(test_image = imgarray)) + expect_identical(data(image(sd)), data(imgarray)) + expect_identical(data(image(sd),2), data(imgarray,2)) + expect_identical(data(image(sd),3), data(imgarray,3)) + expect_identical(image(sd), imgarray) + expect_identical(image(sd, 1), imgarray) +}) + +z <- list(0.1, 0.2) + +for (v in names(z)) { + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write, SpatialDataImage", { + + # create image + set.seed(1) + img <- array(sample(1:255, size = 100*100*3, replace = TRUE), + dim = c(3,100,100)) + + # make image array + imgarray <- SpatialDataImage(img, version = image(sdFormat(v))) + sd <- SpatialData(images = list(test_image = imgarray)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + imgarray2 <- image(sd2) + expect_identical(realize(data(imgarray)), + realize(data(imgarray2))) + expect_equal(meta(imgarray), + meta(imgarray2)) + }) + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write multiscale, SpatialDataImage", { + + # create image + set.seed(1) + img <- array(sample(1:255, size = 100*100*3, replace = TRUE), + dim = c(3,100,100)) + + # make image array + imgarray <- SpatialDataImage(img, scale_factors = c(2,2,2), + version = image(sdFormat(v))) + sd <- SpatialData(images = list(test_image = imgarray)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + imgarray2 <- image(sd2) + expect_identical(realize(data(imgarray, 1)), + realize(data(imgarray2, 1))) + expect_identical(realize(data(imgarray, 2)), + realize(data(imgarray2, 2))) + expect_identical(realize(data(imgarray, 3)), + realize(data(imgarray2, 3))) + expect_equal(meta(imgarray),meta(imgarray2)) + }) +} + test_that("SpatialDataLabel()", { val <- sample(seq_len(12), 20*20, replace=TRUE) mat <- array(val, dim=c(20, 20)) + SpatialDataLabel(mat) # invalid - expect_error(SpatialDataLabel(mat)) expect_error(SpatialDataLabel(mat, 1)) expect_error(SpatialDataLabel(mat, list())) # single scale @@ -82,3 +206,126 @@ test_that("data(),SpatialDataLabel", { expect_error(data(lab, "")) expect_error(data(lab, c(1,2))) }) + +test_that("create,SpatialDataLabel", { + + # create label + set.seed(1) + lbl <- array(sample(0:8L, size = 100*100, replace = TRUE), + dim = c(100,100)) + + # make label array + lblarray <- SpatialDataLabel(lbl) + expect_identical(data(lblarray), lbl) + expect_identical(dim(lblarray),dim(lbl)) + + # coordinate systems + expect_identical(CTname(lblarray), "global") + expect_identical(CTtype(lblarray), "identity") + lblarray_new <- addCT(lblarray, "test", "scale", c(2,2)) + expect_identical(CTname(lblarray_new), c("global", "test")) + expect_identical(CTtype(lblarray_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(labels = list(test_label = lblarray)) + expect_identical(data(label(sd)), data(lblarray)) + expect_identical(label(sd), lblarray) + expect_identical(label(sd, 1), lblarray) +}) + +test_that("create multiscale,SpatialDataLabel", { + + # create label + set.seed(1) + lbl <- array(sample(0:8L, size = 100*100, replace = TRUE), + dim = c(100,100)) + + # make label array + lblarray <- SpatialDataLabel(lbl, scale_factors = c(2,2,2)) + expect_identical(data(lblarray), lbl) + expect_identical(dim(lblarray),dim(lbl)) + + # coordinate systems + expect_identical(CTname(lblarray), "global") + expect_identical(CTtype(lblarray), "identity") + lblarray_new <- addCT(lblarray, "test", "scale", c(2,2)) + expect_identical(CTname(lblarray_new), c("global", "test")) + expect_identical(CTtype(lblarray_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(labels = list(test_label = lblarray)) + expect_identical(data(label(sd)), data(lblarray)) + expect_identical(data(label(sd),2), data(lblarray,2)) + expect_identical(data(label(sd),3), data(lblarray,3)) + expect_identical(label(sd), lblarray) + expect_identical(label(sd, 1), lblarray) +}) + +z <- list(0.1, 0.2) + +for (v in names(z)) { + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write,SpatialDataLabel", { + + # create label + set.seed(1) + lbl <- array(sample(0:8L, size = 100*100, replace = TRUE), + dim = c(100,100)) + + # make label array + lblarray <- SpatialDataLabel(lbl, version = label(sdFormat(v))) + sd <- SpatialData(labels = list(test_label = lblarray)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + lblarray2 <- label(sd2) + expect_identical(realize(data(lblarray)), + realize(data(lblarray2))) + expect_equal(meta(lblarray),meta(lblarray2)) + }) + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write multiscale,SpatialDataLabel", { + + # create label + set.seed(1) + lbl <- array(sample(0:8L, size = 100*100, replace = TRUE), + dim = c(100,100)) + + # make label array + lblarray <- SpatialDataLabel(lbl, scale_factors = c(2,2,2), + version = label(sdFormat(v))) + sd <- SpatialData(labels = list(test_label = lblarray)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + lblarray2 <- label(sd2) + expect_identical(realize(data(lblarray)), + realize(data(lblarray2))) + expect_identical(realize(data(lblarray, 2)), + realize(data(lblarray2, 2))) + expect_identical(realize(data(lblarray, 3)), + realize(data(lblarray2, 3))) + expect_equal(meta(lblarray),meta(lblarray2)) + }) + +} \ No newline at end of file diff --git a/tests/testthat/test-sdattrs.R b/tests/testthat/test-sdattrs.R index 75b44bfb..a18ec5a7 100644 --- a/tests/testthat/test-sdattrs.R +++ b/tests/testthat/test-sdattrs.R @@ -47,5 +47,5 @@ for (v in names(z)) { expect_silent(z <- channels(y <- image(x))) expect_length(z, dim(y)[1]) }) - + } \ No newline at end of file diff --git a/tests/testthat/test-sdframe.R b/tests/testthat/test-sdframe.R index c24c87c5..17960326 100644 --- a/tests/testthat/test-sdframe.R +++ b/tests/testthat/test-sdframe.R @@ -105,3 +105,176 @@ test_that("as.data.frame", { expect_equal(names(y), names(p)) expect_identical(y, as.data.frame(collect(data(p)))) }) + +test_that("create, SpatialDataPoint", { + + # make point frame + df <- example_points() + pf <- SpatialDataPoint(df) + expect_identical(st_coordinates(st_as_sf(data(pf))), + { + dfm <- as.matrix(df) + colnames(dfm) <- c("X", "Y") + dfm + }) + expect_equal(dim(pf), c(100,1)) # geometry column of POINT + expect_identical(names(pf), "geometry") + + # coordinate systems + expect_identical(CTname(pf), "global") + expect_identical(CTtype(pf), "identity") + pf_new <- addCT(pf, "test", "scale", c(2,2)) + expect_identical(CTname(pf_new), c("global", "test")) + expect_identical(CTtype(pf_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(points = list(test_points = pf)) + expect_identical(data(point(sd)), data(pf)) + expect_identical(point(sd), pf) + expect_identical(point(sd, 1), pf) +}) + +test_that("create polygon, SpatialDataShape", { + + # make point frame + df <- example_polygons() + pf <- SpatialDataShape(df) + expect_identical(data(pf), df) + expect_identical(dim(pf),dim(ddbs_collect(df))) + expect_identical(names(pf), colnames(df)) + expect_identical(ddbs_collect(data(pf[1:2,1])), + ddbs_collect(df)[1:2,1]) + + # coordinate systems + expect_identical(CTname(pf), "global") + expect_identical(CTtype(pf), "identity") + pf_new <- addCT(pf, "test", "scale", c(2,2)) + expect_identical(CTname(pf_new), c("global", "test")) + expect_identical(CTtype(pf_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(shapes = list(test_shapes = pf)) + expect_identical(data(shape(sd)), data(pf)) + expect_identical(shape(sd), pf) + expect_identical(shape(sd, 1), pf) +}) + +test_that("create circle, SpatialDataShape", { + + # make point frame + df <- example_circles() + pf <- SpatialDataShape(df) + expect_identical(data(pf), df) + expect_identical(dim(pf),dim(ddbs_collect(df))) + expect_identical(names(pf), colnames(df)) + expect_identical(ddbs_collect(data(pf[1:2,1])), + ddbs_collect(df)[1:2,1]) + + # coordinate systems + expect_identical(CTname(pf), "global") + expect_identical(CTtype(pf), "identity") + pf_new <- addCT(pf, "test", "scale", c(2,2)) + expect_identical(CTname(pf_new), c("global", "test")) + expect_identical(CTtype(pf_new), c("identity", "scale")) + + # make spatial data + sd <- SpatialData(shapes = list(test_shapes = pf)) + expect_identical(data(shape(sd)), data(pf)) + expect_identical(shape(sd), pf) + expect_identical(shape(sd, 1), pf) +}) + +z <- list(0.1, 0.2) + +for (v in z) { + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write, SpatialDataPoint", { + + # make sd data + df <- example_points() + pf <- SpatialDataPoint(df, version = point(sdFormat(v))) + sd <- SpatialData(points = list(test_points = pf)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + pf2 <- point(sd2) + # attr(data(pf), "source_table") is not identical, obviously + expect_equal( + ddbs_collect(data(pf)), + ddbs_collect(data(pf2)) + ) + expect_identical(st_coordinates(st_as_sf(data(pf))), + st_coordinates(st_as_sf(data(pf2)))) + expect_identical(meta(pf),meta(pf2)) + expect_identical(names(pf), names(pf2)) + }) + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write polygon, SpatialDataShape", { + + # make sd data + df <- example_polygons() + pf <- SpatialDataShape(df, version = shape(sdFormat(v))) + sd <- SpatialData(shapes = list(test_shapes = pf)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + pf2 <- shape(sd2) + # TODO: they are not identical, why ? + expect_equal(data(pf) |> collect(), + data(pf2) |> collect()) + expect_identical(meta(pf),meta(pf2)) + expect_identical(names(pf), names(pf2)) + expect_identical(data(pf[1:2, 1]) |> collect(), + data(pf2[1:2,1]) |> collect()) + }) + + td <- tempdir() + zarr.store <- "test.zarr" + zarr.path <- file.path(td, zarr.store) + unlink(zarr.path, recursive = TRUE) + + test_that("write circle, SpatialDataShape", { + + # make sd data + df <- example_circles() + pf <- SpatialDataShape(df, version = shape(sdFormat(v))) + sd <- SpatialData(shapes = list(test_shapes = pf)) + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + pf2 <- shape(sd2) + # TODO: they are not identical, why ? + expect_equal(data(pf) |> collect(), + data(pf2) |> collect()) + expect_identical(meta(pf),meta(pf2)) + expect_identical(names(pf), names(pf2)) + expect_identical(data(pf[1:2, 1]) |> collect(), + data(pf2[1:2,1]) |> collect()) + }) + +} \ No newline at end of file diff --git a/tests/testthat/test-tables.R b/tests/testthat/test-tables.R index ee7a8fc3..61504ff6 100644 --- a/tests/testthat/test-tables.R +++ b/tests/testthat/test-tables.R @@ -1,4 +1,5 @@ require(SingleCellExperiment, quietly=TRUE) +require(anndataR, quietly=TRUE) x <- file.path("extdata", "blobs.zarr") x <- system.file(x, package="SpatialData") @@ -169,3 +170,48 @@ test_that("setTable() fails with invalid inputs", { # Non-existent element expect_error(setTable(x, "non_existent", SingleCellExperiment()), "is not an element of 'x'") }) + +# TODO: update once v3 zarr anndata write support is implemented +v <- 0.1 + +test_that("write, Table (SCE) for shapes", { + + # make sd data + i <- "test_shapes" + df <- example_polygons() + pf <- SpatialDataShape(df, version = shape(sdFormat(v))) + sd <- SpatialData(shapes = setNames(list(pf), i)) + + # create table (SCE) + set.seed(1) + n <- 30 + mat <- matrix(1:(nrow(pf)*n), nrow = n) + sce <- SingleCellExperiment(assays = list(counts = mat)) + + # set table to SpatialData + e <- element(sd, i) + + # set instances and region key manually + int_colData(sce)$instance_id <- instances(e) + colData(sce)$instance_id <- instances(e) + colData(sce)[["region"]] <- i + + # set new table + sd <- setTable(sd, i, name = "test_tables", sce, + version = "0.1") + + # write to location + zarr.path <- tempfile(fileext = ".zarr") + writeSpatialData(sd, path = zarr.path, version = v) + expect_true(dir.exists(zarr.path)) + + # read back and compare + sd2 <- readSpatialData(zarr.path) + expect_true(hasTable(sd2, i)) + expect_contains(meta(table(sd2)), + meta(table(sd))) +}) + +test_that("write, Table (SCE) for labels", { + skip("write, Table (SCE) for labels is not implemented yet!") +}) \ No newline at end of file diff --git a/tests/testthat/test-write.R b/tests/testthat/test-write.R new file mode 100644 index 00000000..e69de29b diff --git a/tests/testthat/test-zarrutils.R b/tests/testthat/test-zarrutils.R new file mode 100644 index 00000000..12b8227f --- /dev/null +++ b/tests/testthat/test-zarrutils.R @@ -0,0 +1,139 @@ +test_that("create zarr/group", { + + # open zarr + store <- tempfile(fileext = ".zarr") + create_zarr(store = store) + expect_true(dir.exists(store)) + expect_true(file.exists(file.path(store, ".zgroup"))) + + # create group one group + create_zarr_group(store = store, name = "group1") + expect_true(dir.exists(file.path(store, "group1"))) + expect_true(file.exists(file.path(store, "group1", ".zgroup"))) + + # create nested two groups + create_zarr_group(store = store, name = "group2/subgroup1") + expect_true(dir.exists(file.path(store, "group2"))) + expect_true(file.exists(file.path(store, "group2", ".zgroup"))) + expect_true(dir.exists(file.path(store, "group2/subgroup1"))) + expect_true(file.exists(file.path(store, "group2/subgroup1", ".zgroup"))) + + # create nested three groups + create_zarr_group(store = store, name = "group3/subgroup1/subsubgroup1") + expect_true(dir.exists(file.path(store, "group3"))) + expect_true(file.exists(file.path(store, "group3", ".zgroup"))) + expect_true(dir.exists(file.path(store, "group3/subgroup1"))) + expect_true(file.exists(file.path(store, "group3/subgroup1", ".zgroup"))) + expect_true(dir.exists(file.path(store, "group3/subgroup1/subsubgroup1"))) + expect_true(file.exists(file.path(store, "group3/subgroup1/subsubgroup1", ".zgroup"))) + + # invalid version string + dir.create(td <- tempfile()) + name <- "test" + expect_error(create_zarr(dir = td, name = name, version = "4"), pattern = "version must be '2' or '3'") +}) + +test_that("create zarr/group v3", { + + # open v3 zarr store + store <- tempfile(fileext = ".zarr") + create_zarr(store = store, version = 3) + expect_true(dir.exists(store)) + expect_true(file.exists(file.path(store, "zarr.json"))) + expect_false(file.exists(file.path(store, ".zgroup"))) + + # check zarr.json exists and attributes are empty + expect_true(file.exists(file.path(store, "zarr.json"))) + expect_equal(Rarr::read_zarr_attributes(store), list()) + + # create a sub-group + create_zarr_group(store = store, name = "images", version = 3) + expect_true(file.exists(file.path(store, "images", "zarr.json"))) + expect_false(file.exists(file.path(store, "images", ".zgroup"))) + + # create nested groups — parent group should also be v3 + create_zarr_group(store = store, name = "points/blobs_points", version = 3) + expect_true(file.exists(file.path(store, "points", "zarr.json"))) + expect_true(file.exists(file.path(store, "points/blobs_points", "zarr.json"))) +}) + + +# create a v2 zarr array for the v2 zattrs tests +dir.create(td <- tempfile()) +path <- file.path(td, "test.zarr") +x <- array(runif(n = 10), dim = c(2, 5)) +Rarr::write_zarr_array( + x = x, zarr_array_path = path, + chunk_dim = c(2, 5), zarr_version = 2L +) + +test_that("read/write zattrs", { + + # add .zattrs to / + zattrs <- list(foo = "foo", bar = "bar") + Rarr::write_zarr_attributes(path, new.zattrs = zattrs) + expect_true(file.exists(file.path(path, ".zattrs"))) + + # check .zattrs + read.zattrs <- Rarr::read_zarr_attributes(path) + expect_equal(read.zattrs, zattrs) + + # add new elements to .zattrs + zattrs.new.elem <- list(foo2 = "foo") + Rarr::write_zarr_attributes(path, new.zattrs = zattrs.new.elem) + read.zattrs <- Rarr::read_zarr_attributes(path) + expect_equal(read.zattrs, c(zattrs,zattrs.new.elem)) + + # overwrite + zattrs.new.elem <- list(foo2 = "foo2") + Rarr::write_zarr_attributes(path, new.zattrs = zattrs.new.elem) + read.zattrs <- Rarr::read_zarr_attributes(path) + zattrs[names(zattrs.new.elem)] <- zattrs.new.elem + expect_equal(read.zattrs, c(zattrs)) + + # overwrite = FALSE + zattrs.new.elem <- list(foo2 = "foo") + Rarr::write_zarr_attributes(path, new.zattrs = zattrs.new.elem, overwrite = FALSE) + read.zattrs <- Rarr::read_zarr_attributes(path) + zattrs[names(zattrs.new.elem)] <- "foo2" + expect_contains(read.zattrs, zattrs) + expect_contains(zattrs, read.zattrs) +}) + +test_that("read/write zattrs v3", { + + # create a v3 zarr group to use as the target path + dir.create(td <- tempfile()) + grp <- file.path(td, "elem") + create_zarr_group(store = td, name = "elem", version = 3) + + # write attributes into zarr.json + zattrs <- list(foo = "foo", bar = "bar") + Rarr::write_zarr_attributes(grp, new.zattrs = zattrs) + expect_true(file.exists(file.path(grp, "zarr.json"))) + expect_false(file.exists(file.path(grp, ".zattrs"))) + + # read back attributes from zarr.json + read.zattrs <- Rarr::read_zarr_attributes(grp) + expect_equal(read.zattrs, zattrs) + + # zarr.json must still exist (zarr_format / node_type preserved by Rarr internally) + expect_true(file.exists(file.path(grp, "zarr.json"))) + + # add new element + Rarr::write_zarr_attributes(grp, new.zattrs = list(baz = "baz")) + read.zattrs <- Rarr::read_zarr_attributes(grp) + expect_equal(read.zattrs, c(zattrs, list(baz = "baz"))) + + # overwrite existing key + Rarr::write_zarr_attributes(grp, new.zattrs = list(foo = "FOO")) + read.zattrs <- Rarr::read_zarr_attributes(grp) + expect_equal(read.zattrs$foo, "FOO") + expect_equal(read.zattrs$bar, "bar") # untouched key preserved + + # overwrite = FALSE should not overwrite existing key + Rarr::write_zarr_attributes(grp, new.zattrs = list(foo = "original"), overwrite = FALSE) + read.zattrs <- Rarr::read_zarr_attributes(grp) + expect_equal(read.zattrs$foo, "FOO") # unchanged + +}) \ No newline at end of file