Consolidate methods and standardize units for calculating saturated vapor pressure deficit

CHANGELOG.md

@@ -11,7 +11,8 @@ section for the next release.
 
 * Add function `clip_and_save_raster_file()` for subsetting rasters to match a polygon of interest (#3537).
 * Add CH4 and N2O to standard_vars in PEcAn.utils
-
+* New function `sat_vapor_pressure()` added for computing saturation vapor pressure from temperature using various methods.
+* 
 ## [1.9.0] - 2025-05-25
 
 ### Added

modules/data.atmosphere/NEWS.md

@@ -2,7 +2,7 @@
 
 * `ERA5_met_process()` can now process ensemble data efficiently in parallel using new option `n_cores`
 * Dependency `ggplot2` is now suggested rather than required. It is used in two vignettes and for optional diagnostic plots from `debias_met_regression`.
-
+* New function `sat_vapor_pressure()` added for computing saturation vapor pressure from temperature using various methods.
 
 # PEcAn.data.atmosphere 1.9.0
 

modules/data.atmosphere/R/metutils.R

@@ -69,31 +69,104 @@ get.vpd <- function(rh, temp) {
   ## calculate saturation vapor pressure
   es <- get.es(temp)
   ## calculate vapor pressure deficit
-  return(((100 - rh)/100) * es)
+  return(((100 - rh) / 100) * es)
 } # get.vpd
 
-##' Calculate saturation vapor pressure
-##'
-##' @title get es
-##' @param temp temperature in degrees C
-##' @return saturation vapor pressure in mb
-##' @export
-##' @author David LeBauer
-##' @examples
-##' temp <- -30:30
-##' plot(temp, get.es(temp))
+#' Saturation vapor pressure
+#'
+#' @md
+#' Compute saturation vapor pressure from temperature using one of the
+#' following methods:
+#' - "Magnus" — Alduchov & Eskridge (1996)
+#' - "ClausiusClapeyron" — FAO-56-style with constant L
+#' - "GoffGratch" — Goff & Gratch (1946; over liquid)
+#'
+#' Each method uses different units internally, users can specify units
+#' for both inputs and outputs, with defaults "degC" and "kPa", respectively.
+#' @param temp numeric vector of temperatures
+#' @param method c("Magnus","ClausiusClapeyron","GoffGratch")
+#' @param temp_units input temperature units ("degC","K","degF"), default "degC"
+#' @param out_units output pressure units ("kPa","hPa","Pa","mb"), default "kPa"
+#' @return numeric vector in `out_units`
+#'
+#' @references
+#' Alduchov & Eskridge (1996) J. Appl. Meteor. 35:601–609.
+#' Allen et al. (1998) FAO-56.
+#' Goff & Gratch (1946) Trans. ASHVE 52:95–122.
+#'
+#' @author David LeBauer
+#' @export
+sat_vapor_pressure <- function(
+    temp,
+    method = c("Magnus", "ClausiusClapeyron", "GoffGratch"),
+    temp_units = "degC",
+    out_units = "kPa") {
+  method <- match.arg(method)
+  # normalize common alias
+  if (tolower(out_units) == "mb") out_units <- "hPa"
+
+  if (method == "Magnus") {
+    # canonical temp: degC; canonical pressure: kPa
+    Tc <- units::ud_convert(temp, temp_units, "degC")
+    es_kPa <- 0.61078 * exp((17.27 * Tc) / (Tc + 237.3))
+    return(units::ud_convert(es_kPa, "kPa", out_units))
+  }
+
+  if (method == "ClausiusClapeyron") {
+    # canonical temp: degC; canonical pressure: hPa
+    Tc <- units::ud_convert(temp, temp_units, "degC")
+    L <- 2.5e6 # J kg^-1
+    Rv <- 461 # J kg^-1 K^-1
+    es_hPa <- 6.11 * exp((L / Rv) * (1 / 273 - 1 / (273 + Tc)))
+    return(units::ud_convert(es_hPa, "hPa", out_units))
+  }
+
+  if (method == "GoffGratch") {
+    # canonical temp: K; canonical pressure: hPa
+    Tk <- units::ud_convert(temp, temp_units, "K")
+    Tst <- 373.15 # K
+    est <- 1013.246 # hPa at steam point
+    lg10 <- function(z) log10(z)
+    log10_es <- -7.90298 * (Tst / Tk - 1) +
+      5.02808 * lg10(Tst / Tk) -
+      1.3816e-7 * (10^(11.344 * (1 - Tk / Tst)) - 1) +
+      8.1328e-3 * (10^(-3.49149 * (Tst / Tk - 1)) - 1) +
+      lg10(est)
+    es_hPa <- 10^log10_es
+    return(units::ud_convert(es_hPa, "hPa", out_units))
+  }
+
+  PEcAn.logger::logger.severe(
+    method,
+    "for converting temperature to saturated vapor pressure not supported"
+  )
+}
+
+# ---- Aliases for backward-compatibility ----
+
+#' @rdname sat_vapor_pressure
+#' @md
+#' @export
 get.es <- function(temp) {
-  return(6.11 * exp((2500000/461) * (1/273 - 1/(273 + temp))))
-} # get.es
+  sat_vapor_pressure(
+    temp = temp,
+    method = "ClausiusClapeyron",
+    temp_units = "degC",
+    out_units = "hPa"
+  )
+}
 
