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
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