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

[Paul Hiemstra]

Hi Mauricio,

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

library(sp)
library(maptools)
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
summary(grd)
# Visualize the grid
spplot(source_poly, sp.layout = list("sp.points", grd))

"grd" can now be used for interpolations.

hth and cheers,
Paul

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
http://intamap.geo.uu.nl/~paul