Support units in UK drainage calculations

Also, address some R CMD check notes

13 files changed, 79 insertions, 44 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 084818b..16f04c2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: pfm
 Type: Package
 Title: Utilities for Pesticide Fate Modelling
 Version: 0.6.4
-Date: 2025-06-24
+Date: 2026-02-06
 Authors@R: c(
   person("Johannes Ranke", email = "johannes.ranke@agroscope.admin.ch",
     role = c("aut", "cre"),
@@ -37,4 +37,4 @@ LazyData: true
 Encoding: UTF-8
 URL: https://pkgdown.jrwb.de/pfm, https://github.com/jranke/pfm, http://jranke.github.io/pfm/
 Roxygen: list(markdown = TRUE, r6 = TRUE)
-RoxygenNote: 7.3.2.9000
+RoxygenNote: 7.3.3
diff --git a/R/PEC_sw_drainage_UK.R b/R/PEC_sw_drainage_UK.R
index 4b3111e..599e15e 100644
--- a/R/PEC_sw_drainage_UK.R
+++ b/R/PEC_sw_drainage_UK.R
@@ -4,11 +4,14 @@
 #' published on the CRC website. Degradation before the start of the drainage period is taken into account if
 #' `latest_application` is specified and the degradation parameters are given either as a `soil_DT50` or a `model`.
 #'
-#' @param rate Application rate in g/ha
+#' @param rate Application rate in g/ha or with a compatible unit specified
+#' with the units package
 #' @param interception The fraction of the application rate that does not reach the soil
-#' @param Koc The sorption coefficient normalised to organic carbon in L/kg
+#' @param Koc The sorption coefficient normalised to organic carbon in L/kg or a unit specified
+#' with the units package
 #' @param latest_application Latest application date, formatted as e.g. "01 July"
-#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used
+#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used, in
+#' days or a time unit specified with the units package
 #' @param model The soil degradation model to be used. Either one of "FOMC",
 #'   "DFOP", "HS", or "IORE", or an mkinmod object
 #' @param model_parms A named numeric vector containing the model parameters
@@ -27,27 +30,41 @@
 #' PEC_sw_drainage_UK(150, Koc = 100)
 #' PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550,
 #'   latest_application = "01 July", soil_DT50 = 200)
-
-PEC_sw_drainage_UK <- function(rate, interception = 0, Koc,
-                                   latest_application = NULL, soil_DT50 = NULL,
-                                   model = NULL, model_parms = NULL)
+PEC_sw_drainage_UK <- function(rate, 
+  interception = 0, Koc,
+  latest_application = NULL, soil_DT50 = NULL,
+  model = NULL, model_parms = NULL)
 {
-  percentage_lost <- SSLRC_mobility_classification(Koc)[[2]]
-  amount_available <- rate * (1 - interception) # g/ha
+  # Set default units if not specified and convert to units used in the calculations
+  if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha")
+  rate_g_ha <- as.numeric(set_units(rate, "g/ha"))
+
+  if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg")
+  Koc_L_kg <- as.numeric(set_units(Koc, "L/kg"))
+
+  if (!missing(soil_DT50)) {
+    if (!inherits(soil_DT50, "units")) {
+      soil_DT50_d <- soil_DT50
+    }
+    soil_DT50_d <- as.numeric(set_units(soil_DT50, "d"))
+  }
+
+  percentage_lost <- SSLRC_mobility_classification(Koc_L_kg)[[2]]
+  amount_available <- rate_g_ha * (1 - interception) # amount in g for 1 ha
 
   if (!missing(latest_application)) {
     lct <- Sys.getlocale("LC_TIME")
     tmp <- Sys.setlocale("LC_TIME", "C")
-    if (latest_application == "29 February") {
+    if (latest_application == "29 February") { # Use a leap year
       ref_year <- 2000
-    } else { ref_year <- 1999}
+    } else { ref_year <- 1999} # Use a non-leap year
     latest <- as.Date(paste(latest_application, ref_year), "%d %b %Y")
     if (is.na(latest)) stop("Please specify the latest application in the format '%d %b', e.g. '01 July'")
     tmp <- Sys.setlocale("LC_TIME", lct)
     degradation_time <- as.numeric(difftime(as.Date(paste0(ref_year,"-10-01")), units = "days", latest))
     if (degradation_time > 0) {
       if (!missing(soil_DT50)) {
-        k = log(2)/soil_DT50
+        k = log(2)/soil_DT50_d
         as.Date(paste(latest_application, "1999"), "%d %B %Y")
 
         amount_available <- amount_available * exp(-k * degradation_time)
@@ -62,6 +79,6 @@ PEC_sw_drainage_UK <- function(rate, interception = 0, Koc,
   } 
 
   volume = 130000 # L/ha
-  PEC = 1e6 * (percentage_lost/100) * amount_available / volume
+  PEC = set_units(1e6 * (percentage_lost/100) * amount_available / volume, "\u00B5g/L")
   return(PEC)
 }
diff --git a/R/PEC_sw_exposit.R b/R/PEC_sw_exposit.R
index 282a1ac..a157189 100644
--- a/R/PEC_sw_exposit.R
+++ b/R/PEC_sw_exposit.R
@@ -52,7 +52,8 @@
 #'
 #' @importFrom units as_units set_units drop_units
 #' @importFrom dplyr across mutate
-#' @param rate The application rate in g/ha
+#' @param rate Application rate in g/ha or with a compatible unit specified
+#' with the units package
 #' @param interception The fraction intercepted by the crop
 #' @param Koc The sorption coefficient to soil organic carbon
 #' @param DT50 The soil half-life in days
diff --git a/R/TOXSWA_cwa.R b/R/TOXSWA_cwa.R
index 310029e..b2f7619 100644
--- a/R/TOXSWA_cwa.R
+++ b/R/TOXSWA_cwa.R
@@ -147,10 +147,10 @@ plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t
 #' and some associated statistics.  like maximum moving window average
 #' concentrations, and dataframes holding the events exceeding specified
 #' thresholds.  Usually, an instance of this class will be generated
-#' by \code{\link{read.TOXSWA_cwa}}.
+#' by [read.TOXSWA_cwa].
 #'
 #' @export
-#' @format An \code{\link{R6Class}} generator object.
+#' @format An [R6::R6Class] generator object.
 #' @field filename Length one character vector holding the filename.
 #' @field basedir Length one character vector holding the directory where the file came from.
 #' @field zipfile If not null, giving the path to the zip file from which the file was read.
diff --git a/R/twa.R b/R/twa.R
index 2018dfa..c756f58 100644
--- a/R/twa.R
+++ b/R/twa.R
@@ -3,9 +3,9 @@
 #' @param x When numeric, this is the half-life to be used for an exponential
 #'   decline. When a character string specifying a parent decline model is given
 #'   e.g. \code{FOMC}, \code{parms} must contain the corresponding parameters.
-#'   If x is an \code{\link{mkinfit}} object, the decline is calculated from this
+#'   If x is an [mkin::mkinfit] object, the decline is calculated from this
 #'   object.
-#' @param ini The initial amount. If x is an \code{\link{mkinfit}} object, and
+#' @param ini The initial amount. If x is an [mkin::mkinfit] object, and
 #'   ini is 'model', the fitted initial concentrations are used. Otherwise, ini
 #'   must be numeric. If it has length one, it is used for the parent and
 #'   initial values of metabolites are zero, otherwise, it must give values for
@@ -246,7 +246,7 @@ twa.one_box <- function(x, window = 21)
 #' \code{\link{plot.one_box}} using the window size for the argument
 #' \code{max_twa}.
 #'
-#' The method working directly on fitted \code{\link{mkinfit}} objects uses the
+#' The method working directly on fitted [mkin::mkinfit]] objects uses the
 #' equations given in the PEC soil section of the FOCUS guidance and is restricted
 #' SFO, FOMC and DFOP models and to the parent compound
 #' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
diff --git a/log/check.log b/log/check.log
index bd2b822..9ed80ce 100644
--- a/log/check.log
+++ b/log/check.log
@@ -1,10 +1,10 @@
 * using log directory ‘/home/jranke/git/pfm/pfm.Rcheck’
-* using R version 4.4.2 (2024-10-31)
+* using R version 4.5.2 (2025-10-31)
 * using platform: x86_64-pc-linux-gnu
 * R was compiled by
-    gcc (Debian 12.2.0-14) 12.2.0
-    GNU Fortran (Debian 12.2.0-14) 12.2.0
-* running under: Debian GNU/Linux 12 (bookworm)
+    gcc (Debian 14.2.0-19) 14.2.0
+    GNU Fortran (Debian 14.2.0-19) 14.2.0
+* running under: Debian GNU/Linux 13 (trixie)
 * using session charset: UTF-8
 * using options ‘--no-tests --as-cran’
 * checking for file ‘pfm/DESCRIPTION’ ... OK
@@ -14,7 +14,7 @@
 * checking CRAN incoming feasibility ... NOTE
 Maintainer: ‘Johannes Ranke <johannes.ranke@agroscope.admin.ch>’
 
-Size of tarball: 8532625 bytes
+Size of tarball: 8533567 bytes
 * checking package namespace information ... OK
 * checking package dependencies ... OK
 * checking if this is a source package ... OK
@@ -24,7 +24,7 @@ Size of tarball: 8532625 bytes
 * checking for portable file names ... OK
 * checking for sufficient/correct file permissions ... OK
 * checking whether package ‘pfm’ can be installed ... OK
-* checking installed package size ... NOTE
+* checking installed package size ... INFO
   installed size is 10.3Mb
   sub-directories of 1Mb or more:
     testdata   9.9Mb
@@ -72,7 +72,7 @@ unable to verify current time
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: 3 NOTEs
+Status: 2 NOTEs
 See
   ‘/home/jranke/git/pfm/pfm.Rcheck/00check.log’
 for details.
diff --git a/log/test.log b/log/test.log
index 815e792..7204076 100644
--- a/log/test.log
+++ b/log/test.log
@@ -9,10 +9,10 @@
 ✔ |          1 | Actual and time weighted average concentrations for SFO kinetics
 ✔ |          9 | FOCUS Step 1 calculations
 ✔ |          8 | FOCUS Steps 12 input files
-✔ |          7 | Read and analyse TOXSWA cwa files [2.8s]
-✔ |         15 | UK drainage PEC calculations
+✔ |          7 | Read and analyse TOXSWA cwa files [2.6s]
+✔ |         17 | UK drainage PEC calculations
 
 ══ Results ═════════════════════════════════════════════════════════════════════════════════════════
 Duration: 3.8 s
 
-[ FAIL 0 | WARN 0 | SKIP 0 | PASS 80 ]
+[ FAIL 0 | WARN 0 | SKIP 0 | PASS 82 ]
diff --git a/man/PEC_sw_drainage_UK.Rd b/man/PEC_sw_drainage_UK.Rd
index 6125c3c..af020c6 100644
--- a/man/PEC_sw_drainage_UK.Rd
+++ b/man/PEC_sw_drainage_UK.Rd
@@ -15,15 +15,18 @@ PEC_sw_drainage_UK(
 )
 }
 \arguments{
-\item{rate}{Application rate in g/ha}
+\item{rate}{Application rate in g/ha or with a compatible unit specified
+with the units package}
 
 \item{interception}{The fraction of the application rate that does not reach the soil}
 
-\item{Koc}{The sorption coefficient normalised to organic carbon in L/kg}
+\item{Koc}{The sorption coefficient normalised to organic carbon in L/kg or a unit specified
+with the units package}
 
 \item{latest_application}{Latest application date, formatted as e.g. "01 July"}
 
-\item{soil_DT50}{Soil degradation half-life, if SFO kinetics are to be used}
+\item{soil_DT50}{Soil degradation half-life, if SFO kinetics are to be used, in
+days or a time unit specified with the units package}
 
 \item{model}{The soil degradation model to be used. Either one of "FOMC",
 "DFOP", "HS", or "IORE", or an mkinmod object}
diff --git a/man/PEC_sw_exposit_runoff.Rd b/man/PEC_sw_exposit_runoff.Rd
index 393ed29..60c9b9f 100644
--- a/man/PEC_sw_exposit_runoff.Rd
+++ b/man/PEC_sw_exposit_runoff.Rd
@@ -21,7 +21,8 @@ PEC_sw_exposit_runoff(
 )
 }
 \arguments{
-\item{rate}{The application rate in g/ha}
+\item{rate}{Application rate in g/ha or with a compatible unit specified
+with the units package}
 
 \item{interception}{The fraction intercepted by the crop}
 
diff --git a/man/TOXSWA_cwa.Rd b/man/TOXSWA_cwa.Rd
index 0923f10..126fed4 100644
--- a/man/TOXSWA_cwa.Rd
+++ b/man/TOXSWA_cwa.Rd
@@ -4,14 +4,14 @@
 \alias{TOXSWA_cwa}
 \title{R6 class for holding TOXSWA water concentration data and associated statistics}
 \format{
-An \code{\link{R6Class}} generator object.
+An \link[R6:R6Class]{R6::R6Class} generator object.
 }
 \description{
 An R6 class for holding TOXSWA water concentration (cwa) data
 and some associated statistics.  like maximum moving window average
 concentrations, and dataframes holding the events exceeding specified
 thresholds.  Usually, an instance of this class will be generated
-by \code{\link{read.TOXSWA_cwa}}.
+by \link{read.TOXSWA_cwa}.
 }
 \examples{
 H_sw_R1_stream  <- read.TOXSWA_cwa("00003s_pa.cwa",
diff --git a/man/max_twa.Rd b/man/max_twa.Rd
index 43ad50e..0fb5c93 100644
--- a/man/max_twa.Rd
+++ b/man/max_twa.Rd
@@ -19,7 +19,7 @@ for finding the maximum. It is therefore recommended to check this using
 \code{max_twa}.
 }
 \details{
-The method working directly on fitted \code{\link{mkinfit}} objects uses the
+The method working directly on fitted \link[mkin:mkinfit]{mkin::mkinfit}] objects uses the
 equations given in the PEC soil section of the FOCUS guidance and is restricted
 SFO, FOMC and DFOP models and to the parent compound
 }
diff --git a/man/one_box.Rd b/man/one_box.Rd
index a60aa39..afa6e10 100644
--- a/man/one_box.Rd
+++ b/man/one_box.Rd
@@ -19,10 +19,10 @@ one_box(x, ini, ..., t_end = 100, res = 0.01)
 \item{x}{When numeric, this is the half-life to be used for an exponential
 decline. When a character string specifying a parent decline model is given
 e.g. \code{FOMC}, \code{parms} must contain the corresponding parameters.
-If x is an \code{\link{mkinfit}} object, the decline is calculated from this
+If x is an \link[mkin:mkinfit]{mkin::mkinfit} object, the decline is calculated from this
 object.}
 
-\item{ini}{The initial amount. If x is an \code{\link{mkinfit}} object, and
+\item{ini}{The initial amount. If x is an \link[mkin:mkinfit]{mkin::mkinfit} object, and
 ini is 'model', the fitted initial concentrations are used. Otherwise, ini
 must be numeric. If it has length one, it is used for the parent and
 initial values of metabolites are zero, otherwise, it must give values for
diff --git a/tests/testthat/test_UK_drainage.R b/tests/testthat/test_UK_drainage.R
index defecbb..75d9c98 100644
--- a/tests/testthat/test_UK_drainage.R
+++ b/tests/testthat/test_UK_drainage.R
@@ -24,15 +24,23 @@ test_that("UK drainflow PECs are correct", {
 
   # This is the first example calculation from the data requirements handbook, where they give
   # 8.07 µg/L as the result (obviously a rounding error).
-  expect_equal(round(PEC_sw_drainage_UK(150, interception = 0, Koc = 100), 4), 8.0769)
+  expect_equal(round(PEC_sw_drainage_UK(150, interception = 0, Koc = 100), 4), 
+    as_units(8.0769, "\u00B5g/L"))
+
+  # Same as above, but using units for the input
+  expect_equal(round(PEC_sw_drainage_UK(as_units(0.15, "kg/ha"), interception = 0, 
+        Koc = as_units(100, "mL/g")), 4), 
+    as_units(8.0769, "\u00B5g/L"))
 
   # This is the second example calculation from the data requirements handbook
-  expect_equal(round(PEC_sw_drainage_UK(90, interception = 0, Koc = 10), 4), 13.1538)
+  expect_equal(round(PEC_sw_drainage_UK(90, interception = 0, Koc = 10), 4),
+    as_units(13.1538, "\u00B5g/L"))
 
   # This is the third example calculation from the data requirements handbook, 
   expect_equal(round(PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550,
-                                        latest_application = "01 July",
-                                        soil_DT50 = 200), 2), 0.84)
+    latest_application = "01 July",
+    soil_DT50 = 200), 2), 
+    as_units(0.84, "\u00B5g/L"))
 
   expect_error(round(PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550,
                                         latest_application = "100 July",
@@ -42,6 +50,11 @@ test_that("UK drainflow PECs are correct", {
                                         latest_application = "29 February",
                                         soil_DT50 = 200), 2))
 
+  expect_silent(round(PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550,
+                                        latest_application = "29 February",
+                                        soil_DT50 = as_units(200, "d")), 2))
+  
+
   # Test that PECsw do not increase if the application is after the beginning
   # of the drainflow period
   expect_equal(