# '
# Uses panel data from "bld/data/" and GIS data from "src/original_data/ to
# create the maps Figure 1 and Supplementary Figure SI8 in the paper.
# Stores both maps to "bld/figures/".
# 
# '

rm(list = ls())
options(scipen = 999) # Turn off scientific notation
options(warn = 0) # 0 to turn on again, -1 to turn off

# Load libraries
library(tidyverse) # More intuitive data wrangling
library(rjson) # Read JSON files in R
library(sf) # GIS package
library(tmap)
library(tmaptools)
library(rnaturalearth)

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

# Load panel data
analysis <- readRDS(paste0("bld/data/analysis_year.RDS"))

# Keep only unique observations
d <- analysis %>%
  distinct(mcode, .keep_all = TRUE)
rm(analysis)

# Create point geometries from capital location
d <- st_as_sf(d, coords = c("mun_lon", "mun_lat"), crs = 4326)

# Load Amazon biome border line
biome_line <- st_read(paste0("src/original_data/biome_line.gpkg"))
biome_line <- st_transform(biome_line, crs = 29101)

# Create buffer around the biome border line.
line_buffer_100km <- st_buffer(biome_line, model$distance_to_border)

# Load Legal Amazon polygon
bla <- st_read(paste0("src/original_data/brazilian_legal_amazon.shp"))

# Load municipality polygons
municipalities <- st_read(paste0("src/original_data/municipalities_2013.shp"))
municipalities <- st_transform(municipalities, crs = "WGS84")
municipalities <- municipalities %>%
  select(code_mn,name_mn) %>%
  rename(mcode = code_mn,
         munResNome = name_mn)%>%
  mutate(mcode = floor(mcode/10))

# Retransform to common CRS
biome_line <- st_transform(biome_line, crs = "WGS84")
bla <- st_transform(bla, crs = "WGS84")

# Create subset geometry of municipalities within 100km of biome border
mun_100km <- subset(d, d$dist2cutoff < 100000)
mun_100km$geometry <- NULL;
mun_100km <- merge(mun_100km, municipalities, by = "mcode")
mun_100km <- st_as_sf(mun_100km)

# Download a background map
bb <- bb(municipalities)
background_map <- tmaptools::read_osm(bb, ext = 1.05, type = "esri-topo", zoom = 6)

# Create the map
tmap_mode(mode ="plot")
main_map <- tm_shape(background_map) +
  tm_rgb() + 
  tm_shape(municipalities) +
  tm_polygons(alpha = 0, border.col ="grey60") +
  tm_shape(mun_100km) +
  tm_polygons("treated.f", alpha = 0.25, palette = c("#386cb0", "#f87f01"), n = 2,  legend.show = FALSE) +
  tm_shape(d) +
  tm_dots("treated", size = .08, title = "Municipality capitals", palette = c("#386cb0", "#f87f01"), n = 2, labels = c("Outside biome, 0", "Inside, 1")) +
  tm_shape(biome_line) +
  tm_lines("LEGEND", palette = "#ef3b2c", lwd = 2, title.col = "Biome border", labels = "") +
  tm_shape(bla) +
  tm_polygons("SPRCLASSE", alpha = 0, border.col = "black", title = "Brazilian Legal Amazon", labels = "") +
  tm_shape(line_buffer_100km) +
  tm_polygons("ID", alpha = 0, border.col = "#984ea3", lwd = 1.5, title = "100km buffer", labels = "") +
  tm_layout(legend.bg.color = "white", frame = F, legend.frame = T, legend.position = c("left", "top")) +
  tm_compass(type = "arrow", position = c("right", "top")) +
  tm_scale_bar(breaks = c(0, 100, 200), position = c("right", "top")) #, text.size = 1)

# view map
main_map

# tmap_save(main_map, filename="bld/figures/figure1.pdf")#, height=8.5, width=11, units="in", dpi=300)

# Create South America minimap
world <- ne_countries(continent = "South America",  scale = "medium", returnclass = "sf")
world_minimap <- tm_shape(world) +
  tm_polygons() +
  tm_shape(st_as_sfc(bb)) +
  tm_borders(col = "brown3", lwd = 2) +
  tm_layout(bg.color = "lightblue") +
  tm_shape(municipalities) +
  tm_polygons(alpha = 0, border.col ="grey60") +
  tm_shape(bla) +
  tm_polygons("SPRCLASSE", alpha = 0, border.col = "black", legend.show = FALSE) +
  tm_shape(biome_line) +
  tm_lines("LEGEND", palette = "#ef3b2c", lwd = 1.5, legend.col.show = FALSE,
           legend.lwd.show = FALSE)

# coordinates for minimap
vp_world = viewport(0.88, 0.13, width = 0.25, height = 0.25)

# save map
tmap_save(main_map, "bld/figures/figure1.pdf",
          insets_tm = world_minimap, insets_vp = vp_world)

# Create segments map
main_map <- tm_shape(background_map) +
  tm_rgb() + #m tm_shape(lonmap, bbox = bhu_bb)
  tm_shape(municipalities) +
  tm_polygons(alpha = 0, border.col ="grey60") +
  tm_shape(d) +
  tm_dots("segments.f", size = .3, palette = "Set2", title = "Border segments") +
  tm_shape(biome_line) +
  tm_lines("LEGEND", palette = "#ef3b2c", lwd = 2, title.col = "Biome border", labels = "") +
  tm_shape(bla) +
  tm_polygons("SPRCLASSE", alpha = 0, border.col = "black", title = "Brazilian Legal Amazon", labels = "") +
  tm_layout(legend.bg.color = "white", frame = F, legend.frame = T, legend.position = c("left", "top")) +
  tm_compass(type = "arrow", position = c("right", "top")) +
  tm_scale_bar(breaks = c(0, 100, 200), position = c("right", "top")) #, text.size = 1)

# view map
main_map

# save map
tmap_save(main_map, "bld/figures/segments.pdf",
          insets_tm = world_minimap, insets_vp = vp_world)