diff --git a/DESCRIPTION b/DESCRIPTION
index f07d184eb2013b489db8a4a9195027486454c4ba..2fbd176abba139123b0471e68aa94e001da912c5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: CropWat
Type: Package
Title: R Implementation of the FAO CropWat Model
-Version: 0.1.0
+Version: 0.1.0.9000
Authors@R: c(
person("David", "Dorchies", , "david.dorchies@inrae.fr", role = c("aut", "cre"),
comment = c(ORCID = "0000-0002-6595-7984")),
@@ -21,7 +21,8 @@ RoxygenNote: 7.3.2
Roxygen: list(markdown = TRUE)
Suggests:
knitr,
- rmarkdown
+ rmarkdown,
+ testthat (>= 3.0.0)
VignetteBuilder: knitr
Imports:
dplyr,
@@ -29,7 +30,9 @@ Imports:
lubridate,
readr,
rlang,
+ stats,
tibble,
tidyr
Depends:
R (>= 2.10)
+Config/testthat/edition: 3
diff --git a/R/CW_create_input.R b/R/CW_create_input.R
index 382b03428a8b78bc74de0b084e182f246aea8d07..9e67d4ed671f0df9a7da85f638c7f0ae5043e0a3 100644
--- a/R/CW_create_input.R
+++ b/R/CW_create_input.R
@@ -1,5 +1,10 @@
#' Create CropWat model input
#'
+#' @details
+#' The crop cycle begins at the first matching sowing date in the simulation
+#' period. So, first time steps of simulation if are not in a crop cycle even
+#' if they occur during the crop cycle.
+#'
#' @inheritParams get_crop_params
#' @inheritParams calc_TAW
#' @param DatesR Simulation period ([vector] of [Date])
@@ -13,9 +18,9 @@
#' # Import example climate dataset
#' data(ZH_3_clim)
#' str(ZH_3_clim)
-#' # Selecting year 2010
+#' # Selecting years 2010 and 2011
#' meteo <-
-#' ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") & ZH_3_clim$Date <= as.Date("2010-12-31"), ]
+#' ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") & ZH_3_clim$Date <= as.Date("2011-12-31"), ]
#' # Create model input
#' cw_input <- CW_create_input("SB2023-soja",
#' DatesR = meteo$Date,
@@ -28,34 +33,48 @@
#' plot(cw_input)
#'
CW_create_input <- function(crop,
- year,
- DatesR = seq(as.Date(paste0(year, "-01-01")), as.Date(paste0(year, "-12-31")), by = "1 day"),
+ DatesR,
ETo,
P,
soil_depth,
AWC,
cp = get_crop_params(crop),
sowing_date = cp$sowing_date) {
- year <- lubridate::year(DatesR[1])
- Kc <- calc_Kc(cp, year = year, sowing_date = sowing_date)
- isCycle <- calc_isCycle(cp, year = year, sowing_date = sowing_date)
- Zr <- calc_root_depth(cp,
- soil_depth,
- year,
- sowing_date)
- TAW <- calc_TAW(cp, AWC, soil_depth, year, sowing_date)
- p <- calc_p(cp, year, sowing_date)
- cw_input <- tibble(
- DatesR = DatesR,
- isCycle = isCycle,
- P = P,
- Kc = Kc,
- Zr = Zr,
- ETc = Kc * ETo,
- p = p,
- TAW = TAW,
- RAW = calc_RAW(TAW, p)
- )
+
+ cycles <- get_yearly_cycles(DatesR)
+
+ l <- lapply(cycles, function(cycle) {
+ dates <- cycle$dates
+ Kc <- calc_Kc(cp, DatesR = dates, sowing_date = sowing_date)
+ isCycle <- calc_isCycle(cp, DatesR = dates, sowing_date = sowing_date)
+ Zr <- calc_root_depth(
+ cp,
+ soil_depth = soil_depth,
+ DatesR = dates,
+ sowing_date = sowing_date
+ )
+ TAW <- calc_TAW(
+ cp,
+ AWC = AWC,
+ soil_depth = soil_depth,
+ DatesR = dates,
+ sowing_date = sowing_date,
+ root_depth = Zr
+ )
+ p <- calc_p(cp, dates, sowing_date)
+ cw_input <- tibble(
+ DatesR = dates,
+ isCycle = isCycle,
+ P = P[cycle$which],
+ Kc = Kc,
+ Zr = Zr,
+ ETc = Kc * ETo[cycle$which],
+ p = p,
+ TAW = TAW,
+ RAW = calc_RAW(TAW, p)
+ )
+ })
+ cw_input <- bind_rows(l)
attr(cw_input, "crop_params") <- cp
attr(cw_input, "soil_params") <- list(soil_depth = soil_depth, AWC = AWC)
class(cw_input) <- c("CW_input", "CW_TS", class(cw_input))
diff --git a/R/CW_model.R b/R/CW_model.R
index d56564e1ae075eb86077ee8787e73fe9c8c19951..351a836886e442c00a5fc9bf6f114a5b43c79e55 100644
--- a/R/CW_model.R
+++ b/R/CW_model.R
@@ -27,13 +27,13 @@
CW_model <- function(X, input, FUN_IRRIG = CW_irrig_Dr) {
Dr <- X$Dr
- if (input$isCycle == 0 && X$isCycle == 1) {
- # TODO Check what does that mean
+ if (!input$isCycle && X$isCycle) {
+ # What happens during the harvest day?
Dr <- Dr * input$TAW / X$TAW - input$P # ???
}
Dr <- max(0, Dr - input$P)
- if (Dr > input$RAW) {
+ if (input$isCycle && Dr > input$RAW) {
Ks <- max(0, (input$TAW - Dr) / (input$TAW - input$RAW))
} else {
Ks <- 1
diff --git a/R/calc_Kc.R b/R/calc_Kc.R
index cc413ddf9a0abb592c21e2166703b8eb82e8a345..202e2b0cde9dad95803e85d319c8b61006736796 100644
--- a/R/calc_Kc.R
+++ b/R/calc_Kc.R
@@ -1,37 +1,60 @@
-#' Compute crop coefficient time series for one crop cycle or calendar year
+#' Compute crop parameters time series from crop parameters pivot points
+#'
+#' These functions compute the following parameters:
+#'
+#' - `Calc_Kc`: crop coefficient $K_c$
+#' - `calc_p`: critical depletion fraction $p=RAW/TAW$
+#' - `calc_root_depth`: root depth in m
+#' - `calc_RAW`: Readily Available Water (RAW) in mm
+#' - `calc_TAW`: Total Available Water (TAW) in mm
#'
#' @param cp Crop parameters (See [get_crop_params])
-#' @param year Year of simulation (Only used for detecting leap year).
+#' @param AWC Available Water Capacity (mm)
+#' @param soil_depth Soil depth (m)
+#' @param root_depth Root depth time series (m). See [calc_root_depth]
+#' @param TAW Total Available Water (mm). See [calc_TAW]
+#' @param p fraction of Readily Available Water time series. See [calc_p]
+#' @param DatesR A [vector] of continuous [Date]
#' If `NULL`, the calculation is limited to the crop cycle of the plant
#' @param sowing_date Sowing date in format "MM-DD"
#'
-#' @return A [vector] of Kc for each day of the year or the crop cycle.
+#' @return A [vector] of the parameter for each day of the crop cycle or the period defined by `DatesR`.
+#'
+#' @details
+#' For `calc_TAW`, parameters `cp`, `year`, and `sowing_date` are useless if `root_depth` is
+#' provided.
+#'
+#'
#' @export
+#' @rdname calc_params
#'
#' @examples
-#' Kc <- calc_Kc(get_crop_params("SB2023-soja"), year = NULL)
+#' # Compute Kc for the crop cycle
+#' Kc <- calc_Kc(get_crop_params("SB2023-soja"))
#' plot(Kc, type = "l")
#'
-#' Kc <- calc_Kc(get_crop_params("SB2023-soja"), year = 2024)
+#' # Compute Kc for a given period (less than one year)
+#' # with default sowing date defined in crop parameters
+#' DatesR <- seq(as.Date("2010-03-01"), as.Date("2010-10-31"), by = "1 day")
+#' Kc <- calc_Kc(get_crop_params("SB2023-soja"), DatesR = DatesR)
#' plot(Kc, type = "l")
#'
+#' # Compute Kc for a given period with user defined sowing date
+#' Kc <- calc_Kc(get_crop_params("FAO-MAIZE"),
+#' DatesR = DatesR,
+#' sowing_date = "05-01")
+#' plot(Kc, type = "l")
calc_Kc <- function(cp,
- year,
+ DatesR = NULL,
sowing_date = cp$sowing_date) {
-
# TODO: stopifnot cp, year, sowing_date format
- kc <- c(
- rep(cp$Kini, cp$Lini),
- cp$Kini + ((0:(cp$Ldev - 1)) / (cp$Ldev)) * (cp$Kmax -
- cp$Kini),
- #pb dernière valeur égale à Kc_mid, 1 j trop tôt
- rep(cp$Kmax, cp$Lmid),
- cp$Kmax + ((1:cp$Lend) / cp$Lend) * (cp$Kend - cp$Kmax)
+ calc_interpolated_param(
+ DatesR,
+ sowing_date,
+ x = c(cp$Lini, cp$Ldev, cp$Lmid, cp$Lend),
+ y = c(cp$Kini, cp$Kmax, cp$Kmax, cp$Kend),
+ yleft = cp$Kini,
+ yright = cp$Kini
)
-
- if (!is.null(year)) {
- kc <- complete_year(kc, year, sowing_date)
- }
- return(kc)
}
diff --git a/R/calc_RAW.R b/R/calc_RAW.R
index 180f2bb905d8c818446973f715e793ff39653568..64ce44ccc65788332c062bc3454307e3138c50f6 100644
--- a/R/calc_RAW.R
+++ b/R/calc_RAW.R
@@ -1,20 +1,5 @@
-#' Calculation of Readily Available Water (RAW)
-#'
-#' @param TAW Total Available Water (mm). See [calc_TAW]
-#' @param p fraction of Readily Available Water time series. See [calc_p]
-#'
-#' @return A [vector] of Readily Available Water.
#' @export
-#'
-#' @examples
-#' cp <- get_crop_params("SB2023-soja")
-#' TAW <- calc_TAW(cp,
-#' AWC = 140,
-#' soil_depth = 1.2,
-#' year = 2024)
-#' p <- calc_p(cp, year = 2024)
-#' RAW <- calc_RAW(TAW, p)
-#' plot(RAW, type = "l")
+#' @rdname calc_params
calc_RAW <- function(TAW, p) {
return(p * TAW)
}
diff --git a/R/calc_TAW.R b/R/calc_TAW.R
index 54e60612fc550ff9ec9f0c794f6eaaaabb82084f..c15c724fe3feac8c824f8fae3a406bd63c5a2126 100644
--- a/R/calc_TAW.R
+++ b/R/calc_TAW.R
@@ -1,29 +1,11 @@
-#' Calculation of Total Available Water (TAW)
-#'
-#' @details
-#' Parameters `cp`, `year`, and `sowing_date` are useless if `root_depth` is
-#' provided.
-#'
-#' @inheritParams calc_root_depth
-#' @param AWC Available Water Capacity (mm)
-#' @param root_depth Root depth time series (m). See [calc_root_depth].
-#'
-#' @return A [vector] of Total Available Water (mm) for each day of the year or the crop cycle.
#' @export
-#'
-#' @examples
-#' TAW <- calc_TAW(get_crop_params("SB2023-soja"),
-#' AWC = 140,
-#' soil_depth = 1.2,
-#' year = 2024)
-#' plot(TAW, type = "l")
-#'
+#' @rdname calc_params
calc_TAW <- function(cp,
AWC,
soil_depth,
- year,
+ DatesR = NULL,
sowing_date = cp$sowing_date,
- root_depth = calc_root_depth(cp, soil_depth, year, sowing_date)) {
+ root_depth = calc_root_depth(cp, soil_depth, DatesR, sowing_date)) {
return(AWC / soil_depth * root_depth)
}
diff --git a/R/calc_isCycle.R b/R/calc_isCycle.R
index 3fbd390e51653ad1a17583ed8b582d901f07b182..ec4482973b1f03009ecc12b8533d3794d295bfc2 100644
--- a/R/calc_isCycle.R
+++ b/R/calc_isCycle.R
@@ -6,12 +6,15 @@
#' @export
#'
#' @examples
-#' isCycle <- calc_isCycle(get_crop_params("SB2023-soja"), year = 2024)
+#' DatesR <- seq(as.Date("2010-03-01"), as.Date("2010-10-31"), by = "1 day")
+#' isCycle <- calc_isCycle(get_crop_params("SB2023-soja"), DatesR)
+#' plot(isCycle)
#'
-calc_isCycle <- function(cp, year, sowing_date = cp$sowing_date) {
- isCycle <- rep(TRUE, cp$Lini + cp$Ldev + cp$Lmid + cp$Lend)
- if (!is.null(year)) {
- isCycle <- complete_year(isCycle, year, sowing_date)
- }
+calc_isCycle <- function(cp, DatesR, sowing_date = cp$sowing_date) {
+ Ltot <- c(cp$Lini, cp$Ldev, cp$Lmid, cp$Lend)
+ l <- get_period(DatesR, sowing_date, Ltot)
+ isCycle <- c(rep(FALSE, l$doy_start),
+ rep(TRUE, sum(Ltot)),
+ rep(FALSE, max(0, length(DatesR) - sum(Ltot) - l$doy_start)))
return(isCycle)
}
diff --git a/R/calc_p.R b/R/calc_p.R
index 1a68c483dc4951f41e05278c881b917e67d763d3..775183ef9e1e303c9d6c3fd9325353b6aecee298 100644
--- a/R/calc_p.R
+++ b/R/calc_p.R
@@ -1,30 +1,14 @@
-#' Compute critical depletion fraction (p) time series for one crop cycle or calendar year
-#'
-#' *p* is the fraction RAW / TAW
-#'
-#' @inheritParams calc_Kc
-#'
-#' @return A [vector] of critical depletion fraction for each day of the year or the crop cycle.
#' @export
-#'
-#' @examples
-#' p <- calc_p(get_crop_params("SB2023-soja"),
-#' year = 2024)
-#' plot(p, type = "l")
-#'
+#' @rdname calc_params
calc_p <- function(cp,
- year,
+ DatesR = NULL,
sowing_date = cp$sowing_date) {
- p <- c(
- cp$p_ini + (0:(cp$Lini + cp$Ldev - 1
- )) * (cp$p_mid - cp$p_ini) / (cp$Lini + cp$Ldev - 1)
- ,
- rep(cp$p_mid, cp$Lmid)
- ,
- cp$p_mid + (1:(cp$Lend)) * (cp$p_end - cp$p_mid) / (cp$Lend)
+ calc_interpolated_param(
+ DatesR,
+ sowing_date,
+ x = c(1, cp$Lini + cp$Ldev - 1, cp$Lmid, cp$Lend),
+ y = c(cp$p_ini, cp$p_mid, cp$p_mid, cp$p_end),
+ yleft = cp$p_ini,
+ yright = cp$p_ini
)
- if (!is.null(year)) {
- p <- complete_year(p, year, sowing_date)
- }
- return(p)
}
diff --git a/R/calc_root_depth.R b/R/calc_root_depth.R
index 88965a9621ce9213b4edb2e66313ac6c67b7d3be..15d8a5a13b3777cd8e779fce2bad62cfb4e2c8f4 100644
--- a/R/calc_root_depth.R
+++ b/R/calc_root_depth.R
@@ -1,30 +1,17 @@
-#' Compute root depth time series for one crop cycle or calendar year
-#'
-#' @inheritParams calc_Kc
-#' @param soil_depth Soil depth (m)
-#'
-#' @return A [vector] of root depth (m) for each day of the year or the crop cycle.
#' @export
-#'
-#' @examples
-#' Zr <- calc_root_depth(get_crop_params("SB2023-soja"),
-#' soil_depth = 1.2,
-#' year = 2024)
-#' plot(Zr, type = "l")
-#'
+#' @rdname calc_params
calc_root_depth <- function(cp,
soil_depth,
- year,
+ DatesR = NULL,
sowing_date = cp$sowing_date) {
- Zr <- c(cp$Zini + (0:(cp$Lini + cp$Ldev - 1) / (cp$Lini + cp$Ldev)) * (cp$Zend - cp$Zini),
- rep(cp$Zend, cp$Lmid + cp$Lend))
-
- if (!is.null(year)) {
- Zr <- complete_year(Zr, year, sowing_date)
- }
-
- # Root depth is limited by soil depth
- Zr[Zr > soil_depth] <- soil_depth
-
- return(Zr)
+ Zini <- min(cp$Zini, soil_depth)
+ Zend <- min(cp$Zend, soil_depth)
+ calc_interpolated_param(
+ DatesR,
+ sowing_date,
+ x = c(1, cp$Lini + cp$Ldev - 1, cp$Lmid + cp$Lend),
+ y = c(Zini, Zend, Zend),
+ yleft = cp$Zini,
+ yright = cp$Zini
+ )
}
diff --git a/R/utils.R b/R/utils.R
index f6f612e72be6f5f3c5d506ac0b9df69c77ae681e..21b3a2aeeda305d122a43499be37b61bb093d96c 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -12,11 +12,42 @@ CW_data <- function(file) {
system.file("extdata", file, package = "CropWat")
}
-complete_year <- function(x, year, sowing_date) {
- # Cycle starting day (DOY)
- doy_start <- julian(as.Date(paste0(year, "-", sowing_date)), origin = as.Date(paste0(year, "-01-01"))) + 1
- nb_days_year <- ifelse(lubridate::leap_year(year), 366, 365)
- x_year <- rep(ifelse(is.logical(x), FALSE, x[1]), nb_days_year)
- x_year[doy_start:(doy_start + length(x) - 1)] <- x
- return(x_year)
+calc_interpolated_param <- function(DatesR, sowing_date, x, y, yleft, yright) {
+ l <- get_period(DatesR, sowing_date, x)
+ stats::approx(x = l$doy_start + cumsum(x), y = y, xout = l$xout, yleft = yleft, yright = yright)$y
+}
+
+get_period <- function(DatesR, sowing_date, x) {
+ if (!is.null(DatesR)) {
+ stopifnot(inherits(DatesR, 'Date'),
+ DatesR[length(DatesR)] <= DatesR[1] + lubridate::years(),
+ !is.na(sowing_date),
+ grepl("^[0-9]{1,2}-[0-9]{1,2}$", sowing_date))
+ year <- lubridate::year(DatesR[1])
+ doy_start <- julian(as.Date(paste0(year, "-", sowing_date)), origin = DatesR[1])
+ xout <- seq_along(DatesR)
+ } else {
+ doy_start <- 0
+ xout <- seq(sum(x))
+ }
+ return(list(doy_start = doy_start, xout = xout))
+}
+
+get_yearly_cycles <- function(DatesR) {
+ years <- unique(lubridate::year(DatesR))
+ start_dates <- rep(DatesR[1], length(years))
+ lubridate::year(start_dates) <- years
+ end_dates <- start_dates - lubridate::days() + lubridate::years()
+ end_dates[length(end_dates)] <- last(DatesR)
+ if (last(end_dates) <= last(start_dates)) {
+ years <- years[- length(years)]
+ }
+ lapply(setNames(seq_along(years), years), function(i) {
+ dates <- seq(start_dates[i], end_dates[i], by = "1 day")
+ list(year = years[i],
+ start = start_dates[i],
+ end = end_dates[i],
+ dates = dates,
+ which = which(DatesR %in% dates))
+ })
}
diff --git a/README.Rmd b/README.Rmd
index 4a1c5a5e4a82bebed246792d3dfe5dda4a699060..ddef6396e2823a25f7757efe029fe195f36a7646 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -17,6 +17,7 @@ knitr::opts_chunk$set(
# CropWat
<!-- badges: start -->
+
<!-- badges: end -->
CropWat is an R Implementation of the FAO CropWat Model.
@@ -43,13 +44,14 @@ irrigated crop.
```{r example}
# Load library
library(CropWat)
+library(dplyr)
# Import example climate dataset
data(ZH_3_clim)
# Selecting year 2010
meteo <- ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") &
- ZH_3_clim$Date <= as.Date("2010-12-31"), ]
+ ZH_3_clim$Date <= as.Date("2014-12-31"), ]
# Formatting model input
cw_input <- CW_create_input("SB2023-soja",
@@ -65,13 +67,13 @@ X <- CW_create_state(cw_input = cw_input)
# Choose an irrigation management
# Example: Fill Soil moisture depletion when half of RAW is reached
-fun_irrig_half_RAW <- CW_irrig_fun_factory(RAW_ratio = 0.5, apply_Dr = TRUE)
+fun_irrig_half_RAW <- CW_irrig_fun_factory(RAW_ratio = 1.25, apply_Dr = TRUE)
# Simulate water balance with an irrigation management
cw_output <- CW_run_simulation(X, cw_input, FUN_IRRIG = fun_irrig_half_RAW)
plot(cw_output)
# Total amount of irrigation applied (mm)
-sum(cw_output$Ir)
+cw_output |> group_by(lubridate::year(DatesR)) |> summarise(irrig_volume = sum(Ir))
```
diff --git a/README.md b/README.md
index 963042f80391391a069d7ab914eda05920231886..6c63728acec783a3b967b031be297d281b28f992 100644
--- a/README.md
+++ b/README.md
@@ -1,75 +1,91 @@
-
-<!-- README.md is generated from README.Rmd. Please edit that file -->
-
-# CropWat
-
-<!-- badges: start -->
-<!-- badges: end -->
-
-CropWat is an R Implementation of the FAO CropWat Model.
-
-This implements the functions describing water balance on an irrigated
-crop as described in FAO publications of the Irrigation and Drainage
-Series, namely, No. 56 “Crop Evapotranspiration - Guidelines for
-computing crop water requirements†and No. 33 titled “Yield response to
-waterâ€.
-
-## Installation
-
-You can install the development version of CropWat like so:
-
-``` r
-# install.packages("remotes")
-remotes::install_git("https://forgemia.inra.fr/umr-g-eau/cropwat.git")
-```
-
-## Example
-
-This is a basic example which shows you how to simulate water balance of
-an irrigated crop.
-
-``` r
-# Load library
-library(CropWat)
-
-# Import example climate dataset
-data(ZH_3_clim)
-
-# Selecting year 2010
-meteo <- ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") &
- ZH_3_clim$Date <= as.Date("2010-12-31"), ]
-
-# Formatting model input
-cw_input <- CW_create_input("SB2023-soja",
- DatesR = meteo$Date,
- ETo = meteo$ETP,
- P = meteo$Ptot,
- soil_depth = 1.2,
- AWC = 140)
-plot(cw_input)
-```
-
-<img src="man/figures/README-example-1.png" width="100%" />
-
-``` r
-
-# Initial state of the model set up
-X <- CW_create_state(cw_input = cw_input)
-
-# Choose an irrigation management
-# Example: Fill Soil moisture depletion when half of RAW is reached
-fun_irrig_half_RAW <- CW_irrig_fun_factory(RAW_ratio = 0.5, apply_Dr = TRUE)
-
-# Simulate water balance with an irrigation management
-cw_output <- CW_run_simulation(X, cw_input, FUN_IRRIG = fun_irrig_half_RAW)
-plot(cw_output)
-```
-
-<img src="man/figures/README-example-2.png" width="100%" />
-
-``` r
-
-# Total amount of irrigation applied (mm)
-sum(cw_output$Ir)
-#> [1] 375.3663
-```
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+# CropWat
+
+<!-- badges: start -->
+<!-- badges: end -->
+
+CropWat is an R Implementation of the FAO CropWat Model.
+
+This implements the functions describing water balance on an irrigated
+crop as described in FAO publications of the Irrigation and Drainage
+Series, namely, No. 56 “Crop Evapotranspiration - Guidelines for
+computing crop water requirements†and No. 33 titled “Yield response to
+waterâ€.
+
+## Installation
+
+You can install the development version of CropWat like so:
+
+``` r
+# install.packages("remotes")
+remotes::install_git("https://forgemia.inra.fr/umr-g-eau/cropwat.git")
+```
+
+## Example
+
+This is a basic example which shows you how to simulate water balance of
+an irrigated crop.
+
+``` r
+# Load library
+library(CropWat)
+library(dplyr)
+#>
+#> Attaching package: 'dplyr'
+#> The following objects are masked from 'package:stats':
+#>
+#> filter, lag
+#> The following objects are masked from 'package:base':
+#>
+#> intersect, setdiff, setequal, union
+
+# Import example climate dataset
+data(ZH_3_clim)
+
+# Selecting year 2010
+meteo <- ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") &
+ ZH_3_clim$Date <= as.Date("2014-12-31"), ]
+
+# Formatting model input
+cw_input <- CW_create_input("SB2023-soja",
+ DatesR = meteo$Date,
+ ETo = meteo$ETP,
+ P = meteo$Ptot,
+ soil_depth = 1.2,
+ AWC = 140)
+plot(cw_input)
+```
+
+<img src="man/figures/README-example-1.png" width="100%" />
+
+``` r
+
+# Initial state of the model set up
+X <- CW_create_state(cw_input = cw_input)
+
+# Choose an irrigation management
+# Example: Fill Soil moisture depletion when half of RAW is reached
+fun_irrig_half_RAW <- CW_irrig_fun_factory(RAW_ratio = 1.25, apply_Dr = TRUE)
+
+# Simulate water balance with an irrigation management
+cw_output <- CW_run_simulation(X, cw_input, FUN_IRRIG = fun_irrig_half_RAW)
+plot(cw_output)
+```
+
+<img src="man/figures/README-example-2.png" width="100%" />
+
+``` r
+
+# Total amount of irrigation applied (mm)
+cw_output |> group_by(lubridate::year(DatesR)) |> summarise(irrig_volume = sum(Ir))
+#> # A tibble: 5 × 2
+#> `lubridate::year(DatesR)` irrig_volume
+#> <dbl> <dbl>
+#> 1 2010 212.
+#> 2 2011 191.
+#> 3 2012 287.
+#> 4 2013 213.
+#> 5 2014 139.
+```
diff --git a/man/CW_create_input.Rd b/man/CW_create_input.Rd
index 9e8fc840a84c31f3823e1a3a2a579fefb7dfd689..7fc559e27846bef0324be705318a5561a489e699 100644
--- a/man/CW_create_input.Rd
+++ b/man/CW_create_input.Rd
@@ -6,9 +6,7 @@
\usage{
CW_create_input(
crop,
- year,
- DatesR = seq(as.Date(paste0(year, "-01-01")), as.Date(paste0(year, "-12-31")), by =
- "1 day"),
+ DatesR,
ETo,
P,
soil_depth,
@@ -20,9 +18,6 @@ CW_create_input(
\arguments{
\item{crop}{The code of the crop}
-\item{year}{Year of simulation (Only used for detecting leap year).
-If \code{NULL}, the calculation is limited to the crop cycle of the plant}
-
\item{DatesR}{Simulation period (\link{vector} of \link{Date})}
\item{ETo}{Potential Evaporation (mm/day)}
@@ -43,13 +38,18 @@ A \link{tibble} containing the input times series.
\description{
Create CropWat model input
}
+\details{
+The crop cycle begins at the first matching sowing date in the simulation
+period. So, first time steps of simulation if are not in a crop cycle even
+if they occur during the crop cycle.
+}
\examples{
# Import example climate dataset
data(ZH_3_clim)
str(ZH_3_clim)
-# Selecting year 2010
+# Selecting years 2010 and 2011
meteo <-
- ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") & ZH_3_clim$Date <= as.Date("2010-12-31"), ]
+ ZH_3_clim[ZH_3_clim$Date >= as.Date("2010-01-01") & ZH_3_clim$Date <= as.Date("2011-12-31"), ]
# Create model input
cw_input <- CW_create_input("SB2023-soja",
DatesR = meteo$Date,
diff --git a/man/calc_Kc.Rd b/man/calc_Kc.Rd
deleted file mode 100644
index 712604f9d78a9b09ee85505d583aff8aa93abbdb..0000000000000000000000000000000000000000
--- a/man/calc_Kc.Rd
+++ /dev/null
@@ -1,30 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/calc_Kc.R
-\name{calc_Kc}
-\alias{calc_Kc}
-\title{Compute crop coefficient time series for one crop cycle or calendar year}
-\usage{
-calc_Kc(cp, year, sowing_date = cp$sowing_date)
-}
-\arguments{
-\item{cp}{Crop parameters (See \link{get_crop_params})}
-
-\item{year}{Year of simulation (Only used for detecting leap year).
-If \code{NULL}, the calculation is limited to the crop cycle of the plant}
-
-\item{sowing_date}{Sowing date in format "MM-DD"}
-}
-\value{
-A \link{vector} of Kc for each day of the year or the crop cycle.
-}
-\description{
-Compute crop coefficient time series for one crop cycle or calendar year
-}
-\examples{
-Kc <- calc_Kc(get_crop_params("SB2023-soja"), year = NULL)
-plot(Kc, type = "l")
-
-Kc <- calc_Kc(get_crop_params("SB2023-soja"), year = 2024)
-plot(Kc, type = "l")
-
-}
diff --git a/man/calc_RAW.Rd b/man/calc_RAW.Rd
deleted file mode 100644
index 4759abf84499a8bf1b323f83789c9f99db0ba0f3..0000000000000000000000000000000000000000
--- a/man/calc_RAW.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/calc_RAW.R
-\name{calc_RAW}
-\alias{calc_RAW}
-\title{Calculation of Readily Available Water (RAW)}
-\usage{
-calc_RAW(TAW, p)
-}
-\arguments{
-\item{TAW}{Total Available Water (mm). See \link{calc_TAW}}
-
-\item{p}{fraction of Readily Available Water time series. See \link{calc_p}}
-}
-\value{
-A \link{vector} of Readily Available Water.
-}
-\description{
-Calculation of Readily Available Water (RAW)
-}
-\examples{
-cp <- get_crop_params("SB2023-soja")
-TAW <- calc_TAW(cp,
- AWC = 140,
- soil_depth = 1.2,
- year = 2024)
-p <- calc_p(cp, year = 2024)
-RAW <- calc_RAW(TAW, p)
-plot(RAW, type = "l")
-}
diff --git a/man/calc_TAW.Rd b/man/calc_TAW.Rd
deleted file mode 100644
index 2fef09401bd4865a7861e59bfaf2b4ec2e0d40e4..0000000000000000000000000000000000000000
--- a/man/calc_TAW.Rd
+++ /dev/null
@@ -1,47 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/calc_TAW.R
-\name{calc_TAW}
-\alias{calc_TAW}
-\title{Calculation of Total Available Water (TAW)}
-\usage{
-calc_TAW(
- cp,
- AWC,
- soil_depth,
- year,
- sowing_date = cp$sowing_date,
- root_depth = calc_root_depth(cp, soil_depth, year, sowing_date)
-)
-}
-\arguments{
-\item{cp}{Crop parameters (See \link{get_crop_params})}
-
-\item{AWC}{Available Water Capacity (mm)}
-
-\item{soil_depth}{Soil depth (m)}
-
-\item{year}{Year of simulation (Only used for detecting leap year).
-If \code{NULL}, the calculation is limited to the crop cycle of the plant}
-
-\item{sowing_date}{Sowing date in format "MM-DD"}
-
-\item{root_depth}{Root depth time series (m). See \link{calc_root_depth}.}
-}
-\value{
-A \link{vector} of Total Available Water (mm) for each day of the year or the crop cycle.
-}
-\description{
-Calculation of Total Available Water (TAW)
-}
-\details{
-Parameters \code{cp}, \code{year}, and \code{sowing_date} are useless if \code{root_depth} is
-provided.
-}
-\examples{
-TAW <- calc_TAW(get_crop_params("SB2023-soja"),
- AWC = 140,
- soil_depth = 1.2,
- year = 2024)
-plot(TAW, type = "l")
-
-}
diff --git a/man/calc_isCycle.Rd b/man/calc_isCycle.Rd
index 53f5642772d812061e618193a38bec38be9653d3..d7b0dafc66751af37fe417978708a3c9afcad290 100644
--- a/man/calc_isCycle.Rd
+++ b/man/calc_isCycle.Rd
@@ -4,12 +4,12 @@
\alias{calc_isCycle}
\title{Compute cycle period extend during the year}
\usage{
-calc_isCycle(cp, year, sowing_date = cp$sowing_date)
+calc_isCycle(cp, DatesR, sowing_date = cp$sowing_date)
}
\arguments{
\item{cp}{Crop parameters (See \link{get_crop_params})}
-\item{year}{Year of simulation (Only used for detecting leap year).
+\item{DatesR}{A \link{vector} of continuous \link{Date}
If \code{NULL}, the calculation is limited to the crop cycle of the plant}
\item{sowing_date}{Sowing date in format "MM-DD"}
@@ -21,6 +21,8 @@ A \link{vector} of \link{logical} of the crop cycle for each day of the year.
Compute cycle period extend during the year
}
\examples{
-isCycle <- calc_isCycle(get_crop_params("SB2023-soja"), year = 2024)
+DatesR <- seq(as.Date("2010-03-01"), as.Date("2010-10-31"), by = "1 day")
+isCycle <- calc_isCycle(get_crop_params("SB2023-soja"), DatesR)
+plot(isCycle)
}
diff --git a/man/calc_p.Rd b/man/calc_p.Rd
deleted file mode 100644
index 9438c02f00e807b3852af399b1e5e45e2280abd8..0000000000000000000000000000000000000000
--- a/man/calc_p.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/calc_p.R
-\name{calc_p}
-\alias{calc_p}
-\title{Compute critical depletion fraction (p) time series for one crop cycle or calendar year}
-\usage{
-calc_p(cp, year, sowing_date = cp$sowing_date)
-}
-\arguments{
-\item{cp}{Crop parameters (See \link{get_crop_params})}
-
-\item{year}{Year of simulation (Only used for detecting leap year).
-If \code{NULL}, the calculation is limited to the crop cycle of the plant}
-
-\item{sowing_date}{Sowing date in format "MM-DD"}
-}
-\value{
-A \link{vector} of critical depletion fraction for each day of the year or the crop cycle.
-}
-\description{
-\emph{p} is the fraction RAW / TAW
-}
-\examples{
-p <- calc_p(get_crop_params("SB2023-soja"),
- year = 2024)
-plot(p, type = "l")
-
-}
diff --git a/man/calc_params.Rd b/man/calc_params.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..ab605bd749c33d1742bed5a4fbd1f86f17e6f394
--- /dev/null
+++ b/man/calc_params.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calc_Kc.R, R/calc_RAW.R, R/calc_TAW.R,
+% R/calc_p.R, R/calc_root_depth.R
+\name{calc_Kc}
+\alias{calc_Kc}
+\alias{calc_RAW}
+\alias{calc_TAW}
+\alias{calc_p}
+\alias{calc_root_depth}
+\title{Compute crop parameters time series from crop parameters pivot points}
+\usage{
+calc_Kc(cp, DatesR = NULL, sowing_date = cp$sowing_date)
+
+calc_RAW(TAW, p)
+
+calc_TAW(
+ cp,
+ AWC,
+ soil_depth,
+ DatesR = NULL,
+ sowing_date = cp$sowing_date,
+ root_depth = calc_root_depth(cp, soil_depth, DatesR, sowing_date)
+)
+
+calc_p(cp, DatesR = NULL, sowing_date = cp$sowing_date)
+
+calc_root_depth(cp, soil_depth, DatesR = NULL, sowing_date = cp$sowing_date)
+}
+\arguments{
+\item{cp}{Crop parameters (See \link{get_crop_params})}
+
+\item{DatesR}{A \link{vector} of continuous \link{Date}
+If \code{NULL}, the calculation is limited to the crop cycle of the plant}
+
+\item{sowing_date}{Sowing date in format "MM-DD"}
+
+\item{TAW}{Total Available Water (mm). See \link{calc_TAW}}
+
+\item{p}{fraction of Readily Available Water time series. See \link{calc_p}}
+
+\item{AWC}{Available Water Capacity (mm)}
+
+\item{soil_depth}{Soil depth (m)}
+
+\item{root_depth}{Root depth time series (m). See \link{calc_root_depth}}
+}
+\value{
+A \link{vector} of the parameter for each day of the crop cycle or the period defined by \code{DatesR}.
+}
+\description{
+These functions compute the following parameters:
+}
+\details{
+\itemize{
+\item \code{Calc_Kc}: crop coefficient $K_c$
+\item \code{calc_p}: critical depletion fraction $p=RAW/TAW$
+\item \code{calc_root_depth}: root depth in m
+\item \code{calc_RAW}: Readily Available Water (RAW) in mm
+\item \code{calc_TAW}: Total Available Water (TAW) in mm
+}
+
+For \code{calc_TAW}, parameters \code{cp}, \code{year}, and \code{sowing_date} are useless if \code{root_depth} is
+provided.
+}
+\examples{
+# Compute Kc for the crop cycle
+Kc <- calc_Kc(get_crop_params("SB2023-soja"))
+plot(Kc, type = "l")
+
+# Compute Kc for a given period (less than one year)
+# with default sowing date defined in crop parameters
+DatesR <- seq(as.Date("2010-03-01"), as.Date("2010-10-31"), by = "1 day")
+Kc <- calc_Kc(get_crop_params("SB2023-soja"), DatesR = DatesR)
+plot(Kc, type = "l")
+
+# Compute Kc for a given period with user defined sowing date
+Kc <- calc_Kc(get_crop_params("FAO-MAIZE"),
+ DatesR = DatesR,
+ sowing_date = "05-01")
+plot(Kc, type = "l")
+}
diff --git a/man/calc_root_depth.Rd b/man/calc_root_depth.Rd
deleted file mode 100644
index 305afe7c6554529a2785cd2557d6e7a5f29128cf..0000000000000000000000000000000000000000
--- a/man/calc_root_depth.Rd
+++ /dev/null
@@ -1,31 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/calc_root_depth.R
-\name{calc_root_depth}
-\alias{calc_root_depth}
-\title{Compute root depth time series for one crop cycle or calendar year}
-\usage{
-calc_root_depth(cp, soil_depth, year, sowing_date = cp$sowing_date)
-}
-\arguments{
-\item{cp}{Crop parameters (See \link{get_crop_params})}
-
-\item{soil_depth}{Soil depth (m)}
-
-\item{year}{Year of simulation (Only used for detecting leap year).
-If \code{NULL}, the calculation is limited to the crop cycle of the plant}
-
-\item{sowing_date}{Sowing date in format "MM-DD"}
-}
-\value{
-A \link{vector} of root depth (m) for each day of the year or the crop cycle.
-}
-\description{
-Compute root depth time series for one crop cycle or calendar year
-}
-\examples{
-Zr <- calc_root_depth(get_crop_params("SB2023-soja"),
- soil_depth = 1.2,
- year = 2024)
-plot(Zr, type = "l")
-
-}
diff --git a/man/figures/README-example-1.png b/man/figures/README-example-1.png
index c6bb85000f3bbf6c9358008cd15ccca8ebc9f82e..d62834a14498b9cdd3eb6d90a38d1534e19ac27f 100644
Binary files a/man/figures/README-example-1.png and b/man/figures/README-example-1.png differ
diff --git a/man/figures/README-example-2.png b/man/figures/README-example-2.png
index 606f7e098a670e0967f2f68d1f06767df333a118..5b4f760f3d13b51e0c1ed62638d01263517166ba 100644
Binary files a/man/figures/README-example-2.png and b/man/figures/README-example-2.png differ
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000000000000000000000000000000000000..80d06924094400190a08075c0107dbce018b4c5e
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(CropWat)
+
+test_check("CropWat")
diff --git a/tests/testthat/test-calc_interpolated_param.R b/tests/testthat/test-calc_interpolated_param.R
new file mode 100644
index 0000000000000000000000000000000000000000..6cca315f600827ce25d3775853da0f2d249616ab
--- /dev/null
+++ b/tests/testthat/test-calc_interpolated_param.R
@@ -0,0 +1,44 @@
+test_that("calc_interpolated_param works", {
+ x <- c(2, 2, 1)
+ y <- c(1, 3, 2)
+ expect_equal(calc_interpolated_param(NULL, NA, x, y, 1, 1),
+ c(1, 1, 2, 3, 2))
+ DatesR <- seq(as.Date("2000-01-01"), as.Date("2000-01-07"), by = "1 day")
+ expect_error(calc_interpolated_param(DatesR, NA, x, y, 1, 1))
+ expect_error(calc_interpolated_param(DatesR, "1-223", x, y, 1, 1))
+ expect_error(calc_interpolated_param(DatesR, "04/04", x, y, 1, 1))
+ expect_error(calc_interpolated_param(seq(as.Date("2000-01-01"), as.Date("2001-01-02"), by = "1 day"), "04-04", x, y, 1, 1))
+ expect_equal(calc_interpolated_param(DatesR, "01-02", x, y, 1, 1),
+ c(1,1,1,2,3,2,1))
+})
+test_that("calc_interpolated_param works over new year's day :)", {
+ cp <- get_crop_params("SB2023-soja")
+ sowing_date = "11-01"
+ x <- c(cp$Lini, cp$Ldev, cp$Lmid, cp$Lend)
+ y <- c(cp$Kini, cp$Kmax, cp$Kmax, cp$Kend)
+ yleft <- cp$Kini
+ yright <- cp$Kini
+ DatesR <- seq(as.Date("2000-10-01"), as.Date("2001-09-30"), by = "1 day")
+ p <- calc_interpolated_param(
+ DatesR,
+ sowing_date = sowing_date,
+ x = x,
+ y = y,
+ yleft = yleft,
+ yright = yright
+ )
+ expect_equal(p[cp$Lini + 31], cp$Kini) # 10-31
+ expect_gt(p[cp$Lini + 32], cp$Kini) # 11-01 = sowing date
+ expect_lt(p[cp$Lini + cp$Ldev + 30], cp$Kmax) # Day before end of dev
+ expect_equal(p[cp$Lini + cp$Ldev + 31], cp$Kmax) # Max plateau
+ DatesR <- seq(as.Date("2001-01-01"), as.Date("2001-10-31"), by = "1 day")
+ p <- calc_interpolated_param(
+ DatesR,
+ sowing_date = sowing_date,
+ x = x,
+ y = y,
+ yleft = yleft,
+ yright = yright
+ )
+ expect_true(all(p == cp$Kini))
+})
diff --git a/tests/testthat/test-get_yearly_cycles.R b/tests/testthat/test-get_yearly_cycles.R
new file mode 100644
index 0000000000000000000000000000000000000000..2947a9f1337556b0cf6327e19e61ffc94d59e62b
--- /dev/null
+++ b/tests/testthat/test-get_yearly_cycles.R
@@ -0,0 +1,46 @@
+test_that("get_yearly_cycles works", {
+ # Regular calendar year
+ DatesR <- seq(as.Date("2019-01-01"), as.Date("2019-12-31"), by = "1 day")
+ expect_equal(get_yearly_cycles(DatesR), list("2019" =
+ list(
+ year = 2019,
+ start = DatesR[1],
+ end = last(DatesR),
+ dates = DatesR,
+ which = seq(365)
+ )))
+ # Calendar leap year
+ DatesR <- seq(as.Date("2020-01-01"), as.Date("2020-12-31"), by = "1 day")
+ expect_equal(get_yearly_cycles(DatesR), list("2020" =
+ list(
+ year = 2020,
+ start = DatesR[1],
+ end = last(DatesR),
+ dates = DatesR,
+ which = seq(366)
+ )))
+ # Not calendar year
+ DatesR <- seq(as.Date("2020-06-01"), as.Date("2021-05-31"), by = "1 day")
+ expect_equal(get_yearly_cycles(DatesR), list("2020" =
+ list(
+ year = 2020,
+ start = DatesR[1],
+ end = last(DatesR),
+ dates = DatesR,
+ which = seq(365)
+ )))
+ # Not complete year into calendar year
+ DatesR <- seq(as.Date("2020-06-01"), as.Date("2020-12-31"), by = "1 day")
+ expect_equal(get_yearly_cycles(DatesR), list("2020" =
+ list(
+ year = 2020,
+ start = DatesR[1],
+ end = last(DatesR),
+ dates = DatesR,
+ which = seq(length(DatesR))
+ )))
+ # 2 years
+ DatesR <- seq(as.Date("2019-01-01"), as.Date("2020-12-31"), by = "1 day")
+ cycles <- get_yearly_cycles(DatesR)
+
+})