NCL: Reading and writing NetCDF-4

NCL supports many NetCDF-4 features, including:

groups
compound data
"vlen" array (variable length arrays)
"string" type
"enum" data type
"opaque" data type

NetCDF-4 allows you to store attributes, variables, and dimensions in hierarchical groups, which look much like a UNIX file system. That is, you have a root group called "/", which then can have any number of groups under itr.

To determine if your NetCDF-4 file has any groups, you can use "ncdump" or "ncl_filedump":

    ncdump -h nc4_out_nc4uvt.nc

If you see something like:

    ...
    group: GROUPNAME {
      dimensions:
      ...
      variables:
      ...
  // group attributes:
                string :Conventions = "None" ;
                string :source_file = "nc4uvt.nc" ;
                string :title = "NCL generated netCDF file" ;
  } // group grp1

then this means the file has a group called GROUPNAME in it.

To use NCL to read any data under this group, first use the special => operator to access the group:

  f = addfile("nc4_out_nc4uvt.nc", "r")
  g = f=>GROUPNAME

The g variable is type "file", so you can then use all the regular "file" operators and functions on it to access data:

  print(getfilevarnames(g))
  V = g->V
  printVarSummary(v)

Below is a sample to read NetCDF file with group in it:

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

 setfileoption("nc", "FileStructure", "Advanced")

 fn = "nc4_out_nc4uvt.nc"

 fi = addfile(fn, "r")

 setfileoption("nc", "Format",  "NetCDF4Classic")

 printVarSummary(fi)
;print(fi)

 time = fi->time
 lev = fi->lev
 lat = fi->lat
 lon = fi->lon

 t = fi->T
 u = fi->U
 v = fi->V

;printVarSummary(t)
;printVarSummary(u)
;printVarSummary(v)

;print("t(0,0,0,0) = " + t(0,0,0,0))
;print("u(0,1,1,1) = " + u(0,1,1,1))
;print("v(0,2,2,2) = " + v(0,2,2,2))

;print("t&lat(0) = " + t&lat(0))
;print("t@units = " + t@units)
;print("u@units = " + u@units)
;print("v@units = " + v@units)

;u1 = fi->U(::2)
;v1 = fi->V(3:121:3)

;===================================================================

 g1 = fi=>/grp1

 print(g1)
 printVarSummary(g1)
;exit

 uf = fi->U

 wks  = gsn_open_wks("x11", "xy")

 res              = True
 res@tiMainString = "Basic plot" 

 plot1 = gsn_csm_contour(wks, uf(0,0,:,:), False)

 ug = g1->U

 plot2 = gsn_csm_contour(wks, ug(0,0,:,:), False)

Write NetCDF file with group:

 setfileoption("nc", "FileStructure", "Advanced")

 fn = "nc4uvt.nc"

 fi = addfile(fn, "r")

 setfileoption("nc", "Format",  "NetCDF4")

;printVarSummary(fi)
;print(fi)

 time = fi->time
 lev = fi->lev
 lat = fi->lat
 lon = fi->lon

 t = fi->T
 u = fi->U
 v = fi->V

;printVarSummary(t)
;printVarSummary(u)
;printVarSummary(v)

;print("t(0,0,0,0) = " + t(0,0,0,0))
;print("u(0,1,1,1) = " + u(0,1,1,1))
;print("v(0,2,2,2) = " + v(0,2,2,2))

;print("t&lat(0) = " + t&lat(0))
;print("t@units = " + t@units)
;print("u@units = " + u@units)
;print("v@units = " + v@units)

;u1 = fi->U(::2)
;v1 = fi->V(3:121:3)

 ntim  = dimsizes(time)	; get dimension sizes
 nlev  = dimsizes(lev)
 nlat  = dimsizes(lat)
 nlon  = dimsizes(lon)

;------------------------------------------------------------
;setfileoption("nc", "FileStructure", "Advanced")

 fon = "nc4_out_" + fn
 system("/bin/rm -f " + fon) ; remove if exists
 fo = addfile(fon, "c")

;===================================================================
; explicitly declare file definition mode. Improve efficiency.
;===================================================================
 setfileoption(fo,"DefineMode",True)

