Load the relevant packages

library(ncdf4);
## Warning: package 'ncdf4' was built under R version 3.5.2
library(ncdf4.helpers);
## Warning: package 'ncdf4.helpers' was built under R version 3.5.2
library(maps);
## Warning: package 'maps' was built under R version 3.5.2
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 3.5.2
library(lattice)
## Warning: package 'lattice' was built under R version 3.5.2
library(latticeExtra)
## Warning: package 'latticeExtra' was built under R version 3.5.2

Clear variables and any open plots before starting.

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

Set working directory

setwd("X:/GY667")

Preliminary program reads 1. NetCDF gridded surface temperature data and 2. calculated global mean surface temperatures Store the script and the data in the same folder. Complete the workshop in this folder. We will work further with these data in the workshop 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]))

You can explore the data by checking the dimensions. Two functions are useful: dim, length e.g. dim(tobs) = 72 36 2014, compare this to length(lon), length(lat), length(ts) to see the data are in an array of dimensions longitude x latitude x time

Flag missing data, always get rid of any missing data etc. missing data flags are defined differently for most datasets

tobs[t<-90] = NA;
tobs[t>490] = NA;
  1. 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))

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

plot(gmst, type = "l", col = "red", lwd = 3, xlab = "Years", ylab = "T anomaly", main = "GMST")
legend("topleft", col='red', lwd = 1,
    legend="GIS TEMP")

From the output there appears to be an increasing trend in global mean temperature from the 1960s onwards.

  1. taking the gridded hadcrut data, stored as tobs, create annually averaged gridded data if you look at dim(tobs) = 72 x 36 x 2018, these are the dimensions of longitude, latitude, and time you can check this by looking at length(ts), length(lat), length(lon)

View the dimensions of the data

dim(tobs)
## [1]   72   36 2014
dim(ts)
## [1] 2014
dim(lat)
## [1] 36
dim(lon)
## [1] 72

3a. visualising the data each time slice is a month of gridded surface data pick out an index, say 1936-08-16 (the dates end in 15 or 16 so visually inspect) or “1896-07-16” or “2017-06-16”

