Workshop

Workshop 1 Global Mean Surface Temps

load necessary packages

library(ncdf4);
library(ncdf4.helpers);
library(maps);
library(RColorBrewer)
library(lattice)
library(latticeExtra)

(not essential) clear variables and any open plots before starting

graphics.off()
rm(list=ls())

Set the working directory where the data is located

setwd("C:/Users/micha/Downloads/R")

Step 1 - Read HadCRUT Data

available from # https://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/gridded_fields/HadCRUT.4.6.0.0.median_netcdf.zip

nc   <- nc_open(file.path(getwd(),"HadCRUT.4.6.0.0.median.nc")); # open netcdf file
tobs <- ncvar_get(nc,"temperature_anomaly"); # load temperature anomaly
ts   <- ncvar_get(nc,"time"); # load the time
lat  <- ncvar_get(nc,"latitude") # load the latitude
lon  <- ncvar_get(nc,"longitude") # load the longitude
# get units, always changeable with time
tunits<-ncatt_get(nc,"time",attname="units")
# tunits = "days since 1850-1-1 00:00:00"
tustr<-strsplit(tunits$value, " ")
# check you are happy with the data time origin
date<-as.character(as.Date(ts,origin=unlist(tustr)[3]))

2. load global mean surface temperature data from giss temp https://data.giss.nasa.gov/gistemp/

gmst = read.csv(file.path(getwd(),'global_temps_monthly.csv'),skip = 1,header=TRUE)
# create an annual average
gmst <- data.frame(Year = gmst$Year, t = rowMeans(gmst[, 2:13], na.rm = T))

plot the global mean temperature data stored in the variable ‘gmst’. Use ‘plot’

plot(gmst,type="l")

3 Visualising the data

year_index <- which(date=="1901-06-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11); # defines the interval of colour plotted
# two lines of code to get nicer colours for the eventual plots
colors <- rev(brewer.pal(10, "RdYlBu"))
rgb.palette <- colorRampPalette(colors)
filled.contour(lon, lat, mapmat, color.palette=rgb.palette, levels=int,
                plot.title=title(main=paste0('HadCRUT Anom. in ', date[year_index]),
                xlab="Latitude",ylab="Longitude", cex.lab=1.5),
                plot.axes={axis(1, cex.axis=1.5);
                axis(2, cex.axis=1.5);
                map('world', add=TRUE);
                grid()},
                key.title=title(main="[oC]"),
                key.axes={axis(4, cex.axis=1.5)})

We want to see data closer to the present. So, we repeat the last step of using the large chunk of code above but replace the date 1901-06-16 to 2017-05-16 as seen below

year_index <- which(date=="2017-05-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11); # defines the interval of colour plotted
# two lines of code to get nicer colours for the eventual plots
colors <- rev(brewer.pal(10, "RdYlBu"))
rgb.palette <- colorRampPalette(colors)
filled.contour(lon, lat, mapmat, color.palette=rgb.palette, levels=int,
               plot.title=title(main=paste0('HadCRUT Anom. in ', date[year_index]),
                                xlab="Latitude",ylab="Longitude", cex.lab=1.5),
               plot.axes={axis(1, cex.axis=1.5);
                 axis(2, cex.axis=1.5);
                 map('world', add=TRUE);
                 grid()},
               key.title=title(main="[oC]"),
               key.axes={axis(4, cex.axis=1.5)})

Repeat the last step of the large chunk of code, but replace the date for more recent data.

year_index <- which(date=="2017-05-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11); # defines the interval of colour plotted
# two lines of code to get nicer colours for the eventual plots
colors <- rev(brewer.pal(10, "RdYlBu"))
rgb.palette <- colorRampPalette(colors)
filled.contour(lon, lat, mapmat, color.palette=rgb.palette, levels=int,
               plot.title=title(main=paste0('HadCRUT Anom. in ', date[year_index]),
                                xlab="Latitude",ylab="Longitude", cex.lab=1.5),
               plot.axes={axis(1, cex.axis=1.5);
                 axis(2, cex.axis=1.5);
                 map('world', add=TRUE);
                 grid()},
               key.title=title(main="[oC]"),
               key.axes={axis(4, cex.axis=1.5)})