initial support for read_mdim() to work with proxy = TRUE

see #559 for motivation

NAMESPACE

@@ -5,6 +5,7 @@ S3method("[",dimension)
 S3method("[",dimensions)
 S3method("[",intervals)
 S3method("[",intervals_list)
+S3method("[",mdim)
 S3method("[",nc_proxy)
 S3method("[",stars)
 S3method("[",stars_proxy)
@@ -104,6 +105,7 @@ S3method(st_as_stars,data.frame)
 S3method(st_as_stars,default)
 S3method(st_as_stars,im)
 S3method(st_as_stars,list)
+S3method(st_as_stars,mdim)
 S3method(st_as_stars,nc_proxy)
 S3method(st_as_stars,ncdfgeom)
 S3method(st_as_stars,sf)
@@ -116,6 +118,7 @@ S3method(st_bbox,dimensions)
 S3method(st_bbox,stars)
 S3method(st_coordinates,dimensions)
 S3method(st_coordinates,stars)
+S3method(st_crop,mdim)
 S3method(st_crop,nc_proxy)
 S3method(st_crop,stars)
 S3method(st_crop,stars_proxy)

NEWS.md

@@ -1,6 +1,6 @@
 # version 0.6-7
 
-* 
+* initial support for `read_mdim()` to work with `proxy = TRUE`; #559 
 
 # version 0.6-6
 

R/mdim.R

