Croparea update

.Rbuildignore

@@ -11,3 +11,4 @@
 ^codecov\.yml$
 ^\.claude$
 ^_pkgdown\.yml$
+^\.positai$

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '3741632'
+ValidationKey: '3928250'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'

.gitignore

@@ -3,3 +3,4 @@
 .RData
 .Ruserdata
 .Rprofile
+.positai

CITATION.cff

@@ -4,8 +4,8 @@ type: software
 title: |-
   mrlandcore: Fundamental MAgPIE preprocessing functions including
       calcLanduseInitialisation
-version: 1.8.1
-date-released: '2026-08-07'
+version: 1.9.0
+date-released: '2026-08-10'
 abstract: Fundamental MAgPIE preprocessing functions including calcLanduseInitialisation.
 authors:
 - family-names: Beier

DESCRIPTION

@@ -2,8 +2,8 @@ Type: Package
 Package: mrlandcore
 Title: Fundamental MAgPIE preprocessing functions including
     calcLanduseInitialisation
-Version: 1.8.1
-Date: 2026-08-07
+Version: 1.9.0
+Date: 2026-08-10
 Authors@R: c(
     person("Felicitas", "Beier", , "beier@pik-potsdam.de", role = c("aut", "cre")),
     person("Kristine", "Karstens", role = "aut"),
@@ -26,7 +26,7 @@ Description: Fundamental MAgPIE preprocessing functions including
 License: LGPL-3
 URL: https://github.com/pik-piam/mrlandcore
 Depends:
-    madrat (>= 2.20.9),
+    madrat (>= 3.37.4),
     magclass (>= 3.17),
     mrdownscale (>= 0.43.0),
     mrdrivers (>= 1.0.0),

R/calcCropareaFAOLUH.R

@@ -0,0 +1,224 @@
+#' @title calcCropareaFAOLUH
+#' @description Returns croparea from LUH (for cellular resolution) and
+#'              FAOSTAT (for country-level resolution).
+#'              Total harvested areas can be lower or higher than arable
+#'              land because of multicropping or fallow land.
+#'              Rice areas are distributed to flooded LUH areas. Additional FAOSTAT
+#'              rice areas are distributed based on country shares.
+#'
+#' @param sectoral   "area_harvested" returns croparea aggregated to FAO products,
+#'                   "ProductionItem" unaggregated ProdSTAT items,
+#'                   "FoodBalanceItem" Food Balance Sheet categories,
+#'                   "kcr" MAgPIE items, and "lpj" LPJmL items
+#' @param physical   if TRUE the sum over all crops agrees with the cropland area per country
+#' @param cellular   if TRUE: calculates cellular MAgPIE crop area for all magpie croptypes.
+#'                   Crop area from LUH3 crop types (c3ann, c4ann, c3per, c4per, cnfx)
+#'                   are mapped to MAgpIE crop types using "FAO2LUH2MAG_croptypes" and doing
+#'                   an intermediate step via harvested areas of FAO weight area within a
+#'                   specific LUH crop type to divide into MAgPIE crop types.
+#' @param irrigation If true: cellular areas are returned separated
+#'                   into irrigated and rainfed (see setup in calcLUH3)
+#'
+#' @return areas of individual crops from FAOSTAT and weight
+#'
+#' @author Ulrich Kreidenweis, Kristine Karstens, Felicitas Beier
+#'
+#' @importFrom utils read.csv
+#' @importFrom magclass setNames getCells collapseDim getItems
+#' @importFrom magpiesets findset addLocation
+#' @importFrom madrat toolAggregate toolGetMapping
+#' @importFrom mrfaocore toolExtrapolateFodder toolFAOcombine
+#' @importFrom utils person
+#' @importFrom withr local_options
+
+calcCropareaFAOLUH <- function(sectoral = "kcr", physical = TRUE, cellular = FALSE, irrigation = FALSE) {
+
+  local_options(magclass_sizeLimit = 1e+10)
+
+  if (!cellular) {
+
+    if (irrigation) stop("Irrigation levels for country based data not yet implemented!")
+
+    #################################
+    ### Croparea on country level ###
+    #################################
+
+    if (!is.null(sectoral) && !(sectoral == "lpj")) {
+
+      cropPrim <- readSource("FAO_online", "CropLive2010")[, , "area_harvested"]
+      # use linear_interpolate
+      fodder   <- readSource("FAO", "Fodder")[, , "area_harvested"]
+      fodder   <- toolExtrapolateFodder(fodder, endyear = max(getYears(cropPrim, as.integer = TRUE)))
+      data     <- toolFAOcombine(cropPrim, fodder) / 10^6 # convert to Mha
+
+      if (sectoral %in% c("FoodBalanceItem", "kcr")) {
+
+        aggregation <- toolGetMapping("FAOitems_online_2010update.csv", type = "sectoral",
+                                      where = "mrfaocore")
+        remove      <- setdiff(getNames(data, dim = 1), aggregation$post2010_ProductionItem)
+        data        <- data[, , remove, invert = TRUE]
+        data        <- toolAggregate(data, rel = aggregation, from = "post2010_ProductionItem",
+                                     to = ifelse(sectoral == "kcr", "k", sectoral),
+                                     dim = 3.1, partrel = TRUE)
+
+        if (sectoral == "kcr") {
+          # add bioenergy with 0 values
+          data <- add_columns(x = data, addnm = c("betr", "begr"), dim = 3.1)
+          data[, , c("betr", "begr")] <- 0
+
+          # remove all non kcr items
+          kcr    <- findset("kcr")
+          remove <- setdiff(getItems(data, dim = 3.1), kcr)
+
+          if (length(remove) > 0) {
+            remainArea <- mean(dimSums(data[, , "area_harvested"][, , remove], dim = 1) /
+                                 dimSums(dimSums(data[, , "area_harvested"], dim = 3), dim = 1))
+            if (remainArea > 0.02) vcat(1, "Aggregation created a 'remaining' category. The area harvested is",
+                                        round(remainArea, digits = 3) * 100, "% of total \n")
+            vcat(2, paste0("Data for the following items removed: ", remove))
+            data <- data[, , kcr]
+          }
+        }
+
+      } else if (sectoral != "ProductionItem") {
+        stop("Sectoral aggregation not supported")
+      }
+
+    } else if (sectoral == "lpj") {
+
+      magCroparea <- calcOutput("CropareaFAOLUH", sectoral = "kcr", physical = physical,
+                                cellular = FALSE, irrigation = FALSE, aggregate = FALSE)
+
+      mag2lpj     <- toolGetMapping(type = "sectoral", name = "MAgPIE_LPJmL.csv",
+                                    where = "mappingfolder")
+      mag2lpj     <- mag2lpj[!(mag2lpj$MAgPIE == "pasture"), ]
+
+      lpjCroparea <- toolAggregate(magCroparea, rel = mag2lpj, from = "MAgPIE", to = "LPJmL", dim = 3.1)
+      data        <- lpjCroparea
+
+    } else {
+      stop("Sectoral aggregation not supported")
+    }
+
+    # use the share of the single crops to calculate their "physical" area
+    if (physical) {
+      # 6620  = (6620|Arable land and Permanent crops or  6620|Cropland)
+      cropland        <- setNames(collapseNames(calcOutput("FAOLand",
+                                                           aggregate = FALSE)[, , "6620", pmatch = TRUE]), "crop")
+      harvestedShare <- data / dimSums(data, dim = 3.1)
+      commonyears     <- intersect(getYears(cropland), getYears(harvestedShare))
+      data            <- collapseNames(cropland[, commonyears, ] * harvestedShare[, commonyears, ])
+    }
+
+    data[is.na(data)] <- 0
+
+  } else {
+    ##################################
+    ### Croparea on cellular level ###
+    ##################################
+
+    if (sectoral == "kcr") {
+      # LUH related data input on cell level
+      luhWeights   <- calcOutput("LUH2MAgPIE", share = "MAGofLUH",
+                                 missing = "fill", rice = "non_flooded", aggregate = FALSE)
+
+      luhCroptypes <- c("c3ann", "c4ann", "c3per", "c4per", "c3nfx")
+
+      luhCroparea  <- calcOutput("LUH3", landuseTypes = "LUH3", aggregate = FALSE,
+                                 irrigation = irrigation, cellular = TRUE)[, , luhCroptypes]
+
+      commonYears  <- intersect(getYears(luhWeights), getYears(luhCroparea))
+      luhWeights   <- luhWeights[, commonYears, ]
+      luhCroparea  <- luhCroparea[, commonYears, ]
+
+      commonCountries <- intersect(getItems(luhWeights, dim = "ISO"), getItems(luhCroparea, dim = "iso"))
+
+      # corrected rice area (in Mha)
+      ricearea <- calcOutput("Ricearea", cellular = TRUE,
+                             share = FALSE, aggregate = FALSE)
+
+      commonYears  <- intersect(getYears(ricearea), getYears(luhCroparea))
+      luhWeights   <- luhWeights[, commonYears, ]
+      luhCroparea  <- luhCroparea[, commonYears, ]
+      ricearea     <- ricearea[, commonYears, ]
+
+      # irrigation
+      if (irrigation == TRUE) {
+        # for check
+        luhCropareaTotal <- dimSums(luhCroparea, dim = 3)
+
+        # calculate irrigation share for rice area correction
+        irrigShr <- new.magpie(cells_and_regions = getCells(luhCroparea),
+                               years = getYears(luhCroparea),
+                               names = getNames(luhCroparea), fill = NA)
+        luhCroparea <- add_columns(luhCroparea, addnm = "total", dim = 3.2, fill = NA)
+        luhCroparea[, , "total"] <-  dimSums(luhCroparea[, , c("rainfed", "irrigated")], dim = 3.2)
+
+        irrigShr[, , "irrigated"] <- collapseNames(ifelse(luhCroparea[, , "total"] > 0,
+                                                          luhCroparea[, , "irrigated"] / luhCroparea[, , "total"], 0))
+        irrigShr[, , "rainfed"]   <- 1 - collapseNames(irrigShr[, , "irrigated"])
+
+        # flooded rice areas
+        floodedRice <- collapseNames(ricearea[, , "flooded"] * irrigShr[, , "c3ann"])
+
+        luhCroparea <- luhCroparea[, , "total", invert = TRUE]
+
+      } else {
+        # for check
+        luhCropareaTotal <- dimSums(luhCroparea, dim = 3)
+
+        # flooded rice areas
+        floodedRice <- collapseNames(ricearea[, , "flooded"])
+
+      }
+
+      # temporarily exclude flooded rice for distribution of other crops and aerobic rice areas
+      luhCroparea[, , "c3ann"] <- luhCroparea[, , "c3ann"] - floodedRice
+
+      # correction of LUH cropareas with FAO country shares
+      luh2mag          <- luhCroparea * luhWeights[commonCountries, , ]
+      magCroparea      <- dimSums(luh2mag, dim = 3.1)
+
+      # total rice area correction
+      magCroparea[, , "rice_pro"] <- magCroparea[, , "rice_pro"] + floodedRice
+
+      # check sums
+      if (any(round(abs(dimSums(magCroparea, dim = 3) - luhCropareaTotal), digits = 6) > 1e-6)) {
+        stop("Sums after rice correction in calcCroparea don't match!")
+      }
+
+      data <- collapseNames(magCroparea)
+
+    } else if (sectoral == "lpj") {
+
+      magCroparea   <- calcOutput("CropareaFAOLUH", sectoral = "kcr", physical = physical,
+                                  cellular = TRUE, irrigation = irrigation, aggregate = FALSE)
+      mag2lpj       <- toolGetMapping(type = "sectoral", name = "MAgPIE_LPJmL.csv",
+                                      where = "mappingfolder")
+      mag2lpj       <- mag2lpj[!(mag2lpj$MAgPIE == "pasture"), ]
+      lpjCroparea   <- toolAggregate(magCroparea, rel = mag2lpj, from = "MAgPIE", to = "LPJmL", dim = "MAG")
+      data          <- lpjCroparea
+
+    } else {
+      stop("Not possible (for now) for the given item set (sectoral)!")
+    }
+
+    if (!physical) {
+
+      multiCropping   <- calcOutput("MulticroppingOld", aggregate = FALSE)
+      commonCountries <- intersect(getItems(multiCropping, dim = "ISO"), getItems(data, dim = "iso"))
+      data            <- data * multiCropping[commonCountries, getYears(data), ]
+    }
+  }
+
+  data <- collapseNames(data)
+
+  # not more precision than 1 ha needed. very small areas can make problems in some weighting scripts
+  data <- round(data, 6)
+
+  return(list(x            = data,
+              weight       = NULL,
+              unit         = "million ha",
+              description  = "croparea",
+              isocountries = !cellular))
+}

R/calcCropareaLandInG.R

@@ -10,7 +10,6 @@
 #'                    agrees with the physical cropland of readLandInG(subtype = physical)
 #' @param cellular    if TRUE: calculates cellular crop area for all magpie croptypes.
 #'                    Option FALSE is not (yet) available.
-#' @param cells       Switch between "magpiecell" (59199) and "lpjcell" (67420)
 #' @param irrigation  If true: cellular areas are returned separated
 #'                    into irrigated and rainfed
 #' @param selectyears extract certain years from the data
@@ -27,7 +26,7 @@
 #' @importFrom mstools toolHoldConstant
 #'
 calcCropareaLandInG <- function(sectoral = "kcr", physical = TRUE, cellular = FALSE,
-                                cells = "magpiecell", irrigation = FALSE, selectyears = "all",
+                                irrigation = FALSE, selectyears = "all",
                                 lpjml = c(natveg = "LPJmL4_for_MAgPIE_44ac93de",
                                           crop = "ggcmi_phase3_nchecks_bft_e511ac58"),
                                 climatetype = "GSWP3-W5E5:historical") {
@@ -50,7 +49,7 @@ calcCropareaLandInG <- function(sectoral = "kcr", physical = TRUE, cellular = FA
 
   ### Calculations ###
   # read in fallow land (for check below)
-  fallow <- calcOutput("FallowLand", aggregate = FALSE)
+  fallow <- calcOutput("FallowLandInG", aggregate = FALSE, cellular = TRUE)
 
   # year selection
   if (any(selectyears == "all")) {
@@ -280,15 +279,6 @@ calcCropareaLandInG <- function(sectoral = "kcr", physical = TRUE, cellular = FA
       # fill missing countries with 0
       cropArea <- toolConditionalReplace(x = toolCountryFill(cropArea),
                                          conditions = "is.na()", replaceby = 0)
-    } else {
-      if (cells == "magpiecell") {
-        cropArea <- toolCoord2Isocell(cropArea)
-      } else if (cells == "lpjcell") {
-        # this is already the format of cropArea
-      } else {
-        stop("This value for the cell parameter is not supported,
-          choose between \"magpiecell\" and \"lpjcell\"")
-      }
     }
     cropAreaList[[y]] <- cropArea
   }
@@ -300,5 +290,5 @@ calcCropareaLandInG <- function(sectoral = "kcr", physical = TRUE, cellular = FA
               weight = NULL,
               description = "Croparea for different croptypes",
               unit = "Mha",
-              isocountries = FALSE))
+              isocountries = !cellular))
 }

