GY667 Assignment 3

library(sp)
library(raster) 
library(ncdf4)
library(rgdal)

FALSE rgdal: version: 1.3-6, (SVN revision 773)
FALSE  Geospatial Data Abstraction Library extensions to R successfully loaded
FALSE  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
FALSE  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
FALSE  GDAL binary built with GEOS: FALSE 
FALSE  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
FALSE  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
FALSE  Linking to sp version: 1.3-1

library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

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

Set the working directory

setwd("~/Dropbox/GY667 Assignment 3")
path <- file.path(getwd(),"Data");

Load in the coastline data for mapping SST and SLP

load("~/Dropbox/GY667 Assignment 3/Data/world.coast.Rdata")

Sea Surface Temperature

Load in the Hadley Sea Surface Temperature Data

nc<-nc_open(file.path(path,'HadISST_sst.nc'))

Etract the latitude and longitude

lat<-ncvar_get(nc, 'latitude') # latitude
lon<-ncvar_get(nc, 'longitude') # logitude

Load time

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

Checking that happy with data time origin

date<-as.character(as.Date(time,origin=unlist(tustr)[3]))

get the data for sea surface temperture (sst)

sst <- ncvar_get(nc, 'sst')

load HadISST_sst.nc

lat <- ncvar_get(nc, 'latitude')
lon <- ncvar_get(nc, 'longitude')
time <- ncvar_get(nc, 'time')

replace missing values with NA, many datasets 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

use annual averages with aggregate and 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=TRUE,dims=2)

two lines of code to get nicer colours for the eventual plots

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

for this we can use ‘levelplot’

levelplot(avsst,col.regions = pal(100));

expand.grid will get lats and lons added to a data frame expanded on each others size, our avsst that we want to plot can be added as a vector, basically expanding what we had

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

not interested in short term variability so let’s make annual averages

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

asst will contain the annual sst, declare it here to have the correct size, this is to speed up computation

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

use a for loop to load the annual data into asst

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

we’ve created annual averages

so the average looks exactly the same

add to the same data frame as earlier for plotting

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

plot

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

remove the global mean from each year, lets do a traditional loop

gmean <- colMeans(asst, na.rm = TRUE, dims=2)
for (k in 1:nyr){
  asst[,,k]<-asst[,,k]-matrix(gmean[k],length(lon),length(lat))
}

lon0 <- -10.5 #
lat0 <- 51.5 #
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]

look at this timeseries:

plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon0, ',Lat=', lat0))

############# #Subpolar Gyre

lon1 <- -35.5
lat1 <- 59.5
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon1, ',Lat=', lat1))

##############

Load in the first indice: “AMO”

nc<-nc_open(file.path(path,'iamo_ersst.nc'))

amv.monthly <- ncvar_get(nc,'AMO') # AMV
amv.time <- ncvar_get(nc, 'time') # time

get units

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

check you are happy with the data time origin

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

let’s create a data frame with these data

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

aggregate annually

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

chop to the same size

dim(asst)

## [1] 360 180 148

yrs[1]

## [1] "1870"

yrs[length(yrs)]

## [1] "2017"

1870-2017

amv.an$year

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

1880-2019, overlapping period is 1880-2017, cut off end of the amv data, this can be done visually

amv.an <- amv.an[1:138,]

similarly, cutt off the asst data, and yrs

yrs  <- yrs[11:148]
asst <- asst[,,11:148]

now the data are the same size, set

sst_ts <- amv.an$amv

you don’t need to modify much below this line

look at this timeseries:

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

lets calculate the correlation in a loop also

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 the same data frame as earlier for plotting

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

Assuming a 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05)

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 , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste0('Correlation of SSTA with AMV')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

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',
          xlim=c(-120,10),ylim=c(0,80),
          main=paste0('Correlation of SSTA with AMV')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.01, col = "black"))

The correlation of the SST with AMV can be seen clearly in the graph above. High SST levels can be seen around the atlantic area.

using LM. Linear regression

lets calculate the correlation in a loop also

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

