Синхронизировать NA среди слоев растрового стека

Я пытаюсь разработать функцию для "синхронизации" NA между слоями растрового стека, т.е. чтобы убедиться, что для любого данного пикселя стека, если один слой имеет NA, все слои должны быть установлены в NA для этого пикселя.,

Это особенно полезно при объединении растров, поступающих из разных источников для моделирования распределения видов, потому что некоторые модели не обрабатывают правильно NA.

Я нашел два способа сделать это, но я не нахожу ни одного из них удовлетворительным. Один из них требует использования функции getValues и, следовательно, не может использоваться для очень больших стеков или компьютеров с низким объемом оперативной памяти. Другой более безопасный для памяти, но намного медленнее. Поэтому я здесь, чтобы спросить, есть ли у кого-нибудь идея улучшить мои попытки?

Вот две возможности:

  1. Использование getValues ​​()

    syncNA1 <- function (x) 
{
  val <- getValues(x)
  NA.pos <- unique(which(is.na(val), arr.ind = T)[, 1])
  val[NA.pos, ] <- NA
  x <- setValues(x, val)
  return(x)
}

  2. Использование calc()

    syncNA2 <- function(y)
{
  calc(y, na.rm = T, fun = function(x, na.rm = na.rm)
  {
    if(any(is.na(x)))
    {
      rep(NA, length(x))
    } else
    {
      x
    }
  })
}

Теперь демонстрация их соответствующего времени вычислений для того же стека:

> system.time(
+ b1 <- syncNA1(a1)
+ )
   user  system elapsed 
   3.04    0.15    3.20 
> system.time(
+ b2 <- syncNA2(a1)
+ )
   user  system elapsed 
   5.89    0.19    6.08

Большое спасибо за вашу помощь,

Борис

Источник

3 ответа

Решение

С именем стекаs"Я бы сначала использовал calc(s, fun = sum) вычислить слой маски, который записывает местоположение всех ячеек со значением NA, по крайней мере, в одном из слоев стека. mask() затем позволит вам применить эту маску к каждому слою в стеке.

Вот пример:

library(raster)

## Construct reproducible data! (Here a raster stack with NA values in each layer) 
m <- raster(ncol=10, nrow=10)
n <- raster(ncol=10, nrow=10)
m[] <- runif(ncell(m))
n[] <- runif(ncell(n)) * 10
m[m < 0.5] <- NA
n[n < 5] <- NA
s <- stack(m,n)

## Synchronize the NA values
s2 <- mask(s, calc(s,fun = sum))

## Check that it worked
plot(s2)

Источник

Я не знаю насчет скорости, но вы можете попробовать преобразовать в массив, загрузить NA и выполнить обратное преобразование. псевдокод:

xarray<-as.array(xstack)
ind.na<-which(is.na(xarray),array.ind=TRUE)
for(j in nrow(ind.na) ) {
    xarray[ind.na[j,1],ind.na[j,2],]<-NA
   }
nastack<-raster(xarray)

Я не проверил правильный выбор индексов там, и я не подтвердил, что я преобразовал обратно в raster stack правильно, но я надеюсь, вы поняли идею.

РЕДАКТИРОВАТЬ: я провел тест времени, с растрами 1000x1000, но в остальном, как Джош создал.

 microbenchmark(josh(s),syncNA1(s),syncNA2(s),times=5)
Unit: milliseconds
       expr       min        lq    median        uq        max
    josh(s)  774.2363  789.1653  800.2511  806.5364   809.9087
 syncNA1(s)  652.3928  659.8327  692.3578  695.8057   743.9123
 syncNA2(s) 7951.3918 8291.7917 8604.2226 8606.3432 10254.4739
 neval
     5
     5
     5
Источник

В итоге я создал гибридную функцию между syncNA1 и решением Джоша. Эта функция защищает память, если на компьютере недостаточно оперативной памяти, но может работать быстрее, если на компьютере достаточно оперативной памяти:

synchroniseNA <- function(x)
{
  if(canProcessInMemory(x, n = 2))
  {
    val <- getValues(x)
    NA.pos <- unique(which(is.na(val), arr.ind = T)[, 1])
    val[NA.pos, ] <- NA
    x <- setValues(x, val)
    return(x)
  } else
  {
    x <- mask(x, calc(x, fun = sum))
    return(x)
  }
}

Тем не менее, я эмпирически определил, что количество оперативной памяти, используемой data.frame в два раза больше растрового файла (для n аргумент canProcessInMemory()), но я не совсем уверен, что я здесь.

Источник
Другие вопросы по тегам