The Satellite Application Facility on Climate Monitoring (CMSAF) generates, archives and distributes widely recognised high-quality satellite-derived products and services relevant for climate monitoring in operational mode. The data is freely accesible here after a registration process.

I have ask them for several files with monthly averages of global solar radiation over the Iberian Peninsula (download). The raster package can easily read these files:

old <- getwd()
##change to your folder...
setwd('/home/oscar/Datos/CMSAF')
listFich <- dir(pattern='2008')
listNC <- lapply(listFich, raster)
stackSIS <- do.call(stack, listNC)
stackSIS <- stackSIS*24 ##from irradiance (W/m2) to irradiation Wh/m2
setwd(old)

It is interesting to know that with the last version of raster the Raster objects include a z slot, able to store information such as a time index. Its value is set with setZ:

idx <- fBTd('prom', year=2008)

SISmm <- setZ(stackSIS, idx)
layerNames(SISmm) <- as.character(idx)

Besides, raster has now an spplot method for Raster objects:

spplot(SISmm, names.attr=layerNames(SISmm))

Let’s improve this plot:

xscale.raster <- function(...){ans <- xscale.components.default(...); ans$top=FALSE; ans}
yscale.raster <- function(...){ans <- yscale.components.default(...); ans$right=FALSE; ans}

myTheme=custom.theme.2(pch=19, cex=0.7,
  region=rev(brewer.pal(9, 'YlOrRd')))
myTheme$strip.background$col='transparent'
myTheme$strip.shingle$col='transparent'
myTheme$strip.border$col='transparent'

p <- spplot(SISmm,
            names.attr=layerNames(SISmm),
            as.table=TRUE,
            par.settings=myTheme,
            between=list(x=0.5, y=0.2),
            scales=list(draw=TRUE),
            xscale.components=xscale.raster,
            yscale.components=yscale.raster)

And now let’s add the administrative borders. This information is available here. It is also accessible with raster::getData although it does not work for me…

library(maptools)
proj <- CRS('+proj=latlon +ellps=WGS84')
old <- getwd()
setwd('/home/oscar/Datos/ESP_adm')##Change to your folder
mapaSHP <- readShapeLines('ESP_adm2.shp', proj4string=proj)
setwd(old)

The layer function from latticeExtra and the sp.lines function from sp add easily this information to our previous plot:

p + layer(sp.lines(mapaSHP, lwd=0.8))

Now it is time for using solaR and calculate the effective irradiance in an inclined plane. First we have to define a suitable function:

foo <- function(x, ...){
  gef <- calcGef(lat=x[1], dataRad=list(G0dm=x[2:13]))
  result <- as.data.frameY(gef)[c('Gefd', 'Befd', 'Defd')]##the results are yearly values
  as.numeric(result)
}

raster::calc will apply this function to each cell of the stack. Previously we need a latitude layer:

latLayer <- init(SISmm, v='y')

Finally we are ready for the combination of raster and solaR. The next piece of code spends a long time: the result, the RasterLayer gefCMSAF, is available here).

gefS <- calc(stack(latLayer, SISmm), foo,
             filename='/home/oscar/Datos/CMSAF/gefCMSAF',##change to your folder
             overwrite=TRUE)
layerNames(gefS)=c('Gefd', 'Befd', 'Defd')##Three layers

spplot(subset(gefS, 'Gefd'), par.settings=myTheme, scales=list(draw=TRUE)) +  layer(sp.lines(mapaSHP))