2017-10-17 1 views
1

NDVI/강수의 래스터 스택에서 픽셀별로 회귀 분석에서 잔차를 추출하려고합니다. 내 스크립트는 내 데이터의 작은 부분으로 실행할 때 작동합니다. 하지만 내 연구 영역 전체를 실행하려고 할 때 : "setValues ​​(out, x) 오류 : 값은 숫자, 정수, 논리 또는 요인이어야합니다."픽셀 회귀에 의한 픽셀의 잔차 추출

lm은 두 기울기를 모두 추출 할 수 있기 때문에 작동합니다. 및 가로 채기. 나는 잔류 물을 추출 할 수 없다.

어떻게 수정 될 수 있는지 알고 싶습니다. 여기

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) 

그리고 내 데이터는 다음과 같습니다 :

여기 내 스크립트입니다 https://1drv.ms/u/s!AhwCgWqhyyDclJRjhh6GtentxFOKwQ

+0

흠, 아마 수식 인수에 문제가있다? 명시 적으로 열을 호출합니다. –

+0

데이터에 대한 링크가 더 이상 사용되지 않습니다. – Borealis

답변

1

함수 첫 번째 레이어 모델을 피팅 피하기 위해 NA 값을 보여줍니다 경우에만 테스트. 다른 레이어에는 NA이있을 수 있습니다. 아시다시피 lmna.action = na.exclude을 추가했기 때문입니다.
NA으로 인해 모델에서 일부 값을 제거하는 경우 잔차의 길이가 비공유 값 만 가질 수 있습니다. 즉, 결과로 생성되는 r 벡터의 레이어 길이는 NA 값에 따라 달라집니다. 그러면 calc은 정의 된 수의 레이어와 스택의 길이가 다른 결과를 결합 할 수 없습니다.
기능을 사용하려면 길이가 r이고 속성 잔차는 비 NA 값으로 지정해야합니다.
사용자가 제공 한 데이터 집합에서 작동하는 다음 함수를 제안합니다. (1) 각 레이어 (완벽한 모델)에서 비교할 값이 두 개인 경우 모델 피팅을 피하십시오. (3) 어떤 이유로 든 모델이 적합하면 try을 추가하면 더 많은 테스트를 위해 쉽게 찾을 수있는 -1e32의 값을 출력합니다.

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) 

Residuals by layer

관련 문제