[R-sig-Geo] creating a geotiff - and a new question

[Monica Pisica]

WOW! thanks so much for the code. That really makes my life easier. I really appreciate it.

Well, i have a new question - and this time it should be more than easy .... i actually feel ashamed, but for whatever reason i get an error.
 
I am reading a multiband image and i want to get tables out of each band.
 
so:
imag <- new("GDALReadOnlyDataset", pca_list)
 
dim(imag)
[1] 1408 1548    6
 
pcab <- list()
 
for(i in 1:6) {
timag<- getRasterTable(imag, band = i)
pcab[[i]] <- timag$band1
}
 
Error: subscript out of bounds
 
i
[1] 2

So it is clear to me that for band = 1 everything is fine but for band = 2 things are not anymore. So what i am doing wrong??
 
Again sorry for this very basic question, and lots of thanks for the code ;-)
 
Monica

----------------------------------------
> Date: Thu, 13 Aug 2009 15:28:51 +0100
> Subject: Re: [R-sig-Geo] creating a geotiff
> From: b.rowlingson at lancaster.ac.uk
> To: pisicandru at hotmail.com
> CC: r-sig-geo at stat.math.ethz.ch
>
> On Thu, Aug 13, 2009 at 3:10 PM, Monica Pisica wrote:
>>
>> Thanks Barry!
>>
>> I will look into it - for now i have my clumsy code and because my dataset is relatively small it does not matter, but for future this might not be the case.
>
> Found it!
>
> writeMgdal <-
> function (dataset, fname, drivername = "GTiff", type =
> "Float32",options=NULL)
> {
> require(rgdal)
> if (nchar(fname) == 0)
> stop("empty file name")
> tds.out <- m2GDAL(dataset = dataset, drivername = drivername,
> type = type, options=options)
> saveDataset(tds.out, fname, options = options)
> invisible(fname)
> }
>
>
> m2GDAL <-
> function (dataset, drivername = "GTiff", type = "Float32", options=NULL)
> {
> #if (is.na(match(type, .GDALDataTypes)))
> # stop(paste("Invalid type:", type, "not in:", paste(.GDALDataTypes,
> # collapse = "\n")))
>
> dataset$z = dataset$z[,ncol(dataset$z):1]
> xcs = diff(range(dataset$x))/ncol(dataset$z)
> ycs = diff(range(dataset$y))/nrow(dataset$z)
> gp = list(cellsize=c(xcs,ycs),
> cells.dim=c(nrow(dataset$z),ncol(dataset$z)),
> cellcentre.offset=c(min(dataset$x),min(dataset$y)))
> cellsize = gp$cellsize
> offset = gp$cellcentre.offset
> dims = gp$cells.dim
> d.drv = new("GDALDriver", drivername)
> nbands = 1
> if (!is.null(options) && !is.character(options))
> stop("options not character")
> tds.out = new("GDALTransientDataset", driver = d.drv, rows = dims[2],
> cols = dims[1], bands = nbands, type = type, options = options,
> handle = NULL)
> gt = c(offset[1] - 0.5 * cellsize[1], cellsize[1], 0, offset[2] +
> (dims[2] - 0.5) * cellsize[2], 0, -cellsize[2])
> .Call("RGDAL_SetGeoTransform", tds.out, gt, PACKAGE = "rgdal")
> p4s <- NA
> if (!is.na(p4s) && nchar(p4s)> 0)
> .Call("RGDAL_SetProject", tds.out, p4s, PACKAGE = "rgdal")
> for (i in 1:nbands) {
> band = as.matrix(dataset$z)
> if (!is.numeric(band))
> stop("Numeric bands required")
> putRasterData(tds.out, band, i)
> }
> tds.out
> }
>
> Usage:
>
>> xyz=list(x=seq(0,1,len=10),y=seq(1,2,len=14),z=matrix(runif(10*14),10,14))
>> image(xyz)
>> writeMgdal(xyz,"xyz.gtiff")
>
> The resulting gtiff file when loaded into Quantum GIS looks the same
> as the image in R, once I've set Qgis to scale the data.
>
> Obviously this only works for single-band data.
>
> Barry
_________________________________________________________________


_sync:082009