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

[Mauricio Zambrano]

Dear Edzer,

Just in case it be useful for somebody else, I solved the problem with
the installation of "rgdal" in R, installing the "libgdal1-dev"
library on my linux system.

Thanks.

Mauricio

2008/5/23 Mauricio Zambrano <hzambran.newsgroups at gmail.com>:
> Sorry for not include the list in my last e-mail, I forgot it. :(.
>
>
>
> 2008/5/23 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
>> 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.
>>>
> Thanks, for clarifying me this.
>
>>> 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.
>
> I tried, but I got an error:
>
>> remove.packages("rgdal")
> Warning in remove.packages("rgdal") :
>  argument 'lib' is missing: using /usr/local/lib/R/site-library
>> install.packages("rgdal",dependencies=TRUE)
> Warning in install.packages("rgdal", dependencies = TRUE) :
>  argument 'lib' is missing: using '/usr/local/lib/R/site-library'
> --- Please select a CRAN mirror for use in this session ---
> Loading Tcl/Tk interface ... done
> probando la URL 'http://cran.miroir-francais.fr/src/contrib/rgdal_0.5-24.tar.gz'
> Content type 'application/x-gzip' length 4426322 bytes (4.2 Mb)
> URL abierta
> ==================================================
> downloaded 4.2 Mb
>
> * Installing *source* package 'rgdal' ...
> gdal-config: gdal-config
> ./configure: line 1311: gdal-config: command not found
>
> The gdal-config script distributed with GDAL could not be found.
> If you have not installed the GDAL libraries, you can
> download the source from  http://www.gdal.org/
> If you have installed the GDAL libraries, then make sure that
> gdal-config is in your path. Try typing gdal-config at a
> shell prompt and see if it runs. If not, use:
>  --configure-args='--with-gdal-config=/usr/local/bin/gdal-config' echo
> with appropriate values for your installation.
>
> ERROR: configuration failed for package 'rgdal'
> ** Removing '/usr/local/lib/R/site-library/rgdal'
>
> The downloaded packages are in
>        /tmp/Rtmphf9aT7/downloaded_packages
> Warning message:
> In install.packages("rgdal", dependencies = TRUE) :
>  installation of package 'rgdal' had non-zero exit status
>
> Do you have any idea about how to solve it ?
>
>>>
>>> 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.
>>>
>
> Sorry, you are right, I got the error when I tried:
> #interpolating with the inverse distance
> pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_catch_nNA, catchment.grid)
>
>>> 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
>> ?
>
> Thanks a lot, I didn't know that function, and looking at the help of
> "coordinates" I didn't find it.
> This allow me to run the interpolation in the right projected system:
>
> pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_catch_nNA, catchment.grid)
>
> Thanks a lot.
>
> Mauricio
>
>> --
>> 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
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>
>>
>
>
>
> --
> hzambran
>
> Linux user #454569 -- Ubuntu user #17469
>



-- 
hzambran

Linux user #454569 -- Ubuntu user #17469