Why the mean value from NCO::ncra is defferent from terra::tapp()?
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)](https://nco.sourceforge.net/).
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][1]
[The plot of terra result][2]
[1]: https://i.stack.imgur.com/SIIZE.png
[2]: https://i.stack.imgur.com/5xkyz.png
