[R-sig-Geo] Including SAGA grid as GDAL-supported format

[Tomislav Hengl]

Dear r-sig-geo,

I would like to initiate the processes of registering the SAGA grid  
format under GDAL (the SAGA 2.0.4 source code is available for  
download here  
http://sourceforge.net/projects/saga-gis/files/SAGA%20-%202.0/SAGA%202.0.4/saga_2.0.4_src.zip/download).

Do I have to follow some formal procedure, or do I have to prepare the  
GDAL driver myself (as explained at  
http://www.gdal.org/gdal_drivertut.html)? Any help/suggestions are  
welcome (apparently it should not be too complicated).

SAGA grid consists of tree types of files:
1. "*.sgrd" - the header file with name, data format, XLL, YLL, rows  
columns, cell size, z-factor and no data value;
2. "*.sdat" - the raw data file;
3. "*.hgrd" - the history file;

Here are some examples hot to read and write the SAGA grids to R:

> library(gstat)
> library(RSAGA)
> library(spatstat)

> data(meuse.grid)
> coordinates(meuse.grid) <- ~x+y
> gridded(meuse.grid) <- TRUE
> proj4string(meuse.grid) = CRS("+init=epsg:28992")

# write to SAGA grid format;
> write.asciigrid(meuse.grid["soil"], "meuse_soil.asc")
> rsaga.esri.to.sgrd(in.grids="meuse_soil.asc",  
> out.sgrd="meuse_soil.sgrd", in.path=getwd())

# read SAGA grid format:
> sgrd <- matrix((unlist(strsplit(readLines(file("meuse_soil.sgrd")),  
> split="\t= "))), ncol=2, byrow=T)
> sgrd
       [,1]              [,2]
  [1,] "NAME"            "meuse_soil"
  [2,] "DESCRIPTION"     "UNIT"
  [3,] "DATAFILE_OFFSET" "0"
  [4,] "DATAFORMAT"      "FLOAT"
  [5,] "BYTEORDER_BIG"   "FALSE"
  [6,] "POSITION_XMIN"   "178460.0000000000"
  [7,] "POSITION_YMIN"   "329620.0000000000"
  [8,] "CELLCOUNT_X"     "78"
  [9,] "CELLCOUNT_Y"     "104"
[10,] "CELLSIZE"        "40.0000000000"
[11,] "Z_FACTOR"        "1.000000"
[12,] "NODATA_VALUE"    "-9999.000000"
[13,] "TOPTOBOTTOM"     "FALSE"

# read the raw data: 4bit, numeric (FLOAT), byte order small;
> sdat <- readBin("meuse_soil.sdat", what="numeric", size=4,  
> n=as.integer(sgrd[8,2])*as.integer(sgrd[9,2]))
> sdat.sp <- as.im(list(x=seq(from=as.integer(sgrd[6,2]),  
> length.out=as.integer(sgrd[8,2]), by=as.integer(sgrd[10,2])),  
> y=seq(from=as.integer(sgrd[7,2]), length.out=as.integer(sgrd[9,2]),  
> by=as.integer(sgrd[10,2])), z=matrix(sdat,  
> nrow=as.integer(sgrd[8,2]), ncol=as.integer(sgrd[9,2]))))
> sdat.sp <- as(sdat.sp, "SpatialGridDataFrame")
# replace the mask value with NA's:
> sdat.sp at data[[1]] <-  
> ifelse(sdat.sp at data[[1]]==as.integer(sgrd[12,2]), NA,  
> sdat.sp at data[[1]])
> spplot(sdat.sp)


Of course, it would be much easier to have this in a single line:

> meuse.grid <- readGDAL("meuse_soil.sgrd")

or

> writeGDAL(meuse.grid["soil"], "meuse_soil.sgrd", "SAGA")



PS: A new version of SAGA has just been released few days ago.


thank you,

T. Hengl
http://home.medewerker.uva.nl/t.hengl/

Connected discussion:  
https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20080130/cd5c3748/attachment.pl