#' @title Day Length Calculator
#' @author Mark Otto (FWS)
#'
#' @description This function calculates the annual daylight hours at a given location.
#' @param Lng Longitude at given location
#' @param Lat Latitude at given location
#' @param Type Annual, monthly, or custom. Defaults to annual
#' @param Labels Optional abels for custom time periods
#' @import maptools
#' @export
#' @examples
#' DayLen(Lng=-65.6342,Lat=38.118352,Type="Annual")
#' DayLen(Lng=-65.6342,Lat=38.118352,Type="Monthly")
#' DayLen(Lng=-75.6342,Lat=38.118352, Type=c("02/29","06/30","09/30"), Labels=c("Fall-Winter","Spring:Breed,Fledge","Summer"))
#' DayLen(Lng=-75.6342,Lat=38.118352,Type=c("02/29","06/30","09/30"), Labels=c("Winter","Spring:Breed,Fledge","Summer","Fall"))
# written by Mark Otto (FWS).
# Definitions for Daylight Hours, Sunrise, and Sunset are from NOAA
DayLen<-function(Lng,Lat,Type="Annual",Labels=NULL){
# Choose a generic leap year to include February 29.
Yr<-2012
# Make a table of days of the year
MoLen<-c(31,28,31,30,31,30,31,31,30,31,30,31)
if(Yr%%4==0)MoLen[2]<-29
DateYr<-data.frame(Yr,Mo=rep(1:12,times=MoLen),
Day=unlist(lapply(MoLen,function(MoLen)1:MoLen))
)
DateYr$Date<-with(DateYr,paste(Yr,Mo,Day,sep="-"))
rownames(DateYr)<-DateYr$Date
# Set the location. Longitude matters. Is it elevation?
coords<-data.frame(
xcoord=Lng,ycoord=Lat
)
# Sunrise, sunset
DateYr$Sunrise<-24*sunriset(crds=SpatialPoints(coords = coords,
proj4string=CRS("+proj=longlat +datum=WGS84")),
dateTime=as.POSIXct(DateYr$Date),# , tz="EET"),
direction="sunrise", POSIXct.out=FALSE
)
DateYr$Sunset<-24*sunriset(crds=SpatialPoints(coords = coords,
proj4string=CRS("+proj=longlat +datum=WGS84")),
dateTime=as.POSIXct(DateYr$Date),# , tz="EET"),
direction="sunset", POSIXct.out=FALSE
)
DateYr$DayLtHr<-with(DateYr,Sunset-Sunrise)
# Leapyear only occurs every four years.
DateYr[paste(Yr,2,29,sep="-"),"DayLtHr"]<-DateYr[paste(Yr,2,29,sep="-"),"DayLtHr"]/4
DateYr$Days<-1
DateYr[paste(Yr,2,29,sep="-"),"Days"]<-0.25
# Annual
if(Type[1]=="Annual"){
DateYr$Factor<-factor("Overall")
BELabel<-c(Overall="1/1-12/31")
# Monthly
} else if(Type[1]=="Monthly"){
cMo<-c(
"Jan","Feb","Mar","Apr","May","Jun",
"Jul","Aug","Sep","Oct","Nov","Dec"
)
DateYr$Factor<-factor(cMo[DateYr$Mo],levels=cMo)
tMoLen<-MoLen
tMoLen[2]<-28.25
BELabel<-paste(1:12,"/1-",1:12,"/",tMoLen,sep="")
names(BELabel)<-cMo
# To define the cut dates by their upper date, include
# the last day of last year as the first cut point and set
# right=T. format the dates so they can be sorted
} else {
Factor<-Type
tFac<-strsplit(Factor,"/")
tFac<-data.frame(t(data.frame(tFac)),stringsAsFactors=FALSE)
tFac<-apply(tFac,2,function(x)sub(" ","0",format(as.numeric(x))))
cCut<-c(
paste(Yr-1,12,31,sep="-"),
apply(tFac,1,function(iMoDay,Yr)
paste(Yr,iMoDay[1],iMoDay[2],sep="-"),Yr=Yr
),
paste(Yr,12,31,sep="-")
)
# Just do this to check then rerun if OK. Where is is.POSIXlt?
tCut<-sapply(cCut,function(iDate)try(as.POSIXct(iDate),TRUE))
# Try did not try hard enough
# sapply(tCut,function(iDate)!inherits(iDate,"try-error"))
ErrIdx<-grep("Error",tCut)
if(length(ErrIdx)>0){
stop(paste("Invalid cut point dates:",
paste(Factor[ErrIdx-1],collapse=",")
))
} else {
tCut<-as.POSIXct(cCut)
}
# Beginning-end labels. Wrap around if the number of labels=number of cuts.
Beg<-as.POSIXct(tCut+24*60*60)
nFactor<-length(tCut)
BELabel<-if(length(Labels)==length(Factor)){
data.frame(Beg=Beg[c(nFactor-1,2:(nFactor-2))],End=tCut[2:(nFactor-1)])
} else {
data.frame(Beg=Beg[-nFactor],End=tCut[-1])
}
rownames(BELabel)<-Labels
BELabel<-apply(BELabel,1,function(iDate){
cDate<-sub("201[123]-","",as.character(iDate))
cDate<-sub(" EST","",cDate)
cDate<-sub("-","/",cDate)
paste(cDate,collapse="-")
})
# If the number of labels is the same as the number cut points,
# the first cut point wraps around to include the end of the year.
tLabels<-if(is.null(Labels)){
1:(length(Factor)+1)
} else {
c(Labels,if(length(Labels)==length(Factor))Labels[1] else c())
}
# Make the cuts, return integers for the index labels to handle the
# wrap around labels.
Idx<-cut(as.POSIXct(DateYr$Date),breaks=tCut,labels=FALSE,right=TRUE)
# Associate the labels with the index and make a factor
DateYr$Factor<-factor(as.character(tLabels[Idx]),
levels=names(BELabel)
)
}
# Summarize the day length and number of days by the factor requested.
DaysLen<-aggregate(cbind(Days,DayLtHr)~Factor,data=DateYr,FUN=sum)
DaysLen$cRange<-BELabel
return(DaysLen)
}
'
DayLen(Lng=-65.6342,Lat=38.118352,Type="Annual")
DayLen(Lng=-65.6342,Lat=38.118352,Type="Monthly")
DayLen(Lng=-75.6342,Lat=38.118352,
Type=c("02/29","06/30","09/30"),
Labels=c("Fall-Winter","Spring:Breed,Fledge","Summer")
)
DayLen(Lng=-75.6342,Lat=38.118352,
Type=c("02/29","06/30","09/30"),
Labels=c("Winter","Spring:Breed,Fledge","Summer","Fall")
)
'
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.