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
(layout fixed)
Line 6: Line 6:
enter in R from GRASS and digit:
enter in R from GRASS and digit:
<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")  
 
#shows the class of "giaciture" (SpatialPointsDataFrame):
>class(giaciture)
 
#get region from GRASS to R:
>G <- gmeta6()
</pre>
</pre>


Line 15: Line 21:


<pre>
<pre>
>grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2),G$south+(G$nsres/2)),cellsize=c(G$ewres, G$nsres),
#create a grid from the region settings of GRASS, it is very important to have
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
square cells, so you can set the region settings of GRASS or you can give directly square dimensions
# directly square dimensions using the values: e.g."cells.dim=c(50,50)":
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)),\
>mask_SG <- SpatialGridDataFrame(grd,data=list(k=rep(1,G$cols*G$rows)),proj4string=CRS(G$proj4)) #create a
cellsize=c(G$ewres, G$nsres), cells.dim=c(G$cols, G$rows))  
SpatialGridDataFrame
 
>class(mask_SG)
#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)
>library(gstat)
>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  
#create variogram, and "IMMERSIONE" here is the our variable, the variable on  
is the ordinary one
#which we have to do the prediction, ~ 1 select the type of kriging, this  
>efitted <- fit.variogram(cvgm,vgm(psill=5000,model="Exp",range=1500,nugget=8000)) #choose the model to fit  
#is the ordinary one:
variogram (here is exponential) and give the estimated parameters of the variogram (partial sill, range and  
>cvgm <- variogram(IMMERSIONE~1,locations=giaciture,width=400,cutoff=4000)  
nugget)
 
>OK_pred <- krige(IMMERSIONE~ 1,locations=giaciture,newdata=mask_SG,model=efitted) # make the kriging prediction
#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
>names(OK_pred) #show the name of variable kriged
>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.
#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)
</pre>
</pre>



Revision as of 12:00, 23 May 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!