# This script is used to create polar plots and pollution roses for the EM27 data
# Load required libraries
library(ggplot2)
library(grid)
library(openair)
library(latticeExtra)
library(dplyr)
library(readr)   # safer than base read.csv

all_rr_df_path <- '/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/Projects/Ratio_paper_2025_v5/Data/ha/all_rolling_regr_df_filtered.csv'

# Use readr::read_csv instead of base read.csv
all_rr_df <- read_csv(all_rr_df_path, show_col_types = FALSE)

# Clean up column names (remove spaces, special chars)
names(all_rr_df) <- make.names(names(all_rr_df), unique = TRUE)

# Inspect what columns exist that match
print(grep("good_group", names(all_rr_df), value = TRUE))

# Parse datetime properly
all_rr_df <- all_rr_df %>%
  mutate(
    dt   = as.POSIXct(dt, format = "%Y-%m-%d %H:%M:%S", tz = "America/Denver"),
    date = as.Date(dt)
  )

[1] "ch4_co2_good_group" "co_co2_good_group" 

# Set the regression label and type
regr_label <- 'co_co2'
x_name = "xco2.ppm._wmean_ex5q_divak"
if (regr_label == 'ch4_co2'){
  y_name = "xch4.ppm._wmean_ex5q_divak"
} else if (regr_label == 'co_co2'){
  y_name = "xco.ppm._wmean_ex5q_divak"
}
regr_type <- 'york'
permil = TRUE

save_plots <- TRUE
show_plots <- FALSE

# Figure out the right "good group" column
good_col <- grep(paste0("^", regr_label, "_good_group$"), names(all_rr_df), value = TRUE)

# Filter safely
frr_df <- all_rr_df %>%
  filter(!is.na(.data[[good_col]]))
if (permil) {
  frr_df <- frr_df %>%
    mutate(across(.cols = grep("slope", names(frr_df), value = TRUE), 
                  .fns = ~ . * 1000))
}
stat_type <- 'slope' # the statistic to plot
var_of_interest <- paste(regr_label,regr_type,stat_type,sep='_') # the variable of interest from the regression label and type

if (save_plots) {
png("/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/Projects/Ratio_paper_2025_v5/Figures/ha/polar_co_co2.png",width=2000,height=1600,res=300)
polarPlot(frr_df,
          pollutant = var_of_interest,
          statistic = 'mean',
          cols = 'viridis',
          limits = c(0, 11),
          angle.scale = 45,
          par.settings = list(fontsize = list(text = 30)),
          units = '',
          key.footer = NULL,
)
dev.off()
}

if (show_plots) {
  polarPlot(frr_df,
            pollutant = var_of_interest,
            statistic = 'mean',
            cols = 'viridis',
            limits = c(0, 11),
            angle.scale = 45,
            par.settings = list(fontsize = list(text = 30)),
            units = '',
            key.footer = NULL,
  )
}

pdf: 2

# Set the regression label and type
regr_label <- 'ch4_co2'
x_name = "xco2.ppm._wmean_ex5q_divak"
if (regr_label == 'ch4_co2'){
  y_name = "xch4.ppm._wmean_ex5q_divak"
} else if (regr_label == 'co_co2'){
  y_name = "xco.ppm._wmean_ex5q_divak"
}
regr_type <- 'york'
permil = TRUE

save_plots <- TRUE
show_plots <- FALSE

# Figure out the right "good group" column
good_col <- grep(paste0("^", regr_label, "_good_group$"), names(all_rr_df), value = TRUE)

# Filter safely
frr_df <- all_rr_df %>%
  filter(!is.na(.data[[good_col]]))
if (permil) {
  frr_df <- frr_df %>%
    mutate(across(.cols = grep("slope", names(frr_df), value = TRUE), 
                  .fns = ~ . * 1000))
}
stat_type <- 'slope' # the statistic to plot
var_of_interest <- paste(regr_label,regr_type,stat_type,sep='_') # the variable of interest from the regression label and type

if (save_plots) {
png("/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/Projects/Ratio_paper_2025_v5/Figures/ha/polar_ch4_co2.png",width=2000,height=1600,res=300)
polarPlot(frr_df,
          pollutant = var_of_interest,
          statistic = 'mean',
          cols = 'viridis',
          limits = c(0, 8.8),
          angle.scale = 45,
          par.settings = list(fontsize = list(text = 30)),
          units = '',
          key.footer = NULL,

)
dev.off()
}

if (show_plots) {
  polarPlot(frr_df,
            pollutant = var_of_interest,
            statistic = 'mean',
            cols = 'viridis',
            limits = c(0, 8.8),
            angle.scale = 45,
            par.settings = list(fontsize = list(text = 30)),
            units = '',
            key.footer = NULL,
  )
}

pdf: 2