Export band names with NetCDF file in R

Question

I am using R to work with multiband raster files and am running into a problem. I have named the bands, but when I export in netcdf format, the band names are being stripped. Anyone know how to specify band names when exporting netcdf files?

For example, if I have a raster object named "ten", and I type names(ten) I get the vector of band names. If I write ten to disk and then reimport using

newraster = brick("ten.nc"), then the names are stripped away. However, brick() is able to import other .nc files fine with the names attached.

The code I am using is

z = writeRaster(ten,filename = a,overwrite=TRUE,) ### Output to ./outdir folder. This is monthly data

Answer 1

It's because NetCDF does not have names for each slice in the 3rd (and higher dims), but raster does. NetCDF has a name for a "variable" (which is the array), but raster has a name for every slice in the variable. (This is the standard mess where we conflate data fields/attributes with dimensions). There's no straightforward way to store these names in NetCDF without messing around - you'd have to create a dim-var to store them essentially, but there's no clear model of using such a thing in NetCDF.

See how the variable gets the name, but when we round-trip the individual "layers" simply get a layer name, and the varname applies to the whole array. It's the same in GDAL, it unwraps the array to multiple bands (layers) - NetCDF has nowhere for the "dsf1" and "dsf2" to go without extra work at both ends (these would be equivalent to R's dimnames for an array).

library(raster)
r0 <- raster(volcano)
r <- brick(r0, r0 * 2)
names(r) <- paste0("dsf", seq_len(nlayers(r)))
writeRaster(r, "a.nc")

class       : RasterBrick 
dimensions  : 87, 61, 5307, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.01639344, 0.01149425  (x, y)
extent      : 0, 1, 4.943962e-17, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /rdsi/PRIVATE/home/mdsumner/Git/nectar/a.nc 
names       : X1, X2 
unknown     : 1, 2 
varname     : variable 

system("gdalinfo a.nc")
Driver: netCDF/Network Common Data Format
Files: a.nc
Size is 61, 87
Coordinate System is `'
Origin = (0.000000000000000,1.000000000000000)
Pixel Size = (0.016393442622951,-0.011494252873563)
Metadata:
  easting#long_name=easting
  easting#units=meter
  NC_GLOBAL#Conventions=CF-1.4
  NC_GLOBAL#created_by=R, packages ncdf and raster (version 2.3-6)
  NC_GLOBAL#date=2014-11-16 07:15:24
  NETCDF_DIM_EXTRA={value}
  NETCDF_DIM_value_DEF={2,4}
  NETCDF_DIM_value_VALUES={1,2}
  northing#long_name=northing
  northing#units=meter
  value#long_name=value
  value#units=unknown
  variable#_FillValue=-3.4e+38
  variable#long_name=variable
  variable#max={195,390}
  variable#min={94,188}
  variable#missing_value=-3.4e+38
Corner Coordinates:
Upper Left  (   0.0000000,   1.0000000) 
Lower Left  (   0.0000000,   0.0000000) 
Upper Right (   1.0000000,   1.0000000) 
Lower Right (   1.0000000,   0.0000000) 
Center      (   0.5000000,   0.5000000) 
Band 1 Block=61x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-3.39999999999999996e+38
  Metadata:
    _FillValue=-3.4e+38
    long_name=variable
    max={195,390}
    min={94,188}
    missing_value=-3.4e+38
    NETCDF_DIM_value=1
    NETCDF_VARNAME=variable
Band 2 Block=61x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-3.39999999999999996e+38
  Metadata:
    _FillValue=-3.4e+38
    long_name=variable
    max={195,390}
    min={94,188}
    missing_value=-3.4e+38
    NETCDF_DIM_value=2
    NETCDF_VARNAME=variable

brick("a.nc")
class       : RasterBrick 
dimensions  : 87, 61, 5307, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.01639344, 0.01149425  (x, y)
extent      : 0, 1, 4.943962e-17, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /rdsi/PRIVATE/home/mdsumner/Git/nectar/a.nc 
names       : X1, X2 
unknown     : 1, 2 
varname     : variable 

And with just one layer, it's clear where the name "dsf" applies - to the entire array:

library(raster)
r <- raster(volcano)
names(r) <- "dsf"
writeRaster(r, "a.nc")

class       : RasterLayer
dimensions  : 87, 61, 5307  (nrow, ncol, ncell)
resolution  : 0.01639344, 0.01149425  (x, y)
extent      : 0, 1, 4.943962e-17, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : a.nc
names       : dsf
zvar        : dsf


system("ncdump -h a.nc")

netcdf a {
dimensions:
    easting = 61 ;
    northing = 87 ;
variables:
    double easting(easting) ;
            easting:units = "meter" ;
            easting:long_name = "easting" ;
    double northing(northing) ;
            northing:units = "meter" ;
            northing:long_name = "northing" ;
    float dsf(northing, easting) ;
            dsf:_FillValue = -3.4e+38 ;
            dsf:missing_value = -3.4e+38 ;
            dsf:long_name = "dsf" ;
            dsf:min = 94. ;
            dsf:max = 195. ;

Without being too critical (since raster is awesome), I think names(raster) is mis-applied, if the raster model was more general it would give only one name no matter whether the grid was 3d, 4d or more, but if you had a different kind of variable - say "temperature" and "pressure" then those would be the names, independent of the dimension of the data. This is a major source of confusion in the murky waters between GIS and the other worlds of modelling.