Standardize VPD treatment between CRU, watch_wfdei, and wfde5

All three cases now compute 'vapr' in Pa within 'ingest_globalfields()'. Previously for CRU there was a variable 'vap' in hPa.
geco-bern · Dec 20, 2024 · 764eb68 · 764eb68
1 parent 6fd1629
commit 764eb68
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 59 deletions.
diff --git a/R/ingest.R b/R/ingest.R
@@ -64,7 +64,7 @@ ingest <- function(
   patm_base <-patm_mean <- month <- tavg <-temp <- temp_fine <-
   tmax <- tmax_fine <- tmin <- tmin_fine <- prec <- prec_fine <-
   days_in_month <- rain <- snow <- srad <- srad_fine <- ppfd <-
-  ppfd_fine <- wind <- wind_fine <- qair <- vap <- vapr <- vapr_fine <-
+  ppfd_fine <- wind <- wind_fine <- qair <- vapr <- vapr_fine <-
   ilon <- data <- yy <- mm <- co2_avg <- year <- . <- bias <- NULL
 
   # Check: all sites are distinct wrt name, lon and lat
@@ -563,25 +563,15 @@ ingest <- function(
         }
 
         # Bias correction for relative humidity (actually vapour pressure): scale
-        if ("vapr" %in% getvars_wc){
+        if ("vapr" %in% getvars_wc){ # i.e. equivalent with: "if (vpd" %in% getvars)"
 
-          # calculate vapour pressure from specific humidity - needed for bias correction with worldclim data
+          # a) correct vapor pressure (vapr, Pa) with worldclim data
           if (source == "watch_wfdei" || source == "wfde5"){
-            # specific humidity (qair, g g-1) is read, convert to vapour pressure (vapr, Pa)
-            ddf <- ddf %>% 
-              rowwise() %>% 
-              dplyr::mutate(
-                vapr = calc_vp(qair = qair,
-                               patm = patm)
-                ) %>% 
-              ungroup()
-
+            # For watch_wfdei, wfde5, vapr has already been computed based on qair and patm 
+            # within `ingest_globalfields()`. Since there is no bias-correction of qair and ptm
+            # it does not need to be recomputed here
           } else if (source == "cru"){
-            # vapour pressure is read from file, convert from hPa to Pa
-            ddf <- ddf %>% 
-              dplyr::mutate(vapr = 1e2 * vap) %>% 
-              dplyr::select(-vap)
-
+            # For cru, vapr has already been unit-transformed based on vap within `ingest_globalfields()`
           }
 
           df_bias <- df_fine %>%
@@ -612,11 +602,8 @@ ingest <- function(
             mutate(vapr = ifelse(is.na(scale), vapr, vapr * scale)) %>%
             dplyr::select(-scale, -month)
 
-        }
-
-        # Calculate vapour pressure deficit from specific humidity
-        if ("vpd" %in% getvars){
-
+          # b) re-calculate vapour pressure deficit VPD from bias-corrected vapor pressure
+          #    (Note: this overwrites the vpd that has been computed within `ingest_globalfields()`)
           if (source == "watch_wfdei" || source == "wfde5"){
             # use daily mean temperature
             ddf <- ddf %>%
@@ -625,7 +612,7 @@ ingest <- function(
                 vapr = calc_vp(
                   qair = qair,
                   patm = patm
-                  ),
+                ),
                 vpd = calc_vpd(eact = vapr, tc = temp)) %>% 
               ungroup()
 
@@ -635,7 +622,7 @@ ingest <- function(
               rowwise() %>%
               dplyr::mutate(
                 vpd = calc_vpd(eact = vapr, tmin = tmin, tmax = tmax)
-                ) %>% 
+              ) %>% 
               ungroup()
           }
         }
@@ -646,12 +633,6 @@ ingest <- function(
 
       }
 
-    } else {
-      if ("vpd" %in% getvars){
-        # For cases where there is no bias correction,
-        # (of sources cru, watch_wfdei, wfde5; but not ndep)
-        # vapour pressure deficit has already been computed within `ingest_globalfields()`
-      }
     }
 
 	} else if (source == "gee"){

diff --git a/R/ingest_bysite.R b/R/ingest_bysite.R
@@ -63,7 +63,7 @@ ingest_bysite <- function(
     patm_base <-patm_mean <- month <- tavg <- temp <- temp_fine <-
     tmax <- tmax_fine <- tmin <- tmin_fine <- prec <- prec_fine <-
     days_in_month <- rain <- snow <- srad <- srad_fine <- ppfd <-
-    ppfd_fine <- wind <- wind_fine <- qair <- vap <- vapr <- vapr_fine <-
+    ppfd_fine <- wind <- wind_fine <- qair <- vapr <- vapr_fine <-
     ilon <- data <- yy <- mm <- co2_avg <- year <- . <- bias <-
     co2 <- lon...1 <- lat...2 <- bottom <- top <- depth <- var <-
     var_wgt <- depth_tot_cm <- NULL
@@ -376,22 +376,15 @@ ingest_bysite <- function(
         }
 
         # Bias correction for relative humidity (actually vapour pressure): scale
-        if ("vapr" %in% getvars_wc){
+        if ("vapr" %in% getvars_wc){ # i.e. equivalent with: "if (vpd" %in% getvars)"
 
-          # calculate vapour pressure from specific humidity - needed for bias correction with worldclim data
-          if (source == "watch_wfdei"){
-            # specific humidity (qair, g g-1) is read, convert to vapour pressure (vapr, Pa)
-            df_tmp <- df_tmp %>% 
-              rowwise() %>% 
-              dplyr::mutate(vapr = calc_vp(qair = qair, patm = patm)) %>% 
-              ungroup()
-
+          # a) correct vapor pressure (vapr, Pa) with worldclim data
+          if (source == "watch_wfdei" || source == "wfde5"){
+            # For watch_wfdei, wfde5, vapr has already been computed based on qair and patm 
+            # within `ingest_globalfields()`. Since there is no bias-correction of qair and ptm
+            # it does not need to be recomputed here
           } else if (source == "cru"){
-            # vapour pressure is read from file, convert from hPa to Pa
-            df_tmp <- df_tmp %>% 
-              dplyr::mutate(vapr = 1e2 * vap) %>% 
-              dplyr::select(-vap)
-
+            # For cru, vapr has already been unit-transformed based on vap within `ingest_globalfields()`
           }
 
           df_bias <- df_fine %>% 
@@ -418,12 +411,9 @@ ingest_bysite <- function(
             left_join(df_bias %>% dplyr::select(month, scale), by = "month") %>% 
             mutate(vapr = ifelse(is.na(scale), vapr, vapr * scale)) %>% 
             dplyr::select(-scale, -month)
-        }      
-
-
-        # Calculate vapour pressure deficit from specific humidity
-        if ("vpd" %in% getvars){
 
+          # b) re-calculate vapour pressure deficit VPD from bias-corrected vapor pressure
+          #    (Note: this overwrites the vpd that has been computed within `ingest_globalfields()`)
           if (source == "watch_wfdei" || source == "wfde5"){
             # use daily mean temperature
             df_tmp <- df_tmp %>%
@@ -438,21 +428,14 @@ ingest_bysite <- function(
               dplyr::mutate(vpd = calc_vpd(eact = vapr, tmin = tmin, tmax = tmax)) %>% 
               ungroup()
           }
-
-        }
+        }      
 
         # keep only required dates
         df_tmp <- df_tmp %>% 
           right_join(ddf_dates, by = c("sitename", "date"))
 
       }
 
-    } else {
-      if ("vpd" %in% getvars){
-        # For cases where there is no bias correction,
-        # (of sources cru, watch_wfdei, wfde5; but not ndep)
-        # vapour pressure deficit has already been computed within `ingest_globalfields()`
-      }
     }
 
   } else if (source == "modis"){

diff --git a/R/ingest_globalfields.R b/R/ingest_globalfields.R
@@ -399,7 +399,8 @@ ingest_globalfields <- function(
     if ("vpd" %in% getvars){
       df_out <- df_out %>% 
         rowwise() %>%
-        mutate(vpd = calc_vpd( eact = 1e2 * vap, tmin = tmin, tmax = tmax )) %>%
+        mutate(vapr = 1e2 * vap) %>% # to go from hPa to Pa
+        mutate(vpd  = calc_vpd( eact = vapr, tmin = tmin, tmax = tmax )) %>%
         ungroup() # undo rowwise()
     }