Add buffer zone stats to plot and make self-contained

.gitignore

@@ -3,5 +3,4 @@
 .Rhistory
 .Rproj.user
 .Ruserdata
-interactive-deforestation-map_files/*
 *.html

interactive-deforestation-map.qmd

@@ -5,18 +5,19 @@ format:
     code-fold: true
     fig-width: 10
     fig-height: 10
+    embed-resources: true
 editor: source
 params:
   wdpa_id: 115772
   buffer: 10000
-  lwd_pa: 4.0
-  col_pa: "#CAE64E"
-  lwd_bf: 3.0
-  col_bf: "#CAE64E"
+  lwd_pa: 8.0
+  col_pa: "#f18e26"
+  lwd_bf: 6.0
+  col_bf: "#032733"
   gfw_version: "GFC-2023-v1.11"
   min_size: 1
   min_cover: 30
-  col_bars: "#032733"
+  stacked: TRUE
 ---
 
 ```{r setup}
@@ -28,7 +29,7 @@ library(leaflet)
 library(ggplot2)
 library(sf)
 stopifnot(length(params$wdpa_id) == 1)
-#params <- list(wdpa_id = 115772, buffer = 10000, gfw_version = "GFC-2023-v1.11", min_size = 1, min_cover = 30, col_pa ="#032733")
+stopifnot(all(c("wdpa_id", "buffer", "lwd_pa", "col_pa", "lwd_bf", "col_bf", "gfw_version", "min_size", "min_cover", "stacked") %in% names(params)))
 ```
 
 This script downloads WDPA data for a specified country and generates a web map displaying the selected protected areas (PAs) including a buffer region around it.
@@ -39,18 +40,20 @@ Please adjust the `params` object in the YAML header to your requirements.
 fetch_wdpaid <- function(id) {
 
   url <- "https://data-gis.unep-wcmc.org/server/rest/services/ProtectedSites/The_World_Database_of_Protected_Areas/FeatureServer/1/query?where=WDPAID+=+%s&geometryType=esriGeometryPolygon&outFields=*&f=geojson"
-  sf::read_sf(sprintf(url, id))
+  read_sf(sprintf(url, id))
 
 }
 
 ```
 
 ```{r process-wdpa}
 #| echo: false
+#| warning: false
 
 wdpa <- fetch_wdpaid(params$wdpa_id)
 if (any(wdpa$status == "Proposed")) wdpa <- subset(wdpa, status != "Proposed")
-buffer <- sf::st_buffer(wdpa, dist = params$buffer)
+buffer <- st_buffer(wdpa, dist = params$buffer)
+st_geometry(buffer) <- sf::st_difference(st_geometry(buffer), st_geometry(wdpa))
  
 ```
 
@@ -78,15 +81,16 @@ basemap_custom <-
   # add legend(s)
   addLayersControl(
     baseGroups = c("Satellite", "CartoDB", "OpenStreetMap", "Topography", "Nightlights"),
-    overlayGroups = c("WDPA Regions", "Buffers", "Forest Cover Loss (2001-2020)", "Regional Primary Forests (2001)"),
+    overlayGroups = c("Protected Area", "Buffer Zone", "Forest Cover Loss (2001-2020)", "Regional Primary Forests (2001)"),
     options = layersControlOptions(collapsed = FALSE)) |>
   # uncheck some layers in layer control
   hideGroup(group = c("Regional Primary Forests (2001)","Labels (PA Names)"))
 
 deforestation_map <- basemap_custom |>
-  addPolygons(data = wdpa, fillOpacity = 0.0, opacity = 0.8, color = params$col_pa, smoothFactor = 0, weight = params$lwd_pa, group = "WDPA Regions") |>
-  addPolygons(data = buffer, fillOpacity = 0.0, opacity = 0.8, color = params$col_bf, smoothFactor = 0, weight = params$lwd_bf, group = "Buffers",
-              dashArray = "10 8")
+  addPolygons(data = buffer, fillOpacity = 0.0, opacity = 0.8, color = params$col_bf, smoothFactor = 0, weight =params$lwd_bf, group = "Buffer Zone",
+              dashArray = "10 8") |>
+  addPolygons(data = wdpa, fillOpacity = 0.0, opacity = 1.0, color = params$col_pa, smoothFactor = 0, weight = params$lwd_pa, group = "Protected Area")
+
 ```
 
 ```{r plot-map}
@@ -98,11 +102,16 @@ deforestation_map <- basemap_custom |>
 ```{r calc_losses}
 #| echo: false
 mapme_options(outdir = NULL)
-wdpa <- get_resources(wdpa, get_gfw_treecover(params$gfw_version), get_gfw_lossyear(params$gfw_version))
-wdpa <- calc_indicators(wdpa, calc_treecover_area(min_size = params$min_size, min_cover = params$min_cover))
-inds <- wdpa[ ,c("wdpaid", "treecover_area")]
+inds <- rbind(wdpa, buffer)
+inds$buffer <- factor(0:1, levels = 0:1, labels = c("Protected Area", "Buffer Zone"))
+inds <- get_resources(inds, get_gfw_treecover(params$gfw_version), get_gfw_lossyear(params$gfw_version))
+inds <- calc_indicators(inds, calc_treecover_area(min_size = params$min_size, min_cover = params$min_cover))
+inds <- inds[ ,c("wdpaid", "buffer", "treecover_area")]
 inds <- portfolio_long(inds, drop_geoms = TRUE)
-inds$losses <- c(0, abs(diff(inds$value)))
+inds$losses <- abs(c(
+  c(0, diff(subset(inds, buffer == "Protected Area")[["value"]])),
+  c(0, diff(subset(inds, buffer == "Buffer Zone")[["value"]]))
+  ))
 ```
 
 
@@ -111,14 +120,16 @@ inds$losses <- c(0, abs(diff(inds$value)))
 ```{r plot_losses}
 
 ggplot(data = inds) +
-  geom_col(aes(x=datetime, y=losses), fill=params$col_bars) +
-  labs(x = "Year", y = "Forest cover losses (in ha)") +
+  geom_col(aes(x=datetime, y=losses, fill=buffer), position=ifelse(params$stacked, "stack", "dodge")) +
+  labs(x = "Year", y = "Forest cover losses (in ha)", fill = "Zone") +
   scale_x_datetime(date_breaks = "1 year", date_labels =  "%Y") +
+  scale_fill_manual(values = c(params$col_pa, params$col_bf)) +
   theme_classic() +
   theme(
     axis.title=element_text(size=16),
     axis.text=element_text(size=14),
-    axis.text.x=element_text(angle=60, hjust=1)
+    axis.text.x=element_text(angle=60, hjust=1),
+    legend.position="bottom"
     )
 
 ```