Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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, lng=lonlatMax)
place <- place[][]\$countryName
p3 <- p2 + layer(panel.text(lonlatMax, lonlatMax, place, cex=0.5))
p3```

We got it!        