In the meantime, I solved the problem using something very similar to
the advice given by Oehler:
# reading the boundary of the catchment
library(maptools)
catchment <- readShapePoly("only_catchment.shp")
catchment.grid <- spsample(catchment, type="regular", cellsize=100)
and it works.
However, if i did:
library(maptools)
#projection string for the UTM ED50, Z30N
p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")
catchment <- readShapePoly("only_catchment.shp", proj4string=p4s)
catchment.grid <- spsample(catchment, type="regular", cellsize=100)
I got an error, because the grid is projected but the coordinates of
the stations are not. How can assign a projection to a given set of
Thanks in advance,
Mauricio
2008/5/23 Edzer Pebesma <[EMAIL PROTECTED]>:
> Good question; the meuse example does not explain how to get meuse.grid.
>
> The asciigrid file you have you should read with readGDAL in package
> rgdal,
> if that goes wrong I'd expect the problem in your data file.
>
> Otherwise, if you have a shapefile with a polygon covering the
> catchment,
> called catch.pol, you'd want to use
>
> sel = overlay(catch.pol, catch.grid)
>
> The grid cells inside the polygon(s) are now found by
>
> catch.sel = catch.grid[!is.na(sel)]
> --
> Edzer
>
> 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
> >
> >
> >
>
