为何NCO::ncra与terra::tapp()计算的月均值存在差异？

我有一份包含2000年全年逐日气象数据的.nc文件，分别通过两种方式计算月均值：

1. 使用Linux下NetCDF Operator（NCO）的ncra命令，输出文件时间戳设为每月中旬，结果看似正常；
2. 使用R包terra的tapp()函数复核。

最终发现二者计算的均值存在明显差异，示例像素的数值对比如下：

完整代码实现 #

加载依赖包与数据准备 #

# load packages
require(pacman)
pacman::p_load(dplyr, tidyverse, terra, RColorBrewer, scales)

wd <- '~/RA/CMIP6/'

base <- 'https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/NEX-GDDP-CMIP6'
mdl <- 'ACCESS-CM2'
prd <- 'historical'
var <- 'tasmax'

# download data
urli <- paste0(base, '/', mdl, '/', prd, '/', 'r1i1p1f1', '/', var, '/', var, '_day_', mdl, '_', prd, '_', 'r1i1p1f1', '_gn_2000.nc')
setwd(wd)
system(paste0("axel -a -n 5 ", urli))
file.name <- paste0(var, '_day_', mdl, '_', prd, '_', 'r1i1p1f1', '_gn_2000.nc')

方式一：NCO ncra计算月均值 #

# ---------- ncra ----------------

# initial date for month
days_start <- seq(as.Date(paste0("2000-01-01")), length = 12, by = "months")
days_end <- seq(as.Date(paste0("2000-02-01")), length = 12, by = "months") - 1

# monthly mean
for(m in 1:12){
  system(paste0('ncra -O -d time,"', days_start[m], '","', days_end[m], '" ', wd, file.name, ' ', wd, 'monthly_', m, '.nc'))
}
# stack monthly mean in one .nc file
system(paste0('ncrcat ', wd, 'monthly_*.nc ', wd, gsub('day', 'month' , file.name)))

# delete intermediate files
system(paste0('rm ', wd, 'monthly_*.nc'))

r_nco <- rast(paste0(wd, gsub('day', 'month' , file.name)))

方式二：terra::tapp计算月均值 #

# ------------ terra::tapp -----------------

r <- rast(paste0(wd, file.name))
indices <- format(time(r), format = "%m") %>%
  as.double()
r_terra <- tapp(r, indices, fun = mean)
time(r_terra) <- seq(as.Date(paste0("2000-01-15")), length = 12, by = "months")

结果可视化与对比 #

# ---------- plot and compare --------------
col <- brewer_pal(palette = "RdYlBu", direction = -1)(5)
plot(r_nco, col = col)
plot(r_terra, col = col)

# find a pixel with value in 2000.01
time(r_nco)
r_nco[[1]][90490]
time(r_terra)
r_terra[[1]][90490]

可视化结果 #

• NCO计算结果图
• terra计算结果图

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