### R Scripts accompanying the following dataset https://doi.org/10.7923/6Z9M-WZ55 
  # this file includes necessary code for kriging to predict foodscapes at the Camas Idaho site.
    # Note: The Camas Idaho site is also referred to as Magic Reservoir or MR within this script, 
    # and within related projects. Camas == Magic Reservoir.

### 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(raster)
library(rgdal)
library(sp)
mr.hab <- readGDAL("Camas_PatchTypeClassification_20130625.tif") #load in Camas Veg raster
mr.hab <- as(mr.hab, "SpatialGridDataFrame")
mr.hab.grid <- as.data.frame(mr.hab)[,c(2,3,1)] #convert to data frame to use in predict
colnames(mr.hab.grid) <- c("Easting", "Northing", "Veg_type")
mr.hab.grid$Veg_type <- factor(mr.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)
mr.hab.grid$Veg_type <- factor(mr.hab.grid$Veg_type, levels(mr.hab.grid$Veg_type)[c(3,2,1)])

# Partition prediction grid so computer doesn't crash
mr.hab.grid2.1 <- mr.hab.grid[1:1000000, ] # 1 million locs seems good
mr.hab.grid2.2 <- mr.hab.grid[1000001:2000000, ]
mr.hab.grid2.3 <- mr.hab.grid[2000001:3000000, ]
mr.hab.grid2.4 <- mr.hab.grid[3000001:4000000, ]
mr.hab.grid2.5 <- mr.hab.grid[4000001:5000000, ]
mr.hab.grid2.6 <- mr.hab.grid[5000001:6000000, ]
mr.hab.grid2.7 <- mr.hab.grid[6000001:7000000, ]
mr.hab.grid2.8 <- mr.hab.grid[7000001:8000000, ]
mr.hab.grid2.9 <- mr.hab.grid[8000001:9041226, ]


#### CP Winter MR ####
cp.win.mr <- data %>%
  filter(!is.na(CP_Win), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, CP_Win)
#  group_by(Patch_ID, Veg_type) %>% 
#  summarize(n(), Northing=mean(Northing), Easting=mean(Easting), CP_Win=mean(CP_Win))

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

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

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

