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