Skip to content

Commit

Permalink
terra based read/write.magpie
Browse files Browse the repository at this point in the history
  • Loading branch information
pfuehrlich-pik committed Jan 11, 2024
1 parent 19fc7d0 commit 2010bb3
Show file tree
Hide file tree
Showing 5 changed files with 56 additions and 81 deletions.
45 changes: 22 additions & 23 deletions R/read.magpie.R
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,8 @@
#' @author Jan Philipp Dietrich, Stephen Bi, Florian Humpenoeder, Pascal Sauer
#' @seealso \code{"\linkS4class{magpie}"}, \code{\link{write.magpie}}
#' @export
read.magpie <- function(file_name, file_folder = "", file_type = NULL, as.array = FALSE, # nolint: object_name_linter.
comment.char = "*", check.names = FALSE, ...) { # nolint: object_name_linter.
read.magpie <- function(file_name, file_folder = "", file_type = NULL, # nolint: object_name_linter, cyclocomp_linter.
as.array = FALSE, comment.char = "*", check.names = FALSE, ...) { # nolint: object_name_linter.

.buildFileName <- function(fileName, fileFolder) {
fileName <- paste0(fileFolder, fileName)
Expand Down Expand Up @@ -111,31 +111,30 @@ read.magpie <- function(file_name, file_folder = "", file_type = NULL, as.array
data = grep(".data", m$metadata$dimtype, fixed = TRUE))
attr(readMagpie, "comment") <- m$comment
} else if (fileType %in% c("asc", "nc", "grd", "tif")) {
if (!requireNamespace("raster", quietly = TRUE)) stop("The package \"raster\" is required!")
if (fileType == "nc") {
if (!requireNamespace("ncdf4", quietly = TRUE)) {
stop("The package \"ncdf4\" is required!")
if (!requireNamespace("terra", quietly = TRUE)) {
stop("The package \"terra\" is required!")
}
nc <- ncdf4::nc_open(fileName)
var <- names(nc[["var"]])
vdim <- vapply(nc[["var"]], function(x) return(x$ndims), integer(1))
var <- var[vdim > 0]
ncdf4::nc_close(nc)
tmp <- list()
for (v in var) {
suppressSpecificWarnings({
warning <- utils::capture.output(tmp[[v]] <- raster::brick(fileName, varname = v, ...))
}, "partial match of 'group' to 'groups'", fixed = TRUE)
if (length(warning) > 0) {
tmp[[v]] <- NULL
next
}
name <- sub("^X([0-9]*)$", "y\\1", names(tmp[[v]]), perl = TRUE)
if (length(name) == 1 && name == "layer") name <- "y0"
names(tmp[[v]]) <- paste0(name, "..", v)

x <- terra::rast(fileName)
if (all(grepl("Time=[0-9]+", names(x)))) {
names(x) <- sub("(.+)_Time=([0-9]+)", "y\\2..\\1", names(x))
} else if (all(grepl("_", names(x)))) {
names(x) <- vapply(names(x), function(n) {
parts <- strsplit(n, "_")[[1]] # e.g. "AFR_3" where 3 means the third entry in terra::time(x)
year <- terra::time(x)[as.integer(parts[2])]
if (is.na(year)) {
year <- as.integer(parts[2])
}
return(paste0("y", year, "..", parts[1]))
}, character(1))
}
readMagpie <- as.magpie(raster::stack(tmp))

readMagpie <- clean_magpie(as.magpie(x))
} else {
if (!requireNamespace("raster", quietly = TRUE)) {
stop("The package \"raster\" is required!")
}
readMagpie <- as.magpie(raster::brick(fileName, ...))
}
} else {
Expand Down
22 changes: 15 additions & 7 deletions R/spatRasterToDataset.R
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,21 @@ spatRasterToDataset <- function(x) {
stop("The package \"terra\" is required!")
}

stopifnot(grepl("^y[0-9]+\\.\\.", names(x)))
varnames <- unique(sub("^y[0-9]+\\.\\.", "", names(x)))
datasets <- lapply(varnames, function(varname) {
spatRaster <- x[paste0("\\.\\.", varname, "$")]
terra::varnames(spatRaster) <- varname
return(spatRaster)
})
if (all(grepl("^y[0-9]+\\.\\.", names(x)))) {
varnames <- unique(sub("^y[0-9]+\\.\\.", "", names(x)))
datasets <- lapply(varnames, function(varname) {
spatRaster <- x[paste0("\\.\\.", varname, "$")]
terra::varnames(spatRaster) <- varname
return(spatRaster)
})
} else {
varnames <- names(x)
datasets <- lapply(varnames, function(varname) {
spatRaster <- x[varname]
terra::varnames(spatRaster) <- varname
return(spatRaster)
})
}
x <- terra::sds(datasets)
names(x) <- varnames
return(x)
Expand Down
64 changes: 16 additions & 48 deletions R/write.magpie.R
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
#' full access for user, read access for group and no acess for anybody else).
#' Set to NULL system defaults will be used. Access codes are identical to the
#' codes used in unix function chmod.
#' @param zname Time variable for the writeRaster function
#' @param zname name of the time variable for raster files like nc, asc, grd and tif
#' @param ... additional arguments passed to specific write functions
#' @note
#'
Expand Down Expand Up @@ -98,7 +98,7 @@ write.magpie <- function(x, # nolint: object_name_linter, cyclocomp_linter.
file_name, file_folder = "", file_type = NULL, # nolint: object_name_linter.
append = FALSE, comment = NULL,
comment.char = "*", # nolint: object_name_linter.
mode = NULL, zname = "Time", ...) {
mode = NULL, zname = "time", ...) {
umask <- Sys.umask()
if (!is.null(mode)) {
umaskMode <- as.character(777 - as.integer(mode))
Expand Down Expand Up @@ -208,56 +208,24 @@ write.magpie <- function(x, # nolint: object_name_linter, cyclocomp_linter.
raster::writeRaster(x, filename = filePath, format = format[file_type], overwrite = TRUE,
zname = zname, zunit = zunit, varname = varname, ...)
} else if (file_type == "nc") {

if (!requireNamespace("ncdf4", quietly = TRUE) || !requireNamespace("raster", quietly = TRUE)) {
stop("The packages \"ncdf4\" and \"raster\" are required!")
}
.sub <- function(rx, name) {
layer <- sub("^.*\\.\\.", "", names(rx))
if (length(unique(layer)) == 1) return(rx)
return(rx[[which(layer == name)]])
}
if (is.magpie(x)) {
varnames <- getItems(x, dim = 3)
zunit <- ifelse(all(isYear(getYears(x))), "years", "")
years <- getYears(x, as.integer = TRUE)
x <- as.RasterBrick(x)
} else if (inherits(x, "RasterBrick")) {
tmp <- names(x)
tmp <- strsplit(tmp, "\\..")
years <- sort(unique(unlist(lapply(tmp, function(x) x[1]))))
varnames <- sort(unique(unlist(lapply(tmp, function(x) x[2]))))
zunit <- ifelse(all(isYear(years)), "years", "")
years <- as.numeric(gsub("y", "", years))
if (!requireNamespace("ncdf4", quietly = TRUE) || !requireNamespace("terra", quietly = TRUE)) {
stop("The packages \"ncdf4\" and \"terra\" are required!")
}
if (is.null(varnames)) varnames <- "Variable"
if (is.null(comment)) {
unit <- "not specified"
} else {
indicators <- sub(":.*$", "", comment)
units <- sub("^.*: ", "", comment)
if (any(grepl("unit", indicators))) {
unit <- units[grep("unit", indicators)]
} else {
unit <- "not specified"
}
spatRasterDataset <- as.SpatRasterDataset(x)

if (zname != "time") {
warning("zname != 'time' is discouraged, as terra will not recognize it as time dimension")
}
# raster is using partial matching resulting in a warning if warnPartialMatchDollar is set
# terra::writeCDF does not set the "axis" attribute for the time dimension, which triggers a warning
suppressSpecificWarnings({
raster::writeRaster(.sub(x, varnames[1]), filename = filePath, format = "CDF", overwrite = TRUE,
compression = 9, zname = zname, zunit = zunit, varname = varnames[1], varunit = unit, ...)
}, "partial match of 'group' to 'groups'", fixed = TRUE)
terra::writeCDF(spatRasterDataset, filePath, overwrite = TRUE, zname = zname, ...)
}, paste0("GDAL Message 1: dimension #0 (", zname, ") is not a Time or Vertical dimension."), fixed = TRUE)

# set the "axis" attribute to "T" for the time dimension to prevent further warnings when reading the file
nc <- ncdf4::nc_open(filePath, write = TRUE)
if (zunit == "years") {
try(ncdf4::ncvar_put(nc, zname, years), silent = TRUE)
}
if (length(varnames) > 1) {
for (i in varnames[-1]) {
suppressSpecificWarnings({
nc <- ncdf4::ncvar_add(nc, ncdf4::ncvar_def(i, unit, nc$dim, compression = 9))
}, "partial match of 'group' to 'groups'", fixed = TRUE)
ncdf4::ncvar_put(nc, i, aperm(as.array(.sub(x, i)), c(2, 1, 3)))
}
# nc will not have time dim when x had no time dimension
if (zname %in% names(nc$dim)) {
ncdf4::ncatt_put(nc, zname, "axis", "T")
}
ncdf4::nc_close(nc)
} else if (file_type == "rds") {
Expand Down
4 changes: 2 additions & 2 deletions man/write.magpie.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion tests/testthat/test-raster.R
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
skip_if_not_installed("terra")

test_that("terra convertion does not alter data", {
test_that("terra conversion does not alter data", {
expect_error(as.SpatRaster(1), "not a magpie object")
for (i in c(0.5, 2)) {
for (j in c(1, 4)) {
Expand Down

0 comments on commit 2010bb3

Please sign in to comment.