Read in data
# Oregon boundaries
orotl <- st_read("orotl.shp", quiet = TRUE)
st_crs(orotl) <- 4326
orotl<- readOGR("orotl.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/jessieshapiro/Documents/Documents/Education/University of Utah/GEOG6000/Labs/lab10/oregon/orotl.shp", layer: "orotl"
## with 36 features
## It has 1 fields
# Oregon temperatures
ortann <- st_read("oregontann.shp", quiet = TRUE)
st_crs(ortann) <- 4326
# Oregon DEM file
orgrid <- st_read("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)
Create a sample variogram
ortann.var = variogram(tann ~ 1, ortann)
#ortann.var = variogram(log(tann) ~ 1, ortann)
plot(ortann.var, plot.numbers=TRUE)
Create a varigoram model
ortann.vgm1 = vgm(psill= 3.5, "Cir", range= 190, nugget=0.5)
plot(ortann.var, ortann.vgm1)
Fit model to sample varigoram
ortann.vgm2 = fit.variogram(ortann.var, ortann.vgm1)
plot(ortann.var, ortann.vgm2)
Final fitted parameters of the variogram
ortann.vgm2
## model psill range
## 1 Nug 0.2076555 0.0000
## 2 Cir 3.7168890 208.1694
Interpolate annual temperature
ppt.pred.ok <- krige(ortann$tann ~ 1,
locations = ortann,
newdata = orgrid.dem,
model = ortann.vgm2,
nmax = 40)
## [using ordinary kriging]
names(ppt.pred.ok)
## [1] "var1.pred" "var1.var"
nclr = 9
plotclr <- brewer.pal(nclr, "Blues")
plot(ppt.pred.ok["var1.pred"], col = plotclr, main = "Oregon Temperature Data (Ordinary Kriging)", reset = FALSE)
plot(orotl, add = TRUE)
Run 5-fold cross-validation
ppt.cv.ok = krige.cv(tann ~ 1, ortann,
ortann.vgm2,
nmax=40,
nfold=5)
head(ppt.cv.ok)
## Simple feature collection with 6 features and 6 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -123.883 ymin: 42.217 xmax: -117.817 ymax: 46.15
## CRS: NA
## var1.pred var1.var observed residual zscore fold
## 1 8.390349 1.1010059 9.6 1.2096507 1.1528297 2
## 2 10.852343 1.2288801 12.5 1.6476567 1.4863189 2
## 3 11.433005 1.0836932 11.1 -0.3330049 -0.3198876 4
## 4 10.743165 1.0352335 10.3 -0.4431648 -0.4355581 1
## 5 8.504996 0.5688224 7.6 -0.9049961 -1.1999368 4
## 6 8.081500 0.5972574 8.4 0.3184998 0.4121244 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(ppt.cv.ok$residual^2))
## [1] 1.234224
R-Squared
cor(ppt.cv.ok$observed, ppt.cv.ok$var1.pred)^2
## [1] 0.5619103