Load the relevant packages

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

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

## 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/Ledi/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/Ledi/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)

Set your working directory.

setwd("C:\\Users\\Ledi\\Documents\\Workshop 3")
getwd()

## [1] "C:/Users/Ledi/Documents/Workshop 3"

Set the path to the data:

path <- file.path(getwd(),"Workshop 3")

load("world.coast.RData")

Convert to the correct projections:

world.coast <- spTransform(world.coast, CRS("+proj=longlat +datum=WGS84"))

Task 1:

Load HadISST_sst.nc:

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

tunits<-ncatt_get(nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
Date<-as.character(as.Date(ts,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 comment out the following two lines, what difference does excluding NA values make? avsst = rowMeans(sst,na.rm=TRUE,dims=2)

avsst = rowMeans(sst,na.rm=FALSE,dims=2) # without sea ice

colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)

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

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) # annual averages from monthly data
}

The above code has created 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')

Removing the global mean from each year, using 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 for near Ireland:

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

The timeseries:

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

Calculating the correlation in a loop:

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)

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

Task 2: Correlation map of a point north of the Gulf Stream:

lon1 <- -35.5
lat1 <- 63.5
sst_ts1<-asst[which(lon==lon1),which(lat==lat1),]

plot(yrs,sst_ts1,type='l',xlab='Year',ylab='SST Anomaly North Gulf Stream',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)
  }
}

Add to the same data frame as earlier for plotting:

grid$corr1 <- as.vector(cmatrix1)

#plot
levelplot(corr1~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))

Workshop 3

Ledi Hoxha

2 April 2019

Task 1:

SST time series for near Ireland:

Task 2: Correlation map of a point north of the Gulf Stream: