## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE ) options(rmarkdown.html_vignette.check_title = FALSE) ## ----setup, include = FALSE--------------------------------------------------- # load packages library(OBIC); library(data.table); library(ggplot2);require(patchwork) setDTthreads(1) # load binnenveld data.table binnenveld <- as.data.table(OBIC::binnenveld) ## ----showbinnenveld----------------------------------------------------------- dim(binnenveld) binnenveld[1] ## ----echo=FALSE, out.width = '85%', out.height = '85%', fig.cap = 'Figure 1. Graphic representation of how measured soil properties are aggregated to scores.'---- # ![](OBIC_score_integratie.png){widht = 25%, height = 20%} knitr::include_graphics('../vignettes/OBIC_score_integratie_2.png') ## ----------------------------------------------------------------------------- # select the relevant columns without management measures and data from Visual Soil Assessment cols <- colnames(binnenveld)[!grepl('_BCS$|^M_',colnames(binnenveld))] # select the first field, a grassland field dt <- binnenveld[ID==1,mget(cols)] # run the obic_field with default management measures and no visual assessment data # test the obic field function via obic_field and give only the final score obic_field(B_SOILTYPE_AGR = dt$B_SOILTYPE_AGR, B_GWL_CLASS = dt$B_GWL_CLASS, B_SC_WENR = dt$B_SC_WENR, B_HELP_WENR = dt$B_HELP_WENR, B_AER_CBS = dt$B_AER_CBS, B_GWL_GLG = dt$B_GWL_GLG, B_GWL_GHG = dt$B_GWL_GHG, B_GWL_ZCRIT = dt$B_GWL_ZCRIT, B_DRAIN = FALSE, B_FERT_NORM_FR = 1, B_LU_BRP = dt$B_LU_BRP, A_SOM_LOI = dt$A_SOM_LOI, A_SAND_MI = dt$A_SAND_MI, A_SILT_MI = dt$A_SILT_MI, A_CLAY_MI = dt$A_CLAY_MI, A_PH_CC = dt$A_PH_CC, A_N_RT = dt$A_N_RT, A_CN_FR = dt$A_CN_FR, A_S_RT = dt$A_S_RT, A_N_PMN = dt$A_N_PMN, A_P_AL = dt$A_P_AL, A_P_CC = dt$A_P_CC, A_P_WA = dt$A_P_WA, A_CEC_CO = dt$A_CEC_CO, A_CA_CO_PO = dt$A_CA_CO_PO, A_MG_CO_PO = dt$A_MG_CO_PO, A_K_CO_PO = dt$A_K_CO_PO, A_K_CC = dt$A_K_CC, A_MG_CC = dt$A_MG_CC, A_MN_CC = dt$A_MN_CC, A_ZN_CC = dt$A_ZN_CC, A_CU_CC = dt$A_CU_CC, output = 'obic_score') # test the obic field function via obic_field_dt and give only the final score obic_field_dt(dt,output = 'obic_score') ## ----results = FALSE---------------------------------------------------------- # run obic_field to retrieve indicators obic_field_dt(dt, output = 'indicators') # run obic_field to retrieve aggregated scores # for chemistry, biology, fysics, management and environment obic_field_dt(dt, output = 'scores') # the default option is to retrieve all output obic_field_dt(dt, output = 'all') ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE, echo=FALSE----------- # make a data.table for sandy and clay soil dt.test <- data.table(B_SOILTYPE_AGR = c(rep('dekzand',10),rep('rivierklei',10)), A_SOM_LOI = c(seq(0.1,10,length.out = 10), seq(0.1,10,length.out = 10)), A_CLAY_MI = c(rep(5,10),rep(25,10))) # estimate bulk density (D_BDS) dt.test[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI)] # plot output ggplot(data = dt.test, aes(x = A_SOM_LOI, y = D_BDS, group = B_SOILTYPE_AGR, color = B_SOILTYPE_AGR)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Bulk density (kg / m3)') + xlab('Soil organic matter content (%)') + theme(legend.position = c(0.8,0.8)) + ggtitle('Estimate bulk density from SOM and clay content') ## ----results = FALSE---------------------------------------------------------- # estimate soil bulk density dt[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI)] # estimate soil sampling depth dt[, D_RD := calc_root_depth(dt$B_LU_BRP)] # estimate Carbon pool dt[, D_OC := calc_organic_carbon(A_SOM_LOI, D_BDS, D_RD)] # estimate age of grassland dt[,D_GA := calc_grass_age(ID, B_LU_BRP)] # estimate crop rotation fraction for sugar beet dt[, D_CP_SUGARBEET := calc_rotation_fraction(ID, B_LU_BRP, crop = "sugarbeet")] ## ----results = FALSE---------------------------------------------------------- # estimate nitrogen supply (kg N / ha) dt[, D_NLV := calc_nlv(B_LU_BRP, B_SOILTYPE_AGR, A_N_RT, A_CN_FR, D_OC, D_BDS, D_GA)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two land uses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate default properties needed to estimate NLV dt.test[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI)] dt.test[, D_RD := calc_root_depth(B_LU_BRP)] dt.test[, D_OC := calc_organic_carbon(A_SOM_LOI, D_BDS, D_RD)] dt.test[, D_GA := calc_grass_age(ID, B_LU_BRP)] # estimate nitrogen supply (kg N / ha) dt.test[, D_NLV := calc_nlv(B_LU_BRP, B_SOILTYPE_AGR, A_N_RT, A_CN_FR, D_OC, D_BDS, D_GA)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged NLV for both soils dt.test1 <- dt.test[,list(D_NLV = mean(D_NLV)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_NLV, x = luse, fill = luse)) + geom_col() + theme_bw() + ylim(0,200)+ ylab('N supplying capacity (kg N / ha)') + xlab('Landuse') + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(legend.position = c(0.2,0.85), plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('N supplying capacity') # estimate NLV for grassland soil over range of N-total levels in soil dt.test2 <- dt.test[ID==1] # overwrite N levels of the soil and the soil texture to sand with a range to illustrate impact of N total on NLV dt.test2[,B_LU_BRP := 265] dt.test2[, A_N_RT := seq(100,5000,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_GA := 5] dt.test2[, D_NLV := calc_nlv(B_LU_BRP, B_SOILTYPE_AGR, A_N_RT, A_CN_FR, D_OC, D_BDS, D_GA)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_N_RT, y = D_NLV)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('N supplying capacity (kg N / ha)') + xlab('Total Nitrogen content (mg / kg)') + ggtitle('Relationship NLV (grasland) and N total') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate phosphate availability index (unitless) dt[, D_PBI := calc_phosphate_availability(B_LU_BRP, A_P_AL, A_P_CC, A_P_WA)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grass and continue maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate PBI dt.test[, D_PBI := calc_phosphate_availability(B_LU_BRP, A_P_AL, A_P_CC, A_P_WA)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged PBI for both soils dt.test1 <- dt.test[,list(D_PBI = mean(D_PBI)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_PBI, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('P availability index (-)') + xlab('Landuse') + ylim(0,7)+ theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('P availability index') # estimate PBI for grassland soil over range of P-CaCl2 levels in soil dt.test2 <- dt.test[ID==1] # overwrite P-CaCl2 levels dt.test2[, A_P_CC := seq(0.5,10,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_PBI := calc_phosphate_availability(B_LU_BRP, A_P_AL, A_P_CC, A_P_WA)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_P_CC, y = D_PBI)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('P availability index (-)') + xlab('Available P-CaCl2 (mg / kg)') + ggtitle('Relationship PBI (grasland) and P-CaCl2') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate potassium availability index (unitless) dt[, D_K := calc_potassium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_K_CC)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate K-availability index dt.test[, D_K := calc_potassium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_K_CC)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged K-availability indices for both soils dt.test1 <- dt.test[,list(D_K = mean(D_K)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_K, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('K availability index (-)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('K availability index') # estimate K-supply for grassland soil over range of K-CaCl2 levels in soil dt.test2 <- dt.test[ID==1] # overwrite K-CaCl2 levels dt.test2[, A_K_CC := seq(5,500,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_K := calc_potassium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_K_CC)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_K_CC, y = D_K)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('K availability index (-)') + xlab('Available K-CaCl2 (mg / kg)') + ggtitle('Relationship K-availability (grasland) and K-CaCl2') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate S supplying capacity dt[, D_SLV := calc_slv(B_LU_BRP, B_SOILTYPE_AGR, B_AER_CBS,A_SOM_LOI,A_S_RT, D_BDS)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate bulk density dt.test[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI)] # estimate Mg-availability index dt.test[,D_SLV := calc_slv(B_LU_BRP, B_SOILTYPE_AGR, B_AER_CBS,A_SOM_LOI,A_S_RT, D_BDS)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged SLV for both soils dt.test1 <- dt.test[,list(D_SLV = mean(D_SLV)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_SLV, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('S supplying capacity (kg S / ha)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('S supplying capacity (kg S / ha)') # estimate SLV for grassland soil over range of total S levels in soil dt.test2 <- dt.test[ID==1] # overwrite A_S_RT levels dt.test2[, A_S_RT := seq(500,5000,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_SLV := calc_slv(B_LU_BRP, B_SOILTYPE_AGR, B_AER_CBS,A_SOM_LOI,A_S_RT, D_BDS)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_S_RT, y = D_SLV)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('S supplying capacity (kg S / ha)') + xlab('total S (mg / kg)') + ggtitle('Relationship SLV (grasland) and total S') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate magnesium availability index (unitless) dt[, D_MG := calc_magnesium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_MG_CC, A_K_CC)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate Mg-availability index dt.test[, D_MG := calc_magnesium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_MG_CC, A_K_CC)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged Mg-availability indices for both soils dt.test1 <- dt.test[,list(D_MG = mean(D_MG)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_MG, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Mg availability index (-)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Mg availability index') # estimate Mg-supply for grassland soil over range of Mg-CaCl2 levels in soil dt.test2 <- dt.test[ID==1] # overwrite Mg-CaCl2 levels dt.test2[, A_MG_CC := seq(5,500,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_MG := calc_magnesium_availability(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, A_CEC_CO, A_K_CO_PO, A_MG_CC, A_K_CC)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_MG_CC, y = D_MG)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('Mg availability index (-)') + xlab('Available Mg-CaCl2 (mg / kg)') + ggtitle('Relationship Mg-availability (grasland) and Mg-CaCl2') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate Cu availability index dt[, D_CU := calc_copper_availability(B_LU_BRP, A_SOM_LOI, A_CLAY_MI, A_K_CC, A_MN_CC, A_CU_CC)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate Cu-availability index dt.test[, D_CU := calc_copper_availability(B_LU_BRP, A_SOM_LOI, A_CLAY_MI, A_K_CC, A_MN_CC, A_CU_CC)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged Cu-availability indices for both soils dt.test1 <- dt.test[,list(D_CU = mean(D_CU)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_CU, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Cu availability index (-)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Cu availability index') # estimate Cu-supply for grassland soil over range of Cu-CaCl2 levels in soil dt.test2 <- dt.test[ID==1] # overwrite Cu-CaCl2 levels dt.test2[, A_CU_CC := seq(5,500,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_CU := calc_copper_availability(B_LU_BRP, A_SOM_LOI, A_CLAY_MI, A_K_CC, A_MN_CC, A_CU_CC)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_CU_CC, y = D_CU)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('Cu availability index (-)') + xlab('Available Cu-CaCl2 (ug / kg)') + ggtitle('Relationship Cu-availability (grasland) and Cu-CaCl2') # plot side by side p1 + p2 ## ----results = FALSE---------------------------------------------------------- # estimate Zn availability index dt[, D_ZN := calc_zinc_availability(B_LU_BRP, B_SOILTYPE_AGR, A_PH_CC, A_ZN_CC)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # estimate Zn-availability index dt.test[, D_ZN := calc_zinc_availability(B_LU_BRP, B_SOILTYPE_AGR, A_PH_CC, A_ZN_CC)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged Zn-availability indices for both soils dt.test1 <- dt.test[,list(D_ZN = mean(D_ZN)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_ZN, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Zn availability index (-)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Zn availability index') # estimate Zn-supply for grassland soil over range of Zn-CaCl2 levels in soil dt.test2 <- dt.test[ID==1] # overwrite Zn-CaCl2 levels dt.test2[, A_ZN_CC := seq(50,5000,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_ZN := calc_zinc_availability(B_LU_BRP, B_SOILTYPE_AGR, A_PH_CC, A_ZN_CC)] # plot output p2 <- ggplot(data = dt.test2, aes(x = A_ZN_CC, y = D_ZN)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('Zn availability index (-)') + xlab('Available Zn-CaCl2 (ug / kg)') + ggtitle('Relationship Zn-availability (grasland) and Zn-CaCl2') # plot side by side p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # estimate distance to required pH # dt[, D_PH_DELTA := calc_ph_delta(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, D_CP_STARCH, # D_CP_POTATO, D_CP_SUGARBEET, D_CP_GRASS, D_CP_MAIS, D_CP_OTHER)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # define a field with two landuses: grassland and maize dt.test <- binnenveld[ID %in% c(1,11)] # Calculate the crop rotation fraction dt.test[, D_CP_STARCH := calc_rotation_fraction(ID, B_LU_BRP, crop = "starch")] dt.test[, D_CP_POTATO := calc_rotation_fraction(ID, B_LU_BRP, crop = "potato")] dt.test[, D_CP_SUGARBEET := calc_rotation_fraction(ID, B_LU_BRP, crop = "sugarbeet")] dt.test[, D_CP_GRASS := calc_rotation_fraction(ID, B_LU_BRP, crop = "grass")] dt.test[, D_CP_MAIS := calc_rotation_fraction(ID, B_LU_BRP, crop = "mais")] dt.test[, D_CP_OTHER := calc_rotation_fraction(ID, B_LU_BRP, crop = "other")] # estimate delta-pH dt.test[, D_PH_DELTA := calc_ph_delta(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, D_CP_STARCH, D_CP_POTATO, D_CP_SUGARBEET, D_CP_GRASS, D_CP_MAIS, D_CP_OTHER)] # add group for figure dt.test[,luse := fifelse(ID==11,'arable','grasland')] # estimate averaged delta-pH for both soils dt.test1 <- dt.test[,list(D_PH_DELTA = mean(D_PH_DELTA)), by = 'luse'] # plot output p1 <- ggplot(data = dt.test1,aes(y = D_PH_DELTA, x = luse, fill = luse)) + geom_col(show.legend = FALSE) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('desired pH-change (-)') + xlab('Landuse') + theme(plot.title = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('desired pH-change (-)') # estimate desired pH-change for grassland soil over range of pH levels in soil dt.test2 <- dt.test[ID==1] # overwrite initial pH dt.test2[, A_PH_CC := seq(3,10,length.out = .N)] dt.test2[, B_SOILTYPE_AGR := 'dekzand'] dt.test2[, D_PH_DELTA := calc_ph_delta(B_LU_BRP, B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI, A_PH_CC, D_CP_STARCH, D_CP_POTATO, D_CP_SUGARBEET, D_CP_GRASS, D_CP_MAIS, D_CP_OTHER)] # plot output p2 <- ggplot(data = dt.test2, aes(y = D_PH_DELTA, x = A_PH_CC)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylab('desired pH-change (-)') + xlab('pH value') + ggtitle('Relationship desired pH-change and pH') # plot side by side p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # estimate importance of CEC supporting crop development # dt[, D_CEC := calc_cec(A_CEC_CO)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE, eval = FALSE---- # # Make a data table with CEC levels between 1 and 100 mmol+/kg # dt <- data.table(A_CEC_CO = seq(1,100,1)) # # # Add D_CEC # dt[,D_CEC := calc_cec(A_CEC_CO)] # # # plot output # p1 <- ggplot(data = dt,aes(y = D_CEC, x = A_CEC_CO)) + # geom_point(show.legend = FALSE) + geom_line()+ theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ # ylab('CEC evaluation') + xlab('Cation exchange capacity') + # theme(plot.title = element_text(size=10), # axis.text = element_text(size=10)) + # ggtitle('') # # p1 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # reformat GWL_CLASS # dt[, B_GWL_CLASS := format_gwt(B_GWL_CLASS)] # # # estimate risk on yield reduction to drought stress or wetness stress # dt[, D_WSI_DS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'droughtstress')] # dt[, D_WSI_WS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'wetnessstress')] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # subset a series of agricultural fields dt.test <- binnenveld[ID %in% c(1,181,8,125,5,11)] # set ID to 1:5 for plotting purpose dt.test[,pID := .GRP, by = ID] # reformat GWL_CLASS dt.test[, B_GWL_CLASS := format_gwt(B_GWL_CLASS)] # estimate delta-pH dt.test[, D_WSI_DS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'droughtstress'),by=ID] dt.test[, D_WSI_WS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'wetnessstress'),by=ID] # estimate averaged moisture stress for the different fields dt.test1 <- dt.test[,list(D_WSI_DS = mean(D_WSI_DS), D_WSI_WS = mean(D_WSI_WS)), by = 'pID'] # melt dt.test1 dt.test1 <- melt(dt.test1,id.vars = 'pID', variable.name = 'stresstype') # plot output p1 <- ggplot(data = dt.test1,aes(y = value, x = pID, fill = stresstype)) + geom_bar(stat='identity') + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylim(0,55)+ theme(legend.position = c(0.3,0.8), plot.title = element_text(size=10), legend.text = element_text(size=10)) + ylab('yield reduction (%)') + xlab('Field') + ggtitle('Soil water stress index') # estimate moisture stress with variable Gt for an arable field dt.test2 <- dt.test[ID==11][1] dt.test2 <- dt.test2[rep(1,4)] dt.test2[,B_GWL_CLASS := c('GtIII','GtIV','GtV','GtVI')] # update WSI dt.test2[, D_WSI_DS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'droughtstress')] dt.test2[, D_WSI_WS := calc_waterstressindex(B_HELP_WENR, B_LU_BRP, B_GWL_CLASS, WSI = 'wetnessstress')] # estimate averaged moisture stress for the different fields dt.test2 <- dt.test2[,list(D_WSI_DS = mean(D_WSI_DS), D_WSI_WS = mean(D_WSI_WS)), by = 'B_GWL_CLASS'] # melt dt.test2 dt.test2 <- melt(dt.test2,id.vars = 'B_GWL_CLASS', variable.name = 'stresstype') # plot output p2 <- ggplot(data = dt.test2,aes(y = value, x = B_GWL_CLASS, fill = stresstype)) + geom_bar(stat='identity') + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylim(0,30)+ theme(legend.position = c(0.3,0.8), plot.title = element_text(size=10), legend.text = element_text(size=10)) + ylab('yield reduction (%)') + xlab('GWL_CLASS') + ggtitle('Soil water stress index') # plot side by side p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # estimate the presence / risk for surface sealing and wind erodibility # dt[, D_SE := calc_sealing_risk(A_SOM_LOI, A_CLAY_MI)] # dt[, D_WE := calc_winderodibility(B_LU_BRP, A_CLAY_MI, A_SILT_MI)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # make a data.table for sandy and clay soil with potatoes dt.test1 <- data.table(B_LU_BRP = rep(859,20), A_CLAY_MI = c(seq(1,25,length.out = 10), seq(1,25,length.out = 10)), A_SOM_LOI = c(rep(3.5,10),rep(3.5,10)), A_SILT_MI = c(rep(5,10),rep(5,10))) # make a data.table for sandy and clay soil with grassland dt.test2 <- data.table(B_LU_BRP = rep(265,20), A_CLAY_MI = c(seq(1,25,length.out = 10), seq(1,25,length.out = 10)), A_SOM_LOI = c(rep(5,10),rep(5,10)), A_SILT_MI = c(rep(5,10),rep(5,10))) # combine both dt.test <- rbind(dt.test1,dt.test2) dt.test[, B_SOILTYPE_AGR := fifelse(A_CLAY_MI < 15,'dekzand','rivierklei')] # estimate the presence / risk for surface sealing and wind erodibility dt.test[, D_SE := calc_sealing_risk(A_SOM_LOI, A_CLAY_MI)] dt.test[, D_WE := calc_winderodibility(B_LU_BRP, A_CLAY_MI, A_SILT_MI)] # set land use dt.test[,luse := fifelse(B_LU_BRP==265,'grass','cropland')] # plot output sealing p1 <- ggplot(data = dt.test, aes(x = A_CLAY_MI, y = D_SE, color = luse,group = luse)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Soil Sealing (-)') + xlab('Clay content (%)') + ylim(4,12)+ theme(legend.position = c(0.4,0.8), plot.title = element_text(size=10), legend.text = element_text(size=10)) + ggtitle('Soil sealing given SOM and clay content') + guides(color = guide_legend(title="land use")) # plot output wind erodibility p2 <- ggplot(data = dt.test, aes(x = A_CLAY_MI, y = D_WE, color = luse)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Wind Erodibility (-)') + xlab('Clay content (%)') + ylim(0,1)+ theme(legend.position = c(0.7,0.8), plot.title = element_text(size=10), legend.text = element_text(size=10)) + ggtitle('Wind erodibility given soil texture and landuse') + guides(color = guide_legend(title="land use")) # plot side by side p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # assess the crumbleability and aggregate stability # dt[, D_CR := calc_crumbleability(A_SOM_LOI, A_CLAY_MI,A_PH_CC)] # dt[, D_AS := calc_aggregatestability(B_SOILTYPE_AGR,A_SOM_LOI,A_K_CO_PO,A_CA_CO_PO,A_MG_CO_PO)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # make a data.table for crumbleability dt.test1 <- data.table(A_CLAY_MI = c(seq(1,25,length.out = 10), seq(1,25,length.out = 10)), A_SOM_LOI = c(rep(3.5,10),rep(3.5,10)), A_PH_CC = rep(6,20)) # make a data.table for aggregate stability dt.test2 <- data.table(B_SOILTYPE_AGR = rep('dekzand', 8), A_SOM_LOI = rep(3.5,8), A_K_CO_PO = seq(28,0,-4), A_CA_CO_PO = seq(30,100,10), A_MG_CO_PO = seq(28,0,-4)) # make second data.table for aggregate stability with clay soil dt.test3 <- data.table(B_SOILTYPE_AGR = rep('rivierklei', 8), A_SOM_LOI = rep(3.5,8), A_K_CO_PO = seq(28,0,-4), A_CA_CO_PO = seq(30,100,10), A_MG_CO_PO = seq(28,0,-4)) dt2 <- rbindlist(list(dt.test2, dt.test3)) # Make a random aggragate stabilty data table with variying cec values dt4 <- data.table(B_SOILTYPE_AGR = rep('rivierklei', 460), A_SOM_LOI = rep(3.5,460), A_CA_CO_PO = rep(seq(50,95,1),10), A_K_CO_PO = c(runif(460, min = 0, max = 10))) dt4 <- dt4[A_K_CO_PO+ A_CA_CO_PO >= 100, A_K_CO_PO := 100 - A_CA_CO_PO] dt4 <- dt4[,A_MG_CO_PO := 100 - A_CA_CO_PO - A_K_CO_PO] # crumbleability and aggregatestability dt.test1[, D_CR := calc_crumbleability(A_SOM_LOI, A_CLAY_MI,A_PH_CC)] dt4 <- dt4[, D_AS := calc_aggregatestability(B_SOILTYPE_AGR,A_SOM_LOI,A_K_CO_PO,A_CA_CO_PO,A_MG_CO_PO)] # plot output crumbleabilty p1 <- ggplot(data = dt.test1, aes(x = A_CLAY_MI, y = D_CR)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Crumbleability') + xlab('Clay content (%)') + ylim(4,12)+ theme(legend.position = c(0.4,0.8), plot.title = element_text(size=10), legend.text = element_text(size=10)) + ggtitle('Soil crumbleability given SOM and clay content') # plot output aggregate stability p4 <- ggplot(data = dt4, aes(x = A_CA_CO_PO, y = D_AS, colour = A_MG_CO_PO)) + geom_point(size=2) + theme_bw() + scale_fill_viridis_d()+ scale_colour_viridis_b()+ ylab('Distance to optimum CEC occupation') + xlab('Calcium ocuppation') + ylim(0,1)+ theme(plot.title = element_text(size=10), legend.position = c(0.35,0.75), legend.text = element_text(size=10), legend.title = element_text(size=10)) + ggtitle('Aggregate stability\ngiven variation in Ca, Mg and K')+ guides(colour = guide_legend(title="Mg occupation (%)")) # plot side by side p1 + p4 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # overwrite soil physical functions for compaction when BCS is available # dt[,D_P_CO := (3 * A_EW_BCS + 3 * A_SC_BCS + 3 * A_RD_BCS - 2 * A_P_BCS - A_RT_BCS)/18] # dt[,D_P_CO := pmax(0, D_P_CO)] # dt[,I_P_CO := fifelse(is.na(D_P_CO),I_P_CO,D_P_CO)] # # # overwrite soil physical functions for aggregate stability when BCS is available # dt[,D_P_CEC := (3 * A_EW_BCS + 3 * A_SS_BCS - A_C_BCS)/12] # dt[,D_P_CEC := pmax(0, D_P_CEC)] # dt[,I_P_CEC := fifelse(is.na(D_P_CEC),I_P_CEC,D_P_CEC)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # make a data table with some BCS values dt <- binnenveld[ID == 1] dt[,year := 1:.N, by=ID] dt <- dt[year==1] # reformat B_SC_WENR dt[, B_SC_WENR := format_soilcompaction(B_SC_WENR)] # calculate soil compaction and aggregate stability dt[, I_P_CO := ind_compaction(B_SC_WENR)] dt[, D_AS := calc_aggregatestability(B_SOILTYPE_AGR,A_SOM_LOI,A_K_CO_PO,A_CA_CO_PO,A_MG_CO_PO)] dt[, I_P_CEC:= ind_aggregatestability(D_AS)] # make three duplicates dt <- dt[rep(1,3)][,id:=.I] # set BCS all from bad quality to good quality # colnames bcscols <- names(dt)[grep('_BCS$',names(dt))] # set BCS dt[id==1, c(bcscols) := 0][,c('A_P_BCS','A_C_BCS','A_RT_BCS') := 2] dt[id==2, c(bcscols) := 1][,c('A_P_BCS','A_C_BCS','A_RT_BCS') := 1] dt[id==3, c(bcscols) := 2][,c('A_P_BCS','A_C_BCS','A_RT_BCS') := 0] # overwrite soil physical functions for compaction when BCS is available dt[,D_P_CO := (3 * A_EW_BCS + 3 * A_SC_BCS + 3 * A_RD_BCS - 2 * A_P_BCS - A_RT_BCS)/18] dt[,D_P_CO := pmax(0.05, D_P_CO)] dt[,I_P_CO2 := fifelse(is.na(D_P_CO),I_P_CO,D_P_CO)] # overwrite soil physical functions for aggregate stability when BCS is available dt[,D_P_CEC := (3 * A_EW_BCS + 3 * A_SS_BCS - A_C_BCS)/12] dt[,D_P_CEC := pmax(0.05, D_P_CEC)] dt[,I_P_CEC2 := fifelse(is.na(D_P_CEC),I_P_CEC,D_P_CEC)] # melt dt2 <- melt(dt[,.(id,I_P_CEC,I_P_CO,I_P_CEC2,I_P_CO2)], id.vars = c('id')) dt2[,treatment := fifelse(grepl('2$',variable),'corrected','uncorrected')] p1 <- ggplot(data = dt2[grepl('CEC',variable)], aes(y = value, x = id-1, fill = treatment)) + geom_bar(stat='identity',position=position_dodge())+ theme_bw() + scale_fill_viridis_d()+ ylab('Soil compaction index') + xlab('Fieldscore VSA')+ theme(legend.position = 'none', plot.title = element_text(size=8), legend.text = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Soil compaction before and after\ncorrection with VSA data\nfor a soil with low compaction') p2 <- ggplot(data = dt2[grepl('CO',variable)], aes(y = value, x = id-1, fill = treatment)) + geom_bar(stat='identity',position=position_dodge())+ theme_bw() + scale_fill_viridis_d()+ ylab('Aggregate stability index') + xlab('Fieldscore VSA')+ theme(legend.position = c(0.3,0.8), plot.title = element_text(size=8), axis.text = element_text(size=10), legend.text = element_text(size=8) ) + ggtitle('Aggregate stability before and\nafter correction withVSA data\nfor a soil with a moderate aggregate stability') # print both figures p1+p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # estimate the plant available water in topsoil # dt[, D_WRI := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'plant available water')] # # # estimate the water holding capacity in topsoil # dt[, D_WRI := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'water holding capacity')] # # # estimate the moisture content of the wilting point in topsoil # dt[, D_WRI := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'wilting point')] # # # estimate the moisture content of the field capacity in topsoil # dt[, D_WRI := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'field capacity')] # # # estimate the saturated permeability in topsoil # dt[, D_WRI := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'Ksat')] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # Create data dt <- data.table(A_CLAY_MI = c(seq(1,36,2.5), rep(5,45)), plot = c(rep('klei',15), rep('zand',15), rep('silt', 15), rep('som',15)), A_SAND_MI = c(rep(20,15), seq(1,36,2.5), rep(20,30)), A_SILT_MI = c(rep(20,30), seq(1,36,2.5), rep(20,15)), A_SOM_LOI = c(rep(2, 45), seq(0.1,5.7,0.4)) ) # Calculate water retention values dt <- dt[, paw := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'plant available water')] dt <- dt[, whc := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'water holding capacity')] dt <- dt[, wp := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'wilting point')] dt <- dt[, fc := calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'field capacity')] dt <- dt[, ksat:= calc_waterretention(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI,type = 'Ksat')] pclay <- ggplot(data = dt[plot == 'klei'], aes(x = A_CLAY_MI, y = paw)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylim(50,80)+ ylab('PAW (mm)') + xlab('Clay content (%)') + ggtitle('Relation between PAW and clay') pcsand <- ggplot(data = dt[plot == 'zand'], aes(x = A_SAND_MI, y = paw)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylim(50,80)+ ylab('PAW (mm)') + xlab('Sand content (%)') + ggtitle('Relation between PAW and sand') pcsilt <- ggplot(data = dt[plot == 'silt'], aes(x = A_SILT_MI, y = paw)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylim(50,80)+ ylab('PAW (mm)') + xlab('Silt content (%)') + ggtitle('Relation between PAW and silt') pcsom <- ggplot(data = dt[plot == 'som'], aes(x = A_SOM_LOI, y = paw)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10))+ ylim(50,80)+ ylab('PAW (mm)') + xlab('Soil organic matter content (%)') + ggtitle('Relation between PAW and SOM') # compose figure in 2 x 2 layout (pclay + pcsand) / (pcsilt + pcsom) ## ----results = FALSE, eval = FALSE-------------------------------------------- # # calculate the index for the potential mineralizable nitrogen pool # dt[, D_PMN := calc_pmn(B_LU_BRP, B_SOILTYPE_AGR, A_N_PMN)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE ,warning=FALSE---- # define test data table with increasing PMN and different land uses dt.test <- data.table(B_LU_BRP = c(rep(265,12), rep(234,12)), B_SOILTYPE_AGR = 'rivierklei', A_N_PMN = rep(seq(0,550,50),2), A_N_PMN2 = rep(seq(0,220,20),2)) dt.test2 <- data.table(B_LU_BRP = c(rep(265,12), rep(234,12)), B_SOILTYPE_AGR = 'dekzand', A_N_PMN = rep(seq(0,550,50),2), A_N_PMN2 = rep(seq(0,220,20),2)) dt.test <- rbindlist(list(dt.test, dt.test2)) # Add unique identifyer dt.test[B_LU_BRP == 265 & B_SOILTYPE_AGR == 'rivierklei',landuse := 'grassland on clay'] dt.test[B_LU_BRP == 265 & B_SOILTYPE_AGR == 'dekzand',landuse := 'grassland on sand'] dt.test[B_LU_BRP == 234 & B_SOILTYPE_AGR == 'rivierklei',landuse := 'arable on clay'] dt.test[B_LU_BRP == 234 & B_SOILTYPE_AGR == 'dekzand',landuse := 'arable on sand'] # Calculate PMN index dt.test <- dt.test[,D_PMN := calc_pmn(B_LU_BRP, B_SOILTYPE_AGR, A_N_PMN)] # plot output p1 <- ggplot(data = dt.test, aes(y = D_PMN, x = A_N_PMN, group = landuse, colour = landuse)) + geom_line() + geom_point(size=3) + ylim(0,800)+ theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=8), legend.title = element_text(size=8), axis.text = element_text(size=10), legend.position = c(0.33,0.75))+ ylab('PMN index') + xlab('Microbial activity') + ggtitle('Relationship between measured PMN and\nadjusted PMN index') # calculate indicator with A_PMN2 (which has lower numbers) dt.test <- dt.test[, I_PMN2 := ind_pmn(calc_pmn(B_LU_BRP, B_SOILTYPE_AGR, A_N_PMN2))] # Plot indicator p2 <- ggplot(data = dt.test[A_N_PMN2<=150], aes(y = I_PMN2, x = A_N_PMN2, group = landuse, colour = landuse)) + geom_line() + geom_point(size=3) + xlim(0,150)+ theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=8), legend.title = element_text(size=8), axis.text = element_text(size=10), legend.position = c(0.6,0.3))+ ylab('PMN indicator score') + xlab('Measured Microbial activity') + ggtitle('Conversion of measured PMN to OBI score') # plot side by side p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # Calculate disease resistance # dt[, I_B_DI := ind_resistance(A_SOM_LOI)] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # Make data table dt <- data.table(A_SOM_LOI = seq(0.1, 6.9,0.4)) # Calculate disease resistance indicator dt[, I_B_DI := ind_resistance(A_SOM_LOI)] # make a plot p1 <- ggplot(data = dt,aes(y = I_B_DI, x = A_SOM_LOI)) + geom_line() + geom_point(size=3) + theme_bw() + scale_fill_viridis_d() + scale_color_viridis_d()+ theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10), legend.position = c(0.75,0.15))+ ylab('Disease resistance indicator score') + xlab('Soil organic matter content (%)') + ggtitle('Relationship between SOM and disease resistance indicator') # plot figure p1 ## ----hidden text nematodes, eval = FALSE, include = FALSE, echo=FALSE--------- # # hide text in unevaluated code block till text is finished # ### Nematodes (WIP) # Nematodes are small animals occurring in all ecotypes on Earth, so also in soil. There is a large variety of nematodes in soils, they vary amongst other things, in size, feeding habits, reproduction speed, and life cycle. Plant parasitic nematodes (PPN) can reproduce in a range of host plants. Some nematodes have a specific host preference while other can reproduce in a large variety of plants. PPN can severely depress crop yields in vulnerable crops. Therefore, farmers can have their fields sampled to analyse which and how many nematodes occur. # # A nematode parameter that is not entered is assumed to be 0. The severity of an infection depends on the number of individuals and differs per parameter. For example: five Pratylenchus fallax (A_RLN_PR_FAL) will barely reduce the indicator score while five Ditylenchus destructor (A_SN_DI_DES) will severly reduce the score. It is assumed that the most severe infection is most limiting for crop production, therefore, the lowest score for a nematode parameter will determine the indicator score. ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ dt <- data.table(nempar = c('A_RLN_PR_FAL', 'A_SN_DI_DES'), nr = 5, nem_ind = c(ind_nematodes(265, A_RLN_PR_FAL=5),ind_nematodes(265, A_SN_DI_DES=5))) ## ----results = FALSE, eval = FALSE-------------------------------------------- # # calculate potential for N leaching to groundwater # dt[,D_NGW := calc_nleach(B_SOILTYPE_AGR, B_LU_BRP, B_GWL_CLASS, D_NLV, B_AER_CBS, leaching_to = "gw")] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # subset a series of agricultural fields dt.test <- binnenveld[ID %in% c(1,2,3,4,5)] dt.test[,year := 1:.N, by = ID] dt.test <- dt.test[year==1] # check B_GWL_class # estimate default properties needed to estimate NLV dt.test[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI),by=ID] dt.test[, D_RD := calc_root_depth(B_LU_BRP),by=ID] dt.test[, D_OC := calc_organic_carbon(A_SOM_LOI, D_BDS, D_RD),by=ID] dt.test[, D_GA := calc_grass_age(ID, B_LU_BRP),by=ID] # estimate nitrogen supply (kg N / ha) dt.test[, D_NLV := calc_nlv(B_LU_BRP, B_SOILTYPE_AGR, A_N_RT, A_CN_FR, D_OC, D_BDS, D_GA),by=ID] # estimate potential for N leaching dt.test[, D_NGW := calc_nleach(B_SOILTYPE_AGR, B_LU_BRP, B_GWL_CLASS, D_NLV, B_AER_CBS, leaching_to = "gw"),by=ID] # estimate potential for N runoff dt.test[, D_NSW := calc_nleach(B_SOILTYPE_AGR, B_LU_BRP, B_GWL_CLASS, D_NLV, B_AER_CBS, leaching_to = "ow"),by=ID] # plot output p1 <- ggplot(data = dt.test, aes(y = D_NGW, x = ID, fill = ID)) + geom_col() + theme_bw() + scale_fill_viridis_b() + theme(legend.position = 'none') + xlab('Field') + ggtitle('Potential for N leaching') p2 <- ggplot(data = dt.test, aes(y = D_NSW, x = ID, fill = ID)) + geom_col() + theme_bw() + scale_fill_viridis_b() + theme(legend.position = 'none') + xlab('Field') + ggtitle('Potential for N runoff') # Make example with different gwt's dt.test2 <- binnenveld[c(1,1,1,1,1)] dt.test2 <- dt.test2[,B_GWL_CLASS := c('GtIII','GtIV', 'GtV', 'GtVI', 'GtVII')] # estimate default properties needed to estimate NLV dt.test2[, D_BDS := calc_bulk_density(B_SOILTYPE_AGR, A_SOM_LOI, A_CLAY_MI)] dt.test2[, D_RD := calc_root_depth(B_LU_BRP)] dt.test2[, D_OC := calc_organic_carbon(A_SOM_LOI, D_BDS, D_RD)] dt.test2[, D_GA := calc_grass_age(ID, B_LU_BRP)] # estimate nitrogen supply (kg N / ha) dt.test2[, D_NLV := calc_nlv(B_LU_BRP, B_SOILTYPE_AGR, A_N_RT, A_CN_FR, D_OC, D_BDS, D_GA)] # estimate potential for N leaching dt.test2[, D_NGW := calc_nleach(B_SOILTYPE_AGR, B_LU_BRP, B_GWL_CLASS, D_NLV, B_AER_CBS, leaching_to = "gw")] # estimate potential for N runoff dt.test2[, D_NSW := calc_nleach(B_SOILTYPE_AGR, B_LU_BRP, B_GWL_CLASS, D_NLV, B_AER_CBS, leaching_to = "ow")] p3 <- ggplot(data = dt.test2, aes(y = D_NGW, x = B_GWL_CLASS, fill = B_GWL_CLASS)) + geom_col() + theme_bw() + scale_fill_viridis_d() + theme(legend.position = 'none') + xlab('B_GWL_CLASS') + ggtitle('Potential for N leaching') p4 <- ggplot(data = dt.test2, aes(y = D_NSW, x = B_GWL_CLASS, fill = B_GWL_CLASS)) + geom_col() + theme_bw() + scale_fill_viridis_d() + theme(legend.position = 'none') + xlab('B_GWL_CLASS') + ggtitle('Potential for N runoff') # plot figures p3 + p4 ## ----fig.width = 7, fig.height = 2.2,fig.fullwidth = TRUE,echo=FALSE, eval = TRUE---- # Table of field properties dt.table <- copy(dt.test2) setnames(dt.table, 'ID', 'field') # estimate N supplying capacity dt.table <- dt.table[,D_NLV := round(D_NLV, 1)] # select only relevant columns dt.table <- dt.table[,.(field, soiltype = B_SOILTYPE_AGR, crop = B_LU_BRP, groundwaterclass = B_GWL_CLASS, nlv = round(D_NLV), regio = B_AER_CBS)] knitr::kable(dt.table, caption = 'Examples of Field properties to calculate N leaching and runoff') ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # create example data dt.test <- data.table(ngw = seq(0,30,2)) dt.test <- dt.test[, I_NGW := ind_nretention(ngw, 'gw')] # plot output p <- ggplot(data = dt.test, aes(y = I_NGW, x = ngw)) + geom_line() + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ xlab('Potential N leaching') # Plot figure p ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # print figures p2 and p5 p1 + p2 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # Calculate organic matter balance # dt[, D_SOM_BAL := calc_sombalance(B_LU_BRP,A_SOM_LOI, A_P_AL, A_P_WA, M_COMPOST, M_GREEN)] # # # add management when input is missing # cols <- c('M_GREEN', 'M_NONBARE', 'M_EARLYCROP','M_COMPOST','M_SLEEPHOSE','M_DRAIN','M_DITCH','M_UNDERSEED', # 'M_LIME', 'M_NONINVTILL', 'M_SSPM', 'M_SOLIDMANURE','M_STRAWRESIDUE','M_MECHWEEDS','M_PESTICIDES_DST') # dt[, c(cols) := add_management(ID,B_LU_BRP, B_SOILTYPE_AGR, # M_GREEN, M_NONBARE, M_EARLYCROP,M_COMPOST,M_SLEEPHOSE,M_DRAIN,M_DITCH,M_UNDERSEED, # M_LIME, M_NONINVTILL, M_SSPM, M_SOLIDMANURE,M_STRAWRESIDUE,M_MECHWEEDS,M_PESTICIDES_DST)] # # # calculate grass age # dt[, D_GA := calc_grass_age(ID, B_LU_BRP)] # # # Calculate the crop rotation fraction # dt[, D_CP_STARCH := calc_rotation_fraction(ID, B_LU_BRP, crop = "starch")] # dt[, D_CP_POTATO := calc_rotation_fraction(ID, B_LU_BRP, crop = "potato")] # dt[, D_CP_SUGARBEET := calc_rotation_fraction(ID, B_LU_BRP, crop = "sugarbeet")] # dt[, D_CP_GRASS := calc_rotation_fraction(ID, B_LU_BRP, crop = "grass")] # dt[, D_CP_MAIS := calc_rotation_fraction(ID, B_LU_BRP, crop = "mais")] # dt[, D_CP_OTHER := calc_rotation_fraction(ID, B_LU_BRP, crop = "other")] # dt[, D_CP_RUST := calc_rotation_fraction(ID, B_LU_BRP, crop = "rustgewas")] # dt[, D_CP_RUSTDEEP := calc_rotation_fraction(ID, B_LU_BRP, crop = "rustgewasdiep")] # # # Calculate series of management actions # dt[, D_MAN := calc_management(A_SOM_LOI,B_LU_BRP, B_SOILTYPE_AGR,B_GWL_CLASS, # D_SOM_BAL,D_CP_GRASS,D_CP_POTATO,D_CP_RUST,D_CP_RUSTDEEP,D_GA, # M_COMPOST,M_GREEN, M_NONBARE, M_EARLYCROP, M_SLEEPHOSE, M_DRAIN, # M_DITCH, M_UNDERSEED, M_LIME, M_NONINVTILL, M_SSPM, # M_SOLIDMANURE,M_STRAWRESIDUE,M_MECHWEEDS,M_PESTICIDES_DST # )] ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # subset binnenveld for first 5 fields dt <- binnenveld[ID %in% c(1:5)] dt <- dt[,ID := as.factor(ID)] # calculate SOM balance for all crops in crop rotation plan dt[, D_SOM_BAL := calc_sombalance(B_LU_BRP,A_SOM_LOI, A_P_AL, A_P_WA, M_COMPOST, M_GREEN)] # calculate the sum per field dt2 <- dt[,list(D_SOM_BAL = sum(D_SOM_BAL)), by = ID] # plot figure p2 <- ggplot(data = dt2,aes(y = D_SOM_BAL, x = ID, fill = ID)) + geom_col() + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('SOM balance') + xlab('Field') + theme(legend.position = c('none'), plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Soil organic matter balance for five fields') # add management when input is missing cols <- c('M_GREEN', 'M_NONBARE', 'M_EARLYCROP','M_COMPOST','M_SLEEPHOSE','M_DRAIN','M_DITCH','M_UNDERSEED', 'M_LIME', 'M_NONINVTILL', 'M_SSPM', 'M_SOLIDMANURE','M_STRAWRESIDUE','M_MECHWEEDS','M_PESTICIDES_DST') dt[, c(cols) := add_management(ID,B_LU_BRP, B_SOILTYPE_AGR, M_GREEN, M_NONBARE,M_EARLYCROP,M_COMPOST, M_SLEEPHOSE,M_DRAIN,M_DITCH,M_UNDERSEED, M_LIME, M_NONINVTILL, M_SSPM,M_SOLIDMANURE, M_STRAWRESIDUE,M_MECHWEEDS,M_PESTICIDES_DST)] # calculate grass age dt[, D_GA := calc_grass_age(ID, B_LU_BRP)] # Calculate the crop rotation fraction dt[, D_CP_STARCH := calc_rotation_fraction(ID, B_LU_BRP, crop = "starch")] dt[, D_CP_POTATO := calc_rotation_fraction(ID, B_LU_BRP, crop = "potato")] dt[, D_CP_SUGARBEET := calc_rotation_fraction(ID, B_LU_BRP, crop = "sugarbeet")] dt[, D_CP_GRASS := calc_rotation_fraction(ID, B_LU_BRP, crop = "grass")] dt[, D_CP_MAIS := calc_rotation_fraction(ID, B_LU_BRP, crop = "mais")] dt[, D_CP_OTHER := calc_rotation_fraction(ID, B_LU_BRP, crop = "other")] dt[, D_CP_RUST := calc_rotation_fraction(ID, B_LU_BRP, crop = "rustgewas")] dt[, D_CP_RUSTDEEP := calc_rotation_fraction(ID, B_LU_BRP, crop = "rustgewasdiep")] # calculate management per field per year dt[, D_MAN := calc_management(A_SOM_LOI,B_LU_BRP, B_SOILTYPE_AGR,B_GWL_CLASS, D_SOM_BAL,D_CP_GRASS,D_CP_POTATO,D_CP_RUST,D_CP_RUSTDEEP,D_GA, M_COMPOST,M_GREEN, M_NONBARE, M_EARLYCROP, M_SLEEPHOSE, M_DRAIN, M_DITCH, M_UNDERSEED, M_LIME, M_NONINVTILL, M_SSPM, M_SOLIDMANURE,M_STRAWRESIDUE,M_MECHWEEDS,M_PESTICIDES_DST )] # calculate mean management score per field dt2 <- dt[,list(D_MAN = mean(D_MAN)), by = ID] # make plot p3 <- ggplot(data = dt2,aes(y = D_MAN, x = ID, fill = ID)) + geom_col() + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Management score') + xlab('Field') + theme(legend.position = c('none'), plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10)) + ggtitle('Management') # print figures p2 + p3 ## ----results = FALSE, eval = FALSE-------------------------------------------- # # evaluate measures # dt.measure <- OBIC::obic_evalmeasure(dt.score, extensive = FALSE) # # # make recommendations of top 3 measures # out.recom <- OBIC::obic_recommendations(dt.measure) ## ----results = FALSE, eval = FALSE-------------------------------------------- # # Be aware these functions below are not given here to be evaluated. # # They illustrate how the different indices can be evaluated (do not execute them here) # # # Calculate indicators for soil chemical functions # dt[, I_C_N := ind_nitrogen(D_NLV, B_LU_BRP)] # dt[, I_C_P := ind_phosphate_availability(D_PBI)] # dt[, I_C_K := ind_potassium(D_K,B_LU_BRP,B_SOILTYPE_AGR,A_SOM_LOI)] # dt[, I_C_MG := ind_magnesium(D_MG, B_LU_BRP, B_SOILTYPE_AGR)] # dt[, I_C_S := ind_sulfur(D_SLV, B_LU_BRP, B_SOILTYPE_AGR, B_AER_CBS)] # dt[, I_C_PH := ind_ph(D_PH_DELTA)] # dt[, I_C_CEC := ind_cec(D_CEC)] # dt[, I_C_CU := ind_copper(D_CU,B_LU_BRP)] # dt[, I_C_ZN := ind_zinc(D_ZN)] # # # Calculate indicators for soil physical functions # dt[, I_P_CR := ind_crumbleability(D_CR, B_LU_BRP)] # dt[, I_P_SE := ind_sealing(D_SE, B_LU_BRP)] # dt[, I_P_DS := ind_waterstressindex(D_WSI_DS)] # dt[, I_P_WS := ind_waterstressindex(D_WSI_WS)] # dt[, I_P_DU := ind_winderodibility(D_WE)] # dt[, I_P_CO := ind_compaction(B_SC_WENR)] # dt[, I_P_WRI := ind_waterretention(D_WRI)] # dt[, I_P_CEC := ind_aggregatestability(D_AS)] # dt[, I_P_WO := ind_workability(D_WO)] # # # Calculate indicators for soil biological functions # dt[, I_B_DI := ind_resistance(A_SOM_LOI)] # dt[, I_B_SF := ind_pmn(D_PMN)] # # # Calculate indicators for environment # dt[, I_E_NGW := ind_nretention(D_NGW, leaching_to = "gw")] # dt[, I_E_NSW := ind_nretention(D_NSW, leaching_to = "ow")] ## ----plot relation function and index values,fig.width = 7, fig.height = 16,fig.fullwidth = TRUE,echo=FALSE---- # create waterstress index relation figure # dt <- data.table(D_WRI = seq(0,100,1)) # dt[, I_P_WRI := ind_waterretention(D_WRI)] # gg <- ggscatter(dt, x = 'D_WRI', y = 'I_P_WRI', title = 'Waterstress evaluation') +xlab('Indicator value') + ylab('Indicator score') +theme_bw() # gg # define a series with numbers reflecting variation in derivatives dt <- data.table(id = 1:202) # add all derivatives dt[,D_PH_DELTA := rep(seq(0,2.5,0.025),2)] dt[,c('D_NLV','D_MG','D_CEC','D_ZN','D_PMN') := rep(seq(0,100,1),2)] dt[,c('D_PBI','D_K','D_SLV') := rep(seq(0,50,0.5),2)] dt[,c('D_CU') := rep(seq(0,20,.2),2)] dt[,c('D_WSI_DS','D_WSI_WS','D_WRI') := rep(seq(0,50,.5),2)] dt[,c('D_NGW','D_NSW') := rep(seq(0,50,.5),2)] dt[,c('D_WE','D_AS') := seq(0,1,length.out = 202)] dt[,D_CR := rep(seq(0,10,0.1),2)] dt[,D_SE := rep(seq(0,50,0.5),2)] dt[,B_LU_BRP := c(rep(265,101),rep(256,101))] dt[,B_SOILTYPE_AGR := rep('dekzand',202)] dt[,A_SOM_LOI := rep(3.5,202)] dt[,D_DUMMY_SOM_LOI := rep(seq(0,35,0.35),2)] dt[,B_AER_CBS := rep("Centraal Veehouderijgebied",202)] dt[,B_SC_WENR := "Matig"] dt[,D_WO := rep(seq(0,1,0.01),2)] # Calculate indicators for soil chemical functions dt[, I_C_N := ind_nitrogen(D_NLV, B_LU_BRP)] dt[, I_C_P := ind_phosphate_availability(D_PBI)] dt[, I_C_K := ind_potassium(D_K,B_LU_BRP,B_SOILTYPE_AGR,A_SOM_LOI)] dt[, I_C_MG := ind_magnesium(D_MG, B_LU_BRP, B_SOILTYPE_AGR)] dt[, I_C_S := ind_sulfur(D_SLV, B_LU_BRP, B_SOILTYPE_AGR, B_AER_CBS)] dt[, I_C_PH := ind_ph(D_PH_DELTA)] dt[, I_C_CEC := ind_cec(D_CEC)] dt[, I_C_CU := ind_copper(D_CU,B_LU_BRP)] dt[, I_C_ZN := ind_zinc(D_ZN)] # Calculate indicators for soil physical functions dt[, I_P_CR := ind_crumbleability(D_CR, B_LU_BRP)] dt[, I_P_SE := ind_sealing(D_SE, B_LU_BRP)] dt[, I_P_DS := ind_waterstressindex(D_WSI_DS)] dt[, I_P_WS := ind_waterstressindex(D_WSI_WS)] dt[, I_P_DU := ind_winderodibility(D_WE)] dt[, I_P_CO := ind_compaction(B_SC_WENR)] dt[, I_P_WRI := ind_waterretention(D_WRI)] dt[, I_P_CEC := ind_aggregatestability(D_AS)] dt[, I_P_WO := ind_workability(D_WO,B_LU_BRP)] # Calculate indicators for soil biological functions dt[, I_B_DI := ind_resistance(D_DUMMY_SOM_LOI)] dt[, I_B_SF := ind_pmn(D_PMN)] # Calculate indicators for environment dt[, I_E_NGW := ind_nretention(D_NGW, leaching_to = "gw")] dt[, I_E_NSW := ind_nretention(D_NSW, leaching_to = "ow")] indcols <- colnames(dt)[grep('^I_',colnames(dt))] cols <- colnames(dt)[grep('^I_|^B_|^A_|^id$',colnames(dt))] dd.melt <- melt.data.table(dt, id.vars = cols,variable.name = 'funct', value.name = 'function_value') dd.melt <- melt.data.table(dd.melt, id.vars = c(cols[!cols%in%indcols],'funct', 'function_value'), measure.vars = indcols, variable.name = 'indicator', value.name = 'index_value') # Reset b_lu_brp to something readable dd.melt <- dd.melt[, B_LU_BRP := as.character(fifelse(B_LU_BRP == 265,'grass','arable'))] # Select correct rows from dd.melt dd.melt <- dd.melt[(funct == 'D_NLV'& indicator == 'I_C_N')| (funct == 'D_PBI'& indicator == 'I_C_P')| (funct == 'D_K'& indicator == 'I_C_K')| (funct == 'D_MG'& indicator == 'I_C_MG')| (funct == 'D_SLV'& indicator == 'I_C_S')| (funct == 'D_PH_DELTA'& indicator == 'I_C_PH')| (funct == 'D_CEC'& indicator == 'I_C_CEC')| (funct == 'D_CU'& indicator == 'I_C_CU')| (funct == 'D_ZN'& indicator == 'I_C_ZN')| (funct == 'D_WSI_DS'& indicator == 'I_P_DS')| (funct == 'D_WSI_WS'& indicator == 'I_P_WS')| (funct == 'D_WRI'& indicator == 'I_P_WRI')| (funct == 'D_CEC'& indicator == 'I_P_CEC')| (funct == 'D_PMN'& indicator == 'I_B_SF')| (funct == 'D_NGW'& indicator == 'I_E_NGW')| (funct == 'D_NSW'& indicator == 'I_E_NSW')| (funct == 'D_CR'& indicator == 'I_P_CR')| (funct == 'D_SE'& indicator == 'I_P_SE')| (funct == 'D_WE'& indicator == 'I_P_DU')| (funct == 'D_DUMMY_SOM_LOI'& indicator == 'I_B_DI')] # make plot ggsp <- ggplot(dd.melt, aes(x= function_value, y= index_value,color=B_LU_BRP)) + geom_line(size = 1.5) + facet_wrap(facets = 'indicator', ncol = 4, scales = 'free_x') + theme_bw() + scale_color_manual(values = c("#440154FF", "#FDE725FF"))+ theme(legend.position = 'bottom') #plot ggsp + labs(caption = 'Relation between function values and their corresponding evaluated index values.') ## ----fig.width = 7, fig.height = 4,fig.fullwidth = TRUE,echo=FALSE------------ # create example data dt.arabl <- data.table(landuse = c(rep('arable on clay',5),rep('arable on peat',5),rep('arable on sand',5),rep('maize',5)), points = c(rep(c(0,4,9,14,18),3),c(0,3,6,9,12)), soiltype = c(rep('rivierklei', 5), rep('veen',5), rep('dekzand',10)), B_LU_BRP = c(rep(234,15),rep(259,5))) dt.gras <- data.table(landuse = c(rep('grass on clay',5),rep('grass on sand',5),rep('grass on peat',5)), points = c(rep(c(0,3,5,8,11),2),c(0,4,9,14,17)), soiltype = c(rep('rivierklei',5),rep('dekzand',5),rep('veen',5))) # calculate indicator dt.gras <- dt.gras[,score := ind_management(points, rep(265,15),soiltype)] dt.arabl <- dt.arabl[,score := ind_management(points, B_LU_BRP,soiltype)] # make plots p.ar <- ggplot(data = dt.arabl, aes(points,score, col = landuse, group = landuse, shape = landuse)) + geom_point(size = 3)+ geom_line() + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Indicator value') + xlab('Soil management points')+theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10), legend.position = c(0.75,0.22)) p.gr <- ggplot(data = dt.gras, aes(points,score, col = landuse, group = landuse, shape = landuse)) + geom_point(size = 3)+ geom_line() + theme_bw() + scale_fill_viridis_d()+ scale_color_viridis_d()+ ylab('Indicator value') + xlab('Soil management points')+theme(plot.title = element_text(size=10), legend.text = element_text(size=10), axis.text = element_text(size=10), legend.position = c(0.75,0.2)) p.ar+p.gr ## ----results = FALSE, eval = FALSE-------------------------------------------- # # load weights.obic (set indicator to zero when not applicable) # w <- as.data.table(OBIC::weight.obic) # # # Add years per field # dt[,year := .I, by = ID] # # # Select all indicators used for scoring # cols <- colnames(dt)[grepl('I_C|I_B|I_P|I_E|I_M|year|crop_cat|SOILT',colnames(dt))] # # # Melt dt and assign main categories for OBI # dt.melt <- melt(dt[,mget(cols)], # id.vars = c('B_SOILTYPE_AGR','crop_category','year'), # variable.name = 'indicator') # # # add categories relevant for aggregating # # C = chemical, P = physics, B = biological, BCS = visual soil assessment # # indicators not used for integrating: IBCS and IM # dt.melt[,cat := tstrsplit(indicator,'_',keep = 2)] # dt.melt[grepl('_BCS$',indicator) & indicator != 'I_BCS', cat := 'IBCS'] # dt.melt[grepl('^I_M_',indicator), cat := 'IM'] # # # Determine number of indicators per category # dt.melt.ncat <- dt.melt[year==1 & !cat %in% c('IBCS','IM')][,list(ncat = .N),by='cat'] # # # add weighing factor to indicator values # dt.melt <- merge(dt.melt,w[,list(crop_category,indicator,weight_nonpeat,weight_peat)], # by = c('crop_category','indicator'), all.x = TRUE) # # # calculate correction factor for indicator values (low values have more impact than high values, a factor 5) # dt.melt[,cf := cf_ind_importance(value)] # # # calculate weighted value for crop category # dt.melt[,value.w := value] # dt.melt[grepl('veen',B_SOILTYPE_AGR) & weight_peat < 0,value.w := -999] # dt.melt[!grepl('veen',B_SOILTYPE_AGR) & weight_nonpeat < 0,value.w := -999] ## ----results=FALSE, eval=FALSE------------------------------------------------ # # subset dt.melt for relevant columns only # out.score <- dt.melt[,list(cat, year, cf, value = value.w)] # # # remove indicator categories that are not used for scoring # out.score <- out.score[!cat %in% c('IBCS','IM','BCS')] # # # calculate weighted average per indicator category # out.score <- out.score[,list(value = sum(cf * pmax(0,value) / sum(cf[value >= 0]))), by = list(cat,year)] # # # for case that a cat has one indicator or one year and has NA # out.score[is.na(value), value := -999] ## ----aggreagate indicators per category, results= FALSE, eval=FALSE----------- # # calculate correction factor per year; recent years are more important # out.score[,cf := log(12 - pmin(10,year))] # # # calculate weighted average per indicator category per year # out.score <- out.score[,list(value = sum(cf * pmax(0,value)/ sum(cf[value >= 0]))), by = cat] ## ----holistic obi score, results = FALSE, eval = FALSE------------------------ # # merge out with number per category # out.score <- merge(out.score,dt.melt.ncat, by='cat') # # # calculate weighing factor depending on number of indicators # out.score[,cf := log(ncat + 1)] # # # calculated final obi score # out.score <- rbind(out.score[,list(cat,value)], # out.score[,list(cat = "T",value = sum(value * cf / sum(cf)))]) ## ----results = FALSE, eval = FALSE-------------------------------------------- # # For example # OBIC::obic_field_dt(binnenveld[ID == 1], output = 'scores') # # OBIC::obic_field_dt(binnenveld[ID == 1], output = 'obic_score') ## ----eval = TRUE, echo=FALSE-------------------------------------------------- OBIC::obic_field_dt(binnenveld[ID == 1], output = 'scores') OBIC::obic_field_dt(binnenveld[ID == 1], output = 'obic_score') ## ----echo=FALSE, fig.height= 6, fig.width=6.8--------------------------------- # Load data dt <- binnenveld dts <- obic_field_dt(dt[ID== unique(dt$ID)[1]]) dts <- dts[,field := unique(dt$ID)[1]] for(i in unique(dt$ID)[2:4]){ dts2 <- obic_field_dt(dt[ID==i]) dts2 <- dts2[,field := i] dts <- rbindlist(list(dts, dts2)) } dts <- dts[, field := as.factor(field)] # g1 <- ggradar(dts[,.(field,S_T_OBI_A, S_C_OBI_A, S_P_OBI_A, S_B_OBI_A, S_M_OBI_A)], values.radar = c(0,0.5,1), # axis.label.size = 4, legend.position = "bottom", legend.title = 'field', # legend.text.size = 12, # grid.label.size = 4, # plot.extent.x.sf = 1.2, # axis.label.offset = 1.2) # pre prepare for lollipop chart cols <- colnames(dts)[grepl('^S_T_|^S_C_|^S_P_|^S_B_|^S_E_|^S_M_|^field', colnames(dts))] dts2 = dts[,mget(cols)] bar <- melt(dts2, id.vars = 'field') bar[grepl('^S_T_', variable), cat := "Total"] bar[grepl('^S_C_', variable), cat := "Chemical"] bar[grepl('^S_B_', variable), cat := "Biological"] bar[grepl('^S_P_', variable), cat := "Physical"] bar[grepl('^S_M_', variable), cat := "Management"] bar[grepl('^S_E_', variable), cat := "Environmental"] # change field for plotting bar <- bar[,field:= paste0('Field ', field)] # make lollipop chart (cleveland dotchart) without ggpubr g2p <- ggplot(bar, aes(x= cat, y= value, group = cat, color = cat)) + geom_col(fill = 'grey', color = NA, width = 0.05) + geom_point(size = 2) + ylab('OBI-score') + xlab('') + theme_bw(12) + theme(legend.position = 'none', panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = 0.5), panel.grid.minor.x = element_blank()) + coord_flip() + facet_wrap(~field) + scale_fill_viridis_d()+ scale_color_viridis_d() g2p ## ----include=FALSE------------------------------------------------------------ knitr::write_bib(c(.packages()), "packages.bib") knitr::write_bib(file = 'packages.bib')