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