diff --git a/DESCRIPTION b/DESCRIPTION index 0674cfc..2d2a2c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: growR Type: Package -Version: 1.3.0.9000 +Version: 1.3.0.9001 Date: 2024-05-23 Authors@R: person( given = "Kevin", diff --git a/NEWS.md b/NEWS.md index 9719e37..2d4665c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,11 @@ ## Added * site_name and run_name are printed to console in `growR_run_loop`. +* PhenologicalAutocut + +## Changed + +* autocut algorithm is now outsourced to its own R6class. # growR 1.3.0 diff --git a/R/autocut.R b/R/autocut.R new file mode 100644 index 0000000..4efb083 --- /dev/null +++ b/R/autocut.R @@ -0,0 +1,436 @@ +#' Autocut +#' +#' @description +#' An algorithm to determine grass cut dates if none are provided. +#' This is an abstract class and not intended for direct use. Instead, use +#' its subclasses that implement a `determine_cut()` method. +#' +#' The expected number of cuts is estimated from management intensity and +#' site altitude based on data for Swiss grasslands by Huguenin et al. +#' +#' @seealso +#' [management_parameters] +#' +#' @references +#' \insertRef{huguenin2017GrundlagenDuengung}{growR} +#' +#' @md +#' @export +Autocut = R6Class( + "Autocut", + public = list( +#' @field MVS The [ModvegeSite] object for which to take the cut decision. + MVS = NULL, +#' @field n_cuts Number of cuts expected for this altitude and management +#' intensity. + n_cuts = NULL, + + #' @description Constructor + #' + #' @param MVS The [ModvegeSite] object for which cuts shall be + #' determined. + #' + initialize = function(MVS) { + self$MVS = MVS + self$n_cuts = self$get_expected_n_cuts(MVS$parameters$ELV, + MVS$management$intensity) + }, + + #' Get number of expected cuts + #' + #' Return the number of expected cuts for a site at a given *elevation* and + #' management *intensity*. + #' + #' This uses data.frame `management_parameters` as a lookup table and + #' interpolates linearly in between the specified values. + #' + #' @param elevation The elevation of the considered site in meters above sea + #' level. + #' @param intensity One of ("high", "middle", "low", "extensive"). Management + #' intensity for considered site. + #' + #' @return Number of expected cuts per season. + #' + get_expected_n_cuts = function(elevation, intensity = "high") { + mask = management_parameters$intensity == intensity + # Find the altitudes below and above given elevation + altitudes = management_parameters[mask, ]$altitude + n_cuts = management_parameters[mask, ]$n_cuts + i0 = which.min(abs(altitudes - elevation)) + alt0 = altitudes[i0] + n0 = n_cuts[i0] + # Check for second closest altitude + i1 = which.min(abs(altitudes[-i0] - elevation)) + alt1 = altitudes[-i0][i1] + n1 = n_cuts[-i0][i1] + # Approximate linearly + n = (n1 - n0) / (alt1 - alt0) * (elevation - alt0) + n0 + return(n) + }, + + #' @description + #' Empty method stub intended for overriding by inheriting subclasses. + #' + #' @param DOY Integer day-of-the-year. + determine_cut = function(DOY) { + stop(paste("`Autocut` is an abstract base class that does not implement the", + "`determine_cut` method. Use a subclass of `Autocut` instead.")) + } +## End of public attributes + ) +) + +#' Autocut based on phenology +#' +#' @description +#' An algorithm to determine grass cut dates if none are provided. +#' This uses empirical data for Switzerland to determine the first and last +#' cut dates of the season from meteorological data. +#' The number of cuts is inferred from Huguenen et al. and these cut events +#' are distributed equally between first and last cut dates. +#' +#' @seealso [management_parameters], [PetersenAutocut] +#' +#' @references +#' \insertRef{huguenin2017GrundlagenDuengung}{growR} +#' +#' @md +#' @export +PhenologicalAutocut = R6Class( + "PhenologicalAutocut", + inherit = Autocut, + public = list( +#' @field cut_DOYs vector containing the integer day-of-year's on which cuts +#' occur. + cut_DOYs = NULL, + + #' @description Constructor + #' + #' Valid cut dates are calculated upon initialization. + #' + #' @param MVS The [ModvegeSite] object for which cuts shall be + #' determined. + #' + #' @seealso [Autocut] + initialize = function(MVS) { + super$initialize(MVS) + + # Cut 200 degree-days before ST2 + target_temperature_sum = self$MVS$parameters$ST2 - 200 + first_cut_DOY = which.min(abs(self$MVS$ST - target_temperature_sum)) + + # Get last day of growth + i0 = 250 + i1 = 350 + threshold = 1 + deltas = diff(self$MVS$ST)[i0:i1] + lower = which(deltas < threshold) + if (length(lower) == 0) { + last_cut_DOY = i1 + } else { + last_cut_DOY = i0 + lower[[1]] + } + + # Time between cuts + integer_n_cuts = round(self$n_cuts) + delta = floor((last_cut_DOY - first_cut_DOY) / (integer_n_cuts - 1)) + self$cut_DOYs = (((1:integer_n_cuts) - 1) * delta) + first_cut_DOY + }, + + #' @description Does a cut occur on *DOY*? + #' + #' Check if *DOY* is in `self$cut_DOYs`. If so, return `TRUE`. Return + #' `FALSE` otherwise. + #' @param DOY Integer day of the year (1-366). + #' @return Boolean + determine_cut = function(DOY) { + if (DOY %in% self$cut_DOYs) { + return(TRUE) + } else { + return(FALSE) + } + } + ) +) + +#' Petersen autocut algorithm +#' +#' @description +#' Simulation routine to realistically predict grass cutting events. +#' This follows an implementation described in Petersen et al. (2021). +#' +#' @details +#' The decision to cut is made based on two criteria. +#' First, it is checked whether a *target biomass* is reached on given +#' DOY. The defined target depends on the DOY and is given through +#' :func:`get_target_biomass`. If said biomass is present, return `TRUE`. +#' +#' Otherwise, it is checked whether a given amount of time has passed +#' since the last cut. Depending on whether this is the first cut of +#' the season or not, the relevant parameters are +#' :int:`last_DOY_for_initial_cut` and :int:`max_cut_period`. +#' If that amount of time has passed, return `TRUE`, otherwise return +#' `FALSE`. +#' +#' The target biomass for a given day is determined following the principles +#' described in Petersen et al. +#' +#' The exact regression for the target biomass is based on Fig. S2 in the +#' supplementary material of Petersen et al. +#' +#' A refinement to expected yield as function of altitude has been +#' implemented according to Table 1a in Huguenin et al. (2017). +#' +#' @references +#' \insertRef{huguenin2017GrundlagenDuengung}{growR} +#' \insertRef{petersen2021DynamicSimulationManagement}{growR} +#' +#' @seealso [PhenologicalAutocut] +#' @md +#' @export +PetersenAutocut = R6Class( + "PetersenAutocut", + inherit = Autocut, + public = list( +#' @field last_DOY_for_initial_cut Start cutting after this DOY, +#' even if yield target is not reached. + last_DOY_for_initial_cut = 150, +#' @field max_cut_period Maximum period to wait between +#' subsequent cuts. + max_cut_period = 55, +#' @field dry_precipitation_limit Maximum amount of allowed +#' precipitation (mm) to consider a day. + dry_precipitation_limit = 1, +#' @field dry_days_before_cut Number of days that shold be dry +#' before a cut is made. + dry_days_before_cut = 1, +#' @field dry_days_after_cut Number of days that shold be dry +#' after a cut is made. + dry_days_after_cut = 2, +#' @field max_cut_delay Number of days a farmer is willing to +#' wait for dry conditions before a cut is made anyways. + max_cut_delay = 5, +#' @field cut_delays Vector to keep track of cut delay times. + cut_delays = c(0), +#' @field dry_window Logical that indicates if DOY at index is +#' considered dry enough to cut. + dry_window = NULL, +#' @field target_biomass Biomass amount that should to be reached +#' by given DOY for a cut to be made. + target_biomass = NULL, +#' @field end_of_cutting_season Determined DOY after which no +#' more cuts are made. + end_of_cutting_season = NULL, + + #' @description Constructor + #' + #' @param MVS The [ModvegeSite] object for which cuts shall be + #' determined. + initialize = function(MVS) { + super$initialize(MVS) + # Precalculate target biomass + self$target_biomass = sapply(1:MVS$days_per_year, + function(day) { + self$get_target_biomass(day, intensity = MVS$management$intensity) + } + ) + # Calculate end of the season and max cut interval. + self$end_of_cutting_season = + get_end_of_cutting_season(MVS$BM_after_cut, MVS$parameters$ELV, + MVS$management$intensity) + delta = self$end_of_cutting_season - self$last_DOY_for_initial_cut + self$max_cut_period = round(delta / self$n_cuts) + # Prepare the intervals of dry-enough cut conditions. + dry_days = MVS$weather$PP < self$dry_precipitation_limit + self$dry_window = logical(MVS$days_per_year) + for (i in 1:MVS$days_per_year) { + i0 = max(i - self$dry_days_before_cut, 1) + i1 = min(i + self$dry_days_after_cut, MVS$days_per_year) + all_dry = all(dry_days[i0:i1]) + self$dry_window[i] = all_dry + } + }, + + #' @description + #' Lookup table returning expected annual gross yields as function of + #' elevation and management intensity. + #' + #' Based on data from Table 1a in + #' Lookup Table of expected yield as functions of height and management + #' intensity after Olivier Huguenin et al. (2017). + #' + #' @param elevation The elevation of the considered site in meters above sea + #' level. + #' @param intensity One of ("high", "middle", "low", "extensive"). Management + #' intensity for considered site. + #' + #' @return Annual gross yield in t / ha (metric tons per hectare). Note that + #' 1 t/ha = 0.1 kg/m^2. + #' + #' @references + #' \insertRef{huguenin2017GrundlagenDuengung}{growR} + get_annual_gross_yield = function(elevation, intensity = "high") { + mask = yield_parameters$intensity == intensity + a = yield_parameters[mask, ]$a + b = yield_parameters[mask, ]$b + return(a + b * max(elevation, 500)) + }, + + #' @description + #' Get target value of biomass on given *DOY*, which determines whether + #' a cut is to occur. + #' + #' The regression for the target biomass is based on Fig. S2 in the + #' supplementary material of + #' Petersen, Krischan, David Kraus, Pierluigi Calanca, Mikhail A. + #' Semenov, Klaus Butterbach-Bahl, and Ralf Kiese. “Dynamic Simulation + #' of Management Events for Assessing Impacts of Climate Change on + #' Pre-Alpine Grassland Productivity.” European Journal of Agronomy + #' 128 (August 1, 2021): 126306. + #' https://doi.org/10.1016/j.eja.2021.126306. + #' + #' A refinement to expected yield as function of altitude has been + #' implemented according to Table 1a in + #' Huguenen-Elie et al. "Düngung von Grasland", Agrarforschung Schweiz, + #' 8, (6), 2017, + #' https://www.agrarforschungschweiz.ch/2017/06/9-duengung-von-grasland-grud-2017/ + #' + #' @param DOY Integer day of the year to consider. + #' @param intensity One of ("high", "middle", "low") specifying + #' management intensity. + #' @return target Biomass (kg / ha) that should be reached on day *DOY* + #' for this management *intensity*. + #' + get_target_biomass = function(DOY, intensity = "high") { + # For DOY < 130, use the value at DOY = 130 + DOY = max(130, DOY) + # The relative fraction is taken from Petersen et al. (Fig. S2). + f = self$get_relative_cut_contribution(DOY) + # The absolute annual yield is approximated by Petersen's equation 19, + # refined according to the source they cite. + # The factor 1000 is to convert from t/ha to kg/ha. + gross_yield = self$get_annual_gross_yield(self$MVS$parameters$ELV, + intensity = intensity) * 1000 + target = f * gross_yield + # The target must not be smaller than the remainder after a cut + return(max(target, self$MVS$BM_after_cut)) + }, + + #' @description + #' Relative cut contribution + #' + #' Get the fraction of the total annual harvested biomass that a cut at given + #' *DOY* is expected to contribute. + #' + #' The regression for the target biomass is based on Fig. S2 in the + #' supplementary material of Petersen et al. (2021). + #' + #' @param DOY Integer representing the day of the year on which a cut occurs. + #' + #' @return The fraction (between 0 and 1) of biomass harvested at the cut at + #' given *DOY* divided by the total annual biomass. + #' + #' @references + #' \insertRef{petersen2021DynamicSimulationManagement}{growR} + get_relative_cut_contribution = function(DOY) { + return((-0.1228 * DOY + 48.96) * 1e-2) + }, + + #' @description + #' Last day of cutting season + #' + #' Estimate the last day on which it still makes sense to cut. This is done + #' by checking at which point the expected target biomass (see + #' `self$get_relative_cut_contribution()`) goes below the minimally harvestable + #' standing biomass. + #' + #' @param min_biomass float A standing biomass below this value cannot even + #' be harvested, + #' @param elevation float Altitude in m.a.s.l. + #' @param intensity string Management intensity. One of "high", "middle", "low" + #' + #' @return float Last (fractional) day of the year on which a cut still makes + #' sense. + #' + #' @seealso `get_relative_cut_contribution()` + #' + #' @export + get_end_of_cutting_season = function(min_biomass, elevation, + intensity = "high") { + m_tot = self$get_annual_gross_yield(elevation, intensity) * 1000 + d_end = (min_biomass - m_tot * 48.96e-2) / (m_tot * -0.1228e-2) + return(d_end) + }, + + #' @description + #' Decide based on simple criteria whether day of year *DOY* would be a + #' good day to cut. + #' + #' This follows an implementation described in + #' Petersen, Krischan, David Kraus, Pierluigi Calanca, Mikhail A. + #' Semenov, Klaus Butterbach-Bahl, and Ralf Kiese. “Dynamic Simulation + #' of Management Events for Assessing Impacts of Climate Change on + #' Pre-Alpine Grassland Productivity.” European Journal of Agronomy + #' 128 (August 1, 2021): 126306. + #' https://doi.org/10.1016/j.eja.2021.126306. + #' + #' The decision to cut is made based on two criteria. + #' First, it is checked whether a *target biomass* is reached on given + #' DOY. The defined target depends on the DOY and is given through + #' :func:`get_target_biomass`. If said biomass is present, return `TRUE`. + #' + #' Otherwise, it is checked whether a given amount of time has passed + #' since the last cut. Depending on whether this is the first cut of + #' the season or not, the relevant parameters are + #' :int:`last_DOY_for_initial_cut` and :int:`max_cut_period`. + #' If that amount of time has passed, return `TRUE`, otherwise return + #' `FALSE`. + #' + #' @param DOY Integer day of the year for which to make a cut decision. + #' @return Boolean `TRUE` if a cut happens on day *DOY*. + #' + #' @seealso `get_target_biomass()` + #' + determine_cut = function(DOY) { + # Get a shortcut + MVS = self$MVS + # Don't cut close to end-of-year + if (DOY > self$end_of_cutting_season) { + return(FALSE) + } + cut_target_reached = FALSE + # Check if target biomass is reached. + if (MVS$BM[DOY] >= self$target_biomass[DOY]) { + cut_target_reached = TRUE + } + # Otherwise, check if enough time has passed to warrant a cut. + n_previous_cuts = length(MVS$cut_DOYs) + if (!cut_target_reached) { + if (n_previous_cuts == 0) { + cut_target_reached = DOY > self$last_DOY_for_initial_cut + } else { + last_cut_DOY = MVS$cut_DOYs[n_previous_cuts] + cut_target_reached = DOY - last_cut_DOY > self$max_cut_period + } + } + # If we're not ready to cut, exit. + if (!cut_target_reached) { + return(FALSE) + } + cut_index = n_previous_cuts + 1 + # Otherwise, check if the weather holds up... + if (self$dry_window[DOY] || + # or if we just cannot wait any longer. + self$cut_delays[cut_index] >= self$max_cut_delay) { + # We're cutting. Prepare cut_delays vector for next cut. + self$cut_delays = c(self$cut_delays, 0) + return(TRUE) + } else { + # Wait another day. + self$cut_delays[cut_index] = self$cut_delays[cut_index] + 1 + return(FALSE) + } + } + ) +) + diff --git a/R/modvegesite.R b/R/modvegesite.R index 37165d4..26dc5ff 100644 --- a/R/modvegesite.R +++ b/R/modvegesite.R @@ -75,7 +75,7 @@ initial_state_variables = list( #' ## Initial conditions #' ```{r child = "vignettes/children/initial_conditions.Rmd"} #' ``` -#' @seealso \link[=autocut]{autocut} +#' @seealso \link[=Autocut]{Autocut} #' #' @references #' \insertRef{jouven2006ModelPredictingDynamics}{growR} @@ -127,37 +127,7 @@ ModvegeSite = R6Class( #' occurred during the growth period, in which case reproductive growth is #' stopped. cut_during_growth_preriod = NULL, -#' @field last_DOY_for_initial_cut [autocut] Start cutting after this DOY, -#' even if yield target is not reached. - last_DOY_for_initial_cut = 150, -#' @field max_cut_period [autocut] Maximum period to wait between -#' subsequent cuts. - max_cut_period = 55, -#' @field dry_precipitation_limit [autocut] Maximum amount of allowed -#' precipitation (mm) to consider a day. - dry_precipitation_limit = 1, -#' @field dry_days_before_cut [autocut] Number of days that shold be dry -#' before a cut is made. - dry_days_before_cut = 1, -#' @field dry_days_after_cut [autocut] Number of days that shold be dry -#' after a cut is made. - dry_days_after_cut = 2, -#' @field max_cut_delay [autocut] Number of days a farmer is willing to -#' wait for dry conditions before a cut is made anyways. - max_cut_delay = 5, -#' @field cut_delays [autocut] Vector to keep track of cut delay times. -#' wait for dry conditions before a cut is made anyways. - cut_delays = c(0), -#' @field dry_window [autocut] Logical that indicates if DOY at index is -#' considered dry enough to cut. - dry_window = NULL, -#' @field target_biomass [autocut] Biomass amount that should to be reached -#' by given DOY for a cut to be made. - target_biomass = NULL, -#' @field end_of_cutting_season [autocut] Determined DOY after which no -#' more cuts are made. - end_of_cutting_season = NULL, -#' @field BM_after_cut [autocut] Amount of biomass that remains after a cut +#' @field BM_after_cut Amount of biomass that remains after a cut #' #' (determined through cut_height and biomass densities BDGV, BDDV, BDGR, #' BDDR). BM_after_cut = NULL, @@ -166,8 +136,11 @@ ModvegeSite = R6Class( weather = NULL, #' @field management A list containing management data as returned by #' [ModvegeEnvironment]'s `get_environment_for_year()` method. If its -#' `is_empty` field is `TRUE`, the [autocut] routine will be employed. +#' `is_empty` field is `TRUE`, the [Autocut] routine will be employed. management = NULL, +#' @field Autocut A subclass of [Autocut]. The algorithm used to determine +#' cut events. + Autocut = NULL, #-Public-methods---------------------------------------------------------------- @@ -255,114 +228,6 @@ ModvegeSite = R6Class( return(M$n_cuts > 0 && (DOY %in% M$cut_DOY)) }, - #' @description - #' Decide based on simple criteria whether day of year *DOY* would be a - #' good day to cut. - #' - #' This follows an implementation described in - #' Petersen, Krischan, David Kraus, Pierluigi Calanca, Mikhail A. - #' Semenov, Klaus Butterbach-Bahl, and Ralf Kiese. “Dynamic Simulation - #' of Management Events for Assessing Impacts of Climate Change on - #' Pre-Alpine Grassland Productivity.” European Journal of Agronomy - #' 128 (August 1, 2021): 126306. - #' https://doi.org/10.1016/j.eja.2021.126306. - #' - #' The decision to cut is made based on two criteria. - #' First, it is checked whether a *target biomass* is reached on given - #' DOY. The defined target depends on the DOY and is given through - #' :func:`get_target_biomass`. If said biomass is present, return `TRUE`. - #' - #' Otherwise, it is checked whether a given amount of time has passed - #' since the last cut. Depending on whether this is the first cut of - #' the season or not, the relevant parameters are - #' :int:`last_DOY_for_initial_cut` and :int:`max_cut_period`. - #' If that amount of time has passed, return `TRUE`, otherwise return - #' `FALSE`. - #' - #' @param DOY Integer day of the year for which to make a cut decision. - #' @return Boolean `TRUE` if a cut happens on day *DOY*. - #' - #' @seealso `get_target_biomass()` - #' - determine_cut_automatically = function(DOY) { - # Don't cut close to end-of-year - if (DOY > self$end_of_cutting_season) { - return(FALSE) - } - cut_target_reached = FALSE - # Check if target biomass is reached. - if (self$BM[DOY] >= self$target_biomass[DOY]) { - cut_target_reached = TRUE - } - # Otherwise, check if enough time has passed to warrant a cut. - n_previous_cuts = length(self$cut_DOYs) - if (!cut_target_reached) { - if (n_previous_cuts == 0) { - cut_target_reached = DOY > self$last_DOY_for_initial_cut - } else { - last_cut_DOY = self$cut_DOYs[n_previous_cuts] - cut_target_reached = DOY - last_cut_DOY > self$max_cut_period - } - } - # If we're not ready to cut, exit. - if (!cut_target_reached) { - return(FALSE) - } - cut_index = n_previous_cuts + 1 - # Otherwise, check if the weather holds up... - if (self$dry_window[DOY] || - # or if we just cannot wait any longer. - self$cut_delays[cut_index] >= self$max_cut_delay) { - # We're cutting. Prepare cut_delays vector for next cut. - self$cut_delays = c(self$cut_delays, 0) - return(TRUE) - } else { - # Wait another day. - self$cut_delays[cut_index] = self$cut_delays[cut_index] + 1 - return(FALSE) - } - }, - - #' @description - #' Get target value of biomass on given *DOY*, which determines whether - #' a cut is to occur. - #' - #' The regression for the target biomass is based on Fig. S2 in the - #' supplementary material of - #' Petersen, Krischan, David Kraus, Pierluigi Calanca, Mikhail A. - #' Semenov, Klaus Butterbach-Bahl, and Ralf Kiese. “Dynamic Simulation - #' of Management Events for Assessing Impacts of Climate Change on - #' Pre-Alpine Grassland Productivity.” European Journal of Agronomy - #' 128 (August 1, 2021): 126306. - #' https://doi.org/10.1016/j.eja.2021.126306. - #' - #' A refinement to expected yield as function of altitude has been - #' implemented according to Table 1a in - #' Huguenen-Elie et al. "Düngung von Grasland", Agrarforschung Schweiz, - #' 8, (6), 2017, - #' https://www.agrarforschungschweiz.ch/2017/06/9-duengung-von-grasland-grud-2017/ - #' - #' @param DOY Integer day of the year to consider. - #' @param intensity One of ("high", "middle", "low") specifying - #' management intensity. - #' @return target Biomass (kg / ha) that should be reached on day *DOY* - #' for this management *intensity*. - #' - get_target_biomass = function(DOY, intensity = "high") { - # For DOY < 130, use the value at DOY = 130 - DOY = max(130, DOY) - # The relative fraction is taken from Petersen et al. (Fig. S2). - f = get_relative_cut_contribution(DOY) - # The absolute annual yield is approximated by Petersen's equation 19, - # refined according to the source they cite. - # The factor 1000 is to convert from t/ha to kg/ha. - gross_yield = get_annual_gross_yield(self$parameters$ELV, - intensity = intensity) * 1000 - target = f * gross_yield - # The target must not be smaller than the remainder after a cut - return(max(target, self$BM_after_cut)) - }, - #-Wrapper-function---------------------------------------------------------- #' @description Carry out a ModVege simulation for one year. @@ -391,34 +256,13 @@ ModvegeSite = R6Class( # Determine whether cuts are to be read from input or to be # determined algorithmically. if (management$is_empty) { - # Precalculate target biomass - self$target_biomass = sapply(1:self$days_per_year, - function(day) { - self$get_target_biomass(day, intensity = management$intensity) - } - ) - # Calculate end of the season and max cut interval. - self$end_of_cutting_season = - get_end_of_cutting_season(self$BM_after_cut, self$parameters$ELV, - management$intensity) - n_cuts = get_expected_n_cuts(self$parameters$ELV, management$intensity) - delta = self$end_of_cutting_season - self$last_DOY_for_initial_cut - self$max_cut_period = round(delta / n_cuts) - # Prepare the intervals of dry-enough cut conditions. - dry_days = weather$PP < self$dry_precipitation_limit - self$dry_window = logical(self$days_per_year) - for (i in 1:self$days_per_year) { - i0 = max(i - self$dry_days_before_cut, 1) - i1 = min(i + self$dry_days_after_cut, self$days_per_year) - all_dry = all(dry_days[i0:i1]) - self$dry_window[i] = all_dry - } +# self$Autocut = PetersenAutocut$new(self) + self$Autocut = PhenologicalAutocut$new(self) # Point the cut determination function to the autocut routine. - self$determine_cut = self$determine_cut_automatically + self$determine_cut = self$Autocut$determine_cut } else { self$determine_cut = self$determine_cut_from_input } - # Loop over days of the year. # :TODO: Avoid if j==1 statements by increasing vector sizes by one. logger("Starting loop over days of the year.", level = TRACE) @@ -733,7 +577,6 @@ ModvegeSite = R6Class( # ceased. self$cut_during_growth_preriod = FALSE self$cut_DOYs = c() - self$cut_delays = c(0) # Water balance and growth self[["SENGV"]] = P$SENGV0 diff --git a/R/run.R b/R/run.R index 38936c0..5ced3b4 100644 --- a/R/run.R +++ b/R/run.R @@ -73,7 +73,7 @@ growR_run_loop = function(modvege_environments, output_dir = "", } #-Store-output------------------------------------------------------------ - results[[run]][[i_year]] = modvege$clone(deep = TRUE) + results[[run]][[i_year]] = modvege$clone(deep = FALSE) #-Adjust-initial-conditions-for-next-year--------------------------------- if (!independent) { diff --git a/R/support_functions.R b/R/support_functions.R index d82e20e..1fbfe3a 100644 --- a/R/support_functions.R +++ b/R/support_functions.R @@ -279,135 +279,6 @@ fC.gs = function(x){1} # function(x){1.5/(1 + .5*x/360)} fC.Tt = function(x){0} # function(x){5*(x - 360)/(720 - 360)} fC.ST = function(x){1} # function(x){1 + .1*(x - 360)/(720 - 360)} -#' Lookup table returning expected annual gross yields as function of -#' elevation and management intensity. -#' -#' Based on data from Table 1a in -#' Lookup Table of expected yield as functions of height and management -#' intensity after Olivier Huguenin et al. (2017). -#' -#' @param elevation The elevation of the considered site in meters above sea -#' level. -#' @param intensity One of ("high", "middle", "low", "extensive"). Management -#' intensity for considered site. -#' -#' @return Annual gross yield in t / ha (metric tons per hectare). Note that -#' 1 t/ha = 0.1 kg/m^2. -#' -#' @md -#' @examples -#' get_annual_gross_yield(1200) -#' get_annual_gross_yield(1200, intensity = "low") -#' -#' @references -#' \insertRef{huguenin2017GrundlagenDuengung}{growR} -#' @md -#' @export -get_annual_gross_yield = function(elevation, intensity = "high") { - mask = yield_parameters$intensity == intensity - a = yield_parameters[mask, ]$a - b = yield_parameters[mask, ]$b - return(a + b * max(elevation, 500)) -# return(1) -} - -#' Get number of expected cuts -#' -#' Return the number of expected cuts for a site at a given *elevation* and -#' management *intensity*. -#' -#' This uses data.frame `management_parameters` as a lookup table and -#' interpolates linearly in between the specified values. -#' -#' @param elevation The elevation of the considered site in meters above sea -#' level. -#' @param intensity One of ("high", "middle", "low", "extensive"). Management -#' intensity for considered site. -#' -#' @return Number of expected cuts per season. -#' -#' @md -#' @examples -#' get_expected_n_cuts(1200) -#' get_expected_n_cuts(1200, intensity = "low") -#' -#' @export -get_expected_n_cuts = function(elevation, intensity = "high") { - mask = management_parameters$intensity == intensity - # Find the altitudes below and above given elevation - altitudes = management_parameters[mask, ]$altitude - n_cuts = management_parameters[mask, ]$n_cuts - i0 = which.min(abs(altitudes - elevation)) - alt0 = altitudes[i0] - n0 = n_cuts[i0] - # Check for second closest altitude - i1 = which.min(abs(altitudes[-i0] - elevation)) - alt1 = altitudes[-i0][i1] - n1 = n_cuts[-i0][i1] - # Approximate linearly - n = (n1 - n0) / (alt1 - alt0) * (elevation - alt0) + n0 - return(n) -} - -#' Relative cut contribution -#' -#' Get the fraction of the total annual harvested biomass that a cut at given -#' *DOY* is expected to contribute. -#' -#' The regression for the target biomass is based on Fig. S2 in the -#' supplementary material of Petersen et al. (2021). -#' -#' @param DOY Integer representing the day of the year on which a cut occurs. -#' -#' @return The fraction (between 0 and 1) of biomass harvested at the cut at -#' given *DOY* divided by the total annual biomass. -#' -#' @md -#' @examples -#' get_relative_cut_contribution(1) -#' get_relative_cut_contribution(150) -#' get_relative_cut_contribution(365) -#' # DOYs larger than 365 are insensible -#' get_relative_cut_contribution(600) -#' -#' @references -#' \insertRef{petersen2021DynamicSimulationManagement}{growR} -#' @md -#' @export -get_relative_cut_contribution = function(DOY) { - return((-0.1228 * DOY + 48.96) * 1e-2) -} - -#' Last day of cutting season -#' -#' Estimate the last day on which it still makes sense to cut. This is done -#' by checking at which point the expected target biomass (see -#' [get_relative_cut_contribution()]) goes below the minimally harvestable -#' standing biomass. -#' -#' @param min_biomass float A standing biomass below this value cannot even -#' be harvested, -#' @param elevation float Altitude in m.a.s.l. -#' @param intensity string Management intensity. One of "high", "middle", "low" -#' -#' @return float Last (fractional) day of the year on which a cut still makes -#' sense. -#' -#' @seealso [get_relative_cut_contribution()] -#' -#' @md -#' @examples -#' get_end_of_cutting_season(50, 1200) -#' get_end_of_cutting_season(50, 1200, intensity = "low") -#' -#' @export -get_end_of_cutting_season = function(min_biomass, elevation, - intensity = "high") { - m_tot = get_annual_gross_yield(elevation, intensity) * 1000 - d_end = (min_biomass - m_tot * 48.96e-2) / (m_tot * -0.1228e-2) - return(d_end) -} - #' Create a weighted temperature sum #' #' A temperature sum is constructed by summing the average daily temperature diff --git a/tests/testthat/test-autocut.R b/tests/testthat/test-autocut.R new file mode 100644 index 0000000..2d67e0e --- /dev/null +++ b/tests/testthat/test-autocut.R @@ -0,0 +1,31 @@ +get_example_environment = function(site = "posieux") { + input_dir = system.file("extdata", package = "growR") + param_file = sprintf("%s_parameters.csv", site) + weather_file = sprintf("%s_weather.txt", site) + management_file = sprintf("%s_management1.txt", site) + E = ModvegeEnvironment$new(paste0(site, "1"), + param_file = param_file, + weather_file = weather_file, + management_file = management_file, + input_dir = input_dir + ) +} + +## Test modvegesite$run when no explicit cut dates are provided. +test_that("modvegesite$run() with autocut", { + E = get_example_environment("posieux") + MV = ModvegeSite$new(E$parameters) + year = E$years[[1]] + E1 = E$get_environment_for_year(year) + management = E1$M + # Replace with empty management + management$is_empty = TRUE + management$cut_DOY = NULL + management$cut_years = NULL + management$n_cuts = NULL + for (intensity in c("high", "middle", "low")) { + management$intensity = intensity + expect_no_error(MV$run(year, E1$W, management)) + } +}) +