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

From: Roger Bivand <Roger.Bivand_at_nhh.no>
Date: Thu, 06 Nov 2008 09:10:29 +0100 (CET)

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

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Thu 06 Nov 2008 - 08:15:45 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 06 Nov 2008 - 09:30:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive