R: How to write CRS information correctly to NetCDF?

Question

As the title says, I'm struggling with writing the CRS of a multi layer raster stack to a NetCDF file using a custom package. This package uses the functions of the ncdf4-package. But the written .nc file doesn't have the CRS information which the source raster had. I've already read this thread and I tried to replicate it, tweaking the WriteNCDF function from above. I tried it with:

ncatt_put(nc, var.name, 'proj4', crs(data), 'double')

with ncas the to be written ncdf file and data as the rasterlayer. This results in the following error:

Error in ncatt_put_inner(idobj$group_id, idobj$id, attname, attval, prec = prec, : ncatt_put: error, passed an attribute with a storage mode not handled. Att name: proj4 Att value: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs Storage mode passed: S4 . Handled types: integer double character

Afterwards I tried it with:

ncatt_put(nc, var.name, 'proj4', crs(data, asText = T), 'text')

This doesn't produce an error but when I open the written .nc file with raster::brick(), the CRS isn't found. So how do I put CRS information correctly into a NetCDF file?

Minimal code to run:

library('ncdf4')
library('raster')
source('WriteNCDF.R')


r <- raster(crs="+proj=utm +zone=32 +datum=WGS84", xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- 1:ncell(r)
s <- stack(r, r/2)
x <- writeRaster(s, filename="test.nc")
crs(x)

z <- WriteNCDF(data.l = list(s$layer.1, s$layer.2), 
               var.name = c('layer1', 'layer2'), 
               var.unit = c('unit', 'unit'),
               file = 'test2.nc',
               utm = T)
zz <- brick('test2.nc', varname = 'layer1')
crs(zz)

The tweaked WriteNCDF.R can be found here.

Answer 1

The test.nc file written via raster::writeRaster uses the CF-1.4 metadata convention for spatial data. If you look at it with ncdump you can see that the CRS is stored as a proj4 attribute of an integer variable called crs:

variables:
    int crs ;
        crs:proj4 = "+proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" ;

The crs variable itself seems to have a missing value - but that's okay, the text is in the attribute above:

data:

 crs = _ ;

 easting = 0.0138888888888889, 0.0416666666666667, 0.0694444444444444, 
    0.0972222222222222, 0.125, 0.152777777777778, 0.180555555555556, 

The WriteNCDF function stores the CRS as an attribute of the layer data. ncdump shows this:

    float layer1(time, north, east) ;
            layer1:units = "unit" ;
            layer1:_FillValue = -9999.f ;
            layer1:scale_factor = 1. ;
            layer1:add_offset = 0. ;
            layer1:proj4 = "+proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" ;

Now I don't think that conforms to the CF-1.4 metadata spec so I doubt any geospatial software is going to read that back when you read the file.

You can either make it conform to the CF-1.4 metadata spec by writing the correct attributes etc (but that might not be possible if your data fundamentally won't fit the spec), or write a counterpart ReadNCDF function that uses the ncdf4 package functions to read these attributes directly after reading the raster grid with raster - for example:

> s2 = nc_open("test2.nc")
> ncatt_get(s2, "layer1","proj4")$value
[1] "+proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

which opens the NetCDF file then gets the "proj4" attribute from variable "layer1".