Daily_Temperatures.Rmd

---
title: "Long-term Trends in Day-Of-Year High and Low Temperatures"
output: html_document
editor_options: 
  chunk_output_type: console
---

Tallahassee is hot and getting hotter. Examples
https://iopscience.iop.org/article/10.1088/2515-7620/ab27cf

Key finding: The most significant warming trends are occurring during summer although the largest trends are during fall and winter.

This increased warming averaged over decades amounts to about 

While a warmer than usual day in December can be experienced as quite pleasant a warmer than usual day in July will be experienced as dreadful.

According to the 6th IPCC assessment, the observed warming in 2010-2019 relative to 1850-1900 is 1.07C. For hot extremes, the evidence is mostly drawn from changes in metrics based on daily maximum temperatures.

But how individuals and communities experience climate change is local.
People of all political leanings highly value listening, openness, and talking about science in the context of society and community. Heslop, C., Dudo, A., and Copple, J. (August, 2021). Communicating science across political divides. Center for Media Engagement. https://mediaengagement.org/research/communicating-science- across-political-divides

Quote: "Throughout the WGI report and unless stated otherwise, uncertainty is quantified using 90% uncertainty intervals. The 90% uncertainty interval, reported in square brackets [x to y], is estimated to have a 90% likelihood of covering the value that is being estimated. The range encompasses the median value and there is an estimated 10% combined likelihood of the value being below the lower end of the range (x) and above its upper end (y). Often the distribution will be considered symmetric about the corresponding best estimate, but this is not always the case."

"For hot extremes, the evidence is mostly drawn from changes in metrics based on daily maximum temperatures" "There has been widespread evidence of human influence on various aspects of temperature extremes, at global, continental, and regional scales. This includes attribution to human influence of observed changes in intensity, frequency, and duration and other relevant characteristics at global and continental scales." But very little work at individual locations. What days of the year are warming the fastest in a particular location. Requires a long-term record

https://rmets.onlinelibrary.wiley.com/doi/10.1002/joc.4688 Diurnal asymmetry to the observed global warming

1. Writing is thinking, 
2. 1 article, 1 question, 1 finding, 
3. it's not a review of the literature, it's an argument about why we need your paper

We compute day-of-year trends in min and max temperature and use descriptive statistics to examine the range, distribution, and extremes of these trends.

What days of the year have warmed the most? What days of the year have cooled the most? 

show the biggest long-term (climate) trends in max/min temperatures? Corollary: what day(s) of the year do not show a normal distribution in max/min temperatures?

Why is this important? Climate change on the scale relevant to our daily lives.

The paper is exploratory. We do not have a hypothesis in mind though we do expect to see more days with upward than downward trends with perhaps a greater number of days with upward trends in low temperatures compared to the number of days with upward trends in high temperatures.

R version 4.1.0 (2021-05-18) -- "Camp Pontanezen"

## Get the required packages
```{r}
library(tidyverse)
library(lubridate)
library(broom)
library(patchwork)
library(ggdist)
```

## Get TLH airport daily weather data

NWS Tallahassee daily weather data:
https://www.ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USW00093805/detail
https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
Column explanations: https://docs.google.com/document/d/1q2WEpXndpMx9lUq-0ON63GojkaOzixrGyYEhI_xkPtw/edit?usp=sharing
AWND = Average daily wind speed (tenths of meters per second) 

