**Omnia sunt Communia! » R-english**, and kindly contributed to R-bloggers)

The solar irradiance incident outside the earth’s atmosphere is called the extraterrestial or extra-atmospheric irradiance. It is derived from the solar constant only with geometric equations. It can be easily calculated with the calcSol function of the solaR package. With this post I will show an example with some packages from the Spatial task view.

We will build a map of the extraterrestial irradiance for a certain day and time. First, we have to compute the mean solar time for a range of longitudes. The function local2Solar will do the task.

hh <- as.POSIXct('2011-05-01 11:00:00', tz='CET') latitude <- seq(70, -70, -1) longitude <- seq(-180, 180, 1) horaLong <- local2Solar(hh, longitude)

Then, we can compute the irradiance for the window defined by *latitude* and *longitude* with *lapply*.

`solList <- lapply(latitude, calcSol, BTi=horaLong)`

Now we have a list of *Sol* objects. Once again we use *lapply* to extract the information we really need.

Bo0List <- lapply(solList, function(x)as.data.frameI(x)$Bo0)

The result is a list of numeric elements, which is converted to a vector with *do.call* and *c*.

Bo0 <- do.call('c', Bo0List)

Finally, we assign the zero value to those *NA* elements, in order to get them black coloured in the map.

`Bo0[is.na(Bo0)] <- 0`

We are now prepared to use the *sp* package. Let’s build a SpatialPixelsDataFrame with our result.

Bo0DF=expand.grid(lon=longitude, lat=latitude) Bo0DF$Bo0 <- c(Bo0)proj <- CRS('+proj=latlon +ellps=WGS84') Bo0SP <- SpatialPixelsDataFrame(points=Bo0DF[,1:2], data=Bo0DF["Bo0"], proj4string=proj)

And we then get the map of extraterrestial irradiance:

paleta=colorRampPalette(rev(brewer.pal('Greys', n=9))) p <- spplot(Bo0SP, scales=list(draw=TRUE), col.regions=paleta, cuts=50)

The (almost) final step: let’s use the layer function of the latticeExtra to add a layer with the countries borders. (This task can also be carried out through the *sp.layout* of spplot)

world <- map("world", plot=FALSE) world_sp <- map2SpatialLines(world, proj4string=proj) p2 <- p+layer(sp.lines(world_sp, lwd=0.5))

This is the real final step. Where is located the maximum value of the irradiance? The geonames package will help to find it.

library(geonames) idxMax <- which.max(Bo0) lonlatMax <- Bo0DF[idxMax,1:2] place <- GNfindNearby(lat=lonlatMax[2], lng=lonlatMax[1]) place <- place[[1]][[1]]$countryName p3 <- p2 + layer(panel.text(lonlatMax[1], lonlatMax[2], place, cex=0.5)) p3

We got it!

The code of this post is available here.

**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...