Adds bar plot and simpler fetching of WDPA IDs

.gitignore

@@ -3,5 +3,8 @@
 .Rhistory
 .Rproj.user
 .Ruserdata
-WDPA/*
-interactive-deforestation-map_files/*
+*.html
+
+/.quarto/
+/docs/
+interactive-deforestation-map_files

interactive-deforestation-map.qmd

@@ -1,79 +1,87 @@
 ---
-title: "Interactive Deforestation Map"
+title: "Interactive visualization of forest cover losses"
 format:
   html:
     code-fold: true
     fig-width: 10
     fig-height: 10
+    embed-resources: true
 editor: source
+params:
+  wdpa_id: 115772 
+  buffer: 10000
+  lwd_pa: 8.0
+  col_pa: "#f18e26"
+  lwd_bf: 6.0
+  col_bf: "#1d7990"
+  gfw_version: "GFC-2023-v1.11"
+  min_size: 1
+  min_cover: 30
+  stacked: true
 ---
 
 ```{r setup}
 #| output: false
-
-library(dplyr)
-library(leaflet)
+#| echo: false
+library(mapme.biodiversity)
 library(leaflet.extras)
+library(leaflet)
+library(ggplot2)
+library(plotly)
 library(sf)
-library(wdpar)
+stopifnot(length(params$wdpa_id) == 1)
+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. You can click on
-the folded `Code` chunks to explore the source code of this report.
+This document can be used to visualize a single protected area from the World Database on Protected Areas (WDPA) on an interactive web map and calculate
+respective forest losses as indicated with Global Forest Watch (GFW).
+
+You can adjust the `params` object in the YAML header to your requirements:
+
+```yaml
+params:
+  wdpa_id: 115772 # ID of the protected area (PA) (see e.g. https://www.protectedplanet.net/115772)
+  buffer: 10000 # Size of the buffer zone in meters
+  lwd_pa: 8.0 # line width of the PA in pixels
+  col_pa: "#f18e26" # color for the PA in the map and barplot
+  lwd_bf: 6.0 # line width of the buffer zone
+  col_bf: "#1d7990" # color for the buffer zone in the map and barplot
+  gfw_version: "GFC-2023-v1.11" # version of GFW to use
+  min_size: 1 # minimum size of forest patches to be included (in ha)
+  min_cover: 30 # minimun percentage of vegetation cover to be considered forest
+  stacked: true # logical, indicating if bars are stacked in the plot
+```
 
-```{r config}
-#| code-fold: false
 
-######## Configuration #########################################################
-################################################################################
+## Deforestation map
 
-# ISO3 code of target country for which to download WDPA data.
-country_iso3 <- "BRA"
+The map displays the selected protected areas (PAs) including a buffer region around it. There are two additional layers included, one shwowing where forest loss
+occured since the year 2000 and the other indicating likley primary forests in the
+year 2001.
 
-# Provide specific WDPA IDs to be displayed here. If left empty, all PAs within the target country will be displayed.
-# Caution: The provided IDs must be located within the selected target country.
-wdpa_ids <- c("115772")
 
-# Buffer around the PAs in m.
-buffer_m <- 10000
+```{r funs}
+#| echo: false
+fetch_wdpaid <- function(id) {
 
-# Basic style options can be set here. For more elaborate style changes, you can edit the leaflet parameters below.
-line_weight_pa <- 4.0
-line_weight_buffer <- 3.0
-line_color_pa = "#CAE64E" # Color provided as hex code
-line_color_buffer = "#CAE64E" # Color provided as hex code
+  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"
+  read_sf(sprintf(url, id))
 
-###############################################################################
-###############################################################################
-```
+}
 
-# Preprocessing WDPA Data
+```
 
 ```{r process-wdpa}
-wdir <- getwd()
-wdpa <- wdpa_fetch(country_iso3, wait = TRUE, download_dir = file.path(wdir, "WDPA"))
-if(!(length(wdpa_ids) == 1 && is.na(wdpa_ids) || all(is.na(wdpa_ids)))) {
-  wdpa <- wdpa |>
-    filter(WDPAID %in% wdpa_ids)
-}
-
-wdpa <- wdpa |>
-  wdpa_clean(retain_status = NULL,
-             erase_overlaps = FALSE,
-             exclude_unesco = FALSE,
-             verbose = FALSE) |>
-  # Drop geometry POINT
-  filter(GEOMETRY_TYPE != "POINT") |>
-  # Remove the PAs that are only proposed, or have geometry type "point"
-  filter(STATUS != "Proposed") |>
-  st_transform(4326)
-
-# Make Buffers around all protected areas
-buffer <- wdpa |>
-  st_buffer(dist = buffer_m)
+#| echo: false
+#| warning: false
+
+wdpa <- fetch_wdpaid(params$wdpa_id)
+if (any(wdpa$status == "Proposed")) wdpa <- subset(wdpa, status != "Proposed")
+buffer <- st_buffer(wdpa, dist = params$buffer)
+st_geometry(buffer) <- sf::st_difference(st_geometry(buffer), st_geometry(wdpa))
+ 
 ```
 
-# Preparing Leaflet Map
 
 ```{r prepare-map}
 basemap_custom <-
@@ -98,19 +106,61 @@ 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-2023)", "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 = line_color_pa, smoothFactor = 0, weight = line_weight_pa, group = "WDPA Regions") |>
-  addPolygons(data = buffer, fillOpacity = 0.0, opacity = 0.8, color = line_color_buffer, smoothFactor = 0, weight = line_weight_buffer, 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")
 
-# Interactive Deforestation Map
+```
 
 ```{r plot-map}
+#| echo: false
   deforestation_map
 ```
+
+
+## Deforestation over time
+
+The barplot shows the amount of forest cover loss occured
+over time. It differentiates between losses that occured
+within the PA and the losses that occured in its respective 
+buffer zone. 
+
+```{r calc_losses}
+#| echo: false
+#| warning: false
+mapme_options(outdir = NULL)
+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 <- round(abs(c(
+  c(0, diff(subset(inds, buffer == "Protected Area")[["value"]])),
+  c(0, diff(subset(inds, buffer == "Buffer Zone")[["value"]]))
+  )))
+```
+
+
+```{r plot_losses}
+
+plt <- ggplot(data = inds) +
+  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=14),
+    axis.text=element_text(size=12),
+    axis.text.x=element_text(angle=60, hjust=1)
+    ) 
+
+ggplotly(plt)
+```