library("timeSeries")
library("reshape")

# http://www.usbr.gov/lc/region/g4000/weekly.pdf
# total capacity ~ 59,710,000 acre-feet
# mead ~ 26,120,000 acre-feet (~43% total)
# powell ~ 24,270,000 acre-feet (~41% total)
# mead + powell ~ 84% system capacity
# mohave ~ 1,815,000 acre-feet
# havasu ~  617,000 acre-feet

# see: LM_AreaCapacityTables2009
# http://www.usbr.gov/lc/region/g4000/LM_AreaCapacityTables2009.pdf
f_mead_elev_vol = "data/mead-elev-volume.txt"
mev_vol = read.table(f_mead_elev_vol,na.strings="9999")

# http://www.usbr.gov/lc/region/g4000/hourly/mead-elv.html
f_mead_elev = "data/mead-elev.txt"
mev = read.table(f_mead_elev, na.strings="9999.99", header=T)
colnames(mev)=c("Year","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

# convert mead elevation data to volume data
mev=melt(mev, id=c("Year")) # matrix to column
mev[,2] = mev[,1]+match(mev[,2],month.abb)/12 - 1/24 # digital yyyy.mo dates
mev = mev[order(mev[,2]),] # reorder the column data by date
elev_idx = round(mev[,3]) - 894 # get resevoir height from base elev
elev_idx[1:5]=NA # the very early dates have negative heights
mvol = cbind(mev[,2],mev_vol[elev_idx,2]) # use list of dated elevations to get matching volumes
mvol = mvol[!is.na(mvol[,2]),]

# http://www.usbr.gov/uc/crsp/GetSiteInfo -> Lake Powell
# default dates, view by selected date -> download data set file -> csv
# http://www.usbr.gov/uc/crsp/GetDateInfo?d0=1719&d1=1792&d2=1862&d3=1872&d4=1928&idCount=5&l=LAKE+POWELL
# http://www.usbr.gov/uc/crsp/GetDataSet?l=LAKE+POWELL&c=1719&strSDate=11-MAR-1963&strEDate=11-JUL-2014
f_pow = "data/acre-feet-powell-20140711.lst"
pvol = read.table(f_pow, na.strings="9999", header=T,sep=",",col.names = c("date","vol"))
pvol = pvol[!is.na(pvol[,2]),] # remove NA
pvol = pvol[4:nrow(pvol),] # remove the tail end of June 1963

# use timeSeries as route to monthy mean data
# beginning date, but this does skew date reported for monthly
pvol[,1] = as.Date(as.numeric(time(pvol[,1])), origin="1963-07-01") 
powell = timeSeries(pvol[,2],pvol[,1])
by = timeSequence(from=start(powell), to = end(powell), by="month")
pvol = aggregate(powell,by,mean)

# 343 is the index which matches the start of Lake Powell data (Jul 63)
# 338 is the index which matches the start of Lake Powell data (Jul 63)
idx = min(length(pvol),nrow(mvol)-338+1)
#lmead = rep(NA,nrow(pvol))
lmead = mvol[338:(338+idx-1),]
lpowell = pvol[1:idx]

# strip na from tail
#naidx = min(sum(!is.na(lpowell)),sum(!is.na(lmead)))
#lpowell = lpowell[1:naidx]
#lmead = lmead[1:naidx]

# get time line
y = as.integer(substr(as.character(time(pvol)),1,4))
yr_idx=length(y)
yrs = rep(NA,length(lpowell))
yrs[1:length(y)] = y
xtcks = y[seq(19,yr_idx,by=60)]

# mead ~ 26,120,000 acre-feet (~43% total)
# powell ~ 24,270,000 acre-feet (~41% total)
sum = rowSums(cbind(lpowell,lmead[,2]),na.rm=T)
subtot = 26.12e6+24.27e6

def.par = par(no.readonly = TRUE) 
graphic_output = TRUE
par(mar=c(5.1,4.1,4.1,3))
if (graphic_output) png("img/mead-powell-volume.png",w=600,h=400)
	plot(100*sum/subtot,type="l",ylab="% combined storage (50 million acre-feet)",xlab="",xaxt="n",xaxs="i",yaxs="i",ylim=c(0,100))
	axis(1, at=seq(7,yr_idx,by=12),labels=rep(" ",length(seq(7,yr_idx,by=12))),las=2)
	axis(1, at=seq(19,yr_idx,by=60),labels=xtcks,las=2)
	lines(100*lmead[,2]/subtot, col="blue")
	lines(100*lpowell/subtot, col="purple")
		title(main="Colorado River Storage\nLake Mead and Lake Powell")
	abline(h=10,col="gray",lty=2)
	abline(h=20,col="gray")
	abline(h=30,col="gray",lty=2)
	abline(h=40,col="gray")
	abline(h=50,col="gray",lty=2)
	abline(h=60,col="gray")
	abline(h=70,col="gray",lty=2)
	abline(h=80,col="gray")
	abline(h=90,col="gray",lty=2)
	abline(h=100,col="gray")
	
	legend("topleft", legend = c("Lake Mead","Lake Powell"), col=c("blue","purple"),
			bty="n", xjust = 1, yjust = 1, cex=0.9, lty=1)
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")

if (graphic_output) dev.off()