Merge pull request #168 from pascal-sauer/master

more memory efficient writeNC, allow specifying whole grid (not just resolution)

Author: pascal-sauer

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '121449200'
+ValidationKey: '121782300'
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'
 - 'Warning: namespace ''.*'' is not available and has been replaced'

CITATION.cff

@@ -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.14.0
-date-released: '2024-02-27'
+version: 6.15.0
+date-released: '2024-03-20'
 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

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: magclass
 Title: Data Class and Tools for Handling Spatial-Temporal Data
-Version: 6.14.0
-Date: 2024-02-27
+Version: 6.15.0
+Date: 2024-03-20
 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"),

R/extend.R

@@ -5,32 +5,54 @@
 #' 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
+#' @param gridDefinition A vector of 5 numeric values: c(xMin, xMax, yMin, yMax, resolution).
+#' Use c(-179.75, 179.75, -89.75, 89.75, 0.5) to extend to a standard 0.5-degree-resolution
+#' lon/lat grid. If NULL, use min/max of coordinates in x and guessResolution.
+#' @param crop If TRUE, discard cells from x which are not in the gridDefinition grid. If FALSE,
+#' throw an error if the coordinates of x are not a subset of the extended coordinates.
+#' @return Magpie object x with dense grid according to gridDefinition, 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)
+extend <- function(x, gridDefinition = NULL, crop = FALSE) {
+  if (is.null(gridDefinition)) {
+    coords <- getCoords(x)
+    firstX <- min(coords$x)
+    lastX <- max(coords$x)
+    firstY <- max(coords$y)
+    lastY <- min(coords$y)
+    res <- guessResolution(coords)
+  } else {
+    stopifnot(length(gridDefinition) == 5)
+    firstX <- gridDefinition[1]
+    lastX <- gridDefinition[2]
+    firstY <- gridDefinition[3]
+    lastY <- gridDefinition[4]
+    res <- gridDefinition[5]
   }
-  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 <- expand.grid(x = seq(firstX, lastX, if (firstX < lastX) res else -res),
+                        y = seq(firstY, lastY, if (firstY < lastY) 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.")
+  if (crop) {
+    x <- x[sparseCoords %in% coords, , ]
+    if (length(x) > 0) {
+      sparseCoords <- paste0(getItems(x, "x", full = TRUE),
+                             ".",
+                             getItems(x, "y", full = TRUE))
+    } else {
+      sparseCoords <- character(0)
+    }
+  } else {
+    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)))]

R/writeNC.R

@@ -3,14 +3,17 @@
 #' @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 ... 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 gridDefinition A vector of 5 numeric values: c(xMin, xMax, yMin, yMax, resolution).
+#' Use c(-179.75, 179.75, -89.75, 89.75, 0.5) to write a standard 0.5-degree-resolution
+#' lon/lat grid. If NULL, use min/max of coordinates in x and guessResolution
 #' @param zname Name of the z dimension in the netCDF file
