This exercise involves investigating the correlation and regression patterns associated with various climate indices. First, let’s load in the required libraries:

library(sp)
library(raster) 
library(ncdf4)
library(rgdal)
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

memory.limit(size=65000)
## [1] 65000

We also need to load a continental data file for plotting purposes. The data was obatined from the Natural Earth website at: http://www.naturalearthdata.com/downloads/110m-physical-vectors/

# Converts into a R Workspace
world.land <- readOGR(dsn = file.path(getwd(),"ne_110m_land"), layer = "ne_110m_land" )
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GY667_2\ne_110m_land", layer: "ne_110m_land"
## with 127 features
## It has 3 fields
# Converts workspace to approapriate projections
world.land <- spTransform(world.land, CRS( "+proj=longlat +datum=WGS84" ) )

Sea Surface Temperature

Now to load the Hadley Centre Sea Surface Temperature dataset, obtained from the Met Office (2019):

nc   <- nc_open(file.path(getwd(),"HadISST_sst.nc"))

#Time
time <- ncvar_get(nc, 'time')

# Latitude and Longitude
lat <- ncvar_get(nc, 'latitude')
lon <- ncvar_get(nc, 'longitude')

# Sea Surface Temperature
sst <- ncvar_get(nc, 'sst')

The SST data may be used to look at monthly and annually averaged sea surface temperature:

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


# Missing Values 
fillvalue <- ncatt_get(nc,"sst","_FillValue")
sst[sst==fillvalue$value] <- NA
missvalue <- ncatt_get(nc,"sst","missing_value")
sst[sst==missvalue$value] <- NA
sst[sst==-1000] <- NA

# Annual Mean
year <- format(as.Date(date, format="%Y-%m-%d"),"%Y")
gmean <- colMeans(sst, na.rm = TRUE, dims=2)
annmean <- aggregate(gmean,by=list(year),FUN=mean,na.rm=TRUE)

avsst = rowMeans(sst,na.rm=FALSE,dims=2)

# Plot Colours
colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)

# Plotting
grid <- expand.grid(x=lon, y=lat)
grid$avsst <- as.vector(avsst)

# Annual Averages
yrs <- annmean$Group.1 
nyr <- length(yrs)

asst <- array(NA,c(dim(lon),dim(lat),nyr)) 
for (k in 1:nyr) {
  asst[,,k] <- rowMeans(sst[,,year==yrs[k]],na.rm=FALSE,dims=2) 
}

# Plotting
grid$an_avsst <- as.vector(rowMeans(asst,na.rm=FALSE,dims=2))

levelplot(an_avsst~x*y, data=grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Annually Averaged SST') + 
  layer(sp.polygons(world.land, fill='white')) 

# Remove the Global Mean 
gmean <- colMeans(asst, na.rm = TRUE, dims=2)

for (k in 1:nyr){
  asst[,,k]<-asst[,,k]-matrix(gmean[k],length(lon),length(lat))
}

Atlantic Multidecadal Oscillation/Variability

The Atlantic Multidecadal Oscillation, or AMO, is a large-scale mode of climate variability affecting the Northern Hemisphere. Simplistically, it is a series of alternating warm and cold phases of sea surface temperature, which occur on a multidecadal scale. This is not to be confused with the North Atlantic Oscillation, a phenomenon which is influenced by the the pressure difference between the Icelandic Low and Azores High. The AMO has been linked to several instances of regional climate variability, including multidecadal rainfall in northeast Brazil and the African Sahel (Knight et al. 2006).

The AMO data used in this exercise was taken fro the KNMI Climate Explorer website (KNMI, 2019). A choice of AMO indices was offered; in the end, the HadSST data was chosen to remain consistent with the global HadISST data used earlier.

i) Sea Surface Temperature

amo_nc   <- nc_open(file.path(getwd(),"amo_hadsst_ts.nc"))

