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

[Mauricio Zambrano]

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

-- 
Mauricio

Linux user #454569 -- Ubuntu user #17469