GRASS and R question

Jonathan,

This is what rgdal was designed to do. Here's an example,

ds1 <- GDAL.open("input.tif")
driver <- new('GDALDriver', 'GTiff')
ds2 <- new('GDALTransientDataset', driver, nrow(ds1), ncol(ds1), type
= "Float64")
for (row in nrow(ds1)) {

  x <- getRasterData(ds1, offset = c(row - 1, 0), region = c(1, ncol(ds1)))
  ndvi <- (x[,2] - x[,1]) / (x[,2] + x[,1])
  putRasterData(ds2, ndvi, offset = c(row - 1, 0))

}
saveDataset(ds2, "result.tif")

This scales to absurd image sizes as only a single scanline is loaded
into memory. (Note: I didn't try this code, but it should be about
right.)

I'm in the process of adding an 'apply' function to rgdal, so pretty
soon you'll be able to do:

ds2 <- blockApply(ds1, function(x) (x[,2]-x[,1]) / (x[,2]+x[,1]))

The apply function will automatically use the optimal block size
supplied by GDAL and so should be really fast on tiled images.

For smaller images, it may be easier to use 'sp' classes. My approach
is to use the rgdal classes for big images and then convert reduced
resolution versions to 'sp' for display.

THK

GRASS and R question

Thread (3 messages)