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
t0 <- as.POSIXct("2014-12-21 10:00:00", tz="UTC") # morning of 
xy <- list(x=-63.60, y=44.65)          # centre of map (Halifax)
D <- 9                                 # 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 use a search interval that should encompass sunrise at this longitude.

1
2
3
4
5
6
library(oce)
sunrise <- uniroot(function(t)
                   sunAngle(t, lat=xy$y, lon=xy$x, useRefraction=TRUE)$altitude,
                   interval=as.numeric(t0 + 3600*c(-5,5)))$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 show the results. The map is done with the OpenStreetMap package.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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)
## Draw lines along which sunrise or sunset can be sighted.
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=60)) {
    lines(cx+c(-1,1)*d*cos(azimuth*pi/180),
          cy+o+c(-1,1)*d*sin(azimuth*pi/180), col='red')
}
mtext(paste("Solstice sunrise at ", format(sunrise-4*3600, "%Y-%m-%d %H:%M")), font=2)

center

Exercises

Resources

This website is written in Jekyll, and the source is available on GitHub.