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

[Roger Bivand]

On Fri, 16 Oct 2009, Tomislav Hengl wrote:

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

Tom,

Writing a read-only driver would be the first step. Once that works, it 
could be extended to copy creation, and finally to direct creation if 
justified. If SAGA I/O has a separate I/O library, it could be shipped 
with the driver as in the pcraster driver, but I don't see one in the 
source tree. Writing a driver for SAGA rasters would permit lots of other 
software to interact with SAGA, not just in R, so would be more robust.

You'd develop on a local copy of the GDAL source, and once done submit a 
file archive of the directory to Frank Warmerdam or other GDAL people. 
Maybe ask on the GDAL developers' list if they have seen questions about 
SAGA rasters before. Are SAGA rasters similar to any existing supported 
formats (Erdas?, IDRISI?, ILWIS?) - could an existing driver be copied and 
modified?

Roger

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

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