#' readin_rhessys_output
#'
#' Pulls output from a single RHESSys run into R
#' @param pre file prefix. Does not include spatial or temporal level of aggrigation
#' ie. "my_watershed" NOT "my_watershed_basin" or "my_watershed_basin.daily"
#' @param b basin
#' @param c canopy
#' @param z basin
#' @param p canopy
#' @param h basin
#' @param g canopy
#' @param f fire
#' @param wy wateryear
#' @param dyn_read dynamically read all the output that exists
#'
#' @export
readin_rhessys_output_old = function(pre, b = 1, c = 0, z = 0, p = 0, h = 0, g = 0, f = 0, wy = 1, dyn_read = FALSE, clean = TRUE) {
# Check for more than 1 line of data in all the output files that meet the prefix
if (dyn_read) {
suf = c("basin.daily", "grow_basin.daily", "hillslope.daily", "grow_hillslope.daily", "zone.daily", "grow_zone.daily",
"patch.daily", "grow_patch.daily", "stratum.daily", "grow_stratum.daily", "fire.daily")
files = paste(pre,"_",suf,sep="")
fexist = file.exists(files)
# check if more than 1 line and set to false if not
fexist[fexist] = sapply(files[fexist],function(x) length(readLines(x,n = 3))) > 1
# set flags based on if files exist and are > 1 line
if (fexist[2] | fexist[4] | fexist[6] | fexist[8] | fexist[10]) {g = 1} else {g = 0}
if (fexist[1]) {b = 1} else {b = 0}
if (fexist[3]) {h = 1} else {h = 0}
if (fexist[5]) {z = 1} else {z = 0}
if (fexist[7]) {p = 1} else {p = 0}
if (fexist[9]) {c = 1} else {c = 0}
if (fexist[11]) {f = 1} else {f = 0}
}
bd = NULL
bdg = NULL
cd = NULL
cdg = NULL
zd = NULL
zdg = NULL
pd = NULL
pdg = NULL
fd = NULL
bd.wy = NULL
bdg.wy = NULL
cd.wy = NULL
cdg.wy = NULL
zd.wy = NULL
zdg.wy = NULL
pd.wy = NULL
pdg.wy = NULL
fd.wy = NULL
bd.wyd = NULL
bdg.wyd = NULL
cd.wyd = NULL
cdg.wyd = NULL
zd.wyd = NULL
zdg.wyd = NULL
pd.wyd = NULL
pdg.wyd = NULL
fd.wyd = NULL
hd = NULL
hdg = NULL
hd.wy = NULL
hdg.wy = NULL
hdg.wyd = NULL
hd.wyd = NULL
if (b == 1) {
nme = sprintf("%s_basin.daily", pre)
bd = read.table(nme, header = T)
bd$et = bd$trans + bd$evap
bd$unfilled_cap = bd$sat_def - bd$unsat_stor - bd$rz_storage
bd = mkdate(bd)
if (wy == 1) {
bd.wy = aggregate(bd, by = list(bd$wy), mean)
bd.wyd = aggregate(bd, by = list(bd$wyd), mean)
tmp = aggregate(bd$snowpack, by = list(bd$wy), max)
bd.wy$pkswe = tmp$x
}
if (g == 1) {
nme = sprintf("%s_grow_basin.daily", pre)
bdg = read.table(nme, header = T)
bdg = mkdate(bdg)
porosity = bd$sat_def/bd$sat_def_z
if (nrow(bd)==nrow(bdg)){
bd$veg_awc = ifelse((bdg$root_depth < bd$sat_def_z),
porosity * (bd$sat_def_z - bdg$root_depth) +
bd$rz_storage, bd$rz_storage)
} else {
print("bd and bdg have a mismatch number of rows. veg_awc not computed.")
}
if (wy == 1) {
bdg.wy = aggregate(bdg, by = list(bdg$wy), mean)
bdg.wyd = aggregate(bdg, by = list(bdg$wyd),
mean)
}
}
}
if (z == 1) {
nme = sprintf("%s_zone.daily", pre)
zd = read.table(nme, header = T)
zd = mkdate(zd)
if (wy == 1) {
zd.wy = aggregate(zd, by = list(zd$wy), mean)
zd.wyd = aggregate(zd, by = list(zd$wyd), mean)
}
if (g == 1) {
nme = sprintf("%s_grow_zone.daily", pre)
zdg = read.table(nme, header = T)
zdg = mkdate(zdg)
if (wy == 1) {
zdg.wy = aggregate(zdg, by = list(zdg$wy), mean)
zdg.wyd = aggregate(zdg, by = list(zdg$wyd),
mean)
}
}
}
if (h == 1) {
nme = sprintf("%s_hillslope.daily", pre)
hd = read.table(nme, header = T)
hd = mkdate(hd)
if (wy == 1) {
hd.wy = aggregate(hd, by = list(hd$wy), mean)
hd.wyd = aggregate(hd, by = list(hd$wyd), mean)
}
if (g == 1) {
nme = sprintf("%s_grow_hillslope.daily", pre)
hdg = read.table(nme, header = T)
hdg = mkdate(hdg)
if (wy == 1) {
hdg.wy = aggregate(hdg, by = list(hdg$wy), mean)
hdg.wyd = aggregate(hdg, by = list(hdg$wyd),
mean)
}
}
}
if (c == 1) {
nme = sprintf("%s_stratum.daily", pre)
cd = read.table(nme, header = T)
nme = sprintf("%s_stratum.daily", pre)
cd = mkdate(cd)
if (wy == 1) {
cd.wy = aggregate(cd, by = list(cd$wy), mean)
cd.wyd = aggregate(cd, by = list(cd$wyd), mean)
}
if (g == 1) {
nme = sprintf("%s_grow_stratum.daily", pre)
cdg = read.table(nme, header = T)
cdg = mkdate(cdg)
cdg$woodc = cdg$live_stemc + cdg$dead_stemc + cdg$live_crootc +
cdg$dead_crootc
cdg$plantc = cdg$woodc + cdg$frootc + cdg$leafc
if (wy == 1) {
cdg.wy = aggregate(cdg, by = list(cdg$wy), mean)
cdg.wyd = aggregate(cdg, by = list(cdg$wyd),
mean)
tmp = aggregate(cdg$cpool, by = list(cdg$wy),
min)
cdg.wy$mincpool = tmp$x
cdg.wy$mincpoolp = tmp$x/cdg.wy$plantc * 100
}
}
}
if (p == 1) {
nme = sprintf("%s_patch.daily", pre)
pd = read.table(nme, header = T)
pd = mkdate(pd)
if (wy == 1) {
pd.wy = aggregate(pd, by = list(pd$wy), mean)
pd.wyd = aggregate(pd, by = list(pd$wyd), mean)
}
if (g == 1) {
nme = sprintf("%s_grow_patch.daily", pre)
pdg = read.table(nme, header = T)
pdg = mkdate(pdg)
if (wy == 1) {
pdg.wy = aggregate(pdg, by = list(pdg$wy), mean)
pdg.wyd = aggregate(pdg, by = list(pdg$wyd),
mean)
}
}
}
if (f == 1) {
nme = sprintf("%s_fire.daily", pre)
fd = read.table(nme, header = T)
fd = mkdate(fd)
if (wy == 1) {
fd.wy = aggregate(fd, by = list(fd$wy), mean)
fd.wyd = aggregate(fd, by = list(fd$wyd), mean)
}
}
a = list(bd = bd, bdg = bdg, bd.wy = bd.wy, bdg.wy = bdg.wy, bd.wyd = bd.wyd,
bdg.wyd = bdg.wyd, fd = fd, fd.wy = fd.wy, fd.wyd = fd.wyd,
hd = hd, hdg = hdg, hd.wy = hd.wy, hdg.wy = hdg.wy, hd.wyd = hd.wyd,
hdg.wyd = hdg.wyd, zd = zd, zdg = zdg, zd.wy = zd.wy, zdg.wy = zdg.wy,
zd.wyd = zd.wyd, zdg.wyd = zdg.wyd, pd = pd, pdg = pdg, pd.wy = pd.wy,
pdg.wy = pdg.wy, pd.wyd = pd.wyd, pdg.wyd = pdg.wyd, cd = cd,
cdg = cdg, cd.wy = cd.wy, cdg.wy = cdg.wy, cd.wyd = cd.wyd,
cdg.wyd = cdg.wyd)
# get rid of empty lists
if (clean) {
b = !sapply(a,is.null)
a = a[b]
}
return(a)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.