add to the same data frame as earlier 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.nao < 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(-15:15)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          xlim=c(-120,10),ylim=c(0,80),
          main=('Regression of SSTA with AMV')) + 
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

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

###Comments The correlation and regression of the SST with AMV can be seen clearly in the graphs above. High SST levels can be seen around the atlantic area, with excessive amounts shown off the coast of Portugal. The regression shows greater significance in the middle of the North Atlantic between Ireland and Canada. The SST influence in the region is dominated by the Atlantic Multidecadal OScillation (AMO) (Gastineau and Frankignoul; 2015). They link the ocean circulation linked withe the atmospheric conditions of AMO to the regression and correlation patterns seen in the maps above.

ENSO

ncenso<-nc_open(file.path(path,'ihadisst1_nino12a.nc'))

timeenso<-ncvar_get(ncenso, 'time') 
tunitsenso<-ncatt_get(ncenso,"time",attname="units")
tustrenso<-strsplit(tunitsenso$value, " ")
dateenso<-as.character(as.Date(timeenso*365.25/12,origin=unlist(tustrenso)[3]))
head(dateenso)

## [1] "1870-01-15" "1870-02-14" "1870-03-16" "1870-04-16" "1870-05-16"
## [6] "1870-06-16"

tail(dateenso)

## [1] "2019-07-17" "2019-08-17" "2019-09-16" "2019-10-17" "2019-11-16"
## [6] "2019-12-17"

get the ENSO data

sstenso <- ncvar_get(ncenso, 'Nino12')

r##eplace missing values with NA

fillvalueenso <- ncatt_get(ncenso,"Nino12","_FillValue")
sstenso[sstenso==fillvalueenso$value] <- NA
missvalueenso <- ncatt_get(ncenso,"Nino12","missing_value")
sstenso[sstenso==missvalueenso$value] <- NA
sstenso[sstenso==3000000000000000000000000000000000] <- NA
summary(sstenso)  

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -2.5295 -0.8158 -0.3026 -0.1794  0.2769  4.3790      11

getting the annual mean of pdo

ensoyear <- format(as.Date(dateenso, format ="%Y-%m-%d"),"%Y")
ensoannmean <- aggregate(sstenso,by=list(ensoyear), FUN= mean, na.rm= TRUE)
head(ensoannmean)

##   Group.1          x
## 1    1870 -1.0840045
## 2    1871 -0.5825608
## 3    1872 -0.9003969
## 4    1873 -0.9435940
## 5    1874 -0.8821288
## 6    1875 -1.0859472

ensoannmean - 1880-2017

colnames(ensoannmean) <- c("Year", "ENSO")    

chop to the same size

dim(asst)

## [1] 360 180 138

yrs[1]

## [1] "1880"

yrs[length(yrs)]

## [1] "2017"

1870-2017

ensoannmean$year

## NULL

1880-2019,overlapping period is 1880-2017, cut off end of the amv data, this can be done visually

ensoannmean <- ensoannmean[1:138,]

similarly, cutt off the asst data, and yrs

yrsenso  <- yrs[1:138]
asstenso <- asst[,,1:138]

now the data are the same size, set

sst_tsenso <- ensoannmean$ENSO

look at this timeseries:

plot(yrsenso,sst_tsenso,type='l',xlab='Year',ylab='SST Anomaly',main='ENSO')

lets calculate the correlation in a loop also

a.matrix <- matrix(NA,length(lon),length(lat))
b.matrix <- matrix(NA,length(lon),length(lat))
                
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
 if (length(asstenso[i,j,][!is.na(asstenso[i,j,])])>2){
   c.matrix[i,j] <- cor(asstenso[i,j,], sst_tsenso)
  p.vals <- cor.test(asstenso[i,j,], sst_tsenso)
 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 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05), 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 , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
 col.regions = pal(100),xlab='Longitude',ylab='Latitude',
 main=paste0('Correlation of SSTA with ENSO')) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

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',
 xlim=c(-80,60),ylim = c(-60,50),
 main=paste0('Correlation of SSTA with ENSO')) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

