Day length calculation

December 21, 2013
By

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

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().

daylength graph

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))

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.

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)