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

[Dylan Beaudette]

On Wednesday 30 January 2008 01:13:55 am 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. 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.sg
>rd;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.sg
>rd;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.
>
>

This is all very interesting, but doesn't the GRASS-R combination already do 
these things- and quite well ? As far as I can tell GRASS can handle the 
massive grid operations, and R+gstat can do the statistical modeling, etc.

But maybe I should check out SAGA again- it wouldn't compile last time...

Thanks for the post,


-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341