5

Taking a raster file of monthly temperature data for multiple years which has a name attached accessible via names(object) in the following format 'Jan.1981', 'Feb.1981' etc (example files for two years that works with code below here - adding all files makes it too big.

Reading in and writing this to NetCDF using the following code:

#Load Packages
library(raster)
library(ncdf4)

#Read in temperature files
r1 <- brick('TavgM_1981.grd')
r2 <- brick('TavgM_1982.grd')

#stack them together 
TempStack = stack(r1, r2)

#set the coordinate system (as it was missing)
crs(TempStack) <- ('+proj=lcc +lat_1=53.5 +lat_2=53.5 +lat_0=46.834 +lon_0=5 +x_0=1488375 +y_0=-203375 +datum=WGS84 +to_meter=2500 +no_defs +ellps=WGS84 +towgs84=0,0,0')

#reproject to get in lat/lon instead of meters
TempStack<-projectRaster(TempStack, crs=CRS("+init=epsg:4326"))

#Extract monthly data names to assign to netCDf later
names <- names(TempStack)

#write the raster file to NetCDF
writeRaster(TempStack, "Temp.nc", overwrite=TRUE, format="CDF",     varname="Temperature", varunit="degC", 
        longname="Temperature -- raster stack to netCDF, monthly average", xname="Longitude",   yname="Latitude", zname='Time', zunit=names)

When I write this to NetCDF and plot the monthly data it is organised from month 1 to month 24, but I want it to have 'Jan 1981', 'Feb 1981' etc.

I thought by adding the zunit argument in writeRaster would work, but it doesn't, the numbers are all still 1-24 instead of Jan, Feb etc.

1 Answer 1

11
+50

There are a couple of misconceptions in your example. First, you should realize that the values in a netcdf dimension must be numeric. They are not just labels for layers, they are actual values of that dimension, and cannot therefore take values like "Jan.1980", which is a string. One way around this is to save your netcdf file and then add the z dimension values to it as a numeric value. Unfortunately that means we can't use date/time variable types either, but must first convert them to numeric equivalents. Here I use the lubridate package to do that.

# first we write the netcdf file to disk
writeRaster(TempStack, "Temp.nc", overwrite=TRUE, 
            format="CDF",     varname="Temperature", varunit="degC", 
            longname="Temperature -- raster stack to netCDF, monthly average", 
            xname="Longitude",   yname="Latitude", zname='Time', zunit='seconds')

# and open a connection to it to make changes.
# note that we use write=TRUE so that we can change it
nc = nc_open('Temp.nc', write = TRUE)

# now convert the strings to numeric values based on their dates
zvals = lubridate::parse_date_time(names, orders = 'm.y', tz = "UTC")
zvals = as.integer(zvals)

# and we can write these numeric dates to the z dimension
ncdf4::ncvar_put(nc, 'Time', zvals)

Haing written the dates to the z dimension like this, we will also need to reverse the process if you want to convert the numeric z values back into raster layer names that look like "Jan.1908" etc. Again, lubridate can help.

ncb = brick('Temp.nc')
zvals = ncvar_get(nc, 'Time')
zvals =  as.POSIXct(zvals, origin = lubridate::origin, tz = "UTC")
znames = paste0(lubridate::month(zvals, label=T), '.', lubridate::year(zvals))
names(ncb) = znames

Let's check that worked:

plot(ncb)

enter image description here

Sign up to request clarification or add additional context in comments.

1 Comment

Thank you for this - I did not realise NetCDF dimensions had to be numeric. Converting to UTC means I can read it with Python too, which is a big help. Thanks!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.