使用terra对SpatRaster像素时间序列应用SPEI::hargreaves的结果差异问题

我在R语言中基于terra包构建了两个带时间维度的SpatRaster对象a（对应Tmin）和b（对应Tmax），希望对每个像素的时间序列应用SPEI::hargreaves函数计算潜在蒸散。自行尝试实现时出现报错，后分别采用raster包的overlay方法和terra包的app方法完成计算：前者耗时超15分钟，后者仅需15秒，但两者输出结果的min、max值差异显著。因数据量超10TB我更倾向terra方案，疑惑为何结果不同，哪个是正确的？

构建带时间维度的SpatRaster #

library(SPEI)
library(terra)

a = array(1:(3*4*12*64),c(3,4,12*64))
a = rast(a)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
terra::time(a)=dates
names(a) <- zoo::as.yearmon(time(a))

b = array(1:(3*4*12*64),c(3,4,12*64))
b = rast(b)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
terra::time(b) <- dates
names(b) <- zoo::as.yearmon(time(b))

自行尝试的报错代码及信息 #

代码 #

library(SPEI)
library(terra)
library(zoo)

har <- function(a, b, lat) {
  SPEI::hargreaves(as.vector(a), as.vector(b), lat,na.rm = TRUE)
} 

lat <- init(rast(a), "y")
PET1 <- terra::lapp(x=c(a, b, lat), fun = Vectorize(har))

报错信息 #

Error: [lapp] cannot use 'fun'. The number of values returned is less 
than the number of input cells.
Perhaps the function is not properly vectorized

raster包实现方案 #

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

har <- function(Tmin, Tmax, lat) {
  SPEI::hargreaves(Tmin, Tmax, lat,na.rm = TRUE)
} 

lat <- init(raster(Tmin), "y")

PET_raster <- raster::overlay(Tmin, Tmax, lat, fun = Vectorize(har))

terra包实现方案（由Robert Hijmans提供） #

library(SPEI)
library(terra)
library(zoo)

lat <- init(rast(Tmin), "y")

r <- c(Tmin, Tmax, lat)
nl <- nlyr(Tmin)
nl2 <- nl + nl

PET_terra <- app(r, \(i) apply(i, 1, 
                       \(j) SPEI::hargreaves(j[1:nl], j[(nl+1):(nl2)], lat=j[nl2+1], verbose=FALSE)))

结果对比 #

raster包输出结果 #

PET_raster 
class      : RasterBrick
dimensions : 68, 104, 7072, 1020  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -98.285, -72.285, 39.875, 56.875  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : r_tmp_2025-01-29_175851.174327_307486_90781.grd
names      :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
min values :  0.0000000,  0.0000000,  0.0000000, 27.9045500, 40.7257577, 50.3362055, 41.2878256, 49.9874998, 47.6142315, 37.7349603,  7.4571957,  0.0000000,  0.0000000,  0.0000000,  0.0000000, ...
max values :   47.11585,   61.65605,   73.30580,  112.98168,  158.29348,  183.51450,  204.36278,  192.53549,  187.70001,  151.89846,   78.30323,   72.70451,   51.80651,   65.65896,   80.28027, ...

terra包输出结果 #

PET_terra 

class       : SpatRaster
dimensions  : 68, 104, 1020  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -98.285, -72.285, 39.875, 56.875  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
source(s)   : memory
names       :    lyr.1,    lyr.2,    lyr.3,    lyr.4,     lyr.5,    lyr.6, ...
min values  :  0.00000,  0.00000,  0.00000, 13.75159,  31.91098,  35.7943, ...
max values  : 17.54078, 28.36545, 49.92122, 89.34686, 138.97809, 171.4730, ...

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