Невозможно рассчитать попиксельную регрессию в R на растровом стеке с удовольствием

Я работаю с растрами, и у меня есть RasterStack с 7n слоями. Я хотел бы рассчитать регрессию по пикселям, используя формулу ниже. Я пытался сделать это с помощью raster::calc, но моя функция не удалась с сообщением:

'Ошибка в lm.fit (x, y, offset = offset, singular.ok = singular.ok, ...): 0 (не-NA) случаев.'

Но все растры в порядке и содержат числа (не только NA), я могу построить график и вычислить общую линейную регрессию по формуле

 cr.sig=lm (raster::as.array(MK_trend.EVI.sig_Only) ~ raster::as.array(stack.pet)+raster::as.array(stack.tmp)+raster::as.array(stack.vap)+raster::as.array(stack.pre)+raster::as.array(stack.wet)+raster::as.array(stack.dtr))

Но когда я складываю слои с

allData = stack(MK_trend.EVI.sig_Only,stack.dtr,stack.wet,stack.pre,stack.vap,stack.tmp,stack.pet)

и попробуйте функцию calc

    # Regression Function, R2
lmFun=function(x){
    x1=as.vector(x);
    if (is.na(x1[1])){
        NA 
    } else {
        m = lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7])
        return(summary(m)$r.squared)
    }
}

Я вижу сообщение об ошибке.
Я новичок в R и программировании, так что, может быть, произошла какая-то глупая ошибка? Буду признателен за любой намек, чтобы обработка работала.


r function regression raster calc
person Rinat S.    schedule 22.11.2017    source источник
comment
Можете ли вы предоставить код, который создает пример данных, воспроизводящих ошибку?   -  person Robert Hijmans    schedule 22.11.2017
comment
@RobertH Дорогой Роберт, мне кажется, что в моем наборе данных есть проблема, поэтому я загружаю его в виде ZIP-архива (1,7 ГБ) dropmefiles.com/PDZL0   -  person Rinat S.    schedule 06.12.2017

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете использовать calc для пиксельной (локальной) регрессии, но ваша формула, похоже, предполагает, что вам нужно что-то еще (глобальная модель).

Если бы регрессия была пиксельной, у вас было бы равное количество значений x и y для каждой ячейки, и вы могли бы использовать calc. См. Примеры ?calc.

Вместо этого у вас есть 1 y (независимая) и 6 x (зависимая) значения переменных для каждой ячейки. Это говорит о том, что вам нужна глобальная модель. Для этого вы можете сделать что-то вроде этого:

library(raster)
# example data
r <- raster(nrow=10, ncol=10)
set.seed(0) 
s <- stack(lapply(1:7, function(i) setValues(r, rnorm(ncell(r), i, 3))))  
x <- values(s)

# model
m <- lm(layer.1 ~ ., data=data.frame(x))

# prediction
p <- predict(s, m)

Это требует загрузки всех значений в память. Если вы не можете этого сделать, вы можете взять большую обычную пробу. См. sampleRegular

И чтобы показать, почему ваш подход не работает:

testFun=function(x1){
    lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7])
}

# first cell
v <- s[1]
v
#      layer.1  layer.2   layer.3  layer.4  layer.5  layer.6  layer.7
#[1,] 4.788863 4.345578 -0.137153 3.626695 3.829971 4.120895 1.936597

m <- testFun(v)
m
#Call:
#lm(formula = x1[1] ~ x1[2] + x1[3] + x1[4] + x1[5] + x1[6] + x1[7])

#Coefficients:
#(Intercept)        x1[2]        x1[3]        x1[4]        x1[5]        x1[6]        x1[7]  
#      4.789           NA           NA           NA           NA           NA           NA  


summary(m)$r.squared
# 0

Хотя я не получаю сообщение об ошибке, о котором вы сообщаете (но все значения R ^ 2 равны нулю).

person Robert Hijmans    schedule 23.11.2017
comment
Уважаемый Роберт, большое спасибо за разъяснения, я попытался выйти из перехвата (return (summary (m) $ coefficients [1]), но возникла та же ошибка, что и раньше. - person Rinat S.; 06.12.2017