Efficient R Raster algebra

Question

In order to calculate reflectance for a multi-spectral radiance image acquired with a linear scanner, I have a 1-row array of white reference values: multi-spectral image A: x rows, y cols, z bands multi-spectral reference image B: 1 row, y cols, z bands

In a simple R code in array form we can set the following example:

a <- array(runif(3*5*6),dim=c(3,5,6)) 
b <- matrix(runif(5*6),nrow=5)

In this example I would just repeat each row of b to make an array as a and then calculate a/b, for which I do:

b2 <- rep(b,nrow(a))
dim(b2) <- c(ncol(a),dim(a)[3],nrow(a))
b3 <- aperm(b2,c(3,1,2))
res <- a/b3

Or using a loop:

res <- a
for(i in 1:nrow(a)){
  res[i,,] <- a[i,,]/b
}

Now the real images are relatively large, x=969, y=640, z=224 and would like to know if there is a relatiely efficient way of doing this using package raster. What I have tried is very slow:

*a and b are read as a raster bricks from envi files
*I calculate b3 as above
*I convert b3 to raster brick, but it is very slow and close to memory problems:

x=969; y=640; z=224
b <- matrix(runif(y*z),nrow=y)
b2 <- rep(b,x)
dim(b2) <- c(y,z,x)
b3 <- aperm(b2,c(3,1,2))
b4 <- brick(b3)
extent(b4) <- c(0,y,0,x)
b4
writeRaster(b4,"b4",overwrite=TRUE)
rm(b2,b3,b4)
b4 <- brick("b4.gri")
b4
res <- overlay(x=RadIma,y=b4,fun=function(x,y){x/y}, filename="res",overwrite=TRUE)

There is probably a way operating with b directly thus avoiding b2,b3 and b4...

I've also tried using merge(), but the merging is impractically slow for the real problem — user2955884
Not entirely sure if I understand what you are trying to do, but does it come down to just dividing two rasters? Also, I don't think your rasters are particularly big, especially not with a modern PC, so most of it should be possible to do in memory without I/O — Val
The problem is not the division, the problem is creating the brick for the denominator by repeating 1 row. — user2955884

user2955884 user2955884 · Accepted Answer · 2018-11-06T08:48:53

I've tried with merge(), but while it works for the reproducible example, it is impractically slow for the real case:

a <- array(runif(3*5*6),dim=c(3,5,6)) 
b <- matrix(runif(5*6),nrow=5)
dim(b) <- c(1,dim(a)[2:3])

a <- brick(a)
extent(a) <- c(0,ncol(a),0,nrow(a))
res(a) <- 1
b <- brick(b)
extent(b) <- c(0,dim(b)[2],0,1)
res(b) <- 1
a
b

v <- vector(length=nrow(a),mode="list")
v[1:nrow(a)] <- list(b)
for(i in 1:nrow(a)) extent(v[[i]]) <- c(0,ncol(a),i-1,i)
b4 <- do.call(merge,args=c(v,filename="b4",format="GTiff",overwrite=TRUE))
b4
dim(a)
dim(b4)
res <- overlay(x=a,y=b4,fun=function(x,y){x/y}, filename="res",overwrite=TRUE)

Therefore this is not really a solution.

Efficient R Raster algebra

1 Answers