GY667 Workshop 3

Load the Libraries

library(sp)
library(raster) 

## Warning: package 'raster' was built under R version 3.5.2

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.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1

library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(reshape2)
library(maps)

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

Setting the Working directory

setwd("~/Dropbox/GY667_Workshop3/data")
path <- file.path(getwd(),"data");

load("~/Dropbox/GY667_Workshop3/data/world.coast.rdata")

Loading in the NCEP data

HAdISST_sst.nc <- nc_open(file.path(getwd(),"HadISST_sst.nc"));

sst

sst <- ncvar_get(HAdISST_sst.nc,"sst"); 

time

time   <- ncvar_get(HAdISST_sst.nc,"time");

latitude

lat  <- ncvar_get(HAdISST_sst.nc,"latitude") 

longitude

lon  <- ncvar_get(HAdISST_sst.nc,"longitude") 

tunits<-ncatt_get(HAdISST_sst.nc,"time",attname="units")
tustr<-strsplit(tunits$value, " ")

checkl data time origin is okay

date<-as.character(as.Date(time,origin=unlist(tustr)[3]))

Replace the missing values

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

use annual averages with aggregate and using colMeans and rowMeans

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)

two lines of code to get nicer colours for the eventual plots

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

‘levelplot’

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

get lats and lons added to a data frame

grid <- expand.grid(x=lon, y=lat)

Addinbg sst as a vector

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

Changing for annual averages

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

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

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

removing globl mean for each year

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

Assigning longitude and latitude

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

looking at the timeseries

plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon0, ',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)
  }
}

add the same data frame to plot

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

changing lat and long

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

timeseries

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

adding the same data frame as before

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