Task 1

Before begining the task the appropiate work directory was loaded.

setwd("/Users/Kennf/Documents/GY667workshop3")

Libraries were loaded and graphics were turned off

library(sp)
library(raster) 
library(ncdf4)
library(rgdal)
## rgdal: version: 1.3-6, (SVN revision 773)
##  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/Kennf/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/Kennf/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 3.5.2
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)
graphics.off()
rm(list=ls())

loading the world.coast.

load("~/GY667workshop3/world.coast.Rdata")

Loading the NCEP data.

Using the following commands, the ncep data was loaded and so the data could be read.

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

once this is completed, the missing values can now be replaced.

fillvalue <- ncatt_get(HAdISST_sst.nc,"sst","_FillValue")
sst[sst==fillvalue$value] <- NA
missvalue <- ncatt_get(HAdISST_sst.nc,"sst","missing_value")
sst[sst==missvalue$value] <- NA
sst[sst==-1000] <- NA

Using the aggregate function to help get the row-means and the col-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)
avsst = rowMeans(sst,na.rm=FALSE,dims=2)

before the plots were to be created. Two lines ofcode were used to give apporaite colurs.

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

the level plot code was added to create the first plot.

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

the next step was adding the latitude and longitude.As well as adding SST as a vector to the data-frame

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

The next line of code was used to change for annual averages.

yrs <- annmean$Group.1 
nyr <- length(yrs)
asst <- array(NA,c(dim(lon),dim(lat),nyr))

next annual data has been added in the asst

for (k in 1:nyr) {
  asst[,,k] <- rowMeans(sst[,,year==yrs[k]],na.rm=FALSE,dims=2)
}

which is then added in to the same data frame so it can be plotted.

grid$an_avsst <- as.vector(rowMeans(asst,na.rm=FALSE,dims=2))

finally the plot can be loaded.

levelplot(an_avsst~x*y, data=grid,col.regions = pal(100),xlab='Longitude',ylab='Latitude',main='Annually Averaged SST') +
 layer(sp.lines(world.coast))

once plotted, the global mean from each year was removed with the following code.

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

the latitude and longitude were added

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

the following code allows the time seires to be viewd

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

The following lines of code are used to plot the second map

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

Task2

Doing the same as task 1 but with cordinates located north of the gulf stream

this was done with changing the latitude and longitude

lon1 <-  -25.5
lat1 <-  56.5
sst_ts<-asst[which(lon==lon1),which(lat==lat1),]

The time seires is then created

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

finally the plot is created

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