1. Why Dates Are Hard
Before diving into code, it helps to understand why dates
cause so much trouble.
The representation problem. “2024-03-15” and
“03/15/24” and “15 March 2024” and “2024-075” (day-of-year notation) all
refer to the same day, but R will interpret each of them differently —
or fail silently and produce NA.
The storage problem. Computers store dates as
numbers (seconds since a reference epoch, or days since a reference
date), but display them as human-readable strings. Converting
between storage and display requires knowing the reference point, the
time zone, and whether daylight saving is in effect.
The sensor problem. Loggers deployed in the field
accumulate clock drift — a few seconds per day that compounds over
months. An iButton that was accurate on deployment day may be 45 minutes
off at retrieval. Merging drifted data with reference data produces
time-offset spikes that look like ecological signal but are
artifacts.
⚠️ Silent NA production is the most dangerous date
bug. as.Date("03/15/2024") returns NA
without any error. If you don’t check for NAs after parsing, you will
silently lose data — and any downstream filter(),
group_by(), or ggplot() call will simply skip
those rows without telling you.
2. Date Types in R
library(tidyverse)
library(lubridate)
R has three built-in date-time classes, and lubridate works with all
of them.
Date
Calendar date only — no time component. Stored as integer days since
1970-01-01.
class(Sys.Date()) # “Date”
POSIXct
Date + time. Stored as seconds since 1970-01-01 00:00:00 UTC. The
workhorse class for sensor data.
class(Sys.time()) # “POSIXct” “POSIXt”
POSIXlt
Date + time stored as a named list (year, mon, mday, …). Flexible for
extraction but slow and list-like — avoid in data frames.
class(strptime(“2024-01-01”,
“%Y-%m-%d”))
# Three ways to represent the same moment
d <- as.Date("2024-03-15") # date only
ct <- as.POSIXct("2024-03-15 09:30:00", tz = "America/Los_Angeles")
lt <- as.POSIXlt("2024-03-15 09:30:00", tz = "America/Los_Angeles")
# What does each store underneath?
unclass(d) # integer: days since 1970-01-01
## [1] 19797
unclass(ct) # numeric: seconds since 1970-01-01 UTC
## [1] 1710520200
## attr(,"tzone")
## [1] "America/Los_Angeles"
str(unclass(lt)) # named list of components
## List of 11
## $ sec : num 0
## $ min : int 30
## $ hour : int 9
## $ mday : int 15
## $ mon : int 2
## $ year : int 124
## $ wday : int 5
## $ yday : int 74
## $ isdst : int 1
## $ zone : chr "PDT"
## $ gmtoff: int NA
## - attr(*, "tzone")= chr "America/Los_Angeles"
## - attr(*, "balanced")= logi TRUE
📌 Use POSIXct in data frames, Date for daily data.
POSIXlt objects are lists — they cannot be stored as a
column in a data frame (tibble). R will silently convert them to
POSIXct if you try. For sensor data with sub-daily
timestamps, always use POSIXct. For data where you only
care about the day (daily weather, survey dates, publication dates),
Date is simpler and uses less memory.
The storage number in practice
# The number underneath tells you the distance from the epoch
as.integer(as.Date("1970-01-01")) # 0 — the epoch
## [1] 0
as.integer(as.Date("2024-03-15")) # days since epoch
## [1] 19797
as.numeric(as.POSIXct("1970-01-01 00:00:00", tz = "UTC")) # 0 — the epoch
## [1] 0
as.numeric(as.POSIXct("2024-03-15 09:30:00", tz = "UTC")) # seconds since epoch
## [1] 1710495000
Understanding this is not just trivia — it explains why
date1 - date2 returns a number of days, why
format() needs a time zone, and why adding
86400 to a POSIXct advances it by exactly one day.
3. Parsing Date Strings
Parsing is converting a character string like
"03/15/2024 09:30" into a Date or
POSIXct object. This is where most date headaches
originate.
3.1 lubridate helper functions
lubridate provides a family of functions named after the order of
date components in the input. These cover the vast majority of
real-world date formats.
ymd()
year-month-day, any separator
ymd(“2024-03-15”)
mdy()
month-day-year (US format)
mdy(“03/15/2024”)
dmy()
day-month-year (European format)
dmy(“15-03-2024”)
ymd_hms()
year-month-day hour:min:sec → POSIXct
ymd_hms(“2024-03-15 09:30:00”)
mdy_hm()
month-day-year hour:min → POSIXct
mdy_hm(“03/15/2024 09:30”)
parse_date_time()
Handles multiple formats in one vector — use when your column mixes
formats.
parse_date_time(x, c(“mdy”, “ymd”))
# These all produce the same Date object
ymd("2024-03-15")
## [1] "2024-03-15"
ymd("20240315") # no separator at all — still works
## [1] "2024-03-15"
mdy("03/15/2024")
## [1] "2024-03-15"
mdy("March 15, 2024") # lubridate understands month names
## [1] "2024-03-15"
dmy("15-Mar-24")
## [1] "2024-03-15"
# Adding a time component → POSIXct (UTC by default)
ymd_hms("2024-03-15 09:30:00")
## [1] "2024-03-15 09:30:00 UTC"
ymd_hms("2024-03-15 09:30:00", tz = "America/Los_Angeles")
## [1] "2024-03-15 09:30:00 PDT"
# A realistic iButton timestamp format: "03/15/24 09:30:00 AM"
mdy_hms("03/15/24 09:30:00 AM", tz = "America/Los_Angeles")
## [1] "2024-03-15 09:30:00 PDT"
💡 Always specify tz when parsing sensor
data. ymd_hms("2024-03-15 09:30:00") returns a UTC
timestamp — even if the sensor was running on Pacific time. If you
forget tz, every downstream plot and join will be off by
7–8 hours without any visible error.
3.3 Diagnosing parse failures
When parsing fails, you get NA. Always check after
parsing.
raw_dates <- c("2024-01-15", "January 20", "not-a-date", "2024-03-05")
parsed <- ymd(raw_dates, quiet = TRUE) # quiet suppresses the warning
# How many failed?
sum(is.na(parsed))
## [1] 2
# Which ones?
raw_dates[is.na(parsed)]
## [1] "January 20" "not-a-date"
⚠️ ymd() without a year produces NA.
"January 20" (no year) fails silently. This is a common
issue when working with historical climate data files that store month
and day in separate columns — join them first with paste()
before parsing.
3.4 Parsing from separate columns
Sensor data and weather station files often store year, month, day,
and sometimes hour as separate columns. Use make_date() and
make_datetime() to combine them.
# Typical weather station export format
weather_raw <- tibble(
YEAR = c(2023L, 2023L, 2023L, 2024L),
MONTH = c(11L, 11L, 12L, 1L),
DAY = c(28L, 29L, 1L, 15L),
HOUR = c(8L, 14L, 0L, 9L),
temp_C = c(12.4, 15.1, 8.9, 11.2)
)
weather_raw %>%
mutate(
date = make_date(YEAR, MONTH, DAY),
datetime = make_datetime(YEAR, MONTH, DAY, HOUR,
tz = "America/Los_Angeles")
)
## # A tibble: 4 × 7
## YEAR MONTH DAY HOUR temp_C date datetime
## <int> <int> <int> <int> <dbl> <date> <dttm>
## 1 2023 11 28 8 12.4 2023-11-28 2023-11-28 08:00:00
## 2 2023 11 29 14 15.1 2023-11-29 2023-11-29 14:00:00
## 3 2023 12 1 0 8.9 2023-12-01 2023-12-01 00:00:00
## 4 2024 1 15 9 11.2 2024-01-15 2024-01-15 09:00:00
5. Rounding and Flooring Timestamps
floor_date(), ceiling_date(), and
round_date() snap a timestamp to the nearest unit boundary.
These are essential for binning irregular sensor data and for joining
datasets with different temporal resolutions.
floor_date()
Round down to the nearest unit. Use for binning: “which 15-min
window does this reading belong to?”
floor_date(ts, “15 minutes”)
ceiling_date()
Round up to the nearest unit boundary.
ceiling_date(ts, “hour”)
round_date()
Round to the nearest unit (up or down). Standard rounding rules
apply.
round_date(ts, “30 minutes”)
# A realistic set of irregular sensor readings
irregular_ts <- ymd_hms(c(
"2024-03-15 09:02:47",
"2024-03-15 09:17:12",
"2024-03-15 09:31:55",
"2024-03-15 09:44:03",
"2024-03-15 10:01:29"
), tz = "America/Los_Angeles")
# Snap to 15-minute bins (floor = assign to start of the bin)
floor_date(irregular_ts, "15 minutes")
## [1] "2024-03-15 09:00:00 PDT" "2024-03-15 09:15:00 PDT"
## [3] "2024-03-15 09:30:00 PDT" "2024-03-15 09:30:00 PDT"
## [5] "2024-03-15 10:00:00 PDT"
# Snap to the hour
floor_date(irregular_ts, "hour")
## [1] "2024-03-15 09:00:00 PDT" "2024-03-15 09:00:00 PDT"
## [3] "2024-03-15 09:00:00 PDT" "2024-03-15 09:00:00 PDT"
## [5] "2024-03-15 10:00:00 PDT"
# Snap to the day
floor_date(irregular_ts, "day")
## [1] "2024-03-15 PDT" "2024-03-15 PDT" "2024-03-15 PDT" "2024-03-15 PDT"
## [5] "2024-03-15 PDT"
Aggregating to regular intervals
The most common workflow: sensor readings arrive at irregular or
sub-minute intervals; you need hourly or daily summaries for
analysis.
# Compute hourly mean PM2.5 from the simulated data
pm_hourly_mean <- pm_hourly %>%
mutate(hour_bin = floor_date(datetime, "hour")) %>%
group_by(site, hour_bin) %>%
summarise(
mean_pm25 = mean(pm25, na.rm = TRUE),
n_obs = n(),
.groups = "drop"
)
# Compute daily mean — same pattern, different unit
pm_daily <- pm_hourly %>%
mutate(day_bin = floor_date(datetime, "day")) %>%
group_by(site, day_bin) %>%
summarise(
mean_pm25 = mean(pm25, na.rm = TRUE),
max_pm25 = max(pm25, na.rm = TRUE),
.groups = "drop"
)
head(pm_daily)
## # A tibble: 6 × 4
## site day_bin mean_pm25 max_pm25
## <chr> <dttm> <dbl> <dbl>
## 1 Claremont 2024-01-01 00:00:00 12.5 29.2
## 2 Claremont 2024-01-02 00:00:00 9.93 17.2
## 3 Claremont 2024-01-03 00:00:00 14.9 23.9
## 4 Claremont 2024-01-04 00:00:00 12.2 19.7
## 5 Claremont 2024-01-05 00:00:00 14.2 26.2
## 6 Claremont 2024-01-06 00:00:00 9.59 22.7
💡 floor_date() is how you join datasets at
different resolutions. If you have 5-minute sensor data and
hourly meteorological data,
mutate(hour_bin = floor_date(datetime, "hour")) on the
sensor data creates a join key that matches the meteorological data’s
timestamp. This avoids the imprecise “nearest timestamp” joins that
introduce subtle alignment errors.
Valid unit strings for floor / ceiling / round
# Any of these work as the unit argument:
"second", "minute", "hour", "day", "week", "month",
"bimonth", "quarter", "season", "halfyear", "year"
# With a multiplier (lubridate >= 1.9):
"15 minutes", "6 hours", "2 weeks", "3 months"
6. Date Arithmetic
6.1 Durations vs. Periods
lubridate distinguishes two ways to add time:
- A Duration is an exact number of seconds.
ddays(1) always means 86 400 seconds.
- A Period is a human-readable unit that respects the
calendar.
days(1) means “tomorrow” — which may be 23, 24,
or 25 hours depending on DST transitions.
dst_day <- ymd_hms("2024-03-10 01:30:00", tz = "America/Los_Angeles")
# March 10, 2024 is the spring-forward DST transition in California
# Adding a Duration: exact seconds (ignores DST)
dst_day + ddays(1) # 25 hours later (crosses the 2-hour skip)
## [1] "2024-03-11 02:30:00 PDT"
# Adding a Period: calendar day (respects DST)
dst_day + days(1) # same clock time the next day
## [1] "2024-03-11 01:30:00 PDT"
⚠️ Use Periods for human-time questions, Durations for
physics. “How long did the sensor run?” → Duration (exact
seconds). “What is tomorrow’s date?” → Period. Mixing them up produces
off-by-one-hour errors in DST-affected time series — a very common bug
in spring air quality data.
6.2 Calculating intervals
# Sensor deployment and retrieval timestamps
deploy <- ymd_hms("2024-01-15 08:00:00", tz = "America/Los_Angeles")
retrieve <- ymd_hms("2024-04-02 16:30:00", tz = "America/Los_Angeles")
# Simple difference (returns a difftime object)
retrieve - deploy
## Time difference of 78.3125 days
# lubridate interval → precise duration
int <- interval(deploy, retrieve)
int_length <- as.duration(int)
int_length
## [1] "6766200s (~11.19 weeks)"
# Convert to convenient units
int_length / ddays(1) # total days (decimal)
## [1] 78.3125
int_length / dweeks(1) # total weeks
## [1] 11.1875
int_length / dhours(1) # total hours
## [1] 1879.5
6.3 Generating regular date sequences
# Daily sequence for a monitoring period
daily_dates <- seq(
from = ymd("2024-01-01"),
to = ymd("2024-03-31"),
by = "1 day"
)
length(daily_dates) # 91 days
## [1] 91
# Hourly sequence
hourly_ts <- seq(
from = ymd_hms("2024-01-01 00:00:00", tz = "America/Los_Angeles"),
to = ymd_hms("2024-01-07 23:00:00", tz = "America/Los_Angeles"),
by = "1 hour"
)
length(hourly_ts) # 168 hours (7 days × 24)
## [1] 168
7. Time Zones
Time zones are the source of more environmental data errors than any
other single issue. Understanding them properly takes about fifteen
minutes; not understanding them costs hours of debugging per
project.
7.1 R’s internal representation
R stores POSIXct timestamps as UTC
internally, always. The tz attribute only affects
how the number is displayed and how arithmetic handles DST.
This has a critical implication: two timestamps that look different on
screen may be identical underneath.
# The same physical moment in two time zones
utc_ts <- ymd_hms("2024-03-15 17:30:00", tz = "UTC")
pst_ts <- with_tz(utc_ts, "America/Los_Angeles")
cat("UTC display: ", format(utc_ts), "\n")
## UTC display: 2024-03-15 17:30:00
cat("PST display: ", format(pst_ts), "\n")
## PST display: 2024-03-15 10:30:00
cat("Same underlying number?", as.numeric(utc_ts) == as.numeric(pst_ts), "\n")
## Same underlying number? TRUE
7.2 with_tz() vs. force_tz()
These two functions look similar but do opposite things.
with_tz() — change the clock on the wall
The underlying moment is the SAME.
Only the display changes.
Use when: you have a correctly stored UTC timestamp and want to display
it in local time.
with_tz(ts,
“America/Los_Angeles”)
→ same instant, Pacific
display
force_tz() — re-label the timestamp
The underlying moment CHANGES.
The wall-clock reading stays the
same.
Use when: your data was recorded in local time but
incorrectly stored as UTC (common sensor bug).
force_tz(ts, “America/Los_Angeles”)
→ same
wall-clock, Pacific storage
# Scenario: iButton stored readings in local time (Pacific),
# but R imported them as UTC. The reading says "09:30" but means
# "09:30 Pacific" not "09:30 UTC" (which would be 01:30 Pacific).
wrong_utc <- ymd_hms("2024-03-15 09:30:00", tz = "UTC") # incorrectly labelled
# with_tz: converts to Pacific — now shows 01:30 (wrong — the sensor really read 09:30)
with_tz(wrong_utc, "America/Los_Angeles")
## [1] "2024-03-15 02:30:00 PDT"
# force_tz: re-labels as Pacific without changing the wall-clock time — shows 09:30 (correct)
corrected <- force_tz(wrong_utc, "America/Los_Angeles")
corrected
## [1] "2024-03-15 09:30:00 PDT"
7.3 Listing valid time zone names
# All valid tz strings on your system
OlsonNames()
# Useful California/Southwest ones
grep("America/Los_Angeles|America/Phoenix|America/Denver|UTC",
OlsonNames(), value = TRUE)
📌 “America/Los_Angeles” not “PST” or “PDT”. Never
use “PST”, “PDT”, “EST”, or “EDT” as time zone strings in R. They are
ambiguous abbreviations, not IANA time zone identifiers, and R may
silently misinterpret or ignore them. Always use the full IANA name from
OlsonNames().
8. Sensor Timestamp Correction
Real-world sensors — iButtons, Raspberry Pis without a
GPS-synchronized clock, Arduino-based loggers — accumulate clock
drift: a systematic offset that grows linearly over the
deployment period. A sensor set correctly at deployment will be several
minutes or hours off at retrieval, depending on the hardware and the
deployment length.
8.2 Diagnosing clock drift
The key diagnostic is to compare the sensor’s timestamps against
known reference events: the time you set the sensor and the time you
retrieved it.
# Recorded times at deployment and retrieval (from your field notebook)
sensor_at_deploy <- ibutton_raw$sensor_time[1]
sensor_at_retrieve <- ibutton_raw$sensor_time[n]
# True times from GPS phone or a reference clock
ref_deploy <- true_deploy
ref_retrieve <- true_retrieve
# Offset at each end
offset_start <- as.numeric(
sensor_at_deploy - ref_deploy, units = "secs")
offset_end <- as.numeric(
sensor_at_retrieve - ref_retrieve, units = "secs")
cat(sprintf("Clock offset at deployment : %+.0f seconds\n", offset_start))
## Clock offset at deployment : +0 seconds
cat(sprintf("Clock offset at retrieval : %+.0f seconds\n", offset_end))
## Clock offset at retrieval : +2820 seconds
cat(sprintf("Total drift over deployment: %+.0f seconds (%.1f min)\n",
offset_end - offset_start,
(offset_end - offset_start) / 60))
## Total drift over deployment: +2820 seconds (47.0 min)
# Visualise: plot sensor time vs true time for the two reference points
ref_check <- tibble(
event = c("Deployment", "Retrieval"),
true_time = as.numeric(c(ref_deploy, ref_retrieve)),
sensor_time = as.numeric(c(sensor_at_deploy, sensor_at_retrieve))
)
ggplot(ref_check, aes(x = true_time, y = sensor_time - true_time,
label = event)) +
geom_point(size = 4, color = "#c8832a") +
geom_line(color = "#c8832a", linewidth = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#4a7c5e") +
geom_text(vjust = -1, size = 3.5, color = "#1e3a2f", fontface = "bold") +
scale_x_continuous(labels = function(x)
format(as.POSIXct(x, origin = "1970-01-01", tz = "America/Los_Angeles"),
"%b %d")) +
labs(
title = "iButton clock drift diagnostic",
subtitle = "Points above zero = sensor running fast",
x = NULL, y = "Sensor offset (seconds)"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(family = "serif", face = "bold",
color = "#1e3a2f"))
8.3 Linear correction
Because clock drift is nearly linear, a simple linear interpolation
between the two known reference points corrects it precisely.
# Linear model: true_time = a + b * sensor_time
# We have two known (sensor, true) pairs
sensor_ref <- as.numeric(c(sensor_at_deploy, sensor_at_retrieve))
true_ref <- as.numeric(c(ref_deploy, ref_retrieve))
drift_model <- lm(true_ref ~ sensor_ref)
coef(drift_model) # intercept and slope
## (Intercept) sensor_ref
## 7.131037e+05 9.995818e-01
# Apply correction to the full dataset
ibutton_corrected <- ibutton_raw %>%
mutate(
sensor_numeric = as.numeric(sensor_time),
true_numeric = predict(drift_model,
newdata = data.frame(sensor_ref = sensor_numeric)),
corrected_time = as.POSIXct(true_numeric,
origin = "1970-01-01",
tz = "America/Los_Angeles")
) %>%
select(sensor_time, corrected_time, temp_C)
# Check: how close is the corrected time to the true time?
max_error_secs <- max(abs(
as.numeric(ibutton_corrected$corrected_time) - as.numeric(true_ts)
))
cat(sprintf("Maximum residual error after correction: %.2f seconds\n",
max_error_secs))
## Maximum residual error after correction: NA seconds
📌 Why linear correction works. Crystal oscillators
— the timekeeping element in most data loggers — drift at a constant
rate determined by temperature and manufacturing tolerance. The drift is
not perfectly constant (temperature changes affect it slightly), but the
linear model removes the dominant component. Residual errors after
linear correction are typically under a few seconds for a 3-month
deployment, which is negligible for hourly or sub-hourly analysis.
8.4 Before and after: checking the correction visually
# Compare sensor vs corrected for a subset in the middle of deployment
mid_window <- ibutton_corrected %>%
filter(between(corrected_time,
ymd_hms("2024-02-15 00:00:00", tz = "America/Los_Angeles"),
ymd_hms("2024-02-17 23:59:59", tz = "America/Los_Angeles")))
pivot_longer(
mid_window,
cols = c(sensor_time, corrected_time),
names_to = "version",
values_to = "timestamp"
) %>%
mutate(version = recode(version,
"sensor_time" = "Raw (drifted)",
"corrected_time" = "Corrected"
)) %>%
ggplot(aes(x = timestamp, y = temp_C, color = version)) +
geom_line(linewidth = 0.8, alpha = 0.9) +
scale_color_manual(
values = c("Raw (drifted)" = "#c8832a", "Corrected" = "#1e3a2f"),
name = NULL
) +
labs(
title = "Drift correction: raw vs. corrected timestamps",
subtitle = "Two-day window — offset visible as horizontal shift in the diurnal curve",
x = NULL, y = "Temperature (°C)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold",
color = "#1e3a2f"),
plot.subtitle = element_text(color = "#7a7a72", size = 10),
legend.position = "top"
)
9. Plotting Dates and Times with ggplot2
ggplot2 recognizes Date and POSIXct objects
and provides specialized scales for each.
9.1 scale_x_date() for daily data
pm_daily %>%
filter(site == "Upland") %>%
ggplot(aes(x = day_bin, y = mean_pm25)) +
geom_line(color = "#4a7c5e", linewidth = 0.7) +
geom_smooth(method = "loess", span = 0.2,
color = "#c8832a", fill = "#fdf8f0",
linewidth = 0.8, se = TRUE) +
# Control axis label frequency and format
scale_x_date(
date_breaks = "2 weeks",
date_labels = "%b %d" # "Jan 15", "Feb 01", etc.
) +
scale_y_continuous(expand = expansion(mult = c(0, 0.08))) +
labs(
title = "Daily mean PM2.5 — Upland, CA (simulated)",
subtitle = "Orange ribbon = LOESS trend with 95% CI",
x = NULL, y = expression(PM[2.5]~(mu*g~m^{-3}))
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold",
color = "#1e3a2f"),
plot.subtitle = element_text(color = "#7a7a72", size = 10),
axis.text.x = element_text(angle = 30, hjust = 1, size = 9)
)
9.2 scale_x_datetime() for sub-daily data
pm_hourly %>%
filter(
site == "Claremont",
between(datetime,
ymd_hms("2024-01-08 00:00:00", tz = "America/Los_Angeles"),
ymd_hms("2024-01-14 23:59:59", tz = "America/Los_Angeles"))
) %>%
ggplot(aes(x = datetime, y = pm25)) +
geom_line(color = "#2d5a3d", linewidth = 0.6, alpha = 0.75) +
scale_x_datetime(
date_breaks = "1 day",
date_labels = "%a\n%b %d" # "Mon\nJan 08" — two lines
) +
labs(
title = "Hourly PM2.5 — Claremont, one week (simulated)",
x = NULL, y = expression(PM[2.5]~(mu*g~m^{-3}))
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold",
color = "#1e3a2f")
)
9.4 Faceting by time period
A powerful pattern for seasonal or diurnal analysis: plot by one time
grouping, facet by another.
# Mean diurnal PM2.5 profile by month
pm_hourly %>%
group_by(site, month, hour) %>%
summarise(mean_pm25 = mean(pm25, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = hour, y = mean_pm25, color = site)) +
geom_line(linewidth = 0.9) +
geom_point(size = 1.8) +
facet_wrap(~ month, nrow = 1) +
scale_color_manual(
values = c("Upland" = "#1e3a2f", "Claremont" = "#c8832a"),
name = NULL
) +
scale_x_continuous(breaks = c(0, 6, 12, 18),
labels = c("midnight", "6am", "noon", "6pm")) +
labs(
title = "Mean diurnal PM2.5 by month — simulated Inland Empire data",
x = "Hour of day", y = expression(PM[2.5]~(mu*g~m^{-3}))
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(family = "serif", face = "bold",
color = "#1e3a2f", size = 12),
legend.position = "top",
axis.text.x = element_text(angle = 40, hjust = 1, size = 8),
strip.text = element_text(face = "bold")
)
