ncdf4-based nc writing for disaggregation.R to avoid memory spike

CHANGELOG.md

@@ -32,6 +32,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **58_peatland** removed realization "on"
 
 ### fixed
+- **scripts** fixed memory spike leading to crashes in disaggregation.R
 - **scripts** Fixed writing of NetCDF files in output/reportMAgPIE2SEALS.R
 - **scripts** Fixed disaggregation.R and disaggregation_LUH2.R to be used with 67k
 - **scripts** bugfix highres.R for bioenergy demand and GHG prices in coupled runs

scripts/output/extra/disaggregation.R

@@ -64,9 +64,49 @@ if (length(map_file) > 1) {
   return(x)
 }
 
+.ncdf4Write <- function(x, filename, unit = "", compression = 2, missval = -9999) {
+  # magclass objects are sparse, fill gaps with NA
+  xy <- getCoords(x)
+  lonCoords <- seq(min(xy$x), max(xy$x), 0.5)
+  latCoords <- seq(max(xy$y), min(xy$y), -0.5)
+
+  coords <- expand.grid(x = lonCoords, y = latCoords)
+  coords <- paste0(coords$x, "|", coords$y)
+  coords <- gsub("\\.", "p", coords)
+  coords <- sub("\\|", ".", coords)
+
+  getItems(x, 1.3) <- NULL
+
+  extended <- new.magpie(coords, getItems(x, 2), getItems(x, 3))
+  extended[getItems(x, 1), , ] <- x
+  x <- extended
+
+  # create netCDF file
+  lonLatTime <- list(ncdf4::ncdim_def("lon", "degrees_east", lonCoords),
+                     ncdf4::ncdim_def("lat", "degrees_north", latCoords),
+                     ncdf4::ncdim_def("time", "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 = lonLatTime,
+                            missval = missval, compression = compression))
+  })
+
+  nc <- ncdf4::nc_create(filename, vars = vars)
+  withr::defer(ncdf4::nc_close(nc))
+
+  ncdf4::ncatt_put(nc, "time", "axis", "T")
+
+  for (vname in getItems(x, 3)) {
+    ncdf4::ncvar_put(nc, vname, x[, , vname])
+  }
+}
+
 .writeDisagg <- function(x, file, comment, message) {
-  write.magpie(.fixCoords(x), file, comment = comment)
-  write.magpie(.fixCoords(x), sub(".mz", ".nc", file), comment = comment, verbose = FALSE)
+  base::message(message)
+  x <- .fixCoords(x)
+  write.magpie(x, file, comment = comment)
+
+  unit <- if (grepl("^unit: ", comment)) sub("^unit: ", "", comment) else ""
+  .ncdf4Write(x, sub(".mz", ".nc", file), unit = unit)
 }
 
 .dissagcrop <- function(gdx, land_hr, map_file) {
@@ -344,6 +384,9 @@ gc()
 land_split_hr <- land_hr[, getYears(area_shr_hr), ]
 area_hr <- area_shr_hr * dimSums(land_split_hr, dim = 3)
 
+rm(area_shr_hr)
+gc()
+
 # replace crop in land_hr in with crop_kfo_rf, crop_kfo_ir, crop_kbe_rf
 # and crop_kbe_ir
 kbe <- c("betr", "begr")
@@ -372,6 +415,8 @@ land_split_hr <- land_split_hr[, , "crop", invert = TRUE]
 # combine land_split_hr with crop_hr.
 land_split_hr <- mbind(crop_hr, fallow, land_split_hr)
 
+rm(crop_kfo_rf, crop_kfo_ir, crop_kbe_rf, crop_kbe_ir, crop_hr, fallow, area_hr)
+
 # split "forestry" into timber plantations, pre-scribed afforestation (NPi/NDC) and endogenous afforestation (CO2 price driven)
 message("Disaggregating forestry")
 farea <- dimSums(landForestry(gdx, level = "cell"), dim = "ac")
@@ -393,6 +438,9 @@ land_split_hr <- land_split_hr[, , "forestry", invert = TRUE]
 # combine land_split_hr with farea_hr
 land_split_hr <- mbind(land_split_hr, farea_hr)
 
+rm(farea, farea_shr, farea_hr)
+gc()
+
 # Write output
 .writeDisagg(land_split_hr, land_hr_split_file,
   comment = "unit: Mha per grid-cell",
@@ -402,6 +450,7 @@ land_split_hr <- mbind(land_split_hr, farea_hr)
   comment = "unit: grid-cell land area fraction",
   message = "Write cropsplit land area share"
 )
+rm(land_split_hr)
 gc()
 
 # --------------------------------
@@ -471,6 +520,9 @@ land_bii_hr <- interpolateAvlCroplandWeighted(
   snv_pol_fader = snv_pol_fader,
   unit = "share"
 )
+
+rm(land_consv_hr, urban_land_hr)
+
 land_bii_hr <- .fixCoords(land_bii_hr)
 
 # Add primary and secondaray other land
@@ -481,12 +533,14 @@ land_bii_hr <- land_bii_hr * side_layers_hr[, , c("forested", "nonforested")]
 
 # Sum over land classes
 bii_hr <- dimSums(land_bii_hr * bii_hr, dim = 3, na.rm = TRUE)
+rm(land_bii_hr)
 
 # Write output
 .writeDisagg(bii_hr, bii_hr_out_file,
   comment = "unitless",
   message = "Write output BII at 0.5°"
 )
+rm(bii_hr)
 gc()
 
 
@@ -524,6 +578,8 @@ out <- peat_hr / dimSums(land_hr[,getYears(peat_hr),], dim = 3)
 out[is.nan(out)] <- 0
 out[is.infinite(out)] <- 0
 
+rm(land_hr, peat_hr)
+
 .writeDisagg(out, peatland_hr_share_out_file,
              comment = "unit: grid-cell land area fraction",
              message = "Write outputs peatland share")

scripts/output/extra/disaggregation_LUH2.R

@@ -490,4 +490,4 @@ gc()
 write.magpie(a, paste0(out_dir, "/LUH2_Yield_Nr.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(a,yield_kr,yield_kr_su)
 gc()
-}
+}