PecanProject · mdietze · Aug 1, 2022 · Jul 27, 2022 · Jul 27, 2022 · Jul 27, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -39,6 +39,11 @@ We are slowly change the license from NCSA opensource to BSD-3 to help with publ
 - when specifying diferent rstudio user, dev setup would mount pecan folder in wrong path.
 - Code cleanup in PEcAn.MA to protect against global namespace pollution (#2965, #2973; @nanu1605) 
 - Fixed vignettes and cleaned up internal warnings in PEcAn.DB (#2966, #2971; @nanu1605).
+- Updated unit conversions throughout PEcAn to use the `units` R package
+  instead of the unmaintained `udunits2`.
+  Note that both `units` and `udunits2` interface with the same underlying
+  compiled code, so the `udunits2` *system library* is still required.
+  (#2989; @nanu1605)
 
 ### Changed
 

diff --git a/Makefile.depends b/Makefile.depends
@@ -27,13 +27,13 @@ $(call depends,models/biocro): | .install/base/logger .install/base/remote .inst
 $(call depends,models/cable): | .install/base/logger .install/base/utils
 $(call depends,models/clm45): | .install/base/logger .install/base/utils
 $(call depends,models/dalec): | .install/base/logger .install/base/remote .install/base/utils
-$(call depends,models/dvmdostem): | .install/base/utils
+$(call depends,models/dvmdostem): | .install/base/logger .install/base/utils
 $(call depends,models/ed): | .install/modules/data.atmosphere .install/modules/data.land .install/base/logger .install/base/remote .install/base/settings .install/base/utils
-$(call depends,models/fates): | .install/base/utils .install/base/logger .install/base/remote
+$(call depends,models/fates): | .install/base/logger .install/base/remote .install/base/utils
 $(call depends,models/gday): | .install/base/utils .install/base/logger .install/base/remote
-$(call depends,models/jules): | .install/base/utils .install/base/logger .install/base/remote
+$(call depends,models/jules): | .install/modules/data.atmosphere .install/base/logger .install/base/remote .install/base/utils
 $(call depends,models/linkages): | .install/base/utils .install/modules/data.atmosphere .install/base/logger .install/base/remote
-$(call depends,models/lpjguess): | .install/base/utils .install/base/logger .install/base/remote
+$(call depends,models/lpjguess): | .install/base/logger .install/base/remote .install/base/utils
 $(call depends,models/maat): | .install/modules/data.atmosphere .install/base/logger .install/base/remote .install/base/settings .install/base/utils
 $(call depends,models/maespa): | .install/modules/data.atmosphere .install/base/logger .install/base/remote .install/base/utils
 $(call depends,models/preles): | .install/base/utils .install/base/logger .install/modules/data.atmosphere .install/base/utils

diff --git a/base/db/DESCRIPTION b/base/db/DESCRIPTION
@@ -58,7 +58,7 @@ Imports:
     rlang,
     tibble,
     tidyr,
-    udunits2,
+    units
 Suggests:
     bit64,
     data.table,

diff --git a/base/db/R/insert.format.vars.R b/base/db/R/insert.format.vars.R
@@ -126,13 +126,15 @@ insert.format.vars <- function(con, format_name, mimetype_id, notes = NULL, head
       u1 <- formats_variables[1,"unit"]
       u2 <- dplyr::tbl(con, "variables") %>% dplyr::select(.data$id, units) %>% dplyr::filter(.data$id %in% !!formats_variables[[1, "variable_id"]]) %>% dplyr::pull(.data$units)
 
-      if(!udunits2::ud.is.parseable(u1)){
+      u1_recognized <- tryCatch(units::as_units(u1), error = function(e) FALSE)
+      if(!u1_recognized){
         PEcAn.logger::logger.error(
-          "Units not parseable. Please enter a unit that is parseable by the udunits library."
+          "Units '", u1,  "' not parseable.",
+          "Please provide a unit that is parseable by the udunits library."
         )
       }
       # Grab the bety units and
-      if(!udunits2::ud.are.convertible(u1, u2)){
+      if(!units::ud_are_convertible(u1, u2)){
         PEcAn.logger::logger.error(
           "Units are not convertable."
         )

diff --git a/base/db/R/query.dplyr.R b/base/db/R/query.dplyr.R
@@ -62,8 +62,8 @@ dplyr.count <- function(df) {
 #' @export
 ncdays2date <- function(time, unit) {
   date    <- lubridate::parse_date_time(unit, c("ymd_hms", "ymd_h", "ymd"))
-  days    <- udunits2::ud.convert(time, unit, paste("days since ", date))
-  seconds <- udunits2::ud.convert(days, "days", "seconds")
+  days    <- PEcAn.utils::ud_convert(time, unit, paste("days since ", date))
+  seconds <- PEcAn.utils::ud_convert(days, "days", "seconds")
   return(as.POSIXct.numeric(seconds, origin = date, tz = "UTC"))
 }  # ncdays2date
 

diff --git a/base/db/R/query.format.vars.R b/base/db/R/query.format.vars.R
@@ -133,7 +133,7 @@ query.format.vars <- function(bety, input.id=NA, format.id=NA, var.ids=NA) {
 
         # This would be a good place to put a test for valid sotrage types. Currently not implemented.
 
-      } else if (udunits2::ud.are.convertible(format$vars$input_units[i], format$vars$pecan_units[i]) == FALSE) {
+      } else if (units::ud_are_convertible(format$vars$input_units[i], format$vars$pecan_units[i]) == FALSE) {
 
         if (PEcAn.utils::misc.are.convertible(format$vars$input_units[i], format$vars$pecan_units[i]) == FALSE) {
           PEcAn.logger::logger.warn("Units not convertible for",format$vars$input_name[i], "with units of",format$vars$input_units[i], ".  Please make sure the varible has units that can be converted to", format$vars$pecan_units[i])

diff --git a/base/db/inst/import-try/96.load.data.R b/base/db/inst/import-try/96.load.data.R
@@ -1,7 +1,6 @@
 # 6. Add TRY data to BETY
 source("common.R")
 load("try.5.RData")
-library(udunits2)
 library(stringr)
 
 setkey(try.dat, ObservationID)
@@ -57,7 +56,7 @@ for(i in 1:nrow(try.entities)){
       site_id = try.sub[j, bety.site.id],
       specie_id = try.sub[j, bety.species.id],
       citation_id = try.sub[j, bety.citation.id],
-      mean = try.sub[j, udunits2::ud.convert(StdValue, UnitName, units)],
+      mean = try.sub[j, PEcAn.utils::ud_convert(StdValue, UnitName, units)],
       n = try.sub[j, Replicates],
       user_id = user_id,
       entity_id = entity,

diff --git a/base/utils/DESCRIPTION b/base/utils/DESCRIPTION
@@ -41,7 +41,7 @@ Imports:
     RCurl,
     rlang,
     stringi,
-    udunits2 (>= 0.11),
+    units,
     XML
 Suggests:
     coda (>= 0.18),

diff --git a/base/utils/NAMESPACE b/base/utils/NAMESPACE
@@ -54,6 +54,7 @@ export(to_ncvar)
 export(trait.lookup)
 export(transformstats)
 export(tryl)
+export(ud_convert)
 export(units_are_equivalent)
 export(vecpaste)
 export(zero.truncate)

diff --git a/base/utils/R/cf2date.R b/base/utils/R/cf2date.R
@@ -18,7 +18,7 @@
 #' cf2datetime(365, "days since 2000-01-01 12:00:00 -05:00")
 cf2datetime <- function(value, unit, tz = "UTC") {
   origin <- "1970-01-01 00:00:00 UTC"
-  ctint <- udunits2::ud.convert(value, unit, paste("seconds since", origin))
+  ctint <- ud_convert(value, unit, paste("seconds since", origin))
   result <- as.POSIXct(ctint, origin = origin)
   attr(result, "tzone") <- tz
   result
@@ -43,7 +43,7 @@ datetime2cf <- function(datetime, unit, ...) {
   }
   origin <- "1970-01-01 00:00:00 UTC"
   raw_ct <- as.numeric(datetime, origin = origin)
-  udunits2::ud.convert(raw_ct, paste("seconds since", origin), unit)
+  ud_convert(raw_ct, paste("seconds since", origin), unit)
 }
 
 #' Extract Julian day from CF or POSIXct date-times

diff --git a/base/utils/R/read.output.R b/base/utils/R/read.output.R
@@ -218,7 +218,7 @@ read.output <- function(runid, outdir,
       }
 
       if (dataframe) {
-        seconds <- udunits2::ud.convert(
+        seconds <- ud_convert(
           nc$dim$time$vals,
           nc$dim$time$units,
           paste("seconds since", run_origin)
@@ -264,7 +264,7 @@ read.output <- function(runid, outdir,
         # Dropping attempt to provide more sensible units because of graph unit errors,
         # issue #792
         # if (v %in% c(cflux, wflux)) {
-        #   newresult <- udunits2::ud.convert(newresult, 'kg m-2 s-1', 'kg ha-1 yr-1')
+        #   newresult <- ud_convert(newresult, 'kg m-2 s-1', 'kg ha-1 yr-1')
         # }
         result[[v]] <- abind::abind(result[[v]], newresult)
       }

diff --git a/base/utils/R/seconds_in_year.R b/base/utils/R/seconds_in_year.R
@@ -11,6 +11,6 @@
 #' @export
 seconds_in_year <- function(year, leap_year = TRUE, ...) {
     diy <- days_in_year(year, leap_year)
-    siy <- udunits2::ud.convert(diy, 'days', 'seconds')
+    siy <- ud_convert(diy, 'days', 'seconds')
     return(siy)
 }
diff --git a/base/utils/R/ud_convert.R b/base/utils/R/ud_convert.R
@@ -0,0 +1,18 @@
+##' Convert units
+##'
+##' Unit conversion to replace the now-unmaintained `udunits2::ud.convert`
+##' @author Chris Black
+##' 
+##' @param x numeric vector
+##' @param u1 string parseable as the units in which `x` is provided
+##' @param u2 string parseable as the units to convert to
+##' 
+##' @return numeric vector with values converted to units in `u2`
+##' @export
+ud_convert <- function(x, u1, u2) {
+  stopifnot(units::ud_are_convertible(u1, u2))
+  x1 <- units::set_units(x, value = u1, mode = "standard")
+  x2 <- units::set_units(x1, value = u2, mode = "standard")
+
+  units::drop_units(x2)
+  } # ud_convert
diff --git a/base/utils/R/units_are_equivalent.R b/base/utils/R/units_are_equivalent.R
@@ -8,6 +8,6 @@
 #' @author Alexey Shiklomanov
 #' @export
 units_are_equivalent <- function(x, y) {
-  x2y <- udunits2::ud.convert(1, x, y)
+  x2y <- ud_convert(1, x, y)
   1 == x2y
 }
diff --git a/base/utils/R/utils.R b/base/utils/R/utils.R
@@ -389,8 +389,8 @@ tabnum <- function(x, n = 3) {
 ##' @export
 ##' @author unknown
 arrhenius.scaling <- function(observed.value, old.temp, new.temp = 25) {
-  new.temp.K <- udunits2::ud.convert(new.temp, "degC", "K")
-  old.temp.K <- udunits2::ud.convert(old.temp, "degC", "K")
+  new.temp.K <- ud_convert(new.temp, "degC", "K")
+  old.temp.K <- ud_convert(old.temp, "degC", "K")
   return(observed.value / exp(3000 * (1 / (new.temp.K) - 1 / (old.temp.K))))
 } # arrhenius.scaling
 #--------------------------------------------------------------------------------------------------#
@@ -558,21 +558,21 @@ misc.convert <- function(x, u1, u2) {
   mmH2O <- 18.01528 # molar mass of H2O, g/mol
 
   if (u1 == "umol C m-2 s-1" & u2 == "kg C m-2 s-1") {
-    val <- udunits2::ud.convert(x, "ug", "kg") * amC
+    val <- ud_convert(x, "ug", "kg") * amC
   } else if (u1 == "kg C m-2 s-1" & u2 == "umol C m-2 s-1") {
-    val <- udunits2::ud.convert(x, "kg", "ug") / amC
+    val <- ud_convert(x, "kg", "ug") / amC
   } else if (u1 == "mol H2O m-2 s-1" & u2 == "kg H2O m-2 s-1") {
-    val <- udunits2::ud.convert(x, "g", "kg") * mmH2O
+    val <- ud_convert(x, "g", "kg") * mmH2O
   } else if (u1 == "kg H2O m-2 s-1" & u2 == "mol H2O m-2 s-1") {
-    val <- udunits2::ud.convert(x, "kg", "g") / mmH2O
+    val <- ud_convert(x, "kg", "g") / mmH2O
   } else if (u1 == "Mg ha-1" & u2 == "kg C m-2") {
-    val <- x * udunits2::ud.convert(1, "Mg", "kg") * udunits2::ud.convert(1, "ha-1", "m-2")
+    val <- x * ud_convert(1, "Mg", "kg") * ud_convert(1, "ha-1", "m-2")
   } else if (u1 == "kg C m-2" & u2 == "Mg ha-1") {
-    val <- x * udunits2::ud.convert(1, "kg", "Mg") * udunits2::ud.convert(1, "m-2", "ha-1")
+    val <- x * ud_convert(1, "kg", "Mg") * ud_convert(1, "m-2", "ha-1")
   } else {
     u1 <- gsub("gC","g*12",u1)
     u2 <- gsub("gC","g*12",u2)
-    val <- udunits2::ud.convert(x,u1,u2)
+    val <- ud_convert(x,u1,u2)
 
 
     #    PEcAn.logger::logger.severe(paste("Unknown units", u1, u2))

diff --git a/base/utils/man/ud_convert.Rd b/base/utils/man/ud_convert.Rd
diff --git a/base/utils/tests/testthat/test-ud_convert.R b/base/utils/tests/testthat/test-ud_convert.R
@@ -0,0 +1,27 @@
+test_that("unit conversions", {
+  expect_equal(ud_convert(4, "in", "cm"), 10.16)
+  expect_equal(ud_convert(23, "degC", "K"), 296.15)
+  expect_equal(ud_convert(6, "days", "seconds"), 518400)
+  expect_equal(ud_convert(32, "hour", "seconds"), 115200)
+  expect_equal(ud_convert(15, "days", "hour"), 360)
+  expect_equal(ud_convert(-7, "ppm", "mol/mol"), -7e-06)
+})
+
+test_that("unit conversion invariants", {
+  expect_equal(ud_convert(1, "g", "g"), 1) 
+  expect_equal(ud_convert(0, "g", "kg"), 0)
+  expect_equal(ud_convert(Inf, "g", "kg"), Inf)
+})
+
+test_that("incompatible units", {
+  expect_error(ud_convert(1, "miles", "grams"))
+  expect_error(ud_convert(1, "radians", "degC"))
+  expect_error(ud_convert(1, "in", "grams"))
+})
+
+test_that("output is type numeric and not class \"units\"", {
+  x <- ud_convert(23, "degC", "K")
+  testthat::expect_failure(expect_s3_class(x, "units"))
+  testthat::expect_type(x, "double")
+
+})
diff --git a/docker/depends/pecan.depends.R b/docker/depends/pecan.depends.R
@@ -132,7 +132,7 @@ wanted <- c(
 'traits',
 'TruncatedNormal',
 'truncnorm',
-'udunits2',
+'units',
 'urltools',
 'utils',
 'vdiffr',

diff --git a/models/basgra/DESCRIPTION b/models/basgra/DESCRIPTION
@@ -13,7 +13,6 @@ Imports:
     PEcAn.utils (>= 1.4.8),
     lubridate,
     ncdf4,
-    udunits2
 Suggests:
     testthat (>= 1.0.2)
 OS_type: unix

diff --git a/models/basgra/R/run_BASGRA.R b/models/basgra/R/run_BASGRA.R
@@ -93,7 +93,7 @@ run_BASGRA <- function(run_met, run_params, site_harvest, site_fertilize, start_
 
         ## convert time to seconds
         sec <- nc$dim$time$vals
-        sec <- udunits2::ud.convert(sec, unlist(strsplit(nc$dim$time$units, " "))[1], "seconds")
+        sec <- PEcAn.utils::ud_convert(sec, unlist(strsplit(nc$dim$time$units, " "))[1], "seconds")
 
         dt <- diff(sec)[1]
         tstep <- round(86400 / dt)
@@ -121,7 +121,7 @@ run_BASGRA <- function(run_met, run_params, site_harvest, site_fertilize, start_
 
         Tair   <- ncdf4::ncvar_get(nc, "air_temperature")  ## in Kelvin
         Tair   <- Tair[ydays %in% simdays]
-        Tair_C <- udunits2::ud.convert(Tair, "K", "degC")
+        Tair_C <- PEcAn.utils::ud_convert(Tair, "K", "degC")
 
 
         #in BASGRA tmin and tmax is only used to calculate the average daily temperature, see environment.f90
@@ -360,36 +360,36 @@ run_BASGRA <- function(run_met, run_params, site_harvest, site_fertilize, start_
     outlist[[length(outlist)+1]]  <- output[thisyear, which(outputNames == "LAI")]  # LAI in (m2 m-2)
 
     CropYield     <- output[thisyear, which(outputNames == "YIELD_POT")] # (g DM m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(CropYield, "g m-2", "kg m-2")  
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(CropYield, "g m-2", "kg m-2")  
 
     clitt         <- output[thisyear, which(outputNames == "CLITT")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(clitt, "g m-2", "kg m-2")  
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(clitt, "g m-2", "kg m-2")  
 
     cstub         <- output[thisyear, which(outputNames == "CSTUB")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(cstub, "g m-2", "kg m-2")  
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(cstub, "g m-2", "kg m-2")  
 
     cst           <- output[thisyear, which(outputNames == "CST")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(cst, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(cst, "g m-2", "kg m-2") 
 
     crt           <- output[thisyear, which(outputNames == "CRT")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(crt, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(crt, "g m-2", "kg m-2") 
 
     cres          <- output[thisyear, which(outputNames == "CRES")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(cres, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(cres, "g m-2", "kg m-2") 
 
     clv           <- output[thisyear, which(outputNames == "CLV")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(clv, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(clv, "g m-2", "kg m-2") 
 
     clvd         <- output[thisyear, which(outputNames == "CLVD")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(clvd, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(clvd, "g m-2", "kg m-2") 
 
     csomf         <- output[thisyear, which(outputNames == "CSOMF")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(csomf, "g m-2", "kg m-2")  
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(csomf, "g m-2", "kg m-2")  
 
     csoms         <- output[thisyear, which(outputNames == "CSOMS")] # (g C m-2)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(csoms, "g m-2", "kg m-2")  
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(csoms, "g m-2", "kg m-2")  
 
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(csomf + csoms, "g m-2", "kg m-2") 
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(csomf + csoms, "g m-2", "kg m-2") 
 
     outlist[[length(outlist)+1]]  <- output[thisyear, which(outputNames == "TILG1")] 
     outlist[[length(outlist)+1]]  <- output[thisyear, which(outputNames == "TILG2")] 
@@ -401,22 +401,22 @@ run_BASGRA <- function(run_met, run_params, site_harvest, site_fertilize, start_
 
     # Soil Respiration in kgC/m2/s
     rsoil         <- output[thisyear, which(outputNames == "Rsoil")] # (g C m-2 d-1)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(rsoil, "g m-2", "kg m-2") / sec_in_day
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(rsoil, "g m-2", "kg m-2") / sec_in_day
 
     # Autotrophic Respiration in kgC/m2/s
     rplantaer     <- output[thisyear, which(outputNames == "RplantAer")] # (g C m-2 d-1)
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(rplantaer, "g m-2", "kg m-2") / sec_in_day
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(rplantaer, "g m-2", "kg m-2") / sec_in_day
 
     # NEE in kgC/m2/s
     # NOTE: According to BASGRA_N documentation: LUEMXQ (used in PHOT calculation) accounts for carbon lost to maintenance respiration, 
     # but not growth respiration. So, photosynthesis rate is gross photosynthesis minus maintenance respiration
     # So this is not really GPP, but it wasn't obvious to add what to get GPP, but I just want NEE for now, so it's OK
     phot          <- output[thisyear, which(outputNames == "PHOT")] # (g C m-2 d-1)
     nee           <- -1.0 * (phot - (rsoil + rplantaer))
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(nee, "g m-2", "kg m-2") / sec_in_day
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(nee, "g m-2", "kg m-2") / sec_in_day
 
     # again this is not technically GPP
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(phot, "g m-2", "kg m-2") / sec_in_day
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(phot, "g m-2", "kg m-2") / sec_in_day
 
     # Qle W/m2
     outlist[[length(outlist)+1]]  <- ( output[thisyear, which(outputNames == "EVAP")] + output[thisyear, which(outputNames == "TRAN")] * 
@@ -425,7 +425,7 @@ run_BASGRA <- function(run_met, run_params, site_harvest, site_fertilize, start_
     # SoilMoist (!!! only liquid water !!!) kg m-2
     # during the groowing season its depth will mainly be equal to the rooting depth, but during winter its depth will be ROOTD-Fdepth
     soilm <- output[thisyear, which(outputNames == "WAL")] # mm
-    outlist[[length(outlist)+1]]  <- udunits2::ud.convert(soilm, "mm", "m") * 1000 # (kg m-3) density of water in soil
+    outlist[[length(outlist)+1]]  <- PEcAn.utils::ud_convert(soilm, "mm", "m") * 1000 # (kg m-3) density of water in soil
 
     # ******************** Declare netCDF dimensions and variables ********************#
     t <- ncdf4::ncdim_def(name = "time",