set.seed(456)

Read in Data:

Oregon boundaries

orotl <- st_read("/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/oregon/orotl.shp", quiet = TRUE)
st_crs(orotl) <- 4326

Oregon temperatures

ortann <- st_read("/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/oregon/oregontann.shp", quiet = TRUE)
st_crs(ortann) <- 4326

Oregon DEM file

orgrid <- st_read("/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/oregon/orgrid.shp", quiet = TRUE) 
st_crs(orgrid) <- 4326
orgrid.dem <- st_rasterize(orgrid, dx = 0.1667, dy = 0.1667)
st_crs(orgrid.dem) <- 4326

Plot

plot(orgrid.dem, reset = FALSE)

plot(ortann["tann"],
     add = TRUE, 
     pch = 16, 
     cex = 1.5,
      legend = TRUE)
plot(orotl,
    add = TRUE, 
     pch = 16, 
     cex = 1.5,
    col = NA,
    fill = NA,
    border = "red")

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

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

Create a variogram model for this data using the vgm() function.

ortann.vgm1 <- vgm(psill = 4, 
                  model = "Cir", 
                  range = 140, 
                  nugget = 0.1)

I have used the values of 4 for sill, 140 for range, and 0.1 for nugget. I used the CIR model as well.

plot(ortann.var, ortann.vgm1, main = "Ortann variogram")

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

ortann.vgm2 = fit.variogram(ortann.var, ortann.vgm1)
plot(ortann.var, ortann.vgm2, main = "Oregon Average Temperature Variogram")

Now use this model to interpolate the annual temperatures using the grid from the DEM, using the krige() function.

ortann.pred <- krige(tann ~ 1, 
                     locations = ortann, 
                     newdata = orgrid.dem, 
                     model = ortann.vgm2, 
                     nmax = 40)
## [using ordinary kriging]
names(ortann.pred)
## [1] "var1.pred" "var1.var"

Produce a map showing the predicted value on the grid and the prediction error.

my.pal <- brewer.pal(9, "Reds")
plot(ortann.pred, col = my.pal, main = "Interpolated Temperature Values")
plot(orotl,
    add = TRUE, 
     pch = 16, 
     cex = 1.5,
    col = NA,
    fill = NA,
    border = "black")

plot(ortann.pred["var1.var"], col = my.pal,
     main = "Interpolated Temperature Prediction Error")

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

ortann.pred2 = krige.cv(tann ~ 1, 
                      locations = ortann, 
                      model = ortann.vgm2, 
                      nmax = 40, 
                      nfold = 5)

head(ortann.pred2)
## 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
## Geodetic CRS:  WGS 84
##   var1.pred    var1.var observed    residual      zscore fold
## 1  8.272909 0.012748689      9.6  1.32709126  11.7535221    5
## 2 10.643671 0.016502693     12.5  1.85632898  14.4503225    5
## 3 10.908796 0.014688052     11.1  0.19120376   1.5776633    3
## 4 10.211644 0.027474641     10.3  0.08835572   0.5330509    5
## 5  8.395192 0.004590634      7.6 -0.79519183 -11.7364176    4
## 6  7.777544 0.004672582      8.4  0.62245567   9.1060489    3
##                  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)

RMSEP

sqrt(mean(ortann.pred2$residual^2))
## [1] 1.220688
sp::bubble(as_Spatial(ortann.pred2)[,"residual"], pch = 16)

R2P

cor(ortann.pred2$observed, ortann.pred2$var1.pred)^2
## [1] 0.5727858
plot(ortann.pred2$var1.pred, ortann.pred2$residual, 
     xlab = 'Predicted Values', 
     ylab = 'Residuals')

abline(h = 0, lty = 2)