Lab12_Exercises

library(ggplot2)

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

library(gstat)

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

library(RColorBrewer)
library(sf)

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

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

library(stars)

## Warning: package 'stars' was built under R version 4.1.2

## Loading required package: abind

library(viridis)

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

## Loading required package: viridisLite

# Oregon Boundaries
orotl <- st_read("C:/Users/taylu/OneDrive/Documents/Geog6000/Datafiles/oregon/oregon/orotl.shp", quiet = TRUE)
st_crs(orotl) <-4326

# Oregon Temperatures
ortann <- st_read("C:/users/taylu/OneDrive/Documents/Geog6000/Datafiles/oregon/oregon/oregontann.shp", quiet = TRUE)
st_crs(ortann) <-4326

# Oregon DEM File
orgrid <- st_read("C:/users/taylu/OneDrive/Documents/Geog6000/Datafiles/oregon/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(orgrid.dem, reset = FALSE)

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

Create a Sample Variogram

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

Create Variogram Model

ortann.vgml <- vgm(psill = 3.5, "Cir", range = 190, nugget = 0.5)
plot(ortann.var, ortann.vgml)

Fit the Model to Sample Variogram

ortann.vgml2 = fit.variogram(ortann.var, ortann.vgml)
plot(ortann.var, ortann.vgml2)

Final Fitted Parameters of Variogram

ortann.vgml2

##   model     psill    range
## 1   Nug 0.2076555   0.0000
## 2   Cir 3.7168890 208.1694

Interpolate the Annual Temperature

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

## [using ordinary kriging]

Predicted Value Map

my.pal <- brewer.pal(8, "Reds")
plot(ortann.pred.ok, col = my.pal, main = "Interpolated Temperature Values")

Predicted Error Map

my.pal1 <- brewer.pal(8, "Blues")
plot(ortann.pred.ok["var1.var"], col = my.pal1, 
     main = "Interpolated Temperature Prediction Error")

Krige.cv() Function with 5 Fold Cross Validation

ortann.cv.ok <- krige.cv(tann ~ 1,
                         locations = ortann,
                         model = ortann.vgml2,
                         nmax = 40,
                         nfold = 5)
head(ortann.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
## Geodetic CRS:  WGS 84
##   var1.pred  var1.var observed   residual     zscore fold
## 1  8.447431 1.0119609      9.6  1.1525694  1.1457378    1
## 2 10.728969 1.2915011     12.5  1.7710306  1.5583990    3
## 3 10.470014 1.5908109     11.1  0.6299860  0.4994841    1
## 4 10.813039 1.0373384     10.3 -0.5130391 -0.5037212    4
## 5  8.629794 0.5698257      7.6 -1.0297938 -1.3642037    5
## 6  7.914157 0.5911681      8.4  0.4858425  0.6318879    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 could be expected when making a prediction

sqrt(mean(ortann.cv.ok$residual^2))

## [1] 1.15334

R2

cor(ortann.cv.ok$observed, ortann.cv.ok$var1.pred)^2

## [1] 0.6172486

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

plot(ortann.cv.ok$var1.pred, ortann.cv.ok$residual,
     xlab = 'Tann Predicted Values',
     ylab = 'Tann Residuals')
abline(h = 0, lty = 2)