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.
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:
- Precipitation amount.
Others you may be interested in are:
- Minimum distance to input observation
- Coordinate reference system
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)
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")