Raster, CMSAF and solaR

(This article was first published on Omnia sunt Communia! » R-english, and kindly contributed to R-bloggers)

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


To leave a comment for the author, please follow the link and comment on his blog: Omnia sunt Communia! » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , , ,

Comments are closed.