R/calcFallow.R

@@ -0,0 +1,31 @@
+#' @title calcFallow
+#' @description Returns fallow land calculated based on physical cropland extent
+#'              and harvested area output of default croparea datasource.
+#'              The resultant areas are scaled to match the land croparea land
+#'              class of calcLanduseInitialisation.
+#' @param cellular TRUE for cellular outputs
+#' @return MAgPIE object containing fallow land in Mha
+#' @author David Hoetten, Felicitas Beier, Benjamin Bodirsky
+#' @seealso
+#' \code{\link{readLandInG}}
+#' @examples
+#' \dontrun{
+#' calcOutput("Fallow")
+#' }
+#' @importFrom magclass dimSums mbind
+#' @importFrom madrat toolConditionalReplace
+#'
+calcFallow <- function(cellular = FALSE) {
+
+  fallow <- calcOutput("Croparea", cellular = cellular,
+                       physical = TRUE, fallow = TRUE,
+                       irrigation = FALSE, aggregate = FALSE)
+  fallow <- fallow[, , "fallow"]
+
+  return(list(x = fallow,
+              weight = NULL,
+              description = "Fallow land",
+              unit = "Mha",
+              isocountries = !cellular))
+
+}

R/calcFallowLand.R → R/calcFallowLandInG.R

@@ -1,4 +1,4 @@
-#' @title calcFallowLand
+#' @title calcFallowLandInG
 #' @description
 #' Calculates fallow land on grid cell level,
 #' based on physical cropland extend and harvested area output
@@ -16,12 +16,12 @@
 #' \code{\link{readLandInG}}
 #' @examples
 #' \dontrun{
-#' calcOutput("FallowLand")
+#' calcOutput("FallowLandInG")
 #' }
 #' @importFrom magclass dimSums mbind
 #' @importFrom madrat toolConditionalReplace
 #'
-calcFallowLand <- function(cellular = TRUE) {
+calcFallowLandInG <- function(cellular = TRUE) {
 
   harvestedArea <- readSource("LandInG", subtype = "harvestedArea")
 
@@ -47,6 +47,6 @@ calcFallowLand <- function(cellular = TRUE) {
               weight = NULL,
               description = "Fallow land",
               unit = "Mha",
-              isocountries = FALSE))
+              isocountries = !cellular))
 
 }