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

[Monica Pisica]

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.Bivand at nhh.no
> To: pisicandru at hotmail.com
> CC: b.rowlingson at lancaster.ac.uk; r-sig-geo at 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.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
>
_________________________________________________________________
Get back to school stuff for them and cashback for you.

B_BackToSchool_Cashback_BTSCashback_1x1