Import the summary-of-the-day data and add columns to the data frame.
```{r}
TLH.df <- read_csv(file = 'Data/TLH_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE)) |>
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN) |>
  filter(Year <= 2019 & Year >= 1973)

Other.df <- read_csv(file = 'Data/LVG_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE)) |>
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN) |>
  filter(Year <= 2020 & Year >= 1943)


# Missing temperatures filled in from NCDC daily summaries https://gis.ncdc.noaa.gov/maps/ncei/summaries/daily

TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 95
TLH.df$MinTemp[TLH.df$Date == "2005-07-08"] <- 72
TLH.df$MinTemp[TLH.df$Date == "2002-10-08"] <- 67 # Quincy

# TLH.df <- Other.df  # Albany, GA shows more cooling days than warming days. Why? Milwaukee shows more warming days than cooling days like Tallahassee. Las Vegas shows more warming days and all warming nights.

TLH.df <- TLH.df |>
  mutate(AvgTemp = (MaxTemp + MinTemp) / 2,
         DTR = MaxTemp - MinTemp)

( AnnualMeanTemperatures <- TLH.df %>%
  group_by(Year) |>
  summarize(AvgHighTemp = mean(MaxTemp, na.rm = TRUE),
            AvgLowTemp = mean(MinTemp, na.rm = TRUE),
            AvgMeanTemp = mean(AvgTemp, na.rm = TRUE),
            AvgDTR = mean(DTR, na.rm = TRUE)) )
```

Diurnal temperature range (DTR) is an index that combines daily high and low temperatures. Since processes that result in high temperature are different than processes that result in low temperature, the metric is not useful for understanding the physics related to long-term trends in day-of-year temperature.

Plot annual mean temperatures.
```{r}
AnnualMeanTemperatures |>
  pivot_longer(cols = c(AvgHighTemp, AvgLowTemp)) |>
  ggplot(mapping = aes(x = Year, y = value, color = name)) +
  scale_y_continuous(limits = c(NA, NA)) +
  scale_color_discrete(guide = "none") +
    geom_point() +
    geom_smooth() +
  theme_minimal() +
  labs(title = "Annual mean daily high and low temperatures (°F)",
       subtitle = "Tallahasee, Florida (1943-2020)",
       x = "", y = "")

AnnualMeanTemperatures |>
  ggplot(mapping = aes(x = Year, y = AvgDTR)) +
  scale_y_continuous(limits = c(NA, NA)) +
    geom_point() +
    geom_smooth() +
  theme_minimal() +
  labs(title = "Annual mean DTR (°F)",
       subtitle = "Tallahasee, Florida (1943-2020)",
       x = "", y = "")
```

Histogram of daily high temperature.
```{r}
labels <- c(paste0(seq(30, 100, by = 10), "°", " F"))

p1 <- TLH.df |> 
ggplot(mapping = aes(x = MaxTemp)) + 
  geom_histogram(binwidth = 1, 
                 mapping = aes(fill = ..count..)) +
  scale_fill_viridis_c(guide = "none") +
  scale_x_continuous(breaks = seq(30, 100, by = 10),
                     labels = labels) +
  theme_minimal() +
  labs(title = "Frequency of daily high temperatures",
       subtitle = "Tallahassee, Florida (1943-2020)",
       x = "", y = "")

arrows <- 
  tibble(
    x1 = c(75),
    x2 = c(89),
    y1 = c(1025), 
    y2 = c(1200)
  )

p1 +
  annotate("text", x = 60, y = 1000, label = "The most common high temperature is 90° F") +
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(.08, "inch")), size = .5,
    color = "gray20", curvature = -.3)
```

Histogram of daily low temperature.
```{r}
labels <- c(paste0(seq(10, 80, by = 10), "°", " F"))

p2 <- TLH.df |> 
ggplot(mapping = aes(x = MinTemp)) + 
  geom_histogram(binwidth = 1, 
                 mapping = aes(fill = ..count..)) +
  scale_fill_viridis_c(guide = "none") +
  scale_x_continuous(breaks = seq(10, 80, by = 10),
                     labels = labels) +
  theme_minimal() +
  labs(title = "Frequency of daily low temperatures",
       subtitle = "Tallahassee, Florida (1943-2020)",
       x = "", y = "")

arrows <- 
  tibble(
    x1 = c(50),
    x2 = c(71),
    y1 = c(1250), 
    y2 = c(1460)
  )

p2 +
  annotate("text", x = 48, y = 1200, label = "The most common low temperature is 72° F") +
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(.08, "inch")), size = .5,
    color = "gray20", curvature = -.3)

labels <- c(paste0(seq(0, 50, by = 10), "°", " F"))
TLH.df |> 
ggplot(mapping = aes(x = DTR)) + 
  geom_histogram(binwidth = 1, 
                 mapping = aes(fill = ..count..)) +
  scale_fill_viridis_c(guide = "none") +
  scale_x_continuous(breaks = seq(0, 50, by = 10),
                     labels = labels) +
  theme_minimal() +
  labs(title = "Frequency of DTR",
       subtitle = "Tallahassee, Florida (1943-2020)",
       x = "", y = "")
```

Statistics of day-of-year max
```{r}
( StatsMax <- 
  TLH.df |>
  group_by(doy) |>
  summarize(Avg = mean(MaxTemp, na.rm = TRUE),
            Highest = max(MaxTemp, na.rm = TRUE),
            Lowest = min(MaxTemp, na.rm = TRUE),
            Q75 = quantile(MaxTemp, prob = .75, na.rm = TRUE),
            Q25 = quantile(MaxTemp, prob = .25, na.rm = TRUE)) )

StatsMax |>
  ggplot(mapping = aes(x = doy, y = Avg ) ) +
  scale_y_continuous(limits = c(NA, NA)) +
  geom_point() + 
  geom_errorbar(mapping = aes(ymin = Lowest, ymax = Highest)) +
  geom_errorbar(mapping = aes(ymin = Q25, ymax = Q75), width = .5) +
  scale_x_continuous(position = "bottom", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal() 
```
            

Non-normality in daily high and low temperatures
```{r}
library(moments)
library(nortest)

TLH.df |>
  group_by(doy) |>
  summarize( AvgByDoY = mean(MaxTemp),
             SDByDoY = sd(MaxTemp),
             SkewByDoY = skewness(MaxTemp) ) |>
  ggplot(mapping = aes(x = doy, y = AvgByDoY)) +
    geom_point() +
    geom_smooth()

normTest <- 
  TLH.df |>
#  filter(Year >= 1990) |>
  group_by(doy) |>
  summarize(pvalueS = shapiro.test(MaxTemp)$p.value,
            pvalueAD = ad.test(MaxTemp)$p.value) |>
  mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
         Month = month(Date))

normTest$doy[normTest$pvalueS > .05]
normTest$doy[normTest$pvalueAD > .05]

normTest |>
  group_by(Month) |>
  summarize(nNormalDaysS = sum(pvalueS > .05),
            nNormalDaysAD = sum(pvalueAD > .05))
  
```

Compute long-term trends in day-of-year high temperatures.
```{r}
( TrendsMax <- 
  TLH.df |> 
    mutate(MaxTemp = (MaxTemp - 32) * 5/9) |>
#  filter(Year >= 1990) |>
    group_by(doy) |>
    do(tidy(lm(MaxTemp ~ Year, data = .))) |>
    filter(term == "Year") |>
    mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
           estimate = estimate * 10,
           std.error = std.error * 10,
           ymax = estimate + std.error,
           ymin = estimate - std.error,
           Daily = "High") )

sum(TrendsMax$estimate > 0)
sum(TrendsMax$estimate < 0)
sum(TrendsMax$estimate > 0) / length(TrendsMax$estimate) * 100

TrendsMax |> 
  arrange(desc(statistic))
```

Histogram of trends
```{r}
library(scales)

TrendsMax |>
  mutate(bins = cut(estimate, seq(-1, 1, by = .1))) |>
  group_by(bins) |>
  ggplot(mapping = aes(x = estimate)) +
    geom_histogram(binwidth = .1, 
                   color = "white",
                   fill = "gray70",
                   center = .05) +
    geom_rug(mapping = aes(color = estimate)) +
    scale_color_gradient2(low = muted("blue"),
                          mid = "white",
                          high = muted("red"),
                          midpoint = 0,
                          guide = 'none') +
    theme_minimal()
```