amo.monthly <- ncvar_get(amo_nc,'AMO') # AMO
amo.time <- ncvar_get(amo_nc, 'time') # time
# get units
tunits<-ncatt_get(amo_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")

# The Date Needs to Be Formatted
amo.date <-as.character(as.Date(amo.time*365.25/12,origin=unlist(tustr)[3]))

# Creating a Data Frame
amo.df <- data.frame("year"=format(as.Date(amo.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(amo.date, format="%Y-%m-%d"),"%m"),
                  "amo"=amo.monthly)

# Annual Aggregation
amo.an <- aggregate(amo~year,amo.df,mean)

# Finding The Overlap
dim(asst)
## [1] 360 180 133
yrs[1]
## [1] "1870"
yrs[length(yrs)]
## [1] "2002"
amo.an$year
##   [1] 1854 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891
##  [15] 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905
##  [29] 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919
##  [43] 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933
##  [57] 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947
##  [71] 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961
##  [85] 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975
##  [99] 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989
## [113] 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003
## [127] 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
## [141] 2018 2019
## 166 Levels: 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 ... 2019
amo.an <- amo.an[2:125,]

yrs  <- yrs[10:133]
asst <- asst[,,10:133]

# SST Timeseries 
sst_ts <- amo.an$amo

# A Quick Look 
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main='AMO')

# Calculating Correlation in a Loop
c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], sst_ts)
    p.vals <- cor.test(asst[i,j,], sst_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

# Adding to the Data Frame
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

# Significance
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

# Plot
levelplot(corr~x*y, data=grid ,
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with AMO')) + 
              layer(sp.polygons(world.land, fill='white')) +
              layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

# LM. Linear regression

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~sst_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01) 
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SSTA with AMO') + 
  layer(sp.polygons(world.land, fill='white')) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

ii) Sea Level Pressure

The sea level pressure data was obtained from a previously prepared nc file name slp.mon.mean.nc:

# Loading the Data
slp_nc   <- nc_open(file.path(getwd(),"slp.mon.mean.nc"))

slp <- ncvar_get(slp_nc, 'slp')
lat<-ncvar_get(slp_nc, 'lat')
lon<-ncvar_get(slp_nc, 'lon') 

# Longitude is in the wrong format (0 to 360 degrees). It must be converted to compatible style (-180 to +180 degrees)
lon <- ifelse(lon>180, lon-360, lon)

# Formatting Time
slp.time<-ncvar_get(slp_nc, 'time')
tunits<-ncatt_get(slp_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
slp.date <-as.character(as.Date(slp.time*24/365,origin=unlist(tustr)[3]))

# Formatting Missing Values
fillvalue <- ncatt_get(slp_nc,"slp","_FillValue")
slp[slp==fillvalue$value] <- NA
missvalue <- ncatt_get(slp_nc,"slp","missing_value")
slp[slp==missvalue$value] <- NA
slp[slp==-1000] <- NA

year <- format(as.Date(slp.date, format="%Y-%m-%d"),"%Y")
gmean <- colMeans(slp, na.rm = TRUE, dims=2)
annmean <- aggregate(gmean,by=list(year),FUN=mean,na.rm=TRUE)

avsst = rowMeans(slp,na.rm=FALSE,dims=2)

colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)

grid <- expand.grid(x=lon, y=lat)

grid$avsst <- as.vector(avsst)

yrs <- annmean$Group.1 
nyr <- length(yrs)

asst <- array(NA,c(dim(lon),dim(lat),nyr))

for (k in 1:nyr) {
  asst[,,k] <- rowMeans(slp[,,year==yrs[k]],na.rm=FALSE,dims=2)
}

grid$an_avsst <- as.vector(rowMeans(asst,na.rm=FALSE,dims=2))

