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

Fri May 23 18:10:56 CEST 2008

On Fri, 23 May 2008, Edzer Pebesma wrote:

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

Yes, it's almost certainly QGIS chewing things up. It looks as though QGIS 
expects the GRASS OGR plugin to be present, or some GRASS dependencies to 
be present - if they are, you need to add the offending GRASS library 
directory to ld.so.conf or similar as root, and refresh ldconfig. On 
Linux, if the above does not work, best remove QGIS and reinstall GDAL 
(from source, not binary, the binaries often have unfulfilled 
dependencies, unfortunately). Then remove and re-install rgdal.

Roger

>>  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
>> > > 
>> > > 
>> > > 
>> > 
>> 
>> 
>> 
>> 
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no