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

[Edzer Pebesma]

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
>>
>>   
>
>