Эффективная растровая алгебра R

Question

Для расчета коэффициента отражения для многоспектрального сияющего изображения Приобретенный с помощью линейного сканера, у меня есть 1-рядный массив эталонных значений белого цвета: мультиспектральное изображение A: x строк, y cols, z полос мульти-спектральный опорное изображение В: 1 ряд, у смещ_по_столбцев, г полоса

В простом R-коде в виде массива мы можем установить следующий пример:

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

В этом примере я просто повторю каждую строку b, чтобы массив как а затем вычислить а / б, для чего я делаю:

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

Или используя цикл:

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

Теперь реальные изображения относительно большие, x = 969, y = 640, z = 224 и хотел бы знать, есть ли относительно эффективный способ сделать это используя пакет растровых. То, что я пробовал, очень медленно:

* a и b читаются как растровые кирпичи из файлов envi
* Я рассчитываю B3, как указано выше
* Я конвертирую b3 в растровый кирпич, но он очень медленный и близок к проблемам с памятью:

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)

Вероятно, есть способ работать с b напрямую, избегая b2, b3 и b4 ...

Answer 1

Я пробовал с merge (), но, хотя он работает для воспроизводимого примера, он нереально медленный для реального случая:

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)

Поэтому это не совсем решение.