# Solstice sunrise-sunset

**Dan Kelley Blog/R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

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.

**leave a comment**for the author, please follow the link and comment on their blog:

**Dan Kelley Blog/R**.

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.