Hi,
a GDAL interface to SAGA grids would of course be very welcome. Right
now the easiest (but indirect) way of importing SAGA grids is with
read.sgrd("mysagagrid")
which uses rsaga.sgrd.to.esri (to a tempfile()) and then read.ascii.grid
In the next RSAGA release I will consider using readBin() in read.sgrd()
in order to read the SAGA grid file directly as you suggested.
Cheers
Alex
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).
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...@data[[1]] <-
ifelse(sdat...@data[[1]]==as.integer(sgrd[12,2]), NA,
sdat...@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@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Alexander Brenning
brenn...@uwaterloo.ca - T +1-519-888-4567 ext 35783
Department of Geography and Environmental Management
University of Waterloo
200 University Ave. W - Waterloo, ON - Canada N2L 3G1
http://www.fes.uwaterloo.ca/geography/faculty/brenning/
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
