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

[Mauricio Zambrano]

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