## ----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 == '-', B_GWL_CLASS := 'I']
# 
#   # 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 == '-', B_GWL_CLASS := 'II']
  
  # 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  for an arable field
  dt.test2 <- dt.test[ID==11][1]
  dt.test2 <- dt.test2[rep(1,4)]
  dt.test2[,B_GWL_CLASS := c('III','IV','V','VI')]
  
  # 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('III','IV', 'V', 'VI', 'VII')]

  # 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')

