To investigate the correlation and regression patterns associated with various climate indices.

Climate indices

  • Atlantic Multidecadal Oscillation/Variability
  • Pacific Decadal Oscillation
  • El Nino/Southern Oscillation

Calculate the correlation and regression patterns associated with each

  • for SST
  • for sea level pressure (SLP)

Begin by loading the relevant packages and clearing variables and plots

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

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

Set relevant working directory and file path.

setwd("F:\\MSc\\GY667 Ocean\\Assignment 2") 

path <- file.path(getwd(),"data")

In order to map the data, need to load a spatial dataframe.

load("F:/MSc/GY667 Ocean/Assignment 2/world.coast.Rdata")

Before calculating correlation and regression patterns, the sea surface temperature and sea level pressure must be loaded and annual averages calculated.

Sea Surface Temperatures

Load the Hadley sea surface temperature data netCDF file.

nc <- nc_open(file.path(getwd(),"HadISST_sst.nc"))
lat <- ncvar_get(nc,"latitude") 
lon <- ncvar_get(nc,"longitude")
sst <- ncvar_get(nc,"sst")
date <- ncvar_get(nc,"time")

Convert the date format to something that R recognises

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

Replace missing values with NA, many datsets have different fillvalues

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

In order to look at global variability, the global mean must be examined. Investigate annual averages with aggregate using colMeans and rowMeans.

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)

Use RColorBrewer to assign a colour palette for plotting this data.

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

The data can be plotted using levelplot. However, levelplot() can only plot data in a dataframe. Hence, in order to plot the SST data with longitude and latitude, the latitude and longitude must be added to a dataframe. expand.grid() can be used to do this and the average SSTs can be added to this dataframe as a vector and plotted using levelplot().

grid <- expand.grid(x=lon, y=lat)
grid$avsst <- as.vector(avsst)
levelplot(avsst~x*y,grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SST'
) + 
  layer(sp.lines(world.coast))

In order to look at annual averages rather than short term variability, annual averages must be examined. An array can be created using a loop that contains the annual SST data calculated from the monthly data.

yrs <- annmean$Group.1  # Extract years for the data
nyr <- length(yrs)      # Determine how many years there are

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

This array contains annual averages i.e. the overall average remains the same as the previous plot. This data can be added to the dataframe formed with longitude, latitude, and monthly averages and plotted using levelplot().

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.lines(world.coast))

Sea level pressure

Load the sea level pressure netCDF data.

slp_nc   <- nc_open(file.path(getwd(),"slp.mon.mean.nc"))
slp_lat <- ncvar_get(slp_nc,"lat") 
slp_lon <- ncvar_get(slp_nc,"lon")
slp <- ncvar_get(slp_nc,"slp")

NOAA measure longitude from 0 to 360 degrees east rather than -180 to 180 as the rest of the data used for this exercise. In order to match the longitudes of all the data, 360 must be subtracted from all values greater than 180 for the sea level pressure data. This will convert these values to agree with the longitude data from the indices and the spatial dataframe.

slp_lon <- ifelse(slp_lon<=180,slp_lon,slp_lon-360)

The time data for sea level pressure must be loaded. Unlike the SST data, the time for SLP is given in hours rather than days. Hence, the time must be converted to days to be recognised by the as.Date() function.

