Solstice sunrise-sunset

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

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
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
D <- D / 111 # deg
Dlon <- D / cos(xy$y * pi / 180)
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)

Results

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.