{"id":10,"date":"2010-10-15T14:04:24","date_gmt":"2010-10-15T18:04:24","guid":{"rendered":"http:\/\/web.colby.edu\/mgimond\/?p=10"},"modified":"2011-10-26T13:58:30","modified_gmt":"2011-10-26T17:58:30","slug":"test","status":"publish","type":"post","link":"https:\/\/web.colby.edu\/mgimond\/2010\/10\/15\/test\/","title":{"rendered":"R: Using matrix algebra to compute a Kriged value"},"content":{"rendered":"

The following example was created in R and was inspired from O’Sullivan and Unwin’s Geographic Information Analysis (1st edition)<\/span><\/strong><\/p>\n

Data for this exercise can be downloaded here<\/a>.<\/span><\/strong><\/p>\n

\n
\r\nlibrary(maptools)\r\n#\r\nS = readShapePoints("points_krig.shp")\r\n# Convert data to spatstat format\r\nSP  = as(S, "SpatialPoints")# This only grabs the X,Y values\r\nSPz = slot(S,"data") # This grabs all attribute values\r\nP   = as(SP, "ppp")  # This creates a spatstat point pattern\r\nmarks(P) = SPz       # This adds the attribute values to P\r\n                     # thus creating a marked point pattern\r\nhist(marks(P)$z)     # plots histogram of attribute z\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Convert data to dataframe\r\nP.frame     = as.data.frame(P)\r\nP.frame.sub = data.frame(x=P.frame$x,y=P.frame$y,z=P.frame$z)\r\n\r\n# Get individual columns\r\nX = P.frame$x\r\nY = P.frame$y\r\nZ = P.frame$z\r\n\r\n# Now sort the vectors to match the order used in O-Sullivan & Unwin\u2019s Book\r\nX = rev(X)\r\nY = rev(Y)\r\nZ = rev(Z)\r\n\r\n# Create distance matrix\r\nA = dist(cbind(X,Y),method="euclidean",diag=T)\r\n\r\n# Replace all distance values by the semivariance at the distance\r\nP.gamma = function(X,range,sill) {P.gamma=sill*(3*X\/(2*range)-0.5* \r\n                                  (X\/range)^3);P.gamma}\r\nA = P.gamma(A,95,100)\r\n\r\n# Augment matrix by one row and one column\r\naugCol = matrix(rep(1,7),nrow=7,ncol=1)\r\naugRow = matrix(c(rep(1,7),0),1,8)\r\nA      = cbind(as.matrix(A),augCol)\r\nA      = rbind(as.matrix(A),augRow)\r\n\r\n# Inverse matrix A\r\nA.inv = solve(A)\r\n# Create distance matrix (single column) between unknown value and all control points\r\n# Unknown value is located at 50,50 in this example\r\nX2 = c(50,X)\r\nY2 = c(50,Y)\r\nb  = as.matrix(as.matrix(dist(cbind(X2,Y2),method="euclidean",diag=T))[,1]) # We only need dist for unk.\r\n# Replace distance values by the semivariance\r\nb = P.gamma(b,95,100)\r\n# Now remove first record and augment by 1\r\nb = as.matrix(rbind(as.matrix(b[-1,]),1))\r\n# Finally, compute the weights\r\nw = A.inv %*% b   # note that all but the last row sum to 1\r\n# The unknown Z value at 50,50 can now be computed from these weights\r\nZ.est = sum(Z * w[-8,]) # equals 10.36\r\n \r\n# The estimated variation is computed from\r\nest.var = sum(w*b)\r\n\r\n# The standard error is\r\nstd.err = est.var^0.5 # equals 5.696\r\n\r\n# To estimate a 95% confidence interval, set limits a 1.96 times std error\r\nupper.range = 1.96 * std.err + Z.est # = 21.53\r\nlower.range = 1.96 * std.err - Z.est # = 0.80\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
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.<\/p>\n","protected":false},"author":1199,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"ngg_post_thumbnail":0,"footnotes":""},"categories":[12807,12808],"tags":[],"_links":{"self":[{"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts\/10"}],"collection":[{"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/users\/1199"}],"replies":[{"embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/comments?post=10"}],"version-history":[{"count":21,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts\/10\/revisions"}],"predecessor-version":[{"id":12,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts\/10\/revisions\/12"}],"wp:attachment":[{"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/media?parent=10"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/categories?post=10"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/tags?post=10"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}