# '
# Creates the fire regression difference in discontinuities graph and the aggregated fires graph.
# Takes the data binned fire dataset stored in "bld/data/" to create
# the Supplementary Figures SI6 and SI7, which is stored as pdf to "bld/figures/".
# 
# '

rm(list = ls())
# load libraries
library(tidyverse)
library(scales)
# disable scientific notation
options(scipen=999)
# load user-written plot themes
source("src/code/ggplot_theme_publication.r")

# load data
d <- readRDS("bld/data/fire_rdid.rds")

# Create binned scatterplot 
binned_scatterplot <- d %>%
  filter(distance >= -250000 & distance <= 250000) %>%
  ggplot(., aes(distance/1000, log(fires), group = time.f, colour = time.f, shape = time.f)) + 
  geom_point(size = .8, alpha = .7) + 
  scale_x_continuous(
    name = "Distance to border (km)",
    labels = scales::comma,
    limits = c(-251,251),
    breaks = seq(-250, 250, 50)) +
  scale_y_continuous(
    name = "log(Number of fires)")+#,
  geom_vline(xintercept = 0) +
  scale_fill_Publication() +
  scale_colour_Publication() +
  theme_Publication() +
  theme(
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1))
  
# view graph
print(binned_scatterplot)

# Add before/after inside/outside LOESS regression lines to graph
p <- binned_scatterplot + 
  geom_smooth(data = subset(d, d$time == 0 & d$treated == 0 & d$distance > -250000 & d$distance < 250000),
              method = loess, formula = "y ~ 1 + x", se = F, colour = "blue", linewidth = 1) +
  geom_smooth(data = subset(d, d$time == 0 & d$treated == 1 & d$distance > -250000 & d$distance < 250000),
              method = loess, formula = "y ~ 1 + x", se = F, colour = "blue", linewidth = 1) +
  geom_smooth(data = subset(d, d$time == 1 & d$treated == 0 & d$distance > -250000 & d$distance < 250000),
              method = loess, formula = "y ~ 1 + x", se = F, colour = "red", linewidth = 1, linetype="dashed") +
  geom_smooth(data = subset(d, d$time == 1 & d$treated == 1 & d$distance > -250000 & d$distance < 250000),
              method = loess, formula = "y ~ 1 + x", se = F, colour = "red", linewidth = 1, linetype="dashed") 

# view graph
print(p)

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

# Create plot of aggregated number of fires within 100km of biome border
p <- d %>%
  filter(distance >= -100000 & distance <= 100000) %>%
  group_by(treated.f, year) %>%
  summarise(fires = sum(fires, na.rm=TRUE)) %>%
  ggplot(., aes(year, fires, group = treated.f, colour = treated.f, shape = treated.f, linetype = treated.f)) +
  geom_line() +
  geom_point() +
  geom_vline(xintercept = 2006) +
  scale_x_continuous(
    name = "Year",
    breaks = 2003:2017) +
  scale_y_continuous(
    name = "Number of fires",
    labels = scales::comma,
    limits = c(0, 65000),
    breaks = seq(0, 65000, 5000)) +
  scale_colour_Publication() +
  scale_fill_Publication() +
  scale_linetype_manual(values = c("solid", "dashed")) +
  theme_Publication() +
  theme(
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1))

# view plot
print(p)

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