lab11

setwd("~/Documents/geog6000/lab11")

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.6.2

library(gstat)

## Warning: package 'gstat' was built under R version 3.6.2

library(RColorBrewer)
library(sf)

## Warning: package 'sf' was built under R version 3.6.2

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(stars)

## Loading required package: abind

library(viridis)

## Warning: package 'viridis' was built under R version 3.6.2

## Loading required package: viridisLite

## Warning: package 'viridisLite' was built under R version 3.6.2

Read in and convert the files

# Oregon boundaries
orotl <- st_read("../datafiles/oregon/orotl.shp", quiet = TRUE)
st_crs(orotl) <- 4326

# Oregon temperatures
ortann <- st_read("../datafiles/oregon/oregontann.shp", quiet = TRUE)
st_crs(ortann) <- 4326

# Oregon DEM file
orgrid <- st_read("../datafiles/oregon/orgrid.shp", quiet = TRUE) 
st_crs(orgrid) <- 4326
orgrid.dem <- st_rasterize(orgrid, dx = 0.1667, dy = 0.1667)

# plot
plot(orgrid.dem, reset = FALSE)

plot(ortann["tann"], 
     add = TRUE, 
     pch = 16, 
     cex = 1.5)

Produce a sample variogram for average annual temperatures in Oregon. (Use the variogram() function)

tann.var <- variogram(tann ~ 1, data = ortann)
plot(tann.var, plot.numbers = TRUE, pch = '+')

tann.var2 <- variogram(tann ~ 1, 
                        data = ortann)

plot(tann.var2, plot.numbers = TRUE, pch = '+')

Create a variogram model for this data using the vgm() function. You will need to choose an appropriate model and initial parameters for the nugget, sill and range. Report the values and model you have used.

modNugget <- -0.1
modRange <- 140
modSill <- 4

tann.vgm1 <- vgm(psill = modSill, 
                  model = "Cir", 
                  range = modRange, 
                  nugget = modNugget)

plot(tann.var, tann.vgm1, main = "Oregon Ltann Variogram")

Use the fit.variogram() function to fit this model to the sample variogram from step a. Produce a plot showing the final variogram.

tann.vgm2 <- fit.variogram(tann.var, tann.vgm1)

plot(tann.var, tann.vgm2, main = "Oregon Ltann Variogram")

tann.vgm2

##   model     psill    range
## 1   Nug 0.2076608   0.0000
## 2   Cir 3.7169254 208.1733

Now use this model to interpolate the annual temperatures using the grid from the DEM, using the krige() function. Produce a map showing the predicted value on the grid and the prediction error.

##ok = ordinary kriging

tann.pred.ok <- krige(tann ~ 1, 
                     locations = ortann, 
                     newdata = orgrid.dem, 
                     model = tann.vgm2, 
                     nmax = 40)

## [using ordinary kriging]

Predicted value map:

my.pal <- brewer.pal(8, "Reds")
plot(tann.pred.ok, col = my.pal, main = "Interpolated temperature values (OK)")

Predicted error map:

my.pal1 <- brewer.pal(8, "Oranges")
plot(tann.pred.ok["var1.var"], col = my.pal1,
     main = "Interpolated temperature prediction error (OK)")

Use the krige.cv() function with 5-fold cross-validation to report the root mean squared error and R2

##cv = cross validation
tann.cv.ok <- krige.cv(tann ~ 1, 
                      locations = ortann, 
                      model = tann.vgm2, 
                      nmax = 40, 
                      nfold = 5)
head(tann.cv.ok)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -123.883 ymin: 42.217 xmax: -117.817 ymax: 46.15
## CRS:           +proj=longlat +datum=WGS84 +no_defs
##   var1.pred  var1.var observed   residual     zscore fold
## 1  8.430224 0.9824586      9.6  1.1697757  1.1801724    2
## 2 10.958861 1.2232092     12.5  1.5411390  1.3934502    1
## 3 11.303051 1.0848295     11.1 -0.2030513 -0.1949508    5
## 4 10.940677 1.1551003     10.3 -0.6406767 -0.5961137    5
## 5  8.519641 0.5690507      7.6 -0.9196409 -1.2191097    3
## 6  8.219295 0.5925382      8.4  0.1807049  0.2347533    2
##                  geometry
## 1 POINT (-120.717 44.917)
## 2   POINT (-120.2 45.717)
## 3 POINT (-122.717 42.217)
## 4  POINT (-123.883 46.15)
## 5 POINT (-117.817 44.833)
## 6 POINT (-117.833 44.783)

Root mean squared error

##average error that might be expected when making a prediction
sqrt(mean(tann.cv.ok$residual^2))

## [1] 1.274448

R2

##amount of variance in the test dataset predicted by the model
cor(tann.cv.ok$observed, tann.cv.ok$var1.pred)^2

## [1] 0.535641

sp::bubble(as_Spatial(tann.cv.ok)[,"residual"], pch = 16)

plot(tann.cv.ok$var1.pred, tann.cv.ok$residual, 
     xlab = 'Tann Predicted Values', 
     ylab = 'Tann Residuals')

abline(h = 0, lty = 2)