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

[Roger Bivand]

On Mon, 4 Feb 2008, Alí Santacruz wrote:

>
>
> Hi,
>
> I am trying to use the Grid Calculator in SAGA, through RSAGA but I have 
> problems with the parameter FORMUL,
>
> I can´t get the function "a/b" works,
>
> the functions "a+b", or "a-b", or "a*b" works, but "a/b" does not,
>
> Any help is appreciated,
>
> see the code below for the errors reported,
>
> thanks,
>
> Ali Santacruz
> M.Sc. Geomatics.
> National University of Colombia,
> Bogota, Colombia
> E-mail: amsantac at unal.edu.co
>
>
> ## this formula does not work:
>> rsaga.geoprocessor("grid_calculus", 1, 
>> list(INPUT="dem.sgrd;brillo.sgrd", RESULT="out2", 
>> FORMUL="a/b"))
>> C:\Users\geomatica\Desktop
>> C:\Users\geomatica\Documents\R\win-library\2.6\RSAGA\saga_vc\saga_cmd.exe
>> grid_calculus 1 -silent -INPUT dem.sgrd;brillo.sgrd -RESULT out2 
>> -FORMUL a\b
           ^^^^

A very smokey smoking gun, isn't it? The interface is converting all 
forward slash to backslash (Windows *does* respect forward slash in file 
path names, so this isn't strictly needed), so you'll have to try to 
escape the value until you see a/b coming through. I'm in doubt whether 
you'll have any luck, though:

...
     if (beep.off & .Platform$OS.type=="windows") {
         command = gsub("/","\\",command,fixed=TRUE)
                        ^^^^^^^^
         batch = c("net stop beep",command)
         batchfilename = paste(tempfile(),".bat",sep="")
         batchfile = file(batchfilename,"wt")
         writeLines(batch,con=batchfile)
         close(batchfile)
         command = batchfilename
     }
...

in R/RSAGA-core.R, which converts everything whether a file path or not.

Roger

PS. Is the RSAGA maintainer subscribed to this list?


> SAGA CMD 2.0
> Copyright (C) 2005 by Olaf Conradunder GNU General Public License (GPL)library path:   C:/Users/geomatica/Documents/R/win-library/2.6/RSAGA/saga_vc/moduleslibrary name:   grid_calculusmodule name:    Grid Calculator
> Load grid: dem.sgrd...
> Load grid: brillo.sgrd...
> Parameters:[Grid system] Grid system: 0.008462; 1506x 1585y; -79.300459x -4.492490y[Grid list] Grids: 2 objects (dem.sgrd, brillo.sgrd))[Grid] Result: new[Text] Formula: a\bError at character #0 of the function: a\b
>
> function not found
> error: executing module [Grid Calculator]
>
> # this formula does not work either:
>> rsaga.geoprocessor("grid_calculus", 1, list(INPUT="dem.sgrd;brillo.sgrd", RESULT="out2", FORMUL="a\b"))
> C:\Users\geomatica\Desktop>C:\Users\geomatica\Documents\R\win-library\2.6\RSAGA\saga_vc\saga_cmd.exe grid_calculus 1 -silent -INPUT dem.sgrd;brillo.sgrd -RESULT out2 -FORMUL
> SAGA CMD 2.0
> Copyright (C) 2005 by Olaf Conradunder GNU General Public License (GPL)library path:   C:/Users/geomatica/Documents/R/win-library/2.6/RSAGA/saga_vc/moduleslibrary name:   grid_calculusmodule name:    Grid Calculator
> Load grid: dem.sgrd...
> Load grid: brillo.sgrd...
> Parameters:[Grid system] Grid system: 0.008462; 1506x 1585y; -79.300459x -4.492490y[Grid list] Grids: 2 objects (dem.sgrd, brillo.sgrd))[Grid] Result: new[Text] Formula: Error at character #0 of the function:
>
> function not found
> error: executing module [Grid Calculator]
>
>
> # this formula works:
>
>> rsaga.geoprocessor("grid_calculus", 1, list(INPUT="dem.sgrd;brillo.sgrd", RESULT="out2", FORMUL="a+b"))C:\Users\geomatica\Desktop>C:\Users\geomatica\Documents\R\win-library\2.6\RSAGA\saga_vc\saga_cmd.exe grid_calculus 1 -silent -INPUT dem.sgrd;brillo.sgrd -RESULT out2 -FORMUL a+b
> SAGA CMD 2.0
> Copyright (C) 2005 by Olaf Conradunder GNU General Public License (GPL)library path:   C:/Users/geomatica/Documents/R/win-library/2.6/RSAGA/saga_vc/moduleslibrary name:   grid_calculusmodule name:    Grid Calculator
> Load grid: dem.sgrd...
> Load grid: brillo.sgrd...
> Parameters:[Grid system] Grid system: 0.008462; 1506x 1585y; -79.300459x -4.492490y[Grid list] Grids: 2 objects (dem.sgrd, brillo.sgrd))[Grid] Result: new[Text] Formula: a+bSave grid: out2...
>
>> rsaga.geoprocessor("grid_calculus", 1, list(INPUT="dem.sgrd;brillo.sgrd", RESULT="out2"))C:\Users\geomatica\Desktop>C:\Users\geomatica\Documents\R\win-library\2.6\RSAGA\saga_vc\saga_cmd.exe grid_calculus 1 -silent -INPUT dem.sgrd;brillo.sgrd -RESULT out2
> SAGA CMD 2.0
> Copyright (C) 2005 by Olaf Conradunder GNU General Public License (GPL)library path:   C:/Users/geomatica/Documents/R/win-library/2.6/RSAGA/saga_vc/moduleslibrary name:   grid_calculusmodule name:    Grid Calculator
> Load grid: dem.sgrd...
> Load grid: brillo.sgrd...
> Parameters:[Grid system] Grid system: 0.008462; 1506x 1585y; -79.300459x -4.492490y[Grid list] Grids: 2 objects (dem.sgrd, brillo.sgrd))[Grid] Result: new[Text] Formula: (a - b) / (a + b)Save grid: out2...
>
>
>> sessionInfo()R version 2.6.1 (2007-11-26) i386-pc-mingw32
> locale:LC_COLLATE=Spanish_Colombia.1252;LC_CTYPE=Spanish_Colombia.1252;LC_MONETARY=Spanish_Colombia.1252;LC_NUMERIC=C;LC_TIME=Spanish_Colombia.1252
> attached base packages:[1] stats     graphics  grDevices utils     datasets  methods   base
> other attached packages:[1] RSAGA_0.9-1
> loaded via a namespace (and not attached):[1] rcompgen_0.1-17 tools_2.6.1
>
>> From: hengl at science.uva.nl> To: r-sig-geo at stat.math.ethz.ch> Date: Wed, 30 Jan 2008 10:13:55 +0100> Subject: [R-sig-Geo] 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> > _______________________________________________> R-sig-Geo mailing list> R-sig-Geo at stat.math.ethz.ch> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> _________________________________________________________________
>
>
> e=wlmailtagline
> 	[[alternative HTML version deleted]]
>
>

-- 
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