Распараллеливание для изменения значений растра

Question

Уважаемое сообщество stackoverflow,

Я прочитал кучу ответов о том, как распараллелить расчет растра. Однако я не нашел ответа, подходящего для моих вопросов. Мой вопрос о распараллеливании для замены значений в растре фреймом данных. Это кажется легким, и я нашел способ сделать это. Тем не менее, у меня есть растр более 5 миллионов ячеек и фрейм данных более 340 000 строк. Таким образом, мой цикл for занимает очень долгое время. Таким образом, я задумался о цикле foreach для выполнения моих расчетов. Моя главная проблема - объединить мои результаты в конце.

Я приведу простой пример (от: http://j.p.rossi.free.fr/rpackages/ecpaysage/TD_Analyse_quanti_paysage.html#(29)), чтобы вы поняли, что я хочу сделать:

library(ECPaysage)
library(raster)
r <- raster(system.file("extdata/r.tif", package="ECPaysage")) 
plot(r)

Это растр землепользования, для которого нам нужна информация о патче. Для этого мы используем пакет SDMTools:

library(SDMTools)
res(r)

matbi <- r
matbi[] <- 0 # let's create the same raster with 0 values
w <- which(r[]==1) 
matbi[w] <- 1 # this creates a binary raster for the specific landuse == 1 

plot(matbi, axes=F, box=F, main="landuse 1 : build") 

Из этого двоичного растра мы хотим создать информацию о патче. Таким образом, мы создаем растр патчей:

matpatch <- ConnCompLabel(matbi)
plot(matpatch, axes=F,box=F, main="patch ID landuse 1 :build") 

Затем мы извлекаем информацию на уровне патча из этого растра патча:

patch <- PatchStat(mat=matpatch, cellsize = res(r)[1], latlon = FALSE)
names(patch)

dim(patch)

patch$patchID
dim(patch) 

Теперь у нас есть количество патчей, содержащихся в растре (каждая ячейка имеет идентификатор патча). И мы хотим извлечь в растре значение патча для каждой ячейки указанного идентификатора патча.

shapeindex <- matpatch

w <- which(shapeindex[]==0) # cells that do not belong to any patch
shapeindex[w] <- NA

Таким образом, чтобы извлечь для каждой ячейки идентификатора патча соответствующий индекс патча, выполняется цикл for:

system.time(for(p in 2:(dim(patch)[1])) {
  w <- which(shapeindex[]==patch$patchID[p])
  shapeindex[w] <- patch$shape.index[p]
})
plot(shape, axes=F,main="shape index - build", box=F, col=rev(terrain.colors(dim(patch)[1]-1))) 

Это прекрасно работает, но с моими данными (растр из более 5 миллионов ячеек и фрейм данных с более чем 340 000 строк) это занимает слишком много времени (и это плохо для воспроизводимости). Одной из моих идей было распараллелить этот цикл с foreach. Тем не менее, моя главная проблема - объединить мои растры в конце. Если бы я распараллелил строки данных, у меня будет 340 000 матриц для объединения. И это может занять то же время, что и цикл for. Таким образом, мне было интересно, есть ли способ (или несколько способов, например, функция), чтобы ускорить процесс и / или легко объединить растры. Спасибо за всю помощь, которую вы можете оказать. Лучший, Adrienne

Answer 1

Спасибо за ваш ответ @ulfelder Я попробовал что-то, как вы сказали. Я создал функцию, использую lapply и комбинирую растры между ними:

Если я возьму предыдущий код перед циклом, это даст:

shape <- matpatch
w <- which(shape[]==0)
shape[w] <- NA
class(shape)
shape
plot(shape)

# create my function for one iteration of the for loop
RastVal <- function(monr,vec,i){
  require(raster)
  w <- which(monr[]==vec$patchID[i])
  x <- which(monr[]!=vec$patchID[i])
  monr[w] <- vec$shape.index[i]
  monr[x] <- 0
  return(monr)
}

# create a list of rasters containing the patch information
RastL <- lapply(X=c(2:(dim(patch)[1])),
                   FUN=function(x)RastVal(shape,patch,x))

# Combining the rasters
new.rast <- RastL[[1]]
for(i in 2:length(RastL)){
  x <- RastL[[i]][] != 0
  new.rast[x] <- RastL[[i]][x]
}
plot(new.rast) # should be the same as the raster shape in the previous code

Однако я не уверен, что это более эффективно ... Может быть, если растр увеличивается, разница во времени между этими двумя способами имеет тенденцию к уменьшению ...

system.time(for(p in 2:(dim(patch)[1])) {
  w <- which(shape[]==patch$patchID[p])
  shape[w] <- patch$shape.index[p]
})

system.time(RastL <- lapply(X=c(2:(dim(patch)[1])),
                   FUN=function(x)RastVal(shape,patch,x)))

new.rast <- RastL[[1]]
system.time(for(i in 2:length(RastL)){
  x <- RastL[[i]][] != 0
  new.rast[x] <- RastL[[i]][x]
})

Более того, создание списка из примерно 340 000 растров может привести к проблемам с памятью ... Еще раз спасибо!

Best, Adrienne