+#' @param progress If TRUE, print progress messages
 #' @author Pascal Sauer
 writeNC <- function(x, filename, unit, ..., compression = 2, missval = NA,
-                    res = NULL, zname = "time") {
+                    gridDefinition = NULL, zname = "time", progress = FALSE) {
   if (!requireNamespace("ncdf4", quietly = TRUE)) {
     stop("The ncdf4 package is required to write netCDF files, please install it.")
   }
@@ -20,10 +23,25 @@ writeNC <- function(x, filename, unit, ..., compression = 2, missval = NA,
     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 (is.null(gridDefinition)) {
+    coords <- getCoords(x)
+    firstX <- min(coords$x)
+    lastX <- max(coords$x)
+    firstY <- max(coords$y)
+    lastY <- min(coords$y)
+    res <- guessResolution(coords)
+  } else {
+    stopifnot(length(gridDefinition) == 5,
+              gridDefinition[1] < gridDefinition[2],
+              gridDefinition[3] < gridDefinition[4])
+    firstX <- gridDefinition[1]
+    lastX <- gridDefinition[2]
+    firstY <- gridDefinition[4]
+    lastY <- gridDefinition[3]
+    res <- gridDefinition[5]
+  }
+  xCoords <- seq(firstX, lastX, res)
+  yCoords <- seq(firstY, lastY, -res)
 
   if (zname != "time") {
     message("terra will not recognize zname != 'time' as time dimension")
@@ -34,10 +52,12 @@ writeNC <- function(x, filename, unit, ..., compression = 2, missval = NA,
   }
 
   # 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)))
+  dimVars <- list(ncdf4::ncdim_def("lon", "degrees_east", xCoords),
+                  ncdf4::ncdim_def("lat", "degrees_north", yCoords))
+  hasTime <- !is.null(getItems(x, 2))
+  if (hasTime) {
+    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,
@@ -47,11 +67,23 @@ writeNC <- function(x, filename, unit, ..., compression = 2, missval = NA,
   nc <- ncdf4::nc_create(filename, vars = vars)
   withr::defer(ncdf4::nc_close(nc))
 
-  if (!is.null(getItems(x, 2)) && zname == "time") {
+  if (hasTime && zname == "time") {
     ncdf4::ncatt_put(nc, "time", "axis", "T")
   }
 
-  for (vname in getItems(x, 3)) {
-    ncdf4::ncvar_put(nc, vname, x[, , vname])
+  messageIf(progress, "Writing ", filename)
+
+  for (i in seq_along(getItems(x, 3))) {
+    vname <- getItems(x, 3)[i]
+    messageIf(progress, i, "/", length(getItems(x, 3)), " Writing ", vname)
+    # memory problems due to the extend happening here? checkout commit
+    # 9e780b5fe00f38b82d20edc245e34605c6eb9c46 to get a chunk-based solution
+    ncdf4::ncvar_put(nc, vname, extend(x[, , vname], gridDefinition = gridDefinition))
+  }
+}
+
+messageIf <- function(condition, ...) {
+  if (condition) {
+    message(...)
   }
 }

README.md

@@ -1,6 +1,6 @@
 # Data Class and Tools for Handling Spatial-Temporal Data
 
-R package **magclass**, version **6.14.0**
+R package **magclass**, version **6.15.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)
 
@@ -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 <https://doi.org/10.5281/zenodo.1158580>, R package version 6.14.0, <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.15.0, <https://github.com/pik-piam/magclass>.
 
 A BibTeX entry for LaTeX users is
 
@@ -65,7 +65,7 @@ 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.14.0},
+  note = {R package version 6.15.0},
   url = {https://github.com/pik-piam/magclass},
   doi = {10.5281/zenodo.1158580},
 }

man/extend.Rd

@@ -2,20 +2,27 @@ test_that("extend works", {
   animal <- maxample("animal")
   # shuffle first dim to ensure order does not matter
   animal <- animal[sample(seq_along(getItems(animal, 1))), , ]
+  animalCoords <- getCoords(animal)
   expect_identical(guessResolution(animal), 0.5)
 
   x <- extend(animal)
-  expectedCoords <- expand.grid(x = seq(-179.75, 179.75, 0.5),
-                                y = seq(89.75, -89.75, -0.5),
+  expectedCoords <- expand.grid(x = seq(min(animalCoords$x), max(animalCoords$x), 0.5),
+                                y = seq(max(animalCoords$y), min(animalCoords$y), -0.5),
                                 KEEP.OUT.ATTRS = FALSE)
   expect_identical(getCoords(x), expectedCoords)
   expect_identical(getSets(x), getSets(animal))
   expect_identical(getItems(x[getItems(animal, 1), , ], 1), getItems(animal, 1))
 
-  expect_error(extend(animal, res = 1),
+  expect_error(extend(animal, gridDefinition = c(3.25, 6.75, 53.25, 49.75, 1)),
                "The coordinates of the input object are not a subset of the extended coordinates")
 
-  x <- extend(animal, xRange = c(-3.25, 16.75), yRange = c(63.25, -59.75), res = 0.25)
+  x <- extend(animal, gridDefinition = c(3.25, 6.75, 53.25, 49.75, 1), crop = TRUE)
+  expectedCoords <- expand.grid(x = seq(3.25, 6.75, 1),
+                                y = seq(53.25, 49.75, -1),
+                                KEEP.OUT.ATTRS = FALSE)
+  expect_identical(getCoords(x), expectedCoords)
+
+  x <- extend(animal, gridDefinition = c(-3.25, 16.75, 63.25, -59.75, 0.25))
 
   expectedCoords <- expand.grid(x = seq(-3.25, 16.75, 0.25),
                                 y = seq(63.25, -59.75, -0.25),
@@ -26,8 +33,8 @@ test_that("extend works", {
   animal2 <- dimOrder(animal, c(3, 4, 2, 1), dim = 1)
   expect_identical(names(dimnames(animal2))[1], "country.cell.y.x")
   x <- extend(animal2)
-  expectedCoords <- expand.grid(x = seq(-179.75, 179.75, 0.5),
-                                y = seq(89.75, -89.75, -0.5),
+  expectedCoords <- expand.grid(x = seq(min(animalCoords$x), max(animalCoords$x), 0.5),
+                                y = seq(max(animalCoords$y), min(animalCoords$y), -0.5),
                                 KEEP.OUT.ATTRS = FALSE)
   expect_identical(getCoords(x), expectedCoords)
 })

man/writeNC.Rd

@@ -78,6 +78,7 @@ test_that("read supports older formats", {
 })
 
 test_that("handling of spatial data works", {
+  skip_if_not_installed("ncdf4")
   td <- withr::local_tempdir()
 
   md <- magclass:::magclassdata$half_deg
@@ -117,6 +118,12 @@ test_that("handling of spatial data works", {
   getItems(anc, dim = 3) <- NULL
   getSets(anc) <- c("x", "y", "d2", "d3")
   expect_identical(anc, a)
+
+  write.magpie(a, file.path(td, "animal2.nc"), gridDefinition = c(3, 7, 49, 54, 0.25))
+  animalRaster <- ncdf4::nc_open(file.path(td, "animal2.nc"))
+  withr::defer(ncdf4::nc_close(animalRaster))
+  expect_identical(as.vector(animalRaster$dim$lon$vals), seq(3, 7, 0.25))
+  expect_identical(as.vector(animalRaster$dim$lat$vals), seq(54, 49, -0.25))
 })
 
 
@@ -162,9 +169,11 @@ test_that("read/write triggers errors and warnings correctly", {
   expect_error(read.magpie(tmpfile), "does not contain a magpie object")
   unlink(tmpfile)
 
-  expect_error(write.magpie(as.magpie(1), "bla.nc"), "No coordinates")
   expect_error(write.magpie(a, "blub.grd"), "no support for multiple variables")
   expect_error(write.magpie(a[, , 1], "blub.asc"), "choose just one")
+
+  skip_if_not_installed("ncdf4")
+  expect_error(write.magpie(as.magpie(1), "bla.nc"), "No coordinates")
 })
 
 test_that("copy.magpie works", {