[gdal-dev] Re: [R-sig-Geo] Using RGDAL to "copy" header information...

[Roger Bivand]

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