r - 在栅格::calc中的tryCatch r

我想编写一个函数(希望)在raster包中的栅格计算器中工作。我想做的是将每个像元值相对于时间 vector 进行回归。有很多这样的例子，但是我想做的是让该方法尝试一种类型的回归(gls，控制AR1残差)，但是如果由于某种原因回归抛出错误(也许没有AR1)然后将其恢复为简单的OLS回归。

我为回归编写了两个函数。一个gls:

# function for calculating the trend, variability, SNR, and residuals for each pixel
## this function will control for AR1 structure in the residuals
funTrAR1 <- function(x, ...) {if (sum(is.na(x)) >= 1) { NA } else {
  mod <- nlme::gls(x ~ Year, na = na.omit, method = "REML", verbose = TRUE,
                   correlation = corAR1(form = ~ Year, fixed = FALSE),
                   control = glsControl(tolerance = 1e-3, msTol = 1e-3, opt = c("nlminb", "optim"),
                                        singular.ok = TRUE, maxIter = 1000, msMaxIter = 1000))
  slope <- mod$coefficients[2]
  names(slope) <- "Trend"
  var <- sd(mod$residuals)
  names(var) <- "Variability"
  snr <- slope/var
  names(snr) <- "SNR"
  residuals <- c(stats::quantile(
    mod$residuals, probs = seq(0,1,0.25),
    na.rm = TRUE, names = TRUE, type = 8),
    base::mean(mod$residuals, na.rm = TRUE))
  names(residuals) <- c("P0", "P25", "P50", "P75", "P100", "AvgResid")
  return(c(slope, var, snr, residuals))}
}

和OLS:

# function for calculating the trend, variability, SNR, and residuals for each pixel
## this function performs simple OLS
funTrOLS <- function(x, ...) {if (sum(is.na(x)) >= 1) { NA } else {
  mod <- lm(x ~ Year, na.action = na.omit)
  slope <- mod$coefficients[2]
  names(slope) <- "TrendOLS"
  var <- sd(mod$residuals)
  names(var) <- "VariabilityOLS"
  snr <- slope/var
  names(snr) <- "SNROLS"
  residuals <- c(stats::quantile(
    mod$residuals, probs = seq(0,1,0.25),
    na.rm = TRUE, names = TRUE, type = 8),
    base::mean(mod$residuals, na.rm = TRUE))
  names(residuals) <- c("P0", "P25", "P50", "P75", "P100", "AvgResid")
  return(c(slope, var, snr, residuals))}
}

我试图将它们包装在一个tryCatch表达式中，该表达式可以传递给raster::calc

xReg <- tryCatch(
  {
    funTrAR1
  },
  error = function(e) {
    ## this should create a text file if a model throws an error
    sink(paste0(inDir, "/Outputs/localOLSErrors.txt"), append = TRUE)
    cat(paste0("Used OLS regression (grid-cell) for model: ", m, ". Scenario: ", t, ". Variable: ", v, ". Realisation/Ensemble: ", r, ". \n"))
    sink()
    ## run the second regression function
    funTrOLS
  }
)

然后将此函数传递给raster::calc，如下所示

cellResults <- calc(rasterStack, fun = xReg)

如果一切正常，将产生一个类似于此输出变量的栅格堆栈



但是，对于我的某些数据集，我正在停止运行所有循环的循环会收到以下错误:

 Error in nlme::gls(x ~ Year, na = na.omit, method = "REML", verbose = TRUE,  :
  false convergence (8)

这直接来自nlme::gls和我希望避免的事情。我以前从未使用过tryCatch(这可能非常明显)，但是如果第一次(AR1)回归失败，有谁知道如何使tryCatch()移至第二个回归函数？

最佳答案

这是另一种编码方式，可能会有所帮助:

xReg <- function(x, ...) {
    r <- try(funTrAR1(x, ...), silent=TRUE)
    # if (class(r) == 'try-error') {
    if (!is.numeric(r)) {  # perhaps a faster test than the one above
        r <- c(funTrOLS(x, ...), 2)
    } else {
        r <- c(r, 1)
    }
    r
}

我添加了一个图层，该图层显示每个单元使用的模型。

你也可以

xReg <- function(x, ...) {
    r <- funTrOLS(x, ...)
    try( r <- funTrAR1(x, ...), silent=TRUE)
    r
}

或两次使用calc，然后再使用cover

xReg1 <- function(x, ...) {
    r <- c(NA, NA, NA, NA)
    try( r <- funTrAR1(x, ...), silent=TRUE)
    r
}
xReg2 <- function(x, ...) {
    funTrOLS(x, ...)
}

a <- calc(rasterStack, xReg1)
b <- calc(rasterStack, xReg2)
d <- cover(a, b)

a将向您显示xReg1失败的位置。

关于r - 在栅格::calc中的tryCatch r，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/49287742/