### using LM. Linear regression ###lets calculate the correlation in a loop also

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(asstenso[i,j,][!is.na(asstenso[i,j,])])>2){
      r.lm <- lm(asstenso[i,j,]~sst_tsenso)
       c.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
    t.matrix[i, j] <- smm$coefficients[8]
   }
 }
}

add to the same data frame as earlier for plotting

grid$reg <- as.vector(c.matrix)
grid$sig <- as.vector(t.matrix)

Assuming a 99% significance level, subset the grid data where pval.nao < 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(-5:5)/10,
       col.regions = pal(100),xlab='Longitude',ylab='Latitude',
        main=('Regression of SSTA With ENSO')) +
        layer(sp.lines(world.coast)) +
        layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

levelplot(reg~x*y, data=grid , at=c(-5:5)/10,
       col.regions = pal(100),xlab='Longitude',ylab='Latitude',
        xlim=c(-80,60),ylim = c(-60,50), 
       main=('Regression of SSTA With ENSO')) +
        layer(sp.lines(world.coast)) +
        layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

##Comments ENSO shows a direct correlation and regression on the African border. While the significance of the patterns observed may not compare well to that of the AMO previouslt, there is still a correlation to be seen from the maps, at South America and the east coast of Africa, and the regression map reiterating this on a smaller scale. The link between SST anomalies and tropical ENSO events in the north Pacific can ve successfully reproduced by forcing the model atmosphere with tropical Pacific SST variations and allowing atmospheric perturbations (lau; 1997). This signifies that if it can be modelled showing the link between ENSO events, influencing SST due to ocean circulation coupling, ocean circulation does have an effect on these patterns.

PDO

ncpdo<-nc_open(file.path(path,'ipdo_ersst.nc'))

timepdo<-ncvar_get(ncpdo, 'time') 
tunitspdo<-ncatt_get(ncpdo,"time",attname="units")
tustrpdo<-strsplit(tunitspdo$value, " ")
datepdo<-as.character(as.Date(timepdo*365.25/12,origin=unlist(tustrpdo)[3]))
head(datepdo)

## [1] "1880-01-15" "1880-02-14" "1880-03-15" "1880-04-15" "1880-05-15"
## [6] "1880-06-15"

tail(datepdo)

## [1] "2019-07-17" "2019-08-16" "2019-09-16" "2019-10-16" "2019-11-16"
## [6] "2019-12-16"

get the PDO data

sstpdo <- ncvar_get(ncpdo, 'index')

replace missing values with NA

fillvaluepdo <- ncatt_get(ncpdo,"index","_FillValue")
sstpdo[sstpdo==fillvaluepdo$value] <- NA
missvaluepdo <- ncatt_get(ncpdo,"index","missing_value")
sstpdo[sstpdo==missvaluepdo$value] <- NA
sstpdo[sstpdo==3000000000000000000000000000000000] <- NA
summary(sstpdo)  

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -3.57816 -0.87640  0.01883  0.00000  0.81825  3.45459        9

getting the annual mean of pdo

pdoyear <- format(as.Date(datepdo, format ="%Y-%m-%d"),"%Y")
pdoannmean <- aggregate(sstpdo,by=list(pdoyear), FUN= mean, na.rm= TRUE)
head(pdoannmean)

##   Group.1          x
## 1    1880 -1.0941062
## 2    1881 -0.1329750
## 3    1882 -1.2800012
## 4    1883 -1.0194251
## 5    1884  0.4898186
## 6    1885  0.8920912

ensoannmean - 1880-2017

colnames(pdoannmean) <- c("Year", "PDO")    

chop to the same size

dim(asst)

## [1] 360 180 138

yrs[1]

## [1] "1880"

yrs[length(yrs)]

## [1] "2017"

1870-2017

pdoannmean$Year

