## R: Using matrix algebra to compute a Kriged value

**The following example was created in R and was inspired from O’Sullivan and Unwin’s Geographic Information Analysis (1st edition)**

**Data for this exercise can be downloaded here.**

library(maptools) # S = readShapePoints("points_krig.shp") # Convert data to spatstat format SP = as(S, "SpatialPoints")# This only grabs the X,Y values SPz = slot(S,"data") # This grabs all attribute values P = as(SP, "ppp") # This creates a spatstat point pattern marks(P) = SPz # This adds the attribute values to P # thus creating a marked point pattern hist(marks(P)$z) # plots histogram of attribute z

# Convert data to dataframe P.frame = as.data.frame(P) P.frame.sub = data.frame(x=P.frame$x,y=P.frame$y,z=P.frame$z) # Get individual columns X = P.frame$x Y = P.frame$y Z = P.frame$z # Now sort the vectors to match the order used in O-Sullivan & Unwinâ€™s Book X = rev(X) Y = rev(Y) Z = rev(Z) # Create distance matrix A = dist(cbind(X,Y),method="euclidean",diag=T) # Replace all distance values by the semivariance at the distance P.gamma = function(X,range,sill) {P.gamma=sill*(3*X/(2*range)-0.5* (X/range)^3);P.gamma} A = P.gamma(A,95,100) # Augment matrix by one row and one column augCol = matrix(rep(1,7),nrow=7,ncol=1) augRow = matrix(c(rep(1,7),0),1,8) A = cbind(as.matrix(A),augCol) A = rbind(as.matrix(A),augRow) # Inverse matrix A A.inv = solve(A) # Create distance matrix (single column) between unknown value and all control points # Unknown value is located at 50,50 in this example X2 = c(50,X) Y2 = c(50,Y) b = as.matrix(as.matrix(dist(cbind(X2,Y2),method="euclidean",diag=T))[,1]) # We only need dist for unk. # Replace distance values by the semivariance b = P.gamma(b,95,100) # Now remove first record and augment by 1 b = as.matrix(rbind(as.matrix(b[-1,]),1)) # Finally, compute the weights w = A.inv %*% b # note that all but the last row sum to 1 # The unknown Z value at 50,50 can now be computed from these weights Z.est = sum(Z * w[-8,]) # equals 10.36 # The estimated variation is computed from est.var = sum(w*b) # The standard error is std.err = est.var^0.5 # equals 5.696 # To estimate a 95% confidence interval, set limits a 1.96 times std error upper.range = 1.96 * std.err + Z.est # = 21.53 lower.range = 1.96 * std.err - Z.est # = 0.80