---
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: