# Solstice sunrise-sunset

December 21, 2013
By

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.