GY667_WS3_18251197

WORKSHOP 3: OCEAN CIRCULATION PATTERNS

load the relevant packages

install.packages(“sp”) # etc.

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

## rgdal: version: 1.3-4, (SVN revision 766)
##  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/Thais/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/Thais/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)

## Warning: package 'maps' was built under R version 3.5.3

library(OpenStreetMap)

(not necessary) clear variables and any open plots before starting

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

Set your working directory, yours will be different

setwd("D:/Documentos D/MSc Climate Change/GY667 - The Ocean and Climate Change/WS3")

Set the path to the data

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

in the directory path store the following data

- sea surface temperature data from Hadley Centre (UK Met Office): https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html HadISST_sst.nc.gz

- coastline data (just for plotting): http://www.naturalearthdata.com/downloads/110m-physical-vectors/ download coastline

read coastline using readOGR as a Spatial Data Frame

part of the rgdal package

higher resolution coastlines available but not needed for this application

world.coast <- readOGR(dsn = file.path(path,“ne_110m_coastline”), layer = “ne_110m_coastline” ) # convert to the correct projections

Task 1:

load HadISST_sst.nc

setwd("D:/Documentos D/MSc Climate Change/GY667 - The Ocean and Climate Change/WS3")
nc   <- nc_open(file.path(getwd(),"HadISST_sst.nc"))

lat <-

lon <-

sst <-

date <-

lat  <- ncvar_get(nc,"latitude") # load the latitude
lon  <- ncvar_get(nc,"longitude") # load the longitude
sst <- ncvar_get(nc,"sst"); # load sea surface temperature
date   <- ncvar_get(nc,"time"); # load the 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

we’re interested in patterns of global variability

so we need to remove the global mean at each step

but first investigate the global mean

use annual averages with aggregate and using colMeans and rowMeans

you could alternatively use apply as in the first workshop

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

year <- format(as.Date(date_nc, format="%Y-%m-%d"),"%Y")
gmean <- colMeans(sst, na.rm = TRUE, dims=2)
annmean <- aggregate(gmean,by=list(year),FUN=mean,na.rm=TRUE)

alternatively comment out the following two lines, what difference does excluding NA values make?

avsst = rowMeans(sst,na.rm=TRUE,dims=2)

avsst_sea_ice = rowMeans(sst,na.rm=TRUE,dims=2) #biased on  sea ice

avsst = rowMeans(sst,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

simple e.g.

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

instantly we have a colorbar and the plot looks nice

however, adding lat and lon is a little more involved as levelplot likes to work with data frames

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

grid <- expand.grid(x=lon, y=lat)

our avsst that we want to plot can be added as a vector, basically 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))

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

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

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

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

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
}

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

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

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

lon0 <- -35.5
lat0 <- 59.5 
sst_ts_subpolar<-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))

lets calculate the correlation in a loop also

can you remove one of these loops? does this speed up computation time?

hints: help(cor)

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

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

lon1 <- -30.5 #
lat1 <- 47.5 #
sst_golf<-asst[which(lon==lon1),which(lat==lat1),]

plot(yrs,sst_golf,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=-30.5, Lat= 47.5'))

cmatrix_golf <- matrix(NA,dim(lon),dim(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    cmatrix_golf[i,j] <- cor(asst[i,j,], sst_golf)
  }
}

grid$corr_golf <- as.vector(cmatrix_golf)

levelplot(corr_golf~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=-30.5,Lat= 47.5')) + 
              layer(sp.lines(world.coast))

Task 3: rather than using “cor” repeat using “lm”. Calculate the slope at each point and plot the pattern

avsst_lm = rowMeans(sst,na.rm=TRUE,dims=2)

grid_lm <- expand.grid(x=lon, y=lat)