# A Quick Look (Pressure in Mb)
levelplot(an_avsst~x*y, data=grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') + 
  layer(sp.polygons(world.land, fill='white')) 

gmean <- colMeans(asst, na.rm = TRUE, dims=2)

for (k in 1:nyr){
  asst[,,k]<-asst[,,k]-matrix(gmean[k],length(lon),length(lat))
}

We’ll repeat the same as before with the AMO code, just substituting the SLP for the SST data:

amo_nc   <- nc_open(file.path(getwd(),"amo_hadsst_ts.nc"))

amo.monthly <- ncvar_get(amo_nc,'AMO')
amo.time <- ncvar_get(amo_nc, 'time')

tunits<-ncatt_get(amo_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")

amo.date <-as.character(as.Date(amo.time*365.25/12,origin=unlist(tustr)[3]))

amo.df <- data.frame("year"=format(as.Date(amo.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(amo.date, format="%Y-%m-%d"),"%m"),
                  "amo"=amo.monthly)

amo.an <- aggregate(amo~year,amo.df,mean)

amo.an <- amo.an[2:101,]
yrs  <- yrs[10:109]
asst <- asst[,,10:109] 

slp_ts <- amo.an$amo

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], slp_ts)
    p.vals <- cor.test(asst[i,j,], slp_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

levelplot(corr~x*y, data=grid ,
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with AMO')) + 
              layer(sp.polygons(world.land, fill='white')) +
              layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~slp_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SLP with AMO') + 
  layer(sp.polygons(world.land, fill='white')) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Analysis

As one might expect, the maps reveal the north Atlantic as having the strongest relationship with the AMO. Indeed, AMO phases have coincided with anomalous climate conditions in this region. For example, the warm AMO phase which occured 1930s-1950s coincided with decreased rainfall in northeast Brazil and increased hurricane formation in the African Sahel (Folland et al., 2001). However, in terms of both SST and SLP projections, there is strong correlation with a number of regions outsude the Atlantic. These include southeast Asia, the southern and southwestern tip of South America and, in the case of sea-level pressure, the southeastern coast of Africa.

Pacific Decadal Oscillation

The Pacific Decadal Oscillation, or PDO, is a pattern of ocean-atmosphere climate variability affecting the mid-latitude Pacific Basin. The PDO usually identified as either warm or cool surface waters in the Pacific Ocean, north of 20°N. The periodicity of the PDO ranges from multiple years to multiple decades (Mantua and Hare, 2002).

Like the AMO data from PART 1, the PDO used in this exercise was taken fro the KNMI Climate Explorer website (KNMI, 2019). A choice of indices was offered; in the end, the HadSST3 data was chosen to remain consistent with the global HadISST data used earlier.

i) Sea Surface Temperature (SST)

# Loading the Data
pdo_nc   <- nc_open(file.path(getwd(),"pdo_hadsst3.nc"))

pdo.monthly <- ncvar_get(pdo_nc,'index')
pdo.time <- ncvar_get(pdo_nc, 'time')

# Formatting Time
tunits<-ncatt_get(pdo_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
pdo.date <-as.character(as.Date(pdo.time*365.25/12,origin=unlist(tustr)[3]))

# Data Frame
pdo.df <- data.frame("year"=format(as.Date(pdo.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(pdo.date, format="%Y-%m-%d"),"%m"),
                  "pdo"=pdo.monthly)

# Annual Aggregation
pdo.an <- aggregate(pdo~year,pdo.df,mean)

# Finding the Overlap
pdo.an <- pdo.an[1:107,]
yrs  <- yrs[15:121]
asst <- asst[,,15:121]

# SST Timeseries 
sst_ts <- pdo.an$pdo

# A Quick Look 
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main='PDO')

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], sst_ts)
    p.vals <- cor.test(asst[i,j,], sst_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

# Add to Plotting Data
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

# Significance
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

# Plot
levelplot(corr~x*y, data=grid , 
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with PDO')) + 
              layer(sp.polygons(world.land, fill='white')) +
              layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

# LM. Linear regression

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~sst_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot
levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SSTA with PDO') + 
  layer(sp.polygons(world.land, fill='white')) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

ii) Sea Level Pressure (SLP)

pdo_nc   <- nc_open(file.path(getwd(),"pdo_hadsst3.nc"))

pdo.monthly <- ncvar_get(pdo_nc,'index')
pdo.time <- ncvar_get(pdo_nc, 'time')

tunits<-ncatt_get(pdo_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
pdo.date <-as.character(as.Date(pdo.time*365.25/12,origin=unlist(tustr)[3]))

pdo.df <- data.frame("year"=format(as.Date(pdo.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(pdo.date, format="%Y-%m-%d"),"%m"),
                  "pdo"=pdo.monthly)

pdo.an <- aggregate(pdo~year,pdo.df,mean)
pdo.an <- pdo.an[1:112,]

slp_ts <- pdo.an$pdo


c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], slp_ts)
    p.vals <- cor.test(asst[i,j,], slp_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)


levelplot(corr~x*y, data=grid , 
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with PDO')) + 
              layer(sp.polygons(world.land, fill='white'))

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~slp_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)


levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SLP with PDO') + 
  layer(sp.polygons(world.land, fill='white'))

Analysis

The north Pacific is recognised as a region of high significance in almost all of the maps. This is to be expected; Mantua and Hare write that, in earlier studies, the climatic fingerprints of the PDO were most visible in this region. Other regions in various parts of the world were recognised, including the south Pacific and Indian Ocean (for SST) and the north Atlantic (for SLP).

El Nino Southern Oscillation

The El Nino Southern Oscillation, also known as ENSO, or just El Nino, refers to the interannual variation in ocean temperatures in the tropical Pacific. The term ‘El Nino’ was originally used to describe a warm current which ran south along the Peruvian and Ecuadorian coasts around Christmastime (the name is Spanish for ‘boy child / Christ’). Now, ‘El Nino’ typically refers to the warm phase of the El Nino Southern Oscillation; that is, a basinwide warming of the tropical Pacific. The opposite phase, known as La Nina (‘little girl’), sees widespread cooling in the tropical pacific. Despite its immediate focus, ENSO has been linked to various climate events around the world (Trenberth, 1997).

The ENSO data used in this exercise was once again taken from the KNMI Climate Explorer website (KNMI, 2019). The site provided a choice of multiple ENSO indices. THe HadSST1 index was chosen to remain consistent with the global ISST data, and with the previous indices. Several other indices, numbered between 1 and 4, were available. These indices are named after ship tracks which cross/have crossed the region. In the end, the Nino3.4a index was chosen. Nino3.4 was created after researchers discovered the key region for ocean-atmosphere interactions lied further west of Nino3 (UCAR, 2019).

i) Sea Surface Temperature

# Loading the Data
enso_nc   <- nc_open(file.path(getwd(),"nino3.4a_hadsst1.nc"))

enso.monthly <- ncvar_get(enso_nc,'Nino3.4')
enso.time <- ncvar_get(enso_nc, 'time')

# Formatting Time
tunits<-ncatt_get(enso_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
enso.date <-as.character(as.Date(enso.time*365.25/12,origin=unlist(tustr)[3]))

# Data Frame
enso.df <- data.frame("year"=format(as.Date(enso.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(enso.date, format="%Y-%m-%d"),"%m"),
                  "enso"=enso.monthly)

# Annual Aggregation
enso.an <- aggregate(enso~year,enso.df,mean)

# Finding the Overlap
enso.an <- enso.an[1:133,]


# SST Timeseries 
sst_ts <- enso.an$enso

# Plot 
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main='ENSO')

# Correlation
c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], sst_ts)
    p.vals <- cor.test(asst[i,j,], sst_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

# Add to Plotting Data
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

# Significance
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

# Plot
levelplot(corr~x*y, data=grid ,
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with ENSO')) + 
              layer(sp.polygons(world.land, fill='white')) +
              layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

# LM. Linear regression

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~sst_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01)
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot
levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SSTA with ENSO') + 
  layer(sp.polygons(world.land, fill='white')) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

ii) Sea Level Pressure

