raster::calc 中的 tryCatch r

Question

我想编写一个函数，该函数将（希望）在raster 包中的光栅计算器中工作。我想要做的是将每个单元格值与时间向量进行回归。这有多个例子，但我想做的是尝试一种回归的方法（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() 移动到第二个回归函数？

Answer 1

这是另一种编码方式，也许会有所帮助：

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 失败的地方。

Answer 2

在阅读了更多内容并查看了@RobertH 的答案后，我编写了一些（非常）丑陋的代码来检查 GLS 是否会失败，如果失败，则改为执行 OLS。我很肯定有一种更好的方法可以做到这一点，但它可以工作并维护我的函数中定义的栅格图层名称，它还会将任何错误导出到 txt 文件。

for (i in 1) {
            j <- tempCentredRas
            cat(paste("Checking to see if gls(AR1) will work for model", m, r,"cell based calculations\n", sep = " "))
            ### This check is particularly annoying as it has to do this for every grid-cell
            ### it therefore has to perform GLS/OLS on every grid cell twice
            ### First to check if it (GLS) will fail, and then again if it does fail (use OLS) or doesn't (use GLS)
            possibleLocalError <- tryCatch(
              raster::calc(j, fun = funTrAR1),
              error = function(err) 
                err
            )
            if (inherits(possibleLocalError, "error")) {
              cat(paste("GLS regression failed for model", m, r, "using OLS instead for cell based results.","\n", sep = " "))
              cellResults <- raster::calc(j, fun = funTrOLS)
            } else {
              cellResults <- raster::calc(j, fun = funTrAR1)
            }
          }