by Mark R Payne
DTU-Aqua, Kgs. Lyngby, Denmark
http://www.staff.dtu.dk/mpay

Fri Feb 9 11:08:40 2018


This notebook describes how to extract pointwise data from a NetCDF file using the raster package. As the working example, we use the HadISST v1.1 dataset, which is available from the UK Met Office’s Hadley Centre here: https://www.metoffice.gov.uk/hadobs/hadisst/

To get started, download the NetCDF version of the file (click on the “Download Data” button, then scroll all the way to the bottom of the next page) and decompress the file (e.g. using WinZip or 7-Zip on windows). You will also need to make sure that you have the “raster”" package installed, and potentially also the “ncdf4” as well, which you can install as follows:

install.packages("raster")
install.packages("ncdf4")

Setting up the raster brick

The raster package has three basic types of data structures:

Raster can also support higher dimensional data, up to around four dimensions, and can interface to high dimensionality data by focussing on one aspect of it at a time.

In this case we want a brick, which we create as follows

library(raster)
## Loading required package: sp
fname <- "~/Documents/common_data/HadISST_sst.nc" #You will need to specify your location here
HadISST.b <- brick(fname)  
## Loading required namespace: ncdf4

Lets have a look and see what it contains:

HadISST.b
## class       : RasterBrick 
## dimensions  : 180, 360, 64800, 1776  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/mpayne/Documents/common_data/HadISST_sst.nc 
## names       : X1870.01.16, X1870.02.14, X1870.03.16, X1870.04.15, X1870.05.16, X1870.06.16, X1870.07.16, X1870.08.16, X1870.09.16, X1870.10.16, X1870.11.16, X1870.12.16, X1871.01.16, X1871.02.15, X1871.03.16, ... 
## Date        : 1870-01-16, 2017-12-16 (min, max)
## varname     : sst

So from this we can see that we have a 1x1 degree resolution grid, that runs from -180 to +180 in the longtitudinal direction, from -90 to 90 in the latitudinal i.e. it’s global. The time steps are stored as “Date” and run from 1870-01-16 to 2017-12-16. We have have a look at the contents using plot

plot(HadISST.b)

which makes for lots of maps - one for each month. We can focus on a single month like so:

plot(HadISST.b[[1]])

But where is the data?

Filtering out values

There is clearly something strange with this - there is no texture in the plot. The reason is because this data set uses -100 to indicate missing values i.e. land or an absence of data. So the first logical thing to do is to conver these to NAs, so that they aren’t plotted. This works much like it would with a standard matrix:

#For simplicity, we'll just work with a single layer
r <- HadISST.b[[1776]]
r[r < -300] <- NA

#Now lets plot it.
plot(r)

That looks a lot better!! Now we can see, for example, that the tropics are hot and the poles are cold. Who would have though?!!

Subsetting

Next step is to focus things down on an area of interest. The “crop()” function is your friend here:

#An extent object defines a box .e.g a region of interest
ROI <- extent(-20,20,45,65)
ROI
## class       : Extent 
## xmin        : -20 
## xmax        : 20 
## ymin        : 45 
## ymax        : 65
#We can use this ROI to crop the raster down to the area that we are interested in - in this case Europe
r.crop <- crop(r,ROI)
plot(r.crop)