{"id":24,"date":"2010-10-15T16:00:26","date_gmt":"2010-10-15T20:00:26","guid":{"rendered":"http:\/\/web.colby.edu\/mgimond\/?p=24"},"modified":"2011-10-27T12:26:05","modified_gmt":"2011-10-27T16:26:05","slug":"24","status":"publish","type":"post","link":"https:\/\/web.colby.edu\/mgimond\/2010\/10\/15\/24\/","title":{"rendered":"R: Kriging"},"content":{"rendered":"
\n

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

\r\nlibrary(maptools)\r\n\r\nS = readShapePoints("points_krig.shp")\r\n\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 add the attribute values to P\r\n#\r\n# thus creating a marked point pattern\r\nhist(marks(P)$z) # plots histogram of attribute z\r\n<\/pre>\n
\"\"<\/span><\/p>\n
\r\n# Convert data to dataframe\r\nP.frame = as.data.frame(P)\r\n# Only keep necessary columns\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# Check plot (note the bullseye effect of the standard interpolation)\r\nplot(density(P,10),axes=T)\r\nplot(P, add=T)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Use geoR package\r\n# Create a trend surface\r\nlibrary(geoR)\r\nP.cloud=variog(coords = cbind(X,Y), data=Z, option=c("cloud"))\r\nplot(P.cloud)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\nP.bin = variog(coords = cbind(X,Y), data=Z, uvec=seq(0,100,l=10))\r\nplot(P.bin)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n#Bins can be displayed with a box-plot for each bin\r\nP.bin1 = variog(coords = cbind(X,Y), data=Z, uvec=seq(0,100,l=5), bin.cloud = T)\r\nplot(P.bin1, bin.cloud=T,main="Classical estimator")\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Theoretical and empirical variograms can be plotted and visually compared\r\nplot(P.bin)\r\nlines.variomodel(cov.model= "spherical", cov.pars = c(100,95), \r\n                 nugget = 0, max.dist=120,lwd=1, col=2)\r\nlines.variomodel(cov.model= "exponential", cov.pars = c(100,95), \r\n                 nugget = 0, max.dist=120,lwd=1, col=3)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# computing a directional variogram\r\nplot(0,0,xlab="Distance",ylab="Semivariogram",xlim=c(0,120),ylim=c(0,1.6),cex.lab=1.2,type="n")\r\ntitle("Directional variogram",cex.main=1.2)\r\nlines(P.bin,pch=1,type="p")\r\nvario.0 = variog(coords=cbind(X,Y), data=Z, lambda=0, max.dist=120, dir=0, tol=pi\/4)\r\nvario.90 = variog(coords=cbind(X,Y), data=Z, lambda=0, max.dist=120, dir=pi\/4, tol=pi\/4)\r\nvario.45 = variog(coords=cbind(X,Y), data=Z, lambda=0, max.dist=120, dir=pi\/8, tol=pi\/4)\r\nvario.120 = variog(coords=cbind(X,Y), data=Z, lambda=0, max.dist=120, dir=3*pi\/8, tol=pi\/4)\r\nlines(vario.0,lty=1,pch=2)\r\nlines(vario.90,lty=2,pch=3)\r\nlines(vario.45,lty=3,pch=4)\r\nlines(vario.120,lty=4,pch=5)\r\nlegend("topleft", legend=c("omnidirectional","N-S","E-W","NW-SE","NE-SW"), pch=c(1,2,3,4,5))\r\n#\r\n# NOTE: pair distances are stored in P.cloud$u (P.cloud is a variogram class)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Using the maximum likelihood technique\r\n# Nugget fit to 0, sill=100, range=95\r\nml = likfit(coord=cbind(X,Y),data=Z, ini = c(100,95), fix.nugget = T) \r\nplot(P.bin1)\r\nlines(ml,col=2)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Predict a single value\r\nplot(cbind(X,Y), xlim=c(0,100), ylim=c(0,100), xlab="Coord X", ylab="Coord Y")\r\nloci =matrix(c(50,50),ncol=2)\r\ntext(loci, as.character(1), col="red")\r\nkc4 = krige.conv(coord=cbind(X,Y),data=Z, locations = loci, krige = krige.control(obj.m = ml))\r\n<\/pre>\n
\r\n> kc4$predict\r\n     data\r\n      10.74597\r\n      data 10.74597\r\n# kc4$predict is the predicted value (e.g. 10.7) and \r\n# kc4$krige.var is the variance (e.g. 9.4)\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Predict a grid\r\npred.grid = expand.grid(seq(0,100,l=51), seq(0,100,l=51)) # Define grid\r\n# kriging calculations\r\nkc = krige.conv(coord=cbind(X,Y),data=Z, loc=pred.grid, krige = krige.control(obj.m = ml))\r\nimage(kc, loc = pred.grid, col=gray(seq(1,0.1,l=30)))\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\n# Using the "fields" library\r\nlibrary(fields)\r\nfit = Krig(cbind(X,Y),Z,theta=95)\r\nsummary( fit) # summary of fit\r\nset.panel( 2,2)\r\nplot(fit) # four diagnostic plots of fit\r\n<\/pre>\n
\"\"<\/a><\/span><\/p>\n
\r\nset.panel()\r\nsurface( fit, type="C") # Display the surface\r\n<\/pre>\n
\"\"<\/a><\/p>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
!! Data for this exercise can be downloaded from 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\/24"}],"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=24"}],"version-history":[{"count":27,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts\/24\/revisions"}],"predecessor-version":[{"id":28,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/posts\/24\/revisions\/28"}],"wp:attachment":[{"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/media?parent=24"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/categories?post=24"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/web.colby.edu\/mgimond\/wp-json\/wp\/v2\/tags?post=24"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}