# Day length calculation

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.

The winter solstice has been on many minds lately. The days are about to start getting longer … but just how fast will they do that?

This post provides R that calculates and graphs day length and its variation, using `uniroot()` to find sunrises and sunsets as indicated by solar altitude that is calculated with `sunAngle()` in the oce package.

The day of the solstice is indicated with vertical lines. All times are in UTC, which is the conventional system for scientific work and the one required by `sunAngle()`. The first step in making the graph shown above is to load the `oce` library and create a function that measures daylength by finding sunrise and sunset times. Note that `uniroot()`, which is used to find times of zero solar altitude, needs lower and upper limits on `t`, and these are calculated by looking back and then forward a half-day. This works well for application to Halifax, but in other timezones other offsets would be needed. Interested readers might want to devised a method based on the longitude, which can be transformed into a timezone.

 `````` 1 2 3 4 5 6 7 8 9 10`````` ```library(oce) daylength <- function(t, lon=-63.60, lat=44.65) { t <- as.numeric(t) alt <- function(t) sunAngle(t, longitude=lon, latitude=lat)\$altitude rise <- uniroot(alt, lower=t-86400/2, upper=t)\$root set <- uniroot(alt, lower=t, upper=t+86400/2)\$root set - rise } ```

Next, use `lappy()` to find the daylength for December, 2013.

 ``````1 2 3`````` ```t0 <- as.POSIXct("2013-12-01 12:00:00", tz="UTC") t <- seq.POSIXt(t0, by="1 day", length.out=1*31) dayLength <- unlist(lapply(t, daylength)) ```

Set up to plot two panels, with narrowed margins.

 ``1`` ```par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0)) ```

Plot daylength in the top panel, indicating the day of the solstice with vertical lines.

 ``````1 2 3 4 5`````` ```plot(t, dayLength/3600, type='o', pch=20, xlab="", ylab="Day length (hours)") grid() solstice <- as.POSIXct("2013-12-21", tz="UTC") abline(v=solstice+c(0, 86400)) ```

Plot daylength difference in the bottom panel, again indicating the solstice.

 ``````1 2 3 4`````` ```plot(t[-1], diff(dayLength), type='o', pch=20, xlab="Day in 2013", ylab="Seconds gained per day") grid() abline(v=solstice+c(0, 86400)) ```

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.