### R Script accompanying the following dataset https://doi.org/10.7923/6Z9M-WZ55 
    # this file includes necessary code for kriging to predict foodscapes at the Cedar Gulch Idaho site.

### Install necessary packages if needed ###
install.packages("dplyr",dependencies = TRUE)
install.packages("gstat",dependencies = TRUE)
install.packages("nlme",dependencies = TRUE)
install.packages("geoR",dependencies = TRUE)
  ### If using Mac OS you may need to reinstall XQuartz in order to 
  ### activate tcltk as a dependency of geoR
  ### https://www.xquartz.org/
install.packages("raster",dependencies = TRUE)
install.packages("rgdal",dependencies = TRUE)
install.packages("sp",dependencies = TRUE)

# Load required packages                 
library(dplyr)
library(gstat)
library(nlme)
library(geoR)
source("VarioScript3.txt")
source("KrigeWrapper3.txt")

# Load data
data <- read.csv("gps-diet-quality-feb2019.csv", header=TRUE)
summary(data)

#### Prediction Grid ####
library(sp)
library(raster)
library(rgdal)

## Cedar Gulch
cg.hab <- readGDAL("CedarGulch_PatchTypeClassification_20150605.tif") #load in Cedar Veg raster
cg.hab <- as(cg.hab, "SpatialGridDataFrame")
gc()
cg.hab.grid <- as.data.frame(cg.hab)[,c(2,3,1)] #convert to data frame to use in predict
colnames(cg.hab.grid) <- c("Easting", "Northing", "Veg_type")
cg.hab.grid$Veg_type <- factor(cg.hab.grid$Veg_type,
                               labels=c("on","off","dwarf")) # 0=on, 1=off, 2=dwarf
# reorder Veg_type so that it matches the order in data (otherwise it thinks on=dwarf and dwarf=on)
cg.hab.grid$Veg_type <- factor(cg.hab.grid$Veg_type, levels(cg.hab.grid$Veg_type)[c(3,2,1)])
gc()

# Partition prediction grid so computer doesn't crash
cg.hab.grid2.1 <- cg.hab.grid[1:1000000, ] # 1 million locs seems good
cg.hab.grid2.2 <- cg.hab.grid[1000001:2000000, ]
cg.hab.grid2.3 <- cg.hab.grid[2000001:3000000, ]
cg.hab.grid2.4 <- cg.hab.grid[3000001:4000000, ]
cg.hab.grid2.5 <- cg.hab.grid[4000001:5000000, ]
cg.hab.grid2.6 <- cg.hab.grid[5000001:6000000, ]
cg.hab.grid2.7 <- cg.hab.grid[6000001:7000000, ]
cg.hab.grid2.8 <- cg.hab.grid[7000001:8000000, ]
cg.hab.grid2.9 <- cg.hab.grid[8000001:9000000, ]
cg.hab.grid2.10 <- cg.hab.grid[9000001:10000000, ]
cg.hab.grid2.11 <- cg.hab.grid[10000001:11000000, ]
cg.hab.grid2.12 <- cg.hab.grid[11000001:12000000, ]
cg.hab.grid2.13 <- cg.hab.grid[12000001:13000000, ]
cg.hab.grid2.14 <- cg.hab.grid[13000001:14000000, ]
cg.hab.grid2.15 <- cg.hab.grid[14000001:15000000, ]
cg.hab.grid2.16 <- cg.hab.grid[15000001:16000000, ]
cg.hab.grid2.17 <- cg.hab.grid[16000001:17000000, ]
cg.hab.grid2.18 <- cg.hab.grid[17000001:18000000, ]
cg.hab.grid2.19 <- cg.hab.grid[18000001:19000000, ]
cg.hab.grid2.20 <- cg.hab.grid[19000001:20000000, ]
cg.hab.grid2.21 <- cg.hab.grid[20000001:21000000, ]
cg.hab.grid2.22 <- cg.hab.grid[21000001:22000000, ]
cg.hab.grid2.23 <- cg.hab.grid[22000001:23000000, ]
cg.hab.grid2.24 <- cg.hab.grid[23000001:24000000, ]
cg.hab.grid2.25 <- cg.hab.grid[24000001:25000000, ]
cg.hab.grid2.26 <- cg.hab.grid[25000001:26000000, ]
cg.hab.grid2.27 <- cg.hab.grid[26000001:27000000, ]
cg.hab.grid2.28 <- cg.hab.grid[27000001:28000000, ]
cg.hab.grid2.29 <- cg.hab.grid[28000001:29000000, ]
cg.hab.grid2.30 <- cg.hab.grid[29000001:30000000, ]
cg.hab.grid2.31 <- cg.hab.grid[30000001:31000000, ]
cg.hab.grid2.32 <- cg.hab.grid[31000001:31466001, ]
gc()

#### CP Winter CG ####
cp.win.cg <- data %>%
  filter(!is.na(CP_Win), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, CP_Win)

