Re: [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.rowling...@lancaster.ac.uk
To: pisican...@hotmail.com
CC: r-sig-geo@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@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.biv...@nhh.no

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo