Lab 12

Set working directory and read files

setwd("/Users/annapeterson/Desktop/Classes/GEOG6000/Lab12")

#Code from Exercise instructions
# Oregon boundaries
orotl <- st_read("/Users/annapeterson/Desktop/Classes/GEOG6000/Lab12/orotl.shp", 
                 quiet = TRUE)
st_crs(orotl) <- 4326

# Oregon temperatures
ortann <- st_read("/Users/annapeterson/Desktop/Classes/GEOG6000/Lab12/oregontann.shp", 
                  quiet = TRUE)
st_crs(ortann) <- 4326

# Oregon DEM file
orgrid <- st_read("/Users/annapeterson/Desktop/Classes/GEOG6000/Lab12/orgrid.shp", 
                  quiet = TRUE) 
st_crs(orgrid) <- 4326
orgrid.dem <- st_rasterize(orgrid, dx = 0.1667, dy = 0.1667)

Create a Sample Variogram

First check the distribution by using a histogram and log-transforming the data

Keeping the regular data because the log-transform does not make it more normally distributed. (If I understand this properly).

Create Variogram Model

Fit Model to Sample Variogram

Kriging

pred_ortann = krige(tann ~ 1, 
                    locations = ortann,
                    newdata = orgrid,
                    model = vgm2_ortann,
                    nmax = 40)

## [using ordinary kriging]

#I could not, for the life of me, get any of the viridis and colorbrewer palettes to work :(
plot(pred_ortann["var1.pred"],
     pal = NULL,
     nbreaks = 11,
     breaks = "pretty",
     pch = 15,
     cex = 1.5,
     main = "Interpolated Values")

plot(pred_ortann["var1.var"],
     pal = NULL,
     nbreaks = 11,
     breaks = "pretty",
     pch = 15,
     cex = 1.5 ,
     main = "Interpolated Prediction Error")

Krige.cv Function with 5-fold cross-validation

ortann_cv_ok = krige.cv(tann ~ 1,
                        locations = ortann,
                        model = vgm2_ortann,
                        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.446901 1.1103573      9.6  1.1530986  1.0942967    1
## 2 10.453338 1.4121705     12.5  2.0466623  1.7222755    1
## 3 11.367931 1.0823681     11.1 -0.2679313 -0.2575349    3
## 4 10.718828 1.0404264     10.3 -0.4188275 -0.4106100    3
## 5  8.511057 0.5676520      7.6 -0.9110572 -1.2092179    5
## 6  8.163135 0.5937375      8.4  0.2368654  0.3074003    4
##                  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_cv_ok$residual^2))

## [1] 1.273185

R²

cor(ortann_cv_ok$observed, ortann_cv_ok$var1.pred)^2

## [1] 0.536004

Using bubbles to see if there’s any pattern to the residual