slp_date <-ncvar_get(slp_nc,"time")
tunits<-ncatt_get(slp_nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
slp_date <- (slp_date)/24    # Is now days since origin rather than hours
slp_date<-as.character(as.Date(slp_date,origin=unlist(tustr)[3]))

Replace the missing SLP values with NA.

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

In order to look at global variability, the global mean must be examined. Investigate annual averages with aggregate using colMeans and rowMeans.

slp_year <- format(as.Date(slp_date, format="%Y-%m-%d"),"%Y")
slp_gmean <- colMeans(slp, na.rm = TRUE, dims=2)
slp_mean <- aggregate(slp_gmean,by=list(slp_year),FUN=mean,na.rm=TRUE)

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

The average SLP can be plotted using levelplot(). As with the SST data, longitude, latitude, and average SLP must be combined in a dataframe to be recognised by levelplot().

grid2 <- expand.grid(x=slp_lon, y=slp_lat)
grid2$avslp <- as.vector(avslp)
levelplot(avslp~x*y,grid2,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SLP'
) + 
  layer(sp.lines(world.coast))

Again, an array containing the annual SLP values can be created using a loop.

slp_yrs <- slp_mean$Group.1 
slp_nyr <- length(slp_yrs)

aslp <- array(NA,c(dim(slp_lon),dim(slp_lat),slp_nyr))
for (k in 1:slp_nyr) {
  aslp[,,k] <- rowMeans(slp[,,slp_year==slp_yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}

The annual averages can be added to the existing data frame and plotted.

grid2$an_avslp <- as.vector(rowMeans(aslp,na.rm=FALSE,dims=2))
#plot
levelplot(an_avslp~x*y, data=grid2,col.regions = pal(100),xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') +
  layer(sp.lines(world.coast))

Climate indices

Monthly time series for the three climate indices can be downloaded from http://climexp.knmi.nl/selectindex.cgi?id=ididsomeone@somewhere

Atlantic Multidecadal Oscillation

Sea surface temperatures

Begin by loading the data.

amo_nc <- nc_open("iamo_hadsst_ts.nc")

amo_time <- ncvar_get(amo_nc, "time")

The date format must be converted to something that R recognises. In this case, the time component of the netCDF file is given in months rather than days. This must be converted to days for the as.Date() formula to recognise it. The origin date can also be extracted from the netCDF file using ncatt_get() to extract the attribute and strsplit() to separate the date from the description.

amounits <- ncatt_get(amo_nc,"time",attname="units")
amoustr <- strsplit(amounits$value, " ")
amo_time <- as.character(as.Date(amo_time*365.25/12,origin=unlist(amoustr)[3]))

The index itself can then be loaded and missing values replaced with NA.

amo <- ncvar_get(amo_nc, 'AMO')

#replace missing values

fillvalue <- ncatt_get(amo_nc,"AMO","_FillValue")
amo[amo==fillvalue$value] <- NA
missvalue <- ncatt_get(amo_nc,"AMO","missing_value")
amo[amo==missvalue$value] <- NA
amo[amo==-1000] <- NA

Annual averages must then be extracted to investigate the global mean.

amo_year <- format(as.Date(amo_time, format="%Y-%m-%d"),"%Y")
amo_mean <- aggregate(amo,by=list(amo_year),FUN=mean,na.rm=TRUE)

The index values and years can then be extracted from the mean data.

amo_yrs <- amo_mean$Group.1  # Years

amo_ts <- amo_mean$x         # Index values

This time series can be plotted simply using plot().

plot(amo_yrs,amo_ts,type="l",xlab="Year",ylab="Annual AMO",main="Annual AMO")

To calculate the correlation and regression patterns for the AMO with SST, the two datasets must run for the same time period. Only the overlapping period for the two datasets will be taken. The two datasets overlap from 1879 - 2017.

amo_ts <- amo_mean$x[26:164]
asst1 <- asst[,,10:148]

Correlation

The correlation can be calculated using a loop and added to a matrix. Significance levels (p-values) of the correlation can also be calculated in the same 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(asst1[i,j,][!is.na(asst1[i,j,])])>2){
      c.matrix[i,j] <- cor(asst1[i,j,], amo_ts)
      p.vals <- cor.test(asst1[i,j,], amo_ts)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

The correlation coefficients and p-values must then be added to the existing data frame for plotting.

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

Assuming a 99% significance level, the data frame can be subsetted where the p-value is less than 0.01. This will extract all data that is significant at the 99% level.

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

This subset can be converted to a spatial points dataframe for plotting with correlation.

sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

The correlation patterns can be plotted using levelplot with the points of significance added as a layer.

levelplot(corr~x*y, data=grid , 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude', main="Correlation of SST with AMO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The strongest significant positive correlation between the AMO and SSTs is evident in the North Atlantic. As the AMO signal denotes patterns of SST variability in the Atlantic, a positive AMO is linked to warmer SSTs in the North Atlantic (McCarthy et al., 2015). However, the relationship between the AMO and SSTs in most of the South Atlantic is weak and not statistically significant. This may be due to the Atlantic Meridional Overturning Circulation (AMOC) transporting energy from the South Atlantic towards the equator, leading to colder SSTs here regardless of AMO phase. Furthermore, a stronger (weaker) AMOC may lead to more (less) transport of warm water northwards, resulting in a positive (negative) AMO.

There is also a significant positive correlation between the AMO and SSTs in the Western Pacific. However, this correlation is weak (~0.2-0.3). In contrast, the correlation between the AMO and SSTs in the Eastern Pacific is significantly negative. Levine et al. (2017) found that AMO phase affects the Walker Circulation over the Pacific Ocean, leading to a relationship between the AMO and ENSO.

Regression

Similar to correlation, the regression patterns and their significance can be calculated using a loop and added to the existing data frame for plotting.

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(asst1[i,j,][!is.na(asst1[i,j,])])>2){
      r.lm <- lm(asst1[i,j,]~amo_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)

Again, points of significance can be subsetted and added to a spatial points dataframe.

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

The regression patterns and significance levels can be plotted using levelplot.

levelplot(reg~x*y, data=grid , at=c(-16:37)/10,
          col.regions = pal(150),xlab='Longitude',ylab='Latitude',main="Linear Regression of SST with AMO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

Again, the AMO can be seen to have the strongest positive effect on SSTs in the North Atlantic, particularly in the North Western Atlantic off the coast of Newfoundland. This area (the Labrador Sea) is a region of deep water formation. The AMOC transports warm water north to this point. As, McCarthy et al. (2015) discuss, there is a widely hypothesised connection between ocean circulation and the AMO as ocean circulation determines ocean heat content. Hence, the strong effect of the AMO on SSTs in the North Western Atlantic may be linked to the AMOC. The AMO also has a weak positive effect on SSTs in the Western Pacific and a weak negative effect on those in the Eastern Pacific, off the west coast of South America. The AMO may have an effect on ENSO. Levine et al. (2017) found that SSTs in the tropical Atlantic, which are affected by the AMO, alter the Walker Circulation over the Pacific, which determines ENSO phase. The above map indicates that this effect is a weak negative one i.e. a positive AMO phase may lead to a La Nina event in the Pacific.

The AMO evidently has the largest impact on SSTs in the North Atlantic but does have a slight effect on SSTs in the Pacific in particular.

Sea level pressure

The AMO data overlaps with the SLP data for a different period than SST. Hence, the AMO time series needs to be reassigned.

amo_yrs <- amo_mean$Group.1
amo_ts <- amo_mean$x

plot(amo_yrs,amo_ts,type="l",xlab="Year",ylab="Annual AMO")

amo_ts <- amo_mean$x[95:165]
aslp1 <- aslp

Correlation

The correlation and significance levels can be calculated using a similar loop as used for the SST patterns. The data can be added to a grid and the correlation and points of significance plotted using levelplot().

c.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
t.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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

grid2$corr <- as.vector(c.matrix)   # Add to data frame
grid2$pval <- as.vector(t.matrix)   # Add to data frame

# Significance at 99% level
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

# plot
levelplot(corr~x*y, data=grid2 , 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude', main="Correlation of SLP with AMO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The AMO has a significantly strong positive relationship with SLP over the Pacific, the south of South America, the South Atlantic, and Greenland. So, when the AMO is in a positive phase, SLP over these regions is higher than when the AMO is in a negative phase. The AMO has a significantly strong negative relationship with SLP over the North Atlantic, Eastern Europe, and the Northern Indian Ocean. When the AMO is in a positive phase, SLP over these regions is lower than when the AMO is in a negative phase.

Air temperatures affect SLP; high air temperatures lead to lower sea level pressure as the upward movement of air over a warm surface decreases the pressure at the surface. As SSTs affect the air temperature above the ocean surface, SSTs also affect SLP. Thus, areas, such as the North Atlantic, that experience higher SSTs during a positive AMO phase will experience lower SLP during a positive AMO phase.

Regression

The same method can be used for linear regression.

r.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
s.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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

grid2$reg <- as.vector(r.matrix)  # Add to data frame
grid2$sig <- as.vector(s.matrix)  # Add to data frame

# Significance at 99% level
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

# plot
levelplot(reg~x*y, data=grid2 , at=c(-90:101)/10,
          col.regions = pal(200),xlab='Longitude',ylab='Latitude',main="Linear Regression of SLP with AMO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The AMO has a significant positive effect on SLP over the the Pacific, South Atlantic and over Greenland. The AMO has a significant negative effect on SLP over the North Atlantic, Eastern Europe, and the Northern Indian Ocean. This may be due to the effect the AMO has on SSTs resulting in the AMO having an inverse effect on SLP over a region. Hence, as a positive AMO leads to warmer SSTs over the North Atlantic, for example, a positive AMO also results in lower than normal SLP over this region. However, as regions such as the South Atlantic experience lower SSTs during a positive AMO, they experience higher SLP during this phase.

Pacific Decadal Oscillation

Sea surface temperatures

Plotting the correlation and regression patterns for SSTs and the Pacific Decadal Oscillation is carried out using the same method used for the AMO.

pdo_nc <- nc_open("ipdo_hadsst3.nc")

pdo_time <- ncvar_get(pdo_nc, "time")

# Convert the date format to something that R recognises

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

# get the data, amo
pdo <- ncvar_get(pdo_nc, 'index')

# replace missing values with NA, many datasets have different fillvalues 
fillvalue <- ncatt_get(pdo_nc,"index","_FillValue")
pdo[pdo==fillvalue$value] <- NA
missvalue <- ncatt_get(pdo_nc,"index","missing_value")
pdo[pdo==missvalue$value] <- NA
pdo[pdo==-1000] <- NA

# use annual averages with aggregate and using colMeans and rowMeans
# you could alternatively use apply as in the first workshop
pdo_year <- format(as.Date(pdo_time, format="%Y-%m-%d"),"%Y")    # Need it for aggregation
pdo_mean <- aggregate(pdo,by=list(pdo_year),FUN=mean,na.rm=TRUE)

# Extract the years from pdo
pdo_yrs <- pdo_mean$Group.1

pdo_ts <- pdo_mean$x

# look at the time series

plot(pdo_yrs,pdo_ts,type="l",xlab="Year",ylab="Annual PDO", main="Annual PDO")

# Need to find overlapping years for amo and sst
# pdo 1899 - 2019
# sst 1870 - 2017

pdo_ts <- pdo_mean$x[16:134]
asst2 <- asst[,,30:148]

Correlation

# Correlation between SST and 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(asst2[i,j,][!is.na(asst2[i,j,])])>2){
      c.matrix[i,j] <- cor(asst2[i,j,], pdo_ts)
      p.vals <- cor.test(asst2[i,j,], pdo_ts)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

# add to the same data frame as earlier for plotting
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01.
# Convert to a spatial points data frame.
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot the correlation patterns
levelplot(corr~x*y, data=grid, 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude', main="Correlation of SST with PDO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The correlation between the PDO and SSTs is strongest in the Pacific, as the PDO denotes patterns of SST anomalies in the Pacific (Newman et al, 2016). There is a very strong, positive, and significant relationship (~0.8) between the PDO and SSTs in the midlatitude Pacific while a very strong, negative, and significant relationship (~-0.8) can be seen in the Western Pacific. During a positive PDO, the Western Pacific is cooler than normal while the Eastern Pacific is warmer, hence the displayed relationship between the PDO and SSTs in the Eastern and Western Pacific.

The PDO is influenced by several factors such as the Aleutian Low and the Kuroshio Current (Newman et al., 2016). Newman et al. (2016) found that the strength and location of the Kuroshio Current impacts the PDO and, hence, SSTs. The above map shows a negative correlation between the PDO and SSTs in the North-Western Pacific i.e. the position of the Kuroshio. The Kuroshio transports warm ocean water from the equator northwards hence, when the PDO is in a positive phase, the correlation patterns calculated here might suggest that the Kuroshio is cooler than normal.

There is also some significant positive correlation between the PDO and SSTs in the Indian Ocean and over some of the Western Atlantic.

Regression

# Linear Regression between SST and PDO

# using 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(asst2[i,j,][!is.na(asst2[i,j,])])>2){
      r.lm <- lm(asst2[i,j,]~pdo_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}


# add to the data frame for plotting
grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01.
# Convert this subset to a spatial points data frame.
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot

levelplot(reg~x*y, data=grid,   at=c(-8:8)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main="Linear Regression of SST with PDO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The regression for the PDO with SSTs displays a similar pattern to that of their correlation patterns. The PDO can be seen to have the greatest impact on SSTs in the Pacific Ocean. It has a positive effect on SSTs in the Eastern Pacific and a negative effect on those in the Western Pacific i.e. a positive PDO leads to warmer SSTs in the Eastern Pacific and cooler SSTs in the Western Pacific. This effect may be due to the strength and position of the Kuroshio current (Newman et al., 2016), with a weaker, cooler (warmer) Kuroshio during a warm (cool) PDO phase. The PDO also displays a signifcant positive effect on SSTs in the Northern Indian Ocean i.e. a positive PDO leads to warmer than normal SSTs in the Indian Ocean.

Sea level pressure

The correlation and regression patterns for SLP and the PDO can be plotted in the same manner.

# Extract the years and time series

pdo_yrs <- pdo_mean$Group.1

pdo_ts <- pdo_mean$x

# plot the time series

plot(pdo_yrs,pdo_ts,type="l",xlab="Year",ylab="Annual PDO")

# Need to find overlapping years for pdo and slp
# pdo 1899 - 2019
# sst 1948 - 2018

pdo_ts <- pdo_mean$x[65:135]
aslp2 <- aslp

Correlation

# Correlation between SLP and PDO

c.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
t.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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

# add to the data frame
grid2$corr <- as.vector(c.matrix)
grid2$pval <- as.vector(t.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01
# Convert this subset to a spatial points data frame.
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)


#plot
levelplot(corr~x*y, data=grid2, 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude', main="Correlation of SLP with PDO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

There is a significant strong negative relationship between the PDO and SLP over the North Pacific with a slighty weaker positive correlation over the South Pacific and Antarctica. There is a significant strong positive relationship between the PDO and SLP over the Central Atlantic, Africa, Southern Asia, and Australasia.

When SLP is lower (higher) than normal over the North Pacific, the PDO is in a warm (cool) phase. The Central Atlantic, Indian Ocean, and Australasia experience higher (lower) SLP than normal when the PDO is in a warm (cool) phase.

Regression

# Linear Regression between SLP and PDO

# using LM. Linear regression

r.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
s.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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


# add to the data frame
grid2$reg <- as.vector(r.matrix)
grid2$sig <- as.vector(s.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01.
# Convert this subset to a spatial points data frame.
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot

levelplot(reg~x*y, data=grid2,   at=c(-24:13)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main="Linear Regression of SLP with PDO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The PDO has a signifcant strong negative effect on SLP over Antarctica. It also has a significant, though weaker effect on SLP over the North and South-Eastern Pacific. The PDO has a significant strong positive effect on SLP over the Central Atlantic, Africa, the Indian Ocean, South Asia, and the Western Pacific. This may be due to the SSTs and associated air temperatures in these regions as determined by the PDO. A positive PDO leads to higher PDOs over the North and South-Eastern Pacific and lower SSTs over the Western Pacific, resulting in lower SLP of the North and South-Eastern Pacifc and higher SLP over the Western Pacific, for example.

El Nino/Southern Oscillation

Sea surface temperature

Plotting the correlation and regression patterns for SSTs and ENSO is again done using the same methods.

enso_nc <- nc_open("ihadisst1_nino3.4a.nc")

enso_time <- ncvar_get(enso_nc, "time")

# Convert the date format to something that R recognises

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

# get the data, amo
enso <- ncvar_get(enso_nc, 'Nino3.4')

# replace missing values with NA, many datasets have different fillvalues 
fillvalue <- ncatt_get(enso_nc,"Nino3.4","_FillValue")
enso[enso==fillvalue$value] <- NA
missvalue <- ncatt_get(enso_nc,"Nino3.4","missing_value")
enso[enso==missvalue$value] <- NA
enso[enso==-1000] <- NA

# use annual averages with aggregate and using colMeans and rowMeans
# you could alternatively use apply as in the first workshop
enso_year <- format(as.Date(enso_time, format="%Y-%m-%d"),"%Y")    # Need it for aggregation
enso_mean <- aggregate(enso,by=list(enso_year),FUN=mean,na.rm=TRUE)

# Extract the years from pdo
enso_yrs <- enso_mean$Group.1

enso_ts <- enso_mean$x

# look at the time series

plot(enso_yrs,enso_ts,type="l",xlab="Year",ylab="Annual ENSO")

# Need to find overlapping years for enso and sst
# enso 1870 - 2019
# sst 1870 - 2017

enso_ts <- enso_mean$x[1:148]
asst3 <- asst

Correlation

# Correlation between SST and 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(asst3[i,j,][!is.na(asst3[i,j,])])>2){
      c.matrix[i,j] <- cor(asst3[i,j,], enso_ts)
      p.vals <- cor.test(asst3[i,j,], enso_ts)
      t.matrix[i, j] <- p.vals$p.value
    }
  }
}

# add to the same data frame as earlier for plotting
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01
# Convert this subset to a spatial points data frame.
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
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="Correlation of SST with ENSO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

There is a very strong positive relationship (~1 i.e. linear relationship) between ENSO and SSTs in the Central and Eastern equatorial Pacific. SSTs in the Eastern Pacific as a whole exhibit a strong and significant positive relationship with ENSO. In contrast, a strong and significant negative relationship is evident in the Western Pacific. During a positive ENSO (El Nino phase), the Walker circulation weakens and upwelling off the west coast of South America decreases (IPCC, 2007). Hence, the thermocline is deeper than during a neutral phase and cold water is not rising to the surface in this region. As a result, a positive ENSO results in warmer than normal SSTs here while the Western Pacific experiences cooler SSTs. The opposite effect is seen during a negative ENSO phase (La Nina phase). Hence, the relationship between SSTs in the Eastern and Western Pacific and ENSO are contrasting. There is also some signifcant positive correlation between ENSO and SSTs in the Indian Ocean and some of the Western Atlantic.

Regression

# Linear Regression between SST and ENSO

# using 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(asst3[i,j,][!is.na(asst3[i,j,])])>2){
      r.lm <- lm(asst3[i,j,]~enso_ts)
      r.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      s.matrix[i, j] <- smm$coefficients[8]
    }
  }
}


# add to the data frame
grid$reg <- as.vector(r.matrix)
grid$sig <- as.vector(s.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01.
# Convert this subset to a spatial points data frame.
sig <- subset(grid[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot

levelplot(reg~x*y, data=grid,   at=c(-10:15)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main="Linear Regression of SST with ENSO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

The regression patterns for ENSO and SSTs are similar to those for the correlation between ENSO and SSTs. ENSO has the strongest effect on SSTs in the Pacific, with a very strong positive effect in the Eastern Pacific and a very strong negative effect in the Western Pacific. The SSTs in the Central and Eastern equatorial Pacific are particularly affected by ENSO, as this is the region where upwelling is decreased (increased) during an El Nino (La Nina) phase and, hence, where ocean water is warmer (cooler) than normal. ENSO also has a significant positive effect on SSTs in the Indian Ocean.

Sea level pressure

The correlation and regression patterns for SLP and ENSO can be plotted in the same manner.

Correlation

# Extract the mean and time series

enso_yrs <- enso_mean$Group.1

enso_ts <- enso_mean$x

# plot the time series

plot(enso_yrs,enso_ts,type="l",xlab="Year",ylab="Annual ENSO")

# Need to find overlapping years for amo and sst
# enso 1870 - 2019
# slp 1948 - 2018

enso_ts <- enso_mean$x[79:149]

# Correlation between SLP and PDO

c.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
t.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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

# add to the same data frame
grid2$corr <- as.vector(c.matrix)
grid2$pval <- as.vector(t.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01
# Convert this subset to a spatial points data frame.
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot
levelplot(corr~x*y, data=grid2, 
          col.regions = pal(100),xlab='Longitude',ylab='Latitude', main="Correlation of SLP with ENSO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

There is a significant strong negative relationship between ENSO and SLP over the Eastern Pacific. There is a significant strong positive relationship between ENSO and SLP over Australasia, and a slightly weaker positive relationship over the Central Atlantic, the Northern Indian Ocean, and Sub-Saharan Africa.

ENSO includes not only variations in SSTs (El Nino) but also variations in the pressure gradient between the Western and Eastern Pacific (Southern Oscillation). When the pressure gradient weakens i.e. when the SLP is higher over the Western Pacific and lower over the Eastern Pacific than normal, warm water is not drawn westwards towards Australia but instead collects in the Central and Eastern Pacific (Diaz et al., 2001). This weakening in the pressure gradient generally coincides with an El Nino phase. Thus, when ENSO is in a positive (negative) phase, SLP over the Eastern Pacific will be lower (higher) than normal while SLP over the Western Pacific will be higher (lower) than normal.

Regression

# Linear Regression between SLP and ENSO

# using LM. Linear regression

r.matrix <- matrix(NA,length(slp_lon),length(slp_lat))
s.matrix <- matrix(NA,length(slp_lon),length(slp_lat))

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


# add to the same data frame
grid2$reg <- as.vector(r.matrix)
grid2$sig <- as.vector(s.matrix)

# Assuming a 99% significance level, subset the grid data where pval < 0.01
# Convert this subset to a spatial points data frame.
sig <- subset(grid2[, c(1, 2, 5, 6)], pval < 0.01) # 1 = lat, 2 = lon, 5 = corr, 6 = pval
sig <- SpatialPointsDataFrame(coords = sig[, c(1, 2)], data = sig)

#plot
levelplot(reg~x*y, data=grid2,   at=c(-23:20)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',main="Linear Regression of SLP with ENSO") + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))

ENSO has a significant negative effect on SLP over the Eastern Pacific. It has a positive effect over Australasia, the Central Atlantic, the Northern Indian Ocean, and Sub-Saharan Africa. ENSO displays a very strong positive effect on SLP over the South Pacific but this is not significant.

An El Nino phase leads to warmer SSTs over the Eastern Pacific than normal and coincides with a negative Southern Oscillation. This leads to lower SLP over the Eastern Pacific than normal. In contrast, SLP over Australia is higher than normal during an El Nino phase.

Bibliography

  • ClimateDataGuide NCAR/UCAR (2018) 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 18 April 2019).
  • Diaz, H.F., Hoerling, M.P. & Eischeid, J.K. (2001) ENSO variability, teleconnections and climate change, International Journal of Climatology, vol. 21, no. 15, pp. 1845-1862.
  • IPCC, 2007: Climate Change 2007: Synthesis Report. Contribution of Working Groups I, II and III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change [Core Writing Team, Pachauri, R.K and Reisinger, A. (eds.)]. IPCC, Geneva, Switzerland, 104 pp
  • Levine, A.F.Z., McPhaden, M.J. & Frierson, D.M.W. (2017) The impact of the AMO on multidecadal ENSO variability: AMO IMPACTS ON ENSO, Geophysical Research Letters, vol. 44, no. 8, pp. 3877-3886.
  • McCarthy, G.D., Haigh, I.D., Hirschi, J.J., Grist, J.P. & Smeed, D.A. (2015) Ocean impact on decadal Atlantic climate variability revealed by sea-level observations, Nature, vol. 521, no. 7553, pp. 508-510.
  • Newman, M., Alexander, M.A., Ault, T.R., Cobb, K.M., Deser, C., Lorenzo, E.D., Mantua, N.J., Miller, A.J., Minobe, S., Nakamura, H., Schneider, N., Vimont, D.J., Phillips, A.S., Scott, J.D. & Smith, C.A. (2016) The Pacific Decadal Oscillation, Revisited, Journal of Climate, vol. 29, no. 12, pp. 4399-4427.