#' @name st_as_stars
#' @param att see \link[raster:factor]{factorValues}; column in the RasterLayer's attribute table
#' @param ignore_file logical; if \code{TRUE}, ignore the Raster object file name
#' @export
st_as_stars.Raster = function(.x, ..., att = 1, ignore_file = FALSE) {
if (!requireNamespace("raster", quietly = TRUE))
stop("package raster required, please install it first") # nocov
#0 360 -90 90
e = as.vector(raster::extent(.x)) # xmin xmax ymin ymax
RasterIO = if ("data" %in% slotNames(.x) && inherits(.x@data, ".SingleLayerData"))
list(bands = .x@data@band)
else
list()
if (!ignore_file) {
file = if ("file" %in% slotNames(.x))
.x@file@name
else if ("filename" %in% slotNames(.x))
.x@filename
else
""
if (file != "") {
r = try(read_stars(file, RasterIO = RasterIO, ...), silent = TRUE)
if (!inherits(r, "try-error")) {
if (is.na(st_crs(r)))
r = st_set_crs(r, st_crs(raster::crs(.x)))
r = fix_dims(r, e)
if (inherits(.x, "RasterLayer"))
names(r) <- names(.x)
return(r)
}
}
}
if (!requireNamespace("sp", quietly = TRUE))
stop("package sp required, please install it first") # nocov
v = raster::values(.x)
dim(v) = dim(.x)[c(2,1,3)]
if (all(raster::is.factor(.x))) {
l = raster::levels(.x)[[1]]$levels
if (length(l) == 0) # get the layer's RAT, column att:
l = raster::factorValues(.x, seq_len(max(v, na.rm = TRUE)), att = att)[[1]]
colors = try(.x@legend@colortable, silent = TRUE)
if (inherits(colors, "try-error") || length(colors) == 0)
colors = NULL
else
colors = colors[-1]
v = structure(v, class = "factor", levels = as.character(l), colors = colors)
# FIXME: should we handle levels for all layers here, or break on multiple different ones?
}
dimensions = list(
x = create_dimension(from = 1, to = dim(v)[1], offset = e[1],
delta = (e[2]-e[1])/dim(v)[1], refsys = st_crs(raster::crs(.x))),
y = create_dimension(from = 1, to = dim(v)[2], offset = e[4],
delta = (e[3]-e[4])/dim(v)[2], refsys = st_crs(raster::crs(.x))))
z = raster::getZ(.x)
dimensions$band = if (is.null(z))
create_dimension(values = names(.x))
else
create_dimension(values = z)
l = if (length(names) > 1)
setNames(list(v), deparse(substitute(.x), 50))
else
setNames(list(v), names(.x)[1])
ret = st_as_stars(l, dimensions = create_dimensions(dimensions, get_raster()))
if (dim(ret)[3] == 1)
adrop(ret, 3)
else
ret
}
#' Coerce stars object into a Raster raster or brick
#'
#' Coerce stars object into a Raster raster or brick
#' @param from object to coerce
#' @name as
#' @rdname coerce-methods
#' @aliases coerce,stars,Raster-method
#' @aliases coerce,stars_proxy,Raster-method
#' @returns RasterLayer or RasterBrick
#' @details If the stars object has more than three dimensions, all dimensions higher than the third will be collapsed into the third dimensions. If the stars object has only an x/y raster but multiple attributes, these are merged first, then put in a raster brick.
setAs("stars", "Raster", function(from) {
if (!requireNamespace("sp", quietly = TRUE))
stop("package sp required, please install it first") # nocov
if (!requireNamespace("raster", quietly = TRUE))
stop("package raster required, please install it first") # nocov
if (!is_regular_grid(from))
stop("only regular rasters can be converted to Raster* objects")
st_as_raster(from, class = "Raster")
})
setAs("stars_proxy", "Raster", function(from) {
if (!requireNamespace("raster", quietly = TRUE))
stop("package raster required, please install it first") # nocov
if (!is_regular_grid(from))
stop("only regular rasters can be converted to Raster* objects")
if (length(attr(from, "call_list"))) {
fname = paste0(tempfile(), ".tif")
write_stars(from, fname)
from = fname
}
raster::brick(unlist(from))
})
get_terra_levels = function(x, colors) {
# create factor levels, as used by stars, from SpatRaster levels in a data.frame
# see https://github.com/r-spatial/stars/pull/484
x = x[order(x[[1]]), ] # sort table on level
missing_labels = x[[2]] == ""
levels = x[[1]]
if (any(levels < 0))
stop("negative IDs in SpatRaster levels not supported")
ex = setdiff(0:max(levels), levels)
exclude = rep(FALSE, max(levels) + 1)
exclude[ex + 1] = TRUE # 0-based vector
list(levels = levels[!missing_labels],
labels = x[[2]][!missing_labels],
exclude = exclude,
colors = colors[!missing_labels])
}
#' @name st_as_stars
#' @param ignore_file logical; if \code{TRUE}, ignore the SpatRaster object file name
#' @param as_attributes logical; if \code{TRUE} and \code{.x} has more than one layer, load these as separate attributes rather than as a band or time dimension (only implemented for the case where \code{ignore_file} is \code{TRUE})
#' @export
st_as_stars.SpatRaster = function(.x, ..., ignore_file = FALSE,
as_attributes = all(terra::is.factor(.x))) {
if (!requireNamespace("terra", quietly = TRUE))
stop("package terra required, please install it first") # nocov
#0 360 -90 90
e = as.vector(terra::ext(.x)) # xmin xmax ymin ymax
src = terra::sources(.x, bands=TRUE)
attr_name = basename(src$source[1])
if (!ignore_file && all(src$source != "")) {
# there can be multiple files, but only the first one is used here.
# perhaps a warning should be given; better would be to iterate over "sid"
# but you might have a situation where some sources are filenames and others are not
if (attr_name == "")
attr_name = "values"
lst = vector("list", length(unique(src$sid)))
for (i in unique(src$sid)) {
file = unique(src$source[src$sid == i])
if (length(file) > 1)
stop("more than one file per sid: giving up; try ignore_file=FALSE")
# RasterIO = if (dim(.x)[3] == 1)
# > 1 would be more sensible?
# But this can only be ignored if the _file_ has 1 band
RasterIO = list(bands = src$bands[src$sid == i])
r = try(read_stars(file, RasterIO = RasterIO, ...), silent = TRUE)
if (! inherits(r, "try-error")) {
if (is.na(st_crs(r)) && terra::crs(.x) != "")
r = st_set_crs(r, st_crs(terra::crs(.x)))
r = fix_dims(r, e)
if (length(unique(src$sid)) > 1 && length(dim(r)) > 2)
r = split(r)
#transfer the layer/band names as well?
# ... = names(.x)[1:(dim(r)[3])]
# perhaps check whether they represent a time dimension (all(!is.na(time(.x))))
} else
stop(paste("error reading", file, "bands", paste0(RasterIO$bands, collapse = " ")))
lst[[i]] = r
}
ret = if (length(lst) > 1)
merge(setNames(do.call(c, lst), names(.x)))
else
r
if (!all(is.na(terra::time(.x))))
ret = st_set_dimensions(ret, 3, values = terra::time(.x), names = "time")
} else { # ignore_file TRUE:
if (terra::nlyr(.x) > 1 && as_attributes) {
ret = do.call(c, lapply(seq_len(terra::nlyr(.x)), function(i) st_as_stars(.x[[i]], ignore_file = TRUE)))
if (!is.null(names(.x)))
names(ret) = names(.x)
return(ret) # RETURNS
}
if (attr_name == "") {
if (all(names(.x) == ""))
attr_name = "values"
else
attr_name = paste(names(.x)[1], collapse = ".")
}
v = terra::values(.x, mat = FALSE)
dimv = dim(v) = setNames(dim(.x)[c(2,1,3)], c("x", "y", "band"))
if (all(terra::is.factor(.x))) {
if (length(terra::levels(.x)) > 1)
warning("ignoring categories/levels for all but first layer")
colors = try(rgb(terra::coltab(.x)[[1]][-1], maxColorValue = 255), silent = TRUE)
if (inherits(colors, "try-error") || length(colors) == 0)
colors = NULL
if (inherits(l <- terra::levels(.x)[[1]], "data.frame"))
l = get_terra_levels(l, colors)
else
stop("terra levels should return a list of data.frame's; pls update terra")
v = structure(factor(as.vector(v), levels = l$levels, labels = l$labels),
dim = dimv, colors = l$colors, exclude = l$exclude)
}
crs = if (terra::crs(.x) == "")
NA_crs_
else
st_crs(terra::crs(.x))
dimensions = list(
x = create_dimension(from = 1, to = dim(v)[1], offset = e[1],
delta = (e[2]-e[1])/dim(v)[1], refsys = crs),
y = create_dimension(from = 1, to = dim(v)[2], offset = e[4],
delta = (e[3]-e[4])/dim(v)[2], refsys = crs))
dimensions$band = create_dimension(values = names(.x))
ret = st_as_stars(setNames(list(v), attr_name), dimensions = create_dimensions(dimensions, get_raster()))
if (dim(ret)[3] == 1)
ret = adrop(ret, 3)
else if (!all(is.na(terra::time(.x))))
ret = st_set_dimensions(ret, 3, values = terra::time(.x), names = "time")
}
ret
}
#' Coerce stars object into a terra SpatRaster
#'
#' Coerce stars object into a terra SpatRaster
#' @param from object to coerce
#' @name as
#' @rdname coerce-methods
#' @aliases coerce,stars,Terra-method
#' @aliases coerce,stars_proxy,Terra-method
#' @returns SpatRaster
#' @details If the stars object has more than three dimensions, all dimensions higher than the third will be collapsed into the third dimensions. If the stars object has only an x/y raster but multiple attributes, these are merged first, then put in a SpatRaster.
setAs("stars", "SpatRaster", function(from) {
if (!requireNamespace("terra", quietly = TRUE))
stop("package terra required, please install it first") # nocov
if (!is_regular_grid(from))
stop("only regular rasters can be converted to SpatRaster objects")
st_as_raster(from, class = "SpatRaster")
})
setAs("stars_proxy", "SpatRaster", function(from) {
if (!requireNamespace("terra", quietly = TRUE))
stop("package terra required, please install it first") # nocov
if (!is_regular_grid(from))
stop("only regular rasters can be converted to SpatRaster objects")
if (length(attr(from, "call_list"))) {
fname = paste0(tempfile(), ".tif")
write_stars(from, fname)
from = fname
}
terra::rast(unlist(from))
})
fix_dims = function(r, e){
e = unname(e)
rdims = attr(r, "dimensions")
dxy = attr(rdims, "raster")$dimensions
dimx = rdims[[dxy[1]]]
dimy = rdims[[dxy[2]]]
xrev = dimx$delta < 0
yrev = dimy$delta < 0
if (xrev) {
dimx$offset = e[2]
dimx$delta = (e[1] - e[2]) / (dimx$to - dimx$from + 1)
} else {
dimx$offset = e[1]
dimx$delta = (e[2] - e[1]) / (dimx$to - dimx$from + 1)
}
if (yrev) {
dimy$offset = e[4]
dimy$delta = (e[3] - e[4]) / (dimy$to - dimy$from + 1)
} else {
dimy$offset = e[3]
dimy$delta = (e[4] - e[3]) / (dimy$to - dimy$from + 1)
}
attr(r, "dimensions")[[dxy[1]]] = dimx
attr(r, "dimensions")[[dxy[2]]] = dimy
r
}
st_as_raster = function(x, class, ...) {
stopifnot(inherits(x, "stars"))
if (is.null(names(x)))
names(x) = make.names(seq_along(x)) # FIXME: why? sdsr_exercises
x_crs = st_crs(x)
x = st_upfront(x) # x/y dimensions first
if (length(dim(x)) > 3) {
warning("folding all higher dimensions into the third dimension") # nocov
x = st_apply(x, 1:2, as.vector) # fortunes::fortune("side effect") # nocov
}
if (length(dim(x)) == 2 && length(x) > 1)
x = merge(x)
d = st_dimensions(x)
if (d[[2]]$delta > 0) { # swap:
ny = dim(x)[2]
d[[2]]$offset = d[[2]]$offset + ny * d[[2]]$delta # top
d[[2]]$delta = -d[[2]]$delta # going down
x[[1]] = if (length(dim(x)) == 2)
x[[1]][,ny:1]
else
x[[1]][,ny:1,]
}
dxy = attr(d, "raster")$dimensions
stopifnot(all(dxy %in% names(d)))
bb = st_bbox(x)
levels = NULL
coltab = vector("list", length(x))
values = if (all(sapply(x, is.factor))) {
ex = attr(x[[1]], "exclude")
if (is.null(ex) || class != "SpatRaster")
structure(merge(x)[[1]], dim = NULL) # return the factor
else {
v = vector("list", length(x))
levels = vector("list", length(x))
for (i in seq_along(v)) {
ex = attr(x[[i]], "exclude")
ix = which(!ex) - 1 # 0-based index
n = as.numeric(structure(x[[i]], dim = NULL)) # factor -> numeric
v[[i]] = ix[n]
levels[[i]] = data.frame(IDs = ix, categories = levels(x[[i]]))
if (!is.null(ct <- attr(x[[i]], "colors"))) {
coltab[[i]] = t(col2rgb(rep("#000000", length(ex)), alpha = TRUE))
coltab[[i]][which(!ex),] = t(col2rgb(ct, alpha = TRUE))
}
}
do.call(c, v)
}
} else {
if (any(sapply(x, is.factor)))
warning("mix of factor and non-factor attributes: all factor levels are ignored")
as.vector(merge(x)[[1]])
}
if (class == "SpatRaster") {
third = setdiff(names(d), dxy)
b = terra::rast(nrows = dim(x)[ dxy[2] ], ncols=dim(x)[ dxy[1] ],
xmin = bb[1], xmax = bb[3], ymin = bb[2], ymax = bb[4],
nlyrs = ifelse(length(dim(x)) == 2, 1, dim(x)[third]),
crs = x_crs$wkt)
terra::values(b) = values
if (!is.null(levels)) {
levels(b) = levels
if (!all(sapply(coltab, is.null)))
for (i in seq_len(terra::nlyr(b))) {
m = coltab[[i]]
terra::coltab(b, layer = i) = cbind(value = seq_len(nrow(m)) - 1, m)
}
}
if (length(dim(x)) != 2) {
z = seq(d[[third]])
if (!any(is.na(z))) {
if (is.character(z)) {
names(b) = z
} else {
names(b) = paste0(third, z)
}
}
} else if (length(x) == terra::nlyr(b))
names(b) = names(x)
} else {
if (length(dim(x)) == 2) {
b = raster::raster(nrows=dim(x)[ dxy[2] ], ncols=dim(x)[ dxy[1] ],
xmn = bb[1], xmx = bb[3], ymn = bb[2], ymx = bb[4],
crs = as(x_crs, "CRS"), vals = values)
} else {
third = setdiff(names(d), dxy)
b = raster::brick(nrows=dim(x)[ dxy[2] ], ncols=dim(x)[ dxy[1] ],
xmn = bb[1], xmx = bb[3], ymn = bb[2], ymx = bb[4], nl = dim(x)[third],
crs = as(x_crs, "CRS"))
raster::values(b) = values
z = seq(d[[third]])
if (!any(is.na(z))) {
if (is.character(z))
names(b) = z
else
b = raster::setZ(b, z)
}
}
}
return(b)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.