FYI: Merging GIS and statistics --- RSAGA

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.

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