Getting started with NetCDF – reading CEH GEAR

Last year CEH released an interpolated precipitation dataset for the UK call GEAR. You can read about this dataset here, or view slides presented at the 2014 BHS Symposium here. In short; these data describe precipitation over Great Britain (yes, excluding Northern Ireland) at a 1 km resolution grid with a daily time step.

I’m using these CEH GEAR data as input to a snow cover and melt model I’m working on. I’ll be using this model to show snow melt potential on a grid basis. A reasonable hurdle has been interacting with the format the GEAR data are archived in. GEAR is stored using NetCDF, with each file representing a different year. Within each file are precipitation and QA data ordered by three dimensions; easting, northing and time. You can read NetCDF data in software like QGIS or using a programming language like Python or R. If you’re a regular reader, you’ll have guessed I’m going to write about accessing NetCDF using R.

First things first, make sure your computer is ready to work with NetCDF, you’ll need to install some software – see this page for a handy script. After that install David Pierce’s ncdf4 package:

install.packages("ncdf4")
library(ncdf4)

You’re now ready to go!

Within the .nc file there are a number of data groups you’ll want to extract, the key ones are:

• Easting
• Northing
• Precipitation amount.

Others you may be interested in are:

• Minimum distance to input observation
• Coordinate reference system
• Latitude
• Longitude
• Time

You can see full details on these by connecting to the NetCDF file and printing a summary:

nc = nc_open("~/dir/dir/CEH_GEAR_daily_GB_1983.nc")
nc

Now you can start to read in these data. To take a single day of observations as a grid slice:

# Read x, y and data
precip = list()
precip$x = ncvar_get(nc, "x")
# Invert northings
precip$y = sort(ncvar_get(nc, "y"))
# Read one day of precipitation data
precip$z = ncvar_get(nc, "rainfall_amount", count=c(-1, -1, 1))
# Flip matrix
precip$z = precip$z[, ncol(precip$z):1]
# Convert to a raster
precip.r = raster(precip)

Note the call to the rainfall variable includes a ‘start’ and ‘count’ command, these detail which dimensions to include (start=where to begin and count=how many), with -1 meaning all and positive integers being an index. The start and count vectors are organised as x, y, time; in this instance we’re reading in the full spatial extent of rainfall, but only for the first time step (1st Jan). You can also see I’ve inverted the northings and flipped the matrix, this is because they’re backwards (or R is reading them backwards).

You can now plot this raster, perform stats on it or whatever else takes your fancy:

plot(precip.r)
# e.g.
hist(precip.r)
boxplot(precip.r)

CEH GEAR precipitation for 1st January 1983.

You might also have a single location you want to extract precipitation for. You can do this through the indexing of rainfall_amount; find the index of the coordinates you want in precip$x and precip$y (remember that you inverted your northings) and subsitute these in the relevant part of your call to rainfall_amount. The following example extracts the precipitation for Ben Nevis (216600, 771200) and plots the result:

# Extract Ben Nevis time series
ben = ncvar_get(nc, "rainfall_amount", start=c(218, 480, 1),
   count=c(1, 1, -1))
# Plot precipitation against day of year
plot(ben, type="h", main="Ben Nevis precipitation (1983)",
   xlab="Day of year", ylab="Daily precipitation (mm)", col="darkblue")

Ben Nevis precipitation for 1983 from CEH GEAR data.