Assignment II
Michael Byrnes
26 April 2019
Preliminaries
Clear The Working station
graphics.off()
rm(list=ls())Load packages necessary for the assessment
library(sp)
library(raster)
library(ncdf4)
library(rgdal)## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/micha/OneDrive/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/micha/OneDrive/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)
library(ncdf4.helpers)Set the directory
setwd("C:/Users/micha/Downloads/WorkshopIII")Set the path to the data
## Set the path to the data
path <- file.path(getwd(),"data");Load in the data
load("world.coast.RData")Question I - Load in Sea Surface Temp (SST) and Sea Level Pressure (SLP) data # Firstly, we are going to work with Sea Surface Temp (SST) using the Had data, previously used in workshop III, using
Load the NetCDF file “HadISST_sst.nc” and convert the date format to something that R recognises.
nc <- nc_open(file.path(getwd(),"HadISST_sst.nc")); # open netcdf file
lat <- ncvar_get(nc,"latitude") # load the latitude
lon <- ncvar_get(nc,"longitude") # load the longitude
ts <- ncvar_get(nc,"time"); # load the time
tunits<-ncatt_get(nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
date<-as.character(as.Date(ts,origin=unlist(tustr)[3]))
sst <- ncvar_get(nc,"sst")Many dataset have different ‘fillvalue’s’, and we want to replace the missing values with NA.
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# We are interested in trends of global variability, therefore we need to remove the global mean at each step. But first, we must investigate the global mean. Use annual averages with aggregate and using colMeans and rowMeans. Firstly format the columns.
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)# Comment out the following lines.
avsst = rowMeans(sst,na.rm=FALSE,dims=2)## The next step is to use some code to enrich the colours for the plots that will be created.
colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)
#We can use the 'levelplot' function, which is found in the lattice package.This will give us a colourbar, making the plot look nice. It is important to add lat & long also, as levelplot works well with data frames.
levelplot(avsst,col.regions = pal(100));
#By expanding the grid, latitudes and longitudes will be added to a data frame expanded on each others size.
grid <- expand.grid(x=lon, y=lat)
# The avsst that we want to plot can be added as a vector, which is essentially expanding what we had.
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))
# The next step is to make annual averages as we are not interested in short term variability.
yrs <- annmean$Group.1
nyr <- length(yrs)
asst <- array(NA,c(dim(lon),dim(lat),nyr))# use a for loop to load the annual data into asst
# note: for loops are a bit old school and I'm sure there is a neat vectorised way to do this-if you know it, let me know
for (k in 1:nyr) {
asst[,,k] <- rowMeans(sst[,,year==yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}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))
}# note: keep the .5 at the end to find a matching lat/lon
lon0 <- -10.5 #
lat0 <- 51.5 #
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]# Check out the timeseries for SST
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon0, ',Lat=', lat0))
# lets calculate the correlation in a loop also
cmatrix <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
cmatrix[i,j] <- cor(asst[i,j,], sst_ts)
}
}
# add to the same data frame as earlier for plotting
grid$corr <- as.vector(cmatrix)#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 Long=', lon0, ',Lat=', lat0)) +
layer(sp.lines(world.coast)) 
Now we are going to examine the Sea Level Pressure (SLP) Data. This was carried out in a very similar fashion to the SST data above.
slp.ncdf <- nc_open(file.path(getwd(),"slp.mon.mean.nc")) # open netcdf file
lat.slp <- ncvar_get(slp.ncdf,"lat") # load the latitude
lon.slp <- ncvar_get(slp.ncdf,"lon") # load the longitude
ts.slp<- ncvar_get(slp.ncdf,"time"); # load the time
change_ts <- ts.slp/24
tunits.slp <-ncatt_get(slp.ncdf,"time",attname="units")
tustr.slp <-strsplit(tunits.slp$value, " ")# check data time origin
slp.date <-as.character(as.Date(change_ts,origin=unlist(tustr.slp)[3]))
slp <- ncvar_get(slp.ncdf,"slp")
slp.lon <- ifelse(lon.slp<=180,lon.slp,lon.slp-360)
# Missing Values Replaced With NA
fillvalue <- ncatt_get(slp.ncdf,"slp","_FillValue")
slp[slp==fillvalue$value] <- NA
missvalue <- ncatt_get(slp.ncdf,"slp","missing_value")
slp[slp==missvalue$value] <- NA
slp[slp==-1000] <- NA#Use the data loaded
slp.year <- format(as.Date(slp.date, format="%Y-%m-%d"),"%Y")
slp.gmean <- colMeans(slp, na.rm = TRUE, dims=2)
slpanmean <- aggregate(slp.gmean,by=list(slp.year),FUN=mean,na.rm=TRUE)
slpav = rowMeans(slp,na.rm=FALSE,dims=2)
colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)
levelplot(slpav,col.regions = pal(100))
grid.2 <- expand.grid(x=slp.lon, y=lat.slp)
grid.2$slpav <- as.vector(slpav)
levelplot(slpav~x*y,grid.2,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Average SLP'
) +
layer(sp.lines(world.coast))
slp.years <- slpanmean$Group.1
slp.nyr <- length(slp.years)
slpa <- array(NA,c(dim(slp.lon),dim(lat.slp),slp.nyr)) # asst possesses the annual sea surface temp.
# Use a loop to increase the speed of the computation
for (k in 1:slp.nyr) {
slpa[,,k] <- rowMeans(slp[,,slp.year==slp.years[k]],na.rm=FALSE,dims=2)
}grid.2$an_slpav <- as.vector(rowMeans(slpa,na.rm=FALSE,dims=2))
#Plot the Annual averages of SLP
levelplot(an_slpav~x*y, data=grid.2,col.regions = pal(100),xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') +
layer(sp.lines(world.coast))
#Take out the Gmean for the years
slp_Gmean <- colMeans(slpa, na.rm = TRUE, dims=2)
for (k in 1:slp.nyr){
slpa[,,k]<-slpa[,,k]-matrix(slp_Gmean[k],length(slp.lon),length(lat.slp))
}#keep the .5 at the end to find a matching lat/lon
lon0 <- -10 #
lat0 <- 50 #
ts.slp<-slpa[which(slp.lon==lon0),which(lat.slp==lat0),]# Plot a timeseries of the data
plot(slp.years,ts.slp,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SST at Long=', lon0, ',Lat=', lat0))
#Calculate the correlation in a loop
cmatrix <- matrix(NA,dim(slp.lon),dim(lat.slp))
for (i in 1:dim(slp.lon)) {
for (j in 1:dim(lat.slp)) {
cmatrix[i,j] <- cor(slpa[i,j,], ts.slp)
}
}# aAdd to the earlier data frame, to help plotting
grid.2$corr <- as.vector(cmatrix)#Plot a correlatiuon of the SST
levelplot(corr~x*y, data=grid.2 , 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 Long=', lon0, ',Lat=', lat0)) +
layer(sp.lines(world.coast)) 
Now SST and SLP have been loaded
Question II - Calculate the Correlation and Regression Patterns associated with the climate indices for both SST and SLP
Now we are going to examine the AMO data, which was one of the climate indices given to examine, found at http://climexp.knmi.nl/selectindex.cgi?id=ididsomeone@somewhere
The Atlantic Ocean has multidecadal warm and cool modes with periods of approximately 30 years. This is similar to the PDO. During warm phases, the AMO is warm in the tropical N. Atlantic and far N. Atlantic and reasonably cool in the central area. During cool phases, the tropical are and Far North Atlantic are cool and the central ocean is warm. (Easterbrook, 2016)
amo <- nc_open(file.path(getwd(),"iamo_ersst.nc")); # open netcdf file
amo_ts <- ncvar_get(amo,"time"); # load the time
amo_units<-ncatt_get(amo,"time",attname="units")
amo_tustr<-strsplit(tunits$value, " ")
amo_date<-as.character(as.Date(amo_ts,origin=unlist(tustr)[3]))
amo_sst <- ncvar_get(amo,"AMO")date<-as.character(as.Date(amo_ts,origin=unlist(tustr)[3]))
AMO <- ncvar_get(amo,"AMO")### Many dataset have different 'fillvalue's', and we want to replace the missing values with NA.
fillvalue <- ncatt_get(amo,"AMO","_FillValue")
AMO[AMO==fillvalue$value] <- NA
missvalue <- ncatt_get(amo,"AMO","missing_value")
AMO[AMO==missvalue$value] <- NA
AMO[AMO==-1000] <- NAAMO_time <- seq.Date(as.Date(unlist(amo_tustr)[3]),length.out = length(amo_ts),by = "month")We want to change the date format of this data
AMO_year <- format(as.Date(AMO_time, format="%Y-%m-%d"),"%Y")
annmean <- aggregate(AMO,by=list(AMO_year),FUN=mean,na.rm=TRUE)# We want to extract mean values and mean years
amo_years <- annmean$Group.1
amo_ts <- annmean$x
## Make a plot for AMO ts
plot(amo_years,amo_ts,xlab="Year",ylab="AMO [C]",col="indianred", type="l",main="AMO ts")
abline(0.0,0.0,type="l", lty=2)## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...):
## graphical parameter "type" is obsoletelegend("bottomright", legend=c("AMO ts"), col=c("indianred"),lty=1:1, cex=0.7)
#The data ranges from 1870-2017. The next step is to extract the dates that happen to overlap for Sea Surface Temp
amo_ts <- annmean$x[21:168]We run a loop, just like we did previously.
cmatrix <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
cmatrix[i,j] <- cor(asst[i,j,], amo_ts)
}
}
grid$corr <- as.vector(cmatrix)We look to plot the correlation of SST and AMO
‘Note’ -> It would plot, but unfortunately it would not knit on RMarkdown. The code is i used can be seen just below
Next, we look at the Sea Level Pressure for AMO Data nad how they correlate
amo.slp <- annmean$x[99:169]
for (i in 1:dim(slp.lon)) {
for (j in 1:dim(lat.slp)) {
cmatrix[i,j] <- cor(slpa[i,j,], amo.slp)
}
}
#plot the correlation of Sea Level Pressure with AMO
levelplot(corr~x*y, data=grid.2 ,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=("SLPA Correlated with AMO")) +
layer(sp.lines(world.coast)) 
The next step is to carry out linear regression of SLP and AMO
The correlation can be run in a loop
matrix_r <- matrix(NA,length(slp.lon),length(lat.slp))
matrix_s <- matrix(NA,length(slp.lon),length(lat.slp))
for (i in 1:dim(slp.lon)) {
for (j in 1:dim(lat.slp)) {
if (length(slpa[i,j,][!is.na(slpa[i,j,])])>2){
lm.r <- lm(slpa[i,j,]~amo.slp)
matrix_r[i,j] <- lm.r$coefficients[2]
smm<-summary(lm.r)
matrix_s[i, j] <- smm$coefficients[8]
}
}
}
# add to the same data frame as earlier for plotting
grid.2$reg <- as.vector(matrix_r)
grid.2$sig <- as.vector(matrix_s)#plot the data of the Sea Level Pressure associated with AMO data
levelplot(reg~x*y, data=grid.2 , at=c(-50:70)/10,
col.regions = pal(200),xlab='Longitude',ylab='Latitude',
main=('SLPA with AMO Regression')) +
layer(sp.lines(world.coast)) 
# layer(sp.points(sig, pch = 20, cex = 0.1, col = "black"))The patters between the correlation and regression of the climate Indice ‘Atlantic Multidecadal Oscillation’ (AMO) appear to be very similar in the context of both SST and SLP. With regards to Sea Level Pressure (SLP) the correlation of SLP the map shows a much more dramtic change, particularly in the mid-latitudes. Although, it is unfortunate that in the regression model there is missing information which is crucial, from a great portion of Asia and North Africa.
Question III - Correlation and Regression Patterns with Pacific Decadal Oscillation (PDO) Data
Load in the pdo data, located at http://climexp.knmi.nl/selectindex.cgi?id=ididsomeone@somewhere
The Pacific Decadal Oscillation (PDO) refers to the cyclical variations in sea surface temps in the Pacific Ocean. The PDO index is defined as the leading principal component of North Pacific monthly sea surface temperature variability. (Easterbrook,2016). Each warm PDO phase lasted about 25-30 years, and vice versa with a cool phase.
pdo_nc <- nc_open(file.path(getwd(),"ipdo.nc")) # open netcdf file
pdo_ts <- ncvar_get(pdo_nc,"time"); # load the time
pdo_tunits <-ncatt_get(pdo_nc,"time",attname="units")
pdo_tustr <-strsplit(pdo_tunits$value, " ")
date<-as.character(as.Date(pdo_ts,origin=unlist(tustr)[3]))
PDO <- ncvar_get(pdo_nc,"PDO")# replace missing values with NA
fillvalue <- ncatt_get(pdo_nc,"PDO","_FillValue")
PDO[PDO==fillvalue$value] <- NA
missvalue <- ncatt_get(pdo_nc,"PDO","missing_value")
PDO[PDO==missvalue$value] <- NA
PDO[PDO==-1000] <- NA
#Extract the Dates
PDO_month <- seq.Date(as.Date(unlist(pdo_tustr)[3]),length.out = length(pdo_ts),by="month")#Change the format of how the data is dated
pdo_yr <- format(as.Date(PDO_month, format="%Y-%m-%d"),"%Y")
pdo_annmean <- aggregate(PDO,by=list(pdo_yr),FUN=mean,na.rm=TRUE)#As seen before, we are going to extract the mean values and years of the PDO data
pdo_years <- pdo_annmean$Group.1
pdo_ts <- pdo_annmean$x#Make a plot to look at the PDO ts
plot(pdo_years,pdo_ts,xlab="Year",ylab="AMO [C]",col="indianred", type="l",main="PDO ts")
abline(0.0,0.0,type="l", lty=2)## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...):
## graphical parameter "type" is obsoletelegend("bottomright", legend=c("PDO hadsst"), col=c("blue"),lty=1:1, cex=0.7)
#The dates run from 1900-2017. Extract the dates that happen to overlap for sst
years_sst <- annmean$x[31:148]
pdoyr <- length(years_sst)
pdo.asst <- array(NA,c(dim(lon),dim(lat),pdoyr))
for (k in 1:pdoyr) {
pdo.asst[,,k] <- rowMeans(sst[,,year[361:1776]==yrs[31:148][k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}grid$an_avsst <- as.vector(rowMeans(pdo.asst,na.rm=FALSE,dims=2))plot the annual avergae of sst for PDO
levelplot(an_avsst~x*y, data=grid,col.regions = pal(100),xlab='Longitude',ylab='Latitude',main='SST Annual Avg') +
layer(sp.lines(world.coast)) 
#Remove the Gmean for all of the years from
#Use a loop as seen before
Pgmean <- colMeans(pdo.asst, na.rm = TRUE, dims=2)
for (k in 1:pdoyr){
pdo.asst[,,k]<-pdo.asst[,,k]-matrix(Pgmean[k],length(lon),length(lat))
}matrix_pdo <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
matrix_pdo[i,j] <- cor(pdo.asst[i,j,], pdo_ts)
}
}
# add to the same data frame as earlier for plotting
grid$corr <- as.vector(matrix_pdo)
#plot
levelplot(corr~x*y, data=grid ,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=('SSTA Correlated with PDO')) +
layer(sp.lines(world.coast)) 
# The next step is to carry out linear regression of Sea Surface Temperature with PDO
matrix_Regpdo1 <- matrix(NA,length(lon),length(lat))
matrix_Regpdo2 <- matrix(NA,length(lon),length(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
if (length(pdo.asst[i,j,][!is.na(pdo.asst[i,j,])])>2){
d.lm <- lm(pdo.asst[i,j,]~pdo_ts)
matrix_Regpdo1[i,j] <- d.lm$coefficients[2]
smm<-summary(d.lm)
matrix_Regpdo2[i, j] <- smm$coefficients[8]
}
}
}
# add to the same data frame, which will be necessary for plotting
grid$reg <- as.vector(matrix_Regpdo1)
grid$sig <- as.vector(matrix_Regpdo2)
#plot The regression of SST with PDO
levelplot(reg~x*y, data=grid , at=c(-5:5)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=('SST Regression with PDO')) +
layer(sp.lines(world.coast)) 
## The Regression and Correlation of SST appear to be very similar from both of the plots above.Às expected, the correlation plot looks to have more dramatic reading to it with its richer colours of red.
We must examine Surface Level Pressure of PDO
Firstly, we will try and examine the correlation between SLP and PDO
pdo_slp <- pdo_annmean$x[49:118]
aslp1 <- slpa[,,1:70]
cmatrix <- matrix(NA,dim(slp.lon),dim(lat.slp))
for (i in 1:dim(slp.lon)) {
for (j in 1:dim(lat.slp)) {
cmatrix[i,j] <- cor(aslp1[i,j,], pdo_slp)
}
}
# Add to the same data frame
grid.2$corr <- as.vector(cmatrix)#Plot the SLPA Correlation with PDO data
levelplot(corr~x*y, data=grid.2 ,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=("SLPA Correlation with PDO")) +
layer(sp.lines(world.coast))
Next we want to examine the regression of SLP and PDO
Matrix_PDO <- matrix(NA,length(slp.lon),length(lat.slp))
Matrix_PDO2 <- matrix(NA,length(slp.lon),length(lat.slp))
for (i in 1:dim(slp.lon)) {
for (j in 1:dim(lat.slp)) {
if (length(aslp1[i,j,][!is.na(aslp1[i,j,])])>2){
r.lm <- lm(aslp1[i,j,]~pdo_slp)
Matrix_PDO[i,j] <- r.lm$coefficients[2]
smm<-summary(r.lm)
Matrix_PDO2[i, j] <- smm$coefficients[8]
}
}
}# add to the same data frame, for plotting later
grid.2$reg <- as.vector(Matrix_PDO)
grid.2$sig <- as.vector(Matrix_PDO2)#Now we plot the regression of SLPA with PDO this time
levelplot(reg~x*y, data=grid.2 , at=c(-20:20)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=('SLPA Regression with PDO')) +
layer(sp.lines(world.coast)) 
# layer(sp.points(sig, pch = 20, cex = 0.1, col = "black")) Similar to SST, the correlation of Sea Level pressure projects much greater change, in comparison to using linear regression of SLP on the Pacific Decadal Oscillation.
Question IV - Correlation and regression patterns of the ENSO data
Investigating the climate indice ‘El Nino Southern Oscillation’ (ENSO) data, found at: http://climexp.knmi.nl/selectindex.cgi?id=ididsomeone@somewhere
El Nino Southern Oscillation (ENSO) occurs at intervals that vary from 2-10 years, where sea surface temperatures are unusually high and the tradewinds are unusually weak over the tropical Pacific Ocean
Firstly, load in the ENSO data
enso_nc <- nc_open(file.path(getwd(),"iersst_nino12a.nc")) # open netcdf file
enso_ts <- ncvar_get(enso_nc,"time"); # load the time
enso_tunits<-ncatt_get(enso_nc,"time",attname="units")
enso_tustr <-strsplit(enso_tunits$value, " ")
date<-as.character(as.Date(enso_ts,origin=unlist(enso_tustr)[3]))
ENSO <- ncvar_get(enso_nc,"Nino12")# replace missing values with NA
fillvalue <- ncatt_get(enso_nc,"Nino12","_FillValue")
ENSO[ENSO==fillvalue$value] <- NA
missvalue <- ncatt_get(enso_nc,"Nino12","missing_value")
ENSO[ENSO==missvalue$value] <- NA
ENSO[ENSO==-1000] <- NA# Extract dates by the month
enso.time <- seq.Date(as.Date(unlist(enso_tustr)[3]),length.out = length(enso_ts),by="month")
#Using loaded data
enso.yr <- format(as.Date(enso.time, format="%Y-%m-%d"),"%Y")
enso.annmean <- aggregate(ENSO,by=list(enso.yr),FUN=mean,na.rm=TRUE)
#Extract the Mean Years and the Mean Values
enso.mean.yrs <- enso.annmean$Group.1
enso.times <- enso.annmean$x#Make a plot of the ENSO ts data
plot(enso.mean.yrs,enso.times,xlab="Year",ylab="ENSO [C]",col="indianred", type="l",main="ENSO times")
abline(0.0,0.0,type="l", lty=2)## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...):
## graphical parameter "type" is obsoletelegend("bottomright", legend=c("ENSO ts"), col=c("indianred"),lty=1:1, cex=0.7)
#The dates that overlap with sst between 1870-2017 must be extracted
enso.times <- enso.annmean$x[1:148]#Do a loop with the enso.times data
cmatrix <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
cmatrix[i,j] <- cor(asst[i,j,], enso.times)
}
}# add to the data frame, which will be used for plotting
grid$corr <- as.vector(cmatrix)
# Plot the Correlation of SSTA with ENSO
levelplot(corr~x*y, data=grid ,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=('SSTA Correlated with ENSO')) +
layer(sp.lines(world.coast)) 
Past records are key in trying to create an accurate projection of future changes in the ocean, such as sea level pressure and sea surface temperature. Abundant physical, geologic evidence from the past provides a record of former periods of recurrent global warming and cooling that were far more intense than recent warming and cooling, according to Easterbrook (2016).
References: