Search code examples
rggplot2rnoaa

Trying to map a value for geom_vline, but is not plotting in the correct place on the x axis with ggplot in R


I am currently trying to generate NOAA tide prediction charts (x = datetime, y = water level) with the dawn/sunrise/dusk/sunset times as vertical lines along the x axis timeline.

The rnoaa package calls the data and gives me the prediction date times in POSIXct. The suncalc library provides me a data frame with each date in the range's sunrise, sunset, etc. in POSIXct format as well.

library(rnoaa)
library(tidyverse)
library(ggplot2)
library(suncalc)
march.tides <- as.data.frame(coops_search(station_name = 8551762, 
                                          begin_date = 20200301, end_date = 20200331, 
                                          datum = "mtl", product = "predictions"))
march.tides <- march.tides %>%
  mutate(yyyy.mm.dd = as.Date(predictions.t))
dates <- unique(march.tides$yyyy.mm.dd)
sunlight.times <- getSunlightTimes(date = seq.Date(as.Date("2020/3/1"), as.Date("2020/3/31"), by = 1), 
                lat = 39.5817, lon = -75.5883, tz = "EST")

I then have a loop that spits out separate plots for each calendar date - which works hunky dory. The vertical lines are drawing on the graph without an error, but are definitely in the wrong spot (sunrise is being drawn around 11am when it should be 06:30).

for (i in seq_along(dates)) {
  plot <- ggplot(subset(march.tides, march.tides$yyyy.mm.dd==dates[i])) +
    aes(x = predictions.t, y = predictions.v) +
    geom_line(size = 1L, colour = "#0c4c8a") +
    theme_bw() + 
    geom_vline(xintercept = sunlight.times$sunrise) +
    geom_vline(xintercept = sunlight.times$sunset) +
    geom_vline(xintercept = sunlight.times$dawn, linetype="dotted") +
    geom_vline(xintercept = sunlight.times$dusk, linetype="dotted") +
    ggtitle(dates[i])
  print(plot)
}

I could alternatively facet the separate dates instead of this looping approach. Even when I subset the data to a single date, the vertical lines still did not draw correctly.

I wondered if maybe the issue was a time zone one. If I try to stick a time zone argument onto the tide prediction data call, I get the error:

Error in if (!repeated && grepl("%[[:xdigit:]]{2}", URL, useBytes = TRUE)) return(URL) : 
  missing value where TRUE/FALSE needed

Solution

  • It looks like you want to use EST as your timezone, so you could include in your conversion of predictions.t.

    I would be explicit in what you want labeled on your xaxis in ggplot using scale_x_datetime, including the timezone.

    library(rnoaa)
    library(tidyverse)
    library(ggplot2)
    library(suncalc)
    library(scales)
    
    march.tides <- as.data.frame(coops_search(station_name = 8551762, 
                                              begin_date = 20200301, end_date = 20200331, 
                                              datum = "mtl", product = "predictions"))
    
    march.tides <- march.tides %>%
      mutate(yyyy.mm.dd = as.Date(predictions.t, tz = "EST"))
    dates <- unique(march.tides$yyyy.mm.dd)
    sunlight.times <- getSunlightTimes(date = seq.Date(as.Date("2020/3/1"), as.Date("2020/3/31"), by = 1), 
                                       lat = 39.5817, lon = -75.5883, tz = "EST")
    
    for (i in seq_along(dates)) {
      plot <- ggplot(subset(march.tides, march.tides$yyyy.mm.dd==dates[i])) +
        aes(x = predictions.t, y = predictions.v) +
        geom_line(size = 1L, colour = "#0c4c8a") +
        theme_bw() + 
        geom_vline(xintercept = sunlight.times$sunrise) +
        geom_vline(xintercept = sunlight.times$sunset) +
        geom_vline(xintercept = sunlight.times$dawn, linetype="dotted") +
        geom_vline(xintercept = sunlight.times$dusk, linetype="dotted") +
        ggtitle(dates[i]) +
        scale_x_datetime(labels = date_format("%b %d %H:%M", tz = "EST"))
      print(plot)
    }
    

    Plot

    tide plot with xaxis date scale with timezone