Plot largest trends.
```{r}
labels <- c(paste0(seq(40, 100, by = 20), "°", " F"))

( p3 <- TLH.df |>
  filter(doy %in% c(26, 350, 192)) |>
  mutate(Day = factor(case_when(Day = doy == 26 ~ "26 January\n Most cooling",
                         Day = doy == 350 ~ "15/16 December\n Most warming",
                         Day = doy == 192 ~ "10/11 July\n Most significant warming"), 
                      levels = c("26 January\n Most cooling", 
                                 "10/11 July\n Most significant warming", "15/16 December\n Most warming"))) |>
  select(Year, MaxTemp, Date, doy, Day) |>
  ggplot(mapping = aes(x = Year, y = MaxTemp)) +
    geom_smooth(method = lm, mapping = aes(color = Day)) +
    geom_point(color = "gray70") +
  facet_wrap(~ Day) +
  scale_x_continuous(breaks = seq(50, 2010, by = 20)) +
  scale_y_continuous(limits = c(35, 105), breaks = seq(40, 100, by = 20), labels = labels) +
  scale_color_manual(values = c("skyblue", "red3", "salmon"), guide = "none") +
  theme_minimal() +
  labs(x = "", y = "") )
```

Plot trends.
```{r}
( p4 <- TrendsMax |>
  ggplot(mapping = aes(x = doy, y = estimate, color = statistic) ) +
  scale_y_continuous(limits = c(-1.5, 1.5)) +
  geom_hline(yintercept = 0, color = "gray70") +
  scale_color_distiller(palette = "RdBu", 
                        limits = c(-1, 1) * max(abs(TrendsMax$statistic)),
                        guide = 'none') +
  geom_point() + 
  geom_errorbar(mapping = aes(ymin = ymin, ymax = ymax)) +
  scale_x_continuous(position = "bottom", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal() +
  labs(title = "Long-term (1943-2020) trends in day-of-year high temperatures (°C/decade), Tallahassee, Florida", 
       subtitle = "Darker colors indicate greater evidence against a no-change hypothesis", 
       x = "", y = "") )
```

Annotate the plot with text and arrows. http://jenrichmond.rbind.io/post/idhtg-how-to-annotate-plots/
```{r}
arrows <- 
  tibble(
    x1 = c(46, 130, 300),
    x2 = c(26.5, 191.5, 349.5),
    y1 = c(-.95, .75, 1.15), 
    y2 = c(-.525, .48, .8)
  )

( p4 <- p4 + 
  annotate("text", x = 46, y = -1, label = "Most cooling: 26 January") +
  annotate("text", x = 130, y = .8, label = "Most significant warming: 10/11 July") +
  annotate("text", x = 300, y = 1.2, label = "Most warming: 15/16 December") +
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(.08, "inch")), size = .5,
    color = "gray50", curvature = .3) )

library(patchwork)

p4/p3
```

Compute long-term trends in day-of-year low temperatures.
```{r}
( TrendsMin <- 
  TLH.df |> 
    mutate(MinTemp = (MinTemp - 32) * 5/9) |>
#  filter(Year >= 1990) |>
    group_by(doy) |>
    do(tidy(lm(MinTemp ~ Year, data = .))) |>
    filter(term == "Year") |>
    mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
           estimate = estimate * 10,
           std.error = std.error * 10,
           ymax = estimate + std.error,
           ymin = estimate - std.error,
           Daily = "Low") )

sum(TrendsMin$estimate > 0)
sum(TrendsMin$estimate < 0)
sum(TrendsMin$estimate > 0) / length(TrendsMin$estimate) * 100

TrendsMin |> 
  arrange(desc(statistic))
```

```{r}
Trends <- rbind(TrendsMax, TrendsMin)
Trends |>
  ggplot(mapping = aes(y = estimate, x = Daily, color = estimate)) +
    stat_slab(color = "gray80", fill = "gray80", scale = .5) +
    geom_point(shape = "_", size = 4) + 
#    stat_halfeye() +
    scale_color_gradient2(low = "blue",
                          mid = "white",
                          high = "red",
                          midpoint = 0,
                          guide = 'none') +
    scale_y_continuous(limits = c(-1.5, 1.5)) +
    theme_minimal() +
    labs(x = "", y = "", title = "Day-of-year trends in high and low temperatures (°C/decade)")
```

