Interactive visualization of forest cover losses

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:

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

Deforestation map

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.

basemap_custom <-
  leaflet() |>
  # add external map providers
  addTiles(group = "OpenStreetMap") |>
  addProviderTiles(providers$Esri.WorldImagery, group="Satellite") |>
  addProviderTiles(providers$CartoDB.Positron, group="CartoDB") |>
  addProviderTiles(providers$Esri.WorldShadedRelief, group="Topography") |>
  addProviderTiles(providers$NASAGIBS.ViirsEarthAtNight2012, group="Nightlights") |>
  addTiles(
    "https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png",
    group = "Forest Cover Loss (2001-2020)",
    attribution = "Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. “High-Resolution Global Maps of 21st-Century Forest Cover Change.” Science 342 (15 November): 850–53. Data available on-line from: http://earthenginepartners.appspot.com/science-2013-global-forest."
  ) |>
  addTiles(
    "https://tiles.globalforestwatch.org/umd_regional_primary_forest_2001/latest/dynamic/{z}/{x}/{y}.png",
    group = "Regional Primary Forests (2001)",
    attribution = "Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. “High-Resolution Global Maps of 21st-Century Forest Cover Change.” Science 342 (15 November): 850–53. Data available on-line from: http://earthenginepartners.appspot.com/science-2013-global-forest."
  ) |>
  addFullscreenControl() |>
  # add legend(s)
  addLayersControl(
    baseGroups = c("Satellite", "CartoDB", "OpenStreetMap", "Topography", "Nightlights"),
    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 = 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")

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.

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)