daylight for tucson ... still needs work

daylight.R

library(maptools)
library(ggplot2)
library(gridExtra)
library(scales)
library(lubridate)

# adapted from http://r.789695.n4.nabble.com/maptools-sunrise-sunset-function-td874148.html
ephemeris <- function(lat, lon, date, span=1, tz="UTC") {
  
  lon.lat <- matrix(c(lon, lat), nrow=1)
  
  # using noon gets us around daylight saving time issues
  day <- as.POSIXct(sprintf("%s 12:00:00", date), tz=tz)
  sequence <- seq(from=day, length.out=span , by="days")
  #browser()
  sunrise <- sunriset(lon.lat, sequence, direction="sunrise", POSIXct.out=TRUE)
  sunset <- sunriset(lon.lat, sequence, direction="sunset", POSIXct.out=TRUE)
  solar_noon <- solarnoon(lon.lat, sequence, POSIXct.out=TRUE)
  
z <-   data.frame(date=as.Date(sunrise$time),
             sunrise=sunrise$time,
             solarnoon=solar_noon$time,
             sunset=sunset$time,
             day_length=as.numeric(sunset$time-sunrise$time)) 
  
}

# for graph #1 y-axis
time_format <- function(hrmn) substr(sprintf("%04d", hrmn),1,2)

# for graph #2 y-axis
pad5 <- function(num) sprintf("%2d", num)


daylight <- function(lat, lon, place, start_date, span=2, tz="UTC", show_solar_noon=FALSE, show_now=TRUE, plot=TRUE) {

  stopifnot(span>=2)
  srss <- ephemeris(lat, lon, start_date, span, tz)
  
  x_label = ""
  srss$sunrise_time <- hour(srss$sunrise) + minute(srss$sunrise)/60
  srss$sunset_time <- hour(srss$sunset) + minute(srss$sunset)/60
  ggplot(srss, aes(x=date)) + 
    geom_ribbon(aes(ymin= sunrise_time, ymax= sunset_time), fill = 'yellow', alpha = 0.7) +
    scale_x_date(expand=c(0,0), date_breaks = '1 month', date_labels = "%b") + 
    scale_y_continuous(limits = c(0,24), expand =c(0,0),breaks = 1:24) + 
    labs(x=x_label, y="time") + 
    theme_bw()+ 
    theme(panel.background=element_rect(fill="#525252"))+ 
    theme(panel.grid=element_line(color='darkblue', size = 0.1))
  
  
  if (plot) gg
  

}
daylight(lat = 32.22, lon = -110.97, 
         place = "Tucson", 
         start_date = "2020-01-01", 
         show_solar_noon=TRUE, 
         show_now=FALSE, 
         span = 365,
         tz="MST",
         plot=TRUE)