@@ -116,10 +116,10 @@ match_raster_dims = function(nms) {
 #' @param variable name of the array to be read; if `"?"`, a list of array names is returned, with group name as list element names.
 #' @param options character; driver specific options regarding the opening (read_mdim) or creation (write_mdim) of the dataset
 #' @param raster names of the raster variables (default: first two dimensions)
-#' @param offset integer; offset for each dimension (pixels) of sub-array to read, defaults to 0 for each dimension(requires sf >= 1.0-9)
+#' @param offset integer; zero-based offset for each dimension (pixels) of sub-array to read, defaults to 0 for each dimension(requires sf >= 1.0-9)
 #' @param count integer; size for each dimension (pixels) of sub-array to read (default: read all); a value of NA will read the corresponding dimension entirely; counts are relative to the step size (requires sf >= 1.0-9)
 #' @param step integer; step size for each dimension (pixels) of sub-array to read; defaults to 1 for each dimension (requires sf >= 1.0-9)
-#' @param proxy logical; return proxy object? (not functional yet)
+#' @param proxy logical; return proxy object?
 #' @param debug logical; print debug info?
 #' @param bounds logical or character: if \code{TRUE} tries to infer from "bounds" attribute; if character, 
 #' named vector of the form \code{c(longitude="lon_bnds", latitude="lat_bnds")} with names dimension names
@@ -132,9 +132,6 @@ read_mdim = function(filename, variable = character(0), ..., options = character
 					 raster = NULL, offset = integer(0), count = integer(0), step = integer(0), proxy = FALSE, 
 					 debug = FALSE, bounds = TRUE, curvilinear = NA) {
 
-	if (proxy)
-		stop("proxy not yet implemented in read_mdim()")
-
 	stopifnot(is.character(filename), is.character(variable), is.character(options))
 	ret = gdal_read_mdim(filename, variable, options, rev(offset), rev(count), rev(step), proxy, debug)
 
@@ -215,33 +212,41 @@ read_mdim = function(filename, variable = character(0), ..., options = character
 		else
 			x
 	}
-	lst = lapply(ret$array_list, function(x) structure(clean_units(x), dim = rev(dim(x))))
 
 	# create return object:
-	st = st_stars(lst, dimensions)
-	if (is.null(ret$srs) || is.na(ret$srs)) {
-		if (missing(curvilinear) || is.character(curvilinear)) { # try curvilinear:
-			xy = raster$dimensions
-			ll = curvilinear
-			if (is.character(curvilinear)) {
-				if (is.null(names(curvilinear)))
-					names(curvilinear) = xy[1:2]
-				ret = try(st_as_stars(st, curvilinear = curvilinear), silent = TRUE)
-				if (inherits(ret, "stars"))
-					st = ret
+	if (proxy) {
+		lst = lapply(ret$array_list, function(x) filename)
+		st = st_stars_proxy(lst, dimensions, NA_value = NA_real_, resolutions = NULL, class = "mdim")
+		if (!is.null(ret$srs))
+			st_set_crs(st, ret$srs)
+		else
+			st
+	} else { 
+		lst = lapply(ret$array_list, function(x) structure(clean_units(x), dim = rev(dim(x))))
+		st = st_stars(lst, dimensions)
+		if (is.null(ret$srs) || is.na(ret$srs)) {
+			if (missing(curvilinear) || is.character(curvilinear)) { # try curvilinear:
+				xy = raster$dimensions
+				ll = curvilinear
+				if (is.character(curvilinear)) {
+					if (is.null(names(curvilinear)))
+						names(curvilinear) = xy[1:2]
+					ret = try(st_as_stars(st, curvilinear = curvilinear), silent = TRUE)
+					if (inherits(ret, "stars"))
+						st = ret
+				}
+				if (!is_curvilinear(st) &&
+					inherits(x <- try(read_mdim(filename, ll[1], curvilinear = FALSE), silent = TRUE), "stars") &&
+					inherits(y <- try(read_mdim(filename, ll[2], curvilinear = FALSE), silent = TRUE), "stars") &&
+						identical(dim(x)[xy], dim(st)[xy]) && identical(dim(y)[xy], dim(st)[xy]))
+					st = st_as_stars(st, curvilinear = setNames(list(x, y), xy))
 			}
-			if (!is_curvilinear(st) &&
-				inherits(x <- try(read_mdim(filename, ll[1], curvilinear = FALSE), silent = TRUE), "stars") &&
-				inherits(y <- try(read_mdim(filename, ll[2], curvilinear = FALSE), silent = TRUE), "stars") &&
-					identical(dim(x)[xy], dim(st)[xy]) && identical(dim(y)[xy], dim(st)[xy]))
-				st = st_as_stars(st, curvilinear = setNames(list(x, y), xy))
-		}
-		st
-	} else
-		st_set_crs(st, ret$srs)
+			st
+		} else
+			st_set_crs(st, ret$srs)
+	}
 }
 
-
 # WRITE helper functions:
 
 add_attr = function(x, at) { # append at to attribute "attrs"
@@ -425,3 +430,37 @@ write_mdim = function(x, filename, driver = detect.driver(filename), ...,
 					as_float = as_float)
 	invisible(x)
 }
+
+#' @export
+st_as_stars.mdim = function(.x, ..., downsample = 0, debug = FALSE,
+		envir = parent.frame()) {
+	d = st_dimensions(.x)
+	step = rep_len(downsample, length(dim(.x))) + 1
+	offset = sapply(d, function(x) x$from) - 1
+	count = (sapply(d, function(x) x$to) - offset) %/% step
+	if (debug)
+		print(data.frame(offset = offset, step = step, count = count))
+	# read:
+	l = vector("list", length(.x))
+	for (i in seq_along(l))
+		l[[i]] = read_mdim(.x[[i]], names(.x)[i], offset = offset, count = count, step = step, proxy = FALSE, ...)
+	ret = do.call(c, l)
+	if (is.na(st_crs(ret)) && !is.na(st_crs(.x))) # as it was set after reading
+		st_crs(ret) = st_crs(.x)
+	# post-process:
+	cl = attr(.x, "call_list")
+	if (!all(downsample == 0))
+		lapply(cl, check_xy_warn, dimensions = st_dimensions(.x))
+	process_call_list(ret, cl, envir = envir, downsample = downsample)
+}
+
+#' @export
+`[.mdim` = function(x, i, j, ...) {
+	structure(NextMethod(), class = class(x))
+}
+
+#' @export
+#' @name st_crop
+st_crop.mdim = function(x, y, ...) {
+	structure(NextMethod(), class = class(x))
+}

R/proxy.R

@@ -54,16 +54,17 @@ plot.stars_proxy = function(x, y, ..., downsample = get_downsample(dim(x))) {
 	plot(st_as_stars(x, downsample = downsample, ...), ..., downsample = 0)
 }
 
-st_stars_proxy = function(x, dimensions, ..., NA_value, resolutions, RasterIO = list(), file_dim = NULL) {
-	stopifnot(!missing(NA_value))
-	stopifnot(!missing(resolutions))
-	stopifnot(length(list(...)) == 0)
-	stopifnot(is.list(x))
-	stopifnot(inherits(dimensions, "dimensions"))
+st_stars_proxy = function(x, dimensions, ..., NA_value, resolutions, RasterIO = list(), file_dim = NULL, class = character(0)) {
+	stopifnot(!missing(NA_value),
+		!missing(resolutions),
+		length(list(...)) == 0,
+		is.list(x),
+		inherits(dimensions, "dimensions"))
+
 	if (length(RasterIO) == 0)
 		RasterIO = NULL
 	structure(x, dimensions = dimensions, NA_value = NA_value, resolutions = resolutions,
-		RasterIO = RasterIO, file_dim = file_dim, class = c("stars_proxy", "stars"))
+		RasterIO = RasterIO, file_dim = file_dim, class = c(class, "stars_proxy", "stars"))
 }
 
 add_resolution = function(lst) {
@@ -336,7 +337,7 @@ st_as_stars.stars_proxy = function(.x, ..., downsample = 0, url = attr(.x, "url"
 		# there are cases where this is not right. Hence:
 		# TODO: only warn when there is a reason to warn.
 		if (!all(downsample == 0))
-			lapply(attr(.x, "call_list"), check_xy_warn, dimensions = st_dimensions(.x))
+			lapply(cl, check_xy_warn, dimensions = st_dimensions(.x))
 		process_call_list(fetch(.x, ..., downsample = downsample), cl, envir = envir, downsample = downsample)
 	}
 }

R/subset.R

@@ -308,8 +308,10 @@ st_crop.stars = function(x, y, ..., crop = TRUE, epsilon = sqrt(.Machine$double.
 		for (i in seq_along(x))
 			x[[i]][mask] = NA
 	}
-	if (normalize[1]) x = st_normalize(x)
-	x
+	if (normalize[1]) 
+		st_normalize(x)
+	else 
+		x
 }
 
 #' @export