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

Paul Hiemstra p.hiemstra at geo.uu.nl
Mon May 26 10:37:10 CEST 2008

Hi Mauricio,

To get a grid based on a shapefile you can use the command "spsample".

source_poly = readShapePoly("/path/to/poly")
# cellsize is in map units (e.g. km), also see "?spsample"
grd = spsample(source_poly, cellsize = c(10e3,10e3), type = "regular")
gridded(grd) = TRUE # Make it a grid
# Visualize the grid
spplot(source_poly, sp.layout = list("sp.points", grd))

"grd" can now be used for interpolations.

hth and cheers,

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

Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: 	+31302535773
Fax:	+31302531145

More information about the R-sig-Geo mailing list