[R-sig-Geo] R, raster & SAGA

[Rainer Hurling]

On 08.12.2010 19:55 (UTC+1), Andy Wilson wrote:
> Hi all,
> I'm new to R (and SAGA) and hope you can help resolve a problem which is
> perplexing me. As part of my project i would like to use the
> ta_morphometry module in the RSAGA library to calculate curvatures for a
> DEM. I've cut back my code to the basics below: read the DEM, write this
> into a SAGA format (I think I need to do this to make use of the SAGA
> functions), call the RSAGA function which outputs the data in SAGA
> format and finally read that data file (curve.sdat) into a raster
> object. It seems to work (no errors thrown), but when I try to do
> anything with the curve raster (including just plot it) I get the error
> message "Error in endrow - nrow : non-numeric argument to binary operator".
>
> This all works fine within the SAGA GUI, including plotting the curve
> data, but I still cant read the output data file into R using raster.
> Instead I have to export from SAGA as another format (eg Surfer) in
> order to be able to read the data into a raster object. However I'd like
> to be able to script the process without resorting to using the GUI to
> change the file format.
>
> Many thanks for your advice,
> Andy Wilson
>
> library (rgdal)
> library (RSAGA)
> library(raster)
>
> r_elev <- raster("C:/Dissertation/DEMs/DEM.IMG")
>
> #set environment for SAGA
> myenv <- rsaga.env(workspace="C:/Dissertation/SAGA")
> #output the DEM raster into a SAGA readable format
> writeRaster(r_elev, filename="C:/Dissertation/SAGA/elev.sgrd",
> format="SAGA", overwrite=TRUE)
> # call the module function using geoprocessor
> rsaga.geoprocessor("ta_morphometry", 0,
> list(ELEVATION="elev.sgrd",SLOPE="slope",ASPECT="aspect",
> CURV="curve",METHOD=5),env=myenv)
> #create a raster from the curve data file
> r_curve <- raster("C:/Dissertation/SAGA/curve.sgrd")

As a workaround until raster 1.7-11 is available the following step 
should work (note that you have to use '.sdat' instead of '.sgrd' with 
rgdal):

r_curve <- raster(readGDAL("C:/Dissertation/SAGA/curve.sdat"))

Hope this helps,
Rainer Hurling


> #plot the curve data
> plot(r_curve)