Defining a grid for interpolations ?

lbrary(maptools)
elevations <- readShapePoints("elevations.shp", proj4string=CRS(as.character(NA)))
library(akima)
elevations.interp <- interp(x=elevations at coords[,1],y=elevations at coords[,2],z=elevations$VALUE,

image(elevations.interp, col=bpy.colors(), asp=1)
library(spatstat)
dem.interp <- as.SpatialGridDataFrame.im(as.im(elevations.interp))
dem.interp at grid

Please note that this "regular" grid is still random, as the first point is
sampled randomly from the first cell. For a fixed, reproducable grid, add
the argument offset = c(0.5, 0.5)
--
Edzer

Paul Hiemstra wrote:

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