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@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.
