---
title: "Evaluation of Forest Cover Development in Protected Areas in Laos"
subtitle: "Analysis of Forest Cover Loss in Four Protected Areas (2000-2024)"
author:
- name: "Johannes Schielein"
affiliations: "KfW Development Bank"
date: "2026-01-12"
date-modified: last-modified
abstract: |
This evaluation report analyzes forest cover development in four protected areas in Laos,
examining forest cover loss dynamics from 2000 to 2024. The analysis utilizes Global Forest
Watch data to assess conservation outcomes and forest stock development, supporting evaluation
of protected area management effectiveness.
categories: [Laos, Protected Areas, Forest Loss, Evaluation, Conservation Outcomes]
bibliography: references.bib
---
```{css}
/*| echo: false */
.map-container {
border: 2px solid #cccccc;
border-radius: 8px;
margin: 0 auto;
}
```
```{r setup}
#| include: false
#| cache: false
# Load libraries
library(tidyverse)
library(sf)
library(leaflet)
library(leaflet.extras)
library(leaflet.extras2)
library(ggsci)
library(scales)
library(htmltools)
library(RColorBrewer)
library(plotly)
library(DT)
options(scipen = 10000)
# Define Laos WDPA IDs for analysis (project areas)
# Note: Nam Ha has two WDPA IDs (71252 = National Protected Area, 555756396 = ASEAN Heritage Park)
# We use only the National Protected Area designation (71252)
target_wdpa_ids <- c(71253, 71252, 555703744, 555703750)
# Define consistent color scheme for project areas
# Using a distinct, appealing color palette
project_colors <- c(
"Nam Kan" = "#1f77b4", # Blue
"Nam Ha" = "#ff7f0e", # Orange
"Hin Nam No" = "#2ca02c", # Green
"Kounxe Nongma" = "#d62728" # Red
)
```
```{r data-loading}
#| include: false
#| cache: false
# ----- Load Protected Areas Data -----
# Load Laos PAs from WDPA file geodatabase
wdpa_lao <- st_read(
"../data/laos/input/wdpa_lao.gdb",
layer = "WDPA_WDOECM_poly_Feb2026_LAO",
quiet = TRUE
)
# Rename columns to standard names
if ("SITE_ID" %in% names(wdpa_lao)) {
wdpa_lao <- wdpa_lao %>%
rename(WDPAID = SITE_ID)
} else if ("WDPAID" %in% names(wdpa_lao)) {
# Already has WDPAID
} else {
# Try to find the ID column
id_col <- grep("ID|id", names(wdpa_lao), value = TRUE)[1]
if (!is.na(id_col)) {
wdpa_lao <- wdpa_lao %>%
rename(WDPAID = !!sym(id_col))
}
}
if ("NAME" %in% names(wdpa_lao) && !"ORIG_NAME" %in% names(wdpa_lao)) {
wdpa_lao <- wdpa_lao %>%
rename(ORIG_NAME = NAME)
}
# Filter for target areas (project areas)
# Note: Exclude Nam Ha ASEAN Heritage Park (555756396), keep only National Protected Area (71252)
target_areas <- wdpa_lao %>%
filter(WDPAID %in% target_wdpa_ids) %>%
filter(!(WDPAID == 555756396)) %>% # Explicitly exclude ASEAN Heritage Park
st_make_valid()
# Transform to WGS84 if needed
if (st_crs(target_areas)$epsg != 4326) {
target_areas <- st_transform(target_areas, 4326)
}
# Get total area from GIS_M_AREA or GIS_AREA (convert from sqkm to ha)
if ("GIS_M_AREA" %in% names(target_areas)) {
target_areas$total_area_ha <- target_areas$GIS_M_AREA * 100
} else if ("GIS_AREA" %in% names(target_areas)) {
target_areas$total_area_ha <- target_areas$GIS_AREA * 100
} else if ("REP_AREA" %in% names(target_areas)) {
target_areas$total_area_ha <- target_areas$REP_AREA * 100
} else {
target_areas$total_area_ha <- as.numeric(st_area(target_areas)) / 10000
}
cat("✓ Loaded", nrow(target_areas), "target protected areas\n")
# ----- Load Forest Cover Data -----
if (file.exists("../data/laos/output/laos_forest_cover_2000_2024.csv")) {
forest_data <- read_csv(
"../data/laos/output/laos_forest_cover_2000_2024.csv",
show_col_types = FALSE
)
# Remove duplicate and cross-border entries:
# - Nam Ha ASEAN Heritage Park (555756396) - duplicate of 71252
# - Nam Poui (555756393) - duplicate of 10132
# - Phou Xieng Thong (555756394) - duplicate of 18893
# - Phong Nha-Ke Bang / Hin Nam No (900883) - cross-border UNESCO site spanning Laos and Vietnam
forest_data <- forest_data %>%
filter(!WDPAID %in% c(555756396, 555756393, 555756394, 900883)) %>%
mutate(WDPAID = as.numeric(WDPAID))
# Separate project and non-project areas
project_data <- forest_data %>% filter(is_project_area == TRUE)
other_data <- forest_data %>% filter(is_project_area == FALSE)
cat("✓ Loaded forest cover data:", nrow(forest_data), "records\n")
cat("✓ Project areas:", length(unique(project_data$WDPAID)), "PAs\n")
cat("✓ Other areas:", length(unique(other_data$WDPAID)), "PAs\n")
} else {
warning("Forest cover data not found. Run process_laos_forest_cover.R first!")
forest_data <- tibble(
WDPAID = integer(),
ORIG_NAME = character(),
year = integer(),
area_ha = numeric(),
loss_ha = numeric(),
cumulative_loss_ha = numeric()
)
}
# ----- Load All Lao PAs for Map and Comparison -----
wdpa_all_lao <- wdpa_lao %>%
st_make_valid()
if (st_crs(wdpa_all_lao)$epsg != 4326) {
wdpa_all_lao <- st_transform(wdpa_all_lao, 4326)
}
# Filter for non-project areas (exclude all project areas)
# Also exclude duplicates and cross-border entries
non_project_areas <- wdpa_all_lao %>%
filter(!WDPAID %in% target_wdpa_ids) %>%
filter(!WDPAID %in% c(555756396, 555756393, 555756394, 900883)) %>%
filter(!grepl("Hin Nam Ho", ORIG_NAME, ignore.case = TRUE)) %>%
filter(!grepl("Hin Nam No", ORIG_NAME, ignore.case = TRUE))
cat("✓ Loaded all Lao PAs:", nrow(wdpa_all_lao), "total\n")
cat("✓ Non-project areas:", nrow(non_project_areas), "\n")
```
## Introduction
This evaluation report analyzes forest cover development in four protected areas in Laos to support the evaluation of conservation outcomes and forest stock development. The analysis focuses on forest cover loss dynamics from 2000 to 2024, utilizing publicly available satellite data to assess trends and patterns in forest conservation.
### Purpose of the Evaluation
The evaluation aims to:
- Assess forest cover development trends in four target protected areas
- Compare deforestation rates in project areas against other national protected areas
- Provide empirical evidence for evaluating conservation outcomes
- Support decision-making for future conservation investments
### Data Sources and Methodology
This analysis utilizes:
- **World Database of Protected Areas (WDPA/IUCN)** - For protected area boundaries and metadata
- **Global Forest Watch (Hansen et al., 2013)** - For forest cover loss detection using GFC-2024-v1.12
- **mapme.biodiversity R package** - For automated data acquisition and processing
**Methodology Notes:**
- Forest cover loss is defined as "a stand-replacement disturbance, or a change from a forest to non-forest state"
- Forest cover is defined as tree cover ≥10% with minimum patch size of 1 hectare
- Analysis period: 2000-2024 (25 years)
- Cumulative deforestation is calculated from the year 2001 onward
::: {.callout-note}
## Project Timeline
The financing agreement for the project was signed in **2015**. However, implementation activities started effectively in **2018** with participatory land use planning. The main implementation phase lasted from **2018 to 2023**, during which village contracts were signed, forest patrols were conducted, and microcredit programs were established.
Note: Hin Nam No National Park was only supported through ancillary measures and is sometimes not listed as a formal project site. Where relevant, averages are shown both with and without Hin Nam No.
:::
## Analyzed Protected Areas
The evaluation focuses on four protected areas in Laos that are part of the evaluation portfolio:
::: {.callout-note}
## Notes on Duplicate and Cross-Border Entries
Several protected areas appear multiple times in the WDPA database. The following duplicates and cross-border entries have been excluded from this analysis:
- **Nam Ha ASEAN Heritage Park** (WDPA ID 555756396) — duplicate of Nam Ha National Protected Area (71252)
- **Nam Poui** (WDPA ID 555756393) — duplicate of Nam Poui (10132)
- **Phou Xieng Thong** (WDPA ID 555756394) — duplicate of Phou Xieng Thong (18893)
- **Phong Nha-Ke Bang / Hin Nam No** (WDPA ID 900883) — cross-border UNESCO World Heritage site spanning Laos and Vietnam (not comparable to national-level PAs)
:::
```{r table-pas}
#| echo: false
#| fig-height: 8
# Get forest cover data for 2000 and latest year
if (nrow(forest_data) > 0) {
forest_2000 <- forest_data %>%
filter(year == 2000) %>%
select(WDPAID, forest_2000 = area_ha)
latest_year <- max(forest_data$year, na.rm = TRUE)
forest_latest <- forest_data %>%
filter(year == latest_year) %>%
select(WDPAID, forest_latest = area_ha)
} else {
forest_2000 <- tibble(WDPAID = numeric(), forest_2000 = numeric())
forest_latest <- tibble(WDPAID = numeric(), forest_latest = numeric())
latest_year <- 2024
}
# Create table data
table_data <- target_areas %>%
st_drop_geometry() %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
left_join(forest_2000, by = "WDPAID") %>%
left_join(forest_latest, by = "WDPAID")
# Add designation and status columns if they exist, otherwise use "N/A"
if ("DESIG_ENG" %in% names(table_data)) {
table_data$Designation <- table_data$DESIG_ENG
} else {
table_data$Designation <- "N/A"
}
if ("STATUS" %in% names(table_data)) {
table_data$Status <- table_data$STATUS
} else {
table_data$Status <- "N/A"
}
# Select and rename columns
table_data <- table_data %>%
mutate(
`Forest {latest_year} (ha)` = forest_latest
) %>%
select(
`Name` = ORIG_NAME,
`WDPA ID` = WDPAID,
`Area (ha)` = total_area_ha,
`Forest 2000 (ha)` = forest_2000,
`Forest {latest_year} (ha)`,
Designation,
Status
) %>%
arrange(`Name`)
# Fix column name
names(table_data)[names(table_data) == "Forest {latest_year} (ha)"] <- paste0("Forest ", latest_year, " (ha)")
table_data %>%
datatable(
filter = 'top',
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
scrollX = TRUE,
pageLength = 10,
buttons = list(
list(extend = 'csv', text = 'Download CSV'),
list(extend = 'excel', text = 'Download Excel')
)
),
caption = "Protected Areas analyzed in this evaluation report"
) %>%
formatRound(c('Area (ha)', 'Forest 2000 (ha)', paste0('Forest ', latest_year, ' (ha)')), 0)
```
## Map
```{r map-laos}
#| echo: false
#| fig-height: 3.5
#| out-width: "100%"
#| classes: map-container
# Prepare data for mapping - calculate forest loss statistics for ALL areas
if (nrow(forest_data) > 0) {
# Calculate absolute and relative loss for ALL areas
forest_stats <- forest_data %>%
group_by(WDPAID, ORIG_NAME) %>%
summarise(
forest_2000 = area_ha[year == 2000],
forest_latest = area_ha[year == max(year, na.rm = TRUE)],
absolute_loss = forest_2000 - forest_latest,
relative_loss = (absolute_loss / forest_2000) * 100,
is_project_area = first(is_project_area),
.groups = "drop"
)
# Join with all areas from WDPA
all_areas_stats <- wdpa_all_lao %>%
st_drop_geometry() %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
left_join(forest_stats, by = "WDPAID")
# Get total area
if ("GIS_M_AREA" %in% names(all_areas_stats)) {
all_areas_stats$total_area_ha <- all_areas_stats$GIS_M_AREA * 100
} else if ("GIS_AREA" %in% names(all_areas_stats)) {
all_areas_stats$total_area_ha <- all_areas_stats$GIS_AREA * 100
} else if ("REP_AREA" %in% names(all_areas_stats)) {
all_areas_stats$total_area_ha <- all_areas_stats$REP_AREA * 100
} else {
# Calculate from geometry
all_areas_with_geom <- wdpa_all_lao %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
left_join(forest_stats, by = "WDPAID")
all_areas_stats$total_area_ha <- as.numeric(st_area(all_areas_with_geom)) / 10000
}
# Categorize losses for all areas
all_areas_stats$relative_loss_categorized <- cut(
all_areas_stats$relative_loss,
breaks = c(-1, 2, 5, 10, 20, 100),
labels = c("0-2%", "2-5%", "5-10%", "10-20%", ">20%")
)
max_abs_loss <- max(all_areas_stats$absolute_loss, na.rm = TRUE)
all_areas_stats$absolute_loss_categorized <- cut(
all_areas_stats$absolute_loss,
breaks = c(-1, 100, 1000, 5000, 10000, max_abs_loss),
labels = c("0-100 ha", "100-1,000 ha", "1,000-5,000 ha", "5,000-10,000 ha", ">10,000 ha")
)
# Get centroids for markers for all areas
all_areas_with_geom <- wdpa_all_lao %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
left_join(forest_stats, by = "WDPAID")
all_centroids <- st_centroid(all_areas_with_geom)
all_coords <- st_coordinates(all_centroids)
all_areas_stats$X <- all_coords[, "X"]
all_areas_stats$Y <- all_coords[, "Y"]
# Create HTML labels for all areas
all_areas_stats$label <- paste(
"<p><b>", all_areas_stats$ORIG_NAME, "</b></br>",
"Total Area:", format(round(all_areas_stats$total_area_ha, 0), big.mark = ","), "ha</br>",
"Forest Cover (2000):", format(round(all_areas_stats$forest_2000, 0), big.mark = ","), "ha</br>",
"Absolute Loss (2000-2024):", format(round(all_areas_stats$absolute_loss, 0), big.mark = ","), "ha</br>",
"Relative Loss (2000-2024):", round(all_areas_stats$relative_loss, 2), "%",
"</p>"
)
# Filter out areas without data
all_areas_stats <- all_areas_stats %>%
filter(!is.na(forest_2000))
} else {
all_areas_stats <- wdpa_all_lao %>%
st_drop_geometry() %>%
mutate(
X = NA, Y = NA,
label = paste0("<p><b>", ORIG_NAME, "</b></p>"),
relative_loss_categorized = "0-2%",
absolute_loss_categorized = "0-100 ha"
)
}
# Prepare label data for all PAs (offset from centroid)
all_pas_for_labels <- bind_rows(
target_areas %>% mutate(pa_type = "Project"),
non_project_areas %>% mutate(pa_type = "Other")
) %>%
filter(!is.na(ORIG_NAME))
# Calculate centroids for labels
if (nrow(all_pas_for_labels) > 0) {
label_centroids <- st_centroid(all_pas_for_labels)
label_coords <- st_coordinates(label_centroids)
all_pas_for_labels$label_lng <- label_coords[, "X"] + 0.15 # Offset east
all_pas_for_labels$label_lat <- label_coords[, "Y"] + 0.08 # Offset north
}
# Color palettes
pal_relative_loss <- colorFactor(
palette = rev(brewer.pal(5, "RdYlGn")),
domain = all_areas_stats$relative_loss_categorized
)
# Create map - build conditionally
map_attribution <- paste(
"Forest Loss: Hansen et al. (2013) |",
"Protected Areas: WDPA/IUCN (2024)"
)
map <- leaflet() %>%
addTiles(group = "OpenStreetMap") %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
addProviderTiles(providers$Esri.WorldShadedRelief, group = "Topography") %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2009&end_year=2013&render_type=true_color",
group = "Forest Cover Loss (2009-2013)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2014&end_year=2018&render_type=true_color",
group = "Forest Cover Loss (2014-2018)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2019&end_year=2023&render_type=true_color",
group = "Forest Cover Loss (2019-2023)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2001&end_year=2024&render_type=true_color",
group = "Forest Cover Loss (2001-2024)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_regional_primary_forest_2001/latest/dynamic/{z}/{x}/{y}.png",
group = "Regional Primary Forests (2001)",
attribution = map_attribution
)
# Add polygon layers
if (nrow(target_areas) > 0 && inherits(target_areas, "sf")) {
# Project PAs outlines only (red outline, transparent fill)
map <- map %>% addPolygons(
data = target_areas,
fillColor = "transparent",
fillOpacity = 0,
color = "red",
weight = 2,
opacity = 1,
group = "Project PAs",
label = ~ORIG_NAME
)
}
if (nrow(non_project_areas) > 0 && inherits(non_project_areas, "sf")) {
map <- map %>% addPolygons(
data = non_project_areas,
fillColor = "#90EE90",
fillOpacity = 0.2,
color = "#006400",
weight = 2,
opacity = 0.9,
group = "Other PAs",
label = ~ORIG_NAME
)
}
# Add circle markers for forest loss statistics (ALL areas)
if (nrow(all_areas_stats) > 0 && !all(is.na(all_areas_stats$X))) {
map <- map %>%
addCircleMarkers(
data = all_areas_stats,
lng = ~X,
lat = ~Y,
color = ~pal_relative_loss(relative_loss_categorized),
group = "Forest Loss Statistics",
radius = ~ifelse(absolute_loss_categorized == ">10,000 ha", 12,
ifelse(absolute_loss_categorized == "5,000-10,000 ha", 8,
ifelse(absolute_loss_categorized == "1,000-5,000 ha", 6,
ifelse(absolute_loss_categorized == "100-1,000 ha", 4, 2)))),
popup = ~label,
stroke = FALSE,
fillOpacity = 0.7
)
}
map <- map %>%
addFullscreenControl() %>%
addLegend(
"bottomright",
data = all_areas_stats,
pal = pal_relative_loss,
values = ~relative_loss_categorized,
title = "Relative Forest Cover Loss<br>(2000-2024)",
opacity = 1
) %>%
addLayersControl(
baseGroups = c("CartoDB", "OpenStreetMap", "Satellite", "Topography"),
overlayGroups = c("Project PAs", "Other PAs", "PA Labels",
"Forest Loss Statistics",
"Forest Cover Loss (2009-2013)",
"Forest Cover Loss (2014-2018)",
"Forest Cover Loss (2019-2023)",
"Forest Cover Loss (2001-2024)",
"Regional Primary Forests (2001)"),
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(c("Forest Cover Loss (2009-2013)", "Forest Cover Loss (2014-2018)",
"Forest Cover Loss (2019-2023)", "Forest Cover Loss (2001-2024)",
"Regional Primary Forests (2001)"))
# Note: Layer opacity control in Leaflet requires custom JavaScript implementation
# The layer control allows toggling layers on/off. For dynamic opacity sliders,
# additional JavaScript would need to be added to manipulate layer opacity.
# Add PA name labels (offset from centroid)
if (nrow(all_pas_for_labels) > 0) {
# Project PA labels (dark green, bold)
project_labels <- all_pas_for_labels %>% filter(pa_type == "Project")
if (nrow(project_labels) > 0) {
map <- map %>% addLabelOnlyMarkers(
data = project_labels,
lng = ~label_lng,
lat = ~label_lat,
label = ~ORIG_NAME,
labelOptions = labelOptions(
noHide = TRUE,
direction = "right",
textOnly = TRUE,
style = list(
"color" = "#006400",
"font-weight" = "bold",
"font-size" = "11px",
"text-shadow" = "1px 1px 2px white, -1px -1px 2px white"
)
),
group = "PA Labels"
)
}
# Other PA labels (dark green, normal)
other_labels <- all_pas_for_labels %>% filter(pa_type == "Other")
if (nrow(other_labels) > 0) {
map <- map %>% addLabelOnlyMarkers(
data = other_labels,
lng = ~label_lng,
lat = ~label_lat,
label = ~ORIG_NAME,
labelOptions = labelOptions(
noHide = TRUE,
direction = "right",
textOnly = TRUE,
style = list(
"color" = "#006400",
"font-weight" = "normal",
"font-size" = "10px",
"text-shadow" = "1px 1px 2px white, -1px -1px 2px white"
)
),
group = "PA Labels"
)
}
}
map
```
## Forest Cover Development (2000-Present)
Forest cover development is quantified using Global Forest Watch data (Hansen et al., 2013). This section presents overall trends in forest cover and loss across the analysis period.
### Summary Statistics
```{r summary-stats}
#| echo: false
if (nrow(forest_data) > 0 && nrow(wdpa_all_lao) > 0) {
# Get total area from WDPA data for ALL areas
# First calculate area from geometry if needed
if (!"GIS_M_AREA" %in% names(wdpa_all_lao) && !"GIS_AREA" %in% names(wdpa_all_lao) && !"REP_AREA" %in% names(wdpa_all_lao)) {
wdpa_all_lao$total_area_ha <- as.numeric(st_area(wdpa_all_lao)) / 10000
}
all_areas_df <- wdpa_all_lao %>%
st_drop_geometry() %>%
mutate(WDPAID = as.numeric(WDPAID))
# Get GIS_M_AREA or GIS_AREA for total area
if ("GIS_M_AREA" %in% names(all_areas_df)) {
all_areas_df$total_area_ha <- all_areas_df$GIS_M_AREA * 100
} else if ("GIS_AREA" %in% names(all_areas_df)) {
all_areas_df$total_area_ha <- all_areas_df$GIS_AREA * 100
} else if ("REP_AREA" %in% names(all_areas_df)) {
all_areas_df$total_area_ha <- all_areas_df$REP_AREA * 100
} else if (!"total_area_ha" %in% names(all_areas_df)) {
# Fallback: use a default or calculate from geometry before dropping
all_areas_df$total_area_ha <- NA_real_
}
# Get forest cover in 2000
forest_2000 <- forest_data %>%
filter(year == 2000) %>%
select(WDPAID, forest_2000 = area_ha) %>%
mutate(WDPAID = as.numeric(WDPAID))
# Calculate cumulative losses for periods
period_losses <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
group_by(WDPAID, ORIG_NAME) %>%
summarise(
cum_loss_2010_2015 = sum(loss_ha[year >= 2010 & year <= 2015], na.rm = TRUE),
cum_loss_2018_2023 = sum(loss_ha[year >= 2018 & year <= 2023], na.rm = TRUE),
cum_loss_total = sum(loss_ha[year >= 2001 & year <= 2024], na.rm = TRUE),
is_project_area = first(is_project_area),
.groups = "drop"
)
# Pivot annual losses to wide format
annual_losses <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2001 & year <= 2024) %>%
select(WDPAID, ORIG_NAME, year, loss_ha) %>%
pivot_wider(
names_from = year,
values_from = loss_ha,
names_prefix = "Loss "
)
# Combine all data
summary_wide <- all_areas_df %>%
select(WDPAID, ORIG_NAME, total_area_ha) %>%
left_join(forest_2000, by = "WDPAID") %>%
left_join(period_losses, by = "WDPAID") %>%
left_join(annual_losses, by = "WDPAID") %>%
mutate(
`Project Area` = ifelse(is_project_area, "Yes", "No"),
`Protected Area Name` = ORIG_NAME,
`WDPA ID` = WDPAID,
`Total Protected Area` = total_area_ha,
`Forest Cover in 2000` = forest_2000,
`Cumulative Loss Pre-Implementation (2010-2015) (ha)` = cum_loss_2010_2015,
`Cumulative Loss Implementation (2018-2023) (ha)` = cum_loss_2018_2023,
`Total Loss (2001-2024) (ha)` = cum_loss_total
) %>%
select(
`Protected Area Name`,
`WDPA ID`,
`Project Area`,
`Total Protected Area`,
`Forest Cover in 2000`,
`Cumulative Loss Pre-Implementation (2010-2015) (ha)`,
`Cumulative Loss Implementation (2018-2023) (ha)`,
`Total Loss (2001-2024) (ha)`,
starts_with("Loss ")
) %>%
arrange(`Project Area`, `Protected Area Name`)
# Get column names for formatting
loss_cols <- names(summary_wide)[grepl("^Loss ", names(summary_wide))]
format_cols <- c('Total Protected Area', 'Forest Cover in 2000',
'Cumulative Loss Pre-Implementation (2010-2015) (ha)',
'Cumulative Loss Implementation (2018-2023) (ha)',
'Total Loss (2001-2024) (ha)',
loss_cols)
summary_wide %>%
datatable(
filter = 'top',
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
scrollX = TRUE,
scrollY = "400px",
pageLength = 10,
buttons = list(
list(extend = 'csv', text = 'Download CSV'),
list(extend = 'excel', text = 'Download Excel')
)
),
caption = "Forest Cover Development Summary (2000-Present) - Wide Format"
) %>%
formatRound(format_cols, 0)
} else {
cat("**Note:** Summary table not yet generated. Run `process_laos_forest_cover.R` to create this data.")
}
```
### Overall Trends
```{r overall-trends}
#| echo: false
#| fig-height: 8
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
# Prepare data for stacked bar charts
# Forest cover by area and year (project areas only for stacking)
forest_by_area <- project_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
select(WDPAID, ORIG_NAME, year, area_ha) %>%
arrange(ORIG_NAME, year)
# Annual losses by area and year (project areas only for stacking)
loss_by_area <- project_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2001) %>%
select(WDPAID, ORIG_NAME, year, loss_ha) %>%
arrange(ORIG_NAME, year)
# Plot 1: Stacked bar chart for forest cover
p1 <- ggplot(forest_by_area, aes(x = year, y = area_ha, fill = ORIG_NAME)) +
geom_bar(stat = "identity", position = "stack") +
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(forest_by_area$area_ha) * 4 * 0.98, label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(forest_by_area$area_ha) * 4 * 0.98, label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_fill_manual(values = project_colors) +
labs(
x = "Year",
y = "Total Forest Cover (ha)",
title = "Total Forest Cover Development (2000-2024)",
subtitle = "Stacked by protected area",
fill = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2000, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom"
)
# Plot 2: Stacked bar chart for annual losses
p2 <- ggplot(loss_by_area, aes(x = year, y = loss_ha, fill = ORIG_NAME)) +
geom_bar(stat = "identity", position = "stack") +
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
scale_fill_manual(values = project_colors) +
labs(
x = "Year",
y = "Annual Forest Loss (ha)",
title = "Annual Forest Cover Loss (2001-2024)",
subtitle = "Stacked by protected area",
fill = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2001, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom"
)
# Make interactive - use separate plots to avoid legend duplication
p1_plotly <- ggplotly(p1, tooltip = c("year", "area_ha", "ORIG_NAME"))
p2_plotly <- ggplotly(p2, tooltip = c("year", "loss_ha", "ORIG_NAME"))
# Remove legend from second plot to avoid duplication
p2_plotly <- p2_plotly %>% layout(showlegend = FALSE)
# Combine with subplot
subplot(
p1_plotly,
p2_plotly,
nrows = 2,
shareX = TRUE,
shareY = FALSE,
margin = 0.05
) %>%
layout(
title = "Forest Cover Development Overview"
)
} else {
cat("**Note:** Forest cover data not yet available. Run `process_laos_forest_cover.R` to generate this data.")
}
```
## Temporal Analysis of Deforestation
This section presents cumulative deforestation trends and forest cover development for project areas and all analyzed protected areas.
### Cumulative Forest Cover Loss (2001-2024) for Project PAs
```{r graph-cumulative-project}
#| echo: false
#| fig-height: 3.5
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
# Calculate cumulative loss from 2001 baseline for project areas only
cumulative_data <- project_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2001) %>%
group_by(WDPAID, ORIG_NAME) %>%
arrange(year) %>%
mutate(
baseline_2001 = area_ha[year == 2001],
cumulative_loss = baseline_2001 - area_ha
) %>%
ungroup()
if (nrow(cumulative_data) > 0) {
p <- ggplot(cumulative_data, aes(
x = year,
y = cumulative_loss,
color = ORIG_NAME,
group = ORIG_NAME,
text = paste(
"Protected Area:", ORIG_NAME,
"\nYear:", year,
"\nCumulative Loss:", format(round(cumulative_loss, 0), big.mark = ","), "ha"
)
)) +
geom_line(linewidth = 1.2) +
geom_point(size = 2) +
# Add vertical line for financing agreement and shaded implementation phase
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(cumulative_data$cumulative_loss, na.rm = TRUE) * 0.95,
label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(cumulative_data$cumulative_loss, na.rm = TRUE) * 0.95,
label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_color_manual(values = project_colors) +
labs(
x = "Year",
y = "Cumulative Forest Loss (ha)",
title = "Cumulative Forest Cover Loss (2001-2024) for Project PAs",
subtitle = "Cumulative loss from 2001 baseline",
color = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2001, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom",
legend.key.width = unit(0.5, "cm"),
legend.text = element_text(size = 8),
legend.title = element_text(size = 9)
) +
guides(color = guide_legend(
nrow = 2,
override.aes = list(size = 1.5),
title.position = "top"
))
# Truncate long names in legend and limit legend width
p_plotly <- ggplotly(p, tooltip = "text")
# Truncate labels
for (i in 1:length(p_plotly$x$data)) {
if (!is.null(p_plotly$x$data[[i]]$name)) {
if (nchar(p_plotly$x$data[[i]]$name) > 25) {
p_plotly$x$data[[i]]$name <- paste0(substr(p_plotly$x$data[[i]]$name, 1, 22), "...")
}
}
}
p_plotly %>% layout(
legend = list(
orientation = "h",
x = 0.5,
xanchor = "center",
y = -0.15,
yanchor = "top",
font = list(size = 9),
xref = "paper",
yref = "paper"
),
margin = list(b = 100)
)
} else {
cat("**Note:** Data not available.")
}
} else {
cat("**Note:** Forest cover data not yet available. Run `process_laos_forest_cover.R` to generate this data.")
}
```
### Cumulative Forest Cover Loss (2001-2024) for All PAs
```{r graph-cumulative-all}
#| echo: false
#| fig-height: 3.5
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
# Combine project and other areas
cumulative_data <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2001) %>%
group_by(WDPAID, ORIG_NAME, is_project_area) %>%
arrange(year) %>%
mutate(
baseline_2001 = area_ha[year == 2001],
cumulative_loss = baseline_2001 - area_ha
) %>%
ungroup()
if (nrow(cumulative_data) > 0) {
# Create color mapping: project areas use project_colors, others use grey
all_areas <- unique(cumulative_data$ORIG_NAME)
color_mapping <- ifelse(all_areas %in% names(project_colors),
project_colors[all_areas],
"grey70")
names(color_mapping) <- all_areas
p <- ggplot(cumulative_data, aes(
x = year,
y = cumulative_loss,
color = ORIG_NAME,
group = ORIG_NAME,
alpha = is_project_area,
text = paste(
"Protected Area:", ORIG_NAME,
"\nYear:", year,
"\nCumulative Loss:", format(round(cumulative_loss, 0), big.mark = ","), "ha",
"\nProject Area:", ifelse(is_project_area, "Yes", "No")
)
)) +
geom_line(linewidth = ifelse(cumulative_data$is_project_area, 1.2, 0.8)) +
geom_point(size = ifelse(cumulative_data$is_project_area, 2, 1.5)) +
# Add vertical line for financing agreement and shaded implementation phase
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(cumulative_data$cumulative_loss, na.rm = TRUE) * 0.95,
label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(cumulative_data$cumulative_loss, na.rm = TRUE) * 0.95,
label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_color_manual(values = color_mapping) +
scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0.6), guide = "none") +
labs(
x = "Year",
y = "Cumulative Forest Loss (ha)",
title = "Cumulative Forest Cover Loss (2001-2024) for All PAs",
subtitle = "Project PAs shown in colors, other PAs in grey",
color = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2001, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom",
legend.key.width = unit(0.5, "cm"),
legend.text = element_text(size = 8),
legend.title = element_text(size = 9)
) +
guides(color = guide_legend(
nrow = 2,
override.aes = list(size = 1.5),
title.position = "top"
))
# Truncate long names in legend and limit legend width
p_plotly <- ggplotly(p, tooltip = "text")
# Truncate labels and adjust legend
for (i in 1:length(p_plotly$x$data)) {
if (!is.null(p_plotly$x$data[[i]]$name)) {
if (nchar(p_plotly$x$data[[i]]$name) > 25) {
p_plotly$x$data[[i]]$name <- paste0(substr(p_plotly$x$data[[i]]$name, 1, 22), "...")
}
}
}
p_plotly %>% layout(
legend = list(
orientation = "h",
x = 0.5,
xanchor = "center",
y = -0.15,
yanchor = "top",
font = list(size = 9),
xref = "paper",
yref = "paper"
),
margin = list(b = 100) # Add bottom margin for legend
)
} else {
cat("**Note:** Data not available.")
}
} else {
cat("**Note:** Forest cover data not yet available. Run `process_laos_forest_cover.R` to generate this data.")
}
```
### Forest Cover (2000-2024) for Project PAs
```{r graph-cover-project}
#| echo: false
#| fig-height: 3.5
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
cover_data <- project_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2000) %>%
select(WDPAID, ORIG_NAME, year, area_ha)
if (nrow(cover_data) > 0) {
p <- ggplot(cover_data, aes(
x = year,
y = area_ha,
color = ORIG_NAME,
group = ORIG_NAME,
text = paste(
"Protected Area:", ORIG_NAME,
"\nYear:", year,
"\nForest Cover:", format(round(area_ha, 0), big.mark = ","), "ha"
)
)) +
geom_line(linewidth = 1.2) +
geom_point(size = 2) +
# Add vertical line for financing agreement and shaded implementation phase
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(cover_data$area_ha, na.rm = TRUE) * 0.98,
label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(cover_data$area_ha, na.rm = TRUE) * 0.98,
label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_color_manual(values = project_colors) +
labs(
x = "Year",
y = "Forest Cover (ha)",
title = "Forest Cover (2000-2024) for Project PAs",
subtitle = "Total forest cover by protected area",
color = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2000, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom",
legend.key.width = unit(0.5, "cm"),
legend.text = element_text(size = 8),
legend.title = element_text(size = 9)
) +
guides(color = guide_legend(
nrow = 2,
override.aes = list(size = 1.5),
title.position = "top"
))
# Truncate long names in legend
p_plotly <- ggplotly(p, tooltip = "text")
for (i in 1:length(p_plotly$x$data)) {
if (!is.null(p_plotly$x$data[[i]]$name)) {
if (nchar(p_plotly$x$data[[i]]$name) > 25) {
p_plotly$x$data[[i]]$name <- paste0(substr(p_plotly$x$data[[i]]$name, 1, 22), "...")
}
}
}
p_plotly %>% layout(
legend = list(
orientation = "h",
x = 0.5,
xanchor = "center",
y = -0.15,
yanchor = "top",
font = list(size = 9)
),
margin = list(b = 100)
)
} else {
cat("**Note:** Data not available.")
}
} else {
cat("**Note:** Forest cover data not yet available. Run `process_laos_forest_cover.R` to generate this data.")
}
```
### Forest Cover (2000-2024) for All PAs
```{r graph-cover-all}
#| echo: false
#| fig-height: 3.5
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
cover_data <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2000) %>%
select(WDPAID, ORIG_NAME, year, area_ha, is_project_area)
if (nrow(cover_data) > 0) {
# Create color mapping: project areas use project_colors, others use grey
all_areas <- unique(cover_data$ORIG_NAME)
color_mapping <- ifelse(all_areas %in% names(project_colors),
project_colors[all_areas],
"grey70")
names(color_mapping) <- all_areas
p <- ggplot(cover_data, aes(
x = year,
y = area_ha,
color = ORIG_NAME,
group = ORIG_NAME,
alpha = is_project_area,
text = paste(
"Protected Area:", ORIG_NAME,
"\nYear:", year,
"\nForest Cover:", format(round(area_ha, 0), big.mark = ","), "ha",
"\nProject Area:", ifelse(is_project_area, "Yes", "No")
)
)) +
geom_line(linewidth = ifelse(cover_data$is_project_area, 1.2, 0.8)) +
geom_point(size = ifelse(cover_data$is_project_area, 2, 1.5)) +
# Add vertical line for financing agreement and shaded implementation phase
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(cover_data$area_ha, na.rm = TRUE) * 0.98,
label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(cover_data$area_ha, na.rm = TRUE) * 0.98,
label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_color_manual(values = color_mapping) +
scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0.6), guide = "none") +
labs(
x = "Year",
y = "Forest Cover (ha)",
title = "Forest Cover (2000-2024) for All PAs",
subtitle = "Project PAs shown in colors, other PAs in grey",
color = "Protected Area"
) +
scale_y_continuous(labels = comma) +
scale_x_continuous(breaks = seq(2000, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom",
legend.key.width = unit(0.5, "cm"),
legend.text = element_text(size = 8),
legend.title = element_text(size = 9)
) +
guides(color = guide_legend(
nrow = 2,
override.aes = list(size = 1.5),
title.position = "top"
))
# Truncate long names in legend and limit legend width
p_plotly <- ggplotly(p, tooltip = "text")
for (i in 1:length(p_plotly$x$data)) {
if (!is.null(p_plotly$x$data[[i]]$name)) {
if (nchar(p_plotly$x$data[[i]]$name) > 25) {
p_plotly$x$data[[i]]$name <- paste0(substr(p_plotly$x$data[[i]]$name, 1, 22), "...")
}
}
}
p_plotly %>% layout(
legend = list(
orientation = "h",
x = 0.5,
xanchor = "center",
y = -0.15,
yanchor = "top",
font = list(size = 9),
xref = "paper",
yref = "paper"
),
margin = list(b = 100)
)
} else {
cat("**Note:** Data not available.")
}
} else {
cat("**Note:** Forest cover data not yet available. Run `process_laos_forest_cover.R` to generate this data.")
}
```
## Key Findings
```{r key-findings-analysis}
#| echo: false
#| fig-height: 8
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
# Calculate detailed statistics for interpretation (all areas)
findings_data <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
group_by(WDPAID, ORIG_NAME, is_project_area) %>%
summarise(
# Overall statistics
forest_2000 = area_ha[year == 2000],
forest_2024 = area_ha[year == max(year, na.rm = TRUE)],
total_loss = forest_2000 - forest_2024,
relative_loss_pct = (total_loss / forest_2000) * 100,
# Pre-implementation and implementation period losses
loss_pre_impl = sum(loss_ha[year >= 2010 & year <= 2015], na.rm = TRUE),
loss_impl = sum(loss_ha[year >= 2018 & year <= 2023], na.rm = TRUE),
# Annual averages for key periods
avg_loss_pre_impl = mean(loss_ha[year >= 2010 & year <= 2015], na.rm = TRUE),
avg_loss_impl = mean(loss_ha[year >= 2018 & year <= 2023], na.rm = TRUE),
avg_loss_overall = mean(loss_ha[year >= 2001 & year <= 2024], na.rm = TRUE),
# Change rate
change_impl_vs_pre = ((avg_loss_impl - avg_loss_pre_impl) / avg_loss_pre_impl) * 100,
.groups = "drop"
) %>%
arrange(desc(total_loss))
# Calculate comparative statistics
project_stats <- findings_data %>% filter(is_project_area == TRUE)
project_stats_excl_hnn <- project_stats %>% filter(ORIG_NAME != "Hin Nam No")
other_stats <- findings_data %>% filter(is_project_area == FALSE)
# Averages
national_avg_relative_loss <- mean(findings_data$relative_loss_pct, na.rm = TRUE)
project_avg_relative_loss <- mean(project_stats$relative_loss_pct, na.rm = TRUE)
project_avg_relative_loss_excl_hnn <- mean(project_stats_excl_hnn$relative_loss_pct, na.rm = TRUE)
other_avg_relative_loss <- mean(other_stats$relative_loss_pct, na.rm = TRUE)
# Period-specific averages
other_avg_pre_impl <- mean(other_stats$avg_loss_pre_impl, na.rm = TRUE)
other_avg_impl <- mean(other_stats$avg_loss_impl, na.rm = TRUE)
project_avg_pre_impl <- mean(project_stats$avg_loss_pre_impl, na.rm = TRUE)
project_avg_impl <- mean(project_stats$avg_loss_impl, na.rm = TRUE)
project_avg_pre_impl_excl_hnn <- mean(project_stats_excl_hnn$avg_loss_pre_impl, na.rm = TRUE)
project_avg_impl_excl_hnn <- mean(project_stats_excl_hnn$avg_loss_impl, na.rm = TRUE)
# Change rates
project_change <- ((project_avg_impl - project_avg_pre_impl) / project_avg_pre_impl) * 100
project_change_excl_hnn <- ((project_avg_impl_excl_hnn - project_avg_pre_impl_excl_hnn) / project_avg_pre_impl_excl_hnn) * 100
other_change <- ((other_avg_impl - other_avg_pre_impl) / other_avg_pre_impl) * 100
# Display findings table
findings_table <- findings_data %>%
mutate(
`Project Area` = ifelse(is_project_area, "Yes", "No")
) %>%
select(
`Protected Area` = ORIG_NAME,
`Project Area`,
`Forest 2000 (ha)` = forest_2000,
`Forest 2024 (ha)` = forest_2024,
`Total Loss (ha)` = total_loss,
`Relative Loss (%)` = relative_loss_pct,
`Avg. Annual Loss Pre-Impl. (2010-2015) (ha)` = avg_loss_pre_impl,
`Avg. Annual Loss Impl. (2018-2023) (ha)` = avg_loss_impl,
`Change Impl. vs Pre-Impl. (%)` = change_impl_vs_pre
) %>%
arrange(`Project Area`, desc(`Total Loss (ha)`))
findings_table %>%
datatable(
options = list(pageLength = 30, dom = 'Bfrtip', scrollX = TRUE,
buttons = list(
list(extend = 'csv', text = 'Download CSV'),
list(extend = 'excel', text = 'Download Excel')
)),
extensions = 'Buttons',
caption = "Detailed Forest Cover Loss Statistics by Protected Area"
) %>%
formatRound(c('Forest 2000 (ha)', 'Forest 2024 (ha)', 'Total Loss (ha)',
'Avg. Annual Loss Pre-Impl. (2010-2015) (ha)',
'Avg. Annual Loss Impl. (2018-2023) (ha)'), 0) %>%
formatRound(c('Relative Loss (%)', 'Change Impl. vs Pre-Impl. (%)'), 1)
}
```
Based on the empirical analysis of forest cover data from 2000 to 2024 for all `r if(nrow(forest_data) > 0) { length(unique(forest_data$WDPAID)) } else { "26" }` protected areas in Laos, the following key findings emerge:
### Overall Forest Cover Development
The analysis reveals that all protected areas in Laos have experienced forest cover loss over the analysis period. The total forest cover loss across all `r if(nrow(forest_data) > 0) { length(unique(forest_data$WDPAID)) } else { "26" }` protected areas from 2000 to 2024 amounts to `r if(nrow(forest_data) > 0) { format(round(sum(forest_data$loss_ha, na.rm = TRUE), 0), big.mark = ",") } else { "data not available" }` hectares.
**Project Areas Performance (all four):**
- The four project areas account for `r if(nrow(forest_data) > 0) { project_total <- sum(project_data$loss_ha, na.rm = TRUE); all_total <- sum(forest_data$loss_ha, na.rm = TRUE); format(round((project_total / all_total) * 100, 1), big.mark = ",") } else { "N/A" }`% of total forest loss across all Lao protected areas
- Average relative loss in project areas (incl. Hin Nam No): `r if(exists("project_avg_relative_loss")) { paste0(round(project_avg_relative_loss, 2), "%") } else { "calculating..." }`
- Average relative loss in project areas (excl. Hin Nam No): `r if(exists("project_avg_relative_loss_excl_hnn")) { paste0(round(project_avg_relative_loss_excl_hnn, 2), "%") } else { "calculating..." }`
- Average relative loss in other protected areas: `r if(exists("other_avg_relative_loss")) { paste0(round(other_avg_relative_loss, 2), "%") } else { "calculating..." }`
::: {.callout-note}
## Note on Hin Nam No
Hin Nam No National Park was only supported through ancillary measures and is sometimes not listed as a formal project site. Its deforestation rate is very low compared to all other national protected areas. Where averages for project areas are reported, values are shown **both with and without Hin Nam No** to provide a more differentiated picture.
:::
### Deforestation Trend Comparison
Comparing average annual deforestation rates in the pre-implementation period (2010–2015) versus the implementation phase (2018–2023):
- **Project areas (incl. Hin Nam No):** average annual loss increased from `r if(exists("project_avg_pre_impl")) { paste0(round(project_avg_pre_impl, 0), " ha/year") } else { "calculating..." }` to `r if(exists("project_avg_impl")) { paste0(round(project_avg_impl, 0), " ha/year") } else { "calculating..." }` (`r if(exists("project_change")) { paste0(ifelse(project_change > 0, "+", ""), round(project_change, 1), "%") } else { "" }`)
- **Project areas (excl. Hin Nam No):** average annual loss increased from `r if(exists("project_avg_pre_impl_excl_hnn")) { paste0(round(project_avg_pre_impl_excl_hnn, 0), " ha/year") } else { "calculating..." }` to `r if(exists("project_avg_impl_excl_hnn")) { paste0(round(project_avg_impl_excl_hnn, 0), " ha/year") } else { "calculating..." }` (`r if(exists("project_change_excl_hnn")) { paste0(ifelse(project_change_excl_hnn > 0, "+", ""), round(project_change_excl_hnn, 1), "%") } else { "" }`)
- **Other (non-project) protected areas:** average annual loss increased from `r if(exists("other_avg_pre_impl")) { paste0(round(other_avg_pre_impl, 0), " ha/year") } else { "calculating..." }` to `r if(exists("other_avg_impl")) { paste0(round(other_avg_impl, 0), " ha/year") } else { "calculating..." }` (`r if(exists("other_change")) { paste0(ifelse(other_change > 0, "+", ""), round(other_change, 1), "%") } else { "" }`)
Annual deforestation rates increased continuously across all protected areas in Laos — including in project areas during the implementation phase. However, the increase in deforestation rates was **markedly steeper in non-project protected areas**.
### Annual Deforestation Rate Trends
The following chart shows the average annual deforestation rate (forest loss as a percentage of the previous year's forest cover) for three groups: all other national protected areas, project areas excluding Hin Nam No, and all project areas including Hin Nam No. Linear trend lines illustrate the long-term trajectory for each group.
```{r graph-deforestation-rate-trends}
#| echo: false
#| fig-height: 5
#| fig-width: 10
#| out-width: "100%"
#| fig-align: center
if (nrow(forest_data) > 0) {
# Calculate annual deforestation rate as % of previous year's forest cover
deforestation_rate <- forest_data %>%
mutate(WDPAID = as.numeric(WDPAID)) %>%
filter(year >= 2001) %>%
group_by(WDPAID, ORIG_NAME, is_project_area) %>%
arrange(year) %>%
mutate(
prev_area = lag(area_ha),
deforestation_rate_pct = (loss_ha / (area_ha + loss_ha)) * 100
) %>%
filter(!is.na(deforestation_rate_pct)) %>%
ungroup()
# Define which PAs are project areas excl. Hin Nam No
deforestation_rate <- deforestation_rate %>%
mutate(
group = case_when(
is_project_area & ORIG_NAME != "Hin Nam No" ~ "Project PAs (excl. Hin Nam No)",
is_project_area & ORIG_NAME == "Hin Nam No" ~ "Hin Nam No only",
TRUE ~ "Other National PAs"
)
)
# Calculate group averages for the three lines:
# 1. Other national PAs
other_avg <- deforestation_rate %>%
filter(group == "Other National PAs") %>%
group_by(year) %>%
summarise(avg_rate = mean(deforestation_rate_pct, na.rm = TRUE), .groups = "drop") %>%
mutate(group = "Other National PAs")
# 2. Project PAs excl. Hin Nam No
project_excl_hnn_avg <- deforestation_rate %>%
filter(group == "Project PAs (excl. Hin Nam No)") %>%
group_by(year) %>%
summarise(avg_rate = mean(deforestation_rate_pct, na.rm = TRUE), .groups = "drop") %>%
mutate(group = "Project PAs (excl. Hin Nam No)")
# 3. All project PAs (incl. Hin Nam No)
project_all_avg <- deforestation_rate %>%
filter(is_project_area) %>%
group_by(year) %>%
summarise(avg_rate = mean(deforestation_rate_pct, na.rm = TRUE), .groups = "drop") %>%
mutate(group = "All Project PAs")
# Combine
trend_data <- bind_rows(other_avg, project_excl_hnn_avg, project_all_avg) %>%
mutate(group = factor(group, levels = c("Other National PAs",
"Project PAs (excl. Hin Nam No)",
"All Project PAs")))
trend_colors <- c(
"Other National PAs" = "#d62728",
"Project PAs (excl. Hin Nam No)" = "#9467bd",
"All Project PAs" = "#2ca02c"
)
trend_linetypes <- c(
"Other National PAs" = "dotted",
"Project PAs (excl. Hin Nam No)" = "longdash",
"All Project PAs" = "dashed"
)
# Compute linear trend lines manually for reliable plotly rendering
trend_lines <- trend_data %>%
group_by(group) %>%
do({
mod <- lm(avg_rate ~ year, data = .)
data.frame(
year = .$year,
fitted = predict(mod, newdata = .),
group = .$group
)
}) %>%
ungroup()
# Lighter colors for trend lines
trend_line_colors <- c(
"Other National PAs" = "#f4a0a0",
"Project PAs (excl. Hin Nam No)" = "#c9b3df",
"All Project PAs" = "#a3d9a3"
)
p <- ggplot() +
# Linear trend lines (drawn first, behind data lines)
geom_line(data = trend_lines, aes(x = year, y = fitted, color = group),
linetype = "dotted", linewidth = 0.7, alpha = 0.5, show.legend = FALSE) +
# Observational data lines (dashed patterns)
geom_line(data = trend_data, aes(
x = year,
y = avg_rate,
color = group,
linetype = group
), linewidth = 1.0) +
# Data point markers along the lines
geom_point(data = trend_data, aes(
x = year,
y = avg_rate,
color = group,
text = paste(
"Group:", group,
"\nYear:", year,
"\nAvg. Deforestation Rate:", round(avg_rate, 3), "%"
)
), size = 2, show.legend = FALSE) +
geom_vline(xintercept = 2015, linetype = "dashed", color = "darkgrey", linewidth = 0.5) +
annotate("rect", xmin = 2017.5, xmax = 2023.5, ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "steelblue") +
annotate("text", x = 2015, y = max(trend_data$avg_rate, na.rm = TRUE) * 0.98,
label = "Financing\nAgreement", color = "darkgrey", size = 2.5, fontface = "italic", hjust = 1.1) +
annotate("text", x = 2020.5, y = max(trend_data$avg_rate, na.rm = TRUE) * 0.98,
label = "Implementation Phase\n(2018-2023)", color = "steelblue", size = 2.5, fontface = "italic", alpha = 0.7) +
scale_color_manual(values = trend_colors) +
scale_linetype_manual(values = trend_linetypes) +
labs(
x = "Year",
y = "Avg. Annual Deforestation Rate (%)",
title = "Average Annual Deforestation Rate by Group (2001-2024)",
subtitle = "Forest loss as % of previous year's forest cover, with linear trend lines",
color = "Group",
linetype = "Group"
) +
scale_y_continuous(labels = function(x) paste0(format(x, nsmall = 2), "%")) +
scale_x_continuous(breaks = seq(2001, 2024, by = 2)) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12, color = "gray50"),
legend.position = "bottom",
legend.key.width = unit(1.5, "cm"),
legend.text = element_text(size = 9),
legend.title = element_text(size = 10)
) +
guides(color = guide_legend(nrow = 2, title.position = "top"),
linetype = guide_legend(nrow = 2, title.position = "top"))
p_plotly <- ggplotly(p, tooltip = "text")
# Post-process traces: lighten trend lines, hide duplicate legend entries
trend_line_colors_map <- c(
"Other National PAs" = "rgba(214, 39, 40, 0.35)",
"Project PAs (excl. Hin Nam No)" = "rgba(148, 103, 189, 0.35)",
"All Project PAs" = "rgba(44, 160, 44, 0.35)"
)
for (i in seq_along(p_plotly$x$data)) {
trace <- p_plotly$x$data[[i]]
trace_name <- trace$name
if (!is.null(trace_name) && grepl("^\\(", trace_name)) {
# Trend lines and marker traces have names like "(Group,1)"
# Extract the group name to set lighter colors for trend lines
group_name <- gsub("^\\((.+),1\\)$", "\\1", trace_name)
if (trace$mode == "lines" && !is.null(trend_line_colors_map[group_name])) {
# This is a trend line — lighten its color
p_plotly$x$data[[i]]$line$color <- trend_line_colors_map[group_name]
}
p_plotly$x$data[[i]]$showlegend <- FALSE
}
}
p_plotly <- p_plotly %>%
layout(
legend = list(
orientation = "h",
x = 0.5,
xanchor = "center",
y = -0.2,
yanchor = "top",
font = list(size = 10)
),
margin = list(b = 120)
)
p_plotly
}
```
The chart shows that while deforestation rates increase across all groups, the trend is substantially steeper for non-project protected areas, suggesting that project interventions may have had a moderating effect on deforestation.
### Area-Specific Patterns
1. **Most Affected Areas**: Among all Lao protected areas, the areas with the highest total forest cover loss are `r if(nrow(findings_data) > 0) { top_areas <- findings_data %>% arrange(desc(total_loss)) %>% head(3); paste(top_areas$ORIG_NAME, collapse = ", ") } else { "to be determined from data" }`. `r if(nrow(findings_data) > 0) { top_areas <- findings_data %>% arrange(desc(total_loss)) %>% head(3); project_in_top <- sum(top_areas$is_project_area); if(project_in_top > 0) { paste0(project_in_top, " of these are project areas.") } else { "None of these are project areas." } } else { "" }`
2. **Relative Loss Comparison**: When considering relative loss (percentage of original forest cover), project areas show `r if(exists("project_avg_relative_loss") && exists("other_avg_relative_loss")) { if(project_avg_relative_loss < other_avg_relative_loss) { paste0("lower relative loss (", round(project_avg_relative_loss, 2), "% vs ", round(other_avg_relative_loss, 2), "% in other PAs)") } else { paste0("higher relative loss (", round(project_avg_relative_loss, 2), "% vs ", round(other_avg_relative_loss, 2), "% in other PAs)") } } else { "varying patterns" }`.
3. **Hin Nam No**: Hin Nam No stands out with a very low deforestation rate compared to all other protected areas. When excluded from project area averages, the remaining three project areas show `r if(exists("project_avg_relative_loss_excl_hnn") && exists("project_avg_relative_loss")) { if(project_avg_relative_loss_excl_hnn > project_avg_relative_loss) { "higher" } else { "lower" } } else { "different" }` relative loss rates.
### Conclusions
The project **could not stop or reverse the overall deforestation trend** — annual deforestation rates continued to increase in project areas, including during the implementation phase (2018–2023). However, the **comparison with non-project protected areas suggests that the project interventions were partially effective**: the increase in deforestation rates was substantially lower in project areas than in comparable non-project areas. This indicates that while absolute deforestation targets were not achieved, the project may have contributed to **avoiding even greater forest loss**.
## Data Sources & References
- **WDPA**: UNEP-WCMC and IUCN (2024), Protected Planet: The World Database on Protected Areas (WDPA)
- **GFW**: Hansen, M. C., et al. (2013). High-Resolution Global Maps of 21st-Century Forest Cover Change. *Science*, 342(15 November): 850–53. Data version: GFC-2024-v1.12.
- **mapme.biodiversity**: R package for automated biodiversity indicator calculation
---
*Report generated: `r Sys.Date()`*