---
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 of 100 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 Data -----
if (nrow(gfw_stats) > 0) {
# Add a unique village identifier (vid) that works for both ESRI and
# geoBoundaries-direct (gb_direct) villages. ESRI villages use their
# numeric objectid; the 4 gb_direct project villages have objectid=NA
# and are identified by their csv_id instead.
gfw_stats <- gfw_stats %>%
mutate(vid = ifelse(!is.na(objectid),
as.character(objectid),
paste0("gb_", csv_id)))
# Compute a CORRECT per-village baseline from the year-2001 observation.
# NOTE: the baseline_2001 column in the CSV is WRONG for 53 non-project
# villages with lgd_villagecode=0 (they share the same group in the
# processing script, so first(area_ha) returns the first village's value
# for all 53). We recompute it here directly from the data.
true_baseline_lookup <- gfw_stats %>%
filter(year == 2001) %>%
select(vid, true_baseline_2001 = area_ha)
# Baseline year data (raw, for aggregate totals)
baseline <- gfw_stats %>% filter(year == 2001)
latest <- gfw_stats %>% filter(year == 2024)
# Project and control subsets
proj_stats <- gfw_stats %>% filter(is_project_village == TRUE)
ctrl_stats <- gfw_stats %>% filter(is_project_village == FALSE)
# Join confidence info for mapping
vdpic_matched <- vdpic %>% filter(confidence != "unmatched" | is.na(confidence))
# Per-village summary (uses correct per-objectid baseline from area_ha[year==2001])
village_summary <- gfw_stats %>%
group_by(vid, 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 — using vid as the unique key and the
# CORRECTED true_baseline_2001 for relative loss calculations.
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(vid, name, csv_name, is_project_village, period, period_pos, area_ha) %>%
pivot_wider(names_from = period_pos, values_from = area_ha) %>%
left_join(true_baseline_lookup, by = "vid") %>%
mutate(
period_loss_ha = start - end,
period_loss_rel = period_loss_ha / true_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 **100 villages** in the Dhalai and North Tripura districts. The project includes agricultural support, land-use zoning, and the designation of conservation set-aside areas. Of the 100 project villages, **83 could be matched to official Census boundary polygons** and are included in this analysis.
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). The project supports **100 villages** in the Dhalai and North Tripura districts; of these, **83 could be matched to official Census boundary polygons** and are included in this analysis. The remaining 17 villages could not be matched due to naming ambiguity, missing Census records, or being sub-hamlets not delineated as separate Census units.
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). Seventeen villages from the field programme could not be matched to any Census polygon — these include sub-hamlets, reserve forest sections, newly settled areas not recorded as separate Census units, and villages where name ambiguity made any reliable match impossible.
### 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 (83 of 100)"
)
}
```
**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)
proj_objectids <- gfw_stats %>%
filter(is_project_village == TRUE & year == 2001 & !is.na(objectid)) %>%
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 project village,
# using the CORRECT per-village 2001 baseline (area_ha at year==2001).
proj_trajectory <- proj_stats %>%
group_by(vid) %>%
mutate(
true_baseline = area_ha[year == 2001],
rel_forest_remaining = area_ha / true_baseline * 100,
cum_rel_loss = (true_baseline - area_ha) / true_baseline * 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),
.groups = "drop"
)
# Period-level data for project villages only
proj_periods <- period_stats %>%
filter(is_project_village == TRUE) %>%
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,
# Label for hover: prefer csv_name (project), fall back to name
village_label = coalesce(csv_name, name)
)
# Build comprehensive download table for all analysed villages
period_wide <- period_stats %>%
select(vid, period, period_loss_ha, period_loss_rel) %>%
mutate(
loss_ha_col = paste0("loss_ha_", gsub("-", "_", period)),
loss_rel_col = paste0("loss_pct_", gsub("-", "_", period))
) %>%
select(vid, loss_ha_col, loss_rel_col, period_loss_ha, period_loss_rel) %>%
pivot_longer(c(loss_ha_col, loss_rel_col), names_to = "col_type", values_to = "col_name") %>%
mutate(value = ifelse(col_type == "loss_ha_col", period_loss_ha, period_loss_rel)) %>%
select(vid, col_name, value) %>%
pivot_wider(names_from = col_name, values_from = value)
download_data <- village_summary %>%
select(vid, name, csv_name, csv_district, csv_block,
confidence, is_project_village, forest_2001, forest_2024, abs_loss, rel_loss) %>%
left_join(period_wide, by = "vid") %>%
mutate(
`Village Name` = coalesce(csv_name, name),
`Project Village` = ifelse(is_project_village, "Yes", "No"),
`Forest 2001 (ha)` = round(forest_2001, 1),
`Forest 2024 (ha)` = round(forest_2024, 1),
`Total Loss (ha)` = round(abs_loss, 1),
`Total Loss (%)` = round(rel_loss, 1),
`Loss 2001-2010 (ha)` = round(loss_ha_2001_2010, 1),
`Loss 2001-2010 (%)` = round(loss_pct_2001_2010, 1),
`Loss 2011-2020 (ha)` = round(loss_ha_2011_2020, 1),
`Loss 2011-2020 (%)` = round(loss_pct_2011_2020, 1),
`Loss 2021-2024 (ha)` = round(loss_ha_2021_2024, 1),
`Loss 2021-2024 (%)` = round(loss_pct_2021_2024, 1),
District = csv_district,
Confidence = confidence
) %>%
arrange(desc(is_project_village), desc(`Total Loss (%)`)) %>%
select(`Village Name`, District, Block = csv_block, `Project Village`, Confidence,
`Forest 2001 (ha)`, `Forest 2024 (ha)`,
`Total Loss (ha)`, `Total Loss (%)`,
`Loss 2001-2010 (ha)`, `Loss 2001-2010 (%)`,
`Loss 2011-2020 (ha)`, `Loss 2011-2020 (%)`,
`Loss 2021-2024 (ha)`, `Loss 2021-2024 (%)`)
# Period-level data for all villages (project + non-project)
# Uses the corrected true_baseline_2001 from period_stats
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")
) %>%
# Remove villages with implausible values: annual loss > 100% is physically
# impossible (a village cannot lose more than 100% of its baseline forest in a year).
# Such values signal data issues (e.g. villages whose geometry changed between
# census vintages, or very small baseline areas causing rounding artefacts).
filter(annual_rel_loss <= 100 & annual_rel_loss >= -10)
}
```
### 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) {
proj_periods$period <- factor(proj_periods$period,
levels = c("2001-2010", "2011-2020", "2021-2024"))
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"
)
period_colors <- c("2001-2010" = "#4575b4", "2011-2020" = "#fc8d59", "2021-2024" = "#d73027")
# Build hover text for individual village points
proj_periods <- proj_periods %>%
mutate(hover_label = paste0(
village_label,
"\nPeriod: ", period,
"\nAnnual loss: ", round(annual_rel_loss, 2), "%/yr"
))
p_box <- ggplot(proj_periods,
aes(x = period, y = annual_rel_loss, fill = period)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(aes(text = hover_label),
width = 0.15, alpha = 0.4, size = 1.5, 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) — hover for village name",
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, tooltip = "text") %>%
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
The plot below compares annualized relative loss between project and non-project villages across the same three periods. Villages with implausible annual loss values (above 100% of the 2001 baseline per year) — arising from a small number of census partition polygons with very small baseline areas — are excluded from this comparison.
```{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"))
# Determine y-axis range from the box whiskers (not outliers)
# Use 1.5*IQR rule to set sensible upper bound
p95 <- quantile(all_periods$annual_rel_loss, 0.95, na.rm = TRUE)
y_max <- ceiling(p95 * 1.2)
p_comp <- ggplot(all_periods, aes(x = period, y = annual_rel_loss, fill = village_type)) +
geom_boxplot(
alpha = 0.7,
outlier.shape = NA, # Hide outliers — they would dominate the scale
position = position_dodge(0.8)
) +
coord_cartesian(ylim = c(0, y_max)) + # Zoom axis without removing data
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 = paste0("Annualized loss as % of 2001 baseline, by period",
"\n(Outliers hidden; y-axis capped at ", y_max, "% for readability)"),
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 = 90),
legend = list(orientation = "h", x = 0.2, y = -0.2)
)
}
```
```{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"
))
}
```
### Forest Cover Loss: All Analysed Villages
The table below covers all 921 analysed villages, sorted by project village first and then by descending total forest cover loss. The **Total Loss (%)** column uses a color bar for quick visual comparison. Use the download buttons to export the full dataset as CSV or Excel.
```{r table-download-all}
#| echo: false
if (nrow(gfw_stats) > 0 && exists("download_data")) {
loss_range <- range(download_data$`Total Loss (%)`, na.rm = TRUE)
download_data %>%
datatable(
filter = 'top',
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
scrollX = TRUE,
pageLength = 10,
buttons = list(
list(extend = 'csv', text = 'Download CSV', filename = 'tripura_village_forest_loss'),
list(extend = 'excel', text = 'Download Excel', filename = 'tripura_village_forest_loss')
)
),
caption = "Forest cover loss by village and period — all 921 analysed villages (project villages first, sorted by total loss)"
) %>%
formatRound(
columns = c('Forest 2001 (ha)', 'Forest 2024 (ha)', 'Total Loss (ha)',
'Loss 2001-2010 (ha)', 'Loss 2011-2020 (ha)', 'Loss 2021-2024 (ha)'),
digits = 0
) %>%
formatStyle(
'Total Loss (%)',
background = styleColorBar(loss_range, '#fc8d59'),
backgroundSize = '98% 88%',
backgroundRepeat = 'no-repeat',
backgroundPosition = 'center'
) %>%
formatStyle(
'Project Village',
backgroundColor = styleEqual(c('Yes', 'No'), c('#e8f4fd', 'white'))
)
}
```
### 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. The **Confidence Overlay** layer shows the village boundaries color-coded by match confidence (as in the first map above) and can be toggled on for reference.
```{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)
# Join village_summary by objectid (numeric); drop NA-objectid gb_direct villages
# from the map since we don't have their polygon centroids in all_villages
village_map_data <- village_centroids %>%
st_drop_geometry() %>%
left_join(
village_summary %>% filter(!is.na(objectid)),
by = "objectid"
) %>%
bind_cols(
st_coordinates(village_centroids) %>% as_tibble() %>% rename(lng = X, lat = Y)
) %>%
filter(!is.na(rel_loss))
# 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")
)
map_proj <- village_map_data %>% filter(is_project_village == TRUE)
map_ctrl <- village_map_data %>% filter(is_project_village == FALSE)
# Build confidence overlay from ESRI polygons (same source as Village Boundaries layer)
# to ensure pixel-perfect alignment. Join confidence from gfw_stats by objectid.
village_conf_info <- gfw_stats %>%
filter(year == 2001, !is.na(objectid), is_project_village == TRUE) %>%
select(objectid, confidence, csv_name, csv_block, match_method) %>%
distinct(objectid, .keep_all = TRUE)
esri_conf_confirmed <- all_villages %>%
inner_join(village_conf_info %>% filter(confidence == "confirmed"),
by = "objectid")
esri_conf_likely <- all_villages %>%
inner_join(village_conf_info %>% filter(confidence == "likely"),
by = "objectid")
esri_conf_uncertain <- all_villages %>%
inner_join(village_conf_info %>% filter(confidence == "uncertain"),
by = "objectid")
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), "%)"
)
pal_loss <- colorFactor(
palette = rev(brewer.pal(5, "RdYlGn")),
domain = village_map_data$rel_loss_cat
)
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
) %>%
# Plain village polygon outlines (all grey, project villages slightly thicker)
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
) %>%
# Confidence overlay (hidden by default) — uses ESRI polygons (same base as
# "Village Boundaries") to ensure pixel-perfect alignment between layers.
addPolygons(
data = non_project,
fillColor = "transparent",
fillOpacity = 0,
color = "darkgrey",
weight = 1,
opacity = 0.5,
group = "Confidence Overlay",
label = ~name
) %>%
addPolygons(
data = esri_conf_confirmed,
fillColor = "transparent",
fillOpacity = 0,
color = "green",
weight = 2,
opacity = 0.9,
group = "Confidence Overlay",
label = ~csv_name
) %>%
addPolygons(
data = esri_conf_likely,
fillColor = "transparent",
fillOpacity = 0,
color = "blue",
weight = 2,
opacity = 0.9,
group = "Confidence Overlay",
label = ~csv_name
) %>%
addPolygons(
data = esri_conf_uncertain,
fillColor = "transparent",
fillOpacity = 0,
color = "orange",
weight = 2,
opacity = 0.9,
group = "Confidence Overlay",
label = ~csv_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", "Confidence Overlay",
"GFW Loss (2001-2010)", "GFW Loss (2011-2020)",
"GFW Loss (2021-2024)", "Primary Forests (2001)"),
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(c("Confidence Overlay",
"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 of 100) | 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(proj_summary$rel_loss > 30, na.rm = TRUE)
n_moderate <- sum(proj_summary$rel_loss > 20 & proj_summary$rel_loss <= 30, na.rm = TRUE)
n_low <- sum(proj_summary$rel_loss <= 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.