Merge pull request #87 from AgroCares/BBWP-67

Bbwp 67

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BBWPC
 Type: Package
 Title: Calculator for BedrijfsBodemWaterPlan (BBWP)
-Version: 2.1.1
+Version: 2.2.1
 Authors@R: c(
     person("Gerard", "Ros", email = "[email protected]", role = c("aut","cre")),
     person("Sven", "Verweij", email = "[email protected]", role = c("aut")),

NEWS.md

@@ -1,3 +1,11 @@
+# BBWPC v2.2.1 2023-12-28
+## Added
+* possibility to falisfy the weighing fraction for evaluation of risks in `wf` function, default set to TRUE, #BBWP-67
+* argument `penalty` added to the functions `bbwp_field_indicators`, `bbwp_field_scores` and `bbwp_main`, #BBWP-67
+
+## Changed
+* function `bbwp_field_indicators` is simplified
+
 # BBWPC v2.1.1 2023-12-27
 ## Added
 * function `bbwp_format_sc_wenr` to ensure that B_SC_WENR is an integer conform format Van den Akker (2006), #BBWP-66

R/bbwp_field_indicators.R

@@ -30,6 +30,7 @@
 #' @param D_WUE_WWRI (numeric) The relative score of soil wetness stress for improved efficiency of water
 #' @param D_WUE_WDRI (numeric) The relative score of drought stress for improved efficiency of water
 #' @param D_WUE_WHC (numeric) The relative score of drought stress for improved efficiency of water
+#' @param penalty (boolean) the option to apply a penalty for high risk BBWP field indicators
 #' 
 #' @import data.table
 #' @import OBIC
@@ -40,10 +41,11 @@ bbwp_field_indicators <- function(D_NGW_SCR,D_NGW_LEA,D_NGW_NLV,
                                   D_NSW_SCR,D_NSW_GWT,D_NSW_RO,D_NSW_SLOPE, D_NSW_WS,D_NSW_NLV,
                                   D_PSW_SCR,D_PSW_GWT,D_PSW_RO,D_PSW_SLOPE,D_PSW_WS,D_PSW_PCC,D_PSW_PSG,D_PSW_PRET,
                                   D_NUE_WRI,D_NUE_PBI,D_NUE_WDRI,D_NUE_NLV,
-                                  D_WUE_WWRI,D_WUE_WDRI,D_WUE_WHC){
+                                  D_WUE_WWRI,D_WUE_WDRI,D_WUE_WHC, penalty = TRUE){
 
   # add visual bindings
   D_RISK_NGW = D_RISK_NSW = D_RISK_PSW = D_RISK_NUE = D_RISK_WB = id = NULL
+  risk_cor = value = group = risk = mcf = WS = SLOPE = NULL
 
   # check length inputs
   arg.length <- max(
@@ -54,86 +56,73 @@ bbwp_field_indicators <- function(D_NGW_SCR,D_NGW_LEA,D_NGW_NLV,
     length(D_WUE_WWRI),length(D_WUE_WDRI),length(D_WUE_WHC)
   )
 
-  # add checks on input
-
-
   # copy input in one data.table
-  dt <- data.table(
-    id = 1:arg.length,
-    D_NGW_SCR = D_NGW_SCR,
-    D_NGW_LEA = D_NGW_LEA,
-    D_NGW_NLV = D_NGW_NLV,
-    D_NSW_SCR = D_NSW_SCR,
-    D_NSW_GWT = D_NSW_GWT,
-    D_NSW_RO = D_NSW_RO,
-    D_NSW_SLOPE = D_NSW_SLOPE,
-    D_NSW_WS = D_NSW_WS,
-    D_NSW_NLV = D_NSW_NLV,
-    D_PSW_SCR = D_PSW_SCR,
-    D_PSW_GWT = D_PSW_GWT,
-    D_PSW_RO = D_PSW_RO,
-    D_PSW_SLOPE = D_PSW_SLOPE,
-    D_PSW_WS = D_PSW_WS,
-    D_PSW_PCC = D_PSW_PCC,
-    D_PSW_PSG = D_PSW_PSG,
-    D_PSW_PRET = D_PSW_PRET,
-    D_NUE_WRI = D_NUE_WRI,
-    D_NUE_PBI = D_NUE_PBI,
-    D_NUE_WDRI = D_NUE_WDRI,
-    D_NUE_NLV = D_NUE_NLV,
-    D_WUE_WWRI = D_WUE_WWRI,
-    D_WUE_WDRI = D_WUE_WDRI,
-    D_WUE_WHC = D_WUE_WHC,
-    D_RISK_NGW = NA_real_,
-    D_RISK_NSW = NA_real_,
-    D_RISK_PSW = NA_real_,
-    D_RISK_NUE = NA_real_,
-    D_RISK_WB = NA_real_
-  )
-
-
-  # integrate all relative field risk indicators into one for indictor for N loss to groundwater
-  dt[, D_RISK_NGW := (wf(D_NGW_SCR) * D_NGW_SCR + 
-                  3 * wf(D_NGW_LEA) * D_NGW_LEA + 
-                  2 * wf(D_NGW_NLV) * D_NGW_NLV) /
-       (wf(D_NGW_SCR) + 3 * wf(D_NGW_LEA) + 2 * wf(D_NGW_NLV))]
-
-  # integrate all relative field risk indicators into one for indictor for N loss to surface water
-  dt[, D_RISK_NSW := (wf(D_NSW_SCR) * D_NSW_SCR + 
-                      wf(D_NSW_GWT) * D_NSW_GWT + 
-                      wf(D_NSW_SLOPE) * D_NSW_SLOPE +
-                      wf(D_NSW_RO) * D_NSW_RO + 
-                      wf(D_NSW_WS) * D_NSW_WS + 
-                  3 * wf(D_NSW_NLV) * D_NSW_NLV ) /
-       (wf(D_NSW_SCR) + wf(D_NSW_GWT) + wf(D_NSW_SLOPE) + wf(D_NSW_RO) + wf(D_NSW_WS) + 3 * wf(D_NSW_NLV))]
-
-  # integrate all relative field risk indicators into one for indictor for P loss to surface water
-  dt[, D_RISK_PSW := (2 * wf(D_PSW_SCR) * D_PSW_SCR + 
-                          wf(D_PSW_GWT) * D_PSW_GWT + 
-                      2 * wf(D_PSW_RO) * D_PSW_RO + 
-                          wf(D_PSW_SLOPE) * D_PSW_SLOPE +
-                      2 * wf(D_PSW_WS) * D_PSW_WS + 
-                          wf(D_PSW_PCC) * D_PSW_PCC + 
-                          wf(D_PSW_PSG) * D_PSW_PSG + 
-                          wf(D_PSW_PRET) * D_PSW_PRET) /
-       (2 * wf(D_PSW_SCR) + wf(D_PSW_GWT) + 2 * wf(D_PSW_RO) + wf(D_PSW_SLOPE) + 2 * wf(D_PSW_WS) + wf(D_PSW_PCC) + wf(D_PSW_PSG) + wf(D_PSW_PRET))]
-
-  # minimize risks when there are no ditches around (wet surrounding fraction < 0.2)
-  dt[D_NSW_WS <= 0.2 & D_PSW_SLOPE < 1,D_RISK_PSW := 0.1]
-  dt[D_NSW_WS <= 0.2 & D_PSW_SLOPE < 1,D_RISK_NSW := 0.1]
-  dt[D_NSW_WS <= 0.1 & D_PSW_SLOPE < 1,D_RISK_PSW := 0.01]
-  dt[D_NSW_WS <= 0.1 & D_PSW_SLOPE < 1,D_RISK_NSW := 0.01]
-
-  # integrate all relative field risk indicators into one for indictor for N and P efficiency of inputs
-  dt[, D_RISK_NUE := (wf(D_NUE_WRI) * D_NUE_WRI + 2 * wf(D_NUE_PBI) * D_NUE_PBI + wf(D_NUE_NLV) * D_NUE_NLV + wf(D_NUE_WDRI) * D_NUE_WDRI) /
-       (wf(D_NUE_WRI) + 2 * wf(D_NUE_PBI) + wf(D_NUE_NLV) + wf(D_NUE_WDRI))]
-
-  # integrate all relative field risk indicators into one for indictor for water retention and efficiency
-  dt[, D_RISK_WB := (wf(D_WUE_WWRI) * D_WUE_WWRI + wf(D_WUE_WDRI) * D_WUE_WDRI + 2 * wf(D_WUE_WHC) * D_WUE_WHC) / (wf(D_WUE_WWRI) + wf(D_WUE_WDRI) + 2 * wf(D_WUE_WHC))]
-
-  # normalise these indicators ???
-
-  setorder(dt, id)
+  dt <- data.table(id = 1:arg.length,
+                   D_NGW_SCR = D_NGW_SCR,
+                   D_NGW_LEA = D_NGW_LEA,
+                   D_NGW_NLV = D_NGW_NLV,
+                   D_NSW_SCR = D_NSW_SCR,
+                   D_NSW_GWT = D_NSW_GWT,
+                   D_NSW_RO = D_NSW_RO,
+                   D_NSW_SLOPE = D_NSW_SLOPE,
+                   D_NSW_WS = D_NSW_WS,
+                   D_NSW_NLV = D_NSW_NLV,
+                   D_PSW_SCR = D_PSW_SCR,
+                   D_PSW_GWT = D_PSW_GWT,
+                   D_PSW_RO = D_PSW_RO,
+                   D_PSW_SLOPE = D_PSW_SLOPE,
+                   D_PSW_WS = D_PSW_WS,
+                   D_PSW_PCC = D_PSW_PCC,
+                   D_PSW_PSG = D_PSW_PSG,
+                   D_PSW_PRET = D_PSW_PRET,
+                   D_NUE_WRI = D_NUE_WRI,
+                   D_NUE_PBI = D_NUE_PBI,
+                   D_NUE_WDRI = D_NUE_WDRI,
+                   D_NUE_NLV = D_NUE_NLV,
+                   D_WUE_WWRI = D_WUE_WWRI,
+                   D_WUE_WDRI = D_WUE_WDRI,
+                   D_WUE_WHC = D_WUE_WHC
+                  )
+
+  # melt the data.table to simplify corrections
+  dt.melt <- data.table::melt(dt, id.vars = 'id',variable.name = 'risk')
+
+  # add correction factor based on risk itself
+  dt.melt[,risk_cor := wf(value,type = "indicators",penalty = penalty)]
+
+  # add groups of risk indicators
+  dt.melt[,group := gsub('_[A-Z]+$','',gsub('D_','',risk))]
+
+  # add manual weighing factor for risks
+  dt.melt[,mcf := 1]
+  dt.melt[group=='NGW' & grepl('_LEA$',risk), mcf := 3]
+  dt.melt[group=='NGW' & grepl('_NLV$',risk), mcf := 2]
+  dt.melt[group=='NSW' & grepl('_NLV$',risk), mcf := 3]
+  dt.melt[group=='PSW' & grepl('_SCR$|_RO$|_WS$',risk), mcf := 2]
+  dt.melt[group=='NUE' & grepl('_PBI$',risk), mcf := 2]
+  dt.melt[group=='WUE' & grepl('_WHC$',risk), mcf := 2]
+
+
+  # minimize risks when there are no ditches around the field (wet surrounding fraction < 0.2)
+
+    # add criteria properties as column (to use as filter)
+    dt.melt[,WS := value[risk=='D_NSW_WS'],by='id']
+    dt.melt[,SLOPE := value[risk=='D_NSW_SLOPE'],by='id']
+
+    # ensure that the final risk after aggregation gets the value 0.1 or 0.01
+    dt.melt[WS <= 0.2 & SLOPE < 1 & group %in% c('NSW','PSW'), c('mcf','risk_cor','value') :=  list(1,1000,0.1)]
+    dt.melt[WS <= 0.1 & SLOPE < 1 & group %in% c('NSW','PSW'), c('mcf','risk_cor','value') :=  list(1,1000,0.01)]
+    dt.melt[,c('WS','SLOPE') := NULL]
+
+  # calculate the mean aggregated risk indicators
+  dt <- dt.melt[,list(risk = sum(risk_cor * value * mcf)/sum(risk_cor * mcf)),by=c('id','group')]
+  dt <- dcast(dt,id~group,value.var='risk')
+
+  # replace output names
+  setnames(dt,old=c('NGW','NSW','NUE','PSW','WUE'),new = c('D_RISK_NGW','D_RISK_NSW','D_RISK_NUE','D_RISK_PSW','D_RISK_WB'))
+
+  # sort output based on id
+  setorder(dt,id)
 
   # extract output
   out <- dt[,mget(c('D_RISK_NGW','D_RISK_NSW','D_RISK_PSW','D_RISK_NUE','D_RISK_WB'))]

R/bbwp_field_scores.R

@@ -23,7 +23,7 @@
 #' @param B_CT_NSW_MAX (numeric) the max critical target for N reduction loss (kg N / ha)
 #' @param measures (data.table) the measures planned / done per fields
 #' @param sector (string) a vector with the farm type given the agricultural sector (options: 'dairy', 'arable', 'tree_nursery', 'bulbs')
-#' 
+#' @param penalty (boolean) the option to apply a penalty for high risk BBWP field indicators 
 #'   
 #' @import data.table
 #'
@@ -32,7 +32,7 @@
 bbwp_field_scores <- function(B_SOILTYPE_AGR, B_GWL_CLASS, A_P_SG, B_SLOPE_DEGREE, B_LU_BBWP,B_AER_CBS,
                               M_DRAIN, D_SA_W, D_RISK_NGW, D_RISK_NSW, D_RISK_PSW, D_RISK_NUE, D_RISK_WB,
                               B_GWP, B_AREA_DROUGHT, B_CT_PSW, B_CT_NSW, 
-                              B_CT_PSW_MAX = 0.5, B_CT_NSW_MAX = 5.0, measures, sector){
+                              B_CT_PSW_MAX = 0.5, B_CT_NSW_MAX = 5.0, measures, sector,penalty = TRUE){
 
   # add visual bindings
   cfngw = cfwb = cfnsw = cfpsw = cfnue = NULL
@@ -182,13 +182,13 @@ bbwp_field_scores <- function(B_SOILTYPE_AGR, B_GWL_CLASS, A_P_SG, B_SLOPE_DEGRE
   dt[,S_BBWP_NUE := 100 * D_OPI_NUE]
   dt[,S_BBWP_WB := 100 * D_OPI_WB]
 
-  dt[,S_BBWP_TOT := (S_BBWP_NGW * wf(S_BBWP_NGW, type="score") + 
-                      S_BBWP_NSW * wf(S_BBWP_NSW, type="score") + 
-                      S_BBWP_PSW * wf(S_BBWP_PSW, type="score") + 
-                      S_BBWP_NUE * wf(S_BBWP_NUE, type="score") + 
-                      S_BBWP_WB * wf(S_BBWP_WB, type="score")) /
-       (wf(S_BBWP_NGW, type="score") + wf(S_BBWP_NSW, type="score") +  wf(S_BBWP_PSW, type="score") +  
-          wf(S_BBWP_NUE, type="score") +  wf(S_BBWP_WB, type="score"))]
+  dt[,S_BBWP_TOT := (S_BBWP_NGW * wf(S_BBWP_NGW, type="score",penalty = penalty) + 
+                      S_BBWP_NSW * wf(S_BBWP_NSW, type="score",penalty = penalty) + 
+                      S_BBWP_PSW * wf(S_BBWP_PSW, type="score",penalty = penalty) + 
+                      S_BBWP_NUE * wf(S_BBWP_NUE, type="score",penalty = penalty) + 
+                      S_BBWP_WB * wf(S_BBWP_WB, type="score",penalty)) /
+       (wf(S_BBWP_NGW, type="score",penalty = penalty) + wf(S_BBWP_NSW, type="score",penalty = penalty) +  wf(S_BBWP_PSW, type="score",penalty = penalty) +  
+          wf(S_BBWP_NUE, type="score",penalty = penalty) +  wf(S_BBWP_WB, type="score",penalty = penalty))]
 
   # order the fields
   setorder(dt, id)

R/bbwp_helpers.R

@@ -33,18 +33,22 @@ cdf_rank <- function(smean,ssd,svalue){
 #' 
 #' @param x The risk or score value to be weighted
 #' @param type Use the weighing function for indicators or score
+#' @param penalty (boolean) the option to apply a penalty for high risk BBWP field indicators
 #' 
 #' @export
-wf <- function(x, type = "indicators") {
+wf <- function(x, type = "indicators", penalty = TRUE) {
 
-  if (type == "indicators") {
+  if (type == "indicators" & penalty == TRUE) {
 
     y <- 1 / (1 - x + 0.2)
 
-  } else if (type == "score") {
+  } else if (type == "score" & penalty == TRUE) {
 
     y <- 1 / (x * 0.01 + 0.2)
 
+  } else if(penalty == FALSE){
+
+    y = 1
   }
 
   return(y)

R/bbwp_main.R

@@ -37,7 +37,8 @@
 #' @param measures (data.table) the measures planned / done per fields
 #' @param sector (string) a vector with the farm type given the agricultural sector (options: options: 'dairy', 'arable', 'tree_nursery', 'bulbs')
 #' @param output (string) a vector specifying the output type of the function. Options: scores, measures 
-#' 
+#' @param penalty (boolean) the option to apply a penalty for high risk BBWP field indicators
+#'  
 #' @details 
 #' B_SLOPE_DEGREE should be used, for backwards compatibility B_SLOPE can still be used. At least one of the must be used, when both are supplied, B_SLOPE is ignored.
 #' LSW is by default a data.table with LSW properties, being calculated from bbwp_lsw_properties. Note that all B_LSW_IDs should be preset in the LSW data.table.
@@ -52,7 +53,7 @@ bbwp <- function(B_SOILTYPE_AGR, B_LU_BBWP,B_GWL_CLASS, B_SC_WENR, B_HELP_WENR,B
                  B_GWP, B_AREA_DROUGHT, B_CT_PSW, B_CT_NSW,B_CT_PSW_MAX = 0.5, B_CT_NSW_MAX = 5.0, 
                  D_SA_W, D_RO_R, B_AREA, 
                  M_DRAIN, B_LSW_ID, LSW = NULL,
-                 measures, sector,output = 'scores'){
+                 measures, sector,output = 'scores',penalty=TRUE){
 
   # add visual binding
   field_id = code = value_min = value_max = NULL
@@ -158,7 +159,8 @@ bbwp <- function(B_SOILTYPE_AGR, B_LU_BBWP,B_GWL_CLASS, B_SC_WENR, B_HELP_WENR,B
                                   D_NUE_NLV = dt$npe_nlv,
                                   D_WUE_WWRI = dt$wue_wwri,
                                   D_WUE_WDRI = dt$wue_wdri,
-                                  D_WUE_WHC = dt$wue_whc
+                                  D_WUE_WHC = dt$wue_whc,
+                                  penalty = penalty
                                 )
 
   # Calculate BBWP field scores
@@ -187,7 +189,8 @@ bbwp <- function(B_SOILTYPE_AGR, B_LU_BBWP,B_GWL_CLASS, B_SC_WENR, B_HELP_WENR,B
                                     B_CT_PSW_MAX = B_CT_PSW_MAX, 
                                     B_CT_NSW_MAX = B_CT_NSW_MAX,
                                     measures = measures,
-                                    sector = sector
+                                    sector = sector,
+                                    penalty = penalty
                                   )
 
   # Calculate the BBWP farm score

tests/testthat/test-bbwp.R

@@ -509,4 +509,5 @@ test <- bbwp(B_SOILTYPE_AGR = c('dekzand', 'loess', 'rivierklei'),
       LSW = NULL
     ))
   })
+