# '
# Uses the panel data stored in "bld/data/" to
# create Supplementary Figure SI5 of the paper and saves the plot as a pdf file to "bld/figures/".
# 
# '

rm(list = ls())
options(scipen = 999)
options(warn = -1) # 0 to turn on again
source(paste0("src/code/ggplot_theme_publication.r")) # Load a user-written ggplot theme and style.

# Load libraries
library(tidyverse) # More intuitive data wrangling
library(scales)
library(Cairo)
library(rjson)

# Load model specification
model <- fromJSON(file = paste0("src/code/specifications.json"))

# Load the data.
analysis <- readRDS(paste0("bld/data/analysis_month.RDS"))

# Drop duplicates and keep only obs at 100km bandwidth
dat <- analysis %>%
  distinct(mcode, yearmonth, .keep_all = TRUE) %>%
  arrange(mcode, yearmonth)

d <- subset(dat, dat$dist2cutoff <= model$distance_to_border & dat$year > 2002)  

p <- d %>%
  arrange(mcode, yearmonth) %>%
  group_by(yearmonth) %>% #, treated) %>%
  summarise(Fires = mean(focos_queimada, na.rm = TRUE),
            `PM2.5`= mean(pm25_ugm3, na.rm = TRUE),
            cases_prevalence = mean(cases_prevalence, na.rm = TRUE),
            death_prevalence = mean(death_prevalence, na.rm = TRUE),
            temperature = mean(temperatura_c, na.rm = TRUE),
            precipitacion = mean(precipitacao_mmdia, na.rm = TRUE),
            pm25_over_limit_relative = mean(pm25_over_limit_relative, na.rm = TRUE),
            humidity = mean(umidade_relativa_percentual, na.rm = TRUE)) %>%
  ungroup() %>% as.data.frame() 

# Convert yearmonth to Date
p$yearmonth <- as.Date(paste0(p$yearmonth, "-01"), format = "%Y-%m-%d")

# Extract month information
p$month <- format(p$yearmonth, "%m")

# Reshape the data for plotting with legend
p_long <- p %>%
  pivot_longer(cols = c(Fires, `PM2.5`), names_to = "variable", values_to = "value")

linetypes <- c("Fires" = "solid", "PM2.5" = "dashed")
  
# Plot the data
pl <- ggplot(data = p_long, aes(x = yearmonth, y = value, color = variable, shape = variable, linetype = variable
                                , group = variable)) +
  geom_line() +
  geom_point(size=.9) +
  
  scale_y_continuous(
    name = "Number of fires",
    limits = c(0,130), 
    breaks = seq(0, 130, 10),
    labels = comma,
    sec.axis = sec_axis(~., name = "PM2.5 concentration (μg/m³)",
                        labels = function(b) { format(b, big.mark = ",", decimal.mark = ".", scientific = FALSE) })
  ) +
  
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  xlab("Year") +
  scale_linetype_manual(values = linetypes) +
  scale_colour_Publication() +
  scale_fill_Publication() +
  theme_Publication() +
  theme(
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(pl) 

# Save it 
ggsave(plot = pl, filename = "bld/figures/monthly.pdf", device = cairo_pdf, height = 9, width = 15, units = "cm")