ошибка raster :: cal c () в calctest (x [1: 5]) с пользовательской функцией - PullRequest
Новая
       3

ошибка raster :: cal c () в calctest (x [1: 5]) с пользовательской функцией

0 голосов
/ 17 июня 2020

Я рассчитываю максимальный климатологический дефицит воды из этого исследования / науки c публикации , используя Rcode, который авторы опубликовали вместе с исследованием, доступным здесь

Это очень простой код, в котором в качестве входных данных используются растры ежемесячных осадков, которые я загрузил из продукта Terraclimate. После вычитания 100 мм / месяц (эвапотранспирация) из каждого растра осадков (месяц) создается стек, и следующая функция используется в cal c () 

wd = stack(month.rainfall)-100 # 100 is the evapotranspiration in mm/month

# MCWD Function
mcwd.f = function(x){
result= as.numeric(x)
for(i in 1:length(result)){
wdn = result[i]
wdn1 = result[i-1]

if(i==1){
  if(wdn>0){ result[i]=0}
  else{result[i]=wdn}
}

if(i!=1){
  cwd = wdn1+wdn
  if( cwd < 0){ result[i]=cwd}
  else{result[i]=0}
   }
 }

return(result)  
 }


# Applying the Function
cwd = calc(wd, fun = mcwd.f)

Однако после того, как я запустил строку при запуске cwd я получаю ошибку Ошибка в .calcTest (x [1: 5], fun, na.rm, forcefun, forceapply): не могу использовать эту функцию

Почему у меня эта ошибка? Как мне это исправить?

Я просмотрел похожие сообщения - 1 , 2 и использовал na.rm = TRUE et c, но я все еще получить ту же ошибку. Я также вижу несколько сообщений, в которых говорится, что cal c () не очень хорошо, но я не знаю, как решить эту проблему или найти альтернативу, поскольку я использую код, опубликованный авторами исследования, которое я копирую.

РЕДАКТИРОВАТЬ - Я использовал вектор из 100 случайных чисел (без НАП), чтобы протестировать функцию mcwd.f, и она работает. Поэтому я думаю, что это связано с наличием пикселей NA, которые, как я думал, должны быть решены, если я использую na.rm = TRUE в cal c (), но это не решает его.

1 Ответ

0 голосов
/ 18 июня 2020

Всегда включайте минимальный, самодостаточный, воспроизводимый пример, когда вы задаете вопрос R. Вы говорите, что ваша функция «работает», но вы не предоставляете входные и ожидаемые выходные данные.

Я получаю это с помощью вашей функции 

mcwd.f(c(-5:5))
# [1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

Вот гораздо более краткая версия вашей функции 

mcwd.f1 <- function(x) {
    x[1] <- min(0, x[1])
    for(i in 2:length(x)){
        x[i] <- min(0, sum(x[c(i-1, i)]))
    }
    return(x)  
}
mcwd.f1(c(-5:5))
# [1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

И вариант, который помогает объяснить следующий 

mcwd.f1a <- function(x) {
    x <- c(0, x)
    for(i in 2:length(x)){
        x[i] <- min(0, sum(x[c(i-1, i)]))
    }
    return(x[-1])  
}

векторизованная версия ( после Габора Гротендика ) 

mcwd.f2 = function(x) {
    f <- function(x, y) min(x + y, 0)
    Reduce(f, x, 0, accumulate = TRUE)[-1]
}
mcwd.f2(c(-5:5))
#[1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

Теперь используйте функцию с RasterStack и calc

Пример данных : 

library(raster)
slogo <- stack(system.file("external/rlogo.grd", package="raster")) 
s <- slogo-100

Используйте cal c: 

x <- calc(s, mcwd.f1)
y <- calc(s, mcwd.f2)

Это работает для обеих функций (а также для вашего mcwd.f). Попробуйте с NA значениями 

s[1000:3000] <- NA 
x <- calc(s, mcwd.f1)
y <- calc(s, mcwd.f2)

Тоже работает. 

x
#class      : RasterBrick 
#dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source     : memory
#names      :  red, green, blue 
#min values : -100,  -200, -300 
#max values :    0,     0,    0

Но не для вашей функции 

y <- calc(s, mcwd.f)
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

И вы можете понять, почему, когда вы делаете 

mcwd.f(c(-5:5,NA))
#Error in if (cwd < 0) { : missing value where TRUE/FALSE needed

Чтобы заставить вашу функцию работать (но вы не должны потому что это неэффективно) вы можете добавить первую строку вроде этой 

mcwd.f = function(x){
    if (any(is.na(x))) return(x*NA)
    ...

В предположении, что при наличии одного NA все значения в ячейке равны или должны стать NA

...