```{r}
labels <- c(paste0(seq(20, 80, by = 20), "°", " F"))

( p5 <- TLH.df |>
  filter(doy %in% c(26, 335, 193)) |>
  mutate(Day = factor(case_when(Day = doy == 26 ~ "26 January\n Most cooling",
                         Day = doy == 335 ~ "30 November/1 December\n Most warming",
                         Day = doy == 193 ~ "11/12 July\n Most significant warming"), 
                      levels = c("26 January\n Most cooling", 
                                 "11/12 July\n Most significant warming", 
                                 "30 November/1 December\n Most warming"))) |>
  select(Year, MinTemp, Date, doy, Day) |>
  ggplot(mapping = aes(x = Year, y = MinTemp)) +
    geom_smooth(method = lm, mapping = aes(color = Day)) +
    geom_point(color = "gray70") +
  facet_wrap(~ Day) +
  scale_x_continuous(breaks = seq(1950, 2010, by = 20)) +
  scale_y_continuous(limits = c(15, 85), breaks = seq(20, 80, by = 20), labels = labels) +
  scale_color_manual(values = c("skyblue", "red3", "salmon"), guide = "none") +
  theme_minimal() +
  labs(x = "", y = "") )

( p6 <- Trends |>
  ggplot(mapping = aes(x = doy, y = estimate, color = statistic) ) +
  scale_y_continuous(limits = c(-1.2, 1.2)) +
  geom_hline(yintercept = 0, color = "gray70") +
  scale_color_distiller(palette = "RdBu", 
                        limits = c(-1, 1) * max(abs(Trends$statistic)),
                        guide = 'none') +
  geom_point() + 
  geom_errorbar(mapping = aes(ymin = ymin, ymax = ymax)) +
  scale_x_continuous(position = "bottom", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal() +
  labs(title = "Long-term (1943-2020) trends in day-of-year low temperatures (°C/decade), Tallahassee, Florida", 
       subtitle = "Darker colors indicate greater evidence against a no-change hypothesis", 
       x = "", y = "") )

arrows <- 
  tibble(
    x1 = c(56, 160, 300),
    x2 = c(26.5, 192.5, 335),
    y1 = c(-1.1, .75, 1.15), 
    y2 = c(-.815, .287, .671)
  )

( p6 <- p6 + 
  annotate("text", x = 66, y = -1.15, label = "Most cooling: 26 January") +
  annotate("text", x = 130, y = .8, label = "Most significant warming: 11/12 July") +
  annotate("text", x = 300, y = 1.2, label = "Most warming: 30 November/1 December") +
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(.08, "inch")), size = .5,
    color = "gray50", curvature = .3) )

p6/p5
```

Personalized weather history. Use daily data from locations where you lived. Start with my weather history. Daily data starting with October 16, 1959 from Milwaukee then from August 1, 1990 from Tallahassee.

Cumulative daily anomalies. 1. Compute the average doy temperature  for each lived location, 2. Subtract average doy from the calendar day temperature conditional on lived location to create a lived anomaly. 3. Start with the birth date and sum the lived anomaly for each day of life.

