“GY667 Workshop 1”

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

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

setwd("/Volumes/Hard Drive/GY667_Workshop1")

nc   <- nc_open(file.path(getwd(),"HadCRUT.4.6.0.0.median.nc")); 
tobs <- ncvar_get(nc,"temperature_anomaly");
ts   <- ncvar_get(nc,"time"); 
lat  <- ncvar_get(nc,"latitude") 
lon  <- ncvar_get(nc,"longitude")

tunits<-ncatt_get(nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
date<-as.character(as.Date(ts,origin=unlist(tustr)[3]))

tobs[t<-90] = NA;
tobs[t>490] = NA;

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

To Do Part 1

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

year_index <- which(date=="1901-06-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11);
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)})

To Do Part 2

year_index <- which(date=="1997-01-16")
mapmat= tobs[,,year_index];
int=seq(-6,6,length.out = 11);
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)})

allyears <- format(as.Date(date, format="%Y-%m-%d"),"%Y") 
years <- unique(allyears)

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= 'Year',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"),
       col= c("red","blue","green"), lty=1,lwd=1, bty="n")

Differences of Means to GMST

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= "blue",main= "GMST",xlab= 'Year',ylab='T anomaly', ylim = c(-.2,.4))
lines(years[31:168], zonal_anomaly,type = 'l',col= 'green')
legend("topleft",
       legend= c("Mean of Meridional Means Anom", "Mean of Zonal Means Anom"),
       col= c("blue","green"), lty=1,lwd=2, bty="n")

The mean of meridional means anomaly better represents the GMST temperature anomaly. There is a wider range of variation in the mean of zonal means anomaly which moves only on the rare occassion remains closer to the GMST temperature.

w<-replicate(length(years),t(replicate(length(lon),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(gmst$Year== "1960")
year_index

## [1] 81

gmst$Year[year_index]

## [1] 1960

plot(gmst$Year[81:168], gmst$t[81:168],type = "l",col= "red",main= "GMST",xlab= 'Year',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","GMST Weighted Mean"),
       col= c("red","blue","green","black"), lty=1,lwd=1, bty="n")

abline(lm(gmst_weighted[111:168]~years[111:168]),col="black",lwd=2)

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

abline(lsfit(years,gmst_weighted))