Обновление растровых значений в циклах foreach в R?(растровое вложение временных рядов NA) - PullRequest
Новая
       43

Обновление растровых значений в циклах foreach в R?(растровое вложение временных рядов NA)

0 голосов
/ 04 июля 2019

Я пытаюсь вписать значения NA в растр временного ряда. Вот воспроизводимый пример моих данных: 

library(raster)
library(rgdal)
library(doParallel)
library(foreach)

r1 <- r2 <- r3 <- r4 <- r5 <- raster(nrow=100, ncol=100)
values(r1) <- runif(ncell(r1))
values(r2) <- runif(ncell(r2))
values(r3) <- runif(ncell(r3))
values(r4) <- runif(ncell(r4))
values(r5) <- runif(ncell(r5))

s <- stack(r1, r2, r3, r4, r5)
time_series <- brick(s)
time_series[1, 30][2] <- NA
time_series[3, 20][3] <- NA
time_series[5, 10][5] <- NA
time_series[8, 40][4] <- NA

Есть пакеты, такие как gapfill, но я нашел их слишком медленными для моей задачи. Я нашел другой метод здесь, как в ответе: https://gis.stackexchange.com/questions/279354/ndvi-time-series-with-missing-values

по: https://gis.stackexchange.com/users/8520/jeffrey-evans

Я хочу преобразовать цикл for в foreach, чтобы рассчитать его для больших изображений. Вот код с циклами for: 

impute.loess <- function(y, x.length = NULL, s = 0.80, 
                         smooth.data = FALSE, ...) {
  if(is.null(x.length)) { x.length = length(y) }
  options(warn = -1)
  x <- 1:x.length
  if (all(is.na(y))) {
    return(y)
  } else {
    p <- loess(y ~ x, span = s, data.frame(x = x, y = y))
    if(smooth.data == TRUE) {
      y <- predict(p, x)
    } else {
      na.idx <- which( is.na(y) )
      if( length(na.idx) > 1 ) {
        y[na.idx] <- predict(p, data.frame(x = na.idx))
      }
    }   
    return(y)
  }
}

time_series_new <- time_series

time_series_new[] <- NA

for (rl in 1:nrow(time_series)) {
  v <- getValues(time_series, rl, 1)
  time_series_new[rl,] <- as.matrix(t(apply(v, MARGIN=1, FUN=impute.loess)))
}

Альтернатива Foreach, которую я пробовал, такова: 

time_series_new2 <- time_series

time_series_new2[] <- NA

cl <- parallel::makeCluster(detectCores()-1)
doParallel::registerDoParallel(cl)

time_series_new2 <- foreach (rl = 1:nrow(time_series),
                             .packages = "raster",
                             .combine = 'rbind') %dopar% {
                        v <- getValues(time_series, rl, 1)
                        time_series_new[rl,] <- as.matrix(t(apply(v,
                                                MARGIN=1, FUN=impute.loess)))
}

parallel::stopCluster(cl)

Однако вот в чем разница: 

> class(time_series_new)
[1] "RasterBrick"
attr(,"package")
[1] "raster"

> class(time_series_new2)
[1] "matrix"

Если я не назначаю цикл foreach объекту, он просто экспортирует результат. Я хочу обновить растровый объект в конце, но не смог найти решение моей проблемы.

Я не смог найти способ установки значений матрицы для растрового объекта - возможно, заданные значения не работали, так как размеры были такими, как: 

> dim(time_series_new)
[1] 100 100   5

> dim(time_series_new2)
[1] 10000     5

Я знаю, что цикл foreach работает по-другому. Есть ли способ обновить объект time_series_new2 внутри цикла foreach, чтобы в конце я мог получить обновленный растровый объект?

EDIT:

setValues ​​() действительно работает! как: 

time_series_new3 <- time_series

time_series_new3[] <- NA #empty raster object

time_series_new3 <- setValues(time_series_new3, time_series_new2) #filled with matrix rendered from foreach loop

> time_series_new3
class      : RasterBrick 
dimensions : 100, 100, 10000, 5  (nrow, ncol, ncell, nlayers)
resolution : 3.6, 1.8  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      :      layer.1,      layer.2,      layer.3,      layer.4,      layer.5 
min values : 1.468023e-04, 3.525158e-04, 9.689084e-05, 5.349121e-05, 4.214607e-05 
max values :    0.9999564,    0.9999854,    0.9997795,    0.9999780,    0.9997880 

> time_series_new2
class      : RasterBrick 
dimensions : 100, 100, 10000, 5  (nrow, ncol, ncell, nlayers)
resolution : 3.6, 1.8  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      :      layer.1,      layer.2,      layer.3,      layer.4,      layer.5 
min values : 1.468023e-04, 3.525158e-04, 9.689084e-05, 5.349121e-05, 4.214607e-05 
max values :    0.9999564,    0.9999854,    0.9997795,    0.9999780,    0.9997880

> all.equal(time_series_new2, time_series_new3)
[1] TRUE

Тем не менее, я хотел бы знать об обновлении в foreach.

1 Ответ

1 голос
/ 05 июля 2019

В цикле foreach вы не можете себе позволить побочный эффект обновления rasterBrick time_series_new.То есть time_series_new знает, что это - тип объекта raster.Комбинация rbind приведет не-data.frames к matrices.Вот так 100 x 100 x 5 становится 10000 x 5.

Я предполагаю, что вы собираетесь parallel вычисления из-за медленной петли for.Если это так, я рекомендую подходить к проблеме по-другому, особенно если пропущенных значений не так много.

Сначала мы можем увидеть, сколько строк на самом деле содержат пропущенные данные: 

missing_dat_rows <- which(is.na(getValues(time_series)) == T, arr.ind = T)[, 1]
missing_dat_rows <- unique(missing_dat_rows)

missing_dat_rows
#[1]  30 220 740 410

Таким образом, вместо того, чтобы перебирать 10 000 результатов, теперь мы можем сосредоточиться на этих 4 результатах.

time_series3 <- time_series
for (mis_row in missing_dat_rows) {
  values(time_series3)[mis_row, ] <- impute.loess(getValues(time_series3)[mis_row, ])
}

К сожалению, мне не удалось заставить функцию impute.loess() возвращать значения для меня.Я сделал пару небольших больших изменений, которые также могут быть полезны, если вы хотите продолжить ваш цикл: 

impute.loess <- function(y, x.length = NULL, s = 0.80, 
                         smooth.data = FALSE, ...) {
  if(is.null(x.length)) { x.length = length(y) }
  options(warn = -1)

  x <- 1:x.length
  if (all(is.na(y))| all(!is.na(y))) { #added the or statement - I don't think we want to do this if there are no missing values.
    return(y)
  } else {
    p <- loess(y ~ x, span = s, data.frame(x = x, y = y))
    if(smooth.data == TRUE) {
      y <- predict(p, x)
    } else {
      na.idx <- which( is.na(y) )
      # if( length(na.idx) > 1 ) { #commented out - I feel as though we should be replacing all NAs
        y[na.idx] <- predict(p, data.frame(x = na.idx))
      # }
    }   
    return(y)
  }
}
...