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)
Exercises
- Try altering
t0
to see how closely this alignment holds over time. - As above, but setting
t0
to an equinox. - Invent some sighting schemes for other times of the day, perhaps based on using local building (viewed at a distance computed from geometry) to check on noontime sun elevation.
Resources
- R source code used here: 2013-12-21-solstice.R.
- Jekyll source code for this blog entry: 2013-12-21-solstice.Rmd