栅格数据集Thornthwaite蒸散计算报错：公式未向量化求助

我帮你梳理下代码里的问题，以及对应的修复方案：

首先看你遇到的错误：

Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : cannot use this formula, probably because it is not vectorized

这个错误主要有几个原因，咱们一步步来修正：

1. 函数定义的语法错误 #

你写的函数定义thornthwaite ET th <- function(Tave, lat) {有语法问题——R里的函数名不能包含空格，得改成合法的名称，比如th_fun。

2. 纬度栅格未正确创建 #

你代码里的lat <- setValues(a, coordinates(tm)[, "y"])中，变量a根本没定义，得先创建一个和温度栅格tm同分辨率的空栅格，再填充纬度值。

3. 确保函数支持向量化输入 #

raster::overlay要求传入的函数能处理向量形式的输入：当你把多层时间序列栅格tm和单层纬度栅格传给overlay时，它会把每个像素的所有时间步温度值作为一个向量，对应像素的纬度作为标量传给你的函数，函数需要返回该像素每个时间步的ET值（也是向量）。

修复后的完整代码 #

library(SPEI)
library(raster)
library(zoo)

# 创建模拟温度栅格数据
tm <- array(1:(3*4*12*64), c(3,4,12*64))
tm <- brick(tm)

# 设置时间序列
dates <- seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
tm <- setZ(tm, dates)
names(tm) <- as.yearmon(getZ(tm))

# 正确创建与tm同分辨率的纬度栅格
lat_raster <- raster(tm)
lat_raster <- setValues(lat_raster, coordinates(tm)[, "y"])

# 定义支持向量化的Thornthwaite计算函数
th_fun <- function(tave_timeseries, lat) {
  # tave_timeseries：单个像素的所有时间步温度值（向量）
  # lat：该像素的纬度（标量）
  SPEI::thornthwaite(tave_timeseries, lat)
}

# 计算每个像素的ET时间序列
out <- raster::overlay(tm, lat_raster, fun = th_fun)

补充说明 #

SPEI::thornthwaite本身就支持处理时间序列向量，所以只要函数能正确接收每个像素的温度序列和纬度值，就能顺利输出对应时间步的ET栅格。运行上述代码后，out就是每个时间步的ET栅格数据集，可以直接用于后续的SPEI计算。

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