Skip to content

FYI: Merging GIS and statistics --- RSAGA

7 messages · Tomislav Hengl, Roger Bivand, Johan Van de Wauw +3 more

Original
Tomislav Hengl
# 
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
Roger Bivand
#
On Wed, 30 Jan 2008, Tomislav Hengl wrote:

            

      
      
      
    

        
No, there are always more missing links, never a silver bullet. SAGA is 
hard to install from source, and it appears that there is not much 
activity on the SF CVS, and few install instructions. It seems to have 
multiple interface dependencies which Linux users will find hard to 
satisfy. It may work as a Windows binary install. The RSAGA maintainer 
does not seem to read this list - my feeling is that SAGA lives in a world 
of its own.

Re. memory problems on Windows, please see the R for Windows FAQ, and 
remember firstly that there is no R for 64-bit windows yet, and that 
Windows memory management (despite much work being done within R to try to 
alleviate things), is inferior to OSX and Linux, because all requests have 
to be met in one chunk. If you want 300MB, Windows will turn you down if 
there isn't a single chunk of that size free, while the alternatives 
collect the chunks they have and hand off addresses making it look like 
continuous memory to the application (if I understand correctly). This 
means that some tasks fail under Windows but succeed under other OS on the 
same hardware.

If SAGA can make a source library for reading and writing its raster file 
formats available to GDAL, and maybe help write a driver, SAGA access 
through rgdal will happen automatically. Lots of other projects do this, 
for example, the PCRaster format is included as source in the GDAL source. 
Were SAGA to split out the raster I/O as a library and provide a copy to 
GDAL, your question would be answered. SAGA does use GDAL for interfacing 
other formats, and would be a "good citizen" if they reciprocated.

I don't see the SAGA format documented (there doesn't seem to be much 
documentation for programmers), but if you know what it is, you can use 
readBin() and writeBin() in R to construct a portable interface.

SAGA may be a good idea, but there are plenty of very viable alternatives.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 30 Jan 2008, Johan Van de Wauw wrote:

        

            

      
      
      
    

        
Please do - or if it is a BIL or similar, just saying which existing GDAL 
driver also works for SAGA would be great. If it is very similar, it might 
be enough to tweak a copy of an existing GDAL format.

            

      
      
      
    

        
That's very helpful. Could you please say whether the coordinate reference 
systems are recognised correctly on both sides, and if you need any extra 
arguments to write (presumably) using functions from the rgdal package?

Roger

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dylan Beaudette
    

    

      
      #
    

  


  

      
On Wednesday 30 January 2008 01:13:55 am Tomislav Hengl wrote:

            

      
      
      
    

        
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,

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Barry Rowlingson
    

    

      
      #
    

  


  

      
Agustin Lobo wrote:

            

      
      
      
    

        
Now this could be done with a GDAL/OGR driver for R-spatial objects 
(stored in .RData files). Then you could use any GIS with GDAL/OGR 
capabilities to make pretty maps with R objects.

  I'm not really au fait with writing and linking GDAL/OGR drivers, but 
the code would have to link with libR and require the sp package... Fun. 
Anyone up for it?

Barry