library(chron); library(gtools) tankstats <- function(dailyrunoff, format = "d/m/y", #data.frame with daily runoff data in mm (runoff.mm.d) and date as format tankvol, tankbegin = tankvol/2, #tank volume in L: volume of water in L at start of run carea, #tank's catchment area in square metres propcon = 1, #proportion of carea that is connected to stormwater initloss = 0, #initial loss in mm (accounted for in dailyrunoff) firstflush = carea*0.2,#first flush bypass volume in L) npeople = 2, #number of people in house toilet = 18.9,#average toilet usage per person per day in L other = rep(0,12), #made up values for car washing wmac1 = 35.31, wmacadd = 23.54, #average washing machine use per day in L for 1 person, #and additional use for each additional person hotwater = 46.9, # average daily hot water use per person in L #(Wilkenfield's estimate of 61.9 minus half the washing machine use outdoor.annual = 12971, #annual outdoor irrigation use in L per 100 sq m of garden garden.area = 350, #area of garden to be irrigated from the tank in sq m outdoor.dist = c(.27,.21,.09,.07,.05,0,0,0,.03,.04,.04,.2),...){ #monthly distribution of garden watering options(chron.origin = c(month = 1, day = 1, year = 1960)) dailyrunoff$date <- chron(as.vector(dailyrunoff$date),format = format) dailyrunoff$month <- month.day.year(dailyrunoff$date)$month ndays <- length(dailyrunoff$date) nyears <- round(ndays/365.25,0) #distribute the number of car wash days into random Saturdays # if(tankvol == 0){ # cat("Note: tank volume is zero. # These results model the effect of a first-flush diverter by itself.") # } # if(as.numeric(carfreq) > 26) # stop("Our model only allows for a maximum of 1 car wash per fortnight (maximum of 26 washes per year)", # call. = FALSE) # saturdays <- which(weekdays(dailyrunoff$date) == "Sat") # carwashdays <- saturdays[order(runif(saturdays))][1:(carfreq*nyears)] daily.outdoor <- outdoor.dist[match(dailyrunoff$month,1:12)] * outdoor.annual * garden.area * 12 /36500 daily.other <- other[match(dailyrunoff$month, 1:12)]/30.4 # if(sum(!is.na(carwashdays))> 0 ) # daily.carwash <- rep(0,ndays) # if(carfreq > 0) # daily.carwash[carwashdays] <- daily.carwash[carwashdays] + carwash annual.outdoor <- sum(daily.outdoor)/(1000*nyears) #I have decided not to include car washing as an outdoor use because it is a small use #and because this water (+soap) is more likely to find its way into the stormwater system #and the purpose of quantifying this is to assume that outdoor use results in loss of N #in unsewered properties if(npeople < 1) stop("Number of people (npeople) must be at least 1 \n", call. = FALSE) toilet <- npeople*toilet hotwater <- npeople*hotwater wmac <- wmac1 + (npeople - 1)*wmacadd x <- cbind(combinations(5, 1, 1:5), array(0, dim = c(5, 5 - 1))) for (i in 2:5) { nc <- dim(combinations(5, i, 1:5))[1] x <- rbind(x, cbind(combinations(5, i, 1:5), array(0, dim = c(nc, 5 - i)))) } usage.desc <- c("garden", "other", "toilet", "laundry", "hot water") cat(usage.desc[x[1,]], file = "temp.txt", sep = ", ") options(warn = -1) usage.descs <- read.table("temp.txt", sep = "#") options(warn = 0) consumption <- daily.outdoor for(i in 2:dim(x)[1]){ uses <- t(array(rep(!is.na( match(1:5,x[i,])),ndays),dim = c(5,ndays))) consumption <- cbind(consumption, (daily.outdoor*uses[,1] + daily.other*uses[,2] + toilet*uses[,3] + wmac*uses[,4] + hotwater*uses[,5])) cat(usage.desc[x[i,]], file = "temp.txt", sep = ", ") options(warn = -1) usage.descs <- rbind(usage.descs, read.table("temp.txt", sep = "#")) options(warn = 0) } usage.descs <- as.vector(usage.descs[,1]) n.uses <- length(usage.descs) tankbegin <- ifelse(length(tankbegin) == 1, rep(tankbegin,n.uses), tankbegin) annual.consumption <- apply(consumption,2,sum)/(1000*nyears) # in kL - #is this correction right?* 12 /36500 annual.outdoor.consumption <- rep(0,31) annual.outdoor.consumption[grep("garden", usage.descs)] <- round(annual.outdoor,0) # if(tankvol < firstflush) # stop("Tank volume is too small", call. = FALSE) #I'm allowing this so that simple first flush diverters are permitted inflow <- apply(cbind(carea*(dailyrunoff$runoff.mm.d - initloss) - firstflush, rep(0,dim(dailyrunoff)[1])),1,max) use <- apply(rbind(consumption[1,], tankbegin), 2, FUN = min) overflow <- apply(rbind(rep(0,n.uses), tankbegin + inflow[1] - use - tankvol), 2,FUN = max) store <- apply(rbind(rep(0,n.uses),tankbegin + inflow[1] - use - overflow), 2, FUN = max) budget <- list(store = store, overflow = overflow, use = use) for(i in 2:dim(dailyrunoff)[1]){ if(i > 2) {store <- budget$store[i-1,]} budget$use <- rbind(budget$use, apply(rbind(consumption[i,], store), 2,FUN = min)) budget$overflow <- rbind(budget$overflow, apply(rbind(rep(0,n.uses), store + inflow[i] - budget$use[i,] - tankvol), 2, FUN = max)) budget$store <- rbind(budget$store, apply(rbind(rep(0,n.uses), store + inflow[i] - budget$use[i,] - budget$overflow[i,]), 2, FUN = max)) } annual.yield.kL <- round(apply(budget$use,2,sum)* 0.36525/dim(dailyrunoff)[1],1) days.runoff <- round(apply(budget$overflow > 0,2,sum)* 365.25/dim(dailyrunoff)[1],0) nat.rof <- 12 #was 14 in original set of data #round(sum(runoffcomplete$runoff.mm.d > 15)*365.25/dim(dailyrunoff)[1],0) #rof = runoff frequency: assume surface runoff occurred from #pre-urban catchment with >15 mm runoff urban.rof <- 121 #was 114 in original set of data #round(sum(runoffcomplete$runoff.mm.d >1)*365.25/dim(dailyrunoff)[1],0) #assume surface runoff occurrs from urban catchment with >1 mm runoff # assume that unconnected surfaces have no existing impact on streams or on N export, so discount #these areas in calculation of stream and N indices (hence *propcon) rof.index <- round(1 - apply(rbind(rep(0,n.uses), (days.runoff - nat.rof)/(urban.rof - nat.rof)),2,max), 2) * carea*propcon/100 urban.rov <- round(sum(dailyrunoff$runoff.mm[dailyrunoff$runoff.mm.d > initloss])* #was > 1 carea*0.36525/dim(dailyrunoff)[1], 0) #annual harvestable volume from carea rov.index <- round(annual.yield.kL/urban.rov,2)*carea/100 #assume firstflush volume is also kept out of the stormwater system N.index <- round((annual.yield.kL + (firstflush*121/1000))/urban.rov,2)*carea*propcon/100 calc.EBI <- function(x){weighted.mean(x,w = c(0.5,0.3,0.2))} list(budget = budget, consumption = consumption, EBstats = data.frame(usage = usage.descs, consumption.kL.y = round(annual.consumption,0), outdoor.consumption = annual.outdoor.consumption, yield.kL.y = annual.yield.kL, days.runoff = days.runoff, stream.index = rof.index, water.index = rov.index, N.index = N.index, EBI = round(apply(cbind(rof.index, propcon*rov.index,rov.index),1,FUN = calc.EBI),2)), uses.matrix = x, parameters = c(tankvol, carea, firstflush, npeople, garden.area)) } #Here's a simple error bar function in R. # produce vertical error bars on a plot-- aliased vbars in anticipation # of the day when I need hbars to plot horizontal error bars.... # # error bars originate at (x, y0) and radiate to y0 + y1 and y0 - y2 # # example (adds +- 1.0 SE to a barplot of means for data vector x): # # y.mean <- mean(x) # mean # y.se <- sqrt(var(x,na.rm=TRUE)/length(x)) # std err # z <- barplot(y.mean,plot=FALSE) # location of bars # barplot(y.mean,...) # the real plot # ebars(z,y.mean,y.mean+y.se,y.mean-y.se) # add error bars ebars <- vbars <- function (x, y0, y1, y2,...) { for(i in 1:length(x)) { if(y0[i]!=0 & !is.na(y0[i]) & y1[i]!=0 & !is.na(y1[i]) & (y0[i]!=y2[i])) { arrows (x[i], y0[i], x[i], y1[i], angle=90, length=0.05,...) } if(y0[i]!=0 & !is.na(y0[i]) & y2[i]!=0 & !is.na(y2[i]) & (y0[i]!=y2[i])) { arrows (x[i], y0[i], x[i], y2[i], angle=90, length=0.05,...) } } }