Added watbal_patch, watbal_patch_mult, and documentation

NAMESPACE

@@ -34,6 +34,8 @@ export(select_parameter_sets)
 export(separate_canopy_output)
 export(tec_repeat)
 export(watbal_basin)
+export(watbal_patch)
+export(watbal_patch_mult)
 export(write_sample_clim)
 importFrom(data.table,"%between%")
 importFrom(data.table,"%like%")

R/watbal_patch.R

@@ -1,87 +1,62 @@
 #' watbal_patch.R
 #'
-#' Water balance for basin daily RHESSys output
+#' Water balance for a single patch daily RHESSys output.
+#'
+#' Before you run the water balance script below, you must have added a "rain" column to your patch daily output
+#' data that contains the patch-specific precipitation. Be careful to add the right precip values if your worldfile
+#' uses multiple base stations, isohyets, or a precip lapse rate.
+#'
+#' This function will add a number of fields to your patch output file, the last of which will be called "watbal".
+#' This is your water balance error. You should see a minor numerical error here even if your water balances (on the order of 10^-6).
+#' If your watbal values are negative then water inputs are less than water outputs, and vice versa.
+#'
+#' If you have multiple patches with only a single stratum per patch, you need to run a different script to loop
+#' through patches individually, call the water balance function, then write the output. For example, if your
+#' patch output is called "rp", your canopy output is called "rc", and you want the water balance output written
+#' to atable called "out", then you would type: "out=pwatbalmult(rp,rc)". This creates a new output table, called
+#' "out" in this example, that contains date fields and a single column for each patch with the water balance error
+#' for that patch (called "P345","P578","P900", etc. where the field label number matches the patch ID). You might
+#' see a minor numerical error here even if your water balances (on the order of 10^-6). If your watbal values are
+#' negative then water inputs are less than wateroutputs and vice versa. Note that this may take a long time to run
+#' if you have many patches; better to run on a single patch or a single hillslope than a full basin.
+#'
+#' See `watbal_patch_mult` for multiple patches.
+#'
+#' @param pd Patch daily rhessys output, read into R with a funciton like `readin_rhessys_output()`
+#' @param cd Canopy daily rhessys output, read into R with a funciton like `readin_rhessys_output()`
+#'
+#' @export
 
