R Markdown

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:/R-3.5.1/library/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:/R-3.5.1/library/rgdal/proj
##  Linking to sp version: 1.3-1
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)
load("D:/Masters in Climate Change/Oceanography/workshop 3/world.coast (1).Rdata")

setting our working directory

setwd("D:/Masters in Climate Change/Oceanography/workshop 3")

Set the path to the data

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

This coastline data was downloaded from the natural earth data website.

Task 1

HadSST <- nc_open(file.path(getwd(),"HadISST_sst.nc"))
lat  <- ncvar_get(HadSST,"latitude")
lon  <- ncvar_get(HadSST,"longitude")
sst <- ncvar_get(HadSST, "sst")
time <- ncvar_get(HadSST, "time")
tunits<-ncatt_get(HadSST,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
date<-as.character(as.Date(time,origin=unlist(tustr)[3]))

Loaded the Hadley sea surface temperature data using the nc_open function. Some of the variables needed to be renamed. These were latitude, longitude, time and sea surface temp.

fillvalue <- ncatt_get(HadSST,"sst","_FillValue")
sst[sst==fillvalue$value] <- NA
missvalue <- ncatt_get(HadSST,"sst","missing_value")
sst[sst==missvalue$value] <- NA
sst[sst==-1000] <- NA
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)
avsst = rowMeans(sst,na.rm=FALSE,dims=2)
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))

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

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

Plotting a timeseries

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

Plotted above is sea surface temperature anomalies for long -10.5 and lat 51.5.

calculating the correlation

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

lon1 <- -44.5
lat1 <- 54.5 
sst_ts<-asst[which(lon==lon1),which(lat==lat1),]

New latitude and longitudes were chosen to graphh SSTA further north of the gulf stream. Longitude is -44.5 and latitude was 54.5. These were gathered from latlong.net.

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

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=', lon1, ',Lat=', lat1)) + 
  layer(sp.lines(world.coast))

Task 3

is.na)grid\(avsst) if(TRUE) skip(grid\)avsst > 2)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.