3

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
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.

  • I think that there must be the capability. I was able to read in netcdf files with the names intact (gis.stackexchange.com/questions/122183/…), but am not able to export them. Maybe this is an issue with writeRaster. – user44796 Nov 17 '14 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.