[R-sig-Geo] Defining a grid for interpolations ?
edzer.pebesma at uni-muenster.de
Mon May 26 13:54:33 CEST 2008
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)
Paul Hiemstra wrote:
> Hi Mauricio,
> To get a grid based on a shapefile you can use the command "spsample".
> 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
> # Visualize the grid
> spplot(source_poly, sp.layout = list("sp.points", grd))
> "grd" can now be used for interpolations.
> hth and cheers,
> 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
>> 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
More information about the R-sig-Geo