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
No edit summary
 
(11 intermediate revisions by 5 users not shown)
Line 1: Line 1:
=== ORDINARY KRIGING IN R WITH GRASS6 DATA ===
-- New development mid-2009:
* GRASS's [[Kriging]] wiki page
* [[v.krige_GSoC_2009]] wxGUI Kriging project by Anne Ghisla
** Anne's {{addonCmd|v.krige}} module
* {{addonCmd|v.autokrige}} module by Mathieu Grelier
--
('''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:
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...
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:
these are the commmands:


enter in R from GRASS and digit:
enter R from the GRASS prompt, and type:
 
<pre>
<pre>
>library(spgrass6)  
library(spgrass6)  
>giaciture <- getSites6sp("giaciture_cat_clean3") #get vector points as SpatialPointsDataFrame
 
>class(giaciture) #shows the class of "giaciture" (SpatialPointsDataFrame)
#get vector points as SpatialPointsDataFrame
>G <- gmeta6() #get region from GRASS to R
#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
</pre>
</pre>


Line 16: Line 35:


<pre>
<pre>
>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 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)"
#create a grid from the region settings of GRASS, it is very important
>mask_SG <- SpatialGridDataFrame(grd,data=list(k=rep(1,G$cols*G$rows)),proj4string=CRS(G$proj4)) #create a SpatialGridDataFrame
# to have square cells, so you can set the region settings of GRASS or
>class(mask_SG)
# you can give directly square dimensions using the values: 
>library(gstat)
# e.g."cells.dim=c(50,50)"
>cvgm <- variogram(IMMERSIONE~1,locations=giaciture,width=400,cutoff=4000) #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
#grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2)
>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)
#                    ,G$south+(G$nsres/2))
>OK_pred <- krige(IMMERSIONE~ 1,locations=giaciture,newdata=mask_SG,model=efitted) # make the kriging prediction
#                    ,cellsize=c(G$ewres, G$nsres)
>names(OK_pred) #show the name of variable kriged
#                    ,cells.dim=c(G$cols, G$rows)) RSB 070405
>writeRast6sp(OK_pred,"OK_pred",zcol="var1.pred",NODATA=-9999) #write a raster file and save it in GRASS, now you can open it from there.
 
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.
</pre>
</pre>


that's all!  
that's all!  


special thanks to Roger Bivand, even ready to lend a hand!
special thanks to Roger Bivand, ever ready to lend a hand!
 
=== More Help ===
 
* A [http://casoilresource.lawr.ucdavis.edu/drupal/node/438 nice working example] for recent versions of R/sp/gstat/spgrass6.
* [[Kriging]]
 
[[Category:Documentation]]
[[Category:R]]
[[Category: Tutorial]]
[[Category: HowTo]]

Latest revision as of 09:49, 9 August 2013

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