## phi = ~100
## relative nugget (nugget / (nugget + sill)) = (1.615/(1.615+1.124)) = ~0.5
cp.win.cg.gls <- gls(CP_Win ~ Easting + Northing + Veg_type, data=cp.win.cg, correlation = corSpher(value=c(100, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(cp.win.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
cp.win.cg.gls <- gls(CP_Win ~ Veg_type, data=cp.win.cg, correlation = corSpher(value=c(100, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(cp.win.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(cp.win.cg.gls)
plot(cp.win.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(cp.win.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(cp.win.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
cp.win.cg.gls2 <- gls(CP_Win ~ Veg_type, data=cp.win.cg) ## non-spatial model
anova(cp.win.cg.gls2, cp.win.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
cp.win.cg.gd <- as.geodata(cp.win.cg[,c(5,4,8)])
cp.win.cg.gd$data <- resid(cp.win.cg.gls, resType="response")[1:256]
#cp.win.cg.gd$residuals <- resid(cp.win.cg.gls, resType="response")[1:138]

## plot the gls residuals
cp.win.cg.vg <- variog(cp.win.cg.gd,max.dist=300)
plot(cp.win.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
cp.win.cg.fit <- variofit(cp.win.cg.vg, ini.cov.pars=c(1.5, 100), nugget=1.5, cov.model="spherical", weights="cressie", messages=F)
summary(cp.win.cg.fit)
lines(cp.win.cg.fit)

## Krige with covariance model and trend model
# Make sure that the "data" has both the CP_Win and Veg_type columns
cp.win.cg.uk2.1 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.1, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.2 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.2, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.3 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.3, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.4 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.4, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.5 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.5, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.6 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.6, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.7 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.7, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.8 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.8, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.9 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.9, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.10 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.10, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.11 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.11, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.12 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.12, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.13 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.13, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.14 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.14, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.15 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.15, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.16 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.16, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.17 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.17, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.18 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.18, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.19 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.19, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.20 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.20, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.21 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.21, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.22 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.22, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.23 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.23, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.24 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.24, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.25 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.25, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.26 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.26, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.27 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.27, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.28 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.28, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.29 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.29, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.30 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.30, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.31 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.31, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()
cp.win.cg.uk2.32 <- krige(data=cp.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.32, cov.mod=cp.win.cg.fit, trend.mod=cp.win.cg.gls)
gc()

# Combine predictions back into the cg.hab.grid data frame
cg.hab.grid$CP_Win <- c(cp.win.cg.uk2.1$predict, cp.win.cg.uk2.2$predict, 
                        cp.win.cg.uk2.3$predict, cp.win.cg.uk2.4$predict, 
                        cp.win.cg.uk2.5$predict, cp.win.cg.uk2.6$predict, 
                        cp.win.cg.uk2.7$predict, cp.win.cg.uk2.8$predict,
                        cp.win.cg.uk2.9$predict, cp.win.cg.uk2.10$predict,
                        cp.win.cg.uk2.11$predict, cp.win.cg.uk2.12$predict, 
                        cp.win.cg.uk2.13$predict, cp.win.cg.uk2.14$predict, 
                        cp.win.cg.uk2.15$predict, cp.win.cg.uk2.16$predict, 
                        cp.win.cg.uk2.17$predict, cp.win.cg.uk2.18$predict,
                        cp.win.cg.uk2.19$predict, cp.win.cg.uk2.20$predict,
                        cp.win.cg.uk2.21$predict, cp.win.cg.uk2.22$predict, 
                        cp.win.cg.uk2.23$predict, cp.win.cg.uk2.24$predict, 
                        cp.win.cg.uk2.25$predict, cp.win.cg.uk2.26$predict, 
                        cp.win.cg.uk2.27$predict, cp.win.cg.uk2.28$predict,
                        cp.win.cg.uk2.29$predict, cp.win.cg.uk2.30$predict,
                        cp.win.cg.uk2.31$predict, cp.win.cg.uk2.32$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
cg.hab.grid$CP_Win_std <- sqrt(c(cp.win.cg.uk2.1$krige.var, cp.win.cg.uk2.2$krige.var, 
                                 cp.win.cg.uk2.3$krige.var, cp.win.cg.uk2.4$krige.var, 
                                 cp.win.cg.uk2.5$krige.var, cp.win.cg.uk2.6$krige.var, 
                                 cp.win.cg.uk2.7$krige.var, cp.win.cg.uk2.8$krige.var,
                                 cp.win.cg.uk2.9$krige.var,cp.win.cg.uk2.10$krige.var,
                                 cp.win.cg.uk2.11$krige.var, cp.win.cg.uk2.12$krige.var, 
                                 cp.win.cg.uk2.13$krige.var, cp.win.cg.uk2.14$krige.var, 
                                 cp.win.cg.uk2.15$krige.var, cp.win.cg.uk2.16$krige.var, 
                                 cp.win.cg.uk2.17$krige.var, cp.win.cg.uk2.18$krige.var,
                                 cp.win.cg.uk2.19$krige.var,cp.win.cg.uk2.20$krige.var,
                                 cp.win.cg.uk2.21$krige.var, cp.win.cg.uk2.22$krige.var, 
                                 cp.win.cg.uk2.23$krige.var, cp.win.cg.uk2.24$krige.var, 
                                 cp.win.cg.uk2.25$krige.var, cp.win.cg.uk2.26$krige.var, 
                                 cp.win.cg.uk2.27$krige.var, cp.win.cg.uk2.28$krige.var,
                                 cp.win.cg.uk2.29$krige.var,cp.win.cg.uk2.30$krige.var,
                                 cp.win.cg.uk2.21$krige.var, cp.win.cg.uk2.32$krige.var))

# # Plot prediction -- note this may take a long time or won't work on some computers
# quilt.plot(cg.hab.grid[,1:2], cg.hab.grid$CP_Win, col=heat.colors(24), add.legend=T)
# points(c.cp.win.cg.gd, pch=19, add=T)

# # Alternate plot with ggplot2, less pixellation and better colormap options
# library(ggplot2)
# ggplot() + geom_raster(data=cg.hab.grid, aes(x=Easting, y=Northing, fill=CP_Win)) + 
#   scale_fill_viridis_c() + coord_quickmap()
# ggplot() + geom_histogram(data=cg.hab.grid, aes(CP_Win, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cg.hab.grid, aes(CP_Win_std, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cp.win.cg, aes(CP_Win, fill=Veg_type))

# Export as GEOTIFF
cg.hab.export <- rasterFromXYZ(cg.hab.grid[,c("Easting","Northing","CP_Win")]) #go from dataframe back to raster
#plot(cg.hab.export) #looks good
crs(cg.hab.export) <- crs(cg.hab) #add back the coordinate reference system

writeRaster(cg.hab.export,'CG_Winter_CP_Kriged_UK.tif',options=c('TFW=YES'))


#### CP Summer CG ####
cp.sum.cg <- data %>%
  filter(!is.na(CP_Sum), Site=="Cedar Gulch", CP_Sum < 25) %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, CP_Sum)

## phi = 120
## relative nugget (nugget / (nugget + sill)) = (1.352/(1.352+0.191)) = 0.88
cp.sum.cg.gls <- gls(CP_Sum ~ Easting + Northing + Veg_type, data=cp.sum.cg, correlation = corSpher(value=c(125, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(cp.sum.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
cp.sum.cg.gls <- gls(CP_Sum ~ Veg_type, data=cp.sum.cg, correlation = corSpher(value=c(125, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(cp.sum.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(cp.sum.cg.gls)
plot(cp.sum.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(cp.sum.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(cp.sum.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
cp.sum.cg.gls2 <- gls(CP_Sum ~ Veg_type, data=cp.sum.cg) ## non-spatial model
anova(cp.sum.cg.gls2, cp.sum.cg.gls)
## the covariance model is significant at p=0.01, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
cp.sum.cg.gd <- as.geodata(cp.sum.cg[,c(5,4,8)])
cp.sum.cg.gd$data <- resid(cp.sum.cg.gls, resType="response")[1:229]

## plot the gls residuals
cp.sum.cg.vg <- variog(cp.sum.cg.gd,max.dist=300)
plot(cp.sum.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
cp.sum.cg.fit <- variofit(cp.sum.cg.vg, ini.cov.pars=c(1.37, 90), nugget=3.4, cov.model="spherical", weights="cressie", messages=F)
summary(cp.sum.cg.fit)
lines(cp.sum.cg.fit)

## Krige with covariance model and trend model
# Make sure that the "data" has both the CP_Sum and Veg_type columns
cp.sum.cg.uk2.1 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.1, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.2 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.2, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.3 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.3, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.4 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.4, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.5 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.5, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.6 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.6, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.7 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.7, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.8 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.8, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.9 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.9, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.10 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.10, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.11 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.11, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.12 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.12, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.13 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.13, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.14 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.14, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.15 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.15, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.16 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.16, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.17 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.17, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.18 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.18, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.19 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.19, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.20 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.20, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.21 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.21, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.22 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.22, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.23 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.23, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.24 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.24, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.25 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.25, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.26 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.26, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.27 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.27, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.28 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.28, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.29 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.29, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.30 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.30, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.31 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.31, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()
cp.sum.cg.uk2.32 <- krige(data=cp.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.32, cov.mod=cp.sum.cg.fit, trend.mod=cp.sum.cg.gls)
gc()

# Combine predictions back into the cg.hab.grid data frame
cg.hab.grid$CP_Sum <- c(cp.sum.cg.uk2.1$predict, cp.sum.cg.uk2.2$predict, 
                        cp.sum.cg.uk2.3$predict, cp.sum.cg.uk2.4$predict, 
                        cp.sum.cg.uk2.5$predict, cp.sum.cg.uk2.6$predict, 
                        cp.sum.cg.uk2.7$predict, cp.sum.cg.uk2.8$predict,
                        cp.sum.cg.uk2.9$predict, cp.sum.cg.uk2.10$predict,
                        cp.sum.cg.uk2.11$predict, cp.sum.cg.uk2.12$predict, 
                        cp.sum.cg.uk2.13$predict, cp.sum.cg.uk2.14$predict, 
                        cp.sum.cg.uk2.15$predict, cp.sum.cg.uk2.16$predict, 
                        cp.sum.cg.uk2.17$predict, cp.sum.cg.uk2.18$predict,
                        cp.sum.cg.uk2.19$predict, cp.sum.cg.uk2.20$predict,
                        cp.sum.cg.uk2.21$predict, cp.sum.cg.uk2.22$predict, 
                        cp.sum.cg.uk2.23$predict, cp.sum.cg.uk2.24$predict, 
                        cp.sum.cg.uk2.25$predict, cp.sum.cg.uk2.26$predict, 
                        cp.sum.cg.uk2.27$predict, cp.sum.cg.uk2.28$predict,
                        cp.sum.cg.uk2.29$predict, cp.sum.cg.uk2.30$predict,
                        cp.sum.cg.uk2.31$predict, cp.sum.cg.uk2.32$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
cg.hab.grid$CP_Sum_std <- sqrt(c(cp.sum.cg.uk2.1$krige.var, cp.sum.cg.uk2.2$krige.var, 
                                 cp.sum.cg.uk2.3$krige.var, cp.sum.cg.uk2.4$krige.var, 
                                 cp.sum.cg.uk2.5$krige.var, cp.sum.cg.uk2.6$krige.var, 
                                 cp.sum.cg.uk2.7$krige.var, cp.sum.cg.uk2.8$krige.var,
                                 cp.sum.cg.uk2.9$krige.var,cp.sum.cg.uk2.10$krige.var,
                                 cp.sum.cg.uk2.11$krige.var, cp.sum.cg.uk2.12$krige.var, 
                                 cp.sum.cg.uk2.13$krige.var, cp.sum.cg.uk2.14$krige.var, 
                                 cp.sum.cg.uk2.15$krige.var, cp.sum.cg.uk2.16$krige.var, 
                                 cp.sum.cg.uk2.17$krige.var, cp.sum.cg.uk2.18$krige.var,
                                 cp.sum.cg.uk2.19$krige.var,cp.sum.cg.uk2.20$krige.var,
                                 cp.sum.cg.uk2.21$krige.var, cp.sum.cg.uk2.22$krige.var, 
                                 cp.sum.cg.uk2.23$krige.var, cp.sum.cg.uk2.24$krige.var, 
                                 cp.sum.cg.uk2.25$krige.var, cp.sum.cg.uk2.26$krige.var, 
                                 cp.sum.cg.uk2.27$krige.var, cp.sum.cg.uk2.28$krige.var,
                                 cp.sum.cg.uk2.29$krige.var,cp.sum.cg.uk2.30$krige.var,
                                 cp.sum.cg.uk2.21$krige.var, cp.sum.cg.uk2.32$krige.var))

# # Plot prediction -- note this may take a long time or won't work on some computers
# quilt.plot(cg.hab.grid[,1:2], cg.hab.grid$CP_Sum, col=heat.colors(24), add.legend=T)
# points(c.cp.sum.cg.gd, pch=19, add=T)

# # Alternate plot with ggplot2, less pixellation and better colormap options
# library(ggplot2)
# ggplot() + geom_raster(data=cg.hab.grid, aes(x=Easting, y=Northing, fill=CP_Sum)) + 
#   scale_fill_viridis_c() + coord_quickmap()
# ggplot() + geom_histogram(data=cg.hab.grid, aes(CP_Sum, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cg.hab.grid, aes(CP_Sum_std, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cp.sum.cg, aes(CP_Sum, fill=Veg_type))

# Export as GEOTIFF
cg.hab.export <- rasterFromXYZ(cg.hab.grid[,c("Easting","Northing","CP_Sum")]) #go from dataframe back to raster
#plot(cg.hab.export) #looks good
crs(cg.hab.export) <- crs(cg.hab) #add back the coordinate reference system

writeRaster(cg.hab.export,'CG_Summer_CP_Kriged_UK.tif',options=c('TFW=YES'))


#### Monoterpenes Winter CG ####
mt.win.cg <- data %>%
  filter(!is.na(MT_Win), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, MT_Win)

## phi = 27
## relative nugget (nugget / (nugget + sill)) = (21630/(21630+32826)) = 0.40
mt.win.cg.gls <- gls(MT_Win ~ Easting + Northing + Veg_type, data=mt.win.cg, correlation = corSpher(value=c(125, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.win.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
mt.win.cg.gls <- gls(MT_Win ~ Veg_type, data=mt.win.cg, correlation = corSpher(value=c(50, 0.75), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.win.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(mt.win.cg.gls)
plot(mt.win.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(mt.win.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(mt.win.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
mt.win.cg.gls2 <- gls(MT_Win ~ Veg_type, data=mt.win.cg) ## non-spatial model
anova(mt.win.cg.gls2, mt.win.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
mt.win.cg.gd <- as.geodata(mt.win.cg[,c(5,4,8)])
mt.win.cg.gd$data <- resid(mt.win.cg.gls, resType="response")[1:273]

## plot the gls residuals
mt.win.cg.vg <- variog(mt.win.cg.gd,max.dist=300)
plot(mt.win.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
mt.win.cg.fit <- variofit(mt.win.cg.vg, ini.cov.pars=c(7000, 10), nugget=20000, cov.model="spherical", weights="cressie", messages=F)
summary(mt.win.cg.fit)
lines(mt.win.cg.fit)

## Krige with covariance model and trend model
# Make sure that the "data" has both the MT_Win and Veg_type columns
mt.win.cg.uk2.1 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.1, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.2 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.2, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.3 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.3, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.4 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.4, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.5 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.5, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.6 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.6, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.7 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.7, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.8 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.8, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.9 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.9, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.10 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.10, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.11 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.11, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.12 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.12, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.13 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.13, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.14 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.14, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.15 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.15, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.16 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.16, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.17 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.17, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.18 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.18, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.19 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.19, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.20 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.20, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.21 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.21, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.22 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.22, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.23 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.23, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.24 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.24, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.25 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.25, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.26 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.26, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.27 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.27, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.28 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.28, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.29 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.29, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.30 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.30, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.31 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.31, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()
mt.win.cg.uk2.32 <- krige(data=mt.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.32, cov.mod=mt.win.cg.fit, trend.mod=mt.win.cg.gls)
gc()

# Combine predictions back into the cg.hab.grid data frame
cg.hab.grid$MT_Win <- c(mt.win.cg.uk2.1$predict, mt.win.cg.uk2.2$predict, 
                        mt.win.cg.uk2.3$predict, mt.win.cg.uk2.4$predict, 
                        mt.win.cg.uk2.5$predict, mt.win.cg.uk2.6$predict, 
                        mt.win.cg.uk2.7$predict, mt.win.cg.uk2.8$predict,
                        mt.win.cg.uk2.9$predict, mt.win.cg.uk2.10$predict,
                        mt.win.cg.uk2.11$predict, mt.win.cg.uk2.12$predict, 
                        mt.win.cg.uk2.13$predict, mt.win.cg.uk2.14$predict, 
                        mt.win.cg.uk2.15$predict, mt.win.cg.uk2.16$predict, 
                        mt.win.cg.uk2.17$predict, mt.win.cg.uk2.18$predict,
                        mt.win.cg.uk2.19$predict, mt.win.cg.uk2.20$predict,
                        mt.win.cg.uk2.21$predict, mt.win.cg.uk2.22$predict, 
                        mt.win.cg.uk2.23$predict, mt.win.cg.uk2.24$predict, 
                        mt.win.cg.uk2.25$predict, mt.win.cg.uk2.26$predict, 
                        mt.win.cg.uk2.27$predict, mt.win.cg.uk2.28$predict,
                        mt.win.cg.uk2.29$predict, mt.win.cg.uk2.30$predict,
                        mt.win.cg.uk2.31$predict, mt.win.cg.uk2.32$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
cg.hab.grid$MT_Win_std <- sqrt(c(mt.win.cg.uk2.1$krige.var, mt.win.cg.uk2.2$krige.var, 
                                 mt.win.cg.uk2.3$krige.var, mt.win.cg.uk2.4$krige.var, 
                                 mt.win.cg.uk2.5$krige.var, mt.win.cg.uk2.6$krige.var, 
                                 mt.win.cg.uk2.7$krige.var, mt.win.cg.uk2.8$krige.var,
                                 mt.win.cg.uk2.9$krige.var,mt.win.cg.uk2.10$krige.var,
                                 mt.win.cg.uk2.11$krige.var, mt.win.cg.uk2.12$krige.var, 
                                 mt.win.cg.uk2.13$krige.var, mt.win.cg.uk2.14$krige.var, 
                                 mt.win.cg.uk2.15$krige.var, mt.win.cg.uk2.16$krige.var, 
                                 mt.win.cg.uk2.17$krige.var, mt.win.cg.uk2.18$krige.var,
                                 mt.win.cg.uk2.19$krige.var,mt.win.cg.uk2.20$krige.var,
                                 mt.win.cg.uk2.21$krige.var, mt.win.cg.uk2.22$krige.var, 
                                 mt.win.cg.uk2.23$krige.var, mt.win.cg.uk2.24$krige.var, 
                                 mt.win.cg.uk2.25$krige.var, mt.win.cg.uk2.26$krige.var, 
                                 mt.win.cg.uk2.27$krige.var, mt.win.cg.uk2.28$krige.var,
                                 mt.win.cg.uk2.29$krige.var,mt.win.cg.uk2.30$krige.var,
                                 mt.win.cg.uk2.21$krige.var, mt.win.cg.uk2.32$krige.var))

# # Plot prediction -- note this may take a long time or won't work on some computers
# quilt.plot(cg.hab.grid[,1:2], cg.hab.grid$MT_Win, col=heat.colors(24), add.legend=T)
# points(c.mt.win.cg.gd, pch=19, add=T)

# # Alternate plot with ggplot2, less pixellation and better colormap options
# library(ggplot2)
# ggplot() + geom_raster(data=cg.hab.grid, aes(x=Easting, y=Northing, fill=MT_Win)) + 
#   scale_fill_viridis_c() + coord_quickmap()
# ggplot() + geom_histogram(data=cg.hab.grid, aes(MT_Win, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cg.hab.grid, aes(MT_Win_std, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=mt.win.cg, aes(MT_Win, fill=Veg_type))

# Export as GEOTIFF
cg.hab.export <- rasterFromXYZ(cg.hab.grid[,c("Easting","Northing","MT_Win")]) #go from dataframe back to raster
#plot(cg.hab.export) #looks good
crs(cg.hab.export) <- crs(cg.hab) #add back the coordinate reference system

writeRaster(cg.hab.export,'CG_Winter_MT_Kriged_UK.tif',options=c('TFW=YES'))
gc()


#### Monoterpenes Summer CG ####
mt.sum.cg <- data %>%
  filter(!is.na(MT_Sum), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, MT_Sum)

## phi = 27
## relative nugget (nugget / (nugget + sill)) = (21630/(21630+32826)) = 0.40
mt.sum.cg.gls <- gls(MT_Sum ~ Easting + Northing + Veg_type, data=mt.sum.cg, correlation = corSpher(value=c(125, 0.5), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.sum.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
mt.sum.cg.gls <- gls(MT_Sum ~ Veg_type, data=mt.sum.cg, correlation = corSpher(value=c(50, 0.75), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.sum.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(mt.sum.cg.gls)
plot(mt.sum.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(mt.sum.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(mt.sum.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
mt.sum.cg.gls2 <- gls(MT_Sum ~ Veg_type, data=mt.sum.cg) ## non-spatial model
anova(mt.sum.cg.gls2, mt.sum.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
mt.sum.cg.gd <- as.geodata(mt.sum.cg[,c(5,4,8)])
mt.sum.cg.gd$data <- resid(mt.sum.cg.gls, resType="response")[1:240]

## plot the gls residuals
mt.sum.cg.vg <- variog(mt.sum.cg.gd,max.dist=300)
plot(mt.sum.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
mt.sum.cg.fit <- variofit(mt.sum.cg.vg, ini.cov.pars=c(10000, 100), nugget=40000, cov.model="spherical", weights="cressie", messages=F)
summary(mt.sum.cg.fit)
lines(mt.sum.cg.fit)

## Krige with covariance model and trend model
# Make sure that the "data" has both the MT_Win and Veg_type columns
mt.sum.cg.uk2.1 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.1, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.2 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.2, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.3 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.3, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.4 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.4, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.5 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.5, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.6 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.6, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.7 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.7, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.8 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.8, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.9 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.9, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.10 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.10, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.11 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.11, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.12 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.12, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.13 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.13, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.14 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.14, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.15 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.15, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.16 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.16, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.17 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.17, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.18 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.18, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.19 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.19, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.20 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.20, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.21 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.21, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.22 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.22, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.23 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.23, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.24 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.24, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.25 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.25, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.26 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.26, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.27 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.27, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.28 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.28, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.29 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.29, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.30 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.30, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.31 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.31, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()
mt.sum.cg.uk2.32 <- krige(data=mt.sum.cg[,c(5,4,8,2)], predict=cg.hab.grid2.32, cov.mod=mt.sum.cg.fit, trend.mod=mt.sum.cg.gls)
gc()

# Combine predictions back into the cg.hab.grid data frame
cg.hab.grid$MT_Sum <- c(mt.sum.cg.uk2.1$predict, mt.sum.cg.uk2.2$predict, 
                        mt.sum.cg.uk2.3$predict, mt.sum.cg.uk2.4$predict, 
                        mt.sum.cg.uk2.5$predict, mt.sum.cg.uk2.6$predict, 
                        mt.sum.cg.uk2.7$predict, mt.sum.cg.uk2.8$predict,
                        mt.sum.cg.uk2.9$predict, mt.sum.cg.uk2.10$predict,
                        mt.sum.cg.uk2.11$predict, mt.sum.cg.uk2.12$predict, 
                        mt.sum.cg.uk2.13$predict, mt.sum.cg.uk2.14$predict, 
                        mt.sum.cg.uk2.15$predict, mt.sum.cg.uk2.16$predict, 
                        mt.sum.cg.uk2.17$predict, mt.sum.cg.uk2.18$predict,
                        mt.sum.cg.uk2.19$predict, mt.sum.cg.uk2.20$predict,
                        mt.sum.cg.uk2.21$predict, mt.sum.cg.uk2.22$predict, 
                        mt.sum.cg.uk2.23$predict, mt.sum.cg.uk2.24$predict, 
                        mt.sum.cg.uk2.25$predict, mt.sum.cg.uk2.26$predict, 
                        mt.sum.cg.uk2.27$predict, mt.sum.cg.uk2.28$predict,
                        mt.sum.cg.uk2.29$predict, mt.sum.cg.uk2.30$predict,
                        mt.sum.cg.uk2.31$predict, mt.sum.cg.uk2.32$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
cg.hab.grid$MT_Sum_std <- sqrt(c(mt.sum.cg.uk2.1$krige.var, mt.sum.cg.uk2.2$krige.var, 
                                 mt.sum.cg.uk2.3$krige.var, mt.sum.cg.uk2.4$krige.var, 
                                 mt.sum.cg.uk2.5$krige.var, mt.sum.cg.uk2.6$krige.var, 
                                 mt.sum.cg.uk2.7$krige.var, mt.sum.cg.uk2.8$krige.var,
                                 mt.sum.cg.uk2.9$krige.var,mt.sum.cg.uk2.10$krige.var,
                                 mt.sum.cg.uk2.11$krige.var, mt.sum.cg.uk2.12$krige.var, 
                                 mt.sum.cg.uk2.13$krige.var, mt.sum.cg.uk2.14$krige.var, 
                                 mt.sum.cg.uk2.15$krige.var, mt.sum.cg.uk2.16$krige.var, 
                                 mt.sum.cg.uk2.17$krige.var, mt.sum.cg.uk2.18$krige.var,
                                 mt.sum.cg.uk2.19$krige.var,mt.sum.cg.uk2.20$krige.var,
                                 mt.sum.cg.uk2.21$krige.var, mt.sum.cg.uk2.22$krige.var, 
                                 mt.sum.cg.uk2.23$krige.var, mt.sum.cg.uk2.24$krige.var, 
                                 mt.sum.cg.uk2.25$krige.var, mt.sum.cg.uk2.26$krige.var, 
                                 mt.sum.cg.uk2.27$krige.var, mt.sum.cg.uk2.28$krige.var,
                                 mt.sum.cg.uk2.29$krige.var,mt.sum.cg.uk2.30$krige.var,
                                 mt.sum.cg.uk2.21$krige.var, mt.sum.cg.uk2.32$krige.var))

# # Plot prediction -- note this may take a long time or won't work on some computers
# quilt.plot(cg.hab.grid[,1:2], cg.hab.grid$MT_Win, col=heat.colors(24), add.legend=T)
# points(c.mt.sum.cg.gd, pch=19, add=T)

# # Alternate plot with ggplot2, less pixellation and better colormap options
# library(ggplot2)
# ggplot() + geom_raster(data=cg.hab.grid, aes(x=Easting, y=Northing, fill=MT_Win)) + 
#   scale_fill_viridis_c() + coord_quickmap()
# ggplot() + geom_histogram(data=cg.hab.grid, aes(MT_Win, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cg.hab.grid, aes(MT_Win_std, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=mt.sum.cg, aes(MT_Win, fill=Veg_type))

# Export as GEOTIFF
cg.hab.export <- rasterFromXYZ(cg.hab.grid[,c("Easting","Northing","MT_Sum")]) #go from dataframe back to raster
#plot(cg.hab.export) #looks good
crs(cg.hab.export) <- crs(cg.hab) #add back the coordinate reference system

writeRaster(cg.hab.export,'CG_Summer_MT_Kriged_UK.tif',options=c('TFW=YES'))
gc()


#### Coumarins Winter CG ####
co.win.cg <- data %>%
  filter(!is.na(Coum_Win), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Coum_Win)

## phi = 33
## relative nugget (nugget / (nugget + sill)) = (0/(0+0.425)) = 0
co.win.cg.gls <- gls(Coum_Win ~ Easting + Northing + Veg_type, data=co.win.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(co.win.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
co.win.cg.gls <- gls(Coum_Win ~ Veg_type, data=co.win.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(co.win.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(co.win.cg.gls)
plot(co.win.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(co.win.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(co.win.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
co.win.cg.gls2 <- gls(Coum_Win ~ Veg_type, data=co.win.cg) ## non-spatial model
anova(co.win.cg.gls2, co.win.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
co.win.cg.gd <- as.geodata(co.win.cg[,c(5,4,8)])
co.win.cg.gd$data <- resid(co.win.cg.gls, resType="response")[1:257]

## plot the gls residuals
co.win.cg.vg <- variog(co.win.cg.gd,max.dist=300)
plot(co.win.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
co.win.cg.fit <- variofit(co.win.cg.vg, ini.cov.pars=c(0.03, 30), nugget=0.03, cov.model="spherical", weights="cressie", messages=F)
summary(co.win.cg.fit)
lines(co.win.cg.fit)

## Krige with covariance model and trend model
# Make sure that the "data" has both the Coum_Win and Veg_type columns
co.win.cg.uk2.1 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.1, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.2 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.2, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.3 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.3, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.4 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.4, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.5 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.5, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.6 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.6, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.7 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.7, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.8 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.8, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.9 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.9, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.10 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.10, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.11 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.11, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.12 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.12, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.13 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.13, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.14 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.14, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.15 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.15, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.16 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.16, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.17 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.17, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.18 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.18, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.19 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.19, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.20 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.20, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.21 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.21, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.22 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.22, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.23 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.23, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.24 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.24, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.25 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.25, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.26 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.26, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.27 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.27, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.28 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.28, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.29 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.29, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.30 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.30, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.31 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.31, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()
co.win.cg.uk2.32 <- krige(data=co.win.cg[,c(5,4,8,2)], predict=cg.hab.grid2.32, cov.mod=co.win.cg.fit, trend.mod=co.win.cg.gls)
gc()

# Combine predictions back into the cg.hab.grid data frame
cg.hab.grid$Coum_Win <- c(co.win.cg.uk2.1$predict, co.win.cg.uk2.2$predict, 
                          co.win.cg.uk2.3$predict, co.win.cg.uk2.4$predict, 
                          co.win.cg.uk2.5$predict, co.win.cg.uk2.6$predict, 
                          co.win.cg.uk2.7$predict, co.win.cg.uk2.8$predict,
                          co.win.cg.uk2.9$predict, co.win.cg.uk2.10$predict,
                          co.win.cg.uk2.11$predict, co.win.cg.uk2.12$predict, 
                          co.win.cg.uk2.13$predict, co.win.cg.uk2.14$predict, 
                          co.win.cg.uk2.15$predict, co.win.cg.uk2.16$predict, 
                          co.win.cg.uk2.17$predict, co.win.cg.uk2.18$predict,
                          co.win.cg.uk2.19$predict, co.win.cg.uk2.20$predict,
                          co.win.cg.uk2.21$predict, co.win.cg.uk2.22$predict, 
                          co.win.cg.uk2.23$predict, co.win.cg.uk2.24$predict, 
                          co.win.cg.uk2.25$predict, co.win.cg.uk2.26$predict, 
                          co.win.cg.uk2.27$predict, co.win.cg.uk2.28$predict,
                          co.win.cg.uk2.29$predict, co.win.cg.uk2.30$predict,
                          co.win.cg.uk2.31$predict, co.win.cg.uk2.32$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
cg.hab.grid$Coum_Win_std <- sqrt(c(co.win.cg.uk2.1$krige.var, co.win.cg.uk2.2$krige.var, 
                                   co.win.cg.uk2.3$krige.var, co.win.cg.uk2.4$krige.var, 
                                   co.win.cg.uk2.5$krige.var, co.win.cg.uk2.6$krige.var, 
                                   co.win.cg.uk2.7$krige.var, co.win.cg.uk2.8$krige.var,
                                   co.win.cg.uk2.9$krige.var,co.win.cg.uk2.10$krige.var,
                                   co.win.cg.uk2.11$krige.var, co.win.cg.uk2.12$krige.var, 
                                   co.win.cg.uk2.13$krige.var, co.win.cg.uk2.14$krige.var, 
                                   co.win.cg.uk2.15$krige.var, co.win.cg.uk2.16$krige.var, 
                                   co.win.cg.uk2.17$krige.var, co.win.cg.uk2.18$krige.var,
                                   co.win.cg.uk2.19$krige.var,co.win.cg.uk2.20$krige.var,
                                   co.win.cg.uk2.21$krige.var, co.win.cg.uk2.22$krige.var, 
                                   co.win.cg.uk2.23$krige.var, co.win.cg.uk2.24$krige.var, 
                                   co.win.cg.uk2.25$krige.var, co.win.cg.uk2.26$krige.var, 
                                   co.win.cg.uk2.27$krige.var, co.win.cg.uk2.28$krige.var,
                                   co.win.cg.uk2.29$krige.var,co.win.cg.uk2.30$krige.var,
                                   co.win.cg.uk2.21$krige.var, co.win.cg.uk2.32$krige.var))

# # Plot prediction -- note this may take a long time or won't work on some computers
# quilt.plot(cg.hab.grid[,1:2], cg.hab.grid$Coum_Win, col=heat.colors(24), add.legend=T)
# points(c.co.win.cg.gd, pch=19, add=T)

# # Alternate plot with ggplot2, less pixellation and better colormap options
# library(ggplot2)
# ggplot() + geom_raster(data=cg.hab.grid, aes(x=Easting, y=Northing, fill=Coum_Win)) + 
#   scale_fill_viridis_c() + coord_quickmap()
# ggplot() + geom_histogram(data=cg.hab.grid, aes(Coum_Win, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=cg.hab.grid, aes(Coum_Win_std, fill=Veg_type), bins=60)
# #ggplot() + geom_histogram(data=co.win.cg, aes(Coum_Win, fill=Veg_type))

# Export as GEOTIFF
cg.hab.export <- rasterFromXYZ(cg.hab.grid[,c("Easting","Northing","Coum_Win")]) #go from dataframe back to raster
#plot(cg.hab.export) #looks good
crs(cg.hab.export) <- crs(cg.hab) #add back the coordinate reference system

writeRaster(cg.hab.export,'CG_Winter_CO_Kriged_UK.tif',options=c('TFW=YES'))

#### Coumarins Summer CG ####
co.sum.cg <- data %>%
  filter(!is.na(Coum_Sum), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Coum_Sum)

## phi = 33
## relative nugget (nugget / (nugget + sill)) = (0/(0+0.425)) = 0
co.sum.cg.gls <- gls(Coum_Sum ~ Easting + Northing + Veg_type, data=co.sum.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(co.sum.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
co.sum.cg.gls <- gls(Coum_Sum ~ Veg_type, data=co.sum.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(co.sum.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(co.sum.cg.gls)
plot(co.sum.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(co.sum.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(co.sum.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
co.sum.cg.gls2 <- gls(Coum_Sum ~ Veg_type, data=co.sum.cg) ## non-spatial model
anova(co.sum.cg.gls2, co.sum.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
co.sum.cg.gd <- as.geodata(co.sum.cg[,c(5,4,8)])
co.sum.cg.gd$data <- resid(co.sum.cg.gls, resType="response")[1:189]

## plot the gls residuals
co.sum.cg.vg <- variog(co.sum.cg.gd,max.dist=300)
plot(co.sum.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
co.sum.cg.fit <- variofit(co.sum.cg.vg, ini.cov.pars=c(0.01, 100), nugget=0.02, cov.model="spherical", weights="cressie", messages=F)
summary(co.sum.cg.fit)
lines(co.sum.cg.fit)

#### Bpinene Winter CG ####
bp.win.cg <- data %>%
  filter(!is.na(Bpine_Win), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Bpine_Win)

## phi = 33
## relative nugget (nugget / (nugget + sill)) = (0/(0+0.425)) = 0
bp.win.cg.gls <- gls(Bpine_Win ~ Easting + Northing + Veg_type, data=bp.win.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(bp.win.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
bp.win.cg.gls <- gls(Bpine_Win ~ Veg_type + Easting, data=bp.win.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(bp.win.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(bp.win.cg.gls)
plot(bp.win.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(bp.win.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(bp.win.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
bp.win.cg.gls2 <- gls(Bpine_Win ~ Veg_type + Easting, data=bp.win.cg) ## non-spatial model
anova(bp.win.cg.gls2, bp.win.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
bp.win.cg.gd <- as.geodata(bp.win.cg[,c(5,4,8)])
bp.win.cg.gd$data <- resid(bp.win.cg.gls, resType="response")[1:278]

## plot the gls residuals
bp.win.cg.vg <- variog(bp.win.cg.gd,max.dist=150)
plot(bp.win.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
bp.win.cg.fit <- variofit(bp.win.cg.vg, ini.cov.pars=c(700, 50), nugget=1800, cov.model="spherical", weights="cressie", messages=F)
summary(bp.win.cg.fit)
lines(bp.win.cg.fit)

#### Bpinene Summer CG ####
bp.sum.cg <- data %>%
  filter(!is.na(Bpine_Sum), Site=="Cedar Gulch") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Bpine_Sum)

## phi = 33
## relative nugget (nugget / (nugget + sill)) = (0/(0+0.425)) = 0
bp.sum.cg.gls <- gls(Bpine_Sum ~ Easting + Northing + Veg_type, data=bp.sum.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(bp.sum.cg.gls)

## Easting and Northing are not important covariates when the residual covariance is considered
## remove it and re-fit the model
bp.sum.cg.gls <- gls(Bpine_Sum ~ Veg_type, data=bp.sum.cg, correlation = corSpher(value=c(70, 0.7), form= ~ Easting + Northing, nugget = TRUE))
summary(bp.sum.cg.gls)

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(bp.sum.cg.gls)
plot(bp.sum.cg.gls, pch=16) ## no heteroscedasticity
plot(Variogram(bp.sum.cg.gls, form= ~Easting + Northing, resType="response"), pch=16) ## variogram is fine
plot(Variogram(bp.sum.cg.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## remaining error shows no correlation structure

## compare the spatial and non-spatial models
bp.sum.cg.gls2 <- gls(Bpine_Sum ~ Veg_type, data=bp.sum.cg) ## non-spatial model
anova(bp.sum.cg.gls2, bp.sum.cg.gls)
## the covariance model is significant at p<0.001, based on a likelihood ratio test

## bring gls residuals into a geodata object so that we can use them to krige in geoR
bp.sum.cg.gd <- as.geodata(bp.sum.cg[,c(5,4,8)])
bp.sum.cg.gd$data <- resid(bp.sum.cg.gls, resType="response")[1:252]

## plot the gls residuals
bp.sum.cg.vg <- variog(bp.sum.cg.gd,max.dist=150)
plot(bp.sum.cg.vg, pch=16)

## fit a variogram to gls residuals using our standard methods
bp.sum.cg.fit <- variofit(bp.sum.cg.vg, ini.cov.pars=c(800, 50), nugget=1500, cov.model="spherical", weights="cressie", messages=F)
summary(bp.sum.cg.fit)
lines(bp.sum.cg.fit)

#### Cross-validation ####

plot.Custom <- function(x, y, ...) {
  .limits <- range(x, y)
  plot(x, y, xlim = .limits, ylim = .limits, ...)
}

## CP Winter CG
cv.cp.win.cg <- xvalid(cp.win.cg.gd, model=cp.win.cg.fit) #LOOCV
summary(cv.cp.win.cg)

pred <- cp.win.cg.gls$fitted + cv.cp.win.cg$predicted
actual <- cp.win.cg.gls$fitted + resid(cp.win.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## CP Summer CG
cv.cp.sum.cg <- xvalid(cp.sum.cg.gd, model=cp.sum.cg.fit) #LOOCV
summary(cv.cp.sum.cg)

pred <- cp.sum.cg.gls$fitted + cv.cp.sum.cg$predicted
actual <- cp.sum.cg.gls$fitted + resid(cp.sum.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## MT Winter CG
cv.mt.win.cg <- xvalid(mt.win.cg.gd, model=mt.win.cg.fit) #LOOCV
summary(cv.mt.win.cg)

pred <- mt.win.cg.gls$fitted + cv.mt.win.cg$predicted
actual <- mt.win.cg.gls$fitted + resid(mt.win.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## MT Summer CG
cv.mt.sum.cg <- xvalid(mt.sum.cg.gd, model=mt.sum.cg.fit) #LOOCV
summary(cv.mt.sum.cg)

pred <- mt.sum.cg.gls$fitted + cv.mt.sum.cg$predicted
actual <- mt.sum.cg.gls$fitted + resid(mt.sum.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## CO Winter CG
cv.co.win.cg <- xvalid(co.win.cg.gd, model=co.win.cg.fit) #LOOCV
summary(cv.co.win.cg)

pred <- co.win.cg.gls$fitted + cv.co.win.cg$predicted
actual <- co.win.cg.gls$fitted + resid(co.win.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## Bpinene Winter CG
cv.bp.win.cg <- xvalid(bp.win.cg.gd, model=bp.win.cg.fit) #LOOCV
summary(cv.bp.win.cg)

pred <- bp.win.cg.gls$fitted + cv.bp.win.cg$predicted
actual <- bp.win.cg.gls$fitted + resid(bp.win.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging

## Bpinene Summer CG
cv.bp.sum.cg <- xvalid(bp.sum.cg.gd, model=bp.sum.cg.fit) #LOOCV
summary(cv.bp.sum.cg)

pred <- bp.sum.cg.gls$fitted + cv.bp.sum.cg$predicted
actual <- bp.sum.cg.gls$fitted + resid(bp.sum.cg.gls, resType="response")[1:NROW(pred)]
summary(lm(actual ~ pred)) #r^2_adj
plot.Custom(actual, pred, pch = 16, xlab="actual", ylab="predicted"); abline(lm(pred ~ actual)); abline(0,1, lty="dashed")
sqrt(sum((lm(actual ~ pred)$residuals)^2)/NROW(lm(actual ~ pred)$residuals)) #RMSE regression kriging


## r2 for GLS only models
summary(lm(cp.win.cg$CP_Win ~ cp.win.cg.gls$fitted))
summary(lm(cp.sum.cg$CP_Sum ~ cp.sum.cg.gls$fitted))
summary(lm(mt.win.cg$MT_Win ~ mt.win.cg.gls$fitted))
summary(lm(mt.sum.cg$MT_Sum ~ mt.sum.cg.gls$fitted))
summary(lm(co.win.cg$Coum_Win ~ co.win.cg.gls$fitted))
summary(lm(co.sum.cg$Coum_Sum ~ co.sum.cg.gls$fitted))