Kd waterfuns

.Rbuildignore

@@ -5,6 +5,7 @@
 ^Meta$
 ^dev$
 ^DEV$
+^data-raw$
 ^\.github$
 ^codecov\.yml$
 ^LICENSE

.github/workflows/check-standard.yaml

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, development]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, development]
 
 name: R-CMD-check
 

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: OBIC
 Type: Package
 Title: Calculate the Open Bodem Index (OBI) Score
-Version: 2.2.1
+Version: 2.3.0
 Authors@R: c(
     person(given = "Sven", 
            family = "Verweij",

DEV/data/column_description.csv

@@ -103,8 +103,12 @@ I_P_CO;indicator;Ondergrondverdichting;Compaction;-;
 I_P_WRI;indicator;Plant beschikbaar water;Water retention;-;
 I_P_CEC;indicator;Aggregaatstabiliteit;Aggregate stability;-;
 I_P_WO;indicator;Bewerkbaarheid;Workability;-;
-I_B_DI;indicator;Ziektewerendheid;Disease / pest resistance;;
-I_B_SF;indicator;Activiteit bodemleven;Soil life activity;;
+I_B_DI;indicator;Ziektewerendheid;Disease / pest resistance;-;
+I_B_SF;indicator;Activiteit bodemleven;Soil life activity;-;
+I_H_GWR;indicator;Doorlatendheid voor water;Unsaturated permeability;-;
+I_H_NGW;indicator;Efficientie van stikstofmeststoffen in relatie tot kwaliteit grondwater;Nitrogen efficiency for ground water;-;
+I_H_NOW;indicator;Efficientie van stikstofmeststoffen in relatie tot kwaliteit oppervlaktewater;Nitrogen efficiency for surface water;-;
+I_H_PEST;indicator;Uitspoelingsrisico gewasbestrijdingsmiddelen naar grondwater;Leaching risk of pesticides to groundwater;-;
 I_B_NEM;indicator;Plant parasitaire nematoden;Plant parasitic nematodes;;
 I_E_NGW;indicator;N  vasthoudend vermogen voor grondwater;Nitrogen retention for ground water;;
 I_E_NSW;indicator;N vasthoudend vermogen voor oppervlaktewater;Nitrogen retention for surface water;;

DEV/data/weight_obic.csv

@@ -29,6 +29,10 @@ I_M_CLIMATE;indicator;akkerbouw;1;1
 I_M_WATERQUALITY;indicator;akkerbouw;1;1
 I_M_BIODIVERSITY;indicator;akkerbouw;1;1
 I_BCS;indicator;akkerbouw;1;1
+I_H_GWR;indicator;akkerbouw;1;1
+I_H_NGW;indicator;akkerbouw;1;1
+I_H_NOW;indicator;akkerbouw;1;1
+I_H_PEST;indicator;akkerbouw;1;1
 M_GREEN;measurement;akkerbouw;1;1
 M_COMPOST;measurement;akkerbouw;1;1
 M_NONBARE;measurement;akkerbouw;1;1
@@ -74,6 +78,10 @@ I_M_CLIMATE;indicator;mais;1;1
 I_M_WATERQUALITY;indicator;mais;1;1
 I_M_BIODIVERSITY;indicator;mais;1;1
 I_BCS;indicator;mais;1;1
+I_H_GWR;indicator;mais;1;1
+I_H_NGW;indicator;mais;1;1
+I_H_NOW;indicator;mais;1;1
+I_H_PEST;indicator;mais;1;1
 M_GREEN;measurement;mais;1;1
 M_COMPOST;measurement;mais;1;1
 M_NONBARE;measurement;mais;1;1
@@ -119,6 +127,10 @@ I_M_CLIMATE;indicator;grasland;1;1
 I_M_WATERQUALITY;indicator;grasland;1;1
 I_M_BIODIVERSITY;indicator;grasland;1;1
 I_BCS;indicator;grasland;1;1
+I_H_GWR;indicator;grasland;1;1
+I_H_NGW;indicator;grasland;1;1
+I_H_NOW;indicator;grasland;1;1
+I_H_PEST;indicator;grasland;-1;-1
 M_GREEN;measurement;grasland;-1;-1
 M_COMPOST;measurement;grasland;1;1
 M_NONBARE;measurement;grasland;-1;-1
@@ -164,6 +176,10 @@ I_M_CLIMATE;indicator;natuur;1;1
 I_M_WATERQUALITY;indicator;natuur;1;1
 I_M_BIODIVERSITY;indicator;natuur;1;1
 I_BCS;indicator;natuur;1;1
+I_H_GWR;indicator;grasland;-1;-1
+I_H_NGW;indicator;grasland;-1;-1
+I_H_NOW;indicator;grasland;-1;-1
+I_H_PEST;indicator;grasland;-1;-1
 M_GREEN;measurement;natuur;-1;-1
 M_COMPOST;measurement;natuur;1;1
 M_NONBARE;measurement;natuur;1;1

DEV/scripts/ppr_crops_obic.R

@@ -2,11 +2,11 @@
 
   # load package
   library(data.table)
-  
+
   # Load crops.obic (copied csv from pandex)
   cr <- fread('DEV/data/crops_obic.csv')
-  
-  # update the csv and Rdata file
+
+  # update the Rdata file
   crops.obic <- copy(cr)
   save(crops.obic, file = 'data/crops_obic.RData')
 

NAMESPACE

@@ -17,6 +17,7 @@ export(calc_n_efficiency)
 export(calc_nleach)
 export(calc_nlv)
 export(calc_organic_carbon)
+export(calc_permeability)
 export(calc_pesticide_leaching)
 export(calc_ph_delta)
 export(calc_phosphate_availability)
@@ -55,11 +56,13 @@ export(ind_nematodes)
 export(ind_nematodes_list)
 export(ind_nitrogen)
 export(ind_nretention)
+export(ind_permeability)
 export(ind_pesticide_leaching)
 export(ind_ph)
 export(ind_phosphate_availability)
 export(ind_pmn)
 export(ind_potassium)
+export(ind_psp)
 export(ind_resistance)
 export(ind_sealing)
 export(ind_sulfur)

NEWS.md

@@ -1,3 +1,27 @@
+# OBIC 2.3.0 2023-10-16
+
+## Added
+* To `crops.obic`, crop_codes: 6632 (zoete aardappelen), 6636 (naakte haver), 6660 (uien, gele, zaai), and 6664 (uien, rode, zaai)
+* Add `calc_permeability` for calculating top soil permeability 
+* Add `ind_permeability` for calculatin the index for top soil permeability
+* Add vignette `obic_water_functions`, describing the functions that evaluate the effect of soil quality on groundwater quality
+
+## Changed
+* Updated N use norms for onions
+* Update calc_nleach by soil and land use dependent default n supply, issue #154
+* Update calculations in the functions `calc_pesticide_leaching`, `ind_pesticide_leaching`, `calc_psp`, `ind_psp`, `ind_gw_recharge`,`calc_n_efficiency` and `ind_n_efficiency`
+* add soil type and land use dependent NLV estimate in calc_nleach
+
+### Fixed
+* Fixes plotting of regime curve in workability vignette #149
+* Extremely high EOS values for some cultivations fixing #162 & #134
+* Updated bibliography of packages in vignettes
+* N use norm for 1935 (maiskolfsilage)
+* N use norms of 1926 and 1927 (agrarisch natuurmengsel and overige akkerbouwgewassen), increased to match "Akkerbouwgewassen, overig" (RVO, 2022)
+
+### Deprecated
+* building and modifying crops.obic from a script in OBIC
+
 # OBIC 2.2.1 2023-10-16
 ### Added
 * CITATION file and additional author data in the package DESCRIPTION
@@ -43,6 +67,15 @@
 * Update the GitHub Actions to the most recent version
 * Update `RoxygenNote` to v7.2.3
 
+# OBIC 2.0.2 2022-05-13
+### Changed
+* Update calculations in the functions `calc_pesticide_leaching`, `ind_pesticide_leaching`, `calc_psp`, `ind_psp`, `ind_gw_recharge`,`calc_n_efficiency` and `ind_n_efficiency` 
+
+# Added
+* Add `calc_permeability` for calculating top soil permeability 
+* Add `ind_permeability` for calculatin the index for top soil permeability
+* Add vignette `obic_water_functions`, describing the functions that evaluate the effect of soil quality on groundwater quality
+
 # OBIC 2.0.1 2022-03-26
 ### Changed
 * The function `ind_sulpher` is renamed to `ind_sulfur` in documentation

R/groundwater_recharge.R

@@ -1,58 +1,144 @@
 #' Calculate groundwater recharge of a soil
 #' 
-#' This function calculates an index score for groundwater storage based on WHC, sealing risk, drainage and subsoil compaction
+#' This function calculates an index score for groundwater storage based on precipitation surplus, infiltration at saturation, sealing risk, drainage and subsoil compaction
 #' 
 #' 
-#' @param D_WRI_WHC (numeric) The value for Water Retention index (WRI) for WHC as calculated by \code{\link{calc_waterretention}}
+#' @param B_LU_BRP (numeric) The crop code from the BRP
 #' @param D_PSP (numeric) The precipitation surplus per crop  calculated by \code{\link{calc_psp}}
+#' @param D_WRI_K (numeric) The value for top soil permeability (cm/d) as calculated by \code{\link{calc_permeability}}
 #' @param I_P_SE (numeric) The indicator value for soil sealing
 #' @param I_P_CO (numeric) The indicator value for occurrence of subsoil compaction
 #' @param B_DRAIN (boolean) Are drains installed to drain the field (options: yes or no)
+#' @param B_GWL_CLASS (character) The groundwater table class
 #' 
 #' @examples 
-#' ind_gw_recharge(D_WRI_WHC = 0.6, D_PSP = 200, I_P_SE = 0.6, I_P_CO = 0.9)
-#' ind_gw_recharge(D_WRI_WHC = 0.8, D_PSP = 400, I_P_SE = 0.4, I_P_CO = 0.2)
+#' ind_gw_recharge(B_LU_BRP = 265,D_PSP = 200, D_WRI_K = 10, I_P_SE = 0.6, I_P_CO = 0.9, 
+#' B_DRAIN = FALSE, B_GWL_CLASS = 'GtV')
+#' ind_gw_recharge(B_LU_BRP = 233, D_PSP = 400, D_WRI_K = 10, I_P_SE = 0.4, I_P_CO = 0.2, 
+#' B_DRAIN = TRUE, B_GWL_CLASS = 'GtII')
 #'
 #' @return 
 #' The evaluated score for the soil function to improve groundwater recharge. A numeric value between 0 and 1.
 #'          
 #' @export
-ind_gw_recharge <- function(D_WRI_WHC, D_PSP, I_P_SE, I_P_CO, B_DRAIN = FALSE){
+ind_gw_recharge <- function(B_LU_BRP, D_PSP, D_WRI_K, I_P_SE, I_P_CO, B_DRAIN, B_GWL_CLASS){
 
-  I_WRI_WHC = I_W_GWS = cf_compaction = cf_drain = D_I_PSP = NULL
+  I_H_GWR = D_I_WRI_K = cf_compaction = cf_drain = D_I_PSP = NULL
 
   # Check inputs
-  arg.length <- max(length(I_P_CO),length(B_DRAIN))
-  checkmate::assert_numeric(D_WRI_WHC, lower = 0, upper = 1000, any.missing = FALSE)
-  checkmate::assert_numeric(D_PSP, lower = 0, upper = 1000, any.missing = FALSE)
-  checkmate::assert_numeric(I_P_CO, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
-  checkmate::assert_numeric(I_P_SE, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
+  arg.length <- max(length(B_LU_BRP),length(D_WRI_K),length(D_PSP),length(I_P_SE),length(I_P_CO),length(B_DRAIN))
+
+
+  checkmate::assert_numeric(B_LU_BRP, any.missing = FALSE, min.len = 1, len = arg.length)
+  checkmate::assert_subset(B_LU_BRP, choices = unique(OBIC::crops.obic$crop_code), empty.ok = FALSE)
+  checkmate::assert_numeric(D_WRI_K, any.missing = FALSE, len = arg.length)
+  checkmate::assert_numeric(D_PSP, any.missing = FALSE, len = arg.length)
+  checkmate::assert_numeric(I_P_SE, any.missing = FALSE, len = arg.length)
+  checkmate::assert_numeric(I_P_CO, any.missing = FALSE, len = arg.length)
   checkmate::assert_logical(B_DRAIN, any.missing = FALSE, len = arg.length)
+  checkmate::assert_subset(B_GWL_CLASS, choices = c('GtI','GtII','GtIII','GtIV','GtV','GtVI','GtVII','GtVIII',
+                                                    'GtIIb','GtIIIb','GtVb'), empty.ok = FALSE)
 
   # import data into table
-  dt <- data.table(D_WRI_WHC = D_WRI_WHC,
+  dt <- data.table(B_LU_BRP = B_LU_BRP,
+                   D_PSP = D_PSP,
+                   D_WRI_K = D_WRI_K,
                    I_P_SE = I_P_SE,
                    I_P_CO = I_P_CO,
-                   B_DRAIN = B_DRAIN)
+                   B_DRAIN = B_DRAIN,
+                   B_GWL_CLASS = B_GWL_CLASS)
 
   # Calculate the indicator value for water holding capacity
-  dt[,I_WRI_WHC := ind_waterretention(D_WRI_WHC,type = 'water holding capacity')]
+  dt[,D_I_WRI_K := ind_permeability(D_WRI_K)]
 
-  # Calculate indicator value for precipitation surpluss
-  dt[,D_I_PSP := evaluate_logistic(D_PSP,0.05,300,2.5)]
+  # Calculate indicator value for precipitation surplus
+  dt[,D_I_PSP := ind_psp(D_PSP,B_LU_BRP)]
 
   # Correct for subsoil compaction or drainage
-  dt[B_DRAIN == TRUE, c('cf_drain','cf_compaction') := list(0.6,1)]
-  dt[B_DRAIN == FALSE, c('cf_drain','cf_compaction') := list(1,fifelse(I_P_CO < 0.5,0.8,1))]
+  dt[B_DRAIN == TRUE & B_GWL_CLASS %in% c('GtIIIb','GtIV'), c('cf_drain','cf_compaction') := list(0.6,1)]
+  dt[B_DRAIN == FALSE | B_GWL_CLASS %in% c('GtI','GtII','GtIII','GtV','GtVI','GtVII','GtVIII','GtIIb','GtVb'),
+     c('cf_drain','cf_compaction') := list(1,fifelse(I_P_CO < 0.5,0.8,1))]
 
   # Calculate aggregated score
-  dt[,I_W_GWS := (0.35 * I_WRI_WHC + 0.35 * D_I_PSP + 0.3 * I_P_SE) * cf_compaction * cf_drain]
+  dt[,I_H_GWR := (0.7 * D_I_PSP + 0.15 * D_I_WRI_K + 0.15 * I_P_SE) * cf_compaction * cf_drain]
 
   # Calculate mean
-  out <- dt[,I_W_GWS]
+  out <- dt[,I_H_GWR]
 
   return(out)
 }
 
 
 
+#' Calculate the permeability of the top soil
+#' 
+#' This function calculates the permeability of the top soil
+#' 
+#' 
+#' @param A_CLAY_MI (numeric) The clay content of the soil (\%)
+#' @param A_SAND_MI (numeric) The sand content of the soil (\%)
+#' @param A_SILT_MI (numeric) The silt content of the soil (\%)
+#' @param A_SOM_LOI (numeric) The organic matter content of the soil (\%)
+#' 
+#'         
+#' @export
+calc_permeability <- function(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI){
+
+  id = thetaS = thetaR = alfa = n = ksat = l = Pleem = mineral = D_WRI_K = NULL
+
+  # Check inputs
+  arg.length <- max(length(A_CLAY_MI), length(A_SAND_MI),length(A_SILT_MI), length(A_SOM_LOI))
+  checkmate::assert_numeric(A_CLAY_MI, lower = 0, upper = 100, any.missing = FALSE)
+  checkmate::assert_numeric(A_SAND_MI, lower = 0, upper = 100, any.missing = FALSE)
+  checkmate::assert_numeric(A_SILT_MI, lower = 0, upper = 100, any.missing = FALSE)
+  checkmate::assert_numeric(A_SOM_LOI, lower = 0, upper = 100, any.missing = FALSE)
+
+  # Collect data in a table
+  dt <- data.table(
+    id = 1:arg.length,
+    A_CLAY_MI = A_CLAY_MI,
+    A_SAND_MI = A_SAND_MI,
+    A_SILT_MI = A_SILT_MI,
+    A_SOM_LOI = A_SOM_LOI,
+    value = NA_real_
+  )
+
+  # Calculate mineral fractions
+  dt[,mineral   := A_CLAY_MI + A_SAND_MI + A_SILT_MI]
+  dt[,A_CLAY_MI := A_CLAY_MI * 100 / mineral]
+  dt[,A_SILT_MI := A_SILT_MI * 100 / mineral]
+  dt[,Pleem     := A_CLAY_MI + A_SILT_MI]
+
+
+  # Calculate unsaturated permeability
+  dt[,  c("Dichtheid", "thetaR", "thetaS", "alfa", "n", "ksat", "l") := pFpara_ptf_Wosten2001(A_CLAY_MI, Pleem, A_SOM_LOI, 150, 1)]
+
+  # Calculate permeabilty
+  dt[,D_WRI_K := ksat * (((1 + alfa * thetaR^n)^(1 - 1/n) - alfa * thetaR^(n - 1))^2) / ((1 + alfa * thetaR^n)^((1 - 1/n) * (l + 2)))]
+
+  # output
+  return(dt[,D_WRI_K])
+
+}
+
+
+
+#' Calculate the indicator score for the permeability of the top soil
+#' 
+#' This function calculates the indicator score for the permeability of the top soil
+#' 
+#' 
+#' @param D_WRI_K (numeric) The value for top soil permeability (cm/d) as calculated by \code{\link{calc_permeability}}
+#' 
+#'         
+#' @export
+ind_permeability <- function(D_WRI_K){
+
+  # Check inputs
+  checkmate::assert_numeric(D_WRI_K, lower = 0, any.missing = FALSE)
+
+  D_I_WRI_K <- evaluate_logistic(D_WRI_K,0.08,50,0.4)
+
+  return(D_I_WRI_K)
+}
+