0

I am trying to extract raster layer names from a netcdf file as previously written from a raster stack. Exporting the raster stack to ncdf works fine. For example:

library(raster)
library(ncdf4)
library(RNetCDF)

#Create some rasters (x3)
r1<-raster(system.file("external/test.grd", package="raster"))
r2<-r1*2
r3<-r2*3

#Stack them
rstack<-stack(r1,r2,r3)

#Give each raster layer a name - in this instance years 2014 to 2016
names(rstack)<-c("2014","2015","2016")

#Write out to netcdf format
writeRaster(rstack, "rstack.nc", overwrite=TRUE, format="CDF", varname="Temperature", varunit="degC", 
        longname="Temperature -- raster stack to netCDF", xname="X", yname="Y",zname="Year",
        zunit="numeric")

However, when I read the ncdf file back into R the Z dimension (i.e. Year) is not retained. E.g:

#Open the new netcdf dataset and look at the Z dimention, i.e. "Year"
data.nc<- open.nc("rstack.nc")
Zdim = var.get.nc(ncfile=data.nc,variable="Year")
print(Zdim)
#[1] 1 2 3

So what we get is the band numbers, i.e. 1,2,3. But what I require is the text defined by Year (e.g. 2014,2015,2016) as defined in:

names(rstack)<-c("2014","2015","2016")

Is it possible to do this?? This problem is not new, refer here: https://gis.stackexchange.com/questions/122167/export-band-names-with-netcdf-file-in-r

There are some convoluted workarounds to get what is required but they seem largely inefficient (i.e. converting the stack to a matrix then manipulating it from here). Just wondering if there is a more elegant way without having to write a large amount of extra code and taking up unnecessary RAM.

Community
  • 1
  • 1
Mat W
  • 35
  • 4

1 Answers1

3

I don't think this is the same issue as that other question. NetCDF variables don't have dimnames, so you can't round-trip a raster stack the way you want to here.

But, the Zdim, is the values on the Z dimension - not the names - and at least I would expect you to setZ(rstack, <zdimvals>) before you write it. I'm not experienced enough with using writeRaster for generating 3D vars, but this seems to work. 

library(raster)
library(ncdf4)
library(RNetCDF)

#Create some rasters (x3)
r1<-raster(system.file("external/test.grd", package="raster"))
r2<-r1*2
r3<-r2*3

#Stack them
rstack<-stack(r1,r2,r3)
rstack <- setZ(rstack, 2014:2016)
#names(rstack)<-c("2014","2015","2016")

#Write out to netcdf format
writeRaster(rstack, "rstack.nc", overwrite=TRUE, format="CDF",     varname="Temperature", varunit="degC", 
        longname="Temperature -- raster stack to netCDF", xname="X",   yname="Y",zname="Year",
        zunit="numeric")

## your ncdf4 code was not right, looked like RNetCDF (which is fine)
data.nc<- nc_open("rstack.nc")
Zdim = ncvar_get(data.nc,varid="Year")
#print(Zdim)  ## now it's numeric
##[1] 2014 2015 2016

You might want to explore the units and metadata you use to store those year values.

Finally, this is very confusing but it comes down to the disconnect between NetCDF and more GIS-y models. I don't know any easy way to understand it but raw brutal experience. NetCDF is very general, and very powerful, but the library itself is very low-level and a bit too simplistic. (It's all slabs and slices, there's no raw "index" abstraction, at least not enough to be useful).

The tools provided by raster are extremely high level, and so less flexible. The only other tool that's anywhere near as high level is Ferret, for whatever reasons heavy use of NetCDF just stays either very NetCDF-focussed, or just dumb arrays with basic metadata. Structured objects like these are rare. Generally NetCDF manipulation is better done with the "nc operator" tools, but it's worth exploring a few options I think, and if you can get R's raster to do what you need you're way ahead.

Python is also very popular for NetCDF, and GDAL - but GDAL suffers from the same "GIS-y" perspective as raster. It's complicated, personally I tend to get by 99% with just R and GDAL, but I do use other tools when I need to.

mdsumner
  • 29,099
  • 6
  • 83
  • 91