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 = > ([EMAIL PROTECTED]"x","max"[EMAIL PROTECTED]"x","min"])*([EMAIL > PROTECTED]"y","max >"[EMAIL PROTECTED]"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([EMAIL PROTECTED])) { > 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([EMAIL PROTECTED])) { > 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([EMAIL PROTECTED])) { > 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 _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
