Here comes the sun

(This article was first published on 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.


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.