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

[Edzer Pebesma]

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