Using RGDAL to "copy" header information...
On Wed, 5 Nov 2008, Jonathan Greenberg wrote:
R-geographers...
I'm trying to solve a problem to implement a line-by-line tiled processing
using RGDAL (read 1 line of an image, process the one line, write one line of
the image to a binary file). Everything except for the final step I'm able
to do using a combination of RGDAL and r-base commands. Below is the basic
structure, the input file "elev" can be any image file GDAL supports -- the
code below just "copies" the image one line at a time:
###
library(rgdal)
infile='elev'
outfile_base='testout'
outfile_ext='.bil'
outfile=paste(outfile_base,outfile_ext,sep='')
outcon <- file(outfile, "wb")
infile_info=GDALinfo(infile)
nl=infile_info[[1]]
ns=infile_info[[2]]
for (row in 1:nl) {
templine <- readGDAL(infile,region.dim=c(1,ns),offset=c(row-1,0))
writeBin(templine[[1]], outcon,size=4)
}
close(outcon)
# Below doesn't work
# writeGDAL(templine,outfile_base,drivername='EHdr',type="Float32")
###
The issue is this: I need to be able to effectively copy the
geographic/header information from the input file (the last line almost
works, but as you can see it sets the line #s to 1, not "ns"), and parse it
over to the output file (which is some form of a flat binary file -- in this
case, an Arc Binary Raster, but an ENVI binary output would also be good) --
ideally I'd like to be able to modify this header, particularly the number
type (e.g. if the input file is an integer, and I'm writing out a floating
point binary, I need that reflected in the output header). I can do this by
*manually* creating the output header via a series of ascii-writes, but I was
wondering if there is a more effective way of doing this using RGDAL that
might apply generically to the header of any binary image file I might write?
Thanks!
It is the internals of writeGDAL() and create2GDAL() without the bands:
data(meuse.grid)
gridded(meuse.grid) <- ~x+y
d.drv <- new("GDALDriver", "EHdr")
gp <- gridparameters(meuse.grid)
cellsize <- gp$cellsize
offset <- gp$cellcentre.offset
dims <- gp$cells.dim
nbands <- 1
tds.out <- new("GDALTransientDataset", driver = d.drv, rows = dims[2],
cols = dims[1], bands = nbands, type = "Float32")
gt <- c(offset[1] - 0.5 * cellsize[1], cellsize[1], 0.0,
offset[2] + (dims[2] -0.5) * cellsize[2], 0.0, -cellsize[2])
.Call("RGDAL_SetGeoTransform", tds.out, gt, PACKAGE = "rgdal")
list.files(tempdir())
fn <- tempfile()
saveDataset(tds.out, fn)
list.files(tempdir())
GDAL.close(tds.out)
list.files(tempdir())
However, the driver does commit the data file too, so you'd need to run
the above code to produce the header so as not to overwrite your output
data. The key thing is to get the input grid parameters right, but you've
got them anyway. If you want the projection information out too, look in
create2GDAL for the .Call() you need. Note that passing unchecked
variables through .Call may crash your session - going through
GridTopology should make sure that they are stored correctly.
Hope this helps,
Roger
--j
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