Defining a grid for interpolations ?
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