Loading the Libraries

library(sp)
library(raster) 
library(ncdf4)
library(rgdal)
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

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

Loading Coastline Data

load("H:/University_Documents/Geography/MSc_Climate_Change/The Ocean and Climate Change/Gerard R Work/Workshop3/world.coast.Rdata") # downloaded from Moodle

# convert to the correct projections
world.coast <- spTransform(world.coast, CRS( "+proj=longlat +datum=WGS84" ) )

Task 1: Loading the NC File

nc   <- nc_open(file.path(getwd(),"HadISST_sst.nc"))
time <- ncvar_get(nc,"time"); # load the time
lat  <- ncvar_get(nc,"latitude") # load the latitude
lon  <- ncvar_get(nc,"longitude") # load the longitude
sst <- ncvar_get(nc, "sst")



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
# patterns of global variability
# remove the global mean at each step
# investigate the global mean
 
# using annual averages with aggregate and using colMeans and rowMeans
# alternatively using apply as in the first workshop
tunits<-ncatt_get(nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
date<-as.character(as.Date(time,origin=unlist(tustr)[3]))

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)
# alternatively commenting out the following two lines
# avsst = rowMeans(sst,na.rm=TRUE,dims=2)
avsst = rowMeans(sst,na.rm=FALSE,dims=2)

Mapping the NC Data

# 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
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, declared here to have the correct size, this is to speed up computation

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

SST Time Series (Ireland)

lon0 <- -10.5 #
lat0 <- 51.5 #
sst_ts <- 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))

SST Time Series (Sub-Polar Gyre)

lon0 <- -35.5 #
lat0 <- 51.5 #
sst_ts <- 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))

Correlation

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

# adding 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: Correlation map using a point north of the Gulf Stream

lon0 <- -60.5 #new coordinates
lat0 <- 45.5 #
sst_ts<-asst[which(lon==lon0),which(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))

Task 3: Repeat using “lm”. Calculate the slope at each point and plot the pattern

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


lmatrix <- matrix(NA,length(lon),length(lat))

for (i in 1:dim(lon)) {
  for (j in 1:dim(lat)) {
    if (length(asst[i,j,][!is.na(asst[i,j,])])>2){
    r.lm <- lm(asst[i,j,]~sst_ts)
    lmatrix[i,j] <- r.lm$coefficients[2]
    } } }


grid$lm <- as.vector(lmatrix)

levelplot(lm~x*y, data=grid , xlim=c(-120,10),ylim=c(0,80), at=c(-20:20)/10,
              col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with Long=',                 lon0, ',Lat=', lat0)) + layer(sp.lines(world.coast))