```{r}
MKE_BirthClimate.df <- 
  read_csv(file = 'Data/MKE_Daily1938-2020.csv') |>
  mutate(Date = as.Date(DATE),
         DoY = yday(Date)) |>
  filter(Date < as.Date("1959-10-16")) |>
  group_by(DoY) |>
  summarize(AvgDoY = mean(TMAX, na.rm = TRUE)) 

TLH_BirthClimate.df <- 
  read_csv(file = 'Data/TLH_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE),
         DoY = yday(Date)) |>
  filter(Date < as.Date("1959-10-16")) |>
  group_by(DoY) |>
  summarize(AvgDoY = mean(TMAX, na.rm = TRUE)) 

MKE.df <-
  read_csv(file = 'Data/MKE_Daily1938-2020.csv') |>
  mutate(Date = as.Date(DATE),
         DoY = yday(Date)) |>
  filter(Date >= as.Date("1959-10-16") & Date < as.Date("1990-08-01")) |>
  select(STATION, Date, DoY, TMAX) |>
  left_join(MKE_BirthClimate.df)

TLH.df <-
  read_csv(file = 'Data/TLH_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE),
         DoY = yday(Date)) |>
  filter(Date >= as.Date("1990-08-01") & Date < as.Date("2021-01-01")) |>
#  filter(Date >= as.Date("1959-10-16") & Date < as.Date("2021-01-01")) |>
  select(STATION, Date, DoY, TMAX) |>
  left_join(TLH_BirthClimate.df)

LocationHistory.df <- rbind(MKE.df, TLH.df) |>
  mutate(Anomaly = TMAX - AvgDoY,
         Anomaly = replace_na(Anomaly, 0),
         AnomalyF = as.factor(Anomaly >= 0))
LocationHistory.df$CumAnom <- cumsum(LocationHistory.df$Anomaly)

ggplot(data = LocationHistory.df,
       mapping = aes(x = Date, y = CumAnom / nrow(LocationHistory.df), color = AnomalyF)) +
#  geom_line(size = 2) +
  geom_point(shape = ".", size = 1) +
  scale_color_manual(values = c("blue", "red"),
                     guide = 'none') +
  theme_minimal() +
  labs(title = "Personal lifetime cumulative temperature anomalies (F/day)",
       subtitle = "Day-of-year averages based on data prior to my birth at locations where I lived",
       x = "", y = "")
```

```{r}
( Trends <- 
    LifeHistory.df |> 
    mutate(MaxTemp = (MaxTemp - 32) * 5/9) |>
    group_by(doy, NAME) |>
    do(tidy(lm(MaxTemp ~ Year, data = .))) |>
    filter(term == "Year") |>
    mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
           estimate = estimate * 10,
           std.error = std.error * 10,
           ymax = estimate + std.error,
           ymin = estimate - std.error) )
sampleSize <- 
    LifeHistory.df |>
    group_by(doy, NAME) |>
    summarize(sampleSize = n()) |>
    pull(sampleSize)

Trends$sampleSize <- sampleSize

sum(Trends$estimate > 0)
sum(Trends$estimate < 0)
sum(Trends$estimate > 0) / length(Trends$estimate) * 100

Trends |>
  ggplot(mapping = aes(y = estimate, x = NAME, color = estimate)) +
    stat_slab(color = "gray80", fill = "gray80", scale = .5) +
    geom_point(shape = "_", size = 4) + 
    scale_color_gradient2(low = "blue",
                          mid = "white",
                          high = "red",
                          midpoint = 0,
                          guide = 'none') +
    scale_y_continuous(limits = c(NA, NA)) +
    theme_minimal() +
    labs(x = "", y = "", title = "Day-of-year trends in high temperatures (°C/decade)")


Trends |>
  group_by(doy) |>
  summarize(AvgTrend = weighted.mean(estimate, sampleSize)) |>
  ggplot(mapping = aes(x = doy, y = AvgTrend, color = AvgTrend)) +
   scale_y_continuous(limits = c(-2.5, 2.5)) +
   geom_hline(yintercept = 0, color = "gray70") +
   scale_color_gradient2(low = "blue",
                         high = "red",
                         midpoint = 0,
                         guide = 'none') +
   geom_point() +
   scale_x_continuous(position = "bottom", 
                      breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal()
  

ggplot(mapping = aes(x = doy, y = estimate, color = statistic) ) +
  scale_y_continuous(limits = c(-1.5, 1.5)) +
  geom_hline(yintercept = 0, color = "gray70") +
  scale_color_distiller(palette = "RdBu", 
                        limits = c(-1, 1) * max(abs(TrendsMax$statistic)),
                        guide = 'none') +
  geom_point() + 
  geom_errorbar(mapping = aes(ymin = ymin, ymax = ymax)) +
  scale_x_continuous(position = "bottom", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal() +
  labs(title = "Long-term (1943-2020) trends in day-of-year high temperatures (°C/decade), Tallahassee, Florida", 
       subtitle = "Darker colors indicate greater evidence against a no-change hypothesis", 
       x = "", y = "") )

```



