extraction des résidus de la régression pixel par pixel

Question

J'essaie d'extraire les résidus d'une régression pixel par pixel sur une pile raster de NDVI/précipitation. Mon script fonctionne quand je l'exécute avec une petite partie de mes données. Mais quand j'essaie d'exécuter toute ma zone d'étude j'obtiens: "Erreur dans setValues ​​(out, x): les valeurs doivent être numériques, entiers, logiques ou factorielles"

Le lm fonctionne, puisque je peux extraire les deux pente et intercepter. Je ne peux tout simplement pas extraire les résidus.

Une idée de la façon dont cela pourrait être résolu?

Voici mon script:

setwd("F:/working folder/test") 
gimms <- list.files(pattern="*ndvi.tif") 
ndvi <- stack(gimms) 
precip <- list.files(pattern="*pre.tif") 
pre <- stack(precip) 
s <- stack(ndvi,pre) 

residualfun = function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[1:6] ~ x[7:12], na.action=na.exclude) 
r <- residuals.lm(m) 
return (r)}} 

res <- calc(s,residualfun) 

Et voici mes données: https://1drv.ms/u/s!AhwCgWqhyyDclJRjhh6GtentxFOKwQ

Answer 1

Votre fonction ne test si la première couche montre NA valeurs pour éviter l'ajustement du modèle. Mais il peut y avoir NA dans d'autres couches. Vous le savez parce que vous avez ajouté na.action = na.exclude dans votre ajustement lm.
Le problème est que si le modèle supprime certaines valeurs à cause de NA, les résidus n'auront que la longueur des valeurs non-NA. Cela signifie que le vecteur r résultant aura des longueurs différentes en fonction de la quantité de valeurs NA dans les couches. Ensuite, calc ne peut pas combiner des résultats de longueurs différentes dans une pile et un nombre défini de couches.
Pour éviter cela, vous devez spécifier la longueur de r dans votre fonction et attribuer des résidus uniquement aux valeurs non-NA.
Je propose la fonction suivante qui fonctionne maintenant sur l'ensemble de données fourni. J'ai ajouté (1) la possibilité de comparer plus de couches de chaque si vous voulez prolonger votre exploration (avec nlayers), (2) évitez d'ajuster le modèle s'il n'y a que deux valeurs à comparer dans chaque couche (modèle parfait), (3) a ajouté un try si pour une raison quelconque le modèle peut s'adapter, cela produira des valeurs de -1e32 facilement trouvables pour d'autres tests.

library(raster) 
setwd("/mnt/Data/Stackoverflow/test") 
gimms <- list.files(pattern="*ndvi.tif") 
ndvi <- stack(gimms) 
precip <- list.files(pattern="*pre.tif") 
pre <- stack(precip) 
s <- stack(ndvi,pre) 

# Number of layers of each 
nlayers <- 6 

residualfun <- function(x) { 
    r <- rep(NA, nlayers) 
    obs <- x[1:nlayers] 
    cov <- x[nlayers + 1:nlayers] 
    # Remove NA values before model 
    x.nona <- which(!is.na(obs) & !is.na(cov)) 
    # If more than 2 points proceed to lm 
    if (length(x.nona) > 2) { 
     m <- NA 
     try(m <- lm(obs[x.nona] ~ cov[x.nona])) 
     # If model worked, calculate residuals 
     if (is(m)[1] == "lm") { 
     r[x.nona] <- residuals.lm(m) 
     } else { 
     # alternate value to find where model did not work 
     r[x.nona] <- -1e32 
     } 
    } 
return(r) 
} 

res <- calc(s, residualfun)