Merge branch 'magpiemodel:develop' into develop

CHANGELOG.md

@@ -33,8 +33,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **58_peatland** removed realization "on"
 
 ### fixed
-- **scripts** Fixed writing of NetCDF files in output/reportMAgPIE2SEALS.R
-- **scripts** Fixed disaggregation.R and disaggregation_LUH2.R to be used with 67k
+- **scripts** disaggregation_LUH2.R which no longer relies on the old raster-based write.magpie for nc files
+- **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
 - **35_natveg** bugfixes ac_est
 - **35_natveg** removed scaling of pm_carbon_density_ac

DESCRIPTION

@@ -21,7 +21,7 @@ Imports:
     lusweave,
     m4fsdp,
     madrat,
-    magclass,
+    magclass (>= 6.14.0),
     magpie4,
     MagpieNCGains,
     magpiesets,

config/scenario_config.csv

@@ -82,7 +82,7 @@ gms$s73_timber_demand_switch;;;;;;;;;;;;;;;;;;;;;;;1;1;1;1;1;;;;;1;1;0;;;;;;;;;;
 gms$s35_forest_damage;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;4;4;4;4;4;;;;;;;
 gms$c32_shock_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;none;002lin2030;004lin2030;008lin2030;016lin2030;;;;;;;
 gms$c35_shock_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;none;002lin2030;004lin2030;008lin2030;016lin2030;;;;;;;
-input['cellular'];;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.96_h12_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.96_h12_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;rev4.96_h12_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.96_h12_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.96_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.96_h12_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;
+input['cellular'];;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.99_h12_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.99_h12_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;rev4.99_h12_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.99_h12_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.99_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.99_h12_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;
 gms$c52_land_carbon_sink_rcp;;nocc;nocc_hist;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;RCP19;RCP26;RCP45;RCP60;RCPBU;RCPBU;
 gms$c57_macc_version;;;;;;;;;;;;;;;;;;;;;;;;;;;;PBL_2022;PBL_2022;PBL_2022;PBL_2022;;;;;;;;;;;;;;;
 gms$c57_macc_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;;Default;Optimistic;Default;Optimistic;;;;;;;;;;;;;;;

config/scenario_fsec.csv

@@ -70,6 +70,6 @@ gms$c73_build_demand;;;;;;;;;;;;;;;;;;;;;;;;50pc;;;;;;;;;;;;;;;;;;;;;;;;;
 input['cellular'];rev4.99_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.99_FSEC_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.99_FSEC_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;;rev4.99_FSEC_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.99_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.99_FSEC_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.99_FSEC_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;
 input['regional'];rev4.99_FSEC_magpie.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['validation'];rev4.99_FSEC_validation.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-input['additional'];additional_data_rev4.47.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+input['additional'];additional_data_rev4.48.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['calibration'];calibration_FSEC_24Mar23.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 magicc_emis_scen;bjoernAR6_C_SSP2-NDC.mif;;;bjoernAR6_C_SSP2-PkBudg900.mif;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;bjoernAR6_C_SSP1-NDC.mif;;;;;;;;;;;;bjoernAR6_C_SSP1-PkBudg900.mif

scripts/output/extra/disaggregation.R

@@ -65,8 +65,10 @@ if (length(map_file) > 1) {
 }
 
 .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)
+  write.magpie(x, sub(".mz", ".nc", file), comment = comment)
 }
 
 .dissagcrop <- function(gdx, land_hr, map_file) {
@@ -344,6 +346,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 +377,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 +400,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 +412,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 +482,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 +495,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 +540,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

@@ -60,6 +60,50 @@ convertLUH2 <- function(x) {
 
 }
 
+# code taken from the old raster-based nc writing in write.magpie
+# https://raw.githubusercontent.com/pik-piam/magclass/19fc7d098fbbea4af26240d6472e7e088356bf57/R/write.magpie.R
+writeRasterBrick <- function(x, filePath, comment = NULL, zname = "Time", ...) {
+  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)]])
+  }
+  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 (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"
+    }
+  }
+  raster::writeRaster(.sub(x, varnames[1]), filename = filePath, format = "CDF", overwrite = TRUE,
+                      compression = 9, zname = zname, zunit = zunit, varname = varnames[1], varunit = unit, ...)
+  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]) {
+      nc <- ncdf4::ncvar_add(nc, ncdf4::ncvar_def(i, unit, nc$dim, compression = 9))
+      ncdf4::ncvar_put(nc, i, aperm(as.array(.sub(x, i)), c(2, 1, 3)))
+    }
+  }
+  ncdf4::nc_close(nc)
+}
+
 withr::local_options(list(magclass_sizeLimit = 1e+12))
 
 
@@ -194,7 +238,7 @@ saveRDS(states,paste0(outputdir,"/states.rds"))
 gc()
 states <- convertLUH2(states)
 gc()
-write.magpie(states, paste0(out_dir, "/LUH2_states.nc"), comment = "unit: fraction of grid-cell area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(states, paste0(out_dir, "/LUH2_states.nc"), comment = "unit: fraction of grid-cell area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(states)
 gc()
 }
@@ -217,7 +261,7 @@ gc()
 if(!file.exists(paste0(out_dir,"/LUH2_protected_area.nc"))){
 b <- convertLUH2(b)
 gc()
-write.magpie(b, paste0(out_dir, "/LUH2_protected_area.nc"), comment = "unit: fraction of grid-cell", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(b, paste0(out_dir, "/LUH2_protected_area.nc"), comment = "unit: fraction of grid-cell", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(b)
 gc()
 }
@@ -245,7 +289,7 @@ if(!is.null(harvested_area_timber(gdx,level = "cell"))) {
   if(!file.exists(paste0(out_dir,"/LUH2_wood_harvest_area.nc"))){
     b <- convertLUH2(b)
     gc()
-    write.magpie(b, paste0(out_dir, "/LUH2_wood_harvest_area.nc"), comment = "unit: fraction of grid-cell area per year", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+    writeRasterBrick(b, paste0(out_dir, "/LUH2_wood_harvest_area.nc"), comment = "unit: fraction of grid-cell area per year", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
     rm(a,b)
     gc()
   }
@@ -271,7 +315,7 @@ if(!is.null(harvested_area_timber(gdx,level = "cell"))) {
   if(!file.exists(paste0(out_dir,"/LUH2_wood_harvest_yields.nc"))){
     b <- convertLUH2(b)
     gc()
-    write.magpie(b, paste0(out_dir, "/LUH2_wood_harvest_yields.nc"), comment = "unit: m3 per ha per year", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+    writeRasterBrick(b, paste0(out_dir, "/LUH2_wood_harvest_yields.nc"), comment = "unit: m3 per ha per year", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
     rm(a,b)
     gc()
   }
@@ -285,7 +329,7 @@ if(!is.null(harvested_area_timber(gdx,level = "cell"))) {
   if(!file.exists(paste0(out_dir,"/LUH2_wood_harvest_biomass_split.nc"))){
     b <- convertLUH2(b)
     gc()
-    write.magpie(b, paste0(out_dir, "/LUH2_wood_harvest_biomass_split.nc"), comment = "unit: fraction of wood harvest biomass", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+    writeRasterBrick(b, paste0(out_dir, "/LUH2_wood_harvest_biomass_split.nc"), comment = "unit: fraction of wood harvest biomass", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
     rm(b)
     gc()
   }
@@ -304,7 +348,7 @@ if (any(abs(d) > 0.1 )) message(paste0("Difference between cluster and grid cell
 if(!file.exists(paste0(out_dir,"/LUH2_irrigation.nc"))){
 irrig_hr_shr <- convertLUH2(irrig_hr_shr)
 gc()
-write.magpie(irrig_hr_shr, paste0(out_dir, "/LUH2_irrigation.nc"), comment = "unit: fraction of crop area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(irrig_hr_shr, paste0(out_dir, "/LUH2_irrigation.nc"), comment = "unit: fraction of crop area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(irrig_hr_shr,d)
 gc()
 }
@@ -322,7 +366,7 @@ if (any(abs(d) > 0.1)) message(paste0("Difference between cluster and grid cell
 if (!file.exists(paste0(out_dir, "/LUH2_flood.nc"))) {
   flooded <- convertLUH2(flooded)
   gc()
-  write.magpie(flooded, paste0(out_dir, "/LUH2_flood.nc"), comment = "unit: flooded fraction of C3 annual crop area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+  writeRasterBrick(flooded, paste0(out_dir, "/LUH2_flood.nc"), comment = "unit: flooded fraction of C3 annual crop area", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
   rm(flooded, d)
   gc()
 }
@@ -341,7 +385,7 @@ getNames(bio_hr_shr,dim=1) <- c("crpbf_c4per","crpbf_c3per")
 if(!file.exists(paste0(out_dir,"/LUH2_bioenergy.nc"))){
 bio_hr_shr <- convertLUH2(bio_hr_shr)
 gc()
-write.magpie(bio_hr_shr, paste0(out_dir, "/LUH2_bioenergy.nc"), comment = "unit: fraction of crop type area occupied by biofuel crops", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(bio_hr_shr, paste0(out_dir, "/LUH2_bioenergy.nc"), comment = "unit: fraction of crop type area occupied by biofuel crops", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(bio_hr_shr,d)
 gc()
 }
@@ -398,19 +442,19 @@ x <- NULL
 if(!file.exists(paste0(out_dir,"/LUH2_Nitrogen_fertilizer.nc"))){
 x <- convertLUH2(clean_magpie(collapseNames(a[,,"fertilizer"],collapsedim = 3.1)))
 gc()
-write.magpie(x, paste0(out_dir, "/LUH2_Nitrogen_fertilizer.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(x, paste0(out_dir, "/LUH2_Nitrogen_fertilizer.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 }
 
 if(!file.exists(paste0(out_dir,"/LUH2_Nitrogen_manure.nc"))){
 x <- convertLUH2(clean_magpie(collapseNames(a[,,"manure"],collapsedim = 3.1)))
 gc()
-write.magpie(x, paste0(out_dir, "/LUH2_Nitrogen_manure.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(x, paste0(out_dir, "/LUH2_Nitrogen_manure.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 }
 
 if(!file.exists(paste0(out_dir,"/LUH2_Nitrogen_surplus.nc"))){
 x <- convertLUH2(clean_magpie(collapseNames(a[,,"surplus"],collapsedim = 3.1)))
 gc()
-write.magpie(x, paste0(out_dir, "/LUH2_Nitrogen_surplus.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(x, paste0(out_dir, "/LUH2_Nitrogen_surplus.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 }
 
 rm(a,x,weight)
@@ -430,7 +474,7 @@ a[,,n]<-dimSums(yield_kr_su[,,n],dim=3)/dimSums(map_LUHMAg_grid[,,n],dim=3)
 if(!file.exists(paste0(out_dir,"/LUH2_Yield_DM.nc"))){
 a <- convertLUH2(a)
 gc()
-write.magpie(a, paste0(out_dir, "/LUH2_Yield_DM.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(a, paste0(out_dir, "/LUH2_Yield_DM.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(a,yield_kr,yield_kr_su)
 gc()
 }
@@ -449,7 +493,7 @@ a[,,n]<-dimSums(yield_kr_su[,,n],dim=3)/dimSums(map_LUHMAg_grid[,,n],dim=3)
 if(!file.exists(paste0(out_dir,"/LUH2_Yield_DM_rainfed.nc"))){
 a <- convertLUH2(a)
 gc()
-write.magpie(a, paste0(out_dir, "/LUH2_Yield_DM_rainfed.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(a, paste0(out_dir, "/LUH2_Yield_DM_rainfed.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(a,yield_kr,yield_kr_su)
 gc()
 }
@@ -468,7 +512,7 @@ a[,,n]<-dimSums(yield_kr_su[,,n],dim=3)/dimSums(map_LUHMAg_grid[,,n],dim=3)
 if(!file.exists(paste0(out_dir,"/LUH2_Yield_DM_irrigated.nc"))){
 a <- convertLUH2(a)
 gc()
-write.magpie(a, paste0(out_dir, "/LUH2_Yield_DM_irrigated.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(a, paste0(out_dir, "/LUH2_Yield_DM_irrigated.nc"), comment = "unit: tDM-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
 rm(a,yield_kr,yield_kr_su)
 gc()
 }
@@ -487,7 +531,7 @@ a[,,n]<-dimSums(yield_kr_su[,,n],dim=3)/dimSums(map_LUHMAg_grid[,,n],dim=3)
 if(!file.exists(paste0(out_dir,"/LUH2_Yield_Nr.nc"))){
 a <- convertLUH2(a)
 gc()
-write.magpie(a, paste0(out_dir, "/LUH2_Yield_Nr.nc"), comment = "unit: kgN-per-ha", datatype = "FLT8S", zname = "time", xname = "lon", yname = "lat")
+writeRasterBrick(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()
 }