enso_nc   <- nc_open(file.path(getwd(),"nino3.4a_hadsst1.nc"))

enso.monthly <- ncvar_get(enso_nc,'Nino3.4')
enso.time <- ncvar_get(enso_nc, 'time')

tunits<-ncatt_get(enso_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
enso.date <-as.character(as.Date(enso.time*365.25/12,origin=unlist(tustr)[3]))

enso.df <- data.frame("year"=format(as.Date(enso.date, format="%Y-%m-%d"),"%Y"),
                  "month"=format(as.Date(enso.date, format="%Y-%m-%d"),"%m"),
                  "enso"=enso.monthly)

enso.an <- aggregate(enso~year,enso.df,mean)

enso.an <- enso.an[39:150,]

slp_ts <- enso.an$enso

c.matrix <- matrix(NA,length(lon),length(lat))
t.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    c.matrix[i,j] <- cor(asst[i,j,], slp_ts)
    p.vals <- cor.test(asst[i,j,], slp_ts)
    t.matrix[i, j] <- p.vals$p.value
    }
  }
}

grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

levelplot(corr~x*y, data=grid , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with ENSO')) + 
              layer(sp.polygons(world.land, fill='white'))

r.matrix <- matrix(NA,length(lon),length(lat))
s.matrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
      r.lm <- lm(asst[i,j,]~slp_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

levelplot(reg~x*y, data=grid , at=c(-30:30)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main='Regression of SLP with ENSO') + 
  layer(sp.polygons(world.land, fill='white'))

Analysis

High significance is detected off the west and southwest coast of South America. One would expect this, as this is the region most synonymous with El Nino and the El Nino Southern Oscillation (Trenberth, 1997). Other regions identified as having high correlation are the mid to north Pacific (for SST) and the southwest Asia (for SLP).

Comparative Analysis

The maps appear to show similarity between the sea level pressure patterns of the AMO and PDO. These two oscillations occur in two diffrent bodies of water, but perhaps the mechanisms of each have a similar effect on sea level pressure. The AMO maps appear to show region of SLP correlation just off the tip of South America. This may be connected to the El Nino Southern Oscillation, in which case, the same hypothesis could be applied. It is also possible that the similarity is caused by the global effects of ENSO.

In terms of both SST and SLP patterns, a connection appears to be present between the ENSO and the Pacific Decadal Oscillation. Established literature would support this statement. The PDO has, in the past, been recognised as a sort of long-lived version of El Nino climate variability (Mantua and Hare, 2002). A previous study, which looked at El Nino-forced variability on the PDO, found that it was dependent on ENSO across all timescales (Newman et al., 2003) That said, the PDO is in fact a unique oscillation pattern; it has an irregular periodicity which is distint from ENSO, and appears to be influenced by a combination of El Nino and overall atmospheric noise.

References

Folland, C. K., A. W. Colman, D. P. Rowell, and M. K. Davey (2001), Predictability of northeast Brazil rainfall and real-time forecast skill, 1987-98, Journal of Climate, 14, 1937-1958.

Knight, J.R., Folland, C.K. and Scaife, A.A. (2006) Climate impacts of the Atlantic multidecadal oscillation. Geophysical Research Letters, 33(17).

KNMI (2019) Monthly climate indices [online]. Available at: http://climexp.knmi.nl/selectindex.cgi?id=ididsomeone@somewhere (accessed 24 Apr 2019).

Mantua, N.J. and Hare, S.R. (2002) The Pacific decadal oscillation. Journal of oceanography, 58(1), 35-44.

Met Office (2019) Hadley Centre observations datasets [online]. Available at: https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html (accessed 24 Apr 2019).

Natural Earth (2019) 1:110m Physical Vectors [online]. Available at: http://www.naturalearthdata.com/downloads/110m-physical-vectors/ (accessed 24 Apr 2019)

Newman, M., Compo, G.P. and Alexander, M.A. (2003) ENSO-forced variability of the Pacific decadal oscillation. Journal of Climate, 16(23), 3853-3857.

Trenberth, K.E. (1997) The definition of el nino. Bulletin of the American Meteorological Society, 78(12), 2771-2778.

University Corporation for Atmospheric Research (2019) Nino SST Indices (Nino 1+2, 3, 3.4, 4; ONI and TNI) [online]. Available at: https://climatedataguide.ucar.edu/climate-data/nino-sst-indices-nino-12-3-34-4-oni-and-tni (accessed 24 Apr 2019).