可以使用以下代码

library(terra)

#List of .nc files in directory
fls <- list.files("/folder/with/files", ".nc$", full.names = TRUE)

#list to SpatRaster object conversion
r_list <- rast(fls)

#For global mean calculation i.e. single mean for the entire SpatRaster
global(r_list, mean, na.rm=TRUE)

#For local mean calculation i.e. pixel-wise mean calculation
mean(r_list, na.rm=TRUE)

terra::mosaic()用于在空间上合并一组光栅。在重叠的单元格中，您可以应用一个函数来确定要将单元格设置为的值，默认值为mean。所以，是的，你的方法会起作用，但效率非常低。
在R中，你可以这样做：

library(abind)

# Get a list of all files to be merged
lf <- list.files("/folder/with/files", "^AQUA.*nc$", full.names = T)

# Read in all the data
d <- apply(lf, function(f) {n <- nc_open(f); x <- ncvar_get(n, "sst"); nc_close(n); x})

# Turn the list of data matrices into an array
arr <- abind(d, along = 3)

# Calculate the mean over the 3rd dimension of the array to create a matrix
sst <- rowMeans(arr, na.rm = T, dims = 2)

如果需要，可以将sst矩阵转换为terraSpatRaster。
这段代码可以满足你的所有需求，非常高效。

内存限制

如果您必须以这种方式合并许多文件，或者使用不同的数据集，其维度更大，则可能会遇到内存问题，这将极大地减慢该过程。如果是这种情况，请使用循环。循环开销与其余代码（主要是阅读数据）相比将完全相形见绌。

# Get a list of all files to be merged
lf <- list.files("/folder/with/files", "^AQUA.*nc$", full.names = T)

# Function to calculate mean SST from multiple MODIS AQUA files
sst <- function(lf) {
  for (ndx in 1:length(lf)) {
    n <- nc_open(lf[[ndx]])
    if (ndx == 1) {
      data <- ncvar_get(n, "sst")
    } else {
      data <- data + ncvar_get(n, "sst")
    }
    nc_close(n)
  }
  data / length(lf)
}

这里的技巧是，没有比两个矩阵更大的内存需求了：data对象和ncvar_get()读入的数据。