From 6d0f2575d9bf8b48d6bc0548af1482a1f9f11009 Mon Sep 17 00:00:00 2001 From: edzer Date: Wed, 24 Jul 2024 17:45:47 +0200 Subject: [PATCH] initial support for `read_mdim()` to work with `proxy = TRUE` see #559 for motivation --- NAMESPACE | 3 ++ NEWS.md | 2 +- R/mdim.R | 93 +++++++++++++++++++++++++++++++++++--------------- R/proxy.R | 17 ++++----- R/subset.R | 6 ++-- man/mdim.Rd | 4 +-- man/st_crop.Rd | 5 ++- 7 files changed, 89 insertions(+), 41 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 7cb378d3..f0733ea8 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 6696abfb..c9ba5460 100644 --- a/NEWS.md +++ b/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 diff --git a/R/mdim.R b/R/mdim.R index b027c891..b42f32d9 100644 --- a/R/mdim.R +++ b/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)) +} diff --git a/R/proxy.R b/R/proxy.R index 2c398b06..8e617a4d 100644 --- a/R/proxy.R +++ b/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) } } diff --git a/R/subset.R b/R/subset.R index 7bbd564d..dae08cf3 100644 --- a/R/subset.R +++ b/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 diff --git a/man/mdim.Rd b/man/mdim.Rd index dfe3e3f4..78876df0 100644 --- a/man/mdim.Rd +++ b/man/mdim.Rd @@ -42,13 +42,13 @@ write_mdim( \item{raster}{names of the raster variables (default: first two dimensions)} -\item{offset}{integer; offset for each dimension (pixels) of sub-array to read, defaults to 0 for each dimension(requires sf >= 1.0-9)} +\item{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)} \item{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)} \item{step}{integer; step size for each dimension (pixels) of sub-array to read; defaults to 1 for each dimension (requires sf >= 1.0-9)} -\item{proxy}{logical; return proxy object? (not functional yet)} +\item{proxy}{logical; return proxy object?} \item{debug}{logical; print debug info?} diff --git a/man/st_crop.Rd b/man/st_crop.Rd index 58e9d3f4..c4bb9adc 100644 --- a/man/st_crop.Rd +++ b/man/st_crop.Rd @@ -1,11 +1,14 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/proxy.R, R/subset.R +% Please edit documentation in R/mdim.R, R/proxy.R, R/subset.R \name{st_crop} \alias{st_crop} +\alias{st_crop.mdim} \alias{st_crop.stars_proxy} \alias{st_crop.stars} \title{crop a stars object} \usage{ +\method{st_crop}{mdim}(x, y, ...) + \method{st_crop}{stars_proxy}( x, y,