Want to share your content on Rbloggers? 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.

t < as.POSIXct("20131221 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 quarterday before the time of the actual solstice.

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.

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.

library(OpenStreetMap)
map < openmap(c(lat=xy$y+D/2, lon=xy$xDlon/2),
c(lat=xy$yD/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).

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

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.
Rbloggers.com offers daily email 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/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.