R: Using matrix algebra to compute a Kriged value

15 October, 2010 (14:04) | R, Spatial analysis | By: Manuel Gimond

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