[R] Days to solstice calculation

[Pascal Oettli]

Hello,

It seems that this kind of calculations are done in package 'insol'.

Regards,
Pascal


On 2 December 2013 15:26, White, William Patrick <white.232 at wright.edu> wrote:
> Hello,
> I've come across a problem in developing a set of custom functions to calculate the number of hours of daylight at a given latitude, and the number of days a date precedes or secedes the summer solstice. I discovered an inconsistency concerning leap years between my derived values and those from the US naval databases. It seems as far as I can figure that my inconsistency arises either in the calculation I used derived from an ecological modeling study in the 90's, in my understanding of the way R itself handles dates, or in my code. I feel like I must be missing something fundamental here and could use a little guidance. The first function returns the hours of daylight given a latitude, and the Julian day of the year (ie Jan 1 = 1 and so on). This appears to be very accurate. The second function takes a given date, extracts the year, determines the number of days in it, and uses the first function to calculate the hours of daylight in each day, and returns the longest or sh!
>  ortest one (Summer or Winter Solstice). But, in the case of leap years and non leap years, the date returned is identical, as is evidenced by Jan 1 in the provided examples being 170 days from summer solstice in both 2008 and 2007. This was not the case, the solstice should vary by one day between these years. Code is provided below and any help is appreciated.
> Patrick
> ps. apologies to you southern ducks your summer and winter solstices are reversed of my code nomenclature. I'm working with a northern dataset.
>
> Daylength <- function(J,L){
> #Amount of daylight
> #Ecological Modeling, volume 80 (1995) pp. 87-95, "A Model Comparison for Daylength as a Function of Latitude and Day of the Year."
> #D = Daylight length
> #L = Latitude in Degrees
> #J = Day of the year (Julian)
> P <- asin(.39795*cos(.2163108 + 2*atan(.9671396*tan(.00860*(J-186)))))
> A <- sin(0.8333*pi/180)+sin(L*pi/180)*sin(P)
> B <- cos(L*pi/180)*cos(P)
> D <- 24 - (24/pi)* acos(A/B)
> return(D)
> }
>
> #Example today and here
> Daylength(2,39.7505)
>
> TillSolstice <- function(date,solstice){
> Yr <- as.POSIXlt(date)$year+1900
> a <- as.Date(paste(as.character(Yr),as.character(rep("-01-01", length(Yr))),sep = ""))
> b <- as.Date(paste(as.character(Yr),as.character(rep("-12-31", length(Yr))),sep = ""))
> Winter <- NA
> Summer <- NA
> for (g in 1: length(a)){
> if(is.na(a[g])== FALSE){
> if(is.na(b[g])== FALSE){
>   cc <- seq.int(a[g],b[g], by = '1 day')
>   d <- length(cc)
>   e <- strptime(cc, "%Y-%m-%d")$yday+2
>   f <- Daylength(e,39.6981478)
>   Winter[g] <- which.min(f)
>   Summer[g] <- which.max(f)
> }
> }
> if(is.na(a[g])== TRUE){
>  Winter[g] <- NA
>   Summer[g] <- NA
> }
> if(is.na(b[g])== TRUE){
>  Winter[g] <- NA
>   Summer[g] <- NA
> }
>
>
> }
> #Days until solstice
> if (solstice =='S'){Countdown <- Summer - (strptime(date, "%Y-%m-%d")$yday+2)}
> if (solstice =='W'){Countdown <- Winter - (strptime(a, "%Y-%m-%d")$yday+2)}
> return(Countdown)
> }
>
> Nonleap <- TillSolstice(seq(as.Date("2007/1/1"), as.Date("2007/12/31"), by = "1 day"), solstice = 'S')
> Leap <- TillSolstice(seq(as.Date("2008/1/1"), as.Date("2008/12/31"), by = "1 day"), solstice = 'S')
> head(Nonleap)
> tail(Nonleap)
> length(Nonleap)
> head(Leap)
> tail(Leap)
> length(Leap)
>
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



-- 
Pascal Oettli
Project Scientist
JAMSTEC
Yokohama, Japan