From b4aceec353266485e364403087d605eeec2a1f12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Darius=20A=2E=20G=C3=B6rgen?= Date: Fri, 26 Jul 2024 10:06:12 +0200 Subject: [PATCH] Adds bar plot and simpler fetching of WDPA IDs --- .gitignore | 10 +- interactive-deforestation-map.qmd | 160 ++++++++++++++++++++---------- 2 files changed, 113 insertions(+), 57 deletions(-) diff --git a/.gitignore b/.gitignore index 327a4f4..dfa6258 100644 --- a/.gitignore +++ b/.gitignore @@ -3,5 +3,11 @@ .Rhistory .Rproj.user .Ruserdata -WDPA/* -interactive-deforestation-map_files/* +*.html +<<<<<<< HEAD +======= + +/.quarto/ +/docs/ +interactive-deforestation-map_files +>>>>>>> 9ca3a4a (Uses plotly for barplots) diff --git a/interactive-deforestation-map.qmd b/interactive-deforestation-map.qmd index 336e250..e318e25 100644 --- a/interactive-deforestation-map.qmd +++ b/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) +```