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

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

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

r median mask na raster
person Louise    schedule 15.11.2018    source источник

Ответы (2)


arrow_upward
0
arrow_downward

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

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)
person Robert Hijmans    schedule 17.11.2018

arrow_upward
0
arrow_downward

Я еще не уверен насчет части со значениями в диапазоне от 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
person niko    schedule 15.11.2018