[R-sig-Geo] creating a geotiff

[Barry Rowlingson]

On Thu, Aug 13, 2009 at 3:10 PM, Monica Pisica<pisicandru at hotmail.com> 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