2

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.

3
  • That last link goes to a 54.6Mb NetCDF file.... Did you mean to point to some R code? – Spacedman Dec 15 '20 at 22:07
  • Why is using writeRaster(s, filename="test.nc") not doing what you want? It saves the data and the CRS. What is WriteNCDF supposed to do extra? – Spacedman Dec 15 '20 at 22:09
  • @Spacedman I fixed the link. The problem I have with writeRaster() is on one hand the limitation of meta-data you can write to the .nc files and on the other hand, how to write a multi-layer raster stack where each layer has the same time dimension. If you write a multi-layer raster stack (e.g. all Sentinel-2 bands) with writeRaster() ncdf option, every band is another data point on the z-axis (time dimension). – Corbjn Dec 15 '20 at 22:43
3

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

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.