diff --git a/NAMESPACE b/NAMESPACE
index cde8cefe..99ccc7c6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -25,6 +25,12 @@ export(readShape)
 export(readSpatialData)
 export(readTable)
 export(select)
+export(writeImage)
+export(writeLabel)
+export(writePoint)
+export(writeShape)
+export(writeSpatialData)
+export(writeTable)
 exportMethods("$")
 exportMethods("$<-")
 exportMethods("[")
@@ -124,17 +130,22 @@ importFrom(BiocGenerics,table)
 importFrom(BiocGenerics,transform)
 importFrom(DBI,dbIsValid)
 importFrom(DelayedArray,DelayedArray)
+importFrom(EBImage,resize)
 importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,sparseVector)
 importFrom(Matrix,summary)
 importFrom(Matrix,t)
 importFrom(RBGL,sp.between)
 importFrom(Rarr,read_zarr_attributes)
+importFrom(Rarr,write_zarr_array)
+importFrom(Rarr,write_zarr_attributes)
+importFrom(Rarr,write_zarr_group)
 importFrom(Rarr,zarr_overview)
 importFrom(S4Vectors,"metadata<-")
 importFrom(S4Vectors,DataFrame)
 importFrom(S4Vectors,SimpleList)
 importFrom(S4Vectors,coolcat)
+importFrom(S4Vectors,isSequence)
 importFrom(S4Vectors,make_zero_col_DFrame)
 importFrom(S4Vectors,metadata)
 importFrom(S4Vectors,setValidity2)
@@ -154,8 +165,10 @@ importFrom(ZarrArray,ZarrArray)
 importFrom(ZarrArray,path)
 importFrom(ZarrArray,type)
 importFrom(anndataR,read_zarr)
+importFrom(anndataR,write_zarr)
 importFrom(dplyr,.data)
 importFrom(dplyr,all_of)
+importFrom(dplyr,bind_cols)
 importFrom(dplyr,coalesce)
 importFrom(dplyr,collect)
 importFrom(dplyr,count)
@@ -176,6 +189,7 @@ importFrom(duckspatial,ddbs_bbox)
 importFrom(duckspatial,ddbs_create_conn)
 importFrom(duckspatial,ddbs_intersects)
 importFrom(duckspatial,ddbs_open_dataset)
+importFrom(duckspatial,ddbs_write_dataset)
 importFrom(duckspatial,ddbs_write_table)
 importFrom(graph,"edgeData<-")
 importFrom(graph,"edgeDataDefaults<-")
