Merge pull request #26 from springgbv/yuki_sulpher

Yuki sulpher
AgroCares · Oct 21, 2019 · d4283e8 · d4283e8
2 parents 2a5018e + 257d146
commit d4283e8
Show file tree

Hide file tree

Showing 6 changed files with 32 additions and 14 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,13 +5,15 @@
 * The uppper limit for `D_BCS` is increased from 40 to 50
 * Switch on crumbleability
 * For `calc_phosphate_availability` the category `arable` is added for the crop categories
+* Changed evaluation of sulphur for arable fields #26
 
 ### Fixed
 * Fixed typo if mais in `ind_managment`
 * Fixed test for winderodibility
 * Use the correct correction factor in `calc_sealing_risk` #19
 * Fix for calculating `I_P_CEC` #24
 * Fix for calculating difficult values in `calc_sombalance` #25
+* Fix for `calc_sbal_arable` where combinations of soil type and region that do not exist in table 6.2 of Handboek Bodem & Bemesting gave a NA #26
 
 ## Version 0.8.0 2019-08-02
 ### Added

diff --git a/R/evaluate.R b/R/evaluate.R
@@ -53,4 +53,4 @@ evaluate_parabolic <- function(x, x.top) {
 
   return(y)
 
-}
+}
diff --git a/R/obic_preprocessing.R b/R/obic_preprocessing.R
@@ -78,7 +78,7 @@ obic_preprocessing <- function(dt) {
   dt[, D_P_DU := calc_winderodibility(A_CLAY_MI, A_SILT_MI, B_LU_BRP)]
 
   # Calculate the sulphur supply
-  dt[, D_SLV := calc_slv(A_OS_GV,A_S_TOT, B_LU_BRP, B_BT_AK, B_LG_CBS,D_BDS)]
+  dt[, D_SLV := calc_slv(A_S_TOT,A_OS_GV,B_LU_BRP, B_BT_AK, B_LG_CBS,D_BDS)]
 
   # Calculate the magnesium index
   dt[, D_MG := calc_magnesium_availability(A_MG_CC,A_PH_CC,A_OS_GV,A_CEC_CO, A_K_CC,

diff --git a/R/sombalance.R b/R/sombalance.R
@@ -12,7 +12,7 @@
 #' @export
 calc_sombalance <- function(A_OS_GV, A_P_PAL, A_P_WA, B_LU_BRP,M_M3, M_M6) {
 
-  c.diss = id = crop_code = crop_name = crop_n = cropinput = mdose = compost = catchcrop = NULL
+  c.diss = crop_code = crop_name = crop_n = cropinput = mdose = compost = catchcrop = NULL
   crop_eos = crop_eos_residue = NULL
 
   # Load in the datasets
@@ -72,12 +72,14 @@ calc_sombalance <- function(A_OS_GV, A_P_PAL, A_P_WA, B_LU_BRP,M_M3, M_M6) {
   dt[crop_n=='akkerbouw' & A_P_WA > 55,mdose := 0.85 * 50 * slurry_EOS_Pratio]
 
   # EOS input via compost to arable soils
-  dt[crop_n == 'akkerbouw', compost := ifelse(M_M3 == 0, 0, 15 * 218 / M_M3) ]
-  dt[crop_n != 'akkerbouw', compost := 0 ]
+  dt[crop_n == 'akkerbouw' & M_M3 == 0, compost := 0]
+  dt[crop_n == 'akkerbouw' & M_M3 != 0, compost :=  15 * 218 / M_M3]
+  dt[crop_n != 'akkerbouw', compost := 0]
 
   # EOS input via catch crops (and mandatory crops)
-  dt[,catchcrop := ifelse(M_M6,850,0)]
-  dt[grepl('mais|aardappel',crop_name),catchcrop := 850]
+  dt[M_M6 == TRUE, catchcrop := 850]
+  dt[M_M6 != TRUE, catchcrop := 0]
+  dt[grepl('mais|aardappel',crop_name), catchcrop := 850]
 
   # calculate simple eos balance (kg EOS / ha / yr)
   dt[,value := cropinput + mdose + compost + catchcrop - c.diss]

diff --git a/R/sulpher.R b/R/sulpher.R
@@ -12,7 +12,7 @@
 #' @import data.table
 #' 
 #' @export
-calc_slv <- function(A_OS_GV,A_S_TOT, B_LU_BRP, B_BT_AK, B_LG_CBS,D_BDS) {
+calc_slv <- function(A_S_TOT, A_OS_GV, B_LU_BRP, B_BT_AK, B_LG_CBS,D_BDS) {
 
   a = c.ass = c.diss = id = crop_code = soiltype = soiltype.n = crop_category = NULL
   minip.a = D_OC = A_CS_RAT = NULL
@@ -87,7 +87,7 @@ calc_slv <- function(A_OS_GV,A_S_TOT, B_LU_BRP, B_BT_AK, B_LG_CBS,D_BDS) {
 
   # Calculate the SLV for nature land
   dt.nature <- dt[crop_category == "natuur"]
-  dt.nature[,value := 1.5 * A_S_TOT * D_BDS]
+  dt.nature[,value := 1.5 * A_S_TOT * 0.001 * D_BDS * 0.001]
 
   # Combine both tables and extract values
   dt <- rbindlist(list(dt.grass, dt.maize,dt.arable,dt.nature), fill = TRUE)
@@ -147,6 +147,8 @@ calc_sbal_arable <- function(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS) {
   dt[,B_LG_CBS := tolower(B_LG_CBS)]
 
   # add cluster variable to be used later (related to soil type and agronomic region)
+  # Remark YF: Not all B_LG_CBS is covered, resulting in NAs for 'clust'.
+  # For example, dekzand soils in Rivierngebied miss a clust value.
   dt[grepl('klei',B_BT_AK) & grepl('bouwh|oldambt',B_LG_CBS), clust := 1]
   dt[grepl('klei',B_BT_AK) & grepl('rivier|zuidwestelijk',B_LG_CBS), clust := 2]
   dt[grepl('klei',B_BT_AK) & grepl('ijsselmeer',B_LG_CBS), clust := 3]
@@ -155,7 +157,9 @@ calc_sbal_arable <- function(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS) {
   dt[B_BT_AK=='veen', clust := 6]
   dt[grepl('dal|zand|xxx',B_BT_AK) & grepl('noord|oldambt',B_LG_CBS), clust := 7]
   dt[grepl('dal|zand|xxx',B_BT_AK) & grepl('oostelijk|centraal|zuidelijk|zuidwest-brabant',B_LG_CBS), clust := 8]
-  dt[B_BT_AK=='loess',clust := 9] 
+  dt[B_BT_AK=='loess',clust := 9]
+  dt[grepl('klei',B_BT_AK) & is.na(clust), clust := 10]
+  dt[grepl('dal|zand|xxx',B_BT_AK) & is.na(clust), clust := 11]
 
   # add crop S requirement classes  
   dt[,cropclass := calc_cropclass(B_LU_BRP,B_BT_AK,nutrient='S')]
@@ -170,6 +174,9 @@ calc_sbal_arable <- function(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS) {
   dt[clust==7, slv_av := 10]
   dt[clust==8, slv_av := 10]
   dt[clust==9, slv_av := 16]
+  # For combinations that are outside table 6.2 of Handboek Bodem & Bemesting the average slv_av per soiltype is used
+  dt[clust==10, slv_av := mean(c(20, 21, 45, 32, 41))]
+  dt[clust==11, slv_av := mean(c(10, 10))]
 
   # estimate required fertilizer dose
   dt[, sfert := 0]
@@ -184,6 +191,13 @@ calc_sbal_arable <- function(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS) {
   dt[cropclass == 'class2' & clust == 9, sfert := 15]
   dt[cropclass == 'class2' & clust == 2, sfert := 10]
   dt[cropclass == 'class3' & clust %in% c(1,7,8,9), sfert := 10]
+  # For combinations that are outside table 6.2 of Handboek Bodem & Bemesting the average sfert per soiltype is used
+  dt[cropclass == 'class1' & clust == 10, sfert := mean(c(50, 25, 10, 15, 10, 0))]
+  dt[cropclass == 'class1' & clust == 11, sfert := mean(c(55, 50))]
+  dt[cropclass == 'class2' & clust == 10, sfert := mean(c(20, 0, 0, 0, 0, 0))]
+  dt[cropclass == 'class2' & clust == 11, sfert := mean(c(25, 20))]
+  dt[cropclass == 'class3' & clust == 10, sfert := mean(c(10, 0, 0, 0, 0, 0))]
+  dt[cropclass == 'class3' & clust == 11, sfert := mean(c(10, 10))]
 
   # total S requirement (kg S / ha)
   dt[,sreq := slv_av + sfert]
@@ -245,10 +259,10 @@ ind_sulpher <- function(D_SLV,B_LU_BRP, B_BT_AK, B_LG_CBS) {
   # Evaluate S availability for arable land  -----
   dt.arable <- dt[crop_category == "akkerbouw"]
   if(nrow(dt.arable)>0){
-  dt.arable[,sbal := calc_sbal_arable(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS)]
-  dt.arable[,value := evaluate_parabolic(sbal, x.top = 0)]
+    dt.arable[,sbal := calc_sbal_arable(D_SLV, B_LU_BRP, B_BT_AK, B_LG_CBS)]
+    dt.arable[,value := evaluate_logistic(sbal, b = 0.5, x0 = -6, 5)]
   }
-  
+
   # Evaluate S availability for maize land -----
   dt.maize <- dt[crop_category == "mais"]
   dt.maize[,value := evaluate_logistic(D_SLV, b = 5, x0 = 5, v = 5)]

diff --git a/tests/testthat/test-sulpher.R b/tests/testthat/test-sulpher.R
@@ -61,7 +61,7 @@ test_that("ind_sulpher works", {
       B_BT_AK = rep('dekzand',10),
       B_LG_CBS = rep('Oostelijk Veehouderijgebied',10)
     ),
-    expected = c(0,0,1,1,1,1,1,1,1,1),
+    expected = c(0.272,0.779,1,1,1,1,1,1,1,1),
     tolerance = 0.01
   )
   # grasland op klei
Original file line number	Diff line number	Diff line change
Expand Up		@@ -53,4 +53,4 @@ evaluate_parabolic <- function(x, x.top) {

		return(y)

		}
		}