How to interpolate point value using kriging method with R and GRASS 6

From GRASS-Wiki
Jump to navigation Jump to search

ORDINARY KRIGING IN R WITH GRASS6 DATA

-- New development mid-2009:

--

(WARNING!! Most of the code quoted here is very out of date, and simply does not work for current R/sp/gstat/spgrass6. Untried suggestions have been edited in, but without a test location, there is no guarantee that they will work! Roger Bivand, 5 April 2007)

Of all the methods we tried this is the most easy and (I suppose) exact too:

You have to have in your library the packages "gstat" and "spgrass6", you can download this last one directly from R using the command "install.packages". In GRASS we have a vector file named "giaciture_cat_clean3" and we want to do a prediction on this data... these are the commmands:

enter R from the GRASS prompt, and type: 

library(spgrass6) 

#get vector points as SpatialPointsDataFrame 
#giaciture <- getSites6sp("giaciture_cat_clean3") RSB 070405

giaciture <- readVECT6("giaciture_cat_clean3") # RSB 070405

class(giaciture) #shows the class of "giaciture" (SpatialPointsDataFrame)
 
# G <- gmeta6() #get region from GRASS to R  RSB 070405

now if you want you can continue to work in R from GRASS or not... 

#create a grid from the region settings of GRASS, it is very important
# to have square cells, so you can set the region settings of GRASS or
# you can give directly square dimensions using the values:  
# e.g."cells.dim=c(50,50)"
#grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2)
#                    ,G$south+(G$nsres/2))
#                    ,cellsize=c(G$ewres, G$nsres)
#                    ,cells.dim=c(G$cols, G$rows)) RSB 070405

grd <- gmeta2grd() # RSB 070405

#create a SpatialGridDataFrame
mask_SG <- SpatialGridDataFrame(grd,
#                                ,data=list(k=rep(1, G$cols*G$rows)) RSB 070405
     data=data.frame(k=rep(1, prod(slot(grd, "cells.dim")))), # RSB 070405
     proj4string=CRS(proj4string(giaciture))) # RSB 070405
# proj4string(giaciture) and proj4string(mask_SG) must agree
class(mask_SG)

library(gstat)

cvgm <- variogram(IMMERSIONE~1, data=giaciture, width=400, cutoff=4000) # RSB 070405
#create variogram, and "IMMERSIONE" 
#here is the our variable, the variable on wich we have to do the prediction,
# ~ 1 select the type of kriging, this is the ordinary one

efitted <- fit.variogram(cvgm, vgm(psill=5000, model="Exp", range=1500, nugget=8000))
# choose the model to fit variogram (here is exponential) and give the
# estimated parameters of the variogram (partial sill, range and nugget)

OK_pred <- krige(IMMERSIONE~1, data=giaciture, newdata=mask_SG, model=efitted) # RSB 070405
# make the kriging prediction

names(OK_pred) #show the name of variable kriged

writeRAST6(OK_pred, "OK_pred", zcol="var1.pred") # RSB 070405
#write a raster file and save it in GRASS, now you can open it from there.

that's all!

special thanks to Roger Bivand, ever ready to lend a hand!

More Help

Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=How_to_interpolate_point_value_using_kriging_method_with_R_and_GRASS_6&oldid=19626"