This is much better ..... although i still don't understand why my code gave an error, although i have to recognize that yours is much more elegant. Again thanks so much for your help, Monica
---------------------------------------- > Date: Thu, 13 Aug 2009 18:20:11 +0200 > From: roger.biv...@nhh.no > To: pisican...@hotmail.com > CC: b.rowling...@lancaster.ac.uk; r-sig-geo@stat.math.ethz.ch > Subject: Re: [R-sig-Geo] creating a geotiff - and a new question > > 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 > _________________________________________________________________ Get back to school stuff for them and cashback for you. B_BackToSchool_Cashback_BTSCashback_1x1 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo