[R-sig-Geo] Defining a grid for interpolations ?

Edzer Pebesma edzer.pebesma at uni-muenster.de
Fri May 23 15:02:58 CEST 2008

Good question; the meuse example does not explain how to get meuse.grid.

The asciigrid file you have you should read with readGDAL in package 
rgdal, if that goes wrong I'd expect the problem in your data file.

Otherwise, if you have a shapefile with a polygon covering the 
catchment, called catch.pol, you'd want to use

sel = overlay(catch.pol, catch.grid)

The grid cells inside the polygon(s) are now found by

catch.sel = catch.grid[!is.na(sel)]

Mauricio Zambrano wrote:
> Dear List,
> My question is about how to define the grid to be used for the
> interpolations, using R 2.7.0 and gstat 0.9-47
> I'm working in a catchment of ~3000 km2, with daily rainfall data of
> several stations (more than 40), and I would like to interpolate daily
> values within all the catchment using ordinary Kriging.
> For defining the grid in which the interpolation will be carried out,
> at the beginning, I tried with
> #setting a grid each 100m vertical and horizontal
> dx <- seq(674400,730700,by=100)
> dy <- seq(4615100,4744400,by=100)
> catch.grid <- expand.grid(dx,dy)
> #setting the names of the columns of the grid
> names(catch.grid) <- c("x","y")
> #setting the coordinates of the grid
> coordinates(catch.grid) <- ~x+y
> #interpolating with the inverse distance
> pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_catch_nNA, catch.grid)
> #plotting the interpolated values
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]")
> but looking at the figure I realized that the interpolations were
> carried out considering all the cells within the squared grid, and not
> only within the catchment.
> After, I tried reading a raster file (each 100m) and using it as the grid,
> DEMM100m <- read.asciigrid("Catch_DEM_c100m")
> but the results that I got seems to be wrong, because I got high
> values where low values were expected and viceversa. (I really
> appreciate any help clarifying me what I'm doing wrong )
> Is there any way to define a grid, starting from a shapefile of the
> catchment boundaries ?. For example, I would like to define something
> similar to the "meuse.grid" dataset.
> Thanks in advance and best regards

More information about the R-sig-Geo mailing list