Plot the GIS GMST estimate:

plot(gmst, type="l", col="red", main="GMST", xlab="Years", ylab="T Anomaly")
legend("topleft",
       legend=c("GIS TEMP"),
       col=c("red"), lty=1, lwd=1, bty="n")

Find and plot in the same way two months that illustrate inhomogeneities of the dataset through time

year_index <- which(date=="1901-06-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11); # defines the interval of colour plotted

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)})

year_index <- which(date=="2017-07-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11); # defines the interval of colour plotted

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)})

Calculate the mean of the zonal means. Plot the GIS GMST, the mean of the meriodional means, and the mean of the zonal means

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)
}

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

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

plot(gmst, type="l", col="red", main="GMST", xlab="Years", ylab="T Anomaly")
lines(years, mean_merid_means, type="l", col="blue")
lines(years, mean_zonal_means, type= "l", col="green")
legend("topleft",
       legend=c("GIS TEMP", "Mean of the Meridional Means", "Mean of the Zonal Means"),
       fill=c("red", "blue", "green"), lty=1, lwd=1, bty="n")

Subtract the meridional and zonal means from the GIS GMST. Plot the anomalies

year_index <- which(years=="1880")
year_index
## [1] 31
years[year_index]
## [1] "1880"
merid_anomaly <- gmst$t - mean_merid_means[31:168]
zonal_anomaly <- gmst$t - mean_zonal_means[31:168]

plot(years[31:168], merid_anomaly, type="l", col="green", main="GMST", xlab="Years", ylab="T Anomaly", ylim = c(-.2,.4))
lines(years[31:168], zonal_anomaly, type="l", col="blue")
legend("topleft",
       legend=c("Mean of the Meridional Mean Anomaly", "Mean of the Zonal Mean Anomaly"),
       fill=c("green", "blue"), lty=1, lwd=1, bty="n")

It would seem the mean of the meridional means is varying less from GMST than the mean of the zonal means. Therfore, the meridional mean should be treated with more confidence than the zonal mean.

Adding a weighted mean and abline

# 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(lat*pi/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)

year_index <- which(years=="1960")
year_index
## [1] 111
years[year_index]
## [1] "1960"
year_index <- which(gmst$Year=="1960")
year_index
## [1] 81
gmst$Year[year_index]
## [1] 1960
plot(gmst$Year[81:138], gmst$t[81:138], type="l", col="red", main="GMST", xlab="Years", ylab="T Anomaly")
lines(years[111:168], mean_merid_means[111:168], type="l", col="blue")
lines(years[111:168], mean_zonal_means[111:168], type= "l", col="green")
lines(years[111:168], gmst_weighted[111:168], type="l", col="black")
legend("topleft",
       legend=c("GIS TEMP", "Mean of the Meridional Means", "Mean of the Zonal Means", "Weighted Mean"),
       fill=c("red", "blue", "green", "black"), lty=1, lwd=1, bty="n")

# 5. trends. considering data only from 1960, what are the rates of global temperature rise?
abline(lm(gmst_weighted[111:168]~years[111:168]),col="red",lwd=2)
abline(lsfit(years[111:168],gmst_weighted[111:168]), col="red")

The line suggest a stedy rise in global temperatures since 1960.