Here comes the sun

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.