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

[Mauricio Zambrano]

Dear Paul and Edzer,

Thanks a lot for your answers.

Before reading the solution proposed by Paul, I had tried with:

# reading the boundary of the catchment
library(maptools)
catchment <- readShapePoly("only_catchment.shp")
catchment.grid <- spsample(catchment, cellsize=100, offset = c(0.5,
0.5)) # reading the boundary of the catchment

 and then I did an IDW interpolation with:

pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo, catchment.grid)

and it works.

However, after reading your solution, I realised that I didn't use the command:

gridded(catchment.grid)

but I got an interpolation anyway.

Where is performed the interpolation when I don't use a gridded grid (
gridded(catchment.grid) ) ?

Best regards

Mauricio

2008/5/26 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
> 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
>>>
>>>
>>
>>
>
>



-- 
hzambran

Linux user #454569 -- Ubuntu user #17469