Why the mean value from NCO::ncra is defferent from terra::tapp()?

Question

I have a .nc file containing daily data for a whole year, and I tried to calculate the mean value for each month.

I tried to us Linux command ncra provided by netCDF Operator (NCO). The result seems nothing wrong, as the time stamp of the output file is the middle date of each month. However, when I use terra::tapp(daily_nc, indices2, fun = mean) to double check my result, it seems ncra get a different mean value

NCO::ncra	terra::tapp
251.8105	251.9165

Here is the complete code I used:

# load pacakges
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')

# ---------- 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 -----------------

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]

The plot of NCO result

The plot of terra result

Answer 1

NCO behaves almost as @Robert Wilson comments above. Unless told otherwise, NCO converts the on-disk type to double-precision, does the arithmetic, then converts back to the input/output type. This could explain the discrepancy if the terra::tapp method treats the data differently.