---
title: "Forest Cover Assessment for Project Villages in Tripura, India"
subtitle: "Analysis of Forest Cover Loss Dynamics (2001-2024)"
author:
- name: "Johannes Schielein"
affiliations: "MAPME Initiative"
date: "2026-03-07"
date-modified: last-modified
abstract: |
This report assesses forest cover loss dynamics in 83 project villages and 838 control
villages across Tripura state, India, using Global Forest Watch data (2001-2024). The
analysis examines whether sufficient forest remains in project villages for planned
conservation zoning and land-use planning activities, and compares loss trends across
three distinct periods.
categories: [India, Tripura, Villages, Forest Loss, Conservation, Land-Use Zoning]
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(scales)
library(htmltools)
library(RColorBrewer)
library(plotly)
library(DT)
options(scipen = 10000)
# ----- Load Forest Cover Data -----
if (file.exists("../data/india/output/india_forest_cover_2001_2024.csv")) {
gfw_stats <- read_csv(
"../data/india/output/india_forest_cover_2001_2024.csv",
show_col_types = FALSE
)
} else {
warning("Forest cover data not found. Run scripts/india/process_forest_cover.R first!")
gfw_stats <- tibble()
}
# ----- Load Village Geometries -----
all_villages <- st_read(
"../data/india/input/tripura_villages_esri.geojson",
quiet = TRUE
) %>% st_make_valid()
vdpic <- st_read(
"../data/india/input/tripura_vdpic_matched.geojson",
quiet = TRUE
) %>% st_make_valid()
# ----- Prepare Project Village Summary -----
if (nrow(gfw_stats) > 0) {
# Baseline year data
baseline <- gfw_stats %>% filter(year == 2001)
latest <- gfw_stats %>% filter(year == 2024)
# Project villages data
proj_stats <- gfw_stats %>% filter(is_project_village == TRUE)
ctrl_stats <- gfw_stats %>% filter(is_project_village == FALSE)
# Join confidence info to village geometries for mapping
vdpic_matched <- vdpic %>% filter(confidence != "unmatched" | is.na(confidence))
# Compute per-village summary
village_summary <- gfw_stats %>%
group_by(objectid, name, csv_name, csv_district, csv_block,
confidence, match_method, is_project_village) %>%
summarise(
forest_2001 = area_ha[year == 2001],
forest_2024 = area_ha[year == 2024],
abs_loss = area_ha[year == 2001] - area_ha[year == 2024],
rel_loss = (area_ha[year == 2001] - area_ha[year == 2024]) / area_ha[year == 2001] * 100,
.groups = "drop"
)
# Period-level stats per village
period_stats <- gfw_stats %>%
filter(year %in% c(2001, 2010, 2011, 2020, 2021, 2024)) %>%
mutate(
period = case_when(
year %in% c(2001, 2010) ~ "2001-2010",
year %in% c(2011, 2020) ~ "2011-2020",
year %in% c(2021, 2024) ~ "2021-2024"
),
period_pos = case_when(
year %in% c(2001, 2011, 2021) ~ "start",
TRUE ~ "end"
)
) %>%
select(objectid, name, csv_name, is_project_village, period, period_pos, area_ha, baseline_2001) %>%
pivot_wider(names_from = period_pos, values_from = area_ha) %>%
mutate(
period_loss_ha = start - end,
period_loss_rel = period_loss_ha / baseline_2001 * 100
)
}
```
## Introduction
This report examines forest cover loss dynamics in villages across Tripura, a small state in northeast India. The analysis serves a specific conservation question: a rural development project with conservation components supports 83 villages in the Dhalai and North Tripura districts. The project includes agricultural support, land-use zoning, and the designation of conservation set-aside areas.
Given the substantial forest cover loss observed across Tripura in recent years, this analysis assesses:
1. **How much forest remains** in the project villages for conservation zoning
2. **How loss has accelerated or decelerated** across three periods (2001--2010, 2011--2020, 2021--2024)
3. **How project villages compare** to the broader trend across all Tripura villages
::: {.callout-note}
## Key Finding
Project villages lost approximately **23% of their 2001 forest cover** by 2024, with loss rates accelerating sharply after 2010. The most recent period (2021--2024) shows the highest annual loss rate, raising urgent questions about the feasibility of conservation zoning in heavily affected villages.
:::
## Data and Methods
### Village Boundary Data
The analysis covers **917 village polygons** from the ESRI India Living Atlas (Census 2021 boundaries) across all 8 districts of Tripura, plus 4 additional villages matched via geoBoundaries (Census 2011 boundaries). Of these, **83 villages** were identified as project villages through a multi-step name-matching process.
The matching linked field programme records (VDPIC village names) to official Census boundary polygons. Because village names in northeast India often have multiple spellings --- reflecting Bengali, tribal, and English transliterations --- the matching used a combination of:
- **Exact name matching** after normalizing spelling (removing punctuation, standardizing whitespace)
- **Fuzzy matching** using string similarity algorithms, with manual review of ambiguous cases
- **Phonetic normalization** specific to Tripura (e.g., the suffix "-cherra"/"-chhara" refers to the same geographic feature --- a stream)
- **Manual overrides** for seven villages where automated matching failed due to significant name divergence (e.g., "Narikel Kunja" = "Kunjaban", both meaning "coconut grove" in Bengali)
Each match was assigned a confidence level: **confirmed** (exact match or high-similarity score), **likely** (moderate similarity), or **uncertain** (lower similarity or parent-village assignment). Nine villages from the field programme could not be matched to any Census polygon --- these are sub-hamlets, reserve forest sections, or settlements not recorded as separate Census units.
### Overview of Project Villages
```{r table-villages}
#| echo: false
if (nrow(gfw_stats) > 0) {
vdpic_table <- vdpic %>%
st_drop_geometry() %>%
filter(confidence != "unmatched") %>%
select(
District = csv_district,
Block = csv_block,
`Village Name` = csv_name,
`Match Method` = match_method,
`Match Score` = match_score,
Confidence = confidence,
`Matched Name` = src_name,
`Census Name` = src_censusname
) %>%
arrange(District, Block, `Village Name`)
vdpic_table %>%
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 = "Project villages matched to Census boundary polygons (90 of 99)"
)
}
```
**Column descriptions:**
- **District / Block**: Administrative location from the field programme records
- **Village Name**: Name used in the field programme (VDPIC records)
- **Match Method**: How the match was established (exact, fuzzy_high, fuzzy_synonym, manual_override, or tripura_phonetic_fuzzy)
- **Match Score**: String similarity score (0--100); blank for exact matches and manual overrides
- **Confidence**: Overall match quality --- confirmed (green on map), likely (blue), or uncertain (orange)
- **Matched Name / Census Name**: Official names from the matched boundary dataset
### Interactive Map of Study Area
```{r map-villages}
#| echo: false
#| fig-height: 6
#| out-width: "100%"
#| classes: map-container
# Prepare confidence-colored vdpic polygons
vdpic_confirmed <- vdpic %>% filter(confidence == "confirmed")
vdpic_likely <- vdpic %>% filter(confidence == "likely")
vdpic_uncertain <- vdpic %>% filter(confidence == "uncertain" & !st_is_empty(geometry))
# Non-project villages (all_villages minus matched ones)
# Get objectids of matched project villages
proj_objectids <- gfw_stats %>%
filter(is_project_village == TRUE & year == 2001) %>%
pull(objectid) %>%
unique()
non_project <- all_villages %>% filter(!objectid %in% proj_objectids)
project_esri <- all_villages %>% filter(objectid %in% proj_objectids)
map_attribution <- "Forest Loss: Hansen et al. (2013) | Village Boundaries: ESRI India Living Atlas"
map1 <- leaflet() %>%
addTiles(group = "OpenStreetMap") %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
addProviderTiles(providers$Esri.WorldShadedRelief, group = "Topography") %>%
# Non-project village polygons
addPolygons(
data = non_project,
fillColor = "grey",
fillOpacity = 0.2,
color = "darkgrey",
weight = 1,
opacity = 0.6,
group = "Non-Project Villages",
label = ~name
) %>%
# Confirmed project villages
addPolygons(
data = vdpic_confirmed,
fillColor = "grey",
fillOpacity = 0.2,
color = "green",
weight = 2.5,
opacity = 1,
group = "Project Villages (Confirmed)",
label = ~csv_name,
popup = ~paste0(
"<b>", csv_name, "</b><br>",
"Block: ", csv_block, "<br>",
"Confidence: ", confidence, "<br>",
"Match: ", match_method
)
) %>%
# Likely project villages
addPolygons(
data = vdpic_likely,
fillColor = "grey",
fillOpacity = 0.2,
color = "blue",
weight = 2.5,
opacity = 1,
group = "Project Villages (Likely)",
label = ~csv_name,
popup = ~paste0(
"<b>", csv_name, "</b><br>",
"Block: ", csv_block, "<br>",
"Confidence: ", confidence, "<br>",
"Match: ", match_method
)
) %>%
# Uncertain project villages
addPolygons(
data = vdpic_uncertain,
fillColor = "grey",
fillOpacity = 0.2,
color = "orange",
weight = 2.5,
opacity = 1,
group = "Project Villages (Uncertain)",
label = ~csv_name,
popup = ~paste0(
"<b>", csv_name, "</b><br>",
"Block: ", csv_block, "<br>",
"Confidence: ", confidence, "<br>",
"Match: ", match_method
)
) %>%
addFullscreenControl() %>%
addLegend(
"bottomright",
colors = c("green", "blue", "orange", "darkgrey"),
labels = c("Project (Confirmed)", "Project (Likely)",
"Project (Uncertain)", "Non-Project"),
title = "Village Classification",
opacity = 1
) %>%
addLayersControl(
baseGroups = c("CartoDB", "OpenStreetMap", "Satellite", "Topography"),
overlayGroups = c("Project Villages (Confirmed)", "Project Villages (Likely)",
"Project Villages (Uncertain)", "Non-Project Villages"),
options = layersControlOptions(collapsed = FALSE)
)
map1
```
### Forest Cover Data
Forest cover statistics were computed using the [mapme.biodiversity](https://github.com/mapme-initiative/mapme.biodiversity) R package with data from the **Global Forest Watch** dataset (Hansen et al., 2013), version GFC-2024-v1.12. The analysis covers the period **2001--2024**.
**Forest definition used:**
- **Canopy cover threshold**: 10% (pixels with at least 10% tree canopy density in the year 2000 are classified as forest)
- **Minimum patch size**: 1 hectare
- **Resolution**: 30 meters (Landsat-derived)
Forest cover loss is defined as a *stand-replacement disturbance* --- a complete removal of tree cover canopy at the 30m pixel level. The GFW dataset detects loss events but does not capture forest degradation (partial canopy reduction) or regrowth.
### Tripura's Forest Ecosystem
Tripura's forests are classified primarily as **tropical semi-evergreen** and **moist deciduous** types, with significant areas of bamboo brakes --- a sub-climax formation resulting from historic shifting cultivation (jhum). The dominant tree species include *Dipterocarpus turbinatus*, *Artocarpus chama*, and *Castanopsis indica*. According to the India State of Forest Report 2023, Tripura maintains approximately **74.7% forest cover** relative to its geographic area, though this has been declining --- the state lost roughly 95 km^2^ of forest cover between the 2021 and 2023 assessment periods.
### Accuracy Considerations
The GFW dataset provides globally consistent forest monitoring but has known limitations relevant to this context:
- **Shifting cultivation**: The dataset cannot distinguish between permanent deforestation and temporary clearing for jhum cultivation followed by regrowth. In areas with active shifting cultivation, GFW may record "loss" events that are part of a rotational cycle rather than permanent conversion.
- **Bamboo and plantation forests**: GFW treats all tree cover uniformly and cannot distinguish natural forests from bamboo stands or plantations. Bamboo regeneration after clearing may not be captured as "regrowth" if canopy height remains below the detection threshold.
- **Small-scale disturbances**: At 30m resolution, selective logging and small clearings below the pixel size may go undetected, potentially underestimating degradation in smallholder landscapes.
- **Cloud cover**: Northeast India's monsoon climate means persistent cloud cover during parts of the year, which can affect satellite image availability and detection accuracy.
These limitations should be considered when interpreting the results. The GFW data is best suited for detecting broad-scale loss trends rather than precise village-level accounting.
## Forest Cover Loss in Project Villages
```{r period-data}
#| include: false
if (nrow(gfw_stats) > 0) {
# Compute annual relative loss trajectory per village (relative to 2001 baseline)
proj_trajectory <- proj_stats %>%
group_by(objectid) %>%
mutate(
rel_forest_remaining = area_ha / baseline_2001 * 100,
cum_rel_loss = (baseline_2001 - area_ha) / baseline_2001 * 100
) %>%
ungroup()
# Aggregate: median and IQR across project villages per year
proj_annual_agg <- proj_trajectory %>%
group_by(year) %>%
summarise(
median_remaining = median(rel_forest_remaining, na.rm = TRUE),
q25_remaining = quantile(rel_forest_remaining, 0.25, na.rm = TRUE),
q75_remaining = quantile(rel_forest_remaining, 0.75, na.rm = TRUE),
mean_remaining = mean(rel_forest_remaining, na.rm = TRUE),
median_cum_loss = median(cum_rel_loss, na.rm = TRUE),
q25_cum_loss = quantile(cum_rel_loss, 0.25, na.rm = TRUE),
q75_cum_loss = quantile(cum_rel_loss, 0.75, na.rm = TRUE),
.groups = "drop"
)
# Period-level data for project villages
proj_periods <- period_stats %>%
filter(is_project_village == TRUE) %>%
mutate(
# Annualized relative loss (% of 2001 baseline per year)
period_years = case_when(
period == "2001-2010" ~ 9,
period == "2011-2020" ~ 9,
period == "2021-2024" ~ 3
),
annual_rel_loss = period_loss_rel / period_years
)
# Same for all villages (for comparison)
all_periods <- period_stats %>%
mutate(
period_years = case_when(
period == "2001-2010" ~ 9,
period == "2011-2020" ~ 9,
period == "2021-2024" ~ 3
),
annual_rel_loss = period_loss_rel / period_years,
village_type = ifelse(is_project_village, "Project Villages", "Non-Project Villages")
)
}
```
### Cumulative Forest Cover Trajectory
The figure below shows how forest cover in project villages has evolved since 2001. Each village's forest area is expressed as a percentage of its 2001 baseline. The shaded ribbon represents the interquartile range (25th--75th percentile) across all 83 project villages, capturing the spread of outcomes.
```{r fig-trajectory}
#| echo: false
#| fig-height: 5
if (nrow(gfw_stats) > 0) {
p_traj <- ggplot(proj_annual_agg, aes(x = year)) +
geom_ribbon(aes(ymin = q25_remaining, ymax = q75_remaining),
fill = "#2166ac", alpha = 0.2) +
geom_line(aes(y = median_remaining), color = "#2166ac", linewidth = 1.2) +
geom_line(aes(y = mean_remaining), color = "#2166ac", linewidth = 0.6, linetype = "dashed") +
# Period boundaries
geom_vline(xintercept = c(2010.5, 2020.5), linetype = "dotted", color = "grey40", linewidth = 0.4) +
annotate("text", x = 2005.5, y = 72, label = "Period 1\n2001-2010",
size = 3, color = "grey40", hjust = 0.5) +
annotate("text", x = 2015.5, y = 72, label = "Period 2\n2011-2020",
size = 3, color = "grey40", hjust = 0.5) +
annotate("text", x = 2022.5, y = 72, label = "Period 3\n2021-2024",
size = 3, color = "grey40", hjust = 0.5) +
scale_y_continuous(labels = function(x) paste0(x, "%"), limits = c(70, 105)) +
scale_x_continuous(breaks = seq(2001, 2024, 2)) +
labs(
title = "Forest Cover Trajectory in Project Villages (2001-2024)",
subtitle = "Median (solid), mean (dashed), interquartile range (shaded) across 83 villages",
x = NULL,
y = "Forest remaining (% of 2001 baseline)"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(color = "grey40", size = 11),
panel.grid.minor = element_blank()
)
ggplotly(p_traj, tooltip = c("x", "y")) %>%
layout(
hovermode = "x unified",
margin = list(t = 80)
)
}
```
### Period Comparison: Loss Acceleration
To assess whether forest loss is accelerating, we compare three periods using **annualized relative loss** --- the percentage of each village's 2001 forest baseline lost per year within each period. This metric normalizes for the different period lengths (9, 9, and 3 years) and makes trends directly comparable.
```{r fig-boxplots}
#| echo: false
#| fig-height: 5.5
if (nrow(gfw_stats) > 0) {
# Order periods chronologically
proj_periods$period <- factor(proj_periods$period,
levels = c("2001-2010", "2011-2020", "2021-2024"))
# Period summary for annotation
period_medians <- proj_periods %>%
group_by(period) %>%
summarise(
median_annual = median(annual_rel_loss, na.rm = TRUE),
mean_annual = mean(annual_rel_loss, na.rm = TRUE),
total_loss = median(period_loss_rel, na.rm = TRUE),
.groups = "drop"
)
# Colors for periods
period_colors <- c("2001-2010" = "#4575b4", "2011-2020" = "#fc8d59", "2021-2024" = "#d73027")
p_box <- ggplot(proj_periods, aes(x = period, y = annual_rel_loss, fill = period)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 2, outlier.alpha = 0.5) +
geom_jitter(width = 0.15, alpha = 0.3, size = 1.2, color = "grey30") +
scale_fill_manual(values = period_colors) +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
labs(
title = "Annualized Forest Cover Loss by Period",
subtitle = "Each point represents one project village (83 villages)",
x = NULL,
y = "Annual loss (% of 2001 baseline per year)"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(color = "grey40", size = 11),
legend.position = "none",
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()
)
ggplotly(p_box) %>%
layout(margin = list(t = 80))
}
```
```{r period-summary-table}
#| echo: false
if (nrow(gfw_stats) > 0) {
period_summary <- proj_periods %>%
group_by(Period = period) %>%
summarise(
`Years` = first(period_years),
`Median Total Loss (% of 2001)` = round(median(period_loss_rel, na.rm = TRUE), 1),
`Mean Total Loss (% of 2001)` = round(mean(period_loss_rel, na.rm = TRUE), 1),
`Median Annual Loss (%/yr)` = round(median(annual_rel_loss, na.rm = TRUE), 2),
`Mean Annual Loss (%/yr)` = round(mean(annual_rel_loss, na.rm = TRUE), 2),
`Max Total Loss (%)` = round(max(period_loss_rel, na.rm = TRUE), 1),
.groups = "drop"
)
period_summary %>%
datatable(
options = list(dom = 't', pageLength = 5, scrollX = TRUE),
caption = "Period-level forest cover loss summary for 83 project villages",
rownames = FALSE
)
}
```
```{r interpretation-periods}
#| echo: false
#| results: asis
if (nrow(gfw_stats) > 0) {
pm <- period_medians
cat(paste0(
"\n**Interpretation:** The data reveals a clear acceleration of forest loss across the three periods. ",
"During 2001--2010, the median project village lost approximately **",
round(pm$median_annual[pm$period == "2001-2010"], 2),
"% of its 2001 forest per year**. ",
"This rate increased substantially in 2011--2020 to **",
round(pm$median_annual[pm$period == "2011-2020"], 2),
"%/year**, and accelerated further in 2021--2024 to **",
round(pm$median_annual[pm$period == "2021-2024"], 2),
"%/year**. ",
"The spread of outcomes (visible in the boxplot) also widens over time, indicating that ",
"some villages are experiencing dramatically faster loss than others. ",
"The most affected village lost ", round(max(proj_periods$period_loss_rel[proj_periods$period == "2021-2024"], na.rm = TRUE), 1),
"% of its 2001 baseline in just the 2021--2024 period alone.\n"
))
}
```
### Distribution of Total Forest Loss (2001--2024)
The histogram below shows how total relative forest loss (2001--2024) is distributed across project villages.
```{r fig-histogram}
#| echo: false
#| fig-height: 4.5
if (nrow(gfw_stats) > 0) {
proj_summary <- village_summary %>% filter(is_project_village == TRUE)
p_hist <- ggplot(proj_summary, aes(x = rel_loss)) +
geom_histogram(binwidth = 5, fill = "#2166ac", color = "white", alpha = 0.8) +
geom_vline(aes(xintercept = median(rel_loss, na.rm = TRUE)),
color = "#d73027", linewidth = 1, linetype = "dashed") +
annotate("text", x = median(proj_summary$rel_loss, na.rm = TRUE) + 1.5,
y = Inf, label = paste0("Median: ", round(median(proj_summary$rel_loss, na.rm = TRUE), 1), "%"),
hjust = 0, vjust = 2, color = "#d73027", fontface = "bold", size = 4) +
scale_x_continuous(labels = function(x) paste0(x, "%")) +
labs(
title = "Distribution of Total Forest Cover Loss in Project Villages",
subtitle = "Relative loss as percentage of 2001 forest baseline (2001-2024)",
x = "Total relative forest cover loss (%)",
y = "Number of villages"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(color = "grey40", size = 11),
panel.grid.minor = element_blank()
)
ggplotly(p_hist) %>%
layout(margin = list(t = 80))
}
```
```{r interpretation-histogram}
#| echo: false
#| results: asis
if (nrow(gfw_stats) > 0) {
ps <- proj_summary
cat(paste0(
"\nThe median project village lost **", round(median(ps$rel_loss, na.rm = TRUE), 1),
"%** of its 2001 forest cover by 2024. ",
round(sum(ps$rel_loss > 30, na.rm = TRUE) / nrow(ps) * 100, 0),
"% of villages lost more than 30% of their forest, while ",
round(sum(ps$rel_loss < 10, na.rm = TRUE) / nrow(ps) * 100, 0),
"% of villages lost less than 10%. ",
"This wide spread suggests that forest loss pressures vary considerably across the project area, ",
"and that conservation zoning may need to be prioritized differently depending on local conditions.\n"
))
}
```
### Comparison: Project Villages vs. All Tripura Villages
```{r fig-comparison}
#| echo: false
#| fig-height: 5.5
if (nrow(gfw_stats) > 0) {
all_periods$period <- factor(all_periods$period,
levels = c("2001-2010", "2011-2020", "2021-2024"))
p_comp <- ggplot(all_periods, aes(x = period, y = annual_rel_loss, fill = village_type)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 1.5,
outlier.alpha = 0.3, position = position_dodge(0.8)) +
scale_fill_manual(values = c("Project Villages" = "#2166ac",
"Non-Project Villages" = "#b2182b")) +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
labs(
title = "Forest Loss Comparison: Project vs. Non-Project Villages",
subtitle = "Annualized loss as % of 2001 baseline, by period",
x = NULL,
y = "Annual loss (% of 2001 baseline per year)",
fill = NULL
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(color = "grey40", size = 11),
legend.position = "bottom",
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()
)
ggplotly(p_comp) %>%
layout(
margin = list(t = 80),
legend = list(orientation = "h", x = 0.2, y = -0.15)
)
}
```
```{r interpretation-comparison}
#| echo: false
#| results: asis
if (nrow(gfw_stats) > 0) {
proj_med <- all_periods %>%
filter(village_type == "Project Villages") %>%
group_by(period) %>%
summarise(med = median(annual_rel_loss, na.rm = TRUE), .groups = "drop")
ctrl_med <- all_periods %>%
filter(village_type == "Non-Project Villages") %>%
group_by(period) %>%
summarise(med = median(annual_rel_loss, na.rm = TRUE), .groups = "drop")
cat(paste0(
"\n**Interpretation:** The comparison shows that project villages and non-project villages across ",
"Tripura follow broadly similar loss trends --- both groups show accelerating loss over time. ",
"In the most recent period (2021--2024), project villages show a median annual loss of **",
round(proj_med$med[proj_med$period == "2021-2024"], 2), "%/year** compared to **",
round(ctrl_med$med[ctrl_med$period == "2021-2024"], 2),
"%/year** for non-project villages. ",
"This suggests that the forest loss pressures are driven by regional factors (road expansion, ",
"population growth, agricultural conversion) rather than village-specific dynamics.\n"
))
}
```
### Top 15 Most Affected Project Villages
```{r table-top-loss}
#| echo: false
if (nrow(gfw_stats) > 0) {
top_loss <- proj_summary %>%
arrange(desc(rel_loss)) %>%
head(15) %>%
mutate(
forest_2001 = round(forest_2001, 0),
forest_2024 = round(forest_2024, 0),
abs_loss = round(abs_loss, 0),
rel_loss = round(rel_loss, 1)
) %>%
select(
`Village` = csv_name,
`Block` = csv_block,
`Forest 2001 (ha)` = forest_2001,
`Forest 2024 (ha)` = forest_2024,
`Loss (ha)` = abs_loss,
`Loss (%)` = rel_loss
)
top_loss %>%
datatable(
options = list(dom = 't', pageLength = 15, scrollX = TRUE),
caption = "Top 15 project villages by relative forest cover loss (2001-2024)",
rownames = FALSE
) %>%
formatStyle(
'Loss (%)',
background = styleColorBar(range(top_loss$`Loss (%)`), '#fc8d59'),
backgroundSize = '98% 88%',
backgroundRepeat = 'no-repeat',
backgroundPosition = 'center'
)
}
```
### Forest Cover Loss Map
The map below shows forest cover loss across all villages. Point symbology encodes both **absolute loss** (circle size) and **relative loss** (circle color). Project villages are shown at full opacity; non-project villages are semi-transparent to allow the project villages to stand out. Toggle the GFW tile layers to see the spatial pattern of loss across three periods.
```{r map-forest-loss}
#| echo: false
#| fig-height: 7
#| out-width: "100%"
#| classes: map-container
if (nrow(gfw_stats) > 0) {
# Prepare loss data with centroids
village_centroids <- all_villages %>%
st_centroid() %>%
st_transform(4326)
village_map_data <- village_centroids %>%
st_drop_geometry() %>%
left_join(village_summary, by = "objectid") %>%
bind_cols(
st_coordinates(village_centroids) %>% as_tibble() %>% rename(lng = X, lat = Y)
) %>%
filter(!is.na(rel_loss))
# Also add gb_direct villages if they exist
gb_direct <- village_summary %>% filter(is.na(objectid) | !objectid %in% all_villages$objectid)
# Categorize losses
village_map_data$rel_loss_cat <- cut(
village_map_data$rel_loss,
breaks = c(-Inf, 5, 10, 20, 30, Inf),
labels = c("0-5%", "5-10%", "10-20%", "20-30%", ">30%")
)
village_map_data$abs_loss_cat <- cut(
village_map_data$abs_loss,
breaks = c(-Inf, 10, 50, 200, 1000, Inf),
labels = c("<10 ha", "10-50 ha", "50-200 ha", "200-1000 ha", ">1000 ha")
)
# Separate project and non-project for different opacity
map_proj <- village_map_data %>% filter(is_project_village == TRUE)
map_ctrl <- village_map_data %>% filter(is_project_village == FALSE)
# Create labels
map_proj$popup <- paste0(
"<b>", map_proj$csv_name, "</b> (Project Village)<br>",
"Block: ", map_proj$csv_block, "<br>",
"Forest 2001: ", format(round(map_proj$forest_2001, 0), big.mark = ","), " ha<br>",
"Forest 2024: ", format(round(map_proj$forest_2024, 0), big.mark = ","), " ha<br>",
"Loss: ", format(round(map_proj$abs_loss, 0), big.mark = ","), " ha (",
round(map_proj$rel_loss, 1), "%)"
)
map_ctrl$popup <- paste0(
"<b>", map_ctrl$name.x, "</b><br>",
"District: ", map_ctrl$district, "<br>",
"Forest 2001: ", format(round(map_ctrl$forest_2001, 0), big.mark = ","), " ha<br>",
"Forest 2024: ", format(round(map_ctrl$forest_2024, 0), big.mark = ","), " ha<br>",
"Loss: ", format(round(map_ctrl$abs_loss, 0), big.mark = ","), " ha (",
round(map_ctrl$rel_loss, 1), "%)"
)
# Color palette
pal_loss <- colorFactor(
palette = rev(brewer.pal(5, "RdYlGn")),
domain = village_map_data$rel_loss_cat
)
# Radius function
get_radius <- function(abs_cat) {
case_when(
abs_cat == ">1000 ha" ~ 10,
abs_cat == "200-1000 ha" ~ 7,
abs_cat == "50-200 ha" ~ 5,
abs_cat == "10-50 ha" ~ 3.5,
TRUE ~ 2
)
}
map2 <- leaflet() %>%
addTiles(group = "OpenStreetMap") %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
addProviderTiles(providers$Esri.WorldShadedRelief, group = "Topography") %>%
# GFW tile layers for three periods
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2001&end_year=2010&render_type=true_color",
group = "GFW Loss (2001-2010)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2011&end_year=2020&render_type=true_color",
group = "GFW Loss (2011-2020)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png?&start_year=2021&end_year=2024&render_type=true_color",
group = "GFW Loss (2021-2024)",
attribution = map_attribution
) %>%
addTiles(
"https://tiles.globalforestwatch.org/umd_regional_primary_forest_2001/latest/dynamic/{z}/{x}/{y}.png",
group = "Primary Forests (2001)",
attribution = map_attribution
) %>%
# Village polygon outlines (project villages highlighted)
addPolygons(
data = non_project,
fillColor = "transparent",
fillOpacity = 0,
color = "darkgrey",
weight = 0.5,
opacity = 0.4,
group = "Village Boundaries",
label = ~name
) %>%
addPolygons(
data = project_esri,
fillColor = "transparent",
fillOpacity = 0,
color = "#2166ac",
weight = 1.5,
opacity = 0.8,
group = "Village Boundaries",
label = ~name
) %>%
# Non-project circle markers (semi-transparent)
addCircleMarkers(
data = map_ctrl,
lng = ~lng, lat = ~lat,
color = ~pal_loss(rel_loss_cat),
radius = ~get_radius(abs_loss_cat),
popup = ~popup,
stroke = FALSE,
fillOpacity = 0.35,
group = "Non-Project Loss"
) %>%
# Project circle markers (full opacity, with dark border)
addCircleMarkers(
data = map_proj,
lng = ~lng, lat = ~lat,
color = ~pal_loss(rel_loss_cat),
radius = ~get_radius(abs_loss_cat),
popup = ~popup,
stroke = TRUE,
weight = 1,
opacity = 0.8,
fillOpacity = 0.85,
group = "Project Village Loss"
) %>%
addFullscreenControl() %>%
addLegend(
"bottomright",
pal = pal_loss,
values = village_map_data$rel_loss_cat,
title = "Relative Loss<br>(2001-2024)",
opacity = 1
) %>%
addLayersControl(
baseGroups = c("CartoDB", "OpenStreetMap", "Satellite", "Topography"),
overlayGroups = c("Project Village Loss", "Non-Project Loss",
"Village Boundaries",
"GFW Loss (2001-2010)", "GFW Loss (2011-2020)",
"GFW Loss (2021-2024)", "Primary Forests (2001)"),
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(c("GFW Loss (2001-2010)", "GFW Loss (2011-2020)",
"GFW Loss (2021-2024)", "Primary Forests (2001)"))
map2
}
```
## Summary and Implications
```{r summary-stats-final}
#| echo: false
#| results: asis
if (nrow(gfw_stats) > 0) {
total_proj_2001 <- sum(proj_stats$area_ha[proj_stats$year == 2001])
total_proj_2024 <- sum(proj_stats$area_ha[proj_stats$year == 2024])
total_proj_loss <- total_proj_2001 - total_proj_2024
total_proj_rel <- total_proj_loss / total_proj_2001 * 100
total_all_2001 <- sum(baseline$area_ha)
total_all_2024 <- sum(latest$area_ha)
total_all_loss <- total_all_2001 - total_all_2024
total_all_rel <- total_all_loss / total_all_2001 * 100
cat(paste0(
"### Key Figures\n\n",
"| Metric | Project Villages (83) | All Villages (921) |\n",
"|--------|----------------------|--------------------|\n",
"| Forest cover 2001 | ", format(round(total_proj_2001, 0), big.mark = ","), " ha | ",
format(round(total_all_2001, 0), big.mark = ","), " ha |\n",
"| Forest cover 2024 | ", format(round(total_proj_2024, 0), big.mark = ","), " ha | ",
format(round(total_all_2024, 0), big.mark = ","), " ha |\n",
"| Total loss | ", format(round(total_proj_loss, 0), big.mark = ","), " ha | ",
format(round(total_all_loss, 0), big.mark = ","), " ha |\n",
"| Relative loss | ", round(total_proj_rel, 1), "% | ",
round(total_all_rel, 1), "% |\n\n"
))
n_severe <- sum(village_summary$rel_loss[village_summary$is_project_village == TRUE] > 30, na.rm = TRUE)
n_moderate <- sum(village_summary$rel_loss[village_summary$is_project_village == TRUE] > 20 &
village_summary$rel_loss[village_summary$is_project_village == TRUE] <= 30, na.rm = TRUE)
n_low <- sum(village_summary$rel_loss[village_summary$is_project_village == TRUE] <= 10, na.rm = TRUE)
cat(paste0(
"### Implications for Conservation Zoning\n\n",
"The accelerating loss trend raises important questions for the project's conservation components:\n\n",
"- **", n_severe, " villages** have lost more than 30% of their 2001 forest cover. In these villages, ",
"the remaining forest may be too fragmented or reduced for effective conservation zoning without ",
"complementary restoration measures.\n",
"- **", n_moderate, " villages** show moderate loss (20--30%). Conservation zoning remains viable but ",
"urgent --- continued loss at current rates would further reduce the forest base within a few years.\n",
"- **", n_low, " villages** retain more than 90% of their 2001 forest and may offer the best prospects ",
"for effective conservation set-asides.\n\n",
"The fact that loss has accelerated in the most recent period (2021--2024) underscores the urgency ",
"of implementing land-use zoning measures. Without intervention, the forest base available for ",
"conservation in many project villages will continue to shrink rapidly.\n"
))
}
```
## Data Sources
| Dataset | Description | Source |
|---------|-------------|--------|
| GFW Treecover / Lossyear | Annual forest cover, version GFC-2024-v1.12 | [Global Forest Watch](https://www.globalforestwatch.org/) |
| ESRI India Living Atlas | Village boundaries (Census 2021), IAB Village 2024 | [ESRI Living Atlas India](https://livingatlas.esri.in/) |
| geoBoundaries ADM5 | Village boundaries (Census 2011) | [geoBoundaries](https://www.geoboundaries.org/) |
| India State of Forest Report 2023 | National forest statistics | [Forest Survey of India](https://fsi.nic.in/) |
## References
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 (6160): 850--53.