## compare the spatial and non-spatial models
cp.win.mr.gls2 <- gls(CP_Win ~ Veg_type, data=cp.win.mr) ## non-spatial model
anova(cp.win.mr.gls2, cp.win.mr.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.mr.gd <- as.geodata(cp.win.mr[,c(5,4,6)])
cp.win.mr.gd$data <- resid(cp.win.mr.gls, resType="response")[1:138]
#cp.win.mr.gd$residuals <- resid(cp.win.mr.gls, resType="response")[1:138]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the CP_Win and Veg_type columns
cp.win.mr.uk2.1 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.2 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.3 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.4 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.5 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.6 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.7 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.8 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()
cp.win.mr.uk2.9 <- krige(data=cp.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cp.win.mr.fit, trend.mod=cp.win.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$CP_Win <- c(cp.win.mr.uk2.1$predict, cp.win.mr.uk2.2$predict, 
                         cp.win.mr.uk2.3$predict, cp.win.mr.uk2.4$predict, 
                         cp.win.mr.uk2.5$predict, cp.win.mr.uk2.6$predict, 
                         cp.win.mr.uk2.7$predict, cp.win.mr.uk2.8$predict,
                         cp.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$CP_Win_std <- sqrt(c(cp.win.mr.uk2.1$krige.var, cp.win.mr.uk2.2$krige.var, 
                        cp.win.mr.uk2.3$krige.var, cp.win.mr.uk2.4$krige.var, 
                        cp.win.mr.uk2.5$krige.var, cp.win.mr.uk2.6$krige.var, 
                        cp.win.mr.uk2.7$krige.var, cp.win.mr.uk2.8$krige.var,
                        cp.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_CP_Kriged_UK.tif',options=c('TFW=YES'))


#### CP Summer MR ####
cp.sum.mr <- data %>%
  filter(!is.na(CP_Sum), Site=="Camas") %>% 
  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.mr.gls <- gls(CP_Sum ~ Easting + Northing + Veg_type, data=cp.sum.mr, correlation = corSpher(value=c(120, 0.88), form= ~ Easting + Northing, nugget = TRUE))
summary(cp.sum.mr.gls)

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

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

## compare the spatial and non-spatial models
cp.sum.mr.gls2 <- gls(CP_Sum ~ Veg_type, data=cp.sum.mr) ## non-spatial model
anova(cp.sum.mr.gls2, cp.sum.mr.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.mr.gd <- as.geodata(cp.sum.mr[,c(5,4,8)])
cp.sum.mr.gd$data <- resid(cp.sum.mr.gls, resType="response")[1:201]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the CP_Sum and Veg_type columns
cp.sum.mr.uk2.1 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.2 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.3 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.4 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.5 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.6 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.7 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.8 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()
cp.sum.mr.uk2.9 <- krige(data=cp.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cp.sum.mr.fit, trend.mod=cp.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$CP_Sum <- c(cp.sum.mr.uk2.1$predict, cp.sum.mr.uk2.2$predict, 
                        cp.sum.mr.uk2.3$predict, cp.sum.mr.uk2.4$predict, 
                        cp.sum.mr.uk2.5$predict, cp.sum.mr.uk2.6$predict, 
                        cp.sum.mr.uk2.7$predict, cp.sum.mr.uk2.8$predict,
                        cp.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$CP_Sum_std <- sqrt(c(cp.sum.mr.uk2.1$krige.var, cp.sum.mr.uk2.2$krige.var, 
                                 cp.sum.mr.uk2.3$krige.var, cp.sum.mr.uk2.4$krige.var, 
                                 cp.sum.mr.uk2.5$krige.var, cp.sum.mr.uk2.6$krige.var, 
                                 cp.sum.mr.uk2.7$krige.var, cp.sum.mr.uk2.8$krige.var,
                                 cp.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_CP_Kriged_UK.tif',options=c('TFW=YES'))


#### Monoterpenes Winter MR ####
mt.win.mr <- data %>%
  filter(!is.na(MT_Win), Site=="Camas") %>% 
  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.mr.gls <- gls(MT_Win ~ Easting + Northing + Veg_type, data=mt.win.mr, correlation = corSpher(value=c(27, 0.40), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.win.mr.gls)

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

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

## compare the spatial and non-spatial models
mt.win.mr.gls2 <- gls(MT_Win ~ Veg_type, data=mt.win.mr) ## non-spatial model
anova(mt.win.mr.gls2, mt.win.mr.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.mr.gd <- as.geodata(mt.win.mr[,c(5,4,8)])
mt.win.mr.gd$data <- resid(mt.win.mr.gls, resType="response")[1:205]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the MT_Win and Veg_type columns
mt.win.mr.uk2.1 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.2 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.3 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.4 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.5 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.6 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.7 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.8 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()
mt.win.mr.uk2.9 <- krige(data=mt.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=mt.win.mr.fit, trend.mod=mt.win.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$MT_Win <- c(mt.win.mr.uk2.1$predict, mt.win.mr.uk2.2$predict, 
                        mt.win.mr.uk2.3$predict, mt.win.mr.uk2.4$predict, 
                        mt.win.mr.uk2.5$predict, mt.win.mr.uk2.6$predict, 
                        mt.win.mr.uk2.7$predict, mt.win.mr.uk2.8$predict,
                        mt.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$MT_Win_std <- sqrt(c(mt.win.mr.uk2.1$krige.var, mt.win.mr.uk2.2$krige.var, 
                                 mt.win.mr.uk2.3$krige.var, mt.win.mr.uk2.4$krige.var, 
                                 mt.win.mr.uk2.5$krige.var, mt.win.mr.uk2.6$krige.var, 
                                 mt.win.mr.uk2.7$krige.var, mt.win.mr.uk2.8$krige.var,
                                 mt.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_MT_Kriged_UK.tif',options=c('TFW=YES'))


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

## phi = 27
## relative nugget (nugget / (nugget + sill)) = (46389/(46389+38592)) = 0.55
mt.sum.mr.gls <- gls(MT_Sum ~ Easting + Northing + Veg_type, data=mt.sum.mr, correlation = corSpher(value=c(27, 0.55), form= ~ Easting + Northing, nugget = TRUE))
summary(mt.sum.mr.gls)

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

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

## compare the spatial and non-spatial models
mt.sum.mr.gls2 <- gls(MT_Sum ~ Veg_type, data=mt.sum.mr) ## non-spatial model
anova(mt.sum.mr.gls2, mt.sum.mr.gls)
## the covariance model is significant at p=0.003, 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.mr.gd <- as.geodata(mt.sum.mr[,c(5,4,8)])
mt.sum.mr.gd$data <- resid(mt.sum.mr.gls, resType="response")[1:195]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the MT_Sum and Veg_type columns
mt.sum.mr.uk2.1 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.2 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.3 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.4 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.5 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.6 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.7 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.8 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()
mt.sum.mr.uk2.9 <- krige(data=mt.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=mt.sum.mr.fit, trend.mod=mt.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$MT_Sum <- c(mt.sum.mr.uk2.1$predict, mt.sum.mr.uk2.2$predict, 
                        mt.sum.mr.uk2.3$predict, mt.sum.mr.uk2.4$predict, 
                        mt.sum.mr.uk2.5$predict, mt.sum.mr.uk2.6$predict, 
                        mt.sum.mr.uk2.7$predict, mt.sum.mr.uk2.8$predict,
                        mt.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$MT_Sum_std <- sqrt(c(mt.sum.mr.uk2.1$krige.var, mt.sum.mr.uk2.2$krige.var, 
                                 mt.sum.mr.uk2.3$krige.var, mt.sum.mr.uk2.4$krige.var, 
                                 mt.sum.mr.uk2.5$krige.var, mt.sum.mr.uk2.6$krige.var, 
                                 mt.sum.mr.uk2.7$krige.var, mt.sum.mr.uk2.8$krige.var,
                                 mt.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_MT_Kriged_UK.tif',options=c('TFW=YES'))


#### Coumarins Winter MR ####
co.win.mr <- data %>%
  filter(!is.na(Coum_Win), Site=="Camas") %>% 
  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.mr.gls <- gls(Coum_Win ~ Easting + Northing + Veg_type, data=co.win.mr, correlation = corSpher(value=c(33, 0.01), form= ~ Easting + Northing, nugget = TRUE))
summary(co.win.mr.gls)

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

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

## compare the spatial and non-spatial models
co.win.mr.gls2 <- gls(Coum_Win ~ Veg_type, data=co.win.mr) ## non-spatial model
anova(co.win.mr.gls2, co.win.mr.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.mr.gd <- as.geodata(co.win.mr[,c(5,4,8)])
co.win.mr.gd$data <- resid(co.win.mr.gls, resType="response")[1:206]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Coum_Win and Veg_type columns
co.win.mr.uk2.1 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.2 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.3 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.4 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.5 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.6 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.7 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.8 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()
co.win.mr.uk2.9 <- krige(data=co.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=co.win.mr.fit, trend.mod=co.win.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Coum_Win <- c(co.win.mr.uk2.1$predict, co.win.mr.uk2.2$predict, 
                        co.win.mr.uk2.3$predict, co.win.mr.uk2.4$predict, 
                        co.win.mr.uk2.5$predict, co.win.mr.uk2.6$predict, 
                        co.win.mr.uk2.7$predict, co.win.mr.uk2.8$predict,
                        co.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Coum_Win_std <- sqrt(c(co.win.mr.uk2.1$krige.var, co.win.mr.uk2.2$krige.var, 
                                 co.win.mr.uk2.3$krige.var, co.win.mr.uk2.4$krige.var, 
                                 co.win.mr.uk2.5$krige.var, co.win.mr.uk2.6$krige.var, 
                                 co.win.mr.uk2.7$krige.var, co.win.mr.uk2.8$krige.var,
                                 co.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_CO_Kriged_UK.tif',options=c('TFW=YES'))


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

## phi = 100
## relative nugget (nugget / (nugget + sill)) = (0.064/(0.064+0.041)) = 0.61
co.sum.mr.gls <- gls(Coum_Sum ~ Easting + Northing + Veg_type, data=co.sum.mr, correlation = corSpher(value=c(100, 0.61), form= ~ Easting + Northing, nugget = TRUE))
summary(co.sum.mr.gls) #p-values are close to 0.05, might be something to watch

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

## double check properties of the residuals to ensure that the assumptions are met
qqnorm(co.sum.mr.gls) # not the best
plot(co.sum.mr.gls, pch=16) ## some heteroscedasticity
plot(Variogram(co.sum.mr.gls, form= ~Easting + Northing, resType="response"), pch=16)
plot(Variogram(co.sum.mr.gls, form= ~Easting + Northing, resType="normalized"), pch=16) ## also a bit off (trails down at the end)

## compare the spatial and non-spatial models
co.sum.mr.gls2 <- gls(Coum_Sum ~ Veg_type, data=co.sum.mr) ## non-spatial model
anova(co.sum.mr.gls2, co.sum.mr.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.mr.gd <- as.geodata(co.sum.mr[,c(5,4,8)])
co.sum.mr.gd$data <- resid(co.sum.mr.gls, resType="response")[1:178]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Coum_Sum and Veg_type columns
co.sum.mr.uk2.1 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.2 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.3 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.4 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.5 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.6 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.7 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.8 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()
co.sum.mr.uk2.9 <- krige(data=co.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=co.sum.mr.fit, trend.mod=co.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Coum_Sum <- c(co.sum.mr.uk2.1$predict, co.sum.mr.uk2.2$predict, 
                        co.sum.mr.uk2.3$predict, co.sum.mr.uk2.4$predict, 
                        co.sum.mr.uk2.5$predict, co.sum.mr.uk2.6$predict, 
                        co.sum.mr.uk2.7$predict, co.sum.mr.uk2.8$predict,
                        co.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Coum_Sum_std <- sqrt(c(co.sum.mr.uk2.1$krige.var, co.sum.mr.uk2.2$krige.var, 
                                 co.sum.mr.uk2.3$krige.var, co.sum.mr.uk2.4$krige.var, 
                                 co.sum.mr.uk2.5$krige.var, co.sum.mr.uk2.6$krige.var, 
                                 co.sum.mr.uk2.7$krige.var, co.sum.mr.uk2.8$krige.var,
                                 co.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_CO_Kriged_UK.tif',options=c('TFW=YES'))


#### Cineole Winter MR ####
cin.win.mr <- data %>%
  filter(!is.na(Cine_Win), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Cine_Win)

## phi = 113
## relative nugget (nugget / (nugget + sill)) = (1.615/(1.615+1.124)) = 0.59
cin.win.mr.gls <- gls(Cine_Win ~ Easting + Northing + Veg_type, data=cin.win.mr, correlation = corSpher(value=c(50, 0.95), form= ~ Easting + Northing, nugget = TRUE))
summary(cin.win.mr.gls)

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

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

## compare the spatial and non-spatial models
cin.win.mr.gls2 <- gls(Cine_Win ~ Veg_type, data=cin.win.mr) ## non-spatial model
anova(cin.win.mr.gls2, cin.win.mr.gls)
## spatial model is worse (adds no value), based on a likelihood ratio test

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

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Cine_Win and Veg_type columns
cin.win.mr.uk2.1 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.2 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.3 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.4 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.5 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.6 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.7 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.8 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()
cin.win.mr.uk2.9 <- krige(data=cin.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cin.win.mr.fit, trend.mod=cin.win.mr.gls2)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Cine_Win <- c(cin.win.mr.uk2.1$predict, cin.win.mr.uk2.2$predict, 
                        cin.win.mr.uk2.3$predict, cin.win.mr.uk2.4$predict, 
                        cin.win.mr.uk2.5$predict, cin.win.mr.uk2.6$predict, 
                        cin.win.mr.uk2.7$predict, cin.win.mr.uk2.8$predict,
                        cin.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Cine_Win_std <- sqrt(c(cin.win.mr.uk2.1$krige.var, cin.win.mr.uk2.2$krige.var, 
                                 cin.win.mr.uk2.3$krige.var, cin.win.mr.uk2.4$krige.var, 
                                 cin.win.mr.uk2.5$krige.var, cin.win.mr.uk2.6$krige.var, 
                                 cin.win.mr.uk2.7$krige.var, cin.win.mr.uk2.8$krige.var,
                                 cin.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_CIN_Kriged_UK.tif',options=c('TFW=YES'))


#### Cineole Summer MR ####
cin.sum.mr <- data %>%
  filter(!is.na(Cine_Sum), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Cine_Sum)

## phi = 25
## relative nugget (nugget / (nugget + sill)) = (3522/(3522+2022)) = 0.64
cin.sum.mr.gls <- gls(Cine_Sum ~ Easting + Northing + Veg_type, data=cin.sum.mr, correlation = corSpher(value=c(25, 0.64), form= ~ Easting + Northing, nugget = TRUE))
summary(cin.sum.mr.gls)

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

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

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

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

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Cine_Sum and Veg_type columns
cin.sum.mr.uk2.1 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.2 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.3 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.4 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.5 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.6 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.7 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.8 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()
cin.sum.mr.uk2.9 <- krige(data=cin.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cin.sum.mr.fit, trend.mod=cin.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Cine_Sum <- c(cin.sum.mr.uk2.1$predict, cin.sum.mr.uk2.2$predict, 
                        cin.sum.mr.uk2.3$predict, cin.sum.mr.uk2.4$predict, 
                        cin.sum.mr.uk2.5$predict, cin.sum.mr.uk2.6$predict, 
                        cin.sum.mr.uk2.7$predict, cin.sum.mr.uk2.8$predict,
                        cin.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Cine_Sum_std <- sqrt(c(cin.sum.mr.uk2.1$krige.var, cin.sum.mr.uk2.2$krige.var, 
                                 cin.sum.mr.uk2.3$krige.var, cin.sum.mr.uk2.4$krige.var, 
                                 cin.sum.mr.uk2.5$krige.var, cin.sum.mr.uk2.6$krige.var, 
                                 cin.sum.mr.uk2.7$krige.var, cin.sum.mr.uk2.8$krige.var,
                                 cin.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_CIN_Kriged_UK.tif',options=c('TFW=YES'))


#### Camphor Winter MR ####
cam.win.mr <- data %>%
  filter(!is.na(Camphor_Win), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Camphor_Win)

## phi = 15
## relative nugget (nugget / (nugget + sill)) = (9554/(9554+7554)) = 0.56
cam.win.mr.gls <- gls(Camphor_Win ~ Easting + Northing + Veg_type, data=cam.win.mr, correlation = corSpher(value=c(15, 0.56), form= ~ Easting + Northing, nugget = TRUE))
summary(cam.win.mr.gls)

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

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

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

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

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Camphor_Win and Veg_type columns
cam.win.mr.uk2.1 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.2 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.3 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.4 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.5 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.6 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.7 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.8 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()
cam.win.mr.uk2.9 <- krige(data=cam.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cam.win.mr.fit, trend.mod=cam.win.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Camphor_Win <- c(cam.win.mr.uk2.1$predict, cam.win.mr.uk2.2$predict, 
                        cam.win.mr.uk2.3$predict, cam.win.mr.uk2.4$predict, 
                        cam.win.mr.uk2.5$predict, cam.win.mr.uk2.6$predict, 
                        cam.win.mr.uk2.7$predict, cam.win.mr.uk2.8$predict,
                        cam.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Camphor_Win_std <- sqrt(c(cam.win.mr.uk2.1$krige.var, cam.win.mr.uk2.2$krige.var, 
                                 cam.win.mr.uk2.3$krige.var, cam.win.mr.uk2.4$krige.var, 
                                 cam.win.mr.uk2.5$krige.var, cam.win.mr.uk2.6$krige.var, 
                                 cam.win.mr.uk2.7$krige.var, cam.win.mr.uk2.8$krige.var,
                                 cam.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_CAM_Kriged_UK.tif',options=c('TFW=YES'))


#### Camphor Summer MR ####
cam.sum.mr <- data %>%
  filter(!is.na(Camphor_Sum), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, Camphor_Sum)

## phi = 20
## relative nugget (nugget / (nugget + sill)) = (13247/(13247+9245)) = 0.59
cam.sum.mr.gls <- gls(Camphor_Sum ~ Easting + Northing + Veg_type, data=cam.sum.mr, correlation = corSpher(value=c(20, 0.59), form= ~ Easting + Northing, nugget = TRUE))
summary(cam.sum.mr.gls)

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

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

## compare the spatial and non-spatial models
cam.sum.mr.gls2 <- gls(Camphor_Sum ~ Veg_type, data=cam.sum.mr) ## non-spatial model
anova(cam.sum.mr.gls2, cam.sum.mr.gls)
## the covariance model is marginally significant at p=0.055, based on a likelihood ratio test

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

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the Camphor_Sum and Veg_type columns
cam.sum.mr.uk2.1 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.2 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.3 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.4 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.5 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.6 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.7 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.8 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()
cam.sum.mr.uk2.9 <- krige(data=cam.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cam.sum.mr.fit, trend.mod=cam.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$Camphor_Sum <- c(cam.sum.mr.uk2.1$predict, cam.sum.mr.uk2.2$predict, 
                        cam.sum.mr.uk2.3$predict, cam.sum.mr.uk2.4$predict, 
                        cam.sum.mr.uk2.5$predict, cam.sum.mr.uk2.6$predict, 
                        cam.sum.mr.uk2.7$predict, cam.sum.mr.uk2.8$predict,
                        cam.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$Camphor_Sum_std <- sqrt(c(cam.sum.mr.uk2.1$krige.var, cam.sum.mr.uk2.2$krige.var, 
                                 cam.sum.mr.uk2.3$krige.var, cam.sum.mr.uk2.4$krige.var, 
                                 cam.sum.mr.uk2.5$krige.var, cam.sum.mr.uk2.6$krige.var, 
                                 cam.sum.mr.uk2.7$krige.var, cam.sum.mr.uk2.8$krige.var,
                                 cam.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_CAM_Kriged_UK.tif',options=c('TFW=YES'))


#### Chem Div Winter MR ####
cd.win.mr <- data %>%
  filter(!is.na(ChemDiv_Win), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, ChemDiv_Win)

## phi = 50
## relative nugget (nugget / (nugget + sill)) = (0.039/(0.039+0.003)) = 0.93
cd.win.mr.gls <- gls(ChemDiv_Win ~ Easting + Northing + Veg_type, data=cd.win.mr, correlation = corSpher(value=c(50, 0.93), form= ~ Easting + Northing, nugget = TRUE))
summary(cd.win.mr.gls)

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

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

## compare the spatial and non-spatial models
cd.win.mr.gls2 <- gls(ChemDiv_Win ~ Veg_type, data=cd.win.mr) ## non-spatial model
anova(cd.win.mr.gls2, cd.win.mr.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
cd.win.mr.gd <- as.geodata(cd.win.mr[,c(5,4,8)])
cd.win.mr.gd$data <- resid(cd.win.mr.gls, resType="response")[1:205]

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the ChemDiv and Veg_type columns
cd.win.mr.uk2.1 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.2 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.3 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.4 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.5 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.6 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.7 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.8 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()
cd.win.mr.uk2.9 <- krige(data=cd.win.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cd.win.mr.fit, trend.mod=cd.win.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$ChemDiv_Win <- c(cd.win.mr.uk2.1$predict, cd.win.mr.uk2.2$predict, 
                        cd.win.mr.uk2.3$predict, cd.win.mr.uk2.4$predict, 
                        cd.win.mr.uk2.5$predict, cd.win.mr.uk2.6$predict, 
                        cd.win.mr.uk2.7$predict, cd.win.mr.uk2.8$predict,
                        cd.win.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$ChemDiv_Win_std <- sqrt(c(cd.win.mr.uk2.1$krige.var, cd.win.mr.uk2.2$krige.var, 
                                 cd.win.mr.uk2.3$krige.var, cd.win.mr.uk2.4$krige.var, 
                                 cd.win.mr.uk2.5$krige.var, cd.win.mr.uk2.6$krige.var, 
                                 cd.win.mr.uk2.7$krige.var, cd.win.mr.uk2.8$krige.var,
                                 cd.win.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Winter_CD_Kriged_UK.tif',options=c('TFW=YES'))


#### Chem Div Summer MR ####
cd.sum.mr <- data %>%
  filter(!is.na(ChemDiv_Sum), Site=="Camas") %>% 
  dplyr::select(GPS_Name, Veg_type, Site, Northing, Easting, Elev, uav_mean_height, ChemDiv_Sum)

## phi = 50
## relative nugget (nugget / (nugget + sill)) = (0.041/(0.041+0.003)) = 0.91
cd.sum.mr.gls <- gls(ChemDiv_Sum ~ Easting + Northing + Veg_type, data=cd.sum.mr, correlation = corSpher(value=c(10, 0.8), form= ~ Easting + Northing, nugget = TRUE))
summary(cd.sum.mr.gls)

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

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

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

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

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

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

## Krige with covariance model and trend model
# Make sure that the "data" has both the ChemDiv_Sum and Veg_type columns
cd.sum.mr.uk2.1 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.1, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.2 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.2, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.3 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.3, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.4 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.4, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.5 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.5, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.6 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.6, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.7 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.7, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.8 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.8, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()
cd.sum.mr.uk2.9 <- krige(data=cd.sum.mr[,c(5,4,8,2)], predict=mr.hab.grid2.9, cov.mod=cd.sum.mr.fit, trend.mod=cd.sum.mr.gls)
gc()

# Combine predictions back into the mr.hab.grid data frame
mr.hab.grid$ChemDiv_Sum <- c(cd.sum.mr.uk2.1$predict, cd.sum.mr.uk2.2$predict, 
                        cd.sum.mr.uk2.3$predict, cd.sum.mr.uk2.4$predict, 
                        cd.sum.mr.uk2.5$predict, cd.sum.mr.uk2.6$predict, 
                        cd.sum.mr.uk2.7$predict, cd.sum.mr.uk2.8$predict,
                        cd.sum.mr.uk2.9$predict)
# Add the square root of Krige variance, aka the standard deviation
# This gives a map version of uncertainty in the kriged prediction
mr.hab.grid$ChemDiv_Sum_std <- sqrt(c(cd.sum.mr.uk2.1$krige.var, cd.sum.mr.uk2.2$krige.var, 
                                 cd.sum.mr.uk2.3$krige.var, cd.sum.mr.uk2.4$krige.var, 
                                 cd.sum.mr.uk2.5$krige.var, cd.sum.mr.uk2.6$krige.var, 
                                 cd.sum.mr.uk2.7$krige.var, cd.sum.mr.uk2.8$krige.var,
                                 cd.sum.mr.uk2.9$krige.var))

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

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

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

writeRaster(mr.hab.export,'MR_Summer_CD_Kriged_UK.tif',options=c('TFW=YES'))

#### Cross-validation ####

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

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

pred <- cp.win.mr.gls$fitted + cv.cp.win.mr$predicted
actual <- cp.win.mr.gls$fitted + resid(cp.win.mr.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 MR
cv.cp.sum.mr <- xvalid(cp.sum.mr.gd, model=cp.sum.mr.fit) #LOOCV
summary(cv.cp.sum.mr)

pred <- cp.sum.mr.gls$fitted + cv.cp.sum.mr$predicted
actual <- cp.sum.mr.gls$fitted + resid(cp.sum.mr.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 MR
cv.mt.win.mr <- xvalid(mt.win.mr.gd, model=mt.win.mr.fit) #LOOCV
summary(cv.mt.win.mr)

pred <- mt.win.mr.gls$fitted + cv.mt.win.mr$predicted
actual <- mt.win.mr.gls$fitted + resid(mt.win.mr.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 MR
cv.mt.sum.mr <- xvalid(mt.sum.mr.gd, model=mt.sum.mr.fit) #LOOCV
summary(cv.mt.sum.mr)

pred <- mt.sum.mr.gls$fitted + cv.mt.sum.mr$predicted
actual <- mt.sum.mr.gls$fitted + resid(mt.sum.mr.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 MR
cv.co.win.mr <- xvalid(co.win.mr.gd, model=co.win.mr.fit) #LOOCV
summary(cv.co.win.mr)

pred <- co.win.mr.gls$fitted + cv.co.win.mr$predicted
actual <- co.win.mr.gls$fitted + resid(co.win.mr.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 Summer MR
cv.co.sum.mr <- xvalid(co.sum.mr.gd, model=co.sum.mr.fit) #LOOCV
summary(cv.co.sum.mr)

pred <- co.sum.mr.gls$fitted + cv.co.sum.mr$predicted
actual <- co.sum.mr.gls$fitted + resid(co.sum.mr.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

## CIN Winter MR
cv.cin.win.mr <- xvalid(cin.win.mr.gd, model=cin.win.mr.fit) #LOOCV
summary(cv.cin.win.mr)

pred <- cin.win.mr.gls$fitted + cv.cin.win.mr$predicted
actual <- cin.win.mr.gls$fitted + resid(cin.win.mr.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

## CIN Summer MR
cv.cin.sum.mr <- xvalid(cin.sum.mr.gd, model=cin.sum.mr.fit) #LOOCV
summary(cv.cin.sum.mr)

pred <- cin.sum.mr.gls$fitted + cv.cin.sum.mr$predicted
actual <- cin.sum.mr.gls$fitted + resid(cin.sum.mr.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

## CAM Winter MR
cv.cam.win.mr <- xvalid(cam.win.mr.gd, model=cam.win.mr.fit) #LOOCV
summary(cv.cam.win.mr)

pred <- cam.win.mr.gls$fitted + cv.cam.win.mr$predicted
actual <- cam.win.mr.gls$fitted + resid(cam.win.mr.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

## CAM Summer MR
cv.cam.sum.mr <- xvalid(cam.sum.mr.gd, model=cam.sum.mr.fit) #LOOCV
summary(cv.cam.sum.mr)

pred <- cam.sum.mr.gls$fitted + cv.cam.sum.mr$predicted
actual <- cam.sum.mr.gls$fitted + resid(cam.sum.mr.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

## CD Winter MR
cv.cd.win.mr <- xvalid(cd.win.mr.gd, model=cd.win.mr.fit) #LOOCV
summary(cv.cd.win.mr)

pred <- cd.win.mr.gls$fitted + cv.cd.win.mr$predicted
actual <- cd.win.mr.gls$fitted + resid(cd.win.mr.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

## CD Summer MR
cv.cd.sum.mr <- xvalid(cd.sum.mr.gd, model=cd.sum.mr.fit) #LOOCV
summary(cv.cd.sum.mr)

pred <- cd.sum.mr.gls$fitted + cv.cd.sum.mr$predicted
actual <- cd.sum.mr.gls$fitted + resid(cd.sum.mr.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.mr$CP_Win ~ cp.win.mr.gls$fitted))
summary(lm(cp.sum.mr$CP_Sum ~ cp.sum.mr.gls$fitted))
summary(lm(mt.win.mr$MT_Win ~ mt.win.mr.gls$fitted))
summary(lm(mt.sum.mr$MT_Sum ~ mt.sum.mr.gls$fitted))
summary(lm(co.win.mr$Coum_Win ~ co.win.mr.gls$fitted))
summary(lm(co.sum.mr$Coum_Sum ~ co.sum.mr.gls$fitted))
summary(lm(cin.win.mr$Cine_Win ~ cin.win.mr.gls2$fitted))
summary(lm(cin.sum.mr$Cine_Sum ~ cin.sum.mr.gls$fitted))
summary(lm(cam.win.mr$Camphor_Win ~ cam.win.mr.gls$fitted))
summary(lm(cam.sum.mr$Camphor_Sum ~ cam.sum.mr.gls$fitted))
summary(lm(cd.win.mr$ChemDiv_Win ~ cd.win.mr.gls$fitted))
summary(lm(cd.sum.mr$ChemDiv_Sum ~ cd.sum.mr.gls$fitted))

# ## load the kriging wrapper function if haven't loaded entire KrigeWrapper3 script
# krige <- function(data, predict, cov.mod, trend.mod="cte"){    # default to Ordinary Kriging
#   data <- cbind(data,data[,1:2])  				# rearrange data to put coordinates in covariates
#   data.gd <- as.geodata(data,covar.col=4:dim(data)[2], covar.names=colnames(data[,4:dim(data)[2]]))  # put into geodata object
#   predict <- cbind(predict[,1:2],data=rep(NA,dim(predict)[1]),predict)  # rearrange predict to put NA in $data and coordinates in covariates
#   predict.gd <- as.geodata(predict, covar.col=4:dim(predict)[2], covar.names=colnames(predict[,4:dim(predict)[2]]), na.action="none")    # put into geodata object
#   
#   if(is.character(trend.mod)) {  # cte, 1st, 2nd for default OK and UK options
#     trend.d <- trend.mod
#     trend.l <- trend.mod
#   }
#   else  {
#     trend.lm <- as.formula(paste("~",as.character(formula(trend.mod))[3]))  # right side of trend formula
#     trend.d <- trend.spatial(trend=trend.lm, geodata=data.gd)		# trend for locations where we have data
#     trend.l <- trend.spatial(trend=trend.lm, geodata=predict.gd)	# trend for prediction locations
#   }
#   krige.mod <- krige.control(cov.model=cov.mod$cov.model, cov.pars=cov.mod$cov.pars, nugget=cov.mod$nugget, trend.d=trend.d, trend.l=trend.l)
#   krige.out<-krige.conv(data.gd, loc=predict.gd$coords, krige=krige.mod)
#   print(krige.out$beta.est)
#   return(krige.out)
# }