GY667_Assignment_2
Natasha Langton
24 April 2019
Assignment 2: Ocean Indice and SST
Climate Indices
Atlantic Multidecadal Oscillation.
Pacific Decadal Oscillation.
El Nino/Souther Oscillation.
Sea Surface Temperature (SLP)
Load the apporpiate libraries for the analysis
library(sp)
library(raster)
library(ncdf4)
library(rgdal)## Warning: package 'rgdal' was built under R version 3.5.3## 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/Natasha/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/Natasha/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)set working directory
setwd("~/College/Msc Climate Change/GY667/Assignment 2")Clear variables and any open plots before starting
graphics.off()
rm(list=ls())
path <- file.path(getwd(),"data")Read coastline using readOGR as a Spatial Data Frame - downloaded from moodle workshop 3 its part of the rgdal package.
world.coast <- readOGR(dsn = file.path(path,"ne_110m_coastline"), layer = "ne_110m_coastline" )## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Natasha\Documents\College\Msc Climate Change\GY667\Assignment 2\data\ne_110m_coastline", layer: "ne_110m_coastline"
## with 134 features
## It has 3 fields
## Integer64 fields read as strings: scalerank#convert to the correct projections
world.coast <- spTransform(world.coast, CRS( "+proj=longlat +datum=WGS84" ) )Need to load in the HADISST_sst.nc, downloaded from Workshop 3
hadisst <- nc_open(file.path(getwd(), "HADISST_sst.nc"))
had.lat <- ncvar_get(hadisst,"latitude")
had.lon <- ncvar_get(hadisst,"longitude")
had.sst <- ncvar_get(hadisst, "sst")
had.date <- ncvar_get(hadisst, "time")Need to get date into the correct format for R to recognise, this will allow for analysis of the data.
had.date <- ncvar_get(hadisst, "time")
had.tunits<-ncatt_get(hadisst,"time",attname="units")
tustr<-strsplit(had.tunits$value, " ")
had.date<-as.character(as.Date(had.date,origin=unlist(tustr)[3]))
head(had.date)## [1] "1870-01-16" "1870-02-14" "1870-03-16" "1870-04-15" "1870-05-16"
## [6] "1870-06-16"Replace missing values with NA. This will allow for a more accurate representation of the data. So that missing data numbers don’t skew the data.
fillvalue <- ncatt_get(hadisst,"sst","_FillValue")
had.sst[had.sst==fillvalue$value] <- NA
missvalue <- ncatt_get(hadisst,"sst","missing_value")
had.sst[had.sst==missvalue$value] <- NA
had.sst[had.sst==-1000] <- NAWe are interested in patterns of global variability, so the global mean needs to be removed at each step. But first investigate the global mean. The HadISST date data is formatted by year, and then the colmeans function is used to calculate the average sst. This is then aggregated together to create the annual mean data for HADISST.
# use annual averages with aggregate and using colMeans and rowMeans
# colmeans is used as the data is 3 dimensional and the output is from the mean of 3 them
had.year <- format(as.Date(had.date, format="%Y-%m-%d"),"%Y")
had.gmean <- colMeans(had.sst, na.rm = TRUE, dims=2)
had.annmean <- aggregate(had.gmean,by=list(had.year),FUN=mean,na.rm=TRUE)The following code is then used to create nicer colours for maps.
colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)Calculate the row means for HADSST data
had.avsst = rowMeans(had.sst,na.rm=FALSE,dims=2)Create a simple plot of the row means for the HadISST data, using the levelplot function.
levelplot(had.avsst,col.regions = pal(100))
However, adding lat and lon is a little more involved as levelplot likes to work with data frames. Expand grid will get latitude and longitude added to a data frame expanded on each others size. Use the average sst data to plot a new map, including the latitude and longitude.
grid <- expand.grid(x=had.lon, y=had.lat)
# our avsst that we want to plot can be added as a vector, basically expanding what we had
grid$had.avsst <- as.vector(had.avsst)
levelplot(had.avsst~x*y,grid,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Average SST'
) +
layer(sp.lines(world.coast))
Make annual averages to assess long term variability. The new variable had.sst will contain the annual sea surface temperature. Declare it here to have the correct size as it will speed up computation time. A loop is then used to load the annual data into average sst.
had.yrs <- had.annmean$Group.1
had.nyr <- length(had.yrs)
had.asst <- array(NA,c(dim(had.lon),dim(had.lat),had.nyr))
for (k in 1:had.nyr) {
had.asst[,,k] <- rowMeans(had.sst[,,had.year==had.yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}Annual averages have been created, so that the average will look the exact same. Need to add to the same data frame as earlier for plotting. And then plot.
grid$had.an_avsst <- as.vector(rowMeans(had.asst,na.rm=FALSE,dims=2))
levelplot(had.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, using a loop.
had.gmean <- colMeans(had.asst, na.rm = TRUE, dims=2)
for (k in 1:had.nyr){
had.asst[,,k]<-had.asst[,,k]-matrix(had.gmean[k],length(had.lon),length(had.lat))
}AMO
For the purpose of this analysis, the climate indice Atlantic Multidecadal Oscillation (AMO), was downloaded from http://climexp.knmi.nl/selectindex.cgi?id=c4c6e72298579658de28c6ec3cd004c4 as a netcdf file. First the data needed to be loaded into R and the date needed to be formatted into a standardised format for R to recognise and to conduct the analysis. The amo is monthly data and therefore needs to be multiplied by the number of days in a year divided by 12. From the output it is clear that the AMO data starts in 1854 and ends in 2019.
amo <- nc_open(file.path(getwd(), "iamo_hadsst_ts.nc"))
hadsst.amo <- ncvar_get(amo, 'AMO')
amo.date <- ncvar_get(amo, "time")
amo.tunits<-ncatt_get(amo,"time",attname="units")
amo.tustr<-strsplit(amo.tunits$value, " ")
amo.date<-as.character(as.Date(amo.date*365.25/12,origin=unlist(amo.tustr)[3]))
head(amo.date)## [1] "1854-01-15" "1854-02-14" "1854-03-16" "1854-04-16" "1854-05-16"
## [6] "1854-06-16"tail(amo.date)## [1] "2019-07-17" "2019-08-17" "2019-09-16" "2019-10-17" "2019-11-16"
## [6] "2019-12-17"The missing data then needs to be replaced with NA for a more accurate analysis. Use the print function to investigate the AMO dataset and find what the missing value number is to be replaced with NA.
print(amo)## File C:/Users/Natasha/Documents/College/Msc Climate Change/GY667/Assignment 2/iamo_hadsst_ts.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float AMO[time]
## long_name: Atlantic Multidecadal Variability index
## units: C
## _FillValue: 3.00000006085843e+33
##
## 1 dimensions:
## time Size:1992
## units: months since 1854-01-15
## standard_name: time
## long_name: time
## axis: T
## calendar: gregorian
##
## 27 global attributes:
## _NCProperties: version=1|netcdflibversion=4.6.0|hdf5libversion=1.10.0
## title: Difference between hadsst_0-60N_0-80W.dat and hadsst_60S-60N.dat
## description: AMO hadsst
## reference: Kennedy J.J., Rayner, N.A., Smith, R.O., Saunby, M. and Parker, D.E. (2011). Reassessing biases and other uncertainties in sea-surface temperature observations since 1850 part 2: biases and homogenisation. in press JGR Atmospheres
## file: hadsst_0-60N_0-80W.dat
## olderfile: HadSST.3.1.1.0.median.nc
## oldertitle: spatial statistic of Monthly 5 degree version of HadSST.3.1.1.0 - Median of realisations
## source: HadSST.3.1.1.0
## supplementary_information: Updates and supplementary information are available from http://www.hadobs.org
## ensemble_members: 100
## institution: KNMI Climate Explorer and Met Office Hadley Centre
## url: https://www.metoffice.gov.uk/hadobs/hadsst3/
## nco: 4.7.2
## geospatial_lat_max: 90
## geospatial_lon_min: -180
## geospatial_lon_units: degrees_east
## geospatial_lon_resolution: 5
## time_coverage_start: 1850-01-15
## time_coverage_end: 2019-01-15
## climexp_url: https://climexp.knmi.nl/select.cgi?hadsst3
## ave_region: lon= -80.000 0.000, lat= 0.000 60.000
## series2_ave_region: lon= -360.000 0.000, lat= -60.000 60.000
## series2_history: user oldenbor 2019-03-21 22:46:36 get_index HadSST.3.1.1.0.median.nc 0 360 -60 60\nThu Mar 21 22:45:26 2019: ncatted -a institution,global,a,c,Met Office Hadley Centre -a url,global,a,c,https://www.metoffice.gov.uk/hadobs/hadsst3/ HadSST.3.1.1.0.median.nc\n18/2/2019 converted to netcdf
## file2: hadsst_60S-60N.dat
## comment:
## history: 2019-04-24 14:08:47 bin/dat2nc data/iamo_hadsst_ts.dat i AMO_hadsst data/iamo_hadsst_ts.nc\nuser oldenbor 2019-03-21 22:46:36 normdiff hadsst_0-60N_0-80W.dat hadsst_60S-60N.dat none none\nuser oldenbor 2019-03-21 22:46:36 get_index HadSST.3.1.1.0.median.nc -80 0 0 60\nThu Mar 21 22:45:26 2019: ncatted -a institution,global,a,c,Met Office Hadley Centre -a url,global,a,c,https://www.metoffice.gov.uk/hadobs/hadsst3/ HadSST.3.1.1.0.median.nc\n18/2/2019 converted to netcdf
## Conventions: CF-1.0amo.fillvalue <- ncatt_get(amo, "AMO", "_FillValue")
hadsst.amo[hadsst.amo==amo.fillvalue$value] <- NA
amo.missvalue <- ncatt_get(amo,"AMO","missing_value")
hadsst.amo[hadsst.amo==amo.missvalue$value] <- NA
hadsst.amo[hadsst.amo==3.00000006085843e+33] <- NATo calculate the annual means for AMO, first get the AMO years and then get the mean of the AMO sst data and aggregate them together. Then rename the columns in the new data frame.
amo.year <- format(as.Date(amo.date, format ="%Y-%m-%d"),"%Y")
amo.annmean <- aggregate(hadsst.amo,by=list(amo.year), FUN= mean, na.rm= TRUE)
colnames(amo.annmean) <- c("Year", "AMO")
colnames(had.annmean) <- c("Year", "HadSST")Merge the HadISST annual mean dataframe and the AMO annual mean dataframe together. Once the two dataframes are merged. Visually inspect the data to see where the two dataframes overlap. Use square brackets to extract the data that overlaps into a new dataframe. The datasets overlap from 1879 to 2017, which is rows 26:164.
had.amo.merge <- merge(had.annmean, amo.annmean, all = TRUE)
had.amo.merge <- had.amo.merge[c(26:164),(1:3)]Similarly need to the the ASST data into the same length. The ASST data ranges from 1870 to 2017. Therefore need to get it to range from 1879 to 2017 so chop off the start of the data.
had.asst.1 <- had.asst[,,10:148]Now the data is in the same size set create a sst timeseries, using the plot function. This shows the variation of sst for AMO. From the plot AMO is relatively stable from 1880 to 1920, and then there is a step increase from 1920 to 1955, followed by a decrease until around 1990 where AMO begins to increase again.
plot(had.amo.merge$Year, had.amo.merge$AMO, type = 'l', xlab = 'Year', ylab = 'SST Anomaly', main = 'AMO')
The correlation between AMO sst and the HADisst average sst is then calcualted in a loop.
c.matrix <- matrix(NA,length(had.lon),length(had.lat))
t.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.1[i,j,][!is.na(had.asst.1[i,j,])])>2){
c.matrix[i,j] <- cor(had.asst.1[i,j,], had.amo.merge$AMO)
p.vals <- cor.test(had.asst.1[i,j,], had.amo.merge$AMO)
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). Then 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 the correlation between SST and AMO.
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SST with AMO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
It is clear from the map there is a significant relationship between SST in the North Atlantic and AMO. Similarly, there appears to be a correlation between AMO and SST in the North Pacific Ocean above Australia. The black dots indicate where there is a significant relationship. From the map one can conclude than positive SST is related to AMO.
Using LM. Linear regression
Calculate the correlation in a loop.
r.matrix <- matrix(NA,length(had.lon),length(had.lat))
s.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.1[i,j,][!is.na(had.asst.1[i,j,])])>2){
r.lm <- lm(had.asst.1[i,j,]~had.amo.merge$AMO)
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 the regression of SST with AMO.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SSTA with AMO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
From the map it is clear that the output from the regression is similar to the correlation of AMO and SST. It evident that there is significant relationship between AMO and the North Atlantic. In the map it appears that positive SST is related to AMO. However, within the North Atlantic there is a blue dot, also known as the subpolar gyre. This region has a negative SST index, but is also strongly correlated to AMO.
PDO
For the purpose of this analysis the correlation and regression of sst and Pacific Decadal Oscillation (PDO) also needs to be investigated. The data for PDO was also downloaded from the KNMI website in netcdf format. The data is read into R using the following code.
pdo <- nc_open(file.path(getwd(), "ipdo_hadsst3.nc"))
had3.pdo <- ncvar_get(pdo, 'index')Retrieve and format the time data for PDO so that it is a standardised format for R to recognise. The data is similarly in monthly format, therefore it needs to be multiplied by the number of days in a year divided by the number of months. It is clear from the output that the data ranges from 1884 to 2019.
pdo.date <- ncvar_get(pdo, 'time')
pdo.tunits<-ncatt_get(pdo,"time",attname="units")
pdo.tustr<-strsplit(pdo.tunits$value, " ")
pdo.date<-as.character(as.Date(pdo.date*365.25/12,origin=unlist(pdo.tustr)[3]))
head(pdo.date, n = 10)## [1] "1884-01-15" "1884-02-14" "1884-03-15" "1884-04-15" "1884-05-15"
## [6] "1884-06-15" "1884-07-15" "1884-08-15" "1884-09-14" "1884-10-14"tail(pdo.date, n = 10)## [1] "2019-03-17" "2019-04-17" "2019-05-17" "2019-06-16" "2019-07-17"
## [6] "2019-08-16" "2019-09-16" "2019-10-16" "2019-11-16" "2019-12-16"The missing values then need to be replaced with NA.
print(pdo)## File C:/Users/Natasha/Documents/College/Msc Climate Change/GY667/Assignment 2/ipdo_hadsst3.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float index[time]
## long_name:
## _FillValue: 3.00000006085843e+33
##
## 1 dimensions:
## time Size:1632
## units: months since 1884-01-15
## standard_name: time
## long_name: time
## axis: T
## calendar: gregorian
##
## 37 global attributes:
## _NCProperties: version=1|netcdflibversion=4.6.0|hdf5libversion=1.10.0
## title: Monthly 5 degree version of HadSST.3.1.1.0 - Median of realisations
## description: PDO HadSST3
## file: hadsst3-tglobal.nc
## source: HadSST.3.1.1.0
## reference: Kennedy J.J., Rayner, N.A., Smith, R.O., Saunby, M. and Parker, D.E. (2011). Reassessing biases and other uncertainties in sea-surface temperature observations since 1850 part 2: biases and homogenisation. in press JGR Atmospheres
## supplementary_information: Updates and supplementary information are available from http://www.hadobs.org
## ensemble_members: 100
## institution: KNMI Climate Explorer and Met Office Hadley Centre
## url: https://www.metoffice.gov.uk/hadobs/hadsst3/
## nco: 4.7.2
## geospatial_lat_min: -90
## geospatial_lat_max: 90
## geospatial_lat_units: degrees_north
## geospatial_lon_min: -180
## geospatial_lon_max: 180
## geospatial_lon_units: degrees_east
## geospatial_lat_resolution: 5
## geospatial_lon_resolution: 5
## time_coverage_start: 1850-01-15
## time_coverage_end: 2019-01-15
## climexp_url: https://climexp.knmi.nl/select.cgi?hadsst3
## series_institution: UK Met Office / Hadley Centre
## series_source_url: http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series//HadCRUT.4.6.0.0.monthly_ns_avg.txt
## series_source: https://www.metoffice.gov.uk/hadobs/hadcrut4/
## series_references: Morice, C. P., J. J. Kennedy, N. A. Rayner, and P. D. Jones (2012), Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 dataset, J. Geophys. Res., 117, D08101, doi:10.1029/2011JD017187.
## series_climexp_url: https://climexp.knmi.nl/getindices.cgi?UKMOData/hadcrut4_ns_avg
## series_history: retrieved Sun Mar 17 12:14:17 MET 2019
## series_title: global mean temperature
## series_file: /home/oldenbor/climexp/UKMOData/hadcrut4_ns_avg.dat
## pattern_history: user oldenbor 2019-04-18 17:52:36 eof ./hadsst3-tglobal.nc 1 normalize varspace mon 1 ave 12 lon1 100 lon2 260 lat1 20 lat2 65 begin 1900 eof1_hadsst.nc\nuser oldenbor 2019-04-18 17:52:35 subfieldseries UKMOData/HadSST.3.1.1.0.median.nc UKMOData/hadcrut4_ns_avg.dat ./hadsst3-tglobal.nc\nThu Mar 21 22:45:26 2019: ncatted -a institution,global,a,c,Met Office Hadley Centre -a url,global,a,c,https://www.metoffice.gov.uk/hadobs/hadsst3/ HadSST.3.1.1.0.median.nc\n18/2/2019 converted to netcdf
## pattern_title: Monthly 5 degree version of HadSST.3.1.1.0 - Median of realisations
## pattern_file: eof1_hadsst.nc
## scale_factor: 3.00000
## comment:
## history: 2019-04-24 13:58:21 bin/dat2nc data/ipdo_hadsst3.dat i PDO_HadSST3 data/ipdo_hadsst3.nc\nuser oldenbor 2019-04-18 17:52:37 scaleseries 3 aap.dat\nuser oldenbor 2019-04-18 17:52:37 patternfield hadsst3-tglobal.nc eof1_hadsst.nc eof1 1\nuser oldenbor 2019-04-18 17:52:35 subfieldseries UKMOData/HadSST.3.1.1.0.median.nc UKMOData/hadcrut4_ns_avg.dat ./hadsst3-tglobal.nc\nThu Mar 21 22:45:26 2019: ncatted -a institution,global,a,c,Met Office Hadley Centre -a url,global,a,c,https://www.metoffice.gov.uk/hadobs/hadsst3/ HadSST.3.1.1.0.median.nc\n18/2/2019 converted to netcdf
## Conventions: CF-1.0pdo.fillvalue <- ncatt_get(pdo, "index", "_FillValue")
had3.pdo[had3.pdo==pdo.fillvalue$value] <- NA
pdo.missvalue <- ncatt_get(pdo,"index","missing_value")
had3.pdo[had3.pdo==pdo.missvalue$value] <- NA
had3.pdo[had3.pdo==3.00000006085843e+33] <- NACreate annual means for the pdo data, format the date by just year and then get the mean of the PDO sst data, these are then aggregated together to create a new dataframe. Rename the columns in the PDO data.
pdo.year <- format(as.Date(pdo.date, format = "%Y-%m-%d"),"%Y")
pdo.annmean <- aggregate(had3.pdo,by=list(pdo.year), FUN= mean, na.rm= TRUE)
head(pdo.annmean)## Group.1 x
## 1 1884 0.4848025
## 2 1885 1.6688299
## 3 1886 -0.2291120
## 4 1887 NaN
## 5 1888 NaN
## 6 1889 NaNcolnames(pdo.annmean) <- c("Year", "PDO")Merge the HADSIIT annual mean data and the pdo annual mean data together. Then visually inspect the merged data to see where they overlap. They overlap from 1899 to 2017, which corresponds to rows 30:148. Then chop off the data doesn’t correspond by using square brackets.
had.pdo.merge <- merge(had.annmean, pdo.annmean, all = TRUE)
had.pdo.merge <- had.pdo.merge[c(30:148),(1:3)]Now the data is in the same size set create a sst timeseries using the plot function. From the output it is clear that there is a lot of variability in PDO, however according to the plot, there appears to be about a 17 year cycle where PDO increases and then is followed by a decrease.
plot(had.pdo.merge$Year, had.pdo.merge$PDO, type = 'l', xlab = 'Year', ylab = 'SST Anomaly', main = 'PDO')
ASST ranges from 1870 to 2017 therefore need to chop to match 1899 to 2017.
had.asst.2 <- had.asst[,,30:148]Calculate the correlation between HADiSST sst and PDO data in a loop.
c.matrix <- matrix(NA,length(had.lon),length(had.lat))
t.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.2[i,j,][!is.na(had.asst.2[i,j,])])>2){
c.matrix[i,j] <- cor(had.asst.2[i,j,], had.pdo.merge$PDO)
p.vals <- cor.test(had.asst.2[i,j,], had.pdo.merge$PDO)
t.matrix[i, j] <- p.vals$p.value
}
}
}Need to 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 the correlation between PDO and sst.
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with PDO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
There is a significant correlation between PDO and SST in the Northern and Southern Pacific Ocean. There appears to be a strong positive correlation between positive SST and PDO along the West coast of Northern and Southern America. There is also a significant relationship with negative SST and PDO in the North and South Pacific Ocean. Similarly there is a positive correlation between SST and PDO in the Indian Ocean.
Using LM. Linear regression
Calculate the correlation in a loop.
r.matrix <- matrix(NA,length(had.lon),length(had.lat))
s.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.2[i,j,][!is.na(had.asst.2[i,j,])])>2){
r.lm <- lm(had.asst.2[i,j,]~had.pdo.merge$PDO)
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 the regression between PDO and SST.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SSTA with PDO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
From mapping the regression of PDO and SST there appears to be a significant relationship along the West coast of North and South America. From the plot there are both negative and positive SST that influence PDO.
ENSO
The data for the climate indice El Nino/ Southern Oscillation (ENSO) is downloaded from KNMI website, in a netcdf format. The data is loaded using the following code.
enso <- nc_open(file.path(getwd(), "ihadisst1_nino4a.nc"))
hadisst.enso <- ncvar_get(enso, 'Nino4')The time variable is then extracted and is formatted to a standardised format. Similar to the two previous climate indices the data was initially in monthly. Therefore it was multiplied by 365.25 and divided by 12 to get it into days. The data runs from 1870 to 2019.
enso.date <- ncvar_get(enso, "time")
enso.tunits <- ncatt_get(enso,"time", attname = "units")
enso.tustr <- strsplit(enso.tunits$value, " ")
enso.date <- as.character(as.Date(enso.date*365.25/12,origin=unlist(enso.tustr)[3]))
head(enso.date, n = 10)## [1] "1870-01-15" "1870-02-14" "1870-03-16" "1870-04-16" "1870-05-16"
## [6] "1870-06-16" "1870-07-16" "1870-08-16" "1870-09-15" "1870-10-15"tail(enso.date, n = 10)## [1] "2019-03-18" "2019-04-17" "2019-05-18" "2019-06-17" "2019-07-17"
## [6] "2019-08-17" "2019-09-16" "2019-10-17" "2019-11-16" "2019-12-17"The missing data is filled with NA’s so that it is a more accurate analysis. The missing value is found by using the print function.
print(enso)## File C:/Users/Natasha/Documents/College/Msc Climate Change/GY667/Assignment 2/ihadisst1_nino4a.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float Nino4[time]
## long_name: HadISST1 Nino4 index
## units: K
## _FillValue: 3.00000006085843e+33
##
## 1 dimensions:
## time Size:1800
## units: months since 1870-01-15
## standard_name: time
## long_name: time
## axis: T
## calendar: gregorian
##
## 21 global attributes:
## _NCProperties: version=1|netcdflibversion=4.6.0|hdf5libversion=1.10.0
## title: spatial statistic of Monthly version of HadISST sea surface temperature component
## description: HadISST 1.1 monthly average sea surface temperature
## file: HadISST_sst.nc
## cdi: Climate Data Interface version ?? (http://mpimet.mpg.de/cdi)
## source: HadISST
## institution: KNMI Climate Explorer and Met Office Hadley Centre
## reference: Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., Kent, E. C., Kaplan, A. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century J. Geophys. Res.Vol. 108, No. D14, 4407 10.1029/2002JD002670
## supplementary_information: Updates and supplementary information will be available from http://www.metoffice.gov.uk/hadobs/hadisst
## comment: \nData restrictions: for academic research use only. Data are Crown copyright see (http://www.opsi.gov.uk/advice/crown-copyright/copyright-guidance/index.htm)
## cdo: Climate Data Operators version 1.9.3 (http://mpimet.mpg.de/cdo)
## geospatial_lat_max: 90
## geospatial_lon_min: -180
## geospatial_lon_units: degrees_east
## geospatial_lon_resolution: 1
## time_coverage_start: 1870-01-15
## time_coverage_end: 2019-01-15
## climexp_url: https://climexp.knmi.nl/select.cgi?hadisst1
## ave_region: lon= -200.000 -150.000, lat= -5.000 5.000
## history: 2019-04-24 14:17:06 bin/dat2nc data/ihadisst1_nino4a.dat i NINO4 data/ihadisst1_nino4a.nc\nuser oldenbor 2019-03-21 22:46:10 get_index HadISST_sst.nc -200 -150 -5 5\nThu Mar 21 22:45:40 2019: cdo -r -f nc4 -z zip -selvar,sst -setctomiss,-1000. HadISST1_sst_tmp.nc HadISST_sst.nc\n12/3/2019 converted to netcdf from pp format
## Conventions: CF-1.0enso.fillvalue <- ncatt_get(enso, "Nino4", "_FillValue")
hadisst.enso[hadisst.enso==enso.fillvalue$value] <- NA
enso.missvalue <- ncatt_get(enso,"Nino4","missing_value")
hadisst.enso[hadisst.enso==enso.missvalue$value] <- NA
hadisst.enso[hadisst.enso==3.00000006085843e+33] <- NAThe annual means for ENSO data are then calculated and the aggregate function is used to create a new dataframe, that contains the annually averaged ENSO data.
enso.year <- format(as.Date(enso.date, format = "%Y-%m-%d"),"%Y")
enso.annmean <- aggregate(hadisst.enso,by=list(enso.year), FUN= mean, na.rm= TRUE)
head(enso.annmean)## Group.1 x
## 1 1870 -0.6953413
## 2 1871 -0.4309183
## 3 1872 -0.6376590
## 4 1873 -0.5882057
## 5 1874 -0.7470818
## 6 1875 -0.7695382Rename column names for annually averaged PDO data
colnames(enso.annmean) <- c("Year", "ENSO")Merge the HADSIIT annually averaged data and the PDO annually averaged data together. Visually inspect to see where the two data frames overlap. They overlap from 1870 to 2017, which is rows 1:148. Use the square brackets to chop the data.
had.enso.merge <- merge(had.annmean, enso.annmean, all = TRUE)
had.enso.merge <- had.enso.merge[c(1:148),(1:3)]Now the data is in the same size set create a sst timeseries for PDO using the plot function. From the output there are no evident cycles, however the variability is clear.
plot(had.enso.merge$Year, had.enso.merge$ENSO, type = 'l', xlab = 'Year', ylab = 'SST Anomaly', main = 'ENSO')
ASST ranges from 1870 to 2017, however the ENSO data also ranges from 1870 to 2017. Therefore do not need to chop the data to match.
had.asst.3 <- had.asstCalculate the correlation between ENSO and HadiSST sst in a loop.
c.matrix <- matrix(NA,length(had.lon),length(had.lat))
t.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.3[i,j,][!is.na(had.asst.3[i,j,])])>2){
c.matrix[i,j] <- cor(had.asst.3[i,j,], had.enso.merge$ENSO)
p.vals <- cor.test(had.asst.3[i,j,], had.enso.merge$ENSO)
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 the correlation between ENSO and sst.
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with ENSO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
From the map ENSO has a strong negative and positive correlation with SST. Along the coast of North and South America there is a strong positive correlation, where positive SST are associated with ENSO. There is also a positive relationship between ENSO and SST in the Indian ocean. In the North and South Pacific there is a strong negative correlation between ENSO and SST.
Using LM. Linear regression
Calculate the correlation in a loop.
r.matrix <- matrix(NA,length(had.lon),length(had.lat))
s.matrix <- matrix(NA,length(had.lon),length(had.lat))
for (i in 1:dim(had.lon)) {
for (j in 1:dim(had.lat)) {
if (length(had.asst.3[i,j,][!is.na(had.asst.3[i,j,])])>2){
r.lm <- lm(had.asst.3[i,j,]~had.enso.merge$ENSO)
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 the regression of SST with ENSO.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SSTA with ENSO', had.lon, ',Lat=', had.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
It is clear from the map that ENSO occurs in regions with both positive and negative SST values. However, ENSO occurs in more areas where SST values are negative, especially in the North and South Pacific. Along the western coast of North and South America, ENSO is associated with positive SST values.
Sea Level Pressure (SLP)
The sea level pressure data was downloaded from Work shop 3 on moodle.
First need to load the data and then extract the variables of interest. The longitude ranges from 0 to 360, however the range needs to be -180 to 180 for the analysis. This is done using the ifelse code, once completed however the data is in the wrong order, so use the rev function.
slp.mon.mean <- nc_open(file.path(getwd(), "slp.mon.mean.nc"))
slp.lat <- ncvar_get(slp.mon.mean,"lat")
slp.lon <- ncvar_get(slp.mon.mean,"lon")
slp.lon <- ifelse(slp.lon>180, slp.lon-360,slp.lon)
slp.lon <- rev(slp.lon)
m.slp <- ncvar_get(slp.mon.mean, "slp")
slp.date <- ncvar_get(slp.mon.mean, "time")The date needs to be in the correct format to conduct the analysis and for R to recognise. Divide by 24 to get the data into days. From the output it is clear that the data starts in 1948 and ends in 2017.
slp.date <- ncvar_get(slp.mon.mean, "time")
slp.tunits<-ncatt_get(slp.mon.mean,"time",attname="units")
slp.tustr<-strsplit(slp.tunits$value, " ")
slp.date<-as.character(as.Date(slp.date/24,origin=unlist(slp.tustr)[3]))
head(slp.date)## [1] "1948-01-01" "1948-02-01" "1948-03-01" "1948-04-01" "1948-05-01"
## [6] "1948-06-01"tail(slp.date)## [1] "2017-11-01" "2017-12-01" "2018-01-01" "2018-02-01" "2018-03-01"
## [6] "2018-04-01"Replace missing values with NA, to find the missing value number the print function was used.
print(slp.mon.mean)## File C:/Users/Natasha/Documents/College/Msc Climate Change/GY667/Assignment 2/slp.mon.mean.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float slp[lon,lat,time]
## long_name: Sea Level Pressure
## valid_range: 870
## valid_range: 1150
## units: millibars
## add_offset: 0
## scale_factor: 1
## missing_value: -9.96920996838687e+36
## precision: 1
## least_significant_digit: 1
## var_desc: Sea Level Pressure
## level_desc: Sea Level
## statistic: Mean
## parent_stat: Other
## dataset: NCEP Reanalysis Derived Products
## actual_range: 955.560852050781
## actual_range: 1082.55822753906
##
## 3 dimensions:
## lat Size:73
## units: degrees_north
## actual_range: 90
## actual_range: -90
## long_name: Latitude
## standard_name: latitude
## axis: Y
## lon Size:144
## units: degrees_east
## long_name: Longitude
## actual_range: 0
## actual_range: 357.5
## standard_name: longitude
## axis: X
## time Size:844 *** is unlimited ***
## long_name: Time
## delta_t: 0000-01-00 00:00:00
## prev_avg_period: 0000-00-01 00:00:00
## standard_name: time
## axis: T
## units: hours since 1800-01-01 00:00:0.0
## actual_range: 1297320
## actual_range: 1913112
##
## 8 global attributes:
## description: Data is from NMC initialized reanalysis
## (4x/day). These are the 0.9950 sigma level values.
## platform: Model
## Conventions: COARDS
## NCO: 20121012
## history: Thu May 4 18:12:35 2000: ncrcat -d time,0,622 /Datasets/ncep.reanalysis.derived/surface/slp.mon.mean.nc ./surface/slp.mon.mean.nc
## Mon Jul 5 23:22:35 1999: ncrcat slp.mon.mean.nc /Datasets/ncep.reanalysis.derived/surface/slp.mon.mean.nc /dm/dmwork/nmc.rean.ingest/combinedMMs/slp.mon.mean.nc
## /home/hoop/crdc/cpreanjuke2farm/cpreanjuke2farm Thu Oct 26 23:42:16 1995 from pre.sig995.85.nc
## created 95/02/06 by Hoop (netCDF2.3)
## Converted to chunked, deflated non-packed NetCDF4 2014/09
## title: monthly mean slp from the NCEP Reanalysis
## References: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
## dataset_title: NCEP-NCAR Reanalysis 1fillvalue <- ncatt_get(slp.mon.mean,"slp","_FillValue")
m.slp[m.slp==fillvalue$value] <- NA
missvalue <- ncatt_get(slp.mon.mean,"slp","missing_value")
m.slp[m.slp==missvalue$value] <- NA
m.slp[m.slp==--9.96920996838687e+36] <- NAThe SLP date data is formatted by year, and then the colmeans function is used to calculate the average SLP. This is then aggregated together to create the annual mean data for SLP.
slp.year <- format(as.Date(slp.date, format="%Y-%m-%d"),"%Y")
slp.gmean <- colMeans(m.slp, na.rm = TRUE, dims=2)
slp.annmean <- aggregate(slp.gmean,by=list(slp.year),FUN=mean,na.rm=TRUE)Row means are then calculated for the SLP data and saved in a new data frame.
slp.avslp = rowMeans(m.slp,na.rm=FALSE,dims=2)The level plot function is used to create a map of average SLP.
levelplot(slp.avslp,col.regions = pal(100))
The grid is then expanded to include the latitude and longitude of the data frame. Then plotted, with an extra layer containing the lines of the world coast.
grid <- expand.grid(x=slp.lon, y=slp.lat)
grid$slp.avslp <- as.vector(slp.avslp)
levelplot(slp.avslp~x*y,grid,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Average SLP'
) +
layer(sp.lines(world.coast))
For the purpose of this analysis annual averages are needed.
slp.yrs <- slp.annmean$Group.1
slp.nyr <- length(slp.yrs)
slp.asst <- array(NA,c(dim(slp.lon),dim(slp.lat),slp.nyr)) A loop is then used to fill the annual data into slp.asst.
for (k in 1:slp.nyr) {
slp.asst[,,k] <- rowMeans(m.slp[,,slp.year==slp.yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}Annual averages have been created. Add this data to the grid used earlier for plotting and plot again.
grid$slp.an_avslp <- as.vector(rowMeans(slp.asst,na.rm=TRUE,dims=2))
levelplot(slp.an_avslp~x*y, data=grid,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Annually Averaged SLP') +
layer(sp.lines(world.coast)) 
The global mean is then removed from each year using a traditional loop.
slp.gmean <- colMeans(slp.asst, na.rm = TRUE, dims=2)
for (k in 1:slp.nyr){
slp.asst[,,k]<-slp.asst[,,k]-matrix(slp.gmean[k],length(slp.lon),length(slp.lat))
}AMO and SLP
The data from AMO was imported earlier and was formatted to a standardised format. Need to get the dataframe with annual means for AMO and compare it to SLP annual means. This can be done using the head or tail function. Need to ensure that SLP and AMO have the same column name for year before merging.
tail(amo.annmean)## Year AMO
## 161 2014 -0.04044478
## 162 2015 -0.15994082
## 163 2016 0.02649527
## 164 2017 0.16055811
## 165 2018 -0.13220485
## 166 2019 -0.10153560tail(slp.annmean)## Group.1 x
## 66 2013 1010.656
## 67 2014 1010.647
## 68 2015 1009.834
## 69 2016 1009.952
## 70 2017 1010.193
## 71 2018 1010.613colnames(amo.annmean) <- c("Year", "AMO")
colnames(slp.annmean) <- c("Year", "SLP")Merge SLP and AMO annual mean data together into one data frame. Visually inspect the dataframe to see where they both overlap and chop off the bits that don’t using square brackets. From visually inspecting the data they both overlap from 1948 to 2018, which corresponds to rows 95:165.
slp.amo.merge <- merge(slp.annmean, amo.annmean, all = TRUE)
slp.amo.merge <- slp.amo.merge[c(95:165),(1:3)]Similarly need to make sure that the annual SLP data is the same length. As SLP ranges from 1948 to 2018 there is no need to alter the data.
slp.asst.1 <- slp.asstNow the data is in the same size set create a SLP timeseries using the plot function. From the output of the plot there is a clear decrease from 1955 to about 1985, and then AMO starts to increase again to about 2012. This may suggest an increasing and decreasing cycle in AMO of about 30 years.
plot(slp.amo.merge$Year, slp.amo.merge$AMO, type = 'l', xlab = 'Year', ylab = 'SLP Anomaly', main = 'AMO')
Calculate the correlation between SLP and AMO in a loop
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(slp.asst.1[i,j,][!is.na(slp.asst.1[i,j,])])>2){
c.matrix[i,j] <- cor(slp.asst.1[i,j,], slp.amo.merge$AMO)
p.vals <- cor.test(slp.asst.1[i,j,], slp.amo.merge$AMO)
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 the correlation between SLP and AMO
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with AMO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
grid$corr <- as.vector(c.matrix)
grid$pval <- as.vector(t.matrix)There is a strong positive correlation between SLP and AMO in the Indian Ocean and the North Pacific Ocean, particularly around Australia, where positive SLP values are associated with positive AMO. There is also a significant relationship between AMO and negative SLP values in the North Atlantic Ocean.
Using LM. Linear regression
Calculate the correlation in a loop.
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(slp.asst.1[i,j,][!is.na(slp.asst.1[i,j,])])>2){
r.lm <- lm(slp.asst.1[i,j,]~slp.amo.merge$AMO)
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 the regression between SLP and AMO.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SLP with AMO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
From the map there are clear areas of significance which are marked by the black dots. In the North Atlantic it is clear that SLP influences AMO. Similarly, around Australia in the Indian Ocean, positive SLP is associated with AMO.
PDO
The PDO data was imported earlier and put into a standardised format. The head function was used to have a quick look at the data. The data overlaps from 1948 to 2018 which corresponds to rows 65 to 135.
head(slp.annmean)## Year SLP
## 1 1948 1012.461
## 2 1949 1012.047
## 3 1950 1012.388
## 4 1951 1012.453
## 5 1952 1012.698
## 6 1953 1012.209head(pdo.annmean)## Year PDO
## 1 1884 0.4848025
## 2 1885 1.6688299
## 3 1886 -0.2291120
## 4 1887 NaN
## 5 1888 NaN
## 6 1889 NaNMerge the SLP annual mean data and the PDO annual mean data together. Visually inspect the merged data and remove the bits that do not overlap by using the square brackets.
slp.pdo.merge <- merge(slp.annmean, pdo.annmean, all = TRUE)
slp.pdo.merge <- slp.pdo.merge[c(65:135),(1:3)]Now the data is in the same size set create a slp timeseries using the plot function. From the output there appears to be a cycle from about 1950 to 1975, where PDO is quite low with some variability. However, from 1975 to about 2000, PDO is relatively high with some variability. Followed by a relatively low period to about 2018.From about 2000 to 2018. This may indicate a cycle of about 20 years where PDO shifts from High to low
plot(slp.pdo.merge$Year, slp.pdo.merge$PDO, type = 'l', xlab = 'Year', ylab = 'SLP Anomaly', main = 'PDO')
SLP data ranges from 1948 to 2018 therefore no need to alter the data to get it on the same timescale
slp.asst.2 <- slp.asstCalculate the correlation between SLP and PDO in a loop.
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(slp.asst.2[i,j,][!is.na(slp.asst.2[i,j,])])>2){
c.matrix[i,j] <- cor(slp.asst.2[i,j,], slp.pdo.merge$PDO)
p.vals <- cor.test(slp.asst.2[i,j,], slp.pdo.merge$PDO)
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 correlation of SLP with PDO using the levelplot function.
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with PDO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
From the output it is clear that there is a strong positive and negative correlation between SLP and PDO. In the North Atlantic, South Atlantic, North Pacific and South Pacific Oceans, positive SLP has a strong correlation with PDO. There is a strong correlation with negative SLP and PDO off the coast of Asia and in the Southern Ocean.
Using LM. Linear regression
Calculate the correlation in a loop.
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(slp.asst.2[i,j,][!is.na(slp.asst.2[i,j,])])>2){
r.lm <- lm(slp.asst.2[i,j,]~slp.pdo.merge$PDO)
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 the regression of PDO with SLP.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SLP with PDO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
The map shows the regression between SLP and PDO. The output highlights that positive SLP in the Pacific and the Atlantic Ocean is associated with PDO. Similarly, negative SLP in the North Pacific Ocean off the coast of Asia is also associated with PDO.
El Nino/ Southern Oscillation (ENSO)
The ENSO data has been previously imported into R and manipulated so that it is in a standardised format. Use the head function to quickly preview the ENSO and SLP annual mean data.
head(slp.annmean)## Year SLP
## 1 1948 1012.461
## 2 1949 1012.047
## 3 1950 1012.388
## 4 1951 1012.453
## 5 1952 1012.698
## 6 1953 1012.209head(enso.annmean)## Year ENSO
## 1 1870 -0.6953413
## 2 1871 -0.4309183
## 3 1872 -0.6376590
## 4 1873 -0.5882057
## 5 1874 -0.7470818
## 6 1875 -0.7695382Merge the SLP and ENSO data together. This will enable visual inspection of the data side by side to see where they overlap. From inspecting the data they overlap from 1948 to 2018 which corresponds with rows 79:149.
slp.enso.merge <- merge(slp.annmean, enso.annmean, all = TRUE)
slp.enso.merge <- slp.enso.merge[c(79:149),(1:3)]Now the data is in the same size set create a SLP timeseries using the plot function.It is clear from the graph that from the years 1948 to 2018, ENSO is quite variable but tends to follow short periods of sharp increase, followed by short periods of sharp decrease.
plot(slp.enso.merge$Year, slp.enso.merge$ENSO, type = 'l', xlab = 'Year', ylab = 'SLP Anomaly', main = 'ENSO')
SLP ranges from 1948 to 2018 therefore as the merged data already matches there is no need for alteration.
slp.asst.3 <- slp.asstCalculate the correlation of SLP with ENSO in a loop.
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(had.asst.3[i,j,][!is.na(slp.asst.3[i,j,])])>2){
c.matrix[i,j] <- cor(slp.asst.3[i,j,], slp.enso.merge$ENSO)
p.vals <- cor.test(slp.asst.3[i,j,], slp.enso.merge$ENSO)
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 the correlation of SLP with ENSO
levelplot(corr~x*y, data=grid , xlim=c(-180,180),ylim=c(-90,90), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SLP with ENSO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.0001, col = "black"))
There is a strong positive correlation between positive SLP and ENSO in the North and South Pacific, and the Atlantic Ocean. There is also a significant relationship between negative SLP in the Indian Ocean and the North Pacific Ocean, particularly around Australia and the coast of Asia.
Using LM. Linear regression
Calculate the correlation in a loop.
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(slp.asst.3[i,j,][!is.na(slp.asst.3[i,j,])])>2){
r.lm <- lm(slp.asst.3[i,j,]~slp.enso.merge$ENSO)
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 the regression between SLP and ENSO.
levelplot(reg~x*y, data=grid, at=c(-30:30)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',
main=paste0('Regression of SLP with ENSO', slp.lon, ',Lat=', slp.lat)) +
layer(sp.lines(world.coast)) +
layer(sp.points(sig, pch = 20, cex = 0.001, col = "black"))
It is clear that there is a significant relationship between positive SLP and ENSO in the Atlantic and the Pacific Ocean. Similarly, there is a significant relationship between ENSO and negative SLP in the Indian Ocean and Northern Pacific, particularly off the coast of Asia and the region around Australia.
AMO
The Atlantic Multidecadal Oscillation (AMO) is a mode of multidecadal climate variability, that switches between warm and cool phases of sea surface temperatures (SST) over large areas of the Northern hemisphere. AMO causes global impacts on the climate such as, rain in the Sahel and hurricanes in the Atlantic. Positive AMO anomalies are associated with decreased mean sea level pressure and influences the climate over the Atlantic and Europe bringing cyclonic pressure anomalies. The correlation between SLP and AMO can be seen above, where negative slp dominates the Atlantic. Certain models highlight that AMO is related to the strength of the overturning circulation in the Atlantic Ocean (AMOC) (Knight, Folland and Scaife, 2006). AMOC is driven by differences in density related to deep-water formation in the North Atlantic. AMOC transports warm water nortwards in the Atlantic, this explains the difference in temperature that can be seen in previous maps. The North Atlantic has higher SST thatn the South Atlantic ocean. AMOC could potentially be slowed down in the future as a result of a reduction in surface ocean density (Rahmstorf et al., 2015). When AMOC is weak or shuts down, the North Atlantic Ocean cools as a result of reduced northward heat transport, this would have drastic impacts on the climate (Clement et al., 2015). In 1970 a large amount of freshwater entered the North Atlantic and caused AMOC to slow down considerably. However, AMOC recovered and is now facing further instability as a result of further freshening from melting ice sheets (Rahmstorf et al., 2015).
PDO
Pacific Decadal Oscillation (PDO) is a multi-decadal pattern of climate variability that affects the Pacific Ocean (Barnard et al., 2015). PDO has strong impacts on the Southern Hemisphere. PDO events occur typically for about 20 to 30 years. During warm PDO phases SST tends to be cool in central North Pacific which coincides with warm SST along the west coast of America. This can be seen in the correlation map between SST and PDO. SLP can enhance winds over the northern subtropical pacific. Warm phases of PDO bring dry periods to Japan and Eastern Australia. Similarly, warm phases of PDO bring wet periods to South East Brazil and western Australia (Mantua and Hare, 2002). The location and strength of the Kuroshio current influences PDO and SST, which can be seen in the map of the correlation between PDO and SST (Schneider and Cornuelle, 2005).
ENSO
El Nino/ Southern Oscillation (ENSO) describes the interannual variability in SST and SLP over the equatorial Pacific (Barnard et al., 2015). ENSO influences the tropical Atlantic region (Clement et al., 2015). ENSO typically varies from about 2 to 7 years. ENSO is a source of interannual climate variations in the Pacific and the tropics (Mantua and Hare, 2002). ENSO is associated with low temperatures in the North Pacific, and with warm anomalies in the eastern Pacific (Schneider and Cornuelle, 2005). When El Nino has positive index numbers, the SST is higher and SLP is lower in the eastern equatorial pacific. However, when El Nino has negative index numbers there is higher SLP and lower SST in the eastern equatorial pacific. This pattern can be seen in the maps for the correlation for SST and SLP with ENSO (Barnard et al., 2015).
Bibliography:
Barnard, P.L., Short, A.D., Harley, M.D., Splinter, K.D., Vitousek, S., Turner, I.L., Allan, J., Banno, M., Bryan, K.R., Doria, A. and Hansen, J.E., 2015. Coastal vulnerability across the Pacific dominated by El Nino/Southern oscillation. Nature Geoscience, 8(10), p.801.
Clement, A., Bellomo, K., Murphy, L.N., Cane, M.A., Mauritsen, T., Rädel, G. and Stevens, B., 2015. The Atlantic Multidecadal Oscillation without a role for ocean circulation. Science, 350(6258), pp.320-324.
Knight, J.R., Folland, C.K. and Scaife, A.A., 2006. Climate impacts of the Atlantic multidecadal oscillation. Geophysical Research Letters, 33(17).
Mantua, N.J. and Hare, S.R., 2002. The Pacific decadal oscillation. Journal of oceanography, 58(1), pp.35-44.
Rahmstorf, S., Box, J.E., Feulner, G., Mann, M.E., Robinson, A., Rutherford, S. and Schaffernicht, E.J., 2015. Exceptional twentieth-century slowdown in Atlantic Ocean overturning circulation. Nature climate change, 5(5), pp.475-480.
Schneider, N. and Cornuelle, B.D., 2005. The forcing of the Pacific decadal oscillation. Journal of Climate, 18(21), pp.4355-4373.