Solstice sunrise-sunset
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Introduction
The sunAngle() function in the oce package provides a handy way to determine sunrise/sunset azimuth angles, and this is used to find alignments for the winter solstice in Halifax, NS.
Methods
First, set up the problem; these may be the only lines that need editing for other places or times.
1 2 3 | t <- as.POSIXct("2013-12-21 17:11:00", tz="UTC") # winter solstice
xy <- list(x=-63.60, y=44.65) # centre of map (Halifax)
D <- 6 # map span in km
|
Next, use sunAngle() from the oce package to find the time of the sunrise and the azimuth at that moment. Here, uniroot() is used to find the time when the altitude is zero (the sun is on the horizon), and we start searching a quarter-day before the time of the actual solstice.
1 2 3 4 5 6 7 | library(oce)
t0 <- t - 86400 / 4
sunrise <- uniroot(function(t)
sunAngle(t, lat=xy$y, lon=xy$x)$altitude,
interval=as.numeric(c(t0, t)))$root
sunrise <- numberAsPOSIXct(sunrise)
azimuth <- 90 - sunAngle(sunrise, lat=xy$y, lon=xy$x)$azimuth
|
The span D is given in kilometres, which we convert to degrees of latitude and longitude.
1 2 | D <- D / 111 # deg Dlon <- D / cos(xy$y * pi / 180) |
Now it is time to start with the mapping, which uses the OpenStreetMap package.
1 2 3 4 5 | library(OpenStreetMap)
map <- openmap(c(lat=xy$y+D/2, lon=xy$x-Dlon/2),
c(lat=xy$y-D/2, lon=xy$x+Dlon/2),
minNumTiles=9)
plot(map)
|
Now, it remains to draw the sunrise/sunset lines. The variables cx and cy are the centre points of the map, in map units (not latitude and longitude).
1 2 3 4 5 6 7 | cx <- mean(par('usr')[1:2])
cy <- mean(par('usr')[3:4])
d <- diff(par('usr')[3:4]) # scales as the map
for (o in d*seq(-1, 1, length.out=30)) {
lines(cx+c(-1,1)*d*cos(azimuth*pi/180),
cy+o+c(-1,1)*d*sin(azimuth*pi/180), col='red')
}
|
Finally, add a title so that plot is self-explanatory.
1 2 | mtext(sprintf("sunrise angles on day of winter solstice (%s)",
format(t, "%Y-%m-%d")), side=3, line=2, cex=1.5)
|
Results
The graph indicates the results; click it to zoom in.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.