---
title: "Assigment: Spatial Autocorrelation"
author: "(Your Name and NetId Here)"
date: "`r Sys.Date()`"
output:
   html_document:
     toc: TRUE
     toc_float: TRUE
     theme: "lumen"
     code_folding: show
     number_section: TRUE
     fig_height: 4
     fig_width: 6
     fig_align: 'center'
---

```{r setup, include=FALSE} 
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

This exercise continues the "Trend Surfaces" assignment. **Read the instructions and notes carefully**.

Task: load the `sf`, `gstat`, `terra` and `nlme` packages.

Task: restore the workspace objects from the previous exercise `ex_TrendSurface_Assignment.Rmd`, using the `load()` function with the `verbose=TRUE` option to see the objects that are loaded. 

Two choices for this:

1. Load the `.Rdata` file you saved as the last step in that assignment.

2. Use the `.Rdata` file from the model answers of that assignment. This is provided with the Assignment on Canvas, and is named `ex_TrendSurface_Assignment_Results.RData`. If you select this option, the object names will be the ones from the model answers. So you should look there for their meaning.

(Note: The reason to save from one file and load here is because when you `knit` an R Markdown (`.Rmd`) file it starts with a clean workspace.)

Task: Add a gridded DEM covering the study area to the workspace, from file `dem_nj_ny_pa_vt_4km.tif`,  using the `rast` method. Name the object `dem.ne.m` (standing for "DEM of the Northeast, elevation in metres").

Note you did this in the previous exercise, but `terra` data structures are not well-preserved with conventional R `save` and `load`, so you need to re-load the original GeoTIFF and convert to a `SpatRaster`.


This assignment follows the methods of "Tutorial: Trend surfaces in R", $\S9$ "Spatial Correlation of trend surface residuals" --$\S12$ "GLS-Regression Kriging", with the difference that here we use the prepared `SpatRaster` object `dem.ne.m` as the prediction grid, whereas in the Tutorial we constructed a grid.

# Spatial correlation of trend surface residuals

Task: Compute the empirical variogram of the 2nd-order trend surface residuals, with the default cutoff.



Q: What is the approximate range of *local* spatial positive autocorrelation?

A: 

Task: Re-compute the empirical variogram of the 2nd-order trend surface residuals, with this range  plus about 10% (to clearly see the sill) as the cutoff (using the `cutoff` optional argument to the `variogram` function), and plot it.



Q:  What are the estimated sill, range, and nugget of this variogram?

A: 

# Generalized Least Squares trend surface

Task: Compute the coefficients of a full second-order trend, using GLS. Initialize the `gls` correlation structure with the estimated parameters from the previous example.

If you identified a nugget variance, specify `nugget=TRUE` in the argument to the `corExp` function inside the `gls` function. This must  be initialized as a **proportion** of the total sill. So the `value` argument must be a **vector** (constructed with the `c` function) of the range parameter (recall: for an exponential model, 1/3 the effective range) and the proportional nugget.



Q: What are the range of spatial correlation and the proportional nugget variance of the exponential model, as estimated by `gls`? Recalling that the range *parameter* (as fit by `gls` or `fit.variogram`) in an exponential model is 1/3 of the effective range of spatial dependence, what is that? How do these compare with the estimates from the OLS trend?


A: 


Task: Compare the coefficients from the GLS and OLS fits (you computed this in the previous exercise), as _absolute_ differences and as _percentages_ of the OLS fit.


Q: How much change is there between the OLS and GLS trend surface coefficients?

A: 

Task: Predict over the grid with the GLS trend.. Compare the predictions here with those from the OLS 2nd-order trend.


Task: Display this map, alongside the OLS trend surface map, with the same colour ramp.


Q: Are there obvious differences in these?

A: 

Task: Compute the difference between the OLS and GLS trend surfaces, and map them. 



Q: Where are the largest differences between the OLS and GLS trend surfaces? Explain why.

A: 

# Ordinary Kriging (OK) interpolation of the residuals

Task: Display the residuals from the GLS trend surface as a post- plot.



Q: Do these appear to be spatially autocorrelated?

A: 

Task: Compute the empirical variogram model residuals from the GLS trend surface model, and model it with an exponential variogram function.



Q: What are the parameters of the fitted variogram? How well do these match the  correlation structure fitted by `gls`?

A: 

Task: Predict the residuals over the grid by OK,using the fitted variogram model.



Task: Display the OK predictions.



Q: Where are the biggest adjustments? Try to explain why.

A: 

Task: Display the OK prediction standard deviations.



Q: Where are the largest and smallest prediction standard deviations? Why?

A: 

  
# GLS-Regression Kriging

Task: Add the OK predictions of the GLS residuals to the prediction grid object, and then add this to the GLS trend surface prediction to obtain a final prediction.



Q: Compare the summary of the GLS-RK predictions with those from the GLS trend surface. Which has a wider range? Which has a higher median and mean?



A: 
