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

[Roger Bivand]

On Thu, 13 Aug 2009, Monica Pisica wrote:

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

If you want to extract the bands, they are:

x <- readGDAL("x.tif")
summary(x)
names(x)

the band1, band2, ... in its data slot, so

x$band1

is the first band (as a vector). If you need it as a matrix, you can 
subset to the single band and coerce:

as(x["band1"], "matrix")

Hope this helps,

Roger

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