Merge pull request magpiemodel#613 from k4rst3ns/f_67k_PR

Fix disaggregation scripts to be used with 67k input

CHANGELOG.md

@@ -27,7 +27,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **inputdata** There was another bug (terra default na.rm changed) in the inputdata that was fixed with rev4.93
 - **inputdata** There was a major bug (related to proj/terra) in the rev4.91 inputdata that was fixed with rev4.92
 - **scripts** Fixed a bug in NPI/NDC calculations leading to missing AD policies when run with 67k
-
+- **scripts** Fixed disaggregation.R and disaggregation_LUH2.R to be used with 67k
 
 ## [4.6.11] - 2023-09-05
 

scripts/output/extra/disaggregation.R

@@ -54,9 +54,14 @@ if (length(map_file) > 1) {
 #  Output functions
 # -----------------------------------------
 
+.fixCoords <- function(x) {
+  getSets(x, fulldim = FALSE)[1] <- "x.y.iso"
+  return(x)
+}
+
 .writeDisagg <- function(x, file, comment, message) {
-  write.magpie(x, file, comment = comment)
-  write.magpie(x, sub(".mz", ".nc", file), comment = comment, verbose = FALSE)
+  write.magpie(.fixCoords(x), file, comment = comment)
+  write.magpie(.fixCoords(x), sub(".mz", ".nc", file), comment = comment, verbose = FALSE)
 }
 
 .dissagcrop <- function(gdx, land_hr, map_file) {
@@ -82,7 +87,8 @@ if (length(map_file) > 1) {
   map <- readRDS(map_file)
 
   # create full time series
-  land_consv_hr <- new.magpie(map[, "cell"], getYears(land_consv_lr), getNames(wdpa_hr))
+  land_consv_hr <- new.magpie(map[, "cell"], getYears(land_consv_lr), getNames(wdpa_hr),
+                             sets = c("x.y.iso", "year", getSets(wdpa_hr, fulldim = FALSE)[3]))
   land_consv_hr[, getYears(land_consv_hr), ] <- wdpa_hr[, nyears(wdpa_hr), ]
   land_consv_hr[, getYears(wdpa_hr), ] <- wdpa_hr
 
@@ -91,7 +97,8 @@ if (length(map_file) > 1) {
       consv_prio_all <- read.magpie(consv_prio_hr_file)
       consv_prio_hr <- new.magpie(
         cells_and_regions = map[, "cell"],
-        names = getNames(consv_prio_all, dim = 2), fill = 0
+        names = getNames(consv_prio_all, dim = 2), fill = 0,
+        sets = c("x.y.iso", "year", "data")
       )
       iso <- readGDX(gdx, "iso")
       consv_iso <- readGDX(gdx, "policy_countries22")
@@ -299,6 +306,7 @@ land_hr <- interpolateAvlCroplandWeighted(
   snv_pol_shr = snv_pol_shr,
   snv_pol_fader = snv_pol_fader
 )
+land_hr <- .fixCoords(land_hr)
 
 # Write output
 .writeDisagg(land_hr, land_hr_out_file,
@@ -458,6 +466,7 @@ land_bii_hr <- interpolateAvlCroplandWeighted(
   snv_pol_fader = snv_pol_fader,
   unit = "share"
 )
+land_bii_hr <- .fixCoords(land_bii_hr)
 
 # Add primary and secondaray other land
 land_bii_hr <- PrimSecdOtherLand(land_bii_hr, land_hr_file)

scripts/output/extra/disaggregation_LUH2.R

@@ -46,15 +46,16 @@ out_dir<-paste0(outputdir,"/disaggregation_LUH2")
 convertLUH2 <- function(x) {
   #interpolate years
   years <- getYears(x,as.integer = TRUE)
+  getSets(x, fulldim = FALSE)[1] <- "x.y.iso"
   x <- toolFillYears(x,seq(range(years)[1],range(years)[2],by=1))
 
-
   for(n in seq(1995,2085,15)){
       x_1<- if(n==1995) as.RasterBrick(x[,n:(n+15),]) else  as.RasterBrick(x[,(n+1):(n+15),])
       x_aux<- if(n==1995) x_1 else stack(x_aux,x_1)
   }
   #re-project raster from 0.5 to 0.25 degree
-  x <- suppressWarnings(projectRaster(x_aux,raster(res=c(0.25,0.25)),method = "ngb"))
+  x <- suppressWarnings(raster::projectRaster(x_aux,raster::raster(res=c(0.25,0.25)),method = "ngb"))
+  crs(x) <- "+proj=utm +zone=1 +datum=WGS84"
   return(x)
 
 }
@@ -489,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()
-}
+}