@@ -214,6 +228,7 @@ importFrom(sf,st_point)
 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 10510bb1..d5e8587d 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -97,3 +97,17 @@ setClass("SpatialDataFrame",
 setClassUnion("SpatialDataElement", c(
     "SpatialDataImage", "SpatialDataLabel", 
     "SpatialDataPoint", "SpatialDataShape"))
+
+.sdFormat <- setClass(
+  Class = "sdFormat",
+  slots = list(
+    version = "character",
+    zarr_version = "integer",
+    ome_version = "character",
+    image = "character",
+    label = "character",
+    point = "character",
+    shape = "character",
+    table = "character"
+  )
+)
diff --git a/R/AllGenerics.R b/R/AllGenerics.R
index 8d7901ea..bce01399 100644
--- a/R/AllGenerics.R
+++ b/R/AllGenerics.R
@@ -106,3 +106,10 @@ setGeneric("datasets", \(x, ...) standardGeneric("datasets"))
 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/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/sdArray.R b/R/sdArray.R
index 937546ae..0832f818 100644
--- a/R/sdArray.R
+++ b/R/sdArray.R
@@ -65,22 +65,61 @@ NULL
 
 #' @export
 #' @rdname SpatialDataArray
+#' @importFrom methods new
 #' @importFrom S4Vectors metadata<-
-SpatialDataImage <- \(data=list(), meta=SpatialDataAttrs(type="image"), metadata=list(), ...) {
-    if (is.array(data)) data <- list(data)
-    x <- .SpatialDataImage(data=data, meta=meta, ...)
-    metadata(x) <- metadata
-    return(x)
+SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(type="image"),
+                             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(type = "image", 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 <- \(data=list(), meta=SpatialDataAttrs(type="label"), metadata=list(), ...) {
-    if (is.array(data)) data <- list(data)
-    x <- .SpatialDataLabel(data=data, meta=meta, ...)
-    metadata(x) <- metadata
-    return(x)
+SpatialDataLabel <- function(data=list(), 
+                             meta=SpatialDataAttrs(type="label"),
+                             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(type = "label", scale_factors = scale_factors) 
+  }
+  x <- .SpatialDataLabel(data=data, meta=meta, ...)
+  metadata(x) <- metadata
+  
+  # update version if provided
+  if(!is.null(version))
+    version(x) <- version
+  return(x)
 }
 
 # utils ----
@@ -111,6 +150,76 @@ setMethod("data_type", "DelayedArray", \(x) {
     return(df$data_type)
 })
 
+# multiscales ----
+
+#' @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!
diff --git a/R/sdAttrs.R b/R/sdAttrs.R
index e70c4caf..1c99958d 100644
--- a/R/sdAttrs.R
+++ b/R/sdAttrs.R
@@ -61,7 +61,7 @@
 #' 
 #' @export
 SpatialDataAttrs <- \(x, type=c("image", "label", "frame"), 
-    trans=NULL, ver="0.4", dim=2, nch=3, ...) 
+    trans=NULL, ver="0.4", dim=2, nch=3, scale_factors = NULL, ...) 
 {
     stopifnot(
         length(dim) == 1, is.numeric(dim), dim %in% seq(2, 4),
@@ -72,17 +72,30 @@ SpatialDataAttrs <- \(x, type=c("image", "label", "frame"),
     ax <- .default_ax(type, dim)
     # transformations:
     ct <- trans %||% .default_ct(ax)
+    # datasets:
+    ds <- .default_ds(.ax_names(ax), scale_factors) 
     # .zattrs list:
     if (type != "frame") {
-        # 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(type == "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)
@@ -101,10 +114,10 @@ SpatialDataAttrs <- \(x, type=c("image", "label", "frame"),
     switch(match.arg(type), 
         # xyzt for points/shapes
         frame={
-            ax <- list(x, y)
+            ax <- list(x$name, y$name)
             if (dim > 2) {
-                ax <- c(ax, list(z))
-                if (dim > 3) ax <- c(ax, list(t))
+                ax <- c(ax, list(z$name))
+                if (dim > 3) ax <- c(ax, list(t$name))
             }
         },
         # tczyx for images/labels
@@ -120,6 +133,15 @@ SpatialDataAttrs <- \(x, type=c("image", "label", "frame"),
     return(ax)
 }
 
+# Internal helper to get axes names
+.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)
@@ -127,6 +149,26 @@ SpatialDataAttrs <- \(x, type=c("image", "label", "frame"),
     list(ct)
 }
 
+# Internal helper to generate datasets
+.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)
@@ -320,3 +362,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 5f24f14c..79332510 100644
--- a/R/sdFrame.R
+++ b/R/sdFrame.R
@@ -123,7 +123,9 @@ NULL
 #' @importFrom methods is
 #' @importFrom sf st_geometry_type
 #' @importFrom S4Vectors metadata<-
-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")
     if (isTRUE(nrow(data) > 0L)) {
         gt <- tryCatch(unique(st_geometry_type(data)), error=\(.) "n/a")
@@ -146,6 +148,10 @@ SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=l
     }
     x <- .SpatialDataPoint(data=data, meta=SpatialDataAttrs(za), ...)
     metadata(x) <- metadata
+    
+    # update version if provided
+    if(!is.null(version))
+      version(x) <- version
     return(x)
 }
 
@@ -153,12 +159,18 @@ SpatialDataPoint <- \(data=NULL, meta=SpatialDataAttrs(type="frame"), metadata=l
 #' @rdname SpatialDataFrame
 #' @importFrom methods is
 #' @importFrom S4Vectors metadata<-
-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")
-    if (!is(data, "duckspatial_df")) 
-        data <- .duck(data, "sdShape")
+    if (!is(data, "duckspatial_df"))
+      data <- .duck(data, "sdShape")
     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/write.R b/R/write.R
new file mode 100644
index 00000000..cfbce0ba
--- /dev/null
+++ b/R/write.R
@@ -0,0 +1,251 @@
+#' @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",
+                             ...) {
+    fmt <- sdFormat(version)
+    zs <- .write_replace_zarr_group(path, 
+                                    "",
+                                    replace, 
+                                    version=zarr_version(fmt))
+
+    # write root-level spatialdata_attrs for v3 
+    # (scverse/spatialdata uses this to pick the read path)
+    if (version == "0.2")
+        Rarr::write_zarr_attributes(zs, new.zattrs=list(
+            spatialdata_attrs=list(version=version),
+            spatialdata_software_version=
+                paste0("SpatialData v", packageVersion("spatialdataR"))))
+
+    # helper for layer writing
+    .writeLayer <- \(l) {
+        s <- substr(l, 1, nchar(l)-1)
+        g <- match.fun(s)
+        f <- match.fun(paste0("write", 
+                              toupper(substr(s, 1, 1)), 
+                              substr(s, 2, nchar(s))))
+        nms <- match.fun(paste0(s, "Names"))(x)
+        lapply(nms, \(.) f(g(x, .), ., path=zs, replace=replace, format=fmt))
+    }
+
+    invisible(lapply(.LAYERS, .writeLayer))
+}
+
+#' @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 <- .write_replace_zarr_group(file.path(path, "points"),
+                                          name, 
+                                          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")
+}
+
+#' @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 <- .write_replace_zarr_group(file.path(path, "shapes"), 
+                                          name,
+                                          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 <- .write_replace_zarr_group(file.path(path, "images"), 
+                                          name,
+                                          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 <- .write_replace_zarr_group(file.path(path, "labels"),
+                                          name,
+                                          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 <- .write_replace_zarr_group(file.path(path, "tables"), 
+                                          name, 
+                                          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")
+}
+
+# utils ----
+
+#' @noRd
+#' @importFrom Rarr write_zarr_group
+.write_replace_zarr_group <- function(path, name, replace, version){
+  ng <- file.path(path, name)
+  
+  # create element parent dir
+  if(!dir.exists(ng))
+    dir.create(ng, recursive = TRUE)
+  
+  # check element dir
+  if(replace){
+    unlink(ng, recursive = TRUE)
+  } else {
+    stop("Directory \"", ng, "\" already exists. ",
+         "Use 'replace=TRUE' to replace it. ",
+         "Its content will be lost!")
+  }
+  
+  # create group
+  write_zarr_group(path, name, zarr_version = version)
+  
+  return(ng)
+}
+
+#' @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)
+}
diff --git a/man/SpatialDataArray.Rd b/man/SpatialDataArray.Rd
index 8b26c4ce..c4e6298a 100644
--- a/man/SpatialDataArray.Rd
+++ b/man/SpatialDataArray.Rd
@@ -20,14 +20,18 @@
 SpatialDataImage(
   data = list(),
   meta = SpatialDataAttrs(type = "image"),
+  version = image(sdFormat(0.1)),
   metadata = list(),
+  scale_factors = NULL,
   ...
 )
 
 SpatialDataLabel(
   data = list(),
   meta = SpatialDataAttrs(type = "label"),
+  version = image(sdFormat(0.1)),
   metadata = list(),
+  scale_factors = NULL,
   ...
 )
 
diff --git a/man/SpatialDataAttrs.Rd b/man/SpatialDataAttrs.Rd
index 839ba4e2..7eef3328 100644
--- a/man/SpatialDataAttrs.Rd
+++ b/man/SpatialDataAttrs.Rd
@@ -43,6 +43,7 @@ SpatialDataAttrs(
   ver = "0.4",
   dim = 2,
   nch = 3,
+  scale_factors = NULL,
   ...
 )
 
diff --git a/man/SpatialDataFrame.Rd b/man/SpatialDataFrame.Rd
index 32709981..18339d3c 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/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 c62eda12..eb22c404 100644
--- a/tests/testthat/test-sdarray.R
+++ b/tests/testthat/test-sdarray.R
@@ -82,3 +82,235 @@ test_that("data(),SpatialDataLabel", {
     expect_error(data(lab, ""))
     expect_error(data(lab, 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)) {
+
+  test_that("write, SpatialDataImage", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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))
+  })
+  
+  test_that("write multiscale, SpatialDataImage", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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("create,SpatialDataLabel", {
+  zarr.path <- tempfile(fileext = ".zarr")
+  
+  # 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)) {
+  
+  test_that("write,SpatialDataLabel", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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))
+  })
+  
+  test_that("write multiscale,SpatialDataLabel", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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))
+  })
+  
+}
diff --git a/tests/testthat/test-sdframe.R b/tests/testthat/test-sdframe.R
index e2851765..754fc15c 100644
--- a/tests/testthat/test-sdframe.R
+++ b/tests/testthat/test-sdframe.R
@@ -105,3 +105,165 @@ 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) {
+  
+  test_that("write, SpatialDataPoint", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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))
+  })
+  
+  test_that("write polygon, SpatialDataShape", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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())
+  })
+  
+  test_that("write circle, SpatialDataShape", {
+    zarr.path <- tempfile(fileext = ".zarr")
+    
+    # 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())
+  })
+  
+}
+