-# PATCH DAILY WATER BALANCE FUNCTIONS
+watbal_patch = function(pd, cd) {
 
-# NOTES:
-# (1) Before you run the water balance script below, you must have
-#     added a "rain" column to your patch daily output data that contains
-#     the patch-specific precipitation. Be careful to add the right precip
-#     values if your worldfile uses multiple base stations, isohyets,
-#     or a precip lapse rate.
-# (2) These scripts will not work if you have multiple stratums.
+  library(tidyverse)
+  pd$watbal.tmp=with(pd,rain+Qin-Qout-trans_sat-trans_unsat-evap-evap_surface-soil_evap)
+  pd$sd=with(pd,sat_def-rz_storage-unsat_stor)
 
-# TO RUN:
-# (1) In order to run this routine, you need a patch daily output file with
-#     a rain field added and a canopy daily output file.
-# (2) Copy and paste into R the water balance and multipatch run
-#     functions below (i.e., all the text below "##### PASTE BELOW #####")
-# (3) If you only have a single patch and a single canopy, you can run
-#     the watbal function as-is.  For example, if your patch output is called
-#     "rp" and your canopy output is called "rc" then type:
-#     rp=pwatbal(rp,rc)
-#     This will add a number of fields to your patch output file, the last
-#     of which will be called "watbal".  This is your water balance error.
-#     You should see a minor numerical error here even if your water balances
-#     (on the order of 10^-6).  If your watbal values are negative then
-#     water inputs are less than water outputs, and vice versa.
-# (4) If you have multiple patches with only a single stratum per patch,
-#     you need to run a different script to loop through patches
-#     individually, call the water balance function, then write the output.
-#     For example, if your patch output is called "rp", your canopy output
-#     is called "rc", and you want the water balance output written to a
-#     table called "out", then you would type:
-#     out=pwatbalmult(rp,rc)
-#     This creates a new output table, called "out" in this example, that
-#     contains date fields and a single column for each patch with the water
-#     balance error for that patch (called "P345","P578","P900", etc. where
-#     the field label number matches the patch ID). You might see a minor
-#     numerical error here even if your water balances (on the order of 10^-6).
-#     If your watbal values are negative then water inputs are less than water
-#     outputs and vice versa.
-#     *** Note that this may take a long time to run if you have many patches;
-#     better to run on a single patch or a single hillslope than a full basin.
+  cd$weighted_snow_stored = cd$snow_stored * cd$covfrac
+  cd$weighted_rain_stored = cd$rain_stored * cd$covfrac
+  tmp = cd %>% group_by(zoneID, hillID, basinID, patchID, day, month, year) %>% summarize(can_snow_stored = sum(weighted_snow_stored), can_rain_stored = sum(weighted_rain_stored) )
 
+  pd = left_join(pd, tmp[,c("basinID","hillID","zoneID","patchID","day","month","year","can_snow_stored","can_rain_stored")])
 
-##### PASTE BELOW #####
-##### WAT BAL FUNCTION #####
-pwatbal = function(qp, qc) {
-  qp$watbal.tmp = with(qp, rain + Qin - Qout - trans_sat - trans_unsat - evap - evap_surface - soil_evap)
-  qp$sd = with(qp, sat_def - rz_storage - unsat_stor)
-  tmp = diff(qp$sd)
-  tmp = c(0, tmp)
-  qp$sddiff = tmp
-  tmp = diff(qp$snow)
-  tmp = c(0, tmp)
-  qp$snodiff = tmp
-  tmp = diff(qp$detention_store)
-  tmp = c(0, tmp)
-  qp$detdiff = tmp
-  tmp = diff(qp$litter.rain_stor)
-  tmp = c(0, tmp)
-  qp$litdiff = tmp
-  tmp = diff(qc$rain_stored + qc$snow_stored)
-  tmp = c(0, tmp)
-  qp$candiff = tmp
-  qp$watbal = with(qp, watbal.tmp + sddiff - snodiff - detdiff - litdiff)
+  pd$can_water_stored = pd$can_rain_stored+pd$can_snow_stored
 
-  return(qp)
+  tmp=diff(pd$sd)
+  tmp=c(0,tmp)
+  pd$sddiff=tmp
+  tmp=diff(pd$snow)
+  tmp=c(0,tmp)
+  pd$snodiff=tmp
+  tmp=diff(pd$detention_store)
+  tmp=c(0,tmp)
+  pd$detdiff=tmp
+  tmp=diff(pd$litter.rain_stor)
+  tmp=c(0,tmp)
+  pd$litdiff=tmp
+  tmp=diff(pd$can_water_stored)
+  tmp=c(0,tmp)
+  pd$candiff=tmp
+  pd$watbal=with(pd,watbal.tmp+sddiff-snodiff-detdiff-litdiff-candiff)
+  pd$watbal[1]=0.0
+  return(pd)
 }
-##### MULTIPATCH RUN FUNCTION #####
-# pwatbalmult = function(rp, rc) {
-#   pids = unique(rp$patchID)
-#   tmp = subset(rp, rp$patchID == pids[1])
-#   wbp = tmp[c("day", "month", "year")]
-#   wbp = mkdate(wbp)
-#   n = ncol(wbp) + 1
-#   for (i in pids) {
-#     tmpp = subset(rp, rp$patchID == i)
-#     tmpc = subset(rc, rc$patchID == i)
-#     tmpp = pwatbal(tmpp, tmpc)
-#     wbp[, n] = tmpp$watbal
-#     names(wbp)[c(n)] = paste("P", i, sep = "")
-#     n = n + 1
-#   }
-#   wbp
-# }

R/watbal_patch_mult.R

@@ -0,0 +1,43 @@
+#' watbal_patch_mult.R
+#'
+#' Water balance for multiple patches daily RHESSys output. 
+#' 
+#' Before you run the water balance script below, you must have added a "rain" column to your patch daily output
+#' data that contains the patch-specific precipitation. Be careful to add the right precip values if your worldfile
+#' uses multiple base stations, isohyets, or a precip lapse rate.
+#'
+#' This function will add a number of fields to your patch output file, the last of which will be called "watbal".
+#' This is your water balance error. You should see a minor numerical error here even if your water balances (on the order of 10^-6).
+#' If your watbal values are negative then water inputs are less than water outputs, and vice versa.
+#'
+#' If you have multiple patches with only a single stratum per patch, you need to run a different script to loop
+#' through patches individually, call the water balance function, then write the output. For example, if your
+#' patch output is called "pd", your canopy output is called "cd", and you want the water balance output written
+#' to atable called "out", then you would type: "out=pwatbalmult(pd,cd)". This creates a new output table, called
+#' "out" in this example, that contains date fields and a single column for each patch with the water balance error
+#' for that patch (called "P345","P578","P900", etc. where the field label number matches the patch ID). You might
+#' see a minor numerical error here even if your water balances (on the order of 10^-6). If your watbal values are
+#' negative then water inputs are less than wateroutputs and vice versa. Note that this may take a long time to run
+#' if you have many patches; better to run on a single patch or a single hillslope than a full basin.
+#'
+#' @param pd Patch daily rhessys output, read into R with a funciton like `readin_rhessys_output()`
+#' @param cd Canopy daily rhessys output, read into R with a funciton like `readin_rhessys_output()`
+#'
+#' @export
+
+watbal_patch_mult = function(pd, cd) {
+  pids=unique(pd$patchID)
+  tmp=subset(pd,pd$patchID==pids[1])
+  wbp=tmp[c("day","month","year")]
+  wbp=mkdate(wbp)
+  n=ncol(wbp)+1
+  for (i in pids) {
+    tmpp=subset(pd,pd$patchID==i)
+    tmpc=subset(cd,cd$patchID==i)
+    tmpp=watbal_patch(tmpp,tmpc)
+    wbp[,n]=tmpp$watbal
+    names(wbp)[c(n)]=paste("P",i,sep="")
+    n=n+1
+  }
+  return(wbp)
+}