Skip to content

Commit

Permalink
Merge pull request #165 from pascal-sauer/master
Browse files Browse the repository at this point in the history
ncdf4-based nc writing (+ extend helper)
  • Loading branch information
pascal-sauer authored Feb 28, 2024
2 parents 3c385d6 + 635295c commit c8ad8a6
Show file tree
Hide file tree
Showing 19 changed files with 338 additions and 74 deletions.
3 changes: 2 additions & 1 deletion .buildlibrary
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
ValidationKey: '121045680'
ValidationKey: '121449200'
AcceptedWarnings:
- 'Warning: package ''.*'' was built under R version'
- 'Warning: namespace ''.*'' is not available and has been replaced'
- is.na\(\) applied to non-\(list or vector\) of type 'closure'
AcceptedNotes: unable to verify current time
AutocreateReadme: yes
allowLinterWarnings: no
enforceVersionUpdate: no
5 changes: 5 additions & 0 deletions .github/workflows/check.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,11 @@ jobs:
shell: Rscript {0}
run: lucode2:::validkey(stopIfInvalid = TRUE)

- name: Verify that lucode2::buildLibrary was successful
if: github.event_name == 'pull_request'
shell: Rscript {0}
run: lucode2:::isVersionUpdated()

- name: Checks
shell: Rscript {0}
run: |
Expand Down
2 changes: 1 addition & 1 deletion .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ repos:
- id: mixed-line-ending

