, , , ,

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 asked them for several files with monthly averages of global solar radiation over the Iberian Peninsula (download). The raster package can easily read these files:

##change to your folder...
old <- setwd('CMSAF')
listFich <- dir(pattern='2008')
stackSIS <- stack(listFich)
stackSIS <- stackSIS*24 ##from irradiance (W/m2) to irradiation Wh/m2

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 <- seq(as.Date('2008-01-15'), as.Date('2008-12-15'), 'month')

SISmm <- setZ(stackSIS, idx)
layerNames(SISmm) <-

The rasterVis package provides a levelplot method for Raster objects:

levelplot(SISmm, layers=1, FUN.margin=median, contour=TRUE)

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…

proj <- CRS(projection(SISmm))
old <- setwd('ESP_adm')##Change to your folder
mapaSHP <- readShapeLines('ESP_adm2.shp', proj4string=proj)

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

levelplot(SISmm) + 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 <-[c('Gefd', 'Befd', 'Defd')]##the results are yearly values

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='CMSAF/gefCMSAF',##change to your folder
layerNames(gefS)=c('Gefd', 'Befd', 'Defd')##Three layers
levelplot(subset(gefS, 'Gefd')) + layer(sp.lines(mapaSHP))