year_index <- which(date=="1915-08-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)})

From the output it is clear that there is not much global data collected from 1915, and Europe tends to be quite cool.

TODO: repeat the same for a different date. different date for each class give them skeleton script, add limits etc… Choose a date thats more recent as it has more data coverage. Use date[2000:2014] and choose a date from the output This will then be compared to the old map

year_index <- which(date=="2016-08-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)})

It is clear from the map that there is greater global coverage. Similarly there is a widespread increase in temperature which is clearly visable.

There are plenty of commands there, the main things you might want to play with are the main [title], xlab, ylab, and int [defined above] we’ll use annual averages

allyears <- format(as.Date(date, format="%Y-%m-%d"),"%Y") # this produces a year for each timestamp i.e. 12 1850s, 12 1851s etc..
years <- unique(allyears) # this extracts the unique year i.e. 1850, 1851, 1852,...

we’re going to loop through each year and average the monthly values together functions used: mean, creates an average, e.g. mean(c(1,2,3)) = 2 apply, useful when working on arrays to tell R what dimension to average over for, for loops are fundamental in programming, e.g. for (i in 1:10) {cat(i,‘’)} will print the value of i to the screen pre allocate a blank array, filling with NA, this is for speed and is good practice

tann <- array(NA,c(dim(tobs)[1:2],length(years))) 
for (i in 1:length(years)){
  tann[,,i]<-apply(tobs[,,allyears==years[i]],c(1,2),mean,na.rm=TRUE)
}

Annually averaged data is now stored in ‘tann’ use dim() function to check its size or dimensions

dim(tann) # The output corresponds to latitude(36) longitude (72) and years (168).
## [1]  72  36 168
  1. creating a global mean surface temperature timeseries 4a. ignorant averaging

dimensions of tann are long x lat x year to get a zonal average want to average so that the output is lat x year (Longitude x years) to get a meridional avarage, want output to be lon x year

merid_means<-apply(tann,c(1,3),mean,na.rm=TRUE)
mean_merid_means <- apply(merid_means,2,mean,na.rm=TRUE)

TO DO: calculate the mean of the zonal means (Latitude x years)

zonal_means<-apply(tann,c(2,3),mean,na.rm=TRUE)
mean_zonal_means <- apply(zonal_means,2,mean,na.rm=TRUE)

note that apply(a_martrix,2,mean) is equivalent to colMeans(a_matrix)

TO DO. plot the GIS GMST, the mean of the meridional means, and mean of zonal means on one plot

plot(gmst, type = "l", col = "red", xlab = "Years", ylab = "T anomaly")
lines(years, mean_zonal_means, type="l", col = "dark green")
lines(years, mean_merid_means, type="l", col = "blue")

legend("topleft",legend = c("GIS TEMP","Mean of Zonal Means","Mean of the Meridional Means"),col =c("red","dark green","blue"),lty=1,lwd = 1,bty="n")

From this graph there is a clear increase in temperature from 1960s onwards.

TO DO: plot the zonal and meridional means as differences relative to the GIS temp get on same timestamp, gis runs from 1880 to 2017; hadcrut from 1850 to 2017

merid_diff <- gmst$t - mean_merid_means [31:168]
zonal_diff <- gmst$t - mean_zonal_means [31:168]


plot(years [31:168], merid_diff, type = "l", col = "blue", main = "GMST", xlab = "years", ylab = "Temp Anomaly", ylim = c(-0.2, 0.4))
lines(years [31:168], zonal_diff, type = "l", col = "red")

legend("topleft",legend = c("Mean of Zonal Means Anom","Mean of the Meridional Means Anom"),col =c("Red","Blue"),lty=1,lwd = 1,bty="n")

4b. area weighted averaging. we will calculate a cosine weighted average

w<-replicate(length(years),t(replicate(length(lon),cos(lat*pi/180))))

test with an array of ones i.e. the answer should be 1 gmst_weighted_merid <- apply(array(1,dim(tann))w,c(1,3),sum,na.rm=TRUE)/sum(cos(latpi/180))

gmst_weighted_merid <- apply(tann*w,c(1,3),sum,na.rm=TRUE)/sum(cos(lat*pi/180))
gmst_weighted <- apply(gmst_weighted_merid,2,mean,na.rm=TRUE)

Create a plot with GIS Temp, Mean of Zonal Means, Mean of the Meridional Means and Weighted mean

plot(gmst, type = "l", col = "red", xlab = "Years", ylab = "T anomaly")
lines(years, mean_zonal_means, type="l", col = "green")
lines(years, mean_merid_means, type="l", col = "blue")
lines(years, gmst_weighted, type ="l", col = "black")

legend("topleft",legend = c("GIS TEMP","Mean of Zonal Means","Mean of the Meridional Means", "Weighted Mean"),col =c("red","green","blue", "black"),lty=1,lwd = 1,bty="n")


# trends. considering data only from 1960, what are the rates of global temperature rise?

abline(lm(gmst_weighted[111:168]~years[111:168]),col="blue",lwd=2)
## Warning in abline(lm(gmst_weighted[111:168] ~ years[111:168]), col =
## "blue", : only using the first two of 58 regression coefficients
abline(lsfit(years[111:168],gmst_weighted [111:168]),col="black",lwd=2)

abline(lm(mean_zonal_means[111:168]~years[111:168]),col="green",lwd=2)
## Warning in abline(lm(mean_zonal_means[111:168] ~ years[111:168]), col =
## "green", : only using the first two of 58 regression coefficients
abline(lsfit(years[111:168],mean_zonal_means [111:168]),col="green",lwd=2)

abline(lm(mean_merid_means[111:168]~years[111:168]),col="blue",lwd=2)
## Warning in abline(lm(mean_merid_means[111:168] ~ years[111:168]), col =
## "blue", : only using the first two of 58 regression coefficients
abline(lsfit(years[111:168],mean_merid_means [111:168]),col="blue",lwd=2)

From the output it is clear that there is an increase from 1960 onwards, which is demonstrated in the graph.