How to interpolate point value using kriging method with R and GRASS 6: Difference between revisions
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: | |||
> | >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 | #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)),\ | |||
>mask_SG <- SpatialGridDataFrame(grd,data=list(k=rep(1,G$cols*G$rows)),proj4string=CRS(G$proj4)) | cellsize=c(G$ewres, G$nsres), cells.dim=c(G$cols, G$rows)) | ||
>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) | ||
here is the our variable, the variable on | #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 | ||
> | #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) | #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: | |||
>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!