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

Wed Jan 30 17:55:58 CET 2008

Good job! This is more or less like the R-GRASS link, but SAGA
does different things, so it's good having this additional tool available.
Anyway, I do not think that these links are the real solution, which
should be being able to display R spatial objects on a geographical
display in which the information could be interactively consulted
and overlayed with other geographic information. So the work is
perhaps more on the GIS side, which should be able to represent
R spatial objects, or provide R with a real geographical display.

Going back and forth with geotif rasters and/or shp vectors soon becomes 
inconvenient.

Agus

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

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster