-
Notifications
You must be signed in to change notification settings - Fork 2
/
bodemconditiescore.R
182 lines (159 loc) · 8.48 KB
/
bodemconditiescore.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#' Calculate the BodemConditieScore
#'
#' This function calculates the BodemConditieScore given input from manual observations made in the field.
#' The individual parameters are scored in three classes: poor (0), neutral (1) or good (2)
#' More information on this test can be found \href{https://mijnbodemconditie.nl/}{here}
#'
#'
#' @param B_LU_BRP (numeric) The crop code from the BRP
#' @param B_SOILTYPE_AGR (character) The agricultural type of soil
#' @param A_SOM_LOI (numeric) The percentage organic matter in the soil (\%)
#' @param D_PH_DELTA (numeric) The pH difference with the optimal pH.
#' @param A_EW_BCS (numeric) The presence of earth worms (score 0-1-2)
#' @param A_SC_BCS (numeric) The presence of compaction of subsoil (score 0-1-2)
#' @param A_GS_BCS (numeric) The presence of waterlogged conditions, gley spots (score 0-1-2)
#' @param A_P_BCS (numeric) The presence / occurrence of water puddles on the land, ponding (score 0-1-2)
#' @param A_C_BCS (numeric) The presence of visible cracks in the top layer (score 0-1-2)
#' @param A_RT_BCS (numeric) The presence of visible tracks / rutting or trampling on the land (score 0-1-2)
#' @param A_RD_BCS (integer) The rooting depth (score 0-1-2)
#' @param A_SS_BCS (integer) The soil structure (score 0-1-2)
#' @param A_CC_BCS (integer) The crop cover on the surface (score 0-1-2)
#' @param type (character) Define output of the function. Options: score (integrated score) and indicator (score per indicator)
#'
#' @references \href{https://mijnbodemconditie.nl/}{mijnbodemconditie.nl}
#'
#' @import data.table
#'
#' @examples
#' calc_bcs(B_LU_BRP = 265, B_SOILTYPE_AGR = 'dekzand', A_SOM_LOI = 3.5, D_PH_DELTA = 0.4,
#' A_EW_BCS = 1, A_SC_BCS = 1, A_GS_BCS = 1, A_P_BCS = 1, A_C_BCS = 1, A_RT_BCS =1, A_RD_BCS = 1,
#' A_SS_BCS = 1, A_CC_BCS = 1)
#'
#' @return
#' A visual soil assessment score derived from field observations driven by organic matter content and soil structure properties. Returns a numeric value.
#' @export
calc_bcs <- function(B_LU_BRP,B_SOILTYPE_AGR,A_SOM_LOI, D_PH_DELTA,
A_EW_BCS = NA, A_SC_BCS = NA, A_GS_BCS = NA, A_P_BCS = NA, A_C_BCS = NA, A_RT_BCS = NA, A_RD_BCS = NA,
A_SS_BCS = NA, A_CC_BCS = NA,
type = 'score'
) {
id = crop_code = crop_n = crop_category = bcs_om = bcs_ph = soiltype = soiltype.n = NULL
# Load in the datasets
crops.obic <- as.data.table(OBIC::crops.obic)
setkey(crops.obic, crop_code)
soils.obic <- as.data.table(OBIC::soils.obic)
setkey(soils.obic, soiltype)
# Check input
arg.length <- max(length(A_EW_BCS), length(A_SC_BCS), length(A_GS_BCS), length(A_P_BCS),
length(A_C_BCS), length(A_RT_BCS), length(A_RD_BCS), length(A_SS_BCS),
length(A_CC_BCS), length(B_LU_BRP))
checkmate::assert_numeric(B_LU_BRP, any.missing = FALSE, min.len = 1, len = arg.length)
checkmate::assert_subset(B_LU_BRP, choices = unique(crops.obic$crop_code), empty.ok = FALSE)
checkmate::assert_character(B_SOILTYPE_AGR, any.missing = FALSE, min.len = 1, len = arg.length)
checkmate::assert_subset(B_SOILTYPE_AGR, choices = unique(soils.obic$soiltype), empty.ok = FALSE)
checkmate::assert_numeric(A_SOM_LOI, lower = 0, upper = 100, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(D_PH_DELTA, lower = 0, upper = 5)
checkmate::assert_numeric(A_EW_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_SC_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_GS_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_P_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_C_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_RT_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_RD_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_SS_BCS,lower = 0, upper = 2)
checkmate::assert_numeric(A_CC_BCS,lower = 0, upper = 2)
# Collect data in a table
dt <- data.table(
id = 1:arg.length,
B_LU_BRP = B_LU_BRP,
B_SOILTYPE_AGR = B_SOILTYPE_AGR,
A_SOM_LOI = A_SOM_LOI,
A_EW_BCS = A_EW_BCS,
A_SC_BCS = A_SC_BCS,
A_GS_BCS = A_GS_BCS,
A_P_BCS = A_P_BCS,
A_C_BCS = A_C_BCS,
A_RT_BCS = A_RT_BCS,
A_RD_BCS= A_RD_BCS,
A_SS_BCS = A_SS_BCS,
A_CC_BCS = A_CC_BCS,
value = NA_real_
)
# merge data.table with crop type and soil type
dt <- merge(dt, crops.obic[, list(crop_code, crop_category)], by.x = "B_LU_BRP", by.y = "crop_code")
dt <- merge(dt, soils.obic[, list(soiltype, soiltype.n)], by.x = "B_SOILTYPE_AGR", by.y = "soiltype")
# calculate the scores of the Visual Soil Assessment Form: these range from 0 (poor) to 2 (good)
# Calculate the score for pH
dt[, bcs_ph := round(ind_ph(D_PH_DELTA) * 2)]
# Calculate the score for organic matter
dt[, bcs_om := 1]
# organic matter class 'low' when the OM content is lower than 30% quantile
dt[crop_category == 'akkerbouw' & soiltype.n == 'klei' & A_SOM_LOI < 2.2, bcs_om := 0]
dt[crop_category == 'akkerbouw' & soiltype.n == 'zand' & A_SOM_LOI < 3.0, bcs_om := 0]
dt[crop_category == 'akkerbouw' & soiltype.n == 'loess' & A_SOM_LOI < 2.4, bcs_om := 0]
dt[crop_category == 'akkerbouw' & soiltype.n == 'veen' & A_SOM_LOI < 7.9, bcs_om := 0]
dt[crop_category == 'grasland' & soiltype.n == 'klei' & A_SOM_LOI < 6.8, bcs_om := 0]
dt[crop_category == 'grasland' & soiltype.n == 'zand' & A_SOM_LOI < 4.6, bcs_om := 0]
dt[crop_category == 'grasland' & soiltype.n == 'loess' & A_SOM_LOI < 5.1, bcs_om := 0]
dt[crop_category == 'grasland' & soiltype.n == 'veen' & A_SOM_LOI < 15.5, bcs_om := 0]
dt[crop_category == 'mais' & soiltype.n == 'klei' & A_SOM_LOI < 3.4, bcs_om := 0]
dt[crop_category == 'mais' & soiltype.n == 'zand' & A_SOM_LOI < 3.4, bcs_om := 0]
dt[crop_category == 'mais' & soiltype.n == 'loess' & A_SOM_LOI < 2.6, bcs_om := 0]
dt[crop_category == 'mais' & soiltype.n == 'veen' & A_SOM_LOI < 8.7, bcs_om := 0]
# organic matter class 'high' when the OM content is higher than 70% quantile
dt[crop_category == 'akkerbouw' & soiltype.n == 'klei' & A_SOM_LOI > 3.8, bcs_om := 2]
dt[crop_category == 'akkerbouw' & soiltype.n == 'zand' & A_SOM_LOI > 4.8, bcs_om := 2]
dt[crop_category == 'akkerbouw' & soiltype.n == 'loess' & A_SOM_LOI > 3.3, bcs_om := 2]
dt[crop_category == 'akkerbouw' & soiltype.n == 'veen' & A_SOM_LOI > 14.6, bcs_om := 2]
dt[crop_category == 'grasland' & soiltype.n == 'klei' & A_SOM_LOI > 12.9, bcs_om := 2]
dt[crop_category == 'grasland' & soiltype.n == 'zand' & A_SOM_LOI > 6.6, bcs_om := 2]
dt[crop_category == 'grasland' & soiltype.n == 'loess' & A_SOM_LOI > 7.7, bcs_om := 2]
dt[crop_category == 'grasland' & soiltype.n == 'veen' & A_SOM_LOI > 28.6, bcs_om := 2]
dt[crop_category == 'mais' & soiltype.n == 'klei' & A_SOM_LOI > 6.2, bcs_om := 2]
dt[crop_category == 'mais' & soiltype.n == 'zand' & A_SOM_LOI > 4.8, bcs_om := 2]
dt[crop_category == 'mais' & soiltype.n == 'loess' & A_SOM_LOI > 3.4, bcs_om := 2]
dt[crop_category == 'mais' & soiltype.n == 'veen' & A_SOM_LOI > 20.1, bcs_om := 2]
# evaluation for nature (not present yet)
dt[crop_category == 'natuur', bcs_om := 2]
# Calculate final score
dt[,value := 2 * A_CC_BCS + 3 * A_RD_BCS + 3 * A_SC_BCS + 3 * A_EW_BCS +
3 * A_SS_BCS +3 * bcs_ph + 3 * bcs_om + 1 * A_GS_BCS -
2 * A_P_BCS - 1 * A_C_BCS - 1 * A_RT_BCS]
# Combine both tables and extract values
setorder(dt, id)
if(type == 'score'){
# select the integrative BCS score
value <- dt[, value]
}else{
# define columns that are given as output of the function when indicators are calculated
cols <- colnames(dt)[grepl('_BCS$',colnames(dt))]
# convert the BCS input to a 0-1 scale
dt[,c(cols) := lapply(.SD,function(x) x * 0.5),.SDcols = cols]
# select the columns as output
value <- dt[,mget(cols)]
}
return(value)
}
#' Calculate the indicator for BodemConditieScore
#'
#' This function calculates the final score for the BodemConditieScore by using the scores calculated by \code{\link{calc_bcs}}
#'
#' @param D_BCS (numeric) The value of BCS calculated by \code{\link{calc_bcs}}
#'
#' @examples
#' ind_bcs(D_BCS = 12)
#' ind_bcs(D_BCS = c(12,18,26,30))
#'
#' @return
#' The evaluated score for the Visual Soil Assessment. A numeric value between 0 and 50.
#'
#' @export
ind_bcs <- function(D_BCS) {
# Check inputs
checkmate::assert_numeric(D_BCS, lower = 0, upper = 50)
# Evaluate the Bodem Conditie Score
# ensure that it does not exceed 1
value <- pmin(D_BCS / 40, 1)
# return output
return(value)
}