Polynomial fit to cCO2

Author: kuadrat

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: growR
 Type: Package
-Version: 1.3.0.9007
+Version: 1.3.0.9008
 Date: 2024-05-23
 Authors@R: person(
     given = "Kevin",

NEWS.md

@@ -9,6 +9,8 @@
 * `growR_run_loop` allows to suppress returning the simulation objects. This 
   is useful to prevent overkill memory usage.
 * CO2 affects plant transpiration, according to `fCO2_transpiration_mod`.
+* Coefficients for polynomial fits to atmospheric CO2 concentration 
+  projections for different RCPs according to IPCC (`cCO2_coefficients`).
 
 ## Changed
 

R/data.R

@@ -126,3 +126,31 @@
 #' @md
 "parameter_scan_example"
 
+#' Coefficients for polynomial fits to atmospheric CO2 concentration scenarios.
+#'
+#' @description
+#' The fifth IPCC assessment report provides data for atmospheric CO2 
+#' concentration under different representative concentration pathway (RCP) 
+#' scenarios on decadal resolution. The values in this matrix represent the 
+#' coefficients of a second order polynomial fit to said data.
+#'
+#' @details
+#' These coefficients can be used to produce functions of the form `f(year) = 
+#' a3*year^3 + a2*year^2 + a1*year + a0` which return the estimated 
+#' atmospheric CO2 concentration in ppm for a given year.
+#' The validity of the polynomial fit is restricted to the years 2000 to 2100.
+#'
+#' @format A data.frame with columns
+#'   \describe{
+#'   \item{RCP}{Name of the RCP.}
+#'   \item{a3}{Coefficient of the quadratic term.}
+#'   \item{a2}{Coefficient of the quadratic term.}
+#'   \item{a1}{Coefficient of the linear term.}
+#'   \item{a0}{Coefficient of the constant term.}
+#'   }
+#'
+#' @seealso [atmospheric_CO2()], [WeatherData()]`$aCO2()`
+#'
+#' @md
+"cCO2_coefficients"
+

R/growR-package.R

@@ -8,7 +8,8 @@
 "_PACKAGE"
 
 globalVariables(c("yield_parameters",
-                  "management_parameters"))
+                  "management_parameters",
+                  "cCO2_coefficients"))
 
 #' Default options introduced by package growR
 #'

R/weather.R

@@ -286,14 +286,13 @@ WeatherData = R6Class(
     #' @description Choose the RCP to be used for the calculation of 
     #' atmospheric CO2 concentration.
     #'
-    #' @param RCP float or string. One of {2.6, 4.5, 8.5}.
+    #' @param RCP float or string. One of 2.6, 4.5, 8.5.
     #'
     set_RCP = function(RCP) {
       RCP = as.character(RCP)
-      allowed_RCPs = private$cCO2_table$RCP
-      if (!RCP %in% allowed_RCPs) {
+      if (!RCP %in% private$allowed_RCPs) {
         message = "[WeatherData]RCP `%s` not recognized. Allowed values:"
-        message = paste(message, paste(allowed_RCPs, collapse = ", "))
+        message = paste(message, paste(private$allowed_RCPs, collapse = ", "))
         warning(sprintf(message, RCP))
       } else {
         logger(sprintf("[WeatherData]Setting RCP to `%s`.", RCP), level = TRACE)
@@ -321,22 +320,21 @@ WeatherData = R6Class(
         if (is.null(RCP)) {
           RCP = private$RCP
         }
-        row = private$cCO2_table[private$cCO2_table$RCP == RCP, ]
-        c2020 = atmospheric_CO2(2020)
-        c2050 = row[["c2050"]]
-        c2100 = row[["c2100"]]
-        if (year < 2050) {
-          return((c2050 - c2020)/(2050 - 2020) * (year - 2020) + c2020)
-        } else {
-          return((c2100 - c2050)/(2100 - 2050) * (year - 2050) + c2050)
-        }
+        RCP_string = paste0("RCP", RCP)
+        row_selector = cCO2_coefficients[["order"]] == 3
+        row_selector = row_selector & cCO2_coefficients[["RCP"]] == RCP_string
+        row = cCO2_coefficients[row_selector, ]
+        res = row[["a3"]]*year^3 + 
+              row[["a2"]]*year^2 + 
+              row[["a1"]]*year + 
+              row[["a0"]]
+        return(res)
       }
     }
   ),
+
   private = list(
     RCP = "2.6",
-    cCO2_table = data.frame(RCP = c("2.6", "4.5", "8.5"),
-                            c2050 = c(450, 500, 550),
-                            c2100 = c(420, 600, 900))
+    allowed_RCPs = c("2.6", "4.5", "8.5")
   )
 )

data-raw/IPCC_CO2_concentration.R

@@ -0,0 +1,3 @@
+library(usethis)
+cCO2_coefficients = read.csv("inst/extdata/IPCC_CO2_concentration.csv")
+usethis::use_data(cCO2_coefficients)

data/cCO2_coefficients.rda

@@ -0,0 +1,9 @@
+order,RCP,a3,a2,a1,a0
+2,RCP2.6,0,-0.0198632946998996,81.8862194078263,-83949.6436918251
+2,RCP4.5,0,-0.0130613032884067,55.3716644757848,-58136.8986124938
+2,RCP6.0,0,0.0165884691095227,-64.9708171781507,63959.1605297345
+2,RCP8.5,0,0.0444312061725033,-176.484106167858,175614.278395903
+3,RCP2.6,0.000145298883194069,-0.913314482967777,1912.91709994921,-1334602.40391603
+3,RCP4.5,-0.00022427439053681,1.36601482116832,-2770.89462030362,1872293.53145368
+3,RCP6.0,-2.7338283163742E-05,0.18469314438927,-409.483036359746,299272.044888644
+3,RCP8.5,-2.09876022881761E-05,0.173485179548328,-440.966165860248,356263.99954439