Use `ggdist` package from Matthew Kay.
```{r}
TLH.df |>
  ggplot(mapping = aes(x = doy, y = MaxTemp)) +
  stat_pointinterval(show_point = FALSE) +
  scale_color_brewer() +
  scale_x_continuous(position = "top", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal()

Trends |>
  ggplot(mapping = aes(y = estimate)) 
```

Note if you use the `GSODR` package. MAX - Maximum temperature reported during the day converted to Celsius to tenths--time of max temp report varies by country and region, so this will sometimes not be the max for the calendar day. Missing = NA;
```{r}
# install.packages("GSODR")
library(GSODR)

near_stations <- nearest_stations(LAT =  30.3965, #42.9476,
                                  LON = -84.3503,
                                  distance = 3)

#data <- get_GSOD(years = 2000:2019, station = "726400-14839")
data <- get_GSOD(years = 1943:2020, station = "722140-93805")
data2 <- get_GSOD(years = 1943:2020, station = "999999-93805")
str(data)

data |> 
  group_by(YDAY) |>
  summarize(AvgMax = mean(MAX, na.rm = TRUE)) |>
ggplot(mapping = aes(x = YDAY, y = AvgMax)) +
  geom_line()

TLH2.df <- data |>
  select(Date = YEARMODA,
         Year = YEAR,
         Month = MONTH,
         doy = YDAY,
         MAX,
         MIN) |>
  mutate(MaxTemp = MAX * 1.8 + 32,
         MinTemp = MIN * 1.8 + 32) |>
  filter(Year <= 2020 & Year >= 1943)

AnnualMeanTemperatures <- TLH2.df %>%
  group_by(Year) |>
  summarize(AvgHighTemp = mean(MaxTemp, na.rm = TRUE),
            AvgLowTemp = mean(MinTemp, na.rm = TRUE)) 
AnnualMeanTemperatures |>
  pivot_longer(cols = c(AvgHighTemp, AvgLowTemp)) |>
  ggplot(mapping = aes(x = Year, y = value, color = name)) +
  scale_y_continuous(limits = c(NA, NA)) +
  scale_color_discrete(guide = "none") +
    geom_point() +
    geom_smooth() +
  theme_minimal() +
  labs(title = "Annual mean daily high and low temperatures (°F)",
       subtitle = "Tallahasee, Florida (1943-2020)",
       x = "", y = "")


TrendsMax <- 
  TLH2.df |> 
    mutate(MaxTemp = (MaxTemp - 32) * 5/9) |>
    group_by(doy) |>
    do(tidy(lm(MaxTemp ~ Year, data = .))) |>
    filter(term == "Year") |>
    mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
           estimate = estimate * 10,
           std.error = std.error * 10,
           ymax = estimate + std.error,
           ymin = estimate - std.error,
           Daily = "High") 

sum(TrendsMax$estimate > 0)
sum(TrendsMax$estimate < 0)
sum(TrendsMax$estimate > 0) / length(TrendsMax$estimate) * 100

TrendsMax |>
  ggplot(mapping = aes(x = doy, y = estimate, color = statistic) ) +
  scale_y_continuous(limits = c(-1.5, 1.5)) +
  geom_hline(yintercept = 0, color = "gray70") +
  scale_color_distiller(palette = "RdBu", 
                        limits = c(-1, 1) * max(abs(TrendsMax$statistic)),
                        guide = 'none') +
  geom_point() + 
  geom_errorbar(mapping = aes(ymin = ymin, ymax = ymax)) +
  scale_x_continuous(position = "bottom", 
                     breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336) + 13,
                     labels = month.name) +
  theme_minimal() +
  labs(title = "Long-term (1943-2020) trends in day-of-year high temperatures (°C/decade), Tallahassee, Florida", 
       subtitle = "Darker colors indicate greater evidence against a no-change hypothesis", 
       x = "", y = "")
```