The first thing we do is load required libraries:

library(ncdf4)
library(raster)
library(rasterVis)
library(maps)
library(maptools)
library(viridis)

1. Use the nc_open function to open the SST file for reading and display meta data.

nc <- nc_open("HadISST_sst.nc")

2. Use the nvar_get function to read in the entire dataset.

sst <- ncvar_get(nc, varid = "sst")

3. Use the dim function to get the dimensions of the data.

dim(sst)
## [1]  360  180 1796

4. Use the start and count arguments to read just the first slice.

sst <- ncvar_get(nc, varid = "sst", start = c(1,1,1), count = c(-1, -1, 1))
dim(sst)
## [1] 360 180

5. Experiment with reading in different slices.

sst <- ncvar_get(nc, varid = "sst", start = c(4,5,9), count = c(100, 8, 4))
dim(sst)
## [1] 100   8   4

6. Challenge: Suppose I want a particular month of data, say November 1997, how do I determine the number (index) of the slice required?

yr_start <- 1870
yr_want <- 1997
whole_years <- yr_want - yr_start
mnth_want <- 11
slice_required <-12 * whole_years + mnth_want
sst <- ncvar_get(nc, varid = "sst", start = c(1,1, slice_required), count = c(-1, -1, 1))

Plotting a time slice

1. Check the geographical range for the data by reading in the latitude and longitude meta data.

long <- ncvar_get(nc, varid = "longitude")
lat <- ncvar_get(nc, varid = "latitude")

2. Plot the Nov 1997 time-slice by converting the array of values to a raster.

#need to get rid of the unrelaistic values which have been put in by the climate model - will give an array of colours.
sst[abs(sst) > 50] <- NA
#need to transpose the array - turn the data a 90 degree angle
sst_t <- t(sst)
#Need to then turn this data into a raster, so it can be plotted. 
  #xmn - means on the x axis, the minimum value, xmx means on the x-axis, the maximum value and same goes for ymn and ymx
sst_raster <- raster(sst_t, xmn = -180, xmx = 180, ymn = -90, ymx = 90,
                     crs = CRS("+proj=longlat + ellps=WGS84 + datum=WGS84"))
#not sure what the second line does. This apparently sets the coordinate reference system with the argument: crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
plot(sst_raster)

3. Plot with a more appropriate colourmap.

plot(sst_raster, col = heat.colors(256, rev = "TRUE"))

plot(sst_raster, col = viridis(256))

4. Label the axes using the xlab and ylab arguments and provide a title using the main argument.

plot(sst_raster,
    col = viridis(256),
    xlab = "Longitude",
    ylab = "Latitude",
    main = "Monthly SST for November 1997")

5. Make a prettier plot using the levelplot function.

Plot <- levelplot(sst_raster, col.regions = heat.colors(256, rev = TRUE), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                  #dont know what it is from here
                   at = seq(from = 0, to = 35, by = 1),
                   colorkey = list(title="°C \n"))
                   #colorkey = list(title="\u00B0C \n"))
Plot

Note, the use of the seq function to supply the at argument. Experiment with seq in the console so you understand what it does. Then experiment with it in the levelplot function to see what impact it has on the map.

Try adding the argument contour = TRUE.

What happens if the margin argument is set to TRUE or removed?

Note, the colorkey argument.

6. Experiment with some of the colourmaps from the viridis palette (requires the viridis library).

Plot <- levelplot(sst_raster, col.regions = inferno(256), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                   at = seq(from = 0, to = 35, by = 1), #= this, i think, makes the map less grainy???
                   colorkey = list(title="°C \n")) # = this line labels the key on the right hand side of the map
Plot

Have a read of: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html

7. Improve the map by setting the intervals at which the axes are annotated.

xticks <- seq(from = -180, to = 180, by = 60) #made as a variable
yticks <- seq(from = -90, to = 90, by = 30) # made as a variable
Plot <- levelplot(sst_raster, col.regions = inferno(256), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                   scales = list(x = list(at = xticks),
                                 y = list(at = yticks)),
                   at = seq(from = 0, to = 32, by = 1.0),
                   colorkey = list(title="°C \n"))
Plot

#the xticks and yticks put basically a scale onto the x and y axis, so from -180 means the x axis starts at -180, and to 180 means it goes to 180. by = 60 means it goes up in intervals of 60. The same for the Y axis
#do need to ask about these

