How to interpolate point value using kriging method with R and GRASS 6: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(layout fixed)
(+cat)
Line 52: Line 52:


special thanks to Roger Bivand, even ready to lend a hand!
special thanks to Roger Bivand, even ready to lend a hand!
[[Category:Documentation]]

Revision as of 19:56, 16 June 2006

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 in R from GRASS and digit:

>library(spgrass6)

#get vector points as SpatialPointsDataFrame:
>giaciture <- getSites6sp("giaciture_cat_clean3") 

#shows the class of "giaciture" (SpatialPointsDataFrame):
>class(giaciture)

#get region from GRASS to R:
>G <- gmeta6()

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)) 

#create a SpatialGridDataFrame:
>mask_SG <- SpatialGridDataFrame(grd,data=list(k=rep(1,G$cols*G$rows)),\
 proj4string=CRS(G$proj4)) >class(mask_SG)
>library(gstat)

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

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

# make the kriging prediction:
>OK_pred <- krige(IMMERSIONE~ 1,locations=giaciture,newdata=mask_SG,model=efitted)
>names(OK_pred) #show the name of variable kriged

#write a raster file and save it in GRASS, now you can open it from there:
>writeRast6sp(OK_pred,"OK_pred",zcol="var1.pred",NODATA=-9999)

that's all!

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