Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add evaluation for nitrogen #2

Merged
merged 10 commits into from
Jul 12, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,13 +1,18 @@
# Generated by roxygen2: do not edit by hand

export(calc_crumbleability)
export(calc_nlv)
export(calc_ph_delta)
export(calc_phosphate_availability)
export(calc_sealing)
export(eval_crumbleability)
export(eval_nitrogen)
export(eval_ph)
export(eval_phosphate_availability)
export(eval_resistance)
export(eval_sealing)
export(eval_structure)
export(evaluate_logistic)
export(evaluate_parabolic)
import(data.table)
importFrom(stats,approxfun)
37 changes: 34 additions & 3 deletions R/evaluate.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,21 @@
#' @param b (numeric) The growth rate
#' @param x0 (numeric) The offset of the x-axis
#' @param v (numeric) Affects the growth rate near the maximum
#' @param increasing (boolean) Should the evaluation increase (\code{TRUE}) with x or decrease (\code{FALSE})?
#'
#' @references \url{https://en.wikipedia.org/wiki/Generalised_logistic_function}
#'
#' @export
evaluate_logistic <- function(x, b, x0, v) {
evaluate_logistic <- function(x, b, x0, v, increasing = TRUE) {

# Settings
A <- 0 # Lower asympote
K <- 1 # Upper asympote
if (increasing) {
A <- 0 # Lower asympote
K <- 1 # Upper asympote
} else {
A <- 1 # Lower asympote
K <- 0 # Upper asympote
}
C <- 1

# General logistic function
Expand All @@ -23,3 +29,28 @@ evaluate_logistic <- function(x, b, x0, v) {
return(y)

}

#' Evaluate using parabolic function with
#'
#' This function evaluates the calculated values from an indicator using a parabolic function. After the optimum is reached the it stays at its plateau.
#'
#' @param x (numeric) The values of a calc function to be converted to an evaluation
#' @param x.top (numeric) The value at which x reaches the plateau
#'
#' @export
evaluate_parabolic <- function(x, x.top) {

# Setting
a <- 1 / x.top^2
b <- x.top

# Calcute the values
y <- 1 - a * (x - b) ^2

# Set plateaus
y <- ifelse(x >= x.top, 1, y)
y <- ifelse(y < 0, 0, y)

return(y)

}
107 changes: 107 additions & 0 deletions R/nitrogen.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
#' Calculate the NLV
#'
#' This function calculates the NLV (nitrogen producing capacity) for the soil
#'
#' @param n.org (numeric) The organic nitrogen content of the soil in g N / kg
#' @param c.org (numeric) The organic carbon content of the soil in kg C / ha
#' @param crop (numeric) The crop code (gewascode) from the BRP
#' @param soiltype (character) The type of soil
#' @param bulk_density (numeric) The bulk density of the soil in kg / m3
#' @param cn_ratio (numeric)
#' @param grass.age (numeric) The age of the grass if present
#'
#' @import data.table
#'
#' @export
calc_nlv <- function(n.org, c.org, crop, soiltype, bulk_density, cn_ratio, grass.age) {

# Load in the datasets
crops.obic <- as.data.table(OBIC::crops.obic)
setkey(crops.obic, crop_code)
soils.obic <- as.data.table(OBIC::soils.obic)
setkey(soils.obic, soiltype)

# Check input
arg.length <- max(length(n.org), length(c.org), length(crop), length(soiltype), length(bulk_density), length(grass.age))
checkmate::assert_numeric(n.org, lower = 0, upper = 30, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(c.org, lower = 0, upper = 1000, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(crop, any.missing = FALSE, min.len = 1, len = arg.length)
checkmate::assert_subset(crop, choices = unique(crops.obic$crop_code), empty.ok = FALSE)
checkmate::assert_character(soiltype, any.missing = FALSE, min.len = 1, len = arg.length)
checkmate::assert_subset(soiltype, choices = unique(soils.obic$soiltype), empty.ok = FALSE)
checkmate::assert_numeric(bulk_density, lower = 0, upper = 1500, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cn_ratio, lower = 0, upper = 20, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(grass.age, lower = 0, upper = 99, len = arg.length)

# Settings
param.a <- 20 # Age of organic matter
param.b <- 2^((14.1 - 9)/ 9) # Temperature correction
param.cn.micro <- 10 # CN ratio of micro organisms
param.t <- 5 / 12 # 5 months a year
param.diss.micro <- 2 # Dissimilation : assimilation ratio of micro organisms

# Collect data in a table
a = c.ass = c.diss = crop_code = id = soiltype.n = crop_n = NULL
dt <- data.table(
id = 1:arg.length,
n.org = n.org,
c.org = c.org,
crop = crop,
soiltype = soiltype,
bulk_density = bulk_density,
cn_ratio = cn_ratio,
grass.age = grass.age,
value.nlv = NA_real_
)
dt <- merge(dt, crops.obic[, list(crop_code, crop_n)], by.x = "crop", by.y = "crop_code")
dt <- merge(dt, soils.obic[, list(soiltype, soiltype.n)], by = "soiltype")

# Calculate NLV for grass
dt.grass <- dt[crop_n == "gras"]
dt.grass[soiltype.n == "zand" & grass.age < 4, a := 30.79]
dt.grass[soiltype.n == "zand" & grass.age >= 4 & grass.age < 7, a := 28.36]
dt.grass[soiltype.n == "zand" & grass.age >= 7 & grass.age < 10, a := 27.78]
dt.grass[soiltype.n == "zand" & grass.age >= 10, a := 26.57]
dt.grass[soiltype.n == "klei" & grass.age < 4, a := 34.25]
dt.grass[soiltype.n == "klei" & grass.age >= 4 & grass.age < 7, a := 31.54]
dt.grass[soiltype.n == "klei" & grass.age >= 7 & grass.age < 10, a := 30.90]
dt.grass[soiltype.n == "klei" & grass.age >= 10, a := 29.56]

dt.grass[soiltype.n == "zand", value.nlv := 78 + a * n.org ^ 1.0046]
dt.grass[soiltype.n == "klei", value.nlv := 31.7 + a * n.org ^ 1.0046]
dt.grass[soiltype.n == "zand" & value.nlv > 200, value.nlv := 200]
SvenVw marked this conversation as resolved.
Show resolved Hide resolved
dt.grass[soiltype.n == "klei" & value.nlv > 250, value.nlv := 250]
dt.grass[soiltype.n == "veen", value.nlv := 250]

# Calculate the NLV for arable land
dt.arable <- dt[crop_n == "akkerbouw"]
dt.arable[, c.diss := c.org * exp(4.7 * ((param.a + param.b * param.t)^-0.6)) - param.a^-0.6]
dt.arable[, c.ass := c.diss / param.diss.micro]
dt.arable[, value.nlv := ((c.diss + c.ass) / cn_ratio) - (c.ass / param.cn.micro)]

# Combine both tables and extract values
dt <- rbindlist(list(dt.grass, dt.arable), fill = TRUE)
setorder(dt, id)
value.nlv <- dt[, value.nlv]

return(value.nlv)
}

#' Evaluate nitrogen
#'
#' This function evaluates the nitrogen content of the soil by using the NLV calculated by \code{\link{calc_nlv}}
#'
#' @param value.nlv (numeric) The value of NLV calculated by \code{\link{calc_nlv}}
#'
#' @export
eval_nitrogen <- function(value.nlv) {

# Check inputs
checkmate::assert_numeric(value.nlv, lower = -30, upper = 250, any.missing = FALSE)

# Evaluate the nitrogen
eval.nitrogen <- OBIC::evaluate_parabolic(value.nlv, x.top = 120)

# return output
return(eval.nitrogen)
}
131 changes: 131 additions & 0 deletions R/ph.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
#' Calculate the difference between pH and optimum
#'
#' This functions calculates the difference between the measured pH and the optimal pH according to the Bemestingsadvies
#'
#' @param ph.cacl2 (numeric) The pH-CaCl2 of the soil
#' @param soiltype (character) The type of soil
#' @param lutum (numeric) The percentage lutum present in the soil
#' @param om (numeric) The organic matter content of soil in percentage
#' @param cp.starch (numeric) The fraction of starch potatoes in the crop plan
#' @param cp.potato (numeric) The fraction of potatoes (excluding starch potatoes) in the crop plan
#' @param cp.sugarbeet (numeric) The fraction of sugar beets in the crop plan
#' @param cp.grass (numeric) The fracgtion of grass in the crop plan
#' @param cp.mais (numeric) The fraction of mais in the crop plan
#' @param cp.other (numeric) The fraction of other crops in the crop plan
#'
#' @references \href{https://www.handboekbodemenbemesting.nl/nl/handboekbodemenbemesting/Handeling/pH-en-bekalking/Advisering-pH-en-bekalking.htm}{Handboek Bodem en Bemesting tabel 5.1, 5.2 en 5.3}
#'
#' @import data.table
#'
#' @export
calc_ph_delta <- function(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, cp.sugarbeet, cp.grass, cp.mais, cp.other) {

# Load in the datasets
soils.obic <- as.data.table(OBIC::soils.obic)
setkey(soils.obic, soiltype)
dt.ph.delta <- as.data.table(OBIC::tbl.ph.delta)

# Check inputs
arg.length <- max(length(ph.cacl2), length(soiltype), length(om), length(lutum), length(cp.starch), length(cp.potato), length(cp.sugarbeet), length(cp.grass), length(cp.mais), length(cp.other))
checkmate::assert_character(soiltype, any.missing = FALSE, len = arg.length)
checkmate::assert_subset(soiltype, choices = unique(soils.obic$soiltype))
checkmate::assert_numeric(om, lower = 0, upper = 100, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(lutum, lower = 0, upper = 100, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.starch, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.potato, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.sugarbeet, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.grass, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.mais, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(cp.other, lower = 0, upper = 1, any.missing = FALSE, len = arg.length)
cp.total <- cp.starch + cp.potato + cp.sugarbeet + cp.grass + cp.mais + cp.other
if (any(cp.total != 1)) {
stop(paste0("The sum of the fraction of cp is not 1, but ", min(cp.total)))
}

# Collect information in table
soiltype.ph = lutum.low = lutum.high = om.low = om.high = potato.low = potato.high = sugarbeet.low = sugarbeet.high = ph.optimum = id = NULL
dt <- data.table(
id = 1:arg.length,
ph.cacl2 = ph.cacl2,
soiltype = soiltype,
om = om,
lutum = lutum,
cp.starch = cp.starch,
cp.potato = cp.potato,
cp.sugarbeet = cp.sugarbeet,
cp.grass = cp.grass,
cp.mais = cp.mais,
cp.other = cp.other,
table = NA_character_
)

# Join soil type used for this function
dt <- merge(dt, soils.obic, by = "soiltype")

# Define which table to be used
dt[soiltype.ph == 1, table := "5.1"]
dt[soiltype.ph == 2, table := "5.3"]
dt[cp.starch > 0.1, table := "5.2"]
dt[cp.grass + cp.mais >= 0.5, table := "5.3"]

dt[, cp.potato := cp.starch + cp.potato]

# Join conditionally the tables with optimum pH to data
dt.53 <- dt[table == "5.3"]
dt.53 <- dt.ph.delta[dt.53, on=list(table == table, lutum.low <= lutum, lutum.high > lutum, om.low <= om, om.high > om)]
dt.512 <- dt[table %in% c("5.1", "5.2")]
dt.512 <- dt.ph.delta[dt.512, on=list(table == table, potato.low <= cp.potato, potato.high > cp.potato, sugarbeet.low <= cp.sugarbeet, sugarbeet.high > cp.sugarbeet, om.low <= om, om.high > om)]
dt <- rbindlist(list(dt.53, dt.512), fill = TRUE)

# Calculate the difference between the measured pH and the optimum pH
dt[, ph.delta := ph.optimum - ph.cacl2]
dt[ph.delta < 0, ph.delta := 0]

# Extract the ph.delta
setorder(dt, id)
ph.delta <- dt[, ph.delta]

return(ph.delta)

}

#' Evaluate the difference between de optimum pH and ph of the soil
#'
#' This function evaluates the pH of the soil by the difference with the optimum pH. This is calculated in \code{\link{calc_ph_delta}}.
#'
#' @param value.ph.delta (numeric) The pH difference with the optimal pH.
#'
#' @export
eval_ph <- function(value.ph.delta) {

# Check inputs
checkmate::assert_numeric(value.ph.delta, lower = 0, upper = 5, any.missing = FALSE)

# Evaluate the pH
eval.ph <- OBIC::evaluate_logistic(x = value.ph.delta, b = 9, x0 = 0.3, v = 0.4, increasing = TRUE)

return(eval.ph)

}

#' Table with optimal pH for different crop plans
#'
#' This table contains the optimal pH for different crop plans and soil types
#'
#' @format A data.frame with 84 rows and 10 columns:
#' \describe{
#' \item{table}{The original table from Hanboek Bodem en Bemesting}
#' \item{lutum.low}{Lower value for lutum}
#' \item{lutum.high}{Upper value for lutum}
#' \item{om.low}{Lower value for organic matter}
#' \item{om.high}{Upper value for organic matter}
#' \item{potato.low}{Lower value for fraction potatoes in crop plan}
#' \item{potato.high}{Upper value for fraction potatoes in crop plan}
#' \item{sugarbeet.low}{Lower value for fraction potatoes in crop plan}
#' \item{sugarbeet.high}{Upper value for fraction potatoes in crop plan}
#' \item{ph.optimum}{The optimal pH for this range}
#' }
#'
#' #' @references \href{https://www.handboekbodemenbemesting.nl/nl/handboekbodemenbemesting/Handeling/pH-en-bekalking/Advisering-pH-en-bekalking.htm}{Handboek Bodem en Bemesting tabel 5.1, 5.2 en 5.3}
#'
"tbl.ph.delta"
15 changes: 14 additions & 1 deletion R/crops_obic.R → R/tables.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,5 +9,18 @@
#' \item{crop_crumbleability}{The category for this crop at crumbleablity}
#' \item{crop_phosphate}{The category for this crop at phosphate availability}
#' \item{crop_sealing}{The category for this crop at soil sealing}
#' \item{crop_n}{The category for this crop at nitrogen}
#' }
"crops.obic"
"crops.obic"

#' Linking table between soils and different functions in OBIC
#'
#' This table helps to link the different crops in the OBIC functions with the crops selected by the user
#'
#' @format A data.frame with 7 rows and 2 columns:
#' \describe{
#' \item{soiltype}{The name of the soil type}
#' \item{soiltype.ph}{The category for this soil at pH}
#' \item{soiltype.n}{The category for this soil at nitrogen}
#' }
"soils.obic"
4 changes: 4 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,8 @@
# Open Bodem Index Calculator (OBIC)

<!-- badges: start -->
[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental)
<!-- badges: end -->

This package contains the functions to calculate the indicators and evaluations used in Open Bodem Index (OBI).
The Open Bodem Index is a tool to evaluate the soil of agricultural fields based on various criteria.
Binary file modified data/crops_obic.RData
Binary file not shown.
Binary file added data/soils_obic.RData
Binary file not shown.
Binary file added data/tbl_ph_delta.RData
Binary file not shown.
26 changes: 26 additions & 0 deletions man/calc_nlv.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading