Large scale weather correlations using R’s raster package

The previous post on this blog used historical wind speed data from ECMWF reanalysis. Reanalysis products are observation-based snapshots of past states of the atmosphere. They provide a physically consistent picture of the past using current models and the historical data.

For example, suppose you want to know the relation between spatially averaged Irish and Scottish wind speeds. The result from ECMWF reanalysis is shown below. Clearly Scottish and Irish wind speeds are correlated. This might have implications for trading of wind power generation between the UK and Ireland, for example.

This type of analysis is easy to do in R using Robert Hijman’s superb raster package. The first step is to download horizontal(U) and vertical(V) wind speed data from ECMWF Data Server for ERA-interim. To decode this file, the wgrib utility or similar must be installed on your system. The following code creates a raster “stack” object consisting of 1342 U-component wind speed raster maps for NW Europe: 
library(raster)
region <- extent(-40,40,0,80) #subset selection -40W 40E 0N 80N
Ucomp.stk <- stack()
for (i in seq(1,2*1342,by=2) ){
system(paste("wgrib ERAdailyWind.grb -d ", i, " -text -o temp.txt",sep=""))
temp <- scan("temp.txt",skip=1)
temp <- t(matrix(temp,240,121))
temp.ras <- rotate(raster(temp,xmn=0,xmx=360,ymn=-90,ymx=90))
temp.ras <- crop(temp.ras,region)
Ucomp.stk <- addLayer(Ucomp.stk,temp.ras)
}
....
wind.stk <- calc(calc(Ucomp.stk,function(x) x*x)+calc(Vcomp.stk,function(x) x*x),sqrt)

The last line calculates absolute wind speeds from Ucomp.stk and Vcomp.stk. rotate() optionally shifts to Greenwich centred maps. The code for Vcomp.stk is the same as Ucomp.stk and omitted.

A rasterized map of Scotland can be created from a level 1 GADM SpatialPolygons object as follows:
UK <- getData("GADM",country="GBR",level=1)
scotland <- SpatialPolygons(UK@polygons[3])
scotland.ras <- rasterize(scotland, wind.stk,getCover=T)
scotland.ras <- scotland.ras/cellStats(scotland.ras, sum)
The option getCover=T in rasterize() means that any reanalysis cells which partially overlap Scotland are included with appropriate weight. The average wind speeds for Scotland are then:

scottish.wind.stk <- scotland.ras * wind.stk
scottish.speeds <- cellStats(scottish.wind.stk, sum)


Irish average wind speeds are found in the same way.

Finally the high density scatterplot shown above was produced using the hexbin package:


library(hexbin)
wind.bin <- hexbin(irish.speeds,scottish.speeds)
count.cols <- colorRampPalette(c("azure2","yellow","orange","red"),space="Lab")
plot(bin, colramp=function(n) count.cols(n), main="12hr Wind Speeds 2009-2010",xlab="Ireland m/s",ylab="Scotland m/s",legend=0)