writeRaster does not preserve names when writing to NetCDF

Hi Barry, Michael,

Thanks for the replies. Maybe I should clarify my point a bit more.  As
Barry says, NetCDF is well suited for storing high-dimensional data. For
example, at the bottom is the header description from a netcdf file written
with writeRaster(). Note that the variable is stored as a 3-D array with
dimensions "value", "easting" and "northing". The values of the eastings and
northings are written appropriately, but the value dimension/variable is
just filled with a numeric value. This dimensions is used by brick() when
reading a netcdf to populate the layer names.

My proposed solution would be to add an argument to writeRaster() called
zvalues (or similar), which would be written to the value
dimension/variable. Initially this could be a numeric vector, (e.g. the
depths that this variable corresponds to: 6, 12, 18, 30, 40, 50...) but I
believe it may be possible to make this a string as well e.g. with NetCDF
v4.

Does that make sense?

Mark



netcdf EP_1992 {
dimensions:
        easting = 23 ;
        northing = 37 ;
        value = UNLIMITED ; // (6 currently)
variables:
        double easting(easting) ;
                easting:units = "meter" ;
                easting:long_name = "easting" ;
        double northing(northing) ;
                northing:units = "meter" ;
                northing:long_name = "northing" ;
        int value(value) ;
                value:units = "unknown" ;
                value:long_name = "value" ;
        float variable(value, northing, easting) ;
                variable:_FillValue = -3.4e+38 ;
                variable:missing_value = -3.4e+38 ;
                variable:long_name = "variable" ;
                variable:projection = "+proj=utm +zone=29 +datum=WGS84
+units=km +ellps=WGS84 +towgs84=0,0,0" ;
                variable:projection_format = "PROJ.4" ;
                variable:min = -2.26894027129976, -4.60561400334956,
-3.66937747634895, -5.97212100926874, -6.12123420979669, -3.44363500443944 ;
                variable:max = 3.11489873186157, 3.39906242513969,
3.3160419018914, 2.74477652780673, 1.8108377756661, 0.70816760511994 ;


// global attributes:
                :Conventions = "CF-1.4" ;
                :created_by = "R, packages ncdf and raster (version 2.2-31)"
;
                :date = "2014-08-11 15:37:25" ;
data:

 easting = 25, 75, 125, 175, 225, 275, 325, 375, 425, 475, 525, 575, 625,
    675, 725, 775, 825, 875, 925, 975, 1025, 1075, 1125 ;

 northing = 6425, 6375, 6325, 6275, 6225, 6175, 6125, 6075, 6025, 5975,
5925,
    5875, 5825, 5775, 5725, 5675, 5625, 5575, 5525, 5475, 5425, 5375, 5325,
    5275, 5225, 5175, 5125, 5075, 5025, 4975, 4925, 4875, 4825, 4775, 4725,
    4675, 4625 ;

 value = 1, 2, 3, 4, 5, 6 ;

}



On 12/08/14 09:53, Barry Rowlingson wrote:

On Tue, Aug 12, 2014 at 12:05 AM, Michael Sumner <mdsumner at gmail.com>
wrote:

Hi,

I don't think that NetCDF has the capacity to store these names in a
natural way. Internally your brick is stored as a NetCDF variable, in
this case a 3d array. To do this it would need to define a 'dimension"
and include a lookup map to the names. A NetCDF dimension variable
must be populated with numbers I believe, so you would need another
variable to store those and then add in features to raster to use
them.

There's at least two blockers here:
- the NetCDF model is not designed for storing tabular data (a RGB
image is a table of records in one sense, not a 3d array of integers)
- gridded data storage usually conflates "dimension" with
"attributes", which is unfortunate.

  Ummmm....

  NetCDF can store just about anything. Its that flexibility that makes
it hard for one piece of software that writes a raster brick totally
compatible with another piece of software that writes a raster brick.
The NetCDF format itself has no convention for how to store it.

  So people have developed conventions!

http://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html
http://www.unidata.ucar.edu/software/netcdf/conventions.html

http://cfconventions.org/Data/cf-convetions/cf-conventions-1.7/build/cf-conventions-multi.html

[yes there's a typo in that last one - "convetions" - that's how it is.]

  You can see which convention a NetCDF file claims to be by dumping it
with ncdump on the command line:

  // global attributes:
         :Conventions = "CF-1.4" ;
         :created_by = "R, packages ncdf and raster (version 2.1-49)" ;
         :date = "2014-08-11 17:17:01" ;

[note this also dumps the data, you might want to pipe it to | more]

The exact convention used by raster and rgdal seems to vary depending
on whether ncdf is loaded when raster saves it, or something. At least
in my blurry morning fumblings I seem to have reached contradictory
examples which need clearing up once the coffee IV hits. The varname
argument to writeRaster doesn't seem to do anything, but I'm probably
looking at the wrong file or have an old package:raster or something.

Anyway, the two main ways of storing multi-layer raster appear to be,
and these are files with CF-1.4/5 convention markers, as a 3d cube
with dimension (Nx,Ny,Nz) or as a set of 2d named variables. The names
of these 2d variables should be settable at write time and recoverable
at read-time. However, If the data is a 3d cube then I don't think the
dimensions can have names, although I've not really understood the CF
standards yet (and that page has links to some awful documentation and
missing links...).

It would be quite simple to add another variable, of characters, to
the NetCDF file, but I suspect keeping to a given NetCDF Convention
would be a good idea assuming it can do it.

End-of-ramble...

Barry

On Tue, Aug 12, 2014 at 12:26 AM, Mark Payne <markpayneatwork at gmail.com>
wrote:

Hi,

I have a problem with writeRaster when writing a brick to NetCDF - the
names of the layers are not preserved. Here is a minimum example to
demonstrate the point:

#Create an arbitrary brick and write it out
b <- brick(system.file("external/rlogo.grd", package="raster"))
fname <- file.path( tempdir(),"test.nc")
dmp <- writeRaster(b,file=fname)

#Now read it back in
a <- brick(fname)

#However, names are not preserved

names(b)

[1] "red"   "green" "blue"

names(a)

[1] "X1" "X2" "X3"

I realise this is a touch tricky, as NetCDF isn't exactly made for
storing
character strings, which is how raster stores the layer names. However,
there are many cases where the string is just a representation of a
number
anyway.... e.g. when working with time or depth levels as the layers in
the
brick. One solution could therefore be to allow the writeRaster function
to
take an argument "zvals", in the same way that it takes "zname" and
"zunit"
when working with netcdf files - this argument would let you specify the
numeric values of the z-axis. Everything else would then work ok I
think,
as brick() picks up the values of the z-axis, converts them to strings
and
assigns them to layer names.

Cheers,

Mark

raster version: 2.2-31
R version 3.0.3
ncdf4 version: 1.10
Platform: Linux Mint 16, 64 bit

         [[alternative HTML version deleted]]

--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner at gmail.com