[R-sig-Geo] Writing a tif using rgdal

[Roger Bivand]

On Tue, 7 Feb 2006, Tim Keitt wrote:

> I seem to remember writing a wrapper for the gdal copy function,
> although none of that has be extensively tested. If you copy an existing
> dataset to a writable dataset, you could then retrieve the data, operate
> on it and write back to the new dataset. It should (?) then have all the
> metadata copied from the original. 

Yes, copyDataset() is there, but won't it expect the type= argument to 
suit the data? My route as below lets you set up a transient dataset and 
then save it to a file, but I hadn't got the:

.Call("RGDAL_SetGeoTransform", dataset, numericvectorof6, PACKAGE="rgdal")

to set the starting point and cell sizes. I'll carry on trying.

> Sorry, I didn't see in your code
> where you were overwriting gdal pointers.
> 

No trouble - I think this hits me when I'm hurrying and I re-use a 
not-closed handle, filling it with another handle. I can't replicate now.

Roger

> THK
> 
> On Tue, 2006-02-07 at 20:19 +0100, Roger Bivand wrote:
> > On Tue, 7 Feb 2006, Andy Bunn wrote:
> > 
> > > Howdy List:
> > > 
> > > I'm reading in a large number of satellite images and analyzing them in R. I
> > > know how to read a tif into R using rgdal and manipulate it. Is there a good
> > > way to write a new tif out using the projection data from the original?
> > > 
> > > So in the example below how can I write out z as a jpg using the metadata
> > > and such from x?
> > > 
> > > -Andy
> > > 
> > > 
> > > library(rgdal)
> > > logo <- system.file("pictures/Rlogo.jpg", package="rgdal")[1]
> > > x <- new("GDALReadOnlyDataset", logo)
> > > displayDataset(x)
> > > y <- getRasterData(x)
> > > z <- rowMeans(y, dims=2)
> > > dim(z)
> > 
> > Taking a deep breath and remembering to do save.image() frequently 
> > [overwriting GDAL handles can give seg.faults]:
> > 
> > library(rgdal)
> > logo <- system.file("pictures/Rlogo.jpg", package="rgdal")[1]
> > x <- new("GDALReadOnlyDataset", logo)
> > displayDataset(x)
> > y <- getRasterData(x)
> > z <- rowMeans(y, dims=2)
> > dim(z)
> > GDAL.close(x)
> > 
> > tif_driver <- new("GDALDriver", "GTiff")
> > getDriverLongName(tif_driver)
> > tif2 <- new("GDALTransientDataset", tif_driver, 77, 101, 1, 'Float64')
> > bnd1 <- putRasterData(tif2, z)
> > displayDataset(tif2)
> > tif_file <- tempfile()
> > saveDataset(tif2, tif_file)
> > GDAL.close(tif2)
> > GDAL.close(tif_driver)
> > tif3 <- GDAL.open(tif_file)
> > displayDataset(tif3)
> > GDAL.close(tif3)
> > 
> > $ gdalinfo zlogo.tif
> > Driver: GTiff/GeoTIFF
> > Size is 101, 77
> > Coordinate System is `'
> > Corner Coordinates:
> > Upper Left  (    0.0,    0.0)
> > Lower Left  (    0.0,   77.0)
> > Upper Right (  101.0,    0.0)
> > Lower Right (  101.0,   77.0)
> > Center      (   50.5,   38.5)
> > Band 1 Block=101x10 Type=Float64, ColorInterp=Gray
> > 
> > (I was playing with filename zlogo.tif, it reads into GRASS OK, but 
> > that's using the same GDAL drivers.)
> > 
> > Well, it works, though it's a bit too exciting dodging the seg.faults 
> > coming from over-writing live object handles while experimenting. For 
> > determined users, deleteDataset() works too (oh yes!). I'll try to wrap it 
> > up to be a bit smoother. No world file yet, but I guess that could just be 
> > copied? Suggestions welcome.
> > 
> > Roger
> > 
> > > 
> > > > version
> > >          _
> > > platform i386-pc-mingw32
> > > arch     i386
> > > os       mingw32
> > > system   i386, mingw32
> > > status
> > > major    2
> > > minor    2.0
> > > year     2005
> > > month    10
> > > day      06
> > > svn rev  35749
> > > language R
> > > 
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > R-sig-Geo at stat.math.ethz.ch
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > 
> > 
> 

-- 
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