[R-sig-Geo] GRASS and R question

[Tim Keitt]

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

On 9/9/06, Jonathan Greenberg <jgreenberg at arc.nasa.gov> wrote:
> I was wondering if I could get some brief feedback on exactly what are the
> capabilities of the sp and spgrass6 -- say I'd like to use R to calculate
> and NDVI image (I know GRASS's mapcalc can do this fine, but I'm thinking of
> some more complicated analyses that aren't easily performed in R and I
> figured if I can figure out NDVI I can figure out the rest of them) -- Does
> sp, spgrass6 and R itself provide both the data ingestion capabilities as
> well as the necessary subsetting capabilities if my file is very large (e.g.
> Can R do a per-pixel or per-line read/process/write routine, or will it
> attempt to load an entire raster into memory)?
>
> If R can deal with import of a multiband image and tiled processing, can
> someone here perhaps post an example of how to accomplish the above task
> (taking a 4 band large file, calculating NDVI and then writing out the NDVI
> to disk) using R?  Thanks!
>
> --j
>
>
> --
> Jonathan A. Greenberg, PhD
> NRC Research Associate
> NASA Ames Research Center
> MS 242-4
> Moffett Field, CA 94035-1000
> Office: 650-604-5896
> Cell: 415-794-5043
> AIM: jgrn307
> MSN: jgrn307 at hotmail.com
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


-- 
Timothy H. Keitt
http://www.keittlab.org/