Defining a grid for interpolations ?

Mauricio, please keep r-sig-geo in the loop.

Mauricio Zambrano wrote:

 Dear Ezder,

 Thanks for your advice about the readGdal and the "overlay" function,
 but why is it better to use "readGdal" than "readasciigrid" ?.

because the latter one is deprecated.

 In any case, when I tried to load "rgdal" I got the following error
 message:

 "Error in fun(...) :
 	GDAL Error 1: libgrass_I.so: no se puede abr?r el archivo de objeto
 compartido: No existe el fichero ? directorio

 Error : .onLoad failed in 'loadNamespace' for 'rgdal'
 Error: package/namespace load failed for 'rgdal'

 library(rgdal)

 Error in fun(...) :
 	GDAL Error 1: libgrass_I.so: no se puede abr?r el archivo de objeto
 compartido: No existe el fichero ? directorio

 Error : .onLoad failed in 'loadNamespace' for 'rgdal'
 Error: package/namespace load failed for 'rgdal'"

 I don't know what is the problem, but yesterday I installed the new
 "qgis" and I don't know if this can be the cause of the error.

It wouldn't surprise me. You may want to try to reinstall rgdal.

 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

Please tell us where exactly you got which error. I suspect you got it a bit 
later.

 points with existing coordinates ?. At the beginning I did:

 #reading data from the file
 meteo<- read.csv2("Meteorological_by_basin.csv")

 #setting the coordinates
 coordinates(meteo) <- ~EAST_ED50 + NORTH_ED50

Did you try
proj4string(meteo) = p4s
?
--
Edzer (not Ezder)

 Thanks in advance,

 Mauricio

 2008/5/23 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:

 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