R/DayLen.R

In josephstatwick/EagleModelCalc: Eagle fatality calculator for wind farms

#' @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")
  )
'