利用NCL对没有时间维度信息的多个nc文件进行合并输出

大呆爱气象呀

      最近下载了国家气象信息中心的陆面再分析数据集CRA40_LAND逐月数据，从1980年1月至2022年2月，数据为nc格式文件，命名格式为CRA40LAND_LAND_YYYYMM_GLB_0P50_MONTH_V1_0_0.nc,为了对文件按照时间维度合并，一开始使用了CDO mergetime 命令，结果由于原始文件中并没有时间维度信息，合并失败，还是选择写NCL程序进行处理，并对各个变量添加时间维度信息。

• 首先查看文件中变量和维度信息
  

    cdo 命令 ：cdo sinfon **.nc

File format : NetCDF4
    -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter name
     1 : unknown  unknown  c instant       1   1     12177   1  F32  : Swnet_tavg   

        ......
    49 : unknown  unknown  c instant       1   1     12177   1  F32  : TotalPrecip_tavg
   Grid coordinates :
     1 : lonlat                   : points=12177 (123x99)
                              lon : 73.75 to 134.75 by 0.5 degrees_east
                              lat : 53.25 to 4.25 by -0.5 degrees_north
   Vertical coordinates :
     1 : surface                  : levels=1
     2 : generic                  : levels=4
     3 : generic                  : levels=4
     4 : generic                  : levels=4
cdo    sinfon: Processed 49 variables [0.34s 49MB].

netCDF命令: ncdump -h **.nc

netcdf CRA40LAND_LAND_198001_GLB_0P50_MONTH_V1_0_0 {
dimensions:
        lon = 123 ;
        lat = 99 ;
        SoilMoist_profiles = 4 ;
        SoilTemp_profiles = 4 ;
        SmLiqFrac_profiles = 4 ;
variables:
        float Swnet_tavg(lat, lon) ;
                Swnet_tavg:standard_name = "surface_net_downward_shortwave_flux" ;
                Swnet_tavg:long_name = "net downward shortwave radiation" ;
                Swnet_tavg:units = "W m-2" ;

    两种命令的结果都显示变量没有时间维度信息！！！！！
2.接下来编写一个NCL文件进行处理
begin
;;由于文件命名格式是统一的，所以可以进行批量读取
files=systemfunc("ls CRA40LAND*.nc")
;;计算文件数
n =dimsizes(files)
;;利用cdo命令读取文件中变量信息，存放在vars中
vars =systemfunc("cdo sinfon CRA40LAND*198001*.nc"print(vars);输出变量信息
;;结果如下：

Variable: vars
Type: string
Total Size: 488 bytes
            61 values
Number of Dimensions: 1
Dimensions and sizes:   [61]
Coordinates:
(0)        File format : NetCDF4
(1)         -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter name
(2)          1 : unknown  unknown  c instant       1   1     12177   1  F32  : Swnet_tavg
(3)          2 : unknown  unknown  c instant       1   1     12177   1  F32  : Lwnet_tavg
(4)          3 : unknown  unknown  c instant       1   1     12177   1  F32  : Qle_tavg
(5)          4 : unknown  unknown  c instant       1   1     12177   1  F32  : Qh_tavg

;;可以看到变量信息以字符格式输出，并且每一字符串有对应行号

;; 创建一个字符串数组，存放变量名称，数组大小为文件中比变量个数
VARNAME = new(49,"string")

;;接下来通过对字符串进行切片提取变量名称，在NCL官网看下str_split函数介绍，很简单

;;你的字符串有什么分隔符，在第二个参数位置写入就行，记得加" "号，对于我的例子，分隔符为空格和冒号，所以我填写的"  :"
do i = 1,49
       str = str_split(vars(i+1), " :") ; 对每一个字符串切片成字符串数组，这里的(i+1)索引是取决于print(vars)结果，跟每个字符串行号对应

       st = str(10)                          ; 切片后数组最后一个元素为变量名
       VARNAME(i-1) = st               ; 存放起来
end do

;;接下来就开始读取变量，这里进行批量读取，使用addfiles 不是addfile

f = addfiles(files, "r")

;; 将相同变量进行拼接，这一命令会在所有变量原有维度上最左侧加一维度，大小为文件数n，后续将这一维度名称改为时间
ListSetType (f, "join")

;;接下来创建新的nc文件用于存放拼接后文件
system("rm CRA40.nc")           ; 清除旧的文件
fout = addfile("CRA40.nc","c") ; 创建一个新的nc文件

;; 定义时间维度信息,这一步根据个人需求来，我是选择"months since 1980-01-01 00:00:0.0"

time = fspan(0,n-1,n)

time@units="months since 1980-01-01 00:00:0.0"time@calendar ="360 day"

;;这一步生成了n*6的时间信息矩阵

utc date = ut calendar(time, 0)

year= tointeger(utc date(:,0));第一维为年，第二维为月，类推

month = tointeger(utc date(:1))

day= tointeger(utc_date(:,2))

hour = tointeger(utc_ date(: 3))

minute = tointeger(utc_ date(:4))

second = utc date(: 5)

date = todouble(cd inv_calendar(year,month,day,hour,minute,second,time@units,0))

date!0 ="time"

date&time = date

date@units= time@units

date@calendar = time@calendar

;依次读取各个变量

do i= 0,48

     x= f[:]->$VARNAME(i)$;;这一步用到了我们之前存放好的变量名称数组

     x!0= "time" ;;修改第一维度名称为时间

     x&time = date ;;写入维度信息

     fout->$VARNAME(i)$=x    ;;写入新的nc文件

     delete(x)     ;;清楚临时变量，不要影响下一个循环

     print("i= " +i)

end do

    fout->time = date

end

3.运行结束后，查看输出的nc文件，时间信息已经被记录下来了。
     

Time coordinate :                             time : 506 steps
     RefTime =  1980-01-01 00:00:00  Units = months  Calendar = 360_day
  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
  1980-01-01 00:00:00  1980-02-01 00:00:00  1980-03-01 00:00:00  1980-04-01 00:00:00
  1980-05-01 00:00:00  1980-06-01 00:00:00  1980-07-01 00:00:00  1980-08-01 00:00:00
  1980-09-01 00:00:00  1980-10-01 00:00:00  1980-11-01 00:00:00  1980-12-01 00:00:00
  1981-01-01 00:00:00  1981-02-01 00:00:00  1981-03-01 00:00:00  1981-04-01 00:00:00

至此完成！