Workshop 3 : Ocean Circulation Patterns

Preliminaries

Load the relevant libraries for this assignment.

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-1
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

Clear variables and any open plots before loading in the data

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

Set the working directory.

setwd("C:/Users/micha/Downloads/Workshop3")

Task 1 - Load NetCDF data & convert date format

Set the path to the data

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

In the directory path, sea surface temperature data is stored. This data is from the Hadley Centre, otherwise known as the UK Met Office. This data can be located at https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html Coastline data, which is used just for plotting can be found at http://www.naturalearthdata.com/downloads/110m-physical-vectors/

Load in the data

load("world.coast.RData")

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))
for (k in 1:nyr) {
  asst[,,k] <- rowMeans(sst[,,year==yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}

What we have done is create annual average, so the average looks just the same. Next is to add to the same data frame as earlier for plotting

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

Remove the global mean from each year

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

in this example, we’re going to take a timeseries of sst at a point lon0, lat0

and look at the patterns of correlation of sst with this

tseries will just be the asst value at a given point lon0, lat0

near ireland lon0 <- -10.5, lat0 <- 51.5

subpolar gyre lon0 <- -35.5, lat0 <- 59.5

try a few of your own

note: keep the .5 at the end to find a matching lat/lon

Next we are going to make a timeseries of sst. 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),]

Look at the time-series

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

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

grid$corr <- as.vector(cmatrix)

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

Task 2: repeat this correlation map using a point north of the Gulf Stream

Change the Lat and Long to North of the Gulf Stream

lon0 <- -60.5
lat0 <- 40.5
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon0, ',Lat=', lat0))

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)
  }
}
grid$corr <- as.vector(cmatrix)
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))