##   [1] "1880" "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889"
##  [11] "1890" "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899"
##  [21] "1900" "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909"
##  [31] "1910" "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919"
##  [41] "1920" "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929"
##  [51] "1930" "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939"
##  [61] "1940" "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949"
##  [71] "1950" "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959"
##  [81] "1960" "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969"
##  [91] "1970" "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979"
## [101] "1980" "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989"
## [111] "1990" "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999"
## [121] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [131] "2010" "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"

1880-2019, overlapping period is 1880-2017, cut off end of the amv data, this can be done visually

pdoannmean <- pdoannmean[1:138,]

similarly, cutt off the asst data, and yrs

yrspdo  <- yrs[1:138]
asstpdo <- asst[,,1:138]

now the data are the same size, set

sst_tspdo <- pdoannmean$PDO

look at this timeseries:

plot(yrspdo,sst_tspdo,type='l',xlab='Year',ylab='SST Anomaly',main='PDO')

lets calculate the correlation in a loop also

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(asstpdo[i,j,][!is.na(asstpdo[i,j,])])>2){
      c.matrix[i,j] <- cor(asstpdo[i,j,], sst_tspdo)
      p.vals <- cor.test(asstpdo[i,j,], sst_tspdo)
      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 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05), 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 , # xlim=c(-120,10),ylim=c(0,80),  # at=c(-1:1),
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=paste0('Correlation of SSTA with PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

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',
          xlim = c(-200,-50),ylim = c(-70,70),
          main=paste0('Correlation of SSTA with PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

using LM. Linear regressionm, lets calculate the correlation in a loop also

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(asstpdo[i,j,][!is.na(asstpdo[i,j,])])>2){
      r.lm <- lm(asstpdo[i,j,]~sst_tspdo)
      c.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      t.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

add to the same data frame as earlier for plotting

grid$reg <- as.vector(c.matrix)
grid$sig <- as.vector(t.matrix)

Assuming a 99% significance level, subset the grid data where pval.nao < 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(-5:5)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=('Regression of SSTA With PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

levelplot(reg~x*y, data=grid , at=c(-5:5)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',       
          xlim = c(-200,-50),ylim = c(-70,70),
          main=('Regression of SSTA With PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

##Comments Much like the AMO, PDO show a high significance level of the patterns shown in the correlation and regression of SST to PDO. There was a high correlation for PDO with SST, certainly off the coast of North and South America, andf then veering off towards the Asian side of the Pacific.(McPhaden and Zhang; 2002) discuss the slowdonw of the meridional circulation in the upper Pacific slowing down since the 1970s. It is causing an upwelling of about 25%. This upwelling has is associated with a rise in SST’s of about 0.8°C. This shows the clear role ocean circulation may have played in producing the patterns shown in the correlationa dn regression maps made of SST with PDO.

Sea Level Pressure

ncslp<-nc_open(file.path(path,'slp.mon.mean.nc'))

lat1<-ncvar_get(ncslp, 'lat') # latitude
lon1 <- ncvar_get(ncslp, 'lon')
lon1<- ifelse( lon1<=180, lon1, lon1-360)

Load time

timeslp<-ncvar_get(ncslp, 'time')
tunitsslp<-ncatt_get(ncslp,"time",attname="units")
tustrslp<-strsplit(tunitsslp$value, " ")     

Checking that happy with data time origin

dateslp<-as.character(as.Date(timeslp/24-2,origin=unlist(tustrslp)[3]))

get the data for sea level pressure (slp)

slp <- ncvar_get(ncslp, 'slp')
# load slp.mon.mean.nc
latslp <- ncvar_get(ncslp, 'lat')
lonslp <- ifelse(lon1<=180,lon1, lon1-360)
timeslp <- ncvar_get(ncslp, 'time')

replace missing values with NA, many datasets have different fillvalues

fillvalue <- ncatt_get(ncslp,"slp","_FillValue")
slp[slp==fillvalue$value] <- NA
missvalue <- ncatt_get(ncslp,"slp","missing_value")
slp[slp==missvalue$value] <- NA
slp[slp==-9.96921e+36] <- NA

use annual averages with aggregate and using colMeans and rowMeans

yearslp <- format(as.Date(dateslp, format="%Y-%m-%d"),"%Y")
gmeanslp <- colMeans(slp, na.rm = TRUE, dims=2)
annmeanslp <- aggregate(gmeanslp,by=list(yearslp),FUN=mean,na.rm=TRUE)

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

two lines of code to get nicer colours for the eventual plots

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

for this we can use ‘levelplot’ in the lattice package

levelplot(avslp,col.regions = pal(100));

expand.grid will get lats and lons added to a data frame expanded on each others size

grid <- expand.grid(x=lon1, y=lat1)
grid$avslp <- as.vector(avslp)

levelplot(avslp~x*y,grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SLP'
) + 
  layer(sp.lines(world.coast))

not interested in short term variability, make annual averages

yrsslp <- annmeanslp$Group.1
nyrslp <- length(yrsslp)

aslp <- array(slp, c(dim(lon1),dim(lat1),nyrslp))  

use a for loop to load the annual data into asst

for (k in 1:nyrslp) {
  aslp[,,k] <- rowMeans(slp[,,yearslp==yrsslp[k]],na.rm=FALSE,dims=1) # annual averages from monthly data
}

note that all we’ve done is create annual averages, so the average looks exactly the same, add to the same data frame as earlier for plotting

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

plot

levelplot(avslp~x*y, data=grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') + 
  layer(sp.lines(world.coast)) 

###remove the global mean from each year ###lets do a traditional loop

gmeanslp <- colMeans(aslp, na.rm = TRUE, dims=2)
for (k in 1:nyrslp){
  aslp[,,k]<-aslp[,,k]-matrix(gmeanslp[k],length(lon1),length(lat1))
}

lon0 <- -10.5 #
lat0 <- 51.5 #
slp_ts<-aslp[which(lon1==lon0),which(lat1==lat0),]
slp_ts <- annmeanslp$x

look at this timeseries:

plot(yrsslp,slp_ts,type='l',xlab='Year',ylab='SLP Anomaly',main=paste0('SLPA at Long=', lon0, ',Lat=', lat0))

PDO

ncpdo<-nc_open(file.path(path,'ipdo_ersst.nc'))

timepdo<-ncvar_get(ncpdo, 'time') 
tunitspdo<-ncatt_get(ncpdo,"time",attname="units")
tustrpdo<-strsplit(tunitspdo$value, " ")
datepdoslp<-as.character(as.Date(timepdo*365.25/12,origin=unlist(tustrpdo)[3]))
head(datepdoslp)

## [1] "1880-01-15" "1880-02-14" "1880-03-15" "1880-04-15" "1880-05-15"
## [6] "1880-06-15"

tail(datepdoslp)

## [1] "2019-07-17" "2019-08-16" "2019-09-16" "2019-10-16" "2019-11-16"
## [6] "2019-12-16"

get the PDO data

slppdo <- ncvar_get(ncpdo, 'index')

replace missing values with NA

fillvaluepdo <- ncatt_get(ncpdo,"index","_FillValue")
slppdo[slppdo==fillvaluepdo$value] <- NA
missvaluepdo <- ncatt_get(ncpdo,"index","missing_value")
slppdo[slppdo==missvaluepdo$value] <- NA
slppdo[slppdo==3000000000000000000000000000000000] <- NA
summary(slppdo)  

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -3.57816 -0.87640  0.01883  0.00000  0.81825  3.45459        9

getting the annual mean of pdo

pdoyearslp <- format(as.Date(datepdoslp, format ="%Y-%m-%d"),"%Y")
pdoannmeanslp <- aggregate(slppdo,by=list(pdoyearslp), FUN= mean, na.rm= TRUE)
head(pdoannmeanslp)

##   Group.1          x
## 1    1880 -1.0941062
## 2    1881 -0.1329750
## 3    1882 -1.2800012
## 4    1883 -1.0194251
## 5    1884  0.4898186
## 6    1885  0.8920912

ensoannmean - 1880-2017

colnames(pdoannmeanslp) <- c("Year", "PDO")    

chop to the same size

dim(aslp)

## [1] 144  73  72

yrsslp[1]

## [1] "1947"

yrsslp[length(yrsslp)]

## [1] "2018"

1870-2017

pdoannmeanslp$Year

##   [1] "1880" "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889"
##  [11] "1890" "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899"
##  [21] "1900" "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909"
##  [31] "1910" "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919"
##  [41] "1920" "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929"
##  [51] "1930" "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939"
##  [61] "1940" "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949"
##  [71] "1950" "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959"
##  [81] "1960" "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969"
##  [91] "1970" "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979"
## [101] "1980" "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989"
## [111] "1990" "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999"
## [121] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [131] "2010" "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"

1880-2019, overlapping period is 1880-2017, cut off end of the amv data, this can be done visually

pdoannmeanslp1 <- pdoannmeanslp[67:138,]

similarly, cutt off the aslp data, and yrs

yrspdoslp  <- yrsslp[1:72]
aslppdo <- aslp[,,1:72]

now the data are the same size, set

slp_tspdo <- pdoannmeanslp1$PDO

look at this timeseries:

plot(yrspdoslp,slp_tspdo,type='l',xlab='Year',ylab='SST Anomaly',main='PDO SLP Time Series')

lets calculate the correlation in a loop also

c.matrix <- matrix(NA,length(lon1),length(lat1))
t.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslppdo[i,j,][!is.na(aslppdo[i,j,])])>2){
      c.matrix[i,j] <- cor(aslppdo[i,j,], slp_tspdo)
      p.vals <- cor.test(aslppdo[i,j,], slp_tspdo)
      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 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05), 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 , # 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 PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

###using LM. Linear regression ###lets calculate the correlation in a loop also

m.matrix <- matrix(NA,length(lon1),length(lat1))
n.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslppdo[i,j,][!is.na(aslppdo[i,j,])])>2){
      r.lm <- lm(aslppdo[i,j,]~slp_tspdo)
      m.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      n.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

add to the same data frame as earlier for plotting

grid$reg <- as.vector(m.matrix)
grid$sig <- as.vector(n.matrix)

Assuming a 99% significance level, subset the grid data where pval.nao < 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(-2:2)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=('Regression of SLP With PDO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

##Comment I am not able to comment on the maps made bny any of the incdes due to a small error that i cannot find or seem to fix when coding the maps. It has restricted the reading of the maps, thus patterns cannot be seen or discussed unfortunately. It is not for a lack of tryingm, as I have been working on this a long time.

ENSO

ncenso<-nc_open(file.path(path,'ihadisst1_nino12a.nc'))

timeenso<-ncvar_get(ncenso, 'time') 
tunitsenso<-ncatt_get(ncenso,"time",attname="units")
tustrenso<-strsplit(tunitsenso$value, " ")
dateensoslp<-as.character(as.Date(timeenso*365.25/12,origin=unlist(tustrenso)[3]))
head(dateensoslp)

## [1] "1870-01-15" "1870-02-14" "1870-03-16" "1870-04-16" "1870-05-16"
## [6] "1870-06-16"

tail(dateensoslp)

## [1] "2019-07-17" "2019-08-17" "2019-09-16" "2019-10-17" "2019-11-16"
## [6] "2019-12-17"

get the ENSO data

slpenso <- ncvar_get(ncenso, 'Nino12')

replace missing values with NA

fillvalueenso <- ncatt_get(ncenso,"Nino12","_FillValue")
slpenso[slpenso==fillvalueenso$value] <- NA
missvalueenso <- ncatt_get(ncenso,"Nino12","missing_value")
slpenso[slpenso==missvalueenso$value] <- NA
slpenso[slpenso==3000000000000000000000000000000000] <- NA
summary(slpenso)  

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -2.5295 -0.8158 -0.3026 -0.1794  0.2769  4.3790      11

getting the annual mean of enso

ensoyearslp <- format(as.Date(dateensoslp, format ="%Y-%m-%d"),"%Y")
ensoannmeanslp <- aggregate(slpenso,by=list(ensoyearslp), FUN= mean, na.rm= TRUE)
head(ensoannmeanslp)

##   Group.1          x
## 1    1870 -1.0840045
## 2    1871 -0.5825608
## 3    1872 -0.9003969
## 4    1873 -0.9435940
## 5    1874 -0.8821288
## 6    1875 -1.0859472

ensoannmean - 1880-2017

colnames(ensoannmeanslp) <- c("Year", "ENSO")    

chop to the same size

dim(aslp)

## [1] 144  73  72

yrsslp[1]

## [1] "1947"

yrsslp[length(yrsslp)]

## [1] "2018"

1870-2017

ensoannmeanslp$Year

##   [1] "1870" "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879"
##  [11] "1880" "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889"
##  [21] "1890" "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899"
##  [31] "1900" "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909"
##  [41] "1910" "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919"
##  [51] "1920" "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929"
##  [61] "1930" "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939"
##  [71] "1940" "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949"
##  [81] "1950" "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959"
##  [91] "1960" "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969"
## [101] "1970" "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979"
## [111] "1980" "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989"
## [121] "1990" "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999"
## [131] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [141] "2010" "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"

1880-2019, overlapping period is 1880-2017, cut off end of the amv data, this can be done visually

ensoannmeanslp1 <- ensoannmeanslp[78:138,]

similarly, cutt off the asst data, and yrs

yrsensoslp  <- yrsslp[1:61]
aslpenso <- aslp[,,1:61]

now the data are the same size, set

slp_tsenso <- ensoannmeanslp1$ENSO

look at this timeseries:

plot(yrsensoslp,slp_tsenso,type='l',xlab='Year',ylab='SST Anomaly',main='ENSO SLP Time Series')

lets calculate the correlation in a loop also

c.matrix <- matrix(NA,length(lon1),length(lat1))
t.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslpenso[i,j,][!is.na(aslpenso[i,j,])])>2){
      c.matrix[i,j] <- cor(aslpenso[i,j,], slp_tsenso)
      p.vals <- cor.test(aslpenso[i,j,], slp_tsenso)
      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 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05), 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 , # 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.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

### using LM. Linear regression ###lets calculate the correlation in a loop also

c.matrix <- matrix(NA,length(lon1),length(lat1))
t.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslpenso[i,j,][!is.na(aslpenso[i,j,])])>2){
      r.lm <- lm(aslpenso[i,j,]~slp_tsenso)
      c.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      t.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

add to the same data frame as earlier for plotting

grid$reg <- as.vector(m.matrix)
grid$sig <- as.vector(n.matrix)

Assuming a 99% significance level, subset the grid data where pval.nao < 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(-2:2)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=('Regression of SLP With ENSO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

AMO

ncamo<-nc_open(file.path(path,'iamo_ersst.nc'))

timeamo<-ncvar_get(ncamo, 'time') 
tunitsamo<-ncatt_get(ncamo,"time",attname="units")
tustramo<-strsplit(tunitsamo$value, " ")
dateamoslp<-as.character(as.Date(timeamo*365.25/12,origin=unlist(tustramo)[3]))
head(dateamoslp)

## [1] "1880-01-15" "1880-02-14" "1880-03-15" "1880-04-15" "1880-05-15"
## [6] "1880-06-15"

tail(dateamoslp)

## [1] "2019-07-17" "2019-08-16" "2019-09-16" "2019-10-16" "2019-11-16"
## [6] "2019-12-16"

get the AMO data

slpamo <- ncvar_get(ncamo, 'AMO')

replace missing values with NA

fillvalueamo <- ncatt_get(ncamo,"AMO","_FillValue")
slpamo[slpamo==fillvalueamo$value] <- NA
missvalueamo <- ncatt_get(ncamo,"AMO","missing_value")
slpamo[slpamo==missvalueamo$value] <- NA
slpamo[slpamo==3e+33] <- NA
summary(slpamo)  

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## -0.733522 -0.169803 -0.004612  0.000000  0.168052  0.851506         9

getting the annual mean of AMO

amoyearslp <- format(as.Date(dateamoslp, format ="%Y-%m-%d"),"%Y")
amoannmeanslp <- aggregate(slpamo,by=list(amoyearslp), FUN= mean, na.rm= TRUE)
head(amoannmeanslp)

##   Group.1           x
## 1    1880  0.26576526
## 2    1881 -0.02418830
## 3    1882  0.06306856
## 4    1883  0.01907216
## 5    1884 -0.08529102
## 6    1885  0.04853270

amooannmean - 1880-2017

colnames(amoannmeanslp) <- c("Year", "AMO")    

chop to the same size

dim(aslp)

## [1] 144  73  72

yrsslp[1]

## [1] "1947"

yrsslp[length(yrsslp)]

## [1] "2018"

1870-2017

amoannmeanslp$Year

##   [1] "1880" "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889"
##  [11] "1890" "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899"
##  [21] "1900" "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909"
##  [31] "1910" "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919"
##  [41] "1920" "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929"
##  [51] "1930" "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939"
##  [61] "1940" "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949"
##  [71] "1950" "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959"
##  [81] "1960" "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969"
##  [91] "1970" "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979"
## [101] "1980" "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989"
## [111] "1990" "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999"
## [121] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [131] "2010" "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"

1947-2019, overlapping period is 1947-2017, cut off end of the amv data, this can be done visually

amoannmeanslp1 <- amoannmeanslp[68:138,]

similarly, cutt off the asst data, and yrs

yrsamoslp  <- yrsslp[1:71]
aslpamo <- aslp[,,1:71]

now the data are the same size, set

slp_tsamo <- amoannmeanslp1$AMO

look at this timeseries:

plot(yrsamoslp,slp_tsamo,type='l',xlab='Year',ylab='SLP Anomaly',main='AMO SLP Time Series')

lets calculate the correlation in a loop also

c.matrix <- matrix(NA,length(lon1),length(lat1))
t.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslpamo[i,j,][!is.na(aslpamo[i,j,])])>2){
      c.matrix[i,j] <- cor(aslpamo[i,j,], slp_tsamo)
      p.vals <- cor.test(aslpamo[i,j,], slp_tsamo)
      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 5% significance level, subset the grid data where pval < 0.01. (if this were 95% CI, then pval<0.05), 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 , # 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 AMO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

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',
          xlim=c(-120,10),ylim=c(0,80),  
          main=paste0('Correlation of SLP with AMO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

###using LM. Linear regression

c.matrix <- matrix(NA,length(lon1),length(lat1))
t.matrix <- matrix(NA,length(lon1),length(lat1))

for (i in 1:dim(lon1)) {
  for (j in 1:dim(lat1)) {
    if (length(aslpamo[i,j,][!is.na(aslpamo[i,j,])])>2){
      r.lm <- lm(aslpamo[i,j,]~slp_tsamo)
      c.matrix[i,j] <- r.lm$coefficients[2]
      smm<-summary(r.lm)
      t.matrix[i, j] <- smm$coefficients[8]
    }
  }
}

add to the same data frame as earlier for plotting

grid$reg <- as.vector(c.matrix)
grid$sig <- as.vector(t.matrix)

Assuming a 99% significance level, subset the grid data where pval.nao < 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(-2:2)/10,
          col.regions = pal(100),xlab='Longitude',ylab='Latitude',
          main=('Regression of SLP With AMO')) +
  layer(sp.lines(world.coast)) +
  layer(sp.points(sig, pch = 20, cex = 0.005, col = "black"))

Bibliography

Gastineau, G. and Frankignoul, C., 2015. Influence of the North Atlantic SST variability on the atmospheric circulation during the twentieth century. Journal of Climate, 28(4), pp.1396-1416.

Lau, N.C., 1997. Interactions between global SST anomalies and the midlatitude atmospheric circulation. Bulletin of the American Meteorological Society, 78(1), pp.21-34.

McPhaden, M.J. and Zhang, D., 2002. Slowdown of the meridional overturning circulation in the upper Pacific Ocean. Nature, 415(6872), p.603.