Данные вашего примера (немного более краткие)
library(raster)
b <- brick(nrows=100, ncols=100, nl=5)
r.start <- raster(b)
nc <- ncell(b)
set.seed(88)
values(b) <- cbind(runif(nc, min = -10, max = 10), runif(nc, min = -20, max = 10),
runif(ncell(b), min = -30, max = 10), runif(nc, min = -40, max = 10), runif(nc, min = -50, max = 10))
values(r.start) <- round(runif(nc, min = 1, max = 2))
r.end <- r.start + 3
Теперь вы можете сделать:
s <- stack(r.start, r.end, b)
x <- calc(s, fun=function(x) sum(x[(x[1]:x[2])+2]))
И это, кажется, работает
a <- s[1]
a
# layer.1.1 layer.2.1 layer.1.2 layer.2.2 layer.3 layer.4 layer.5
#[1,] 2 5 -1.789974 2.640807 4.431439 -23.09203 -5.688119
fun <- function(x) sum(x[(x[1]:x[2])+2])
fun(a)
#[1] -21.70791
x[1]
#[1] -21.70791
calc для растровых объектов то же, что apply для матриц. (именно поэтому он называется app в terra.
Для начала нужно написать функцию, которая делает то, что вы хотите с вектором.
x <- 1:10
test1 <- function(start, end, values) {
mean(values[start:end])
}
test1(2, 5, x)
test1(5, 8, x)
calc принимает только один аргумент, поэтому такая функция
test2 <- function(values) {
# the +2 to skip the first two elements in the computation
start <- values[1] + 2
end <- values[2] + 2
mean(values[start:end])
}
test2(c(2, 5, x))
test2(c(5, 8, x))
И более краткая версия
test3 <- function(v) {
mean(v[ (v[1]:v[2])+2 ] )
}
test3(c(2, 5, x))
#[1] 3.5
test3(c(5, 8, x))
#[1] 6.5
Второе добавление (и напоминание всегда проверять со значениями NA!). test3 прерывается, когда один из индексов (начало и конец) равен NA (это нормально, если остальные NA)
test3(c(NA, 5, x))
#Error in v[1]:v[2] : NA/NaN argument
Итак, нам нужна функция, которая ловит эти
test4 <- function(v) {
if (any(is.na(v[1:2]))) {
NA
} else {
mean(v[ (v[1]:v[2])+2 ] )
}
}
test4(c(NA, 5, x))
#[1] NA
test4(c(1, 5, x))
#[1] 3
Обычно "start" и "end" оба будут NA одновременно, поэтому более простая версия, которая также должна работать, может быть
test5 <- function(v) {
if (is.na(v[1])) {
NA
} else {
mean(v[ (v[1]:v[2])+2 ] )
}
}
. Этот подход с calc может будьте медленнее, поскольку он превращает RasterBrick в RasterStack с 365 + 2 слоями. Это значительно замедляет чтение данных. Таким образом, вы можете попробовать этот подход вместо overlay (здесь снова используя sum)
f <- function(i, v) {
j <- !is.na(i[,1])
r <- rep(NA, nrow(i))
x <- cbind(i[j,,drop=FALSE], v[j,,drop=FALSE])
r[j] <- apply(x, 1, function(y) sum(y[ (y[1]:y[2])+2 ] ))
r
}
cal <-stack(r.start, r.end)
x <- overlay(cal, b, fun= f, recycle=FALSE)
Вы можете ускорить алгоритм, написав его на Rcpp / C ++
library(Rcpp)
cppFunction('std::vector<double> gtemp(NumericMatrix cal, NumericMatrix wth) {
std::vector<double> out(cal.nrow(), NAN);
for (int i=0; i<cal.nrow(); i++) {
if (!std::isnan(cal(i,0))){
NumericVector v = wth(i,_);
size_t start = cal(i,0)-1;
size_t end = cal(i,1);
out[i] = std::accumulate(v.begin()+start, v.begin()+end, 0.0);
}
}
return out;
}')
x <- overlay(cal, b, fun=gtemp, recycle=FALSE)