在 R 中创建/连接多维 NetCDF

Question

将此作为来源：How to concatenate monthly TRMM netCDF files into a single netCDF file using NCO or R on windows 7?

install.packages("ncdf4")
library(ncdf4)
install.packages("abind")
library(abind)
install.packages("RNetCDF")
library(RNetCDF)
install.packages("ncdf.tools")
library(ncdf.tools)
filenames=read.csv('TRMM.filenames.csv',head=F) 
filenames=as.character(filenames[,1]) 
n.lon=4 
n.lat=7 
NA.matrix=matrix(rep(NA,n.lon*n.lat),nrow=n.lon) 
prcp=array(NA.matrix,c(n.lon,n.lat,1)) 


for (i in 1:length(filenames)){ncdata=nc_open(filenames[i])
+ nc=ncvar_get(ncdata,"precipitation") prcp=abind(prcp,nc)}

prcp=prcp[,,-1] 
dim(prcp) 
saveRDS(prcp,'TRMM.all.rds')

我能够创建一个 rds 文件。但是，我真的很想将它保存为 nc 文件。我尝试通过以下方式创建一个具有 12 个步骤（每个月一个）时间维度的新 netCDF 文件：

dimx <- ncdim_def( "Lon", "degreesE", as.double(-90:-87))
dimy <- ncdim_def( "Lat", "degreesN", as.double(14:16))
dimTime <- ncdim_def( "Time", "months", 1:12, unlim=TRUE )
dimlist<-list(dimx,dimy,dimTime)
precip.ncvar<- ncvar_def("Precip", "mm/hr", dimlist, -1, longname="Precipitation", prec="float")
precip.nccreate<- nc_create( "precip.nccreate.nc", precip.ncvar, force_v4=FALSE, verbose=FALSE )
nc_close(precip.nccreate)

现在的挑战是将降水数据添加到每个月。
在第一个脚本之后，我尝试使用 ncvar_put 函数但没有成功。

filenames1=read.csv('TRMM.filenames.csv',head=F) 
filenames1=as.character(filenames1[,1]) 

for (i in 1:length(filenames1)){ncdata1=nc_open(filenames1[i])
nc1=ncvar_get(ncdata1,"precipitation") 
prcp1=abind(prcp1,nc1)}

n.lon1=4 
n.lat1=7

data2d<-(4*7)

for (i in 1:length(filenames1))
  ncvar_put( precip.nccreate, precip.ncvar, data2d, start=c(1), count=c(1) )

precip.nccreate<- nc_create( "precip.nccreate.nc", precip.ncvar, force_v4=FALSE, verbose=FALSE )

我得到的

ncvar_put 错误（precip.nccreate, precip.ncvar, data2d, start = c(1), : 找不到对象“precip.nccreate”

nc_create 中的错误（“precip.nccreate.nc”，precip.ncvar，force_v4 = FALSE，：找不到对象 'precip.ncvar'

无论如何，我想我只是想找到一种简单的方法将多个 netcdf 文件连接到一个 netcdf 中。

谢谢

Answer 1

由于时间维度是无限的，可以使用NCO的ncrcat命令，例如，

ncrcat in1.nc in2.nc ... out.nc
ncrcat in*.nc out.nc

Answer 2

TRMM 文件没有时间变量或维度。我无法获得 NCO 工具来处理这种情况，但我不是专家。我的解决方案是使用 python 创建文件的副本，同时根据输入文件名添加时间维度和时间值。然后按照@CharlieZender 的建议使用 ncrcat 。抱歉，脚本太长了。

#! /usr/bin/env python
import os
import glob
import netCDF4
import datetime

# get all the TRMM files with file names like, add a time dimension and time variable value from the file name.
# E.G. 3B43.19980101.7.HDF.nc, 3B43.19980201.7.HDF.nc, 3B43.20160701.7.HDF.nc
# Creates a list of monthly files which can be concatenated by NCO Tools ncrcat command.

# assumes all the 3B43.*.nc files are in the subdir files/
# creates out_3B43.*.nc in the current directory
infiles = glob.iglob( "files/*.nc" )

# add time dimension to these gridded variables. we only care about precipiation.
#wanted_vars = ['precipitation', 'relativeError', 'gaugeRelativeWeighting']
wanted_vars = ['precipitation']

for infile in sorted(infiles):
  try:
    nci = netCDF4.Dataset(infile, 'r')
  except RuntimeError:
    print "ERROR: Could not open " + infile
    exit()

  ofile = 'out_' + os.path.basename(infile)

  try:
    nco = netCDF4.Dataset(ofile,'w',clobber=True, format=nci.data_model)
  except RuntimeError:
    print "ERROR: Could not open " + ofile
    exit()

  # copy global attributes
  in_g_attdict = nci.__dict__
  # copy dimensions
  nco.setncatts(in_g_attdict) 
  for dname, d in nci.dimensions.iteritems():
    if d.isunlimited():
      d_size = None
    else:
      d_size = len(d)
    nco.createDimension(dname, d_size)
  # Add a time dimension, as unlimited, None means unlimited.
  nco.createDimension('time', None)

  # Get YYYYMMDD from file name you could use dy = 15 to reflect that these are monthly files.
  file_date = infile.split('.')[1]
  yr = file_date[0:4]
  mo = file_date[4:6]
  dy = file_date[6:]
  odt = datetime.datetime(year=int(yr), month=int(mo), day=int(dy))

  # create the time variable. Note: the time variable must go first to be the outer or record dimension for ncrcat
  time_v = nco.createVariable('time', 'int32', ('time'))
  time_v.setncattr('standard_name', 'time') 
  time_v.setncattr('long_name', 'reference time of data field') 
  time_v.setncattr('axis', 'T') 
  # for time_urnits, I've used seconds since the epoch. But I'm guessing that 'days since XXX' should work as well. But harded to
  time_v.setncattr('units', 'seconds since 1970-01-01 00:00:00 UTC')
  # calculate the number.
  ntime = netCDF4.date2num(odt, time_v.units)
  time_v[:] = ntime

  # Copy variables, skip the wanted_vars, for now
  for iname, ivar in nci.variables.iteritems():
    if iname in wanted_vars:
      continue
    ovar = nco.createVariable(iname, ivar.dtype, ivar.dimensions)
    iattdict = ivar.__dict__
    ovar.setncatts(iattdict) 
    ovar[:] = ivar[:]
  # now copy the wanted 2-D gridded variables, adding a time dimension, as axis 0
  for data_var_name in wanted_vars:
    ivar = nci.variables[data_var_name]
    # original precipation variable is dimensioned by nlon,nlat
    ovar = nco.createVariable(data_var_name, ivar.dtype, ('time', 'nlon', 'nlat') )
    iattdict = ivar.__dict__
    ovar.setncatts(iattdict) 
    ovar[0,:,:] = ivar[:]

  nci.close()
  nco.close()

exit()

这将在当前目录中留下 out_xxxxx.nc 文件，并且可以使用 ncrcat。

ncrcat -H -h out_3B43.19980101.7.HDF.nc out_3B43.19980201.7.HDF.nc final.nc

Answer 3

时间上的合并是通过这种方式与 CDO 完成的：

cdo mergetime in1.nc in2.nc ... in12.nc out.nc

我经常使用 cdo 来操作 TRMM 文件。