Use the seq function to create a set of suitable tick marks for the x-axis and similarly for the y-axis.

Use this form: xticks <- seq(from =, to =, by =)

We give these to the levelplot function in the following form:

scales = list(x = list(at = xticks), y = list(at = yticks))

list is a function for grouping together different (heterogeneous) data objects.

8. Improve the map by fine control of the colour key.

cticks <- seq(from = 0, to = 32, by = 4)
Plot <- levelplot(sst_raster, col.regions = inferno(256), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                   scales = list(x = list(at = xticks),
                                 y = list(at = yticks)),
                   at = seq(from = 0, to = 32, by = 1.0),
                   colorkey = list(title="°C \n",
                                   labels = list (at = cticks),
                                   width = 2.0))
Plot

Use the seq function as before to create a set of suitable tick marks for the colour key.

Add the argument labels = list(at = cticks) to the colorkey argument list.

Adjust the width of the colour key by adding the argument width to the colorkey argument list. Start with a width of 2.0

9. Make the colourmap continuous (smooth) by adjusting the at argument of levelplot

cticks <- seq(from = 0, to = 32, by = 4)
Plot <- levelplot(sst_raster, col.regions = inferno(10240), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                   scales = list(x = list(at = xticks),
                                 y = list(at = yticks)),
                   at = seq(from = 0, to = 32, by = 0.01),
                   colorkey = list(title="°C \n",
                                   labels = list (at = cticks),
                                   width = 1.0))
Plot

10. Challenge: Rotate the map so that the Pacific is in the centre of the map.

Here we will use the extent and rotate functions.

sst_raster
## class      : RasterLayer 
## dimensions : 180, 360, 64800  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -1.8, 30.41759  (min, max)
extent(sst_raster) <- c(0, 360, -90, 90)
sst_raster
## class      : RasterLayer 
## dimensions : 180, 360, 64800  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -1.8, 30.41759  (min, max)
sst_raster_2 <- rotate(sst_raster)
extent(sst_raster_2) <- c(0, 360, -90, 90)
xticks <- seq(from = 0, to = 360, by = 60) #made as a variable
yticks <- seq(from = -90, to = 90, by = 30) 
cticks <- seq(from = 0, to = 32, by = 4)
Plot <- levelplot(sst_raster_2, col.regions = inferno(10240), margin = FALSE,
                   xlab = "Longitude", ylab = "Latitude",
                   main = "Monthly SST for Nov 1997",
                   scales = list(x = list(at = xticks),
                                 y = list(at = yticks)),
                   at = seq(from = 0, to = 32, by = 0.01),
                   colorkey = list(title="°C \n",
                                   labels = list (at = cticks),
                                   width = 1.0))
Plot

Steps:

  1. First, use the extent function to incorrectly set the longitude range to 0-360.

  2. Next, use the rotate function to rotate the map by 180 degrees. It will “think” it is rotating from 0-360 to -180 - 180. But in fact it will be doing to opposite because of the first step.

  3. Finally, we use the extent function again to correctly set the rotated longitude range to 0-360.

  4. You will need to update the x-axis tick marks accordingly when plotting the rotated data.

11. Extra challenge: Add a map of country outlines to the SST map.

Steps:

  1. First, use the map function to create a map of world outlines.

  2. Next, use the map2SpatialLines function to turn this map into an SP object.

  3. Finally, use the layer and sp.lines functions to add this to the map.

#world_outlines <- map("world", plot = FALSE, wrap = c(0, 360))
#world_outlines_sp <- map2SpatialLines(world_outlines, proj4string = CRS("+proj=longlat + ellps=WGS84 + datum=WGS84"))
#world_outlines <- map(map("world", plot = FALSE, wrap = c(0, 360)))
#Plot <- levelplot(sst_raster_2, col.regions = inferno(10240), margin = FALSE,
                  # xlab = "Longitude", ylab = "Latitude",
                   #main = "Monthly SST for Nov 1997",
                   #scales = list(x = list(at = xticks),
                                # y = list(at = yticks)),
                  # at = seq(from = 0, to = 32, by = 0.01),
                   #colorkey = list(title="°C \n", 
                                   #labels = list(at = cticks), 
                                  # width = 1.0))
#Plot+layer(sp.lines(World_outlines_sp, col = "black", lwd = 0.5))