Solstice sunrise-sunset

December 21, 2013
By

(This article was first published on Dan Kelley Blog/R, and kindly contributed to R-bloggers)

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.

graph

To 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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)