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

(12 intermediate revisions by 5 users not shown) | |||

Line 1: | Line 1: | ||

− | ==ORDINARY KRIGING IN R WITH GRASS6 DATA== | + | === 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 | + | enter R from the GRASS prompt, and type: |

+ | |||

<pre> | <pre> | ||

− | + | 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 | ||

</pre> | </pre> | ||

Line 17: | Line 35: | ||

<pre> | <pre> | ||

− | + | #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. | ||

</pre> | </pre> | ||

that's all! | that's all! | ||

− | special thanks to Roger Bivand, | + | 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 02:49, 9 August 2013

### 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 v.krige module

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

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

- A nice working example for recent versions of R/sp/gstat/spgrass6.
- Kriging