R语言两时间序列栅格逐像元线性回归calc函数报错排查

分析背景 #

• 待处理数据：两套空间分辨率、时间范围完全一致的年度格网栅格数据集（NPP、Tmax），均为RasterBrick对象，行列数321×401，共128721个像元，各包含37个年度图层，分辨率0.0375°×0.0375°，空间范围为xmin=32.98125、xmax=48.01875、ymin=2.98125、ymax=15.01875，坐标系为WGS84（+proj=longlat +datum=WGS84 +no_defs），源文件分别为Tmax_ANNUAL_RES01.tif、NPP_ANNUAL_RES01.tif。
• 分析目标：开展逐像元简单线性回归分析，以Tmax为自变量预测NPP，提取每个像元对应的回归斜率。

已运行代码 #

首先将两个栅格文件合并为74层的RasterStack对象：

stack_ANNUAL <- stack(list.files(pattern="*.tif", full.names=TRUE))

随后编写自定义函数逐像元拟合线性模型提取斜率：若像元第一个值为NA则返回NA，否则以前37层（NPP年度序列）为因变量、后37层（Tmax年度序列）为自变量拟合线性模型，提取斜率系数，通过calc函数生成全区域斜率栅格：

fun=function(x) { if (is.na(x[1])){ NA } else { lm(x[1:37] ~ x[38:74])$coefficients[2] }}
slope <- calc(stack_ANNUAL, fun)

触发报错信息 #

Error in (function (classes, fdef, table)  : 
unable to find an inherited method for function ‘writeValues’ for signature 
‘"RasterBrick", "numeric"’

报错原因 #

calc函数要求自定义函数对所有像元返回的结果类型、长度完全统一。原代码存在两个核心问题：

1. NA判断分支返回的是R默认的逻辑型NA（数据类型为logical），而正常拟合模型返回的斜率是长度为1的数值型（numeric）结果，两类结果混杂导致输出数据类型不统一，raster包无法将混杂类型的结果写入RasterBrick对象，最终触发writeValues方法匹配失败的报错。
2. 未对序列其余位置的缺失值做判断，若像元时间序列中间存在NA值，lm拟合会自动剔除缺失值，可能出现返回结果长度异常的问题。

修复方案 #

将分支返回的缺失值替换为数值型NA_real_，保证所有分支返回值类型统一为长度1的数值向量，同时增加全序列缺失值判断，避免拟合异常：

fun <- function(x) {
  # 首个值为NA直接返回数值型缺失值
  if (is.na(x[1])) return(NA_real_)
  # 分别提取NPP、Tmax时间序列
  npp_ts <- x[1:37]
  tmax_ts <- x[38:74]
  # 序列任意位置存在缺失值直接返回NA
  if (any(is.na(npp_ts)) | any(is.na(tmax_ts))) return(NA_real_)
  # 拟合线性模型返回斜率
  lm(npp_ts ~ tmax_ts)$coefficients[2]
}
slope <- calc(stack_ANNUAL, fun)

注意：list.files默认按文件名字母顺序读取文件，运行前需确认合并后的RasterStack图层顺序，保证前37层为NPP、后37层为Tmax，避免自变量和因变量顺序错配导致结果错误。

内容的提问来源于stack exchange，提问作者Addisu D