;setfileoption(fo,"CompressionLevel", 6)
;setfileoption(fo,"CacheSize", 3200000)
;setfileoption(fo,"CacheNelems", 1027)
;setfileoption(fo,"CachePreemption", 0.25)

; create global attributes of the file
;===================================================================
 fAtt               = True            ; assign file attributes
 fAtt@title         = "NCL generated netCDF file"
 fAtt@source_file   = fn
 fAtt@Conventions   = "None"
;fAtt@creation_date = systemfunc ("date")
 fileattdef(fo, fAtt)      ; copy file attributes

;===================================================================
; predefine the coordinate variables and their dimensionality
; Note: to get an UNLIMITED record dimension, we set the dimensionality
; to -1 (or the actual size) and set the dimension name to True.
;===================================================================
 dimNames = (/"time", "lev", "lat", "lon"/)
 dimSizes = (/ 1 ,   nlev,  nlat,  nlon /)
 dimUnlim = (/ True , False, False, False/)
 filedimdef(fo, dimNames, dimSizes, dimUnlim)

;===================================================================

 mtim = 1
 mlat = nlat/2
 mlon = nlon/2
 mlev = nlev/2
 chunkSizes = (/ mtim, mlev, mlat, mlon /)
;dimUnlim(0) = False
 filechunkdimdef(fo,dimNames,chunkSizes,dimUnlim)

;===================================================================

 grpnames = (/"grp1", "group2", "g3"/)

 filegrpdef(fo, grpnames)

; predefine the the dimensionality of the variables to be written out
;===================================================================
; Here we are using NCL functions to facilitate defining
; each variable's dimension name(s) and type.
;===================================================================
 filevardef(fo, "time", typeof(time), getvardims(time))
 filevarattdef(fo,"time", time)                   ; copy time attributes
 fo->time   = (/time/)

 filevardef(fo, "lev", typeof(lev), getvardims(lev) )
 filevarattdef(fo,"lev", lev)                     ; copy lev attributes
 fo->lev    = (/lev/)

 filevardef(fo, "lat", typeof(lat), getvardims(lat))
 filevarattdef(fo,"lat", lat)                     ; copy lat attributes
 fo->lat    = (/lat/)

 filevardef(fo, "lon", typeof(lon), getvardims(lon))
 filevarattdef(fo,"lon", lon)                     ; copy lon attributes
 fo->lon    = (/lon/)

 filevardef(fo, "T", typeof(t), getvardims(t))
 filevarattdef(fo,"T", t)                         ; copy T attributes
 filevarchunkdef(fo, "T", chunkSizes)
 filevarcompressleveldef(fo, "T", 2)
 fo->T      = (/t/)
;print(fo)
;printVarSummary(t)
;exit

 filevardef(fo, "U", typeof(u), getvardims(u))
 filevarattdef(fo,"U", u)                         ; copy U attributes
;filevarchunkdef(fo, "U", chunkSizes)
;filevarcompressleveldef(fo, "U", 4)
 fo->U      = (/u/)

 filevardef(fo, "V", typeof(v), getvardims(v))
 filevarattdef(fo,"V", v)                         ; copy V attributes
;filevarchunkdef(fo, "V", chunkSizes)
;filevarcompressleveldef(fo, "V", 6)
;filevarchunkcachedef(fo, "V", 3200000, 1027, 0.75)
 fo->V      = (/v/)

;===================================================================

 g1 = fo=>/grp1

 print(g1)

 fileattdef(g1, fAtt)
 filedimdef(g1, dimNames, dimSizes, dimUnlim)
 filechunkdimdef(g1,dimNames,chunkSizes,dimUnlim)

 filevardef(g1, "time", typeof(time), getvardims(time))
 filevarattdef(g1,"time", time)                   ; copy time attributes
 g1->time   = (/time/)

 filevardef(g1, "lev", typeof(lev), getvardims(lev) )
 filevarattdef(g1,"lev", lev)                     ; copy lev attributes
 g1->lev    = (/lev/)

 filevardef(g1, "lat", typeof(lat), getvardims(lat))
 filevarattdef(g1,"lat", lat)                     ; copy lat attributes
 g1->lat    = (/lat/)

 filevardef(g1, "lon", typeof(lon), getvardims(lon))
 filevarattdef(g1,"lon", lon)                     ; copy lon attributes
 g1->lon    = (/lon/)

 filevardef(g1, "T", typeof(t), getvardims(t))
 filevarattdef(g1,"T", t)                         ; copy T attributes
 filevarchunkdef(g1, "T", chunkSizes)
 filevarcompressleveldef(g1, "T", 2)
 g1->T      = (/t/)

 filevardef(g1, "U", typeof(u), getvardims(u))
 filevarattdef(g1,"U", u)                         ; copy U attributes
