## national lake temperature and air temperature models
## 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", "doc.result",
                                "therm.sav",
                               "domean", "tempmean","tbot",
                                "do.hyp.fail", "do.hyp.pass",
                               "yday.x", "st.nla2012",
                               "year", "tmax",
                               "index.site.depth", "tmean.pt",
                               "elevation", "index.lat.dd",
                               "index.lon.dd", "tmin.pt")]

    ## calculated site mean values
    varlist <- c("index.site.depth", "tmean.pt", "elevation",
                 "index.lat.dd", "index.lon.dd", "tmin.pt")
    varlist2 <- c("Depth", "Temp", "Elev", "Lat", "Lon",
                  "TempMin")
    matout <- matrix(NA, ncol = length(varlist),
                     nrow = length(levels(dat.merge$site.id.c)))
    for (i in 1:length(varlist)) {
        matout[,i] <- as.vector(tapply(dat.merge[, varlist[i]],
                                       dat.merge$site.id.c, mean, na.rm = T))
    }
    dfmean <- data.frame(levels(dat.merge$site.id.c), matout)
    names(dfmean) <- c("site.id.c", varlist2)

    ## drop 2 missing temperatures
    incvec <- !is.na(dfmean$Temp)
    dfmean <- dfmean[incvec,]
    print(summary(dfmean))

    ## model for minimum air temperature
    ## fill in missing TempMin with GAM
    dftemp <- na.omit(dfmean[, c("Lat", "Lon", "Elev", "TempMin")])
    modmintemp <- gam(TempMin ~ s(Lat, Lon, k = 100) + Elev, data = dftemp)
    print(summary(modmintemp))

    predout <- predict(modmintemp, dfmean[, c("Lat", "Lon", "Elev")])
    dfmean$TempMin.old <- dfmean$TempMin
    dfmean$TempMin <- predout

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

    ## 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 = dfmean)
    print(summary(modairtemp))
    predout <- predict(modairtemp)
    save(modairtemp, file = "modairtemp.rda")

    dat.merge<-merge(dat.merge, dfmean, by = "site.id.c")
    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 in lake temperature excluded
    incvec <- dat.merge$tmax < 5
    dat.merge <- dat.merge[!incvec,]
    incvec <- dat.merge$tmax >  35
    dat.merge <- dat.merge[!incvec,]

    ## model for maximum lake temperature as a function of location,
    ## elevation, and day of the year
    mod2 <- gam(tmax ~ s(yday.x, k =7) + s(Lat, Lon, k = 30) + Elev,
                data = dat.merge)
    print(summary(mod2))

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

    save(mod2, file = "tempmod.rda")

    ## show map of location effect
    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
    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(mod2, 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(mod2, 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()

    return()

}

tempmod(dat.merge.all)
