[R-sig-Geo] FYI: Merging GIS and statistics --- RSAGA

Wed Jan 30 10:59:23 CET 2008

On Wed, 30 Jan 2008, Tomislav Hengl wrote:

>
>
> Dear all,
>
> I just started running analysis with the RSAGA package 
> (http://cran.r-project.org/src/contrib/Descriptions/RSAGA.html), i.e. 
> the R scripting link to SAGA GIS (by Olaf Conrad and colleagues, over 
> 120 modules), that was suggested to me by Paulo van Breugel and I think 
> that this could really be the missing link between statistics and GIS.

No, there are always more missing links, never a silver bullet. SAGA is 
hard to install from source, and it appears that there is not much 
activity on the SF CVS, and few install instructions. It seems to have 
multiple interface dependencies which Linux users will find hard to 
satisfy. It may work as a Windows binary install. The RSAGA maintainer 
does not seem to read this list - my feeling is that SAGA lives in a world 
of its own.

Re. memory problems on Windows, please see the R for Windows FAQ, and 
remember firstly that there is no R for 64-bit windows yet, and that 
Windows memory management (despite much work being done within R to try to 
alleviate things), is inferior to OSX and Linux, because all requests have 
to be met in one chunk. If you want 300MB, Windows will turn you down if 
there isn't a single chunk of that size free, while the alternatives 
collect the chunks they have and hand off addresses making it look like 
continuous memory to the application (if I understand correctly). This 
means that some tasks fail under Windows but succeed under other OS on the 
same hardware.

If SAGA can make a source library for reading and writing its raster file 
formats available to GDAL, and maybe help write a driver, SAGA access 
through rgdal will happen automatically. Lots of other projects do this, 
for example, the PCRaster format is included as source in the GDAL source. 
Were SAGA to split out the raster I/O as a library and provide a copy to 
GDAL, your question would be answered. SAGA does use GDAL for interfacing 
other formats, and would be a "good citizen" if they reciprocated.

I don't see the SAGA format documented (there doesn't seem to be much 
documentation for programmers), but if you know what it is, you can use 
readBin() and writeBin() in R to construct a portable interface.

SAGA may be a good idea, but there are plenty of very viable alternatives.

Roger