;filevarchunkdef(g1, "U", chunkSizes)
;filevarcompressleveldef(g1, "U", 4)
 g1->U      = (/u/)

 filevardef(g1, "V", typeof(v), getvardims(v))
 filevarattdef(g1,"V", v)                         ; copy V attributes
;filevarchunkdef(g1, "V", chunkSizes)
;filevarcompressleveldef(g1, "V", 6)
;filevarchunkcachedef(g1, "V", 3200000, 1027, 0.75)
 g1->V      = (/v/)

;===================================================================
 setfileoption(fo,"DefineMode",False)
;===================================================================
; output only the data values since the dimensionality and such have
; been predefined. The "(/", "/)" syntax tells NCL to only output the
; data values to the predefined locations on the file.
;====================================================================

 printVarSummary(g1)
;print(g1)
;exit

 printVarSummary(fo)
;print(fo)

;to = fo->T
;uo = fo->U
;vo = fo->V

;printVarSummary(to)
;printVarSummary(uo)
;printVarSummary(vo)

;print("to(0,0,0,0) = " + to(0,0,0,0))
;print("uo(0,1,1,1) = " + uo(0,1,1,1))
;print("vo(0,2,2,2) = " + vo(0,2,2,2))

;print("to&lat(0) = " + to&lat(0))
;print("to@units = " + to@units)
;print("uo@units = " + uo@units)
;print("vo@units = " + vo@units)

Read and write string

Write string:

 setfileoption("nc", "FileStructure", "Advanced")

 fn = "./tst_strings.nc"

 f = addfile(fn, "r")

;print(f)
;printVarSummary(f)

 strs = f->universal_declaration_of_human_rights

 print(strs)

We can write the above read string:

;===================================================================

 setfileoption("nc", "FileStructure", "Advanced")
 setfileoption("nc", "Format",  "NetCDF4")
 fn = "ncl_wrt_string.nc"
 system("/bin/rm -f " + fn)
 fo = addfile(fn, "c")

;setfileoption(fo,"DefineMode",True)

 fAtt               = True            ; assign file attributes
 fAtt@title         = "NCL generated netCDF 4 file"
 fAtt@source_file   = fn
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc("date")

;print(fAtt)

 fileattdef(fo, fAtt)

;===================================================================
 nstrs    = dimsizes(strs)
 dimNames = (/"Sentence"/)
 dimSizes = (/ nstrs    /)
 dimUnlim = (/ False    /)
 filedimdef(fo, dimNames, dimSizes, dimUnlim)

 strsdims = getvardims(strs)

 strsdims = "Sentence"

 filevardef(fo, "nc4_string", typeof(strs), strsdims)

 fo->nc4_string = (/strs/)

Read and write compound

If command:

ncdump -h ncl_wrt_comp_dat.nc

produces something like:

types:
  compound COMPOUNDNAME {
    ......
  };
dimensions:
        ......
variables:
        COMPOUNDNAME vn(...)
        ......

Then this file contains compound data. where COMPOUNDNAME is name of this compound data. Right below COMPOUNDNAME, there are few lines with type and names, where these names are called compound data components. To accesee these components, NCL uses "." to access those components, as demostrated below.

Read compound in NetCDF:

 setfileoption("nc", "FileStructure", "Advanced")

 fn = "comp_dat.nc"
;fn = "ncl_wrt_comp_dat.nc"

 f = addfile(fn, "r")

;print(f)
 printVarSummary(f)

 sf_id = f->/dimension_data.starfleet_id

;print(sf_id)
 printVarSummary(sf_id)

 abili = f->/dimension_data.abilities

;print(abili)
 printVarSummary(abili)

Write the above read compound in NetCDF:

