[R-sig-Geo] GRASS and R question

[Roger Bivand]

On Sat, 9 Sep 2006, Tim Keitt wrote:

> 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.)

For reading, if the GDAL GRASS plugin is available, you can read the GRASS 
raster as a single band as Tim suggests in rgdal using offset= and 
region.dim=. The plugin doesn't have a write functionality, though.

If spgrass6 is a feasible route, you'd use g.region to imitate the 
subwindow choice, perhaps setting up sequences of regions and storing them 
first - I guess you could use loops and paste() in R to construct the 
regions if you like. Output would also be for the regions, probably to be 
r.patch-ed together later. Now both reading from and writing to GRASS use 
r.*.bin and temporary files.

> 
> 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.
> 

This will be a useful avenue if you are outside GRASS, I'm not sure if the 
plugin GRASS drive is blocked/tiled. 

Roger

> 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
> >
> 
> 
> 

-- 
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