# temperature model
# 5.31.2019

tempmod <- function(dat.merge) {

    require(mgcv)
    require(maps)
    require(mapproj)
    require(akima)

    # drop littoral duplicates
    dat.merge <- subset(dat.merge, sample.type == "MICX")

    dat.merge.sav <- dat.merge
    # select only variables that we need
    dat.merge <- dat.merge[, c("site.id.c",
                                "therm.sav",
                               "domean", "tempmean","tbot",
                               "doint", "do.hyp.fail", "do.hyp.pass",
                               "yday.x", "st.nla2012",
                               "year", "tmax", "index.lat.dd", "tmean.pt",
                               "tmin.pt",
                 "index.lon.dd", "elevation", "area.ha", "index.site.depth")]

    sitemean <- function(x, site, varname) {
        xm <- tapply(x, site, mean, na.rm = T)
        return(xm)
    }
    ## compute mean values of lake characteristics and merge
    varlist <- c( "index.lat.dd", "tmean.pt", "tmin.pt",
                 "index.lon.dd", "elevation", "area.ha", "index.site.depth")
    matout <- matrix(NA, ncol = length(varlist),
                     nrow = length(levels(dat.merge$site.id.c)))
    for (i in 1:length(varlist)) {
        matout[,i] <- sitemean(dat.merge[, varlist[i]],
                               dat.merge$site.id.c)
    }
    dat.mean <- data.frame(levels(dat.merge$site.id.c), matout)
    names(dat.mean) <- c("site.id.c", "Lat", "Temp", "TempMin",
                         "Lon", "Elev", "Area", "Depth")
    dat.mean$LakeRatio <- (dat.mean$Area*10000)^0.25/dat.mean$Depth

    dat.merge<-merge(dat.merge, dat.mean[, c("site.id.c", "Depth", "LakeRatio",
                                             "Temp", "Elev", "Lat", "Lon",
                                             "TempMin")],
                     by = "site.id.c")
    print(nrow(dat.merge))


    incvec <- !is.na(dat.merge$tmax) & ! is.na(dat.merge$yday.x) &
        ! is.na(dat.merge$Lat) & ! is.na(dat.merge$Elev) &
            ! is.na(dat.merge$Temp)

    dat.merge <- dat.merge[incvec,]

    ## two outliers excluded
    incvec <- dat.merge$tmax < 5
    dat.merge <- dat.merge[!incvec,]
    incvec <- dat.merge$tmax >  35
    dat.merge <- dat.merge[!incvec,]

    ## function to calculate RMS error
    rmserr <- function(x,y) sqrt(sum((x-y)^2)/length(x))

    ## model for mean air temperature as function of location and elevation
    ## this model can be used with the parameters b[1] and b[2]
    ## from the model to predict the first day of stratification for
    ## an arbitrary location
    modairtemp <- gam(Temp ~ s(Lat, Lon, k = 30) + Elev, data = dat.merge)
    print(summary(modairtemp))

    ## plot predicted vs. observed
    predout <- predict(modairtemp)
    dev.new()
    par(mar = c(4,4,1,1), mgp = c(2.3,1,0), bty = "l", mfrow = c(1,2),
        bty = "l")
    plot(predout, dat.merge$Temp, xlab = "Predicted mean air temp",
         ylab = "Observed mean air temp", pch = 21, col = "grey39",
         bg = "white")
    abline(0,1)
    cat("Mean air temp RMS err:", rmserr(predout, dat.merge$Temp), "\n")
    save(modairtemp, file = "modairtemp.rda")

    ## model for minimum air temperature
    df1 <- na.omit(dat.mean[, c("Lon", "Lat", "Elev", "TempMin")])
    modmintemp <- gam(TempMin ~ s(Lon, Lat, k = 100) + Elev, data = df1)
    print(summary(modmintemp))
    predout <- predict(modmintemp, type = "response")

    plot(predout, df1$TempMin, xlab ="Predicted min air temp",
         ylab = "Observed min air temp", pch = 21, col = "grey39",
         bg = "white")
    abline(0,1)
    cat("Min air temp RMS err:", rmserr(predout, df1$TempMin), "\n")

    save(modmintemp, file = "modmintemp.rda")

    ## model for lake surface water temperature
    mod.watertemp <- gam(tmax ~ s(yday.x, k =7) +
                             s(Lat, Lon, k = 30) + Elev,
                         data = dat.merge)
    print(summary(mod.watertemp))
    predout <- predict(mod.watertemp)
    dev.new()
    par(mar = c(4,4,1,1), mgp = c(2.3,1,0), bty="l")
    plot(predout, dat.merge$tmax, pch = 21, col = "grey39", bg = "white",
         xlab = "Predicted surface layer temperature",
         ylab = "Observed surface layer temperature")
    abline(0,1)
    cat("Water temperature RMS error:", rmserr(predout, dat.merge$tmax), "\n")

    save(mod.watertemp, file = "mod.watertemp.rda")

    new.data <- dat.merge
    new.data$Elev <- mean(new.data$Elev)
    new.data$yday.x <- mean(new.data$yday.x)
    predout <- predict(mod.watertemp, new.data)

    ## plot map of location effect (Fig 14)
    require(maps)
    require(akima)
    png(width = 5, height = 3, pointsize = 8, units = "in",
        file = "laketempmap.png", res = 600)
    par(mar = c(1,1,1,1))
    map("state", proj = "albers", par = c(30,40), col = "grey")
    pout <- mapproject(dat.merge$Lon, dat.merge$Lat, proj = "")
    contour(interp(pout$x, pout$y, predout, duplicate = "mean"),add = T)
    dev.off()

    ## plot two other predictors (Fig 15)
    png(width = 6, height = 2.5, pointsize = 8, units = "in",
        res = 600, file ="temppreds.png")
    par(mar = c(4,4,1,1), mfrow = c(1,2), mgp = c(2.3,1,0), bty = "l")
    xnew <- seq(min(dat.merge$yday.x), max(dat.merge$yday.x), length = 100)
    new.data <- data.frame(yday.x = xnew)
    new.data$Elev <- mean(dat.merge$Elev)
    new.data$Lat <- mean(dat.merge$Lat)
    new.data$Lon <- mean(dat.merge$Lon)
    predout <- predict(mod.watertemp, new.data, se.fit = T)
    up <- predout$fit + 2*predout$se.fit
    dn <- predout$fit - 2*predout$se.fit
    plot(xnew, predout$fit, ylim = range(c(up, dn)),
         type = "n", xlab = "Day of year", ylab = "Lake temperature")
    polygon(c(xnew, rev(xnew)), c(up, rev(dn)), col = "grey",
            border = NA)
    lines(xnew, predout$fit)

    xnew <- seq(min(dat.merge$Elev), max(dat.merge$Elev), length = 100)
    new.data <- data.frame(Elev = xnew)
    new.data$yday.x <- mean(dat.merge$yday.x)
    new.data$Lat <- mean(dat.merge$Lat)
    new.data$Lon <- mean(dat.merge$Lon)
    predout <- predict(mod.watertemp, new.data, se.fit = T)
    up <- predout$fit + 2*predout$se.fit
    dn <- predout$fit - 2*predout$se.fit
    plot(xnew, predout$fit, ylim = range(c(up, dn)),
         type = "n", xlab = "Elevation (m)", ylab = "Lake temperature")
    polygon(c(xnew, rev(xnew)), c(up, rev(dn)), col = "grey",
            border = NA)
    lines(xnew, predout$fit)
    dev.off()

    cc <- coef(mod.watertemp)
    ## find day of year for maximum temperature
    xnew <- seq(min(dat.merge$yday.x), max(dat.merge$yday.x), by = 1)
    new.data <- data.frame(yday.x =xnew)
    new.data$Lat <- mean(dat.merge$Lat)
    new.data$Lon <- mean(dat.merge$Lon)
    new.data$Elev <- mean(dat.merge$Elev)
    predout <- predict(mod.watertemp, new.data)
    ip <- which(predout == max(predout))
    daymax <- xnew[ip]

    predout.terms <- predict(mod.watertemp, type = "terms")
    dat.merge$loceff <- predout.terms[,3] # contribution of location for each NLA sample

    predout.terms <- predict(mod.watertemp, new.data, type = "terms")

    ## extract yday curve for decreasing temperature time of year
    daycurve <- predout.terms[ip:nrow(predout.terms), 2]

    ## find critical day of the year for each lat and elev
    tcrit <- 24
    y <- tcrit - dat.merge$loceff - cc["Elev"]*dat.merge$Elev -
        cc["(Intercept)"] # adjust critical temperature with location
    # elevation and intercept effects

    ## find day that corresponds with critical temp
    daycrit <- rep(NA, times = nrow(dat.merge))
    for (i in 1:nrow(dat.merge)) {
        daycrit[i] <- approx(daycurve, xnew[ip:nrow(predout.terms)],  y[i])$y
    }
    daycrit <- round(daycrit)

    ## plot temperature release day for critical temperature of 24 (Fig 16)
    ## Adjust temperatures for other critical temperatures
    png(width = 5, height = 3, pointsize = 8, units = "in", res = 600,
        file = "yday24C.png")
    par(mar = c(1,1,1,1))
    map("state", proj = "albers", par = c(30,40), col = "grey")
    pout <- mapproject(dat.merge$Lon, dat.merge$Lat, proj = "")

    require(akima)
    incvec <- !is.na(daycrit)
    points(pout$x[!incvec], pout$y[!incvec], pch = ".")
    points(pout$x[incvec], pout$y[incvec], pch = 21, col = "grey", bg = "white")
    contour(interp(pout$x[incvec], pout$y[incvec], daycrit[incvec],
                   duplicate = "mean"),
            add = T)
    dev.off()

    return()

}

tempmod(dat.merge.all)