;-----------------------------------------------------
 setfileoption("nc", "Format",  "NetCDF4")
 setfileoption("nc", "FileStructure", "Advanced")

 fon = "ncl_wrt_comp_dat.nc"
 system("/bin/rm -f " + fon)
 fo = addfile(fon, "c")

 fAtt               = True            ; assign file attributes
 fAtt@title         = "NCL generated netCDF file with compound"
 fAtt@source_file   = fon
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc("date")

;print(fAtt)

 fileattdef(fo, fAtt)

;===================================================================
 nvals   = 100
 dimNames = (/"starfleet"/)
 dimSizes = (/ nvals   /)
 dimUnlim = (/ False    /)
 filedimdef(fo, dimNames, dimSizes, dimUnlim)

 var_name = "dimension_data"
 mem_name = (/"starfleet_id", "abilities"/)
 mem_type = (/"integer", "integer"/)
 mem_size = (/1, 7/)
 filecompounddef(fo, "compound_var", var_name, dimNames, \
                 mem_name, mem_type, mem_size)

 complist = new(dimSizes, list)
 printVarSummary(complist)

 do n = 0, nvals - 1
    ListPush(complist[n], (/sf_id(n)/))
    s_abi = abili(n,:)
    ListPush(complist[n], (/s_abi/))

    delete(s_abi)
 end do

 filewritecompound(fo, "compound_var", var_name, mem_name, complist)

Read and write enum

read enum in NetCDF:

 setfileoption("nc", "FileStructure", "Advanced")

 fn = "../tst/tst_enums.nc"

 fi = addfile(fn, "r")

 print(fi)
;printVarSummary(fi)

 v = fi->primary_cloud

 enum_name = v@enum_name
 enum_value = v@enum_value

 print(v)

 print(enum_name)
 print(enum_value)

 nv = dimsizes(v)
 
 print(nv)

 tc = dimsizes(enum_value)

;do n = 0, nv - 1
;   m = ind(v(n) .eq. enum_value)
;   print("v{key, value} = {" + enum_value(m) + "," + enum_name(m) + "}")
;end do

write the above read data out NetCDF:

;===================================================================

 setfileoption("nc", "FileStructure", "Advanced")
 setfileoption("nc", "Format",  "NetCDF4")
 fn = "ncl_wrt_enum.nc"
 system("/bin/rm -f " + fn)
 fo = addfile(fn, "c")

 fAtt               = True            ; assign file attributes
 fAtt@title         = "NCL generated netCDF file with enum"
 fAtt@source_file   = fn
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc("date")

;print(fAtt)

 fileattdef(fo, fAtt)

;===================================================================
 nvals   = dimsizes(v)
 dimNames = (/"cloud_types"/)
 dimSizes = (/ nvals   /)
 dimUnlim = (/ False    /)
 filedimdef(fo, dimNames, dimSizes, dimUnlim)

 var_name = "cloud_observations"
 fileenumdef(fo, "enum_var", var_name, dimNames, enum_name, enum_value)

 fo->$var_name$ = (/v/)

Read and write opaque

read opaque data:

setfileoption("nc", "FileStructure", "Advanced")

 fn = "../tst/tst_opaques.nc"

 f = addfile(fn, "r")

 print(f)
;printVarSummary(f)

 v = f->var_defined_by_netcdf_user

 print(v)

write above read opaque data out:

;===================================================================

 setfileoption("nc", "FileStructure", "Advanced")
 setfileoption("nc", "Format",  "NetCDF4")
 fn = "ncl_wrt_opaque.nc"
 system("/bin/rm -f " + fn)
 fo = addfile(fn, "c")

 fAtt               = True            ; assign file attributes
 fAtt@title         = "NCL generated netCDF file with enum"
 fAtt@source_file   = fn
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc("date")

;print(fAtt)

 fileattdef(fo, fAtt)

;===================================================================
 nvals   = nv(0)
 dimNames = (/"opaque_base_size"/)
 dimSizes = (/ nvals   /)
 dimUnlim = (/ False    /)
 filedimdef(fo, dimNames, dimSizes, dimUnlim)

 var_name = "opaque_variable"
 fileopaquedef(fo, "opaque_type", var_name, var_size, dimNames)

 fo->$var_name$ = (/v/)