Р: Годовой составной растр, основанный на медиане: Как получить индекс исходного слоя для каждого пикселя?

Question

У меня есть список растровых стэков, которые являются временными сериями (более 300 слоев) для разных групп. Временные ряды нерегулярны, и поэтому я хочу создавать ежегодные композиты, основанные на медиане ndvi. Таким образом, из доступных изображений за один год выбирается пиксель со средним значением ndvi (близким к нему). Для других групп я хочу создавать ежегодные композиты на основе этого композита ndvi. Я пытаюсь создать маску, в которой каждый пиксель имеет значение индекса изображения, используемого для медианного композита ndvi. Я буду применять это к другим группам, поэтому у меня есть «одинаковый» годовой композит для каждой группы в год.

К сожалению, я застрял в создании маски. Я создал несколько фиктивных данных и каким-то образом он возвращает два индексатора (один со значениями 1-3, другой 2-4), в то время как я ожидал один (со значениями 1-4). Также моя функция не может обрабатывать значения NA, и добавление na.rm в функцию calc не решает эту проблему. Мне интересно, что мне нужно настроить, чтобы получить один выходной слой со значениями от 1 до 4, и как разрешить функции 'which' работать с NA.

#dummy data:
r <- raster(nrow=5, ncol=5)
set.seed(20181114)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))
names(s) <- c("X2001", "X2002a", "X2002b", "X2002c", "X2002d")
#s$X2002a[2] <- NA 

AnnualMask <- function(ts, year){
  year <- as.character(year)
  ts_year <- subset(ts, (grep(year, names(ts))))
  indexraster <- calc(ts_year, function(x){
    medianval <- median(x, na.rm = TRUE)
    which(abs(x - medianval) == min(abs(x - medianval)))
           })
  return(indexraster)
}

mask2002 <- AnnualMask(s, 2002)
plot(mask2002)

Answer 1

Спасибо за предоставленный пример данных:

library(raster)
r <- raster(nrow=5, ncol=5)
set.seed(20181114)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))
names(s) <- c("X2001", "X2002a", "X2002b", "X2002c", "X2002d")
s[[2]][1:3] <- NA   

Вот две функции, которые делают то же самое. # 1, возможно, проще, но, вероятно, менее эффективно.

annualIndex1 <- function(ts, year) {
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))
    medianval <- calc(ts_year, median, na.rm = TRUE)
    dif <- abs(ts_year - medianval)
    which.min(dif)
}
x1 <- annualIndex1(s, 2002)

Этот использует calc для объединения нескольких шагов

annualIndex2 <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}  
x2 <- annualIndex2(s, 2002)

Результаты одинаковы (и нет NA)

all(values(x1) == values(x2))
#[1] TRUE

Однако, если вы создадите ячейку NA во всех слоях

s[5] <- NA 

функция 2 не работает. Это можно настроить следующим образом

annualIndex2b <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        if (all(is.na(x))) return( NA )
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}

x2b <- annualIndex2b(s, 2002)

Вы используете термин «маска», но на самом деле вы создаете индекс, можете ли вы использовать соответствующие ячейки с strackSelect следующим образом:

ts_year <- subset(ts, (grep(year, names(ts))))
v <- stackSelect(ts_year, x2)

Answer 2

Я пока не уверен насчет части со значениями в диапазоне 1-4, однако здесь есть кое-что, касающееся части which и работы с NA:

# here is some dummy data
x <- c(3,4,5,NA,1,-2-45,2,54)

# now here is the which part as found in your function
medianval <- median(x, na.rm = TRUE)
which(abs(x - medianval) == min(abs(x - medianval)))
# returning in this case
integer(0)

# now let's add na.rm = T in the min function
medianval <- median(x, na.rm = TRUE)
which(abs(x - medianval) == min(abs(x - medianval), na.rm = T))
# returning
[1] 1