[R-sig-Geo] Tiled processing...

[Alexander Brenning]

Hi,

maybe the RSAGA package has the solution to your problem; there are 
actually two ways of applying functions to grids in RSAGA, (1) by 
row-by-row processing (special form of tiles), or (2) using the SAGA 
binaries.

1) local.function

local.function and focal.function are very flexible tools, but they are 
slow because they are written in R. They work with ASCII grids (see e.g. 
write.ascii.grid in RSAGA, or the appropriate GDAL export functions that 
work with your GeoTIFFs).

In your case:

local.function("ingrid", varnames = "outgrid",
      fun = function(x) x + 1000)

(or use focal.function with radius = 0). This will of course work with 
much more general R functions (even with predict methods if you look at 
grid.predict and multi.focal.function).

The RSAGA package technically depends on Windows (because most of its 
functions use SAGA GIS Windows binaries) but the local.function and 
focal.function are (supposed to be) platform-independent; I can provide 
you with the source code if you work on a non-Windows system and can't 
get it from CRAN.

2) rsaga.grid.calculus

Another approach, which involves SAGA itself and therefore (currently) 
depends on Windows, uses SAGA's grid calculator. SAGA is able to process 
large grids efficiently. E.g.:

# first convert the ASCII grid to a SAGA grid (.sgrd):
rsaga.esri.to.sgrd("ingrid")
rsaga.grid.calculus("ingrid", "outgrid", formula = "a+1000")
   # 'ingrid' is treated as 'a' in the formula
rsaga.sgrd.to.esri("outgrid", prec = 2)

I hope this helps...

Cheers
  Alex


Jonathan Greenberg wrote:
> I've recently got back into using R to perform spatial analyses, and I'm 
> trying to figure out how to perform "true" tiled processing, e.g. 
> controlled reading of subsets of an input file, performing a function on 
> this subset, and writing the output, subset by subset, to an output file 
> and, finally, setting up the appropriate "header" info (the metadata).
> 
> Using suggestions Tim Keitt and Roger Bivand gave me some years back, I 
> wrote this simple code (all it should do is take an elevation geotiff 
> and add 1000 to the elevations, writing the output):
> 
> ***
> 
> library(rgdal)
> ds1 <- GDAL.open("elev.tif")
> driver <- new('GDALDriver', 'GTiff')
> ds2 <- new('GDALTransientDataset', driver, nrow(ds1), ncol(ds1), type = 
> "Float32")
>  for (row in nrow(ds1)) {
>    row
>    x <- getRasterData(ds1, offset = c(row - 1, 0), region = c(1, 
> ncol(ds1)))
>    elevnew <- x+1000
>    putRasterData(ds2, elevnew, offset = c(row - 1, 0))
> 
>  }
> saveDataset(ds2, "out.tif")
> closeDataset(ds2)
> closeDataset(ds1)
> 
> ***
> 
> Couple of issues/questions:
> 1) The code doesn't actually seem to work -- I get a tiny, unreadable 
> TIFF out of the back of this algorithm -- what is wrong?
> 2) I'm actually unclear about exactly what this is doing -- am I really 
> just creating a large in-memory matrix (ds2) and then writing the entire 
> output image out to disk, or is this somehow writing line-by-line?  If 
> the former is true, how do I modify this to write within the loop, 
> rather than all-at-once?
> 
> Thanks!
> 
> --j
> 
> 

-- 
Alexander Brenning
brenning at uwaterloo.ca - T +1-519-888-4567 ext 35783
Department of Geography and Environmental Management
University of Waterloo
200 University Ave. W - Waterloo, ON - Canada N2L 3G1
http://www.fes.uwaterloo.ca/geography/faculty/brenning/