Как повторить коды, меняющие переменные в последовательности в R

Question

Это код, который я хочу повторить

A_1981 <- Base[1:12]]
B <- sum(A_1981)
MFI_1981 <- sum(A_1981^2)/B

База - растровый кирпич

A_1981 - на год

MFI_1981 - окончательный результат

Так что я должен продолжить со следующего года

A_1982 <- Base[13:24]]
B <- sum(A_1982)
MFI_1982 <- sum(A_1982^2)/B

Чтобы повторить тот же код, я думаю заменить значения только в именах:

a <- seq(1,421,by=12)
b <- seq(12,432,by=12)
c <- seq(1981,2016, by=1)

И сделать это в последовательностив течение следующего третьего года будет что-то вроде этого

A_a[3] <- Base[[b[3]:c[3]]
B <- sum(A_a[3])
MFI_a[3] <- sum(A_[3]^2)/B

Должен быть какой-то способ с for или сделать функцию.Но понятия не имею, с чего начать.

Answer 1

Я думаю, вы ищете что-то вроде этого

Пример данных (48 слоев, то есть 4 "года")

library(raster)
f <- system.file("external/rlogo.grd", package="raster")
Base <- stack(rep(f, 4*4))

Подход 1

f <- function(year) {
    start <- (year-1981) * 12 + 1
    A <- Base[[start:(start+11)]]
    sum(A^2)/sum(A)
}

mfi <- lapply(1981:1984, f)
MFI <- stack(mfi)

Подход 2

for (year in 1981:1984) {
    start <- (year-1981) * 12 + 1
    A <- Base[[start:(start+11)]]
    mfi <- sum(A^2)/sum(A)
    writeRaster(mfi, paste0(year, ".tif"))
}

s <- stack(paste0(1981:1984, ".tif"))

Подход 3, с mapply, как в ответе Руй Баррадаса, но с фиксированным значением, когда Base является RasterBrick (и также включая последний год)

n <- nlayers(Base)
a <- seq(1, n, by = 12)
mfi <- mapply(function(i, j) sum(Base[[i:j]]^2)/sum(Base[[i:j]]), a, a+11)
s <- stack(mfi)

Answer 2

Следующее делает то, что вы хотите, используя mapply и создает только один объект в .GlobalEnv, который я назвал MFI.

Я начинаю с создания вектора Base, так как вы не опубликовали пример набора данных.

set.seed(2469)    # Make the results reproducible

n <- 432
Base <- sample(100, n, TRUE)

step <- 12
b <- seq(1 + step, n, by = step)
a <- seq(1, n - step, by = step)

MFI <- mapply(function(i, j) sum(Base[i:j]^2)/sum(Base[i:j]), a, b)


head(MFI)
#[1] 63.66472 70.54014 67.60567 53.15550 58.71111 65.37008

Другой способ - использовать Map, как предлагает @Parfait в своемкомментарий.

obj <- Map(function(i, j) sum(Base[i:j]^2)/sum(Base[i:j]), a, b)
names(obj) <- paste("MFI", 1980 + seq_along(obj), sep = "_")

obj$MFI_1981
#[1] 63.66472

Обратите внимание, что length(obj) равно 35, и поэтому последний obj равен obj$MFI_2015, а не MFI_2016, как сказано в вопросе.Эту проблему легко решить, набрав n <- 444 в самом начале кода.