GY667 Workshop 3

graphics.off()
rm(list=ls())
library(sp)
## Warning: package 'sp' was built under R version 3.5.3
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/Jordan/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/Jordan/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
path <- file.path(getwd(),"data")

Task 1 - Loading the ‘HadISST_sst.nc’ data

nc <- nc_open(file.path(getwd(),"data","HadISST_sst.nc"));
lat  <- ncvar_get(nc,"latitude") 
lon  <- ncvar_get(nc,"longitude") 
sst <- ncvar_get(nc,"sst")
date <- ncvar_get(nc,"time")
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

In order to have the date readable in R, the steps from Workshop 1 were used:

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

Below, the formatted time data could be used to extract a year, as well as global and annual means:

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)

Using rowMeans() without the NAs removed gives a very different picture of average SSTs.

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

Choosing the FALSE option as the TRUE option will bias area of sea ice, as true will be seasonally biased with regards to the presence of sea ice.

colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)
levelplot(avsst,col.regions = pal(100));

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

Making annual averages

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

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

Looking timeseries in locations = #1 Ireland

lon0 <- -10.5 #
lat0 <- 51.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))

Creating a correlation matrix

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

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

First, a point in the north atlantic above the gulf stream was identified: off the south-eastern coast of Iceland.

grid1 <- expand.grid(x=lon, y=lat)
grid1$avsst <- as.vector(avsst)
grid1$an_avsst <- as.vector(rowMeans(asst,na.rm=FALSE,dims=2))
lon1 <- -59.5
lat1 <- 45.5
sst_ts1<-asst[which(lon==lon1),which(lat==lat1),]
plot(yrs,sst_ts1,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon1, ',Lat=', lat1))

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

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    cmatrix1[i,j] <- cor(asst[i,j,], sst_ts1)
  }
}
grid1$corr <- as.vector(cmatrix1)
levelplot(corr~x*y, data=grid1, 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=', lon1, ',Lat=', lat1)) + 
layer(sp.lines(world.coast))

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

avsst1 = rowMeans(sst,na.rm=TRUE,dims=2) 
grid2 <- expand.grid(x=lon, y=lat)
grid2$avsst1 <- as.vector(avsst)
levelplot(avsst1~x*y,grid,col.regions = pal(100),
          xlab='Longitude',ylab='Latitude',main='Average SST'
) + 
  layer(sp.lines(world.coast))

Replicating the annual averages

yrs1 <- annmean$Group.1 
nyr1 <- length(yrs)