From 1532a4625779bac013a87c962daf35d925cea502 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Tue, 9 Jul 2019 11:18:03 +0200 Subject: [PATCH 01/10] First checks for function to calculate ph_delta --- R/ph.R | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 R/ph.R diff --git a/R/ph.R b/R/ph.R new file mode 100644 index 00000000..9dac745f --- /dev/null +++ b/R/ph.R @@ -0,0 +1,40 @@ +#' 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 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.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, cp.starch, cp.potato, cp.sugarbeet, cp.mais, cp.other) { + + # Check inputs + arg.length <- max(length(ph.cacl2), length(cp.starch), length(cp.potato), length(cp.sugarbeet), length(cp.mais), length(cp.other)) + checkmate::assert_numeric(ph.cacl2, lower = 0, upper = 14, 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.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.mais + cp.other + if (cp.total != 1) { + stop(paste0("The sum of the fraction of cp is not 1, but ", cp.total)) + } + + # Use tables for starch potatoes if present in crop plan + if (cp.starch > 0.1) { + pr.starch <- TRUE + } else { + pr.starch <- FALSE + } + + +} \ No newline at end of file From f9914e40479a3b1b9922ed503e39f5e08d412ed0 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Tue, 9 Jul 2019 11:37:30 +0200 Subject: [PATCH 02/10] Adds linking table for soiltypes --- R/{crops_obic.R => tables.R} | 13 ++++++++++++- data/soils_obic.RData | Bin 0 -> 219 bytes man/crops.obic.Rd | 2 +- man/soils.obic.Rd | 18 ++++++++++++++++++ 4 files changed, 31 insertions(+), 2 deletions(-) rename R/{crops_obic.R => tables.R} (60%) create mode 100644 data/soils_obic.RData create mode 100644 man/soils.obic.Rd diff --git a/R/crops_obic.R b/R/tables.R similarity index 60% rename from R/crops_obic.R rename to R/tables.R index 92da4063..30240fa1 100644 --- a/R/crops_obic.R +++ b/R/tables.R @@ -10,4 +10,15 @@ #' \item{crop_phosphate}{The category for this crop at phosphate availability} #' \item{crop_sealing}{The category for this crop at soil sealing} #' } -"crops.obic" \ No newline at end of file +"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} +#' } +"soils.obic" \ No newline at end of file diff --git a/data/soils_obic.RData b/data/soils_obic.RData new file mode 100644 index 0000000000000000000000000000000000000000..64d35dc69c341436d7778a0bb026ff701d85cd6f GIT binary patch literal 219 zcmV<103`n(iwFP!0000027ORV4uUWc9a~VrPfWZA3ooDxPvD9jDKaXtQbGwvRvy`{ z8*!>_Vp!N@X5M?9*U5CZjaOkD0ssTZ!%}H&@${mw%8!`*+)ut*(uyK;mIkqK_(%6!(V!$SH?ubIJ+7`F|EufkZQ~c> zg4!;?tT;9+*0&NZ9~mB%r=B;c)HQPcnepE#bzS%D57|qh-4JfYVV%k(J6Ohvxev}n VEuuq34- Date: Tue, 9 Jul 2019 16:19:32 +0200 Subject: [PATCH 03/10] Adds function to calculate delta ph --- NAMESPACE | 1 + R/ph.R | 98 +++++++++++++++++++++++++++++++++------ data/tbl_ph_delta.RData | Bin 0 -> 709 bytes man/calc_ph_delta.Rd | 36 ++++++++++++++ man/tbl.ph.delta.Rd | 28 +++++++++++ tests/testthat/test-ph.R | 34 ++++++++++++++ 6 files changed, 184 insertions(+), 13 deletions(-) create mode 100644 data/tbl_ph_delta.RData create mode 100644 man/calc_ph_delta.Rd create mode 100644 man/tbl.ph.delta.Rd create mode 100644 tests/testthat/test-ph.R diff --git a/NAMESPACE b/NAMESPACE index 64a488fa..66c93848 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(calc_crumbleability) +export(calc_ph_delta) export(calc_phosphate_availability) export(calc_sealing) export(eval_crumbleability) diff --git a/R/ph.R b/R/ph.R index 9dac745f..f206996c 100644 --- a/R/ph.R +++ b/R/ph.R @@ -3,9 +3,13 @@ #' 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 #' @@ -14,27 +18,95 @@ #' @import data.table #' #' @export -calc_ph_delta <- function(ph.cacl2, cp.starch, cp.potato, cp.sugarbeet, cp.mais, cp.other) { +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(cp.starch), length(cp.potato), length(cp.sugarbeet), length(cp.mais), length(cp.other)) - checkmate::assert_numeric(ph.cacl2, lower = 0, upper = 14, any.missing = FALSE, len = arg.length) + 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.mais + cp.other - if (cp.total != 1) { - stop(paste0("The sum of the fraction of cp is not 1, but ", cp.total)) + 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))) } - # Use tables for starch potatoes if present in crop plan - if (cp.starch > 0.1) { - pr.starch <- TRUE - } else { - pr.starch <- FALSE - } + # Collect information in table + soil.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[soil.ph == 1, table := "5.1"] + dt[soil.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 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) + +} -} \ No newline at end of file +#' 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" diff --git a/data/tbl_ph_delta.RData b/data/tbl_ph_delta.RData new file mode 100644 index 0000000000000000000000000000000000000000..ebc6d205e8db540ac4bc342dbe09a5640dd8ecf8 GIT binary patch literal 709 zcmV;$0y_O4iwFP!000002JMcNuX6}N*me;bzo?hlSpA;A&0~#U@8(4 zx{#!)bm5}pqAp17z|IUK8y|s{xlBk*e26YcOax!)Y8VbyoqO{Z4(Lyst15W2{ISzK>mlYQXq z+N!PE>P%Me{lnF=txj4cTB2yrHeB2-!`wE^vSB_n%twa#_}D0qPNVz62g{2gVV|(; zn4w+AXY-kBlo!2D{&$jco^U0vfveyh@GiImI^ca!o%fhW*1z?CF8-yztT|&LAt_PonZg%a8A(@?OiYyMv-7 z)1R19`PLW|9l8E#OXXzy$LF`tceG^2$30)N`-J5sVC59^T}N&o;rgfCM_B*dI^}41 zAK=DMcAv1k1f-Ou@(JsCN~VpUdVJhI!7+XQcAv1Ed8VD`?C*VArX3%*&P|Vh!t;c? z5uZQY+wnQeEz+R4YtXMjz7^-?u0k*W|0tl}N_@_7Z^CXFc8P!Ji?CaVJart%(*FkR zk~}7j1G@s`71+r_?1 zQ%R*RAM`6iX7fP2B(9QQ(hfPXQ=Yf(y-WLbg`B|H^CCyA*#go0rPatNL{aIlmY7T8FZ7&7M>rPJuju*>@)v?km@EJQk!^ER literal 0 HcmV?d00001 diff --git a/man/calc_ph_delta.Rd b/man/calc_ph_delta.Rd new file mode 100644 index 00000000..45df455b --- /dev/null +++ b/man/calc_ph_delta.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ph.R +\name{calc_ph_delta} +\alias{calc_ph_delta} +\title{Calculate the difference between pH and optimum} +\usage{ +calc_ph_delta(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, + cp.sugarbeet, cp.grass, cp.mais, cp.other) +} +\arguments{ +\item{ph.cacl2}{(numeric) The pH-CaCl2 of the soil} + +\item{soiltype}{(character) The type of soil} + +\item{om}{(numeric) The organic matter content of soil in percentage} + +\item{lutum}{(numeric) The percentage lutum present in the soil} + +\item{cp.starch}{(numeric) The fraction of starch potatoes in the crop plan} + +\item{cp.potato}{(numeric) The fraction of potatoes (excluding starch potatoes) in the crop plan} + +\item{cp.sugarbeet}{(numeric) The fraction of sugar beets in the crop plan} + +\item{cp.grass}{(numeric) The fracgtion of grass in the crop plan} + +\item{cp.mais}{(numeric) The fraction of mais in the crop plan} + +\item{cp.other}{(numeric) The fraction of other crops in the crop plan} +} +\description{ +This functions calculates the difference between the measured pH and the optimal pH according to the Bemestingsadvies +} +\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} +} diff --git a/man/tbl.ph.delta.Rd b/man/tbl.ph.delta.Rd new file mode 100644 index 00000000..8a27eda5 --- /dev/null +++ b/man/tbl.ph.delta.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ph.R +\docType{data} +\name{tbl.ph.delta} +\alias{tbl.ph.delta} +\title{Table with optimal pH for different crop plans} +\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}} +\usage{ +tbl.ph.delta +} +\description{ +This table contains the optimal pH for different crop plans and soil types +} +\keyword{datasets} diff --git a/tests/testthat/test-ph.R b/tests/testthat/test-ph.R new file mode 100644 index 00000000..58d17c2e --- /dev/null +++ b/tests/testthat/test-ph.R @@ -0,0 +1,34 @@ +test_that("calc_ph_delta works", { + expect_equal( + calc_ph_delta( + ph.cacl2 = 6, + soiltype = "klei", + lutum = 20, + om = 5, + cp.starch = 0, + cp.potato = 0.3, + cp.grass = 0, + cp.mais = 0.2, + cp.sugarbeet = 0.2, + cp.other = 0.3 + ), + expected = 0.3, + tolerance = 0.001 + ) + expect_equal( + calc_ph_delta( + ph.cacl2 = c(6, 5.8, 6.2, 5.6), + soiltype = c("klei", "veen", "veen", "zavel"), + lutum = c(20, 5, 8, 12), + om = c(5, 20, 23, 8), + cp.starch = c(0, 0.2, 0, 0.4), + cp.potato = c(0.3, 0.15, 0, 0), + cp.grass = c(0, 0.45, 0.7, 0), + cp.mais = c(0.1, 0.2, 0.15, 0.1), + cp.sugarbeet = c(0.2, 0, 0.15, 0), + cp.other = c(0.4, 0, 0, 0.5) + ), + expected = c(0.3, 0, 0, 0), + tolerance = 0.001 + ) +}) \ No newline at end of file From 55a3e1bdfb83b495037f5a28eb8c59520ff8f83e Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Tue, 9 Jul 2019 16:41:14 +0200 Subject: [PATCH 04/10] Adds evaluation for pH --- NAMESPACE | 1 + R/evaluate.R | 12 +++++++++--- R/ph.R | 21 ++++++++++++++++++++- man/eval_ph.Rd | 14 ++++++++++++++ man/evaluate_logistic.Rd | 4 +++- tests/testthat/test-ph.R | 10 ++++++++++ 6 files changed, 57 insertions(+), 5 deletions(-) create mode 100644 man/eval_ph.Rd diff --git a/NAMESPACE b/NAMESPACE index 66c93848..55c639ff 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ export(calc_ph_delta) export(calc_phosphate_availability) export(calc_sealing) export(eval_crumbleability) +export(eval_ph) export(eval_phosphate_availability) export(eval_resistance) export(eval_sealing) diff --git a/R/evaluate.R b/R/evaluate.R index 0ccaace2..e3efcbe4 100644 --- a/R/evaluate.R +++ b/R/evaluate.R @@ -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 diff --git a/R/ph.R b/R/ph.R index f206996c..532ff918 100644 --- a/R/ph.R +++ b/R/ph.R @@ -70,7 +70,7 @@ calc_ph_delta <- function(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, c dt[, cp.potato := cp.starch + cp.potato] - # Join the tables with optimum pH to data + # 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")] @@ -89,6 +89,25 @@ calc_ph_delta <- function(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, c } +#' 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 diff --git a/man/eval_ph.Rd b/man/eval_ph.Rd new file mode 100644 index 00000000..80438bae --- /dev/null +++ b/man/eval_ph.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ph.R +\name{eval_ph} +\alias{eval_ph} +\title{Evaluate the difference between de optimum pH and ph of the soil} +\usage{ +eval_ph(value.ph.delta) +} +\arguments{ +\item{value.ph.delta}{(numeric) The pH difference with the optimal pH.} +} +\description{ +This function evaluates the pH of the soil by the difference with the optimum pH. This is calculated in \code{\link{calc_ph_delta}}. +} diff --git a/man/evaluate_logistic.Rd b/man/evaluate_logistic.Rd index 90c1ba14..1285b7b7 100644 --- a/man/evaluate_logistic.Rd +++ b/man/evaluate_logistic.Rd @@ -4,7 +4,7 @@ \alias{evaluate_logistic} \title{Evaluate using the general logistice function} \usage{ -evaluate_logistic(x, b, x0, v) +evaluate_logistic(x, b, x0, v, increasing = TRUE) } \arguments{ \item{x}{(numeric) The values of a calc function to be converted to an evaluation} @@ -14,6 +14,8 @@ evaluate_logistic(x, b, x0, v) \item{x0}{(numeric) The offset of the x-axis} \item{v}{(numeric) Affects the growth rate near the maximum} + +\item{increasing}{(boolean) Should the evaluation increase (\code{TRUE}) with x or decrease (\code{FALSE})?} } \description{ This function evaluates the calculated values from an indicator using a general logistic function diff --git a/tests/testthat/test-ph.R b/tests/testthat/test-ph.R index 58d17c2e..64e3d1e1 100644 --- a/tests/testthat/test-ph.R +++ b/tests/testthat/test-ph.R @@ -31,4 +31,14 @@ test_that("calc_ph_delta works", { expected = c(0.3, 0, 0, 0), tolerance = 0.001 ) +}) + +test_that("eval_ works", { + expect_equal( + eval_ph( + value.ph.delta = seq(from = 0, to = 1.5, by = 0.2) + ), + expected = c(0.0009951581, 0.0449194368, 0.4261837504, 0.8499235247, 0.9727585651, 0.9954239501, 0.9992415551, 0.9998745743), + tolerance = 0.001 + ) }) \ No newline at end of file From f6a1aa176bd77c0f22e9fd310ee5f11a260ae235 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Wed, 10 Jul 2019 10:17:41 +0200 Subject: [PATCH 05/10] Add first test for calc_nlv --- R/tables.R | 2 ++ data/crops_obic.RData | Bin 5262 -> 5305 bytes data/soils_obic.RData | Bin 219 -> 233 bytes tests/testthat/test-nitrogen.R | 26 ++++++++++++++++++++++++++ 4 files changed, 28 insertions(+) create mode 100644 tests/testthat/test-nitrogen.R diff --git a/R/tables.R b/R/tables.R index 30240fa1..cdb60864 100644 --- a/R/tables.R +++ b/R/tables.R @@ -9,6 +9,7 @@ #' \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" @@ -20,5 +21,6 @@ #' \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" \ No newline at end of file diff --git a/data/crops_obic.RData b/data/crops_obic.RData index 39f4f020c6cdcb00745d83adc38c335430e3607c..911b84e4828aa31261e74ec2eca58900c1274693 100644 GIT binary patch literal 5305 zcmYL{byO4pv-VLW6lnqJUO<|aW??B&VCjYhr8}1HkWNYIhNY2|MiHc&C6<=%4q;iY z-}{~S-Z{^i-!spd|7Pf8a54T<*hf&OSQ-VL)r;1jouH3jYsnfy$yD*AL5sw z$D6eMfv|RuwB2XuTf}YJ9Ht5x%v?ofWoL6vjSR zmF*_WAb$Qj3)6#ZZc@YMsf}jt%Q;;@GvttI3F;BH0f|Y{&rP3i%Tx527{`=b3;t73 z$m-QF`esQeR5lP{^hsB()ibJ<+nCRw^vKooMTxH(l6=&7eI;h;3I!=Dl)9$=(`l z$#5s=;n9@FYMy;^B9wLOmTKu@mRagZF_az2mm`nG=fX#9~w&S-+|KETM1c9 z5t_qYF+3Msmb%s>X|44`W-&)6h}&(&Fii(c?Fc6S*U4mGvD4Lp>DKVnHr{-HqgjlV z#Uh-J%A!vFQEXiLB3zcLybPW;xSV8Puy|qsgRF7^l@gp(-*$yNLml^sC0~=Jh1q?T zAGFSN$}gNCUA7D$YSFm1*+?QnC<=ZX?i$3`hmuI!Q90yN7D2QM3{SD-p{}jVfF?@w zF};YX_E*BPTKId}Z6V3cE%E_IUt<@m3NCF^T6TaR=)P25@=IIoOp7s}F-J);H~i9_ z0b6^tW-n#Gv_V?fA$R$$vr z^rpZ$3d>YW69n!k1(q6al)sg^y}@9Z_-Yb>9QxmvdJSl8QR&V@R-A_R#H#L;nOj_*AfU z^1Ld+S?rjzJh9_`_wyrods*=}HbjMcy{nL6l03#$lLzvS${!Nr;mqT#OmbJ%y3v-dOyXmPXjj;9GpN zq!427Mg<33+nQdxL(lYI_6gKiawMm>Pes#B>Nlq(s@7eHnex}b;ry8_FS$iw{&pHh z&c)0c$EV@3sl%o6F7Pobcj8+YxPOU3VuQc{M?4bMkoGFtG`DP*%yBN{8k&B8tOusnXEN8l!CS=3RAIW7 zoLqJ0Bb_pjDE`%0gvcA&noz!>jVTS$8D$V)BffuKTLD zETAvV4Tkx=xp8b1!5=kJ|9xo+rJv}vb$ie}E|1xKUP?6URYh<;pqZZ2aS7nm91^3e zjWFGRrPFWXQBwDNy@8hr@(E@-OBVN0SqSH?%ZnnmodJNKgyg++7GRw0eA^|$9ZU+h zvEula&EG_tSjs1Ql^L5!cAR>!;WJhA0MSTS=}t0dW*ZRkBXqJ&?eT$?4lW6E+Zka= zX6LJZvf9-zh*s*Uw|+XE$5HgDl464Ppw$n1YjVB-h>L$@OX;StrufEZE4^lsjg?0dU~Q+9%9Fxh|9uV>ZtD&~#*uvfQB@NwI!;e* zolH)M^Zp&)*$7u6LQ3fwXl9Oy+wqDx_89#`8)?@DD9LW*MPxl4WmBifaV*0mAc+yq z9MRPlOGS#>;SEiOu5g^(QW`K~O?Qp}$Auk!Q*m^%?AQBbWDDE;EFn&yeo1o45iKUs zuVYc@o`J9}H7L$o;ab(AQcufQ)sU1M7V}0&f3vQ%QSOpal4y`u#2JsvAKkaQksn&P zt6%Emont@Y#(`|mShPFsyLqK8j0k!&@gqIM=C#-hibv9f(tEYp=n1)oXU1VIBY!Qp zrj%~{jOx8#u_dMIn@YLY9Cz*2xq#3g6D0bg@2%`JBvKbF3*v7o>gC$3a;LV9N!rXx zncgQD=<2pAQ5HD89VJVxsO?=i5JhUGSi4pUI%LBcY#RjNIaY>_Ec+?>rFst|L~UH6 zHzO`Q!cwd&`7FkUT#Dh*JqBcs4W23P3~x-zvK$yQAKX<(fJ$Lb3b&f?QjYFog@82@ z9nnqq?@eDM`eg2df0SrxtJdW@OjUdUsJhewbu=3kEO5>^Y#3Kev7E}X9W)_l9ZJ|N zAQ!W9!Jk~Me|5>Zfx|ECONB@qbG{7 z!dZJsIf_nOja%4#5VF5}Gk{#Pp&wbZYkgRg+s!+}W=#HX$Z|doL7t!(8Ld`N-O>OQ zp!x#+5?X6p-sGghv7*in^M!9C+qhXYP?)zYWo<8LnJ4b_R$Dl~6aV+Zt8+EiNFmkh zM9aOEUARsGQ{3@nWn;LrS+AZ-$=#=lg9AXsdJ=OKPoT;P&EVyMnaKWoY?I;hg~aMS z>Z2~^MPt-qZuSsw1_>bB$9#1>{_-Rdg11*ASwA|%S4%IL4`fCKYPD;Y13wk}vEcyv z0DpExG_AOR@xKPUrhmU=mMzx_XJsSV0?@&+LifSMXFksCdnJ_x)qc-uW5uaGcKc^%GQ9N)B;e1-V{3JAJt$+1>k)_$_qq60e=zY zdhvzC=%E4ETtlr|DLb1_m;zp*gn62$?!K+8e2{ke?2t|L3wPt{YmM*sjidlVJKmO9 zq2wbyZt-ZXFCwdm8s3)K&0q40pc2_Z!3+;!ckhX1Jo`q?dQO~cOu{$G-{QQcfwH?aMCGF5==<8M+0h`HPVJ4!ibLvu4oL>&Zb~`^3^nNaZ1b# zK_5KET!7vxkMpLjbOU`X#tn;d%2bE6G&Gq-BpQUS)-^-RdL3Ve ztqMc=a%M9P=#9LvKfd8@Jzv);yJE(`mV$-@FW}PCnY`(|IiJ#-hJRvXA2Iu$h?67g zx{mWz{$a7ewFX9`Z>3{L6U#IO`u7ESB}T7bBAtfUkTraD2%jLi0`t|whp_x ziY)o2MI(+C9v3e>PUvJ7lUZ;a87?HbVwcUjC@U|>NWI_UJ^fR$ULYp9quk9a6P?f} zm~qF`R%k??mQi?OtWEwkVjRKIf7y~KSpmPRY-phsVS75oT8+(t6B*T}T+WW6murO2 z!Cw3xVUp8#`^GDhw%4S&V!*8!$$s8!6fZELRTwfqBiQ;`Q2OyrN|F zD{S5y5s72h)NDvO(b5nkn=?n+CgW-b(mqQtftct{v@KjjTL!Fe ztO!22)z=CEnRtc>-PmpVpA3JLJFFf1<6+8Akjpw$UC2z<}mnPox9m1ZsPCJN1vRXLf2 zR}}(tubHko-3WaI)W22>D>zQGoY|;?YEDxqVPZL4G56*2NG0!mweZUPPs#yGjzf?-JhqE5Z_hye~A%jVANlWGPfXy-xkSTRns`oK64t{LYR{Ivoy9wIK@q~ zSMI>4Gs@Ib+t7iaN@+{Wa`Xr3;;*}3xLF&N3kLt(>lMa{r7VtI&sZiYCuDfGiW*J) z`w*a+g`A{ZN3wrY#tK|YB)jgGdZcE;)XevcEmU5`vMwE}WSoD(mE<(A&bDe1Uasm; z=;>`vADGU4Uf0$QASiNrr;KyuUchyZ34WD$ebe@FwQ!AWO2FOCN4=%VN21D? zJm-hI+2V1D$QwfOdnxsWR&${+lbtiEs1`*V*$c2_G0(R*n|hd~%+rowG#`8F(BM!e7&06I^p& zy9ZA!=6{!V`aCv!l_$5Wj6si4nB$5a89;;6LLe7Nck5er<0g7{CK*5RNE8V(`6^5B zIE{-rKnl*qIDfWA|B>d&-SZ#cvD(NtZGyYP0t|NHVF5_%)pz$on|Bs7CV|gzyH@4z zPtXdX!ZlA0z?YsI_zasJPtd~XXORJpC_$q08o8TtNdF2?aQ zNsLgGlPvl>13KKhn_7SARZpLqT)=V^^5qSkDKS{3m`SxS%SE zpo#uFSMaXIgZ!F4!)2o+PBsqDlY;M8c{s+HIzbZ~ZCGwWl0-hxNBOlAOsQvRVxl{v zN8p;mkM3tO-3NhqsMoiI;9cJblb}3|M_Sy6+DEoi&9%TMO&~N|PZlY7ehY92yO3PN z1kdjP4qZvFTWv5t%{p_B(Os}ebizd4B~yV{#0Cd?|-| z`AByA{B$)a5)QrnKoKw%mglz4wrmIfzu|v#S_V`xp448CYFi{5v^SXtRpLvF8@ShF z(sufl2eoxg$#4jx_Gra#UkZmh{^I%0>^l+IZV|xaI~aJLhhlpi4D6bhB7Jz|%yNKO z25{t&T50%?dw$*#@n0@K7Z?03UXqLw4IIn6*!c%^z2uJF)??nz68T?ObhgyaxJmpD zuc-w6Jx%^&m4(7;H{1Q74|OTLeb~Rl_{4rRhwdh4P6wUVz2)@R@IR0mbnHeO>!8-Q!VJoZ2w5UCUdz2{xCSNwJCLEp7bezwlds zaZba}JK)dCJAii{(!uZ_NI10MF^{(>%{PbTc8Am-3E9Dn_>>PPmWFGbjlT-|V{NU} z1}xB%wMaLZyXXjTPtSEOQ8E76OYMp*YLV^dN)~!cK^$5n>365JCtc6d{BVLTH5$ zLI|M<@q-`y$Pa$-1MPCE?pL>OcUSlA&90@8ZA;^GtInx&s!p9()xBdczWVxeYF=n_ppSK+Jt0Et>dah!k4Ld9}f7)3NYGLzIk)tLV_=_VjMas;pCh2zBXCabP` zP8hW*i@p$LMV_E@DYPIZ*kv9wq)At)xmxCR76b|Xkr}sWuddnLq}c?OwRqw}wr&@g zk&8xUaXCE`8}XDi!)TD_FgLUuxd}yX1>9~3!h!ZnCfYS2d{!zmCI`^p{RaoMQ)Mo7 zL~Q&uhgCyv6>V<0qHL8p%`sUAkvN<*8Dj86|G|C7oNABBQo@5XQ91Ify^1`c5m_7# z&T?KWl0stBX!pu(m*&Nyx@Njx_mVK(_jTRZs575dmJ?pym)06euc~vu8f#j&l!k4$ z7|HO*&wKTFM$96PetxI0O!XB6SpQ{Co%UXIY?AD*QpNaj*~wrf>*sCWT=cvb{fZfSOY zQq=^Mk@Ubv6eb{Y7TDxQ-FS=elTxV_#>y(l$!VtYQnKwPI!}bElk^fYQl~HHq=!5z zMK|K5R=BSQc@(HCI%u6B0_jzh1Gvn9HJYW7a1Rl#^3+;AZmKX&CYUT!87ImhkadFY zxp$VFKXJNoYPm(10!fm9W_?j82*)-gpf)N??oZ7|r4m_YNcKg6Rtm9lwp(zbz~&{3 zBDxz!lyxj2gZmooRmf7tc6c^=D+)~8a~zm~8}jgl3UvkCGCSMTo#?E!JanPrza7w4 z8AILEZR{{n?%_7Msj%96fviFJ+}0K81rbrU3qwOyaXN*|tf?;07Ua3_zpdTjRaS9q z)*I`s&R~C3(;1e2M&pKA_?}l~meSHM7HA7tfio*dj^pwQ<*;+TOFZ9)E^E2>UvwJW z?Xdd?F3LtgrpRn6Aa&}Sz35X#MFJ>0-w=!lxUWjA0&|?mWh$RjMt5ob?Uz zh}PGk_P*C>F|@r}rLWPUEON;EBV7pXdyPU#UNjQ^4^W3oew7!Tw#C*_?d&CP2f|Ud zZpjoiFmLU@B;!KoTn=86u|%hLoy!iOS7m42dGZp~H@0%@4AXw7O=RNrmEa0*brD7>^Rkf{I>M_p@F(yXn5l7l1 zLOQ~8IY0+DiAoiB`&3)HM4NyPi^U`Z?!U~c6%mnnw6o%#%WN%)BL^>Y3JRW#r48J%!Gfe1MLTIp{vk7t- z13LJPGM#Y?QO&wQ%F|)`&^LT}{+t!gnve(Gqz;Pv4N_0Vm|;X7jY$_PbfgotZUG@9 z=O=gWj+h6&$voxdmDEPuKgV1;l=IZ?J~_v{s9>+U($RA&TFT=RMa6n5emy?N{JOm= z$5FK?w~fw!p|wi8mIQ=3;-HEh!F5P^Dd{8AOP{Oc)e)W9hB#J@OBQ5WZPgu{bL5gH ztQm*DMVmkxqjhSd(nzXenC~(cb|V(J-%x=xnKuxr6X5~fSGU_ics%Y*>F)5HYkQNH zC|Qyp_bWUp5awJ_a;T`%^#rw@fHH>MUqh+WiO$NKPe^i#NmcGl$G~vhGSFm}P;&_K zm<)-VrCLN$ks}5sZNKRgDI+f}AkXduSs78icukMOY)0rwdjT94BJ(gHEIb||BPgPgB zPK$fzxt&fXPtH^6lX;hR1XKo_I}Ow@)^>KDQkBY~t}QQ2sH|2M(V8O{M>RW_PdZXO zG#`1kDoqu2UZx^Hxgol5K5B5^P4eCqCX0lE@ZEyQTQ&PFt%VqIk@lsW=o6h);rs1+ z?EbhjaP8H}PP)z;&HJv9vi;!62)35(Q_b6ndF|~ zDiIf{*06F*Jfs4KB|-Np9h9=)Asp3YreLaZLj}{Ts*0xh)8XHhCq8*wGG$dbg}daA zD<)b~RkEtAae<<1kXP-UIah|E`?3(V=CDU;p5Z-LG@rywo4&CEWA9l*7J)JGQ{3(Z$K4=eAnT>zV>F_`_qx_3Vr=|Y>`ez@}mQz zm^3m=%w9>HI_n-;WO?T2O`15>MrGRWA}iXffx&h)BJ=FK$f`V)Vb#h}Z%CybrN^=| zNga^*qPo~Bldm@9HJo%yeT^sYH`K+KJhsb(su4U+trJ#7-kLQ?ip<>Vrb0TmC6j3f z6A1I}mM}@LE?_OSICTq~l3Ti4a+CYOw6u#RX;do2SBO8L%HH^_%=)Qh%F11w`6|)M zTbSgi&F+#I;dTjn?XwCps`^UI&9Tes4q12!FX_lgNm_g)`Tw0=r1{q78rA`57mXiC-yhALMR#ag+*(?4_gNQTVnl5plNJut1 z9$T`>u7zut;??wCX`Jbz7lk>eLYXk6jz^ZNl>d=?98>N)CTC8RME$hO%|KRW*!O@HSMfjC{%1=`638JFVkr=Sr&Ji zix6RV(2nCe74q^ZEI)@xOpM}ekprM)+BOEp{UAi5sz0VT`7a^yNKL>i{S#!shkS@!(c zbwQT@*b^j*lGL2CtLK^YVk;RDpQ-y8Kv}hgQvARaaTbcLe*c>kIN^yf5Uf zeoac#lRR}@U0C7GL`}$R%s%gHCi{^ND%{R?49a>XbbbK!ASqVoIHXC#w~t1;!LMOf zq^mj}oth=R{Pc(Np3@FIhgVp?uNV_d<_nh$3AbzsC%%2LFyBQtCG6+3i9sV>N+Odm{kqBF0gR^&YU ze0+H(jm+Y7Xhn(il1t1)?IJ41UlPb~X{(aT5;ZEkJbCi9a{MhB;W`wYd@*-JI(U%> zRMZ*1f9wpxk986og7LGR6vG7N{Ws*NjGb)8RbQ3(Cod3h)a9p36=G2Bem^3N;py#@FSfHwixW7s+kXaSxCJPY^);3I&K0X`YSJ7El)%V>jc z6LL3o`w+(E0H49Q%YgGS+=9F7`^$hY#&G)s@OdSM?}G2UkP~>I z^JBP<{%>R6%|8JA4~+Yt7~TQDZR~gJtAH<~e+~M89@inK^>sjukD!+(c zu?hMn_%~r!zw{~SUIWzo3O#RP-#0Y0Vz`qB-J5_#%o}`n&~E-Q;D141n1_2ahPU>D zz64kST*dlfKRXWizaPWK$1$uyehvIK&tsmDvw?kVV&2=Z6R#P=Xe@@|M}Xk#K~5d| zxw9MNHvoSXazXzc$lrbk<21@;0lVK=z`WqcO{{MN<2Im&CfYlI>zGdy`{80; zF6Qw8)_)u8f9Dg7*Z$@L!oJ>tJ?>!Mx4>^3>)F)y^i4qcd+RD7*0KG&fbfsjA|T>R z3-;22U9_;@Ev)kn?BLFR$b}tlU|t)TPYZU@g1%ZI|@v++?3L#%WA*JHSYxWE1TG5j5@^EUK-2XTB8@;Ar84|EOi zZD5_PzX!y6TK@p}k1^cBdRmB&EuD8XT!0+#Z$jP{{9+4s-b8+A!jD?8-xl&!6ZX@D zpEl9I4cJ1yXuS#eZ-Acw{yQMzbn6!ZZ$dwiyZI^le;E*Vwf-^MzXAxmZM_SKc(IOk zhR|08{&fren9qlQ0*L*rzXE#DZT}kZe*ySskc0TtYR0hkehhDoK`(y=_`fmy7Ub<9 zPPM?d1$$}z1|aNt`?t{65dN`)xZi}|ZE1Obhc@fg`2oQN{=O)_Fe+&9=Y5gKz>|p=5HN-wg zn12NRH4E?+%-aAdLBK0k)B8XHG00=c6GON*%v9_jh;u(qvz4{5B=PE*K%*i7h$8;s5NShTBFt% zyLC4z@6Xfve)L-UJk4v=8ns5Pk++XpqgJ=&-u=w}#8nu3WTb=(2*&dhG zI4%&!_Mmg0{GZ9*keH5+Li=BAJ`#8MQ7!LqmbV={)c-j5L51=vBJe6K`K8kRsyE0s6V;CjIsO?|<(*-IFaiEYOqPrPmeaoi~d5SL%0B3yj%A(d?a04g$Rj{pDw diff --git a/data/soils_obic.RData b/data/soils_obic.RData index 64d35dc69c341436d7778a0bb026ff701d85cd6f..b866c77364ad432074a92407184da1466a7a690d 100644 GIT binary patch literal 233 zcmVP?xG=)#Yn4hJs5spLumv)vWR5nl z74Aujl(1ADAR+8TKD4*7u|TC(37(eIpG%X-S}BZ1T&+xhR>n6di;-6_ixr; zL;S++`UZb*%FB5sv=wBQ9Gdl2*Fw+&C=ZHbd&3QtcRG$rbqT2J)O^w%U0v7DuJlfF j-IeZFv5bSc_Vp!N@X5M?9*U5CZjaOkD0ssTZ!%}H&@${mw%8!`*+)ut*(uyK;mIkqK_(%6!(V!$SH?ubIJ+7`F|EufkZQ~c> zg4!;?tT;9+*0&NZ9~mB%r=B;c)HQPcnepE#bzS%D57|qh-4JfYVV%k(J6Ohvxev}n VEuuq34- Date: Wed, 10 Jul 2019 11:20:43 +0200 Subject: [PATCH 06/10] Adds function to calculate nlv --- NAMESPACE | 1 + R/nitrogen.R | 86 ++++++++++++++++++++++++++++++++++ R/ph.R | 6 +-- man/calc_nlv.Rd | 26 ++++++++++ man/crops.obic.Rd | 1 + man/soils.obic.Rd | 1 + tests/testthat/test-nitrogen.R | 6 ++- 7 files changed, 122 insertions(+), 5 deletions(-) create mode 100644 R/nitrogen.R create mode 100644 man/calc_nlv.Rd diff --git a/NAMESPACE b/NAMESPACE index 55c639ff..1d103d4a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # 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) diff --git a/R/nitrogen.R b/R/nitrogen.R new file mode 100644 index 00000000..54240268 --- /dev/null +++ b/R/nitrogen.R @@ -0,0 +1,86 @@ +#' 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 = 100, any.missing = FALSE, len = arg.length) + checkmate::assert_numeric(c.org, lower = 0, upper = 500, 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 = 10, 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 == "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) +} diff --git a/R/ph.R b/R/ph.R index 532ff918..f8c3ca84 100644 --- a/R/ph.R +++ b/R/ph.R @@ -43,7 +43,7 @@ calc_ph_delta <- function(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, c } # Collect information in table - soil.ph = lutum.low = lutum.high = om.low = om.high = potato.low = potato.high = sugarbeet.low = sugarbeet.high = ph.optimum = id = NULL + 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, @@ -63,8 +63,8 @@ calc_ph_delta <- function(ph.cacl2, soiltype, om, lutum, cp.starch, cp.potato, c dt <- merge(dt, soils.obic, by = "soiltype") # Define which table to be used - dt[soil.ph == 1, table := "5.1"] - dt[soil.ph == 2, table := "5.3"] + 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"] diff --git a/man/calc_nlv.Rd b/man/calc_nlv.Rd new file mode 100644 index 00000000..8130cbb4 --- /dev/null +++ b/man/calc_nlv.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/nitrogen.R +\name{calc_nlv} +\alias{calc_nlv} +\title{Calculate the NLV} +\usage{ +calc_nlv(n.org, c.org, crop, soiltype, bulk_density, cn_ratio, grass.age) +} +\arguments{ +\item{n.org}{(numeric) The organic nitrogen content of the soil in g N / kg} + +\item{c.org}{(numeric) The organic carbon content of the soil in kg C / ha} + +\item{crop}{(numeric) The crop code (gewascode) from the BRP} + +\item{soiltype}{(character) The type of soil} + +\item{bulk_density}{(numeric) The bulk density of the soil in kg / m3} + +\item{cn_ratio}{(numeric)} + +\item{grass.age}{(numeric) The age of the grass if present} +} +\description{ +This function calculates the NLV (nitrogen producing capacity) for the soil +} diff --git a/man/crops.obic.Rd b/man/crops.obic.Rd index 3b17615c..82d1e7fb 100644 --- a/man/crops.obic.Rd +++ b/man/crops.obic.Rd @@ -11,6 +11,7 @@ \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} }} \usage{ crops.obic diff --git a/man/soils.obic.Rd b/man/soils.obic.Rd index ac90780c..84cd9291 100644 --- a/man/soils.obic.Rd +++ b/man/soils.obic.Rd @@ -8,6 +8,7 @@ \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} }} \usage{ soils.obic diff --git a/tests/testthat/test-nitrogen.R b/tests/testthat/test-nitrogen.R index 6e5e5cd8..1ba9fdf4 100644 --- a/tests/testthat/test-nitrogen.R +++ b/tests/testthat/test-nitrogen.R @@ -6,9 +6,10 @@ test_that("calc_nlv works", { crop = 265, soiltype = "klei", bulk_density = 1300, + cn_ratio = 4, grass.age = 2 ), - expected = 0.3, + expected = 204.22254, tolerance = 0.001 ) expect_equal( @@ -18,9 +19,10 @@ test_that("calc_nlv works", { crop = c(265, 265, 235, 235), soiltype = c("klei", "veen", "veen", "zavel"), bulk_density = c(1300, 800, 850, 1100), + cn_ratio = c(4, 4, 4, 4), grass.age = c(2, 8, 0, 0) ), - expected = c(0.3, 0, 0, 0), + expected = c(204.22254, 250, 209.43603, 52.31861), tolerance = 0.001 ) }) \ No newline at end of file From 65f163736f9b3960b10b74d1cd99a5c05acd19a8 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Wed, 10 Jul 2019 11:26:40 +0200 Subject: [PATCH 07/10] Fix for maximum nlv for zand and klei --- R/nitrogen.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/nitrogen.R b/R/nitrogen.R index 54240268..1e0a84b6 100644 --- a/R/nitrogen.R +++ b/R/nitrogen.R @@ -69,6 +69,8 @@ calc_nlv <- function(n.org, c.org, crop, soiltype, bulk_density, cn_ratio, grass 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] + 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 From abfcf16fc182bcafad77acfedb714a5e1c5a0c79 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Wed, 10 Jul 2019 12:06:47 +0200 Subject: [PATCH 08/10] Adds evaluation for nitrogen --- NAMESPACE | 2 ++ R/evaluate.R | 25 +++++++++++++++++++++++++ R/nitrogen.R | 19 +++++++++++++++++++ man/eval_nitrogen.Rd | 14 ++++++++++++++ man/evaluate_parabolic.Rd | 16 ++++++++++++++++ tests/testthat/test-nitrogen.R | 17 +++++++++++++++++ 6 files changed, 93 insertions(+) create mode 100644 man/eval_nitrogen.Rd create mode 100644 man/evaluate_parabolic.Rd diff --git a/NAMESPACE b/NAMESPACE index 1d103d4a..6cc9405c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,11 +6,13 @@ 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) diff --git a/R/evaluate.R b/R/evaluate.R index e3efcbe4..1874479d 100644 --- a/R/evaluate.R +++ b/R/evaluate.R @@ -29,3 +29,28 @@ evaluate_logistic <- function(x, b, x0, v, increasing = TRUE) { 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) + +} \ No newline at end of file diff --git a/R/nitrogen.R b/R/nitrogen.R index 1e0a84b6..dd53fb52 100644 --- a/R/nitrogen.R +++ b/R/nitrogen.R @@ -86,3 +86,22 @@ calc_nlv <- function(n.org, c.org, crop, soiltype, bulk_density, cn_ratio, grass 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) +} diff --git a/man/eval_nitrogen.Rd b/man/eval_nitrogen.Rd new file mode 100644 index 00000000..482610da --- /dev/null +++ b/man/eval_nitrogen.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/nitrogen.R +\name{eval_nitrogen} +\alias{eval_nitrogen} +\title{Evaluate nitrogen} +\usage{ +eval_nitrogen(value.nlv) +} +\arguments{ +\item{value.nlv}{(numeric) The value of NLV calculated by \code{\link{calc_nlv}}} +} +\description{ +This function evaluates the nitrogen content of the soil by using the NLV calculated by \code{\link{calc_nlv}} +} diff --git a/man/evaluate_parabolic.Rd b/man/evaluate_parabolic.Rd new file mode 100644 index 00000000..beb45c37 --- /dev/null +++ b/man/evaluate_parabolic.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/evaluate.R +\name{evaluate_parabolic} +\alias{evaluate_parabolic} +\title{Evaluate using parabolic function with} +\usage{ +evaluate_parabolic(x, x.top) +} +\arguments{ +\item{x}{(numeric) The values of a calc function to be converted to an evaluation} + +\item{x.top}{(numeric) The value at which x reaches the plateau} +} +\description{ +This function evaluates the calculated values from an indicator using a parabolic function. After the optimum is reached the it stays at its plateau. +} diff --git a/tests/testthat/test-nitrogen.R b/tests/testthat/test-nitrogen.R index 1ba9fdf4..44a06aaf 100644 --- a/tests/testthat/test-nitrogen.R +++ b/tests/testthat/test-nitrogen.R @@ -25,4 +25,21 @@ test_that("calc_nlv works", { expected = c(204.22254, 250, 209.43603, 52.31861), tolerance = 0.001 ) +}) + +test_that("eval_nitrogen works", { + expect_equal( + eval_nitrogen( + value.nlv = 120 + ), + expected = 1, + tolerance = 0.001 + ) + expect_equal( + eval_nitrogen( + value.nlv = seq(from = -30, to = 250, by = 50) + ), + expected = c(0, 0.3055556, 0.8263889, 1, 1, 1), + tolerance = 0.001 + ) }) \ No newline at end of file From adfb9fdfbdb081a71b680971c9540a691cb41abc Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Fri, 12 Jul 2019 09:24:51 +0200 Subject: [PATCH 09/10] Add badge to readme --- README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/README.md b/README.md index 7b865a15..63db08b7 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,8 @@ # Open Bodem Index Calculator (OBIC) + +[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) + + 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. From 0abb80583f8c1350cdce37ffb031f6ade00ef648 Mon Sep 17 00:00:00 2001 From: Sven <37927107+SvenVw@users.noreply.github.com> Date: Fri, 12 Jul 2019 09:41:20 +0200 Subject: [PATCH 10/10] Improved checks for n.org and c.org --- R/nitrogen.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/nitrogen.R b/R/nitrogen.R index dd53fb52..1bca2987 100644 --- a/R/nitrogen.R +++ b/R/nitrogen.R @@ -23,14 +23,14 @@ calc_nlv <- function(n.org, c.org, crop, soiltype, bulk_density, cn_ratio, grass # 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 = 100, any.missing = FALSE, len = arg.length) - checkmate::assert_numeric(c.org, lower = 0, upper = 500, any.missing = FALSE, len = arg.length) + 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 = 10, 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