## ---- include=FALSE----------------------------------------------------------- library(knitr) opts_chunk$set(fig.width = 5, fig.height = 5, fig.cap = '', collapse = TRUE) library(rgdal) library(raster) library(gstat) library(fields) ## ---- krig0------------------------------------------------------------------- if (!require("rspatial")) remotes::install_github('rspatial/rspatial') ## ---- krig1------------------------------------------------------------------- library(rspatial) a <- sp_data('alberta.csv') m <- lm(z ~ x + y, data=a) summary(m) plot(a[,2:3], xlim=c(0,60), ylim=c(0,60), las=1, pch=20, yaxs="i", xaxs="i") text(a[,2:3], labels=a$z, pos=4) # make the contour lines x <- seq(0, 60, 1) y <- seq(0, 60, 1) # all combinations of x and y xy <- data.frame(expand.grid(x=x,y=y)) z <- predict(m, xy) z <- matrix(z, 61, 61) contour(x, y, z, add=TRUE, labcex=1.25) ## ----------------------------------------------------------------------------- library(raster) dp <- pointDistance(a[,2:3], lonlat=FALSE) dim(a) dim(dp) dp[1:5, 1:5] diag(dp) <- NA ## ----------------------------------------------------------------------------- dist(a$z[1:3]) dz <- dist(a$z) ## ----------------------------------------------------------------------------- dp <- as.dist(dp) ## ---- krig5------------------------------------------------------------------- # semivariance semivar <- dz^2 / 2 plot(dp, semivar, xlim=c(0, 80), ylim=c(0,220), xlab=c('Distance between locations'), ylab=c('Semivariance'), pch=20, axes=FALSE, xaxs="i") axis(1, at=seq(0,80,10)) axis(2, las=1) ## ---- krig10------------------------------------------------------------------ # choose a bin width (in spatial distance) binwidth <- 8 # assign a lag (bin number) to each record lag <- floor(dp/binwidth) + 1 # average value for each lag lsv <- tapply(semivar, lag, mean) # compute the average distance for each lag dlag <- tapply(dp, lag, mean) plot(dlag, lsv, pch=20, axes=FALSE, xlab='Distance', ylab='Semivariance', xlim=c(0,80)) axis(1, at=seq(0,80,10)) axis(2, las=1)