> My experiences so far are very positive --- especially if you work with 
> large grids, because SAGA is quite fast for calculations. Here are some 
> examples from Geomorphometry / Digital Soil Mapping:
>
> 0. Getting started:
>
> **************************************************************************** 
> # Download the SAGA 2.0.1 binaries 
> (http://sourceforge.net/projects/saga-gis/) and unzip them to a local 
> directory e.g. "C:/Progra~1/saga_vc"; # Start R and install the RSAGA 
> package; # load the library and set the directory where the SAGA 
> binaries sit:
>
> library(RSAGA)
> rsaga.env(path="C:/Progra~1/saga_vc")
>
> # To get the exact names of parameters look for a name in the "/modules" directory and then use:
>
> rsaga.get.modules("geostatistics_kriging")
> rsaga.get.usage("geostatistics_kriging", 2)
>
> ****************************************************************************
>
> 1. Error propagation and geomorphometry (both can be run via R now):
>
> ****************************************************************************
>
> # Import the point measurements of heights to generate a DEM:
>
> elevations <- read.delim("elevations.txt") coordinates(elevations)=~X+Y
> spplot(elevations)
>
> # Import the grid definition:
>
> gridmaps = readGDAL("SMU1.asc")
> gridmaps$SMU1 = gridmaps$band1
>
> # Derive area in km^2:
>
> maparea =
> (gridmaps at bbox["x","max"]-gridmaps at bbox["x","min"])*(gridmaps at bbox["y","max"]-gridmaps at bbox["y","min
> "])/1e+06
>
> # Fit a variogram for elevations and produce 50 realizations of a DEM using Sequential Gaussian
> Simulations:
>
> elevations.or = variogram(Z~1, elevations)
> elevations.ovgm = fit.variogram(elevations.or, vgm(1, "Sph", 1000, 1))
> plot(elevations.or, elevations.ovgm, plot.nu=F, pch="+")
>
> DEM.sim = krige(Z~1, elevations, gridmaps, elevations.ovgm, nmax=40, nsim=50)
>
> # Visualize the simulated DEMs in R:
>
> for (i in 1:length(DEM.sim at data)) {
>      image(as.image.SpatialGridDataFrame(DEM.sim[i]), col=terrain.colors(16), asp=1) }
>
> # Write the simulated DEMs in ArcInfo ASCII format:
>
> for (i in 1:length(DEM.sim at data)) {
>      write.asciigrid(DEM.sim[i], c(paste("DEM",as.character(i),".asc",sep="")))
> }
>
> # Now, derive SLOPE maps in SAGA 50 times:
> # ESRI wrapper is used to get the maps directly in ArcInfo ASCII format;
>
> for (i in 1:length(DEM.sim at data)) {
>   rsaga.esri.wrapper(rsaga.slope, method="poly2zevenbergen",
> in.dem=c(paste("DEM",as.character(i),sep="")), out.slope=c(paste("SLOPE",as.character(i),sep="")),
> prec=3, condensed.res=FALSE, intern=FALSE, show.output.on.console=FALSE) }
>
> # Optional: generate a DEM using the Thin Plate Spline (local) interpolation in SAGA:
>
> writeOGR(elevations, "elevations.shp", "elevations", "ESRI Shapefile")
>
> rsaga.get.usage("grid_spline", 1) rsaga.geoprocessor(lib="grid_spline", module=1,
> param=list(GRID="DEMtps.sgrd", SHAPES="elevations.shp", FIELD=1, RADIUS=sqrt(maparea)*1000/3,
> SELECT=1, MAXPOINTS=30, TARGET=2, GRID_GRID="DEM1.sgrd")) rsaga.sgrd.to.esri(in.sgrds="DEMtps.sgrd",
> out.grids="DEMtps.asc", out.path="D:/GEOSTAT/maps/RSAGA", prec=1)
>
>
> ****************************************************************************
>
> 2. Spatial interpolation
> Especially suitable for large maps (R+gstat often fail due to memory limit problems):
>
> ****************************************************************************
> # Export the predictors to SAGA format:
>
> predict.list = gl(n=9, k=1,
> labels=c("DEM","SLOPE","PLANC","TWI","SINS","SMU1","SMU3","SMU4","SMU9"))
> rsaga.esri.to.sgrd(in.grids=levels(predict.list),
> out.sgrds=set.file.extension(levels(predict.list),".sgrd"), in.path="D:/GEOSTAT/maps/RSAGA")
>
> # predict values in SAGA using only regression model:
>
> rsaga.get.usage("geostatistics_grid", 4) rsaga.geoprocessor(lib="geostatistics_grid", module=4,
> param=list(GRIDS="DEM.sgrd;SLOPE.sgrd;PLANC.sgrd;TWI.sgrd;SINS.sgrd;SMU1.sgrd;SMU3.sgrd;SMU4.sgrd;SM
> U9.sgrd", SHAPES="baranja.shp", ATTRIBUTE=0, TABLE="regout.dbf", RESIDUAL="solum_res.shp",
> REGRESSION="SOLUM_reg.sgrd", INTERPOL=0))
>
> # Ordinary kriging:
>
> rsaga.get.usage("geostatistics_kriging", 1) rsaga.geoprocessor(lib="geostatistics_kriging",
> module=1, param=list(GRID="SOLUM_ok.sgrd", VARIANCE="SOLUM_okvar.sgrd", SHAPES="baranja.shp",
> FIELD=0, MODEL=1, NUGGET=0, SILL=200, RANGE=500, TARGET=2, GRID_GRID="SLOPE.sgrd"))
>
> # Regression-kriging:
>
> rsaga.get.usage("geostatistics_kriging", 3) rsaga.geoprocessor(lib="geostatistics_kriging",
> module=3,
> param=list(GRIDS="DEM.sgrd;SLOPE.sgrd;PLANC.sgrd;TWI.sgrd;SINS.sgrd;SMU1.sgrd;SMU3.sgrd;SMU4.sgrd;SM
> U9.sgrd", GRID="SOLUM_rk.sgrd", SHAPES="baranja.shp", FIELD=0, MODEL=1, NUGGET=0, SILL=200,
> RANGE=500, INTERPOL=0))
> # Does not work yet. Possibly a bug in the saga_cmd.exe?
>
> ****************************************************************************
>
> The complete script and datasets are available at:
>
> http://spatial-analyst.net/GRK/examplesRSAGA.zip   (400 KB)
>
> So the only real problem is the import/export from R to SAGA, which I guess could be solved very
> easily if the next version of rgdal would support SAGA format.
>
>
> Tom Hengl
> http://spatial-analyst.net
>
> _______________________________________________
> 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