-## TODO: merge SatVapPress with get.es; add option to choose method
+#' @rdname sat_vapor_pressure
+#' @md
+#' @export
 SatVapPres <- function(T) {
-  # /estimates saturation vapor pressure (kPa) Goff-Gratch 1946 /input: T = absolute temperature
-  T_st <- 373.15  ##steam temperature (K)
-  e_st <- 1013.25  ##/saturation vapor pressure at steam temp (hPa)
-  return(0.1 * exp(-7.90298 * (T_st/T - 1) + 5.02808 * log(T_st/T) - 1.3816e-07 * (10^(11.344 * (1 - T/T_st)) -
-    1) + 0.0081328 * (10^(-3.49149 * (T_st/T - 1)) - 1) + log(e_st)))
-} # SatVapPres
+  sat_vapor_pressure(
+    temp = T,
+    method = "GoffGratch",
+    temp_units = "K",
+    out_units = "kPa"
+  )
+}
 
 
 ##' Calculate RH from temperature and dewpoint

modules/data.atmosphere/tests/testthat/test.metutils.R

@@ -31,4 +31,52 @@ test_that("get.rh RH from dewpoint",{
   expect_equal(getrhtest(25, 10), 38.82, tolerance = 0.2)
   expect_equal(getrhtest(0, -5), 69, tolerance = 0.2)
 })
+
+test_that("different methods of sat_vapor_pressure work correctly", {
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "Magnus"),
+    c(0.286, 1.228),
+    tolerance = 0.001
+  )
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "ClausiusClapeyron"),
+    c(0.287, 1.233), 
+    tolerance = 0.001
+  )
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "GoffGratch"),
+    c(0.286, 1.227),
+    tolerance = 0.001
+  )
+})
+
+test_that("sat_vapor_pressure works with different units", {
+  expect_equal(
+    sat_vapor_pressure(283.15,
+      method = "GoffGratch",
+      temp_units = "K",
+      out_units = "mb"
+    ),
+    12.27,
+    tolerance = 0.01
+  )
 
+    expect_equal(
+      sat_vapor_pressure(283.15,
+        method = "ClausiusClapeyron",
+        temp_units = "K",
+        out_units = "kPa"
+      ),
+      1.227,
+      tolerance = 0.01
+    ),
+    expect_equal(
+      sat_vapor_pressure(283.15,
+        method = "Magnus",
+        temp_units = "K",
+        out_units = "Pa"
+      ),
+      1227,
+      tolerance = 1
+    )
+})