This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(tidyverse)
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

lake_ice_noNA <- read_csv("douglaslake_southfishtailbay_icecover_1931-2024.csv",) %>%
  drop_na()
## Rows: 105 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): observation
## date (1): date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
lake_ice_noNA
write_csv(lake_ice_noNA,"douglaslake_southfishtailbay_icecover_1931-2024.csv")
# Convert the 'date' column to a Date object
lake_ice <- lake_ice_noNA %>%
  mutate(date = ymd(date))

# Create a lagged version of the 'date' and 'observation' column and calculate the duration
lake_ice <- lake_ice %>%
  arrange(date) %>%
  mutate(lead_observation = lead(observation),
         lead_date = lead(date)) %>%
  filter(observation == "ice in", lead_observation == "ice out") %>%
  mutate(duration = lead_date - date)

# Identify winter by subtracting a year from the 'date' if 'date' is in Jan or Feb
lake_ice <- lake_ice %>%
  mutate(winter = ifelse(month(date) %in% 1:2, year(date) - 1, year(date)))

# Create a column for the winter season string representation
lake_ice <- lake_ice %>%
  mutate(winter = paste(winter, winter + 1, sep="-"))

# Select the relevant columns
new_lake_ice <- lake_ice %>%
  select(winter, date, lead_date, duration)

# Rename the columns to be more informative
new_lake_ice <- rename(new_lake_ice,
                       `Ice On Date` = date,
                       `Ice Off Date` = lead_date,
                       `Duration` = duration)

# View the new table
print(new_lake_ice)
## # A tibble: 50 × 4
##    winter    `Ice On Date` `Ice Off Date` Duration
##    <chr>     <date>        <date>         <drtn>  
##  1 1933-1934 1933-11-16    1934-04-28     163 days
##  2 1950-1951 1950-11-24    1951-04-13     140 days
##  3 1973-1974 1973-12-12    1974-04-21     130 days
##  4 1974-1975 1974-12-13    1975-05-01     139 days
##  5 1975-1976 1975-12-07    1976-04-14     129 days
##  6 1979-1980 1979-12-02    1980-04-21     141 days
##  7 1980-1981 1980-12-03    1981-04-01     119 days
##  8 1981-1982 1981-12-11    1982-04-30     140 days
##  9 1982-1983 1982-12-12    1983-04-14     123 days
## 10 1983-1984 1983-12-03    1984-04-13     132 days
## # ℹ 40 more rows
library(ggplot2)  # Part of tidyverse, for plotting

# Convert Duration to numeric in case it's not already.
new_lake_ice <- new_lake_ice %>%
  mutate(Duration = as.numeric(Duration))

# Create a column with just the starting year of the winter for plotting
new_lake_ice <- new_lake_ice %>%
  mutate(Year = as.numeric(sub("-.*", "", winter)))

# Use ggplot to create a scatter plot with Year and Duration
ggplot(new_lake_ice, aes(x = Year, y = Duration)) +
  geom_point() +  # Add points representing each observation
  geom_smooth(method = "lm", se = FALSE, color = "blue") +  # Add a linear trend line, no confidence interval shading
  labs(x = "Year", y = "Ice Duration (days)", title = "Trend of Lake Ice Duration Over Years") +
  theme_minimal()  # Use the minimal theme for a cleaner look
## `geom_smooth()` using formula = 'y ~ x'

library(ggplot2)
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.3.3
# Modify the data to fix November as the start of the period
new_lake_ice <- new_lake_ice %>%
  mutate(`Ice On Fixed` = as.Date(ifelse(month(`Ice On Date`) < 6, 
                                         format(`Ice On Date`, "2001-%m-%d"), 
                                         format(`Ice On Date`, "2000-%m-%d"))),
         `Ice Off Fixed` = as.Date(ifelse(month(`Ice Off Date`) < 6, 
                                          format(`Ice Off Date`, "2001-%m-%d"), 
                                          format(`Ice Off Date`, "2000-%m-%d"))))

# Custom colors from the user
ice_on_color <- "#1b9e77"
ice_off_color <- "#d95f02"
duration_color <- "#7570b3"

# Create the 'Ice On' and 'Ice Off' dates plot
ice_dates_plot <- ggplot(new_lake_ice, aes(x = Year)) +
  geom_point(aes(y = `Ice On Fixed`), color = ice_on_color) +
  geom_point(aes(y = `Ice Off Fixed`), color = ice_off_color) +
  geom_smooth(aes(y = `Ice On Fixed`), method = "lm", se = FALSE, color = ice_on_color) +
  geom_smooth(aes(y = `Ice Off Fixed`), method = "lm", se = FALSE, color = ice_off_color) +
  scale_y_date(date_breaks = "1 month", labels = scales::date_format("%B")) +
  labs(title = "Douglas Lake ice on, off and duration") +
  theme_minimal() +
  theme(plot.margin = margin(b = 0),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none")

# Annotation position calculation
last_year <- max(new_lake_ice$Year)
latest_ice_on <- max(new_lake_ice$`Ice On Fixed`, na.rm = TRUE)
earliest_ice_off <- min(new_lake_ice$`Ice Off Fixed`, na.rm = TRUE)

# Annotate the Ice On and Ice Off trend lines
ice_dates_plot <- ice_dates_plot +
  annotate("text", x = last_year, y = latest_ice_on, label = "Ice On", hjust = 13.3, vjust = 2, 
           color = ice_on_color, size = 3.5, angle = 0) +
  annotate("text", x = last_year, y = earliest_ice_off, label = "Ice Off", hjust = 13.6, vjust = -0.5, 
           color = ice_off_color, size = 3.5, angle = 0)

# Modify the 'Duration' plot to include linear trend lines and annotations
shortest_duration <- min(new_lake_ice$Duration, na.rm = TRUE)
duration_plot <- ggplot(new_lake_ice, aes(x = Year, y = Duration)) +
  geom_point(color = duration_color) +
  geom_smooth(method = "lm", se = FALSE, color = duration_color) +
  annotate("text", x = last_year, y = shortest_duration, label = "Ice Duration", hjust = 7.3, vjust = -9, 
           color = duration_color, size = 3.5, angle = 0) +
  labs(x = "Year", y = "Ice Duration (days)") +
  theme_minimal()

# Combine the plots
combined_plot <- (ice_dates_plot + plot_layout(heights = c(1))) / 
                 (duration_plot + plot_layout(heights = c(1)))

# Display the combined plot
combined_plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

ggsave(combined_plot, 
       filename = "UMBS_ice_cover_plot.jpeg",
       device = "jpeg",
       height = 4, width = 6, units = "in")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'