- repo: https://github.com/lorenzwalthert/precommit
rev: v0.3.2.9025
rev: v0.4.0
hooks:
- id: parsable-R
- id: deps-in-desc
Expand Down
4 changes: 2 additions & 2 deletions CITATION.cff
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@ cff-version: 1.2.0
message: If you use this software, please cite it using the metadata from this file.
type: software
title: 'magclass: Data Class and Tools for Handling Spatial-Temporal Data'
version: 6.13.2
date-released: '2024-01-18'
version: 6.14.0
date-released: '2024-02-27'
abstract: Data class for increased interoperability working with spatial-temporal
data together with corresponding functions and methods (conversions, basic calculations
and basic data manipulation). The class distinguishes between spatial, temporal
Expand Down
6 changes: 3 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Type: Package
Package: magclass
Title: Data Class and Tools for Handling Spatial-Temporal Data
Version: 6.13.2
Date: 2024-01-18
Version: 6.14.0
Date: 2024-02-27
Authors@R: c(
person("Jan Philipp", "Dietrich", , "dietrich@pik-potsdam.de", role = c("aut", "cre")),
person("Benjamin Leon", "Bodirsky", , "bodirsky@pik-potsdam.de", role = "aut"),
Expand Down Expand Up @@ -62,4 +62,4 @@ VignetteBuilder:
knitr
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.3.0
RoxygenNote: 7.3.1
2 changes: 2 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ export(dimExists)
export(dimOrder)
export(dimReduce)
export(dimSums)
export(extend)
export(fulldim)
export(getCPR)
export(getCells)
Expand All @@ -46,6 +47,7 @@ export(getRegionList)
export(getRegions)
export(getSets)
export(getYears)
export(guessResolution)
export(hasCoords)
export(hasSets)
export(is.magpie)
Expand Down
11 changes: 1 addition & 10 deletions R/as.RasterBrick.R
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,6 @@ as.RasterBrick <- function(x, res = NULL) { # nolint
if (!is.magpie(x)) stop("Input is not a magpie object")
if (!requireNamespace("raster", quietly = TRUE)) stop("The package \"raster\" is required!")

.guessRes <- function(xy) {
.tmp <- function(x) {
return(suppressWarnings(min(diff(sort(unique(x))))))
}
guess <- min(.tmp(xy[[1]]), .tmp(xy[[2]]))
if (is.infinite(guess)) guess <- 0.5
return(guess)
}

if (!hasCoords(x) && dimExists(1.2, x)) {
items <- getItems(x, dim = 1.2)
if (length(items) == 59199 && all(items == seq_len(59199))) {
Expand All @@ -44,7 +35,7 @@ as.RasterBrick <- function(x, res = NULL) { # nolint
}
}
xy <- getCoords(x)
if (is.null(res)) res <- .guessRes(xy)
if (is.null(res)) res <- guessResolution(xy)
out <- raster::brick(ncols = 360 / res, nrows = 180 / res, nl = nyears(x) * ndata(x))
m <- wrap(as.array(x), list(1, 2:3), sep = "..")
layerNames <- colnames(m)
Expand Down
13 changes: 1 addition & 12 deletions R/as.SpatRaster.R
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,6 @@ as.SpatRaster <- function(x, res = NULL) { # nolint: object_name_linter.
if (!is.magpie(x)) stop("Input is not a magpie object")
if (!requireNamespace("terra", quietly = TRUE)) stop("The package \"terra\" is required!")

.guessRes <- function(xy) {
.tmp <- function(x) {
return(min(diff(sort(unique(x)))))
}
guess <- min(.tmp(xy[[1]]), .tmp(xy[[2]]))
# use 0.5deg as guess if it cannot be determined otherwise as this is
# the default spatial resolution in the magpie universe.
if (is.infinite(guess)) guess <- 0.5
return(guess)
}

if (!hasCoords(x) && dimExists(1.2, x)) {
items <- getItems(x, dim = 1.2)
if (length(items) == 59199 && all(items == seq_len(59199))) {
Expand All @@ -46,7 +35,7 @@ as.SpatRaster <- function(x, res = NULL) { # nolint: object_name_linter.
}
}
xy <- getCoords(x)
if (is.null(res)) res <- .guessRes(xy)
if (is.null(res)) res <- guessResolution(xy)

m <- wrap(as.array(x), list(1, 2:3), sep = "..")
xyz <- cbind(xy, m)
Expand Down
59 changes: 59 additions & 0 deletions R/extend.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
#' extend
#'
#' Extend a magpie object to a dense grid based on the given xRange, yRange and resolution.
#' This is e.g. required when writing netCDF files. Extending a sparse magpie object to a dense grid
#' requires much more memory, so use with caution.
#'
#' @param x A magpie object
#' @param xRange Range of x coordinates, in ascending order when writing to netCDF
#' @param yRange Range of y coordinates, in descending order when writing to netCDF
#' @param res Resolution of the data, if not provided it will be guessed
#' @return An extended magpie object with the same data as x and gaps filled with NA
#' @author Pascal Sauer
#' @export
extend <- function(x,
xRange = c(-179.75, 179.75),
yRange = c(89.75, -89.75),
res = NULL) {
stopifnot(length(xRange) == 2, length(yRange) == 2)
if (is.null(res)) {
res <- guessResolution(x)
}
coords <- expand.grid(x = seq(xRange[1], xRange[2], if (xRange[1] < xRange[2]) res else -res),
y = seq(yRange[1], yRange[2], if (yRange[1] < yRange[2]) res else -res))
coords <- paste0(coords$x, "|", coords$y)
coords <- gsub("\\.", "p", coords)
coords <- sub("\\|", ".", coords)

sparseCoords <- paste0(getItems(x, "x", full = TRUE),
".",
getItems(x, "y", full = TRUE))
if (!all(sparseCoords %in% coords)) {
stop("The coordinates of the input object are not a subset of the extended coordinates. ",
"Try changing res, xRange or yRange.")
}

subdims1 <- getSets(x)[grep("^d1\\.", names(getSets(x)))]
notCoords <- subdims1[!subdims1 %in% c("x", "y")]

dim1 <- paste0(coords, paste(rep(".NA", length(notCoords)), collapse = ""))

extended <- new.magpie(dim1, getItems(x, 2), getItems(x, 3))

# dimOrder subdims x and y to where they are in x
perm <- match(getSets(x), c("x", "y"))[seq_along(subdims1)]
perm[is.na(perm)] <- seq_along(notCoords) + 2
extended <- dimOrder(extended, perm, dim = 1)

getSets(extended) <- getSets(x)

posInCoords <- match(sparseCoords, coords)

# fill subdims for dim 1 other than x and y (e.g. country code) if they were set in x
dimnames(extended)[[1]][posInCoords] <- getItems(x, 1)

# fill extended with data from x where available
extended[posInCoords, , ] <- x

return(extended)
}
25 changes: 25 additions & 0 deletions R/guessResolution.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
#' guessResolution
#'
#' Guess the resolution of the given magpie object/coordinates by looking at the minimum
#' difference between unique sorted values. Fall back to 0.5 if guess is infinite.
#'
#' @param x A magpie object or the coordinates of a magpie object (the result
#' of \code{\link{getCoords}})
#' @return The guessed resolution of the data
#' @author Jan Philipp Dietrich, Pascal Sauer
#' @export
guessResolution <- function(x) {
if (is.magpie(x)) {
xy <- getCoords(x)
} else {
xy <- x
}
.tmp <- function(x) {
return(min(diff(sort(unique(x)))))
}
guess <- min(.tmp(xy[[1]]), .tmp(xy[[2]]))
# use 0.5deg as guess if it cannot be determined otherwise as this is
# the default spatial resolution in the magpie universe.
if (is.infinite(guess)) guess <- 0.5
return(guess)
}
55 changes: 28 additions & 27 deletions R/write.magpie.R
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,12 @@
#' structure filename_year_datacolumn.asc
#' \item "tif" is the GEOtiff format for gridded data.
#' \item "grd" is the native raster format for gridded data.
#' \item "nc" is the netCDF format for gridded data.
#' \item "nc" is the netCDF format for gridded data. Check ?magclass::writeNC
#' for more details on how nc files are written.
#' }
#'
#' @param x a magclass object. An exception is that formats written via the raster package
#' (currently "nc", "asc", "grd" and "tif") also accept RasterBrick objects which have
#' (currently "asc", "grd" and "tif") also accept RasterBrick objects which have
#' been previously created from a magclass object via as.RasterBrick)
#' @param file_name file name including file ending (wildcards are supported).
#' Optionally also the full path can be specified here (instead of splitting it
Expand All @@ -52,7 +53,8 @@
#' existing data can be combined with the new data using the mbind function
#' @param comment Vector of strings: Optional comment giving additional
#' information about the data. If different to NULL this will overwrite the
#' content of attr(x,"comment")
#' content of attr(x,"comment"). For nc files the unit can be passed
#' with e.g. 'comment = "unit: kg"'.
#' @param comment.char character: a character vector of length one containing a
#' single character or an empty string. Use "" to turn off the interpretation
#' of comments altogether.
Expand All @@ -62,7 +64,8 @@
#' Set to NULL system defaults will be used. Access codes are identical to the
#' codes used in unix function chmod.
#' @param zname name of the time variable for raster files like nc, asc, grd and tif
#' @param ... additional arguments passed to specific write functions
#' @param ... additional arguments passed to specific write functions.
#' Check ?magclass::writeNC for available arguments when writing nc files.
#' @note
#'
#' The binary MAgPIE formats .m and .mz have the following content/structure
Expand Down Expand Up @@ -114,8 +117,8 @@ write.magpie <- function(x, # nolint: object_name_linter, cyclocomp_linter.
if (is.null(file_type)) {
file_type <- tail(strsplit(file_name, "\\.")[[1]], 1) # nolint: object_name_linter.
}
if (inherits(x, "RasterBrick") && !(file_type %in% c("nc", "asc", "grd", "tif"))) {
stop("RasterBrick format is only allowed for file types: nc, asc, grd and tif")
if (inherits(x, "RasterBrick") && !(file_type %in% c("asc", "grd", "tif"))) {
stop("RasterBrick format is only allowed for file types: asc, grd and tif")
}
if (!file_folder == "") {
filePath <- paste(file_folder, file_name, sep = "/")
Expand Down Expand Up @@ -208,29 +211,27 @@ 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("terra", quietly = TRUE)) {
stop("The packages \"ncdf4\" and \"terra\" are required!")
}
if (is.null(getItems(x, 3))) {
getItems(x, 3) <- "Variable"
}
spatRasterDataset <- as.SpatRasterDataset(x)

if (zname != "time") {
warning("zname != 'time' is discouraged, as terra will not recognize it as time dimension")
if (!hasCoords(x) && dimExists(1.2, x)) {
items <- getItems(x, dim = 1.2)
if (length(items) == 59199 && all(items == seq_len(59199))) {
# special treatment for data with 59199 cells as this
# is the originally used 0.5deg resolution in magclass
# before it had been generalized
getCoords(x) <- magclassdata$half_deg[c("lon", "lat")]
}
}
# terra::writeCDF does not set the "axis" attribute for the time dimension, which triggers a warning
suppressSpecificWarnings({
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)
# 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")
if (is.null(comment)) {
unit <- ""
} else {
indicators <- sub(":.*$", "", comment)
units <- sub("^.*: ", "", comment)
if (any(grepl("unit", indicators))) {
unit <- units[grep("unit", indicators)]
} else {
unit <- ""
}
}
ncdf4::nc_close(nc)
writeNC(x = x, filename = file_name, unit = unit, zname = zname, ...)
} else if (file_type == "rds") {
saveRDS(object = x, file = filePath, ...)
} else if (file_type == "cs3" || file_type == "cs3r") {
Expand Down
57 changes: 57 additions & 0 deletions R/writeNC.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
#' Write a magpie object to a netCDF file
#'
#' @param x A magpie object
#' @param filename Name of the netCDF file to write
#' @param unit Unit of the data, to omit pass "" (empty string)
#' @param ... For future expansion.
#' @param compression Level of compression to use (1-9), NA for no compression
#' @param missval The value that encodes NA in the resulting netCDF file
#' @param res Resolution of the data, if not provided it will be guessed
#' @param zname Name of the z dimension in the netCDF file
#' @author Pascal Sauer
writeNC <- function(x, filename, unit, ..., compression = 2, missval = NA,
res = NULL, zname = "time") {
if (!requireNamespace("ncdf4", quietly = TRUE)) {
stop("The ncdf4 package is required to write netCDF files, please install it.")
}
# fail immediately if arguments are not set
stopifnot(is.character(filename), is.character(unit))
if (!(...length() == 0 || all(...names() == "verbose"))) {
stop("Unknown argument passed to writeNC: ", paste(...names(), collapse = ", "))
}

# magclass objects are sparse, fill gaps with NA
coords <- getCoords(x)
x <- extend(x, xRange = c(min(coords$x), max(coords$x)), yRange = c(max(coords$y), min(coords$y)), res = res)
coords <- getCoords(x)

if (zname != "time") {
message("terra will not recognize zname != 'time' as time dimension")
}

if (is.null(getItems(x, 3))) {
getItems(x, 3) <- sub("\\.nc$", "", basename(filename))
}

# create netCDF file
dimVars <- list(ncdf4::ncdim_def("lon", "degrees_east", unique(coords$x)),
ncdf4::ncdim_def("lat", "degrees_north", unique(coords$y)))
if (!is.null(getItems(x, 2))) {
dimVars <- c(dimVars, list(ncdf4::ncdim_def(zname, "years since 0", getYears(x, as.integer = TRUE), unlim = TRUE)))
}
vars <- lapply(getItems(x, 3), function(vname) {
return(ncdf4::ncvar_def(vname, units = unit, dim = dimVars,
missval = missval, compression = compression))
})

nc <- ncdf4::nc_create(filename, vars = vars)
withr::defer(ncdf4::nc_close(nc))

if (!is.null(getItems(x, 2)) && zname == "time") {
ncdf4::ncatt_put(nc, "time", "axis", "T")
}

for (vname in getItems(x, 3)) {
ncdf4::ncvar_put(nc, vname, x[, , vname])
}
}
8 changes: 4 additions & 4 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Data Class and Tools for Handling Spatial-Temporal Data

R package **magclass**, version **6.13.2**
R package **magclass**, version **6.14.0**

[![CRAN status](https://www.r-pkg.org/badges/version/magclass)](https://cran.r-project.org/package=magclass) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1158580.svg)](https://doi.org/10.5281/zenodo.1158580) [![R build status](https://github.com/pik-piam/magclass/workflows/check/badge.svg)](https://github.com/pik-piam/magclass/actions) [![codecov](https://codecov.io/gh/pik-piam/magclass/branch/master/graph/badge.svg)](https://app.codecov.io/gh/pik-piam/magclass) [![r-universe](https://pik-piam.r-universe.dev/badges/magclass)](https://pik-piam.r-universe.dev/builds)

Expand Down Expand Up @@ -56,7 +56,7 @@ In case of questions / problems please contact Jan Philipp Dietrich <dietrich@pi

To cite package **magclass** in publications use:

Dietrich J, Bodirsky B, Bonsch M, Humpenoeder F, Bi S, Karstens K, Leip D, Sauer P (2024). _magclass: Data Class and Tools for Handling Spatial-Temporal Data_. doi: 10.5281/zenodo.1158580 (URL: https://doi.org/10.5281/zenodo.1158580), R package version 6.13.2, <URL: https://github.com/pik-piam/magclass>.
Dietrich J, Bodirsky B, Bonsch M, Humpenoeder F, Bi S, Karstens K, Leip D, Sauer P (2024). _magclass: Data Class and Tools for Handling Spatial-Temporal Data_. doi:10.5281/zenodo.1158580 <https://doi.org/10.5281/zenodo.1158580>, R package version 6.14.0, <https://github.com/pik-piam/magclass>.

A BibTeX entry for LaTeX users is

Expand All @@ -65,8 +65,8 @@ A BibTeX entry for LaTeX users is
title = {magclass: Data Class and Tools for Handling Spatial-Temporal Data},
author = {Jan Philipp Dietrich and Benjamin Leon Bodirsky and Markus Bonsch and Florian Humpenoeder and Stephen Bi and Kristine Karstens and Debbora Leip and Pascal Sauer},
year = {2024},
note = {R package version 6.13.2},
doi = {10.5281/zenodo.1158580},
note = {R package version 6.14.0},
url = {https://github.com/pik-piam/magclass},
doi = {10.5281/zenodo.1158580},
}
```
Loading

0 comments on commit c8ad8a6

Please sign in to comment.