## Data source: http://mesonet.agron.iastate.edu/request/coop/fe.phtml
## First read in data and process

precip <- read.table("../data/changeme.txt", sep = ",", header=TRUE)

## We have a variable called day, but it is not in a format R would understand
precip$date <- as.Date(precip$day, format = "%Y/%m/%d")

precip$doy <- as.numeric(format(precip$date, "%j"))

head(precip)

summary(precip)

library(lattice)

precip$year <- as.numeric(format(precip$date, "%Y"))

precip2 <- precip[precip$year > 1991,]

xyplot(precip ~ doy | factor(year), data = precip2,
       panel = function(x,y,...){
         panel.xyplot(x,y, cex=0.5, pch = 16, ...)
         panel.abline(h = 3, col = "red", lwd = 2)
       })

## Create a rule with 3 days with more than 3 in

flood <- data.frame(Date = as.Date(precip2$date), precip.sum = NA)
flood$doy <- as.numeric(format(flood$Date, "%j"))
flood$year <- as.numeric(format(flood$Date, "%Y"))


for(i in 1:I(dim(precip2)[1]-2)){

  sum.3.precip <- sum(precip2[i:I(i+2),"precip"])

  if(is.na(sum.3.precip)){
    next
  }else{
    if(sum.3.precip > 3){
      flood[i:I(i+2),2] <- sum.3.precip
    }
  }
}

pdf("../figs/precip-sum.pdf")
xyplot(precip.sum ~ doy | factor(year) , data = flood,
       panel = function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(h = 3, lwd = 2, col = "red")
       })
dev.off()
    
## Predicted flooding in 1993, 1997, 1998, 1999, 2001,
##                       2002, 2003, 2004, 2005, 2006,
##                       2007, 2008, 2009, 2010


## Read in the gauge data
## http://waterdata.usgs.gov/ia/nwis/dv/?site_no=05471000&agency_cd=USGS&amp;referred_module=sw

gage <- read.table("../data/gage.txt", skip = 28, fill = TRUE)

## Accoring to the site over 9 feet is a flood

names(gage) <- c("agency", "site", "date", "height", "status")

gage$Date <- as.Date(gage$date, format="%Y-%m-%d")

gage$year <- format(gage$Date, format = "%Y")
gage$doy <- as.numeric(format(gage$Date, format = "%j"))

head(gage)

last.day <- as.Date("2010-10-02")
gage2 <- gage[gage$year > 1991 & gage$Date < last.day,]

pdf("../figs/gauge-height.pdf")
xyplot(height ~ doy | year , data = gage2,
       panel = function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(h = 9, lwd = 2, col = "red")
       })
dev.off()

obs <- na.omit(as.numeric(unique(gage2[gage2$height > 9,"year"])))

## Flooding 
years <- 1992:2010

obs <- ifelse(years %in% obs, 1, 0)

## Predicted

prd <- unique(flood[flood$precip.sum > 2.5,"year"])

pred <- ifelse(years %in% prd, 1, 0)

pred.obs <- data.frame(year = 1992:2010,
                  obs = obs, pred = pred)

## Selecting Days with predicted floods

flood$pred <- ifelse(!is.na(flood$precip.sum),1,NA)

gage2$obs <- ifelse(gage2$height > 9, 1, NA)


pdf("../figs/obs-pred.pdf")
xyplot(pred + gage2$obs ~ doy | factor(year), data = flood,
       pch = c(20,1), cex=1.5,
       col = c("black", "red"),
       ylab = "Agreement predicted and observed",
       xlab = "Day of the year",
       key = list(text = list(c("pred","obs")),
         pch = c(20,1),
         col = c("black","red"),
         points = TRUE,
         lty = 0,
         lines = FALSE))
dev.off()