R help - Nov 2008 - Using RGDAL to "copy" header information...

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!

--j


-- 

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307

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