R code used to create ‘levelplot’ of wind speed over the North Atlantic during AMO+ and AMO-.

Load relevant packages and set working directory

library(sp)
## Warning: package 'sp' was built under R version 3.4.3
library(raster)
## Warning: package 'raster' was built under R version 3.4.3
library(rasterVis)
## Warning: package 'rasterVis' was built under R version 3.4.4
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.4.3
## Loading required package: latticeExtra
## Warning: package 'latticeExtra' was built under R version 3.4.3
## Loading required package: RColorBrewer
library(ncdf4)
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.4.3
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/James/Documents/R/win-library/3.4/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/James/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-7
library(RColorBrewer)
library(lattice)
library(latticeExtra)
library(OpenStreetMap)
## Warning: package 'OpenStreetMap' was built under R version 3.4.4
library(reshape2)
## Warning: package 'reshape2' was built under R version 3.4.3
library(utils)
library(colorRamps)
library(maptools)
## Warning: package 'maptools' was built under R version 3.4.3
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
setwd("G:/Thesis data")
getwd()
## [1] "G:/Thesis data"
path <- file.path(getwd());

read coastline using readOGR as a Spatial Data Frame

world.countries <- readOGR(dsn = file.path(path,"ne_110m_admin_0_countries"), layer = "ne_110m_admin_0_countries" )
## OGR data source with driver: ESRI Shapefile 
## Source: "G:/Thesis data/ne_110m_admin_0_countries", layer: "ne_110m_admin_0_countries"
## with 177 features
## It has 94 fields
## Integer64 fields read as strings:  POP_EST NE_ID

convert to the correct projections

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

Open the u and v wind files in working directory

nc_u <-nc_open(file.path(path,'uwnd.sig995.mon.mean.nc'))
nc_v <-nc_open(file.path(path,'vwnd.sig995.mon.mean.nc'))

get the data, uwnd = monthly mean eastward wind

uwind <- ncvar_get(nc_u, 'uwnd')

get the data, vwnd = monthly mean northward wind

vwind <- ncvar_get(nc_v, 'vwnd')

Getting the dimensions of u and v wind datasets

dim(uwind) # = 180   91 1968
## [1]  180   91 1968
dim(vwind) # = 180   91 1968
## [1]  180   91 1968

Dimensions are the same in both arrays. 180 refers to longitudes (0-358 degrees in 2 degree increments) 91 refers to latitudes (-90-+90 degrees in 2 degree increments) 1968 refers to how many months are each in each array i.e 1850-2014 equals 164 years, multiplied by 12 months per year = 1,968 monts in total.

Replace missing values with NA

uwind[uwind==-9.99999] <- NA
vwind[vwind==-9.99999] <- NA
#^^Missing values in the v and u wind datasets are -9.99999, and so replaced with NA

Getting mean windspeed from u and v wind

windspeed <- sqrt(uwind^2 + vwind^2)
dim(windspeed) #180   91 1968
## [1]  180   91 1968

Extracting remaining variables; lat, lon and time

lat<-ncvar_get(nc_u, 'lat') # latitude, -90  90
lon<-ncvar_get(nc_v, 'lon') # logitude, 0 358
#^^ Note that lat and lon values can be extracted from either 'nc_u' or 'nc_v' as they share same dimensions.

load time

time<-ncvar_get(nc_u, 'time') # time
##^^ Same story with time, same dimensions in both u and v wind datasets, so no need to create two 
##separate values named 'time'.

get units

tunits<-ncatt_get(nc_u,"time",attname="units")
#^^This tells us that the time units used in both u and v wind datasets are as follows: 'hours since 01-01-1800 00:00:0.0'

tustr<-strsplit(tunits$value, " ")
date<-as.character(as.Date(time/24,origin=unlist(tustr)[3]))
#'date' now contains 1968 points (one for each month in the dataset) with each point being on the first day of each month.

windspeed is on a 0 to 360 grid. This can be seen by viewing lon, values go from 0 to 358 degrees.

lon_split = which(lon==180) # find the date line 
lon = c(lon[(lon_split+1):length(lon)]-360,lon[1:lon_split]) # rearrange the longitudes
windspeed = windspeed[c((lon_split+1):length(lon),1:lon_split),,] # and the data
#View(lon)
#View(lat)

Data Manipulation

Subsetting the data

AMO_pos_series <- windspeed[ , ,889:1380] #Monthly mean windspeed, 1925-1965
AMO_neg_series <- windspeed[ , ,1381:1740] #Monthly mean windspeed, 1966-1995

Getting annual means

Array of 1968 months

year <- format(as.Date(date, format="%Y-%m-%d"),"%Y")

Trim year to match extent of AMO positive and negative series.

year_pos <- year[889:1380]
year_neg <- year[1381:1740]

average the data

avwindspeed_AMOpos = rowMeans(AMO_pos_series,na.rm = FALSE,dims = 2)
avwindspeed_AMOneg = rowMeans(AMO_neg_series,na.rm = FALSE,dims = 2)

We’re not interested in short term variability so let’s make annual averages

yrs_pos <- unique(year_pos) 
yrs_neg <- unique(year_neg)

nyr_pos <- length(yrs_pos)
nyr_neg <- length(yrs_neg)

awindspeed will contain the annual windspeed

awindspeed_pos <- array(NA,c(length(lon),length(lat),nyr_pos))
awindspeed_neg <- array(NA,c(length(lon),length(lat),nyr_neg))

use a for loop to load the annual data into awindspeed

AMO positive phase first

for (k in 1:nyr_pos) {
  awindspeed_pos[,,k] <- rowMeans(AMO_pos_series[,,year_pos==yrs_pos[k]],na.rm=FALSE,dims=2) 
}

Then AMO negative phase

for (k in 1:nyr_neg) {
  awindspeed_neg[,,k] <- rowMeans(AMO_neg_series[,,year_neg==yrs_neg[k]],na.rm=FALSE,dims=2) 
}

‘expand.grid’ will get lats and lons added to a data frame expanded on each others size

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

our an_avwindspeed that we want to plot can be added as a vector

grid1$an_avwindspeed_pos <- as.vector(rowMeans(awindspeed_pos,na.rm = FALSE,dims = 2))
grid2$an_avwindspeed_neg <- as.vector(rowMeans(awindspeed_neg,na.rm = FALSE,dims = 2))

Plotting

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

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

AMO positive series (1925-1965)

Plot

levelplot(an_avwindspeed_pos~x*y, data=grid1,col.regions = pal(100),at=c(0:10),xlab='Longitude',
          ylab='Latitude', xlim=c(-120,10),ylim=c(0,80), 
          main='Annual Average Windspeed 1925-1965') + layer(sp.lines(world.countries))

AMO negative series (1966-1995)

Plot

levelplot(an_avwindspeed_neg~x*y, data=grid2,col.regions = pal(100),at=c(0:10),xlab='Longitude',
          ylab='Latitude', xlim=c(-120,10),ylim=c(0,80), 
          main='Annual Average Windspeed 1966-1995') + layer(sp.lines(world.countries))