2  Fuel Stratification

Overview

This chapter describes the classification of burned areas into IPCC vegetation strata using MODIS Land Cover (MCD12Q1) and auxiliary datasets. The classification determines which emission factors and fuel consumption values apply to each burned pixel, directly impacting greenhouse gas emission estimates.

Environment Setup

easypackages::packages(
  "bslib",
  "cols4all", "covr", "cowplot",
  "dendextend", "digest","DiagrammeR","dtwclust", "downlit",
  "e1071", "exactextractr","elevatr",
  "FNN", "future", "forestdata",
  "gdalcubes", "gdalUtilities", "geojsonsf", "geos", "ggplot2", "ggstats",
  "ggspatial", "ggmap", "ggplotify", "ggpubr", "ggrepel", "giscoR",
  "hdf5r", "httr", "httr2", "htmltools",
  "jsonlite",
  "kohonen", "knitr",
  "leaflet.providers", "leafem", "libgeos","luz","lwgeom", "leaflet", "leafgl",
  "mapedit", "mapview", "maptiles", "methods", "mgcv",
  "ncdf4", "nnet",
  "openxlsx", "parallel", "plotly",
  "randomForest", "rasterVis", "raster", "Rcpp", "RcppArmadillo",
  "RcppCensSpatial","rayshader", "RcppEigen", "RcppParallel",
  "RColorBrewer", "reactable", "rgl", "rsconnect","RStoolbox", "rts", "reticulate",
  "s2", "sf", "scales", "sits","spdep", "stars", "stringr","supercells",
  "terra", "testthat", "tidyverse", "tidyterra","tools",
  "tmap", "tmaptools", "terrainr",
  "xgboost",
  prompt = F)

# Disable spherical/obloid geometry 
sf::sf_use_s2(use_s2 = F)

# Assign working directory to current location
knitr::opts_knit$set(root.dir = here::here()) 

# Point environment to your configured python
reticulate::use_python("/opt/local/bin/python", required=T)

# Initialise Earth Engine
ee <- reticulate::import("ee")
ee$Initialize(project = "murphys-deforisk")
cat("GEE initialised\n")
NA GEE initialised

# Helper Function to visualize EE tile in tmap/leaflet
ee_tile_url <- function(ee_image, vis_params) {
  ee_image$getMapId(vis_params)$tile_fetcher$url_format}

# EE FeatureCollection => sf
ee_to_sf <- function(fc) {
  geojsonsf::geojson_sf(
    jsonlite::toJSON(fc$getInfo(), auto_unbox = TRUE))}

---

2.1 IPCC Vegetation Strata

FAOSTAT implements three primary fire categories aligned with IPCC 2006 Guidelines (IPCC, 2006):

1. Forest Fires:

• Humid Tropical Forest: Fires in tropical moist ecological zones
• Other Forest: All remaining forest types (boreal, temperate, tropical dry)

2. Savanna Fires (Grassland):

• Grassland
• Savanna
• Woody Savanna
• Open Shrubland
• Closed Shrubland

3. Fires in Organic Soils:

• Peatland fires (Histosols classification)
• Independent of overlying vegetation

2.1.1 Classification Decision Tree

# Conceptual decision tree (not executable code)
#
# For each burned pixel:
#   ├─ Is pixel on Histosol?
#   │   └─ YES → Organic Soil Fire
#   │   └─ NO → Continue
#   ├─ Is land cover = Forest (Classes 1-5)?
#   │   ├─ YES → Is pixel in Tropical Moist zone?
#   │   │   ├─ YES → Humid Tropical Forest Fire
#   │   │   └─ NO → Other Forest Fire
#   │   └─ NO → Continue
#   ├─ Is land cover = Grassland/Savanna/Shrubland (Classes 6-10)?
#   │   └─ YES → Savanna Fire
#   └─ ELSE → Not classified (excluded from analysis)

---

2.2 GEE Initialisation

# Point environment to your configured python
reticulate::use_python("/opt/local/bin/python", required=T)

# Initialise Earth Engine
ee <- reticulate::import("ee")
ee$Initialize(project = "murphys-deforisk")
cat("GEE initialised\n")
## GEE initialised

# Helper: visualize EE tile in tmap/leaflet
ee_tile_url <- function(ee_image, vis_params) {
  ee_image$getMapId(vis_params)$tile_fetcher$url_format
}

# Helper: EE FeatureCollection => sf
ee_to_sf <- function(fc) {
  geojsonsf::geojson_sf(
    jsonlite::toJSON(fc$getInfo(), auto_unbox = TRUE))
}

# Load AOI boundaries from FAO GAUL
aoi_country_ee <- ee$FeatureCollection("FAO/GAUL/2015/level0")$
  filter(ee$Filter$eq("ADM0_NAME", "Honduras"))
aoi_states_ee  <- ee$FeatureCollection("FAO/GAUL/2015/level1")$
  filter(ee$Filter$eq("ADM0_NAME", "Honduras"))
aoi_country_sf <- ee_to_sf(aoi_country_ee)
aoi_states_sf  <- ee_to_sf(aoi_states_ee)

cat("AOI loaded:", aoi_country_ee$size()$getInfo(), "country,",
    aoi_states_ee$size()$getInfo(), "states\n")
## AOI loaded: 1 country, 19 states

---

2.3 MODIS Land Cover: MCD12Q1

2.3.1 IGBP Classification Scheme

2.3.1.1 Product Overview

• GEE Asset: MODIS/061/MCD12Q1
• Spatial Resolution: 500m
• Temporal Resolution: Annual (2001-present)
• Classification System: Land Cover Type 1 - IGBP (International Geosphere-Biosphere Programme)
• Algorithm: Supervised classification using Random Forest trained on reference sites

2.3.2 IGBP Burnable Classes

Forest Classes (1-5):

Class	Name	IPCC Category	Fuel Consumption (t ha⁻¹)
1	Evergreen Needleleaf Forests	Other Forest	41.0 (mean)
2	Evergreen Broadleaf Forests	Humid Tropical Forest	119.6 (mean)
3	Deciduous Needleleaf Forests	Other Forest	41.0 (mean)
4	Deciduous Broadleaf Forests	Other/Tropical (zone-dependent)	69.4-119.6
5	Mixed Forests	Other Forest	50.4 (mean)

Grassland/Savanna Classes (6-10):

Class	Name	IPCC Category	Fuel Consumption (t ha⁻¹)
6	Closed Shrublands	Savanna (Shrubland)	14.3 (mean)
7	Open Shrublands	Savanna (Shrubland)	14.3 (mean)
8	Woody Savannas	Savanna (Woody)	2.6-4.6 (season)
9	Savannas	Savanna	2.1-10.0 (season)
10	Grasslands	Savanna (Grassland)	2.1-10.0 (season)

Excluded Classes:

Class	Name	Reason for Exclusion
0	Water	Not burnable
11	Permanent Wetlands	Handled separately if Histosol
12	Croplands	Agricultural fires (separate category)
13	Urban and Built-up	Fire prevention/suppression active
14	Cropland/Natural Vegetation Mosaics	Ambiguous, often excluded
15	Permanent Snow and Ice	Not burnable
16	Barren	Insufficient fuel
17	Unclassified	No data

2.3.3 Earth Engine: MCD12Q1

# Load land cover collection
landCover <- ee$ImageCollection("MODIS/061/MCD12Q1")

# Select IGBP classification (Type 1) for 2020 (matching Ch1 burn year)
lc_2020 <- landCover$
  select("LC_Type1")$
  filterDate("2020-01-01", "2020-12-31")$
  first()$
  clip(aoi_country_ee$geometry())

cat("Land cover bands:", paste(lc_2020$bandNames()$getInfo(),
    collapse = ", "), "\n")
NA Land cover bands: LC_Type1

# IGBP palette
lc_vis <- list(
  min = 0L,
  max = 17L,
  palette = c(
    '05450a', '086a10', '54a708', '78d203', '009900',  # Forests 1-5
    'c6b044', 'dcd159', 'dade48', 'fbff13', 'b6ff05',  # Savanna/Grass 6-10
    '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c', '69fff8',  # Other 11-15
    'f9ffa4', '1c0dff'                                 # Barren/Unclass 16-17
  )
)

tmap::tmap_mode("view")
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(ee_tile_url(lc_2020, lc_vis),group = "IGBP Land Cover 2020") +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("#05450a", "#086a10", "#54a708", "#78d203", "#009900",
             "#c6b044", "#dcd159", "#dade48", "#fbff13", "#b6ff05",
             "#27ff87", "#c24f44", "#a5a5a5", "#ff6d4c", "#69fff8",
             "#f9ffa4", "#1c0dff"),
    labels = c("1 Evergreen Needleleaf", "2 Evergreen Broadleaf",
               "3 Deciduous Needleleaf", "4 Deciduous Broadleaf",
               "5 Mixed Forests", "6 Closed Shrublands",
               "7 Open Shrublands", "8 Woody Savannas",
               "9 Savannas", "10 Grasslands",
               "11 Permanent Wetlands", "12 Croplands",
               "13 Urban", "14 Cropland/Veg Mosaic",
               "15 Snow/Ice", "16 Barren", "17 Unclassified"),
    title = "IGBP Land Cover (2020)") +
    tmap::tm_basemap("Esri.WorldImagery")

2.3.4 Land Cover Validation and Uncertainty

Known Classification Issues:

1. Woody Savanna vs. Open Forest: Confusion in transition zones
2. Grassland vs. Cropland: Agricultural areas misclassified
3. Seasonal Deciduous Forests: May shift between forest/savanna
4. Temporal Lag: Land use changes not immediately reflected

Regional Accuracy (Sulla-Menashe et al., 2018):

Region	Overall Accuracy	Forest Accuracy	Savanna Accuracy
Africa	73%	68%	78%
South America	77%	82%	71%
Europe	79%	86%	65%
Asia	74%	75%	73%
North America	81%	88%	74%

2.3.5 Land Cover Class Distribution

# IGBP class labels
igbp_labels <- c(
  "0"  = "Water",
  "1"  = "Evergreen Needleleaf",
  "2"  = "Evergreen Broadleaf",
  "3"  = "Deciduous Needleleaf",
  "4"  = "Deciduous Broadleaf",
  "5"  = "Mixed Forests",
  "6"  = "Closed Shrublands",
  "7"  = "Open Shrublands",
  "8"  = "Woody Savannas",
  "9"  = "Savannas",
  "10" = "Grasslands",
  "11" = "Permanent Wetlands",
  "12" = "Croplands",
  "13" = "Urban",
  "14" = "Cropland/Veg Mosaic",
  "15" = "Snow/Ice",
  "16" = "Barren",
  "17" = "Unclassified"
)

# Frequency histogram of IGBP classes
lc_hist <- lc_2020$reduceRegion(
  reducer  = ee$Reducer$frequencyHistogram(),
  geometry = aoi_country_ee$geometry(),
  scale    = 500L,
  maxPixels = 1e9
)$getInfo()

lc_counts <- lc_hist[["LC_Type1"]]
lc_df <- data.frame(
  class = as.integer(names(lc_counts)),
  pixels = as.integer(unlist(lc_counts))
)
lc_df$label <- igbp_labels[as.character(lc_df$class)]
lc_df$area_ha <- round(lc_df$pixels * 25, 1)  # 500m pixel ~ 25 ha
lc_df <- lc_df[order(-lc_df$pixels), ]

cat("Honduras 2020 — IGBP Land Cover Distribution\n")
NA Honduras 2020 — IGBP Land Cover Distribution
for (i in seq_len(nrow(lc_df))) {
  cat(sprintf("  Class %2d %-25s: %7d pixels (%9.0f ha)\n",
    lc_df$class[i], lc_df$label[i],
    lc_df$pixels[i], lc_df$area_ha[i]))
}
NA   Class  8 Woody Savannas           :  197605 pixels (  4940125 ha)
NA   Class  2 Evergreen Broadleaf      :  109633 pixels (  2740825 ha)
NA   Class  9 Savannas                 :   94455 pixels (  2361375 ha)
NA   Class 10 Grasslands               :   25527 pixels (   638175 ha)
NA   Class 12 Croplands                :    9225 pixels (   230625 ha)
NA   Class 11 Permanent Wetlands       :    4498 pixels (   112450 ha)
NA   Class 17 Unclassified             :    3402 pixels (    85050 ha)
NA   Class  4 Deciduous Broadleaf      :    2390 pixels (    59750 ha)
NA   Class 13 Urban                    :    1281 pixels (    32025 ha)
NA   Class  5 Mixed Forests            :    1076 pixels (    26900 ha)
NA   Class 14 Cropland/Veg Mosaic      :     879 pixels (    21975 ha)
NA   Class  6 Closed Shrublands        :      84 pixels (     2100 ha)
NA   Class 16 Barren                   :      18 pixels (      450 ha)
NA   Class  1 Evergreen Needleleaf     :      11 pixels (      275 ha)
NA   Class  7 Open Shrublands          :       2 pixels (       50 ha)

---

2.4 Forest Fire Classification

2.4.1 Forest Type Identification

# Define forest classes (IGBP 1-5)
forest_mask <- lc_2020$gte(1L)$And(lc_2020$lte(5L))$rename("forest_mask")

# Count forest pixels
forest_count <- forest_mask$selfMask()$reduceRegion(
  reducer   = ee$Reducer$count(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Forest pixels in Honduras:", forest_count[["forest_mask"]], "\n")
## Forest pixels in Honduras: 113127
cat("Forest area (approx):",
    round(as.integer(forest_count[["forest_mask"]]) * 25 / 1e3, 1),
    "thousand ha\n")
## Forest area (approx): 2828 thousand ha

# Visualize forest mask
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(ee_tile_url(forest_mask$selfMask()$clip(aoi_country_ee$geometry()),
    list(palette = "086a10", min = 1L, max = 1L)), group = "Forest (IGBP 1-5)") +
  tmap::tm_add_legend(
    type = "polygons",
    fill = "#086a10",
    labels = "Forest (IGBP 1-5)",
    title = "Forest Cover (2020)") +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.4.2 Humid Tropical Forest Delineation

2.4.2.1 FAO-FRA Global Ecological Zones

Dataset: FAO-FRA Global Ecological Zones 2020 Update

Tropical Moist Zones (for Humid Tropical Forest):

• Tropical rainforest
• Tropical moist deciduous forest
• Tropical mountain systems (in moist zones)

Implementation:

# Climate-based proxy for tropical moist zones
# Load WorldClim data
bio <- ee$Image("WORLDCLIM/V1/BIO")
precip <- bio$select("bio12")  # Annual precipitation (mm)
temp   <- bio$select("bio01")  # Annual mean temp (°C * 10)

# Define tropical moist criteria
# Tropical: Mean annual temp > 18°C (WorldClim stores as °C * 10)
# Moist: Annual precip > 1500mm
tropical <- temp$gt(180L)
moist    <- precip$gt(1500L)
humid_tropical_mask <- tropical$And(moist)

# Combine with forest mask
humid_tropical_forest <- forest_mask$And(humid_tropical_mask)
other_forest          <- forest_mask$And(humid_tropical_mask$Not())

# Count pixels
htf_count <- humid_tropical_forest$selfMask()$reduceRegion(
  reducer   = ee$Reducer$count(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

of_count <- other_forest$selfMask()$reduceRegion(
  reducer   = ee$Reducer$count(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Honduras — Forest Type Split\n")
## Honduras — Forest Type Split
cat("  Humid Tropical Forest:", htf_count[["forest_mask"]] %||% 0, "pixels\n")
##   Humid Tropical Forest: 89648 pixels
cat("  Other Forest:         ", of_count[["forest_mask"]] %||% 0, "pixels\n")
##   Other Forest:          23174 pixels

# Visualize
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      humid_tropical_forest$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "darkgreen", min = 1L, max = 1L)
    ), group = "Humid Tropical Forest"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      other_forest$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "orange", min = 1L, max = 1L)
    ), group = "Other Forest"
  ) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("darkgreen", "orange"),
    labels = c("Humid Tropical Forest", "Other Forest"),
    title = "Forest Type Classification"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.4.3 Forest Fire Area Calculation

# Load burned area from Chapter 1 (2020)
mcd64 <- ee$ImageCollection("MODIS/061/MCD64A1")
burned_2020 <- mcd64$
  filterBounds(aoi_country_ee$geometry())$
  filterDate("2020-01-01", "2020-12-31")

# FAOSTAT quality filter (matching Ch1)
faostat_filter <- function(image) {
  unc  <- image$select("Uncertainty")
  burn <- image$select("BurnDate")
  mask <- unc$lt(20L)$And(burn$gt(0L))
  image$updateMask(mask)
}

# Annual burned mask for 2020
annual_burned_2020 <- burned_2020$
  map(faostat_filter)$
  select("BurnDate")$max()$gt(0L)$selfMask()

# Overlay with vegetation masks
htf_fires_2020    <- annual_burned_2020$updateMask(humid_tropical_forest)
other_fires_2020  <- annual_burned_2020$updateMask(other_forest)

# Calculate areas (ha)
pixelArea_ha <- ee$Image$pixelArea()$divide(10000L)

htf_area <- htf_fires_2020$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

of_area <- other_fires_2020$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

cat("Honduras 2020 — Forest Fire Areas\n")
## Honduras 2020 — Forest Fire Areas
cat("  Humid Tropical Forest fires:",
    round(htf_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Humid Tropical Forest fires: 124750 ha
cat("  Other Forest fires:         ",
    round(of_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Other Forest fires:          43816 ha

# Visualize
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      htf_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "darkred", min = 0L, max = 1L)
    ), group = "Humid Tropical Forest Fires"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      other_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "orange", min = 0L, max = 1L)
    ), group = "Other Forest Fires"
  ) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("darkred", "orange"),
    labels = c("Humid Tropical Forest Fires", "Other Forest Fires"),
    title = "Forest Fire Areas (2020)"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

---

2.5 Savanna Fire Classification

2.5.1 Grassland and Woody Vegetation Types

# Savanna/grassland classes (IGBP 6-10)
savanna_mask <- lc_2020$gte(6L)$And(lc_2020$lte(10L))

# Individual sub-type masks
grassland_mask     <- lc_2020$eq(10L)
savanna_only_mask  <- lc_2020$eq(9L)
woody_savanna_mask <- lc_2020$eq(8L)
open_shrub_mask    <- lc_2020$eq(7L)
closed_shrub_mask  <- lc_2020$eq(6L)

# Count pixels by sub-type
savanna_stats <- ee$Image$cat(
  grassland_mask$rename("grassland"),
  savanna_only_mask$rename("savanna"),
  woody_savanna_mask$rename("woody_savanna"),
  open_shrub_mask$rename("open_shrub"),
  closed_shrub_mask$rename("closed_shrub")
)$reduceRegion(
  reducer   = ee$Reducer$sum(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Honduras 2020 — Savanna/Grassland Sub-Types\n")
## Honduras 2020 — Savanna/Grassland Sub-Types
cat("  Grasslands:       ", savanna_stats$grassland, "pixels\n")
##   Grasslands:        25527 pixels
cat("  Savannas:         ", savanna_stats$savanna, "pixels\n")
##   Savannas:          94455 pixels
cat("  Woody Savannas:   ", savanna_stats$woody_savanna, "pixels\n")
##   Woody Savannas:    197605 pixels
cat("  Open Shrublands:  ", savanna_stats$open_shrub, "pixels\n")
##   Open Shrublands:   2 pixels
cat("  Closed Shrublands:", savanna_stats$closed_shrub, "pixels\n")
##   Closed Shrublands: 84 pixels

# Overlay with burned area
savanna_fires_2020 <- annual_burned_2020$updateMask(savanna_mask)

sav_area <- savanna_fires_2020$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

cat("  Total savanna fire area:", round(sav_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Total savanna fire area: 559731 ha

# Visualize
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      savanna_mask$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "fbff13", min = 1L, max = 1L)
    ), group = "Savanna/Grassland (IGBP 6-10)"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      savanna_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "red", min = 0L, max = 1L)
    ), group = "Savanna Fires 2020"
  ) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("#fbff13", "red"),
    labels = c("Savanna/Grassland (IGBP 6-10)", "Savanna Fires 2020"),
    title = "Savanna Classification & Fires"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.5.2 Climate Zone Stratification

Rationale: Fuel consumption varies by climate (IPCC Table 2.4)

IPCC Climate Zones:

• Boreal (Cold)
• Temperate (Warm/Cool)
• Tropical (Wet/Moist/Dry/Montane)

# Climate zones from WorldClim (already loaded)
# Tropical: MAT > 18°C
tropical_zone  <- temp$gt(180L)
# Boreal: MAT < 0°C
boreal_zone    <- temp$lt(0L)
# Temperate: Everything else
temperate_zone <- temp$gte(0L)$And(temp$lte(180L))

# Climate zone pixel counts within Honduras
zone_stats <- ee$Image$cat(
  tropical_zone$rename("tropical"),
  temperate_zone$rename("temperate"),
  boreal_zone$rename("boreal")
)$reduceRegion(
  reducer   = ee$Reducer$sum(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Honduras — Climate Zone Distribution\n")
## Honduras — Climate Zone Distribution
cat("  Tropical pixels: ", zone_stats$tropical, "\n")
##   Tropical pixels:  449368
cat("  Temperate pixels:", zone_stats$temperate, "\n")
##   Temperate pixels: 15854
cat("  Boreal pixels:   ", zone_stats$boreal, "\n")
##   Boreal pixels:    0

# Apply climate stratification to savanna fires
savanna_tropical  <- savanna_fires_2020$updateMask(tropical_zone)
savanna_temperate <- savanna_fires_2020$updateMask(temperate_zone)

sav_trop_area <- savanna_tropical$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

sav_temp_area <- savanna_temperate$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

cat("Honduras 2020 — Savanna Fires by Climate Zone\n")
## Honduras 2020 — Savanna Fires by Climate Zone
cat("  Tropical savanna fires: ",
    round(sav_trop_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Tropical savanna fires:  556845 ha
cat("  Temperate savanna fires:",
    round(sav_temp_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Temperate savanna fires: 2838 ha

# Visualize climate-stratified fires
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      tropical_zone$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "ff6600", min = 1L, max = 1L)
    ), group = "Tropical Zone"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      temperate_zone$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "3399ff", min = 1L, max = 1L)
    ), group = "Temperate Zone"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.5.3 Seasonal Burn Patterns

Importance: Fuel consumption differs between early/late dry season (IPCC Table 2.4)

• Early dry season: Lower fuel consumption (2.6 t ha⁻¹ savanna woodland)
• Late dry season: Higher fuel consumption (4.6 t ha⁻¹ savanna woodland)

# Extract burn month from burn date
extractBurnMonth <- function(mcd64_image) {
  burn_date <- mcd64_image$select("BurnDate")

  # Get image month from system:time_start
  img_date <- ee$Date(mcd64_image$get("system:time_start"))
  month    <- img_date$get("month")

  # Create month band (only where burned)
  month_band <- ee$Image$constant(month)$
    updateMask(burn_date$gt(0L))$
    rename("burn_month")$
    toInt()

  burn_date$addBands(month_band)
}

# Apply to 2020 collection
mcd64_with_months <- burned_2020$
  map(faostat_filter)$
  map(extractBurnMonth)

# Composite: latest month if pixel burned multiple times
month_composite <- mcd64_with_months$
  select("burn_month")$
  max()

# Honduras dry season: roughly Feb-May (Central America)
# Early dry: Feb-Mar (months 2-3)
# Late dry:  Apr-May (months 4-5)
early_season <- month_composite$gte(2L)$And(month_composite$lte(3L))
late_season  <- month_composite$gte(4L)$And(month_composite$lte(5L))

# Calculate areas for savanna fires by season
savanna_early <- savanna_fires_2020$updateMask(early_season)
savanna_late  <- savanna_fires_2020$updateMask(late_season)

early_area <- savanna_early$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

late_area <- savanna_late$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

cat("Honduras 2020 — Seasonal Savanna Fire Areas\n")
## Honduras 2020 — Seasonal Savanna Fire Areas
cat("  Early dry season (Feb-Mar):",
    round(early_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Early dry season (Feb-Mar): 32670 ha
cat("  Late dry season (Apr-May): ",
    round(late_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Late dry season (Apr-May):  523122 ha

---

2.6 Peatland Fire Classification

2.6.1 Histosols Identification

• Dataset: Harmonized World Soil Database (HWSD v1.2)
• Histosols: Organic soils with > 20-30% organic matter by weight
• Detection Challenge: Fires in organic soils independent of surface land cover

# Proxy approach: Permanent wetlands (IGBP 11) in lowland tropics
wetlands_mask <- lc_2020$eq(11L)

# Elevation threshold for lowlands
elevation <- ee$Image("USGS/SRTMGL1_003")
lowlands  <- elevation$lt(100L)

# Organic soil proxy = wetlands in lowlands
organic_soil_mask <- wetlands_mask$And(lowlands)

# Count in Honduras
org_count <- organic_soil_mask$selfMask()$reduceRegion(
  reducer   = ee$Reducer$count(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Honduras — Organic Soil Proxy\n")
## Honduras — Organic Soil Proxy
cat("  Lowland wetland pixels:", org_count[["LC_Type1"]] %||% 0, "\n")
##   Lowland wetland pixels: 4463

# Overlay with burned area
organic_fires_2020 <- annual_burned_2020$
  updateMask(organic_soil_mask)

org_area <- organic_fires_2020$
  multiply(pixelArea_ha)$
  reduceRegion(
    reducer   = ee$Reducer$sum(),
    geometry  = aoi_country_ee$geometry(),
    scale     = 500L,
    maxPixels = 1e9
  )$getInfo()

cat("  Organic soil fire area:", round(org_area[["BurnDate"]] %||% 0, 1), "ha\n")
##   Organic soil fire area: 1612 ha

2.6.2 Regional Constraints

FAOSTAT approach (per methodology guide):

“Due to the current lack of evidence for other regions of the world, we set to 0 all the country values outside Southeastern Asia.”

Southeast Asian Countries (non-zero organic soil fires):

• Indonesia
• Malaysia
• Brunei Darussalam

# For Honduras: organic soil fires = 0 per FAOSTAT methodology
# (Honduras is outside SE Asia)
cat("Honduras — FAOSTAT Organic Soil Fires\n")
## Honduras — FAOSTAT Organic Soil Fires
cat("  Per FAOSTAT methodology, organic soil fires are set to 0\n")
##   Per FAOSTAT methodology, organic soil fires are set to 0
cat("  for all countries outside Southeast Asia.\n")
##   for all countries outside Southeast Asia.
cat("  Honduras organic soil fire area: 0 ha (FAOSTAT convention)\n")
##   Honduras organic soil fire area: 0 ha (FAOSTAT convention)

2.6.3 Uncertainty in Organic Soil Classification

High uncertainty sources (IPCC 2014 Wetlands Supplement):

1. Detection limitations:

• Smoldering fires difficult to detect via optical remote sensing
• Dense canopy obscures ground fires

2. Soil classification uncertainty:

• HWSD spatial resolution (~1km) coarser than fire detection (500m)
• Organic soil depth variable within pixels
• Drainage status unknown (drained vs. undrained)

3. Fire depth uncertainty:

• Surface burn vs. deep peat burn
• Fuel consumption highly variable: 41 ± 1.4 t ha⁻¹ (IPCC Table 2.6)

FAOSTAT approach: Conservative reporting (SE Asia only) until more validation data available

---

2.7 Multi-Source Validation

2.7.1 GLAD Forest Height Validation

# GLAD canopy height in GEE 
glad_height <- ee$Image("UMD/hansen/global_forest_change_2024_v1_12")$
  select("treecover2000")$
  clip(aoi_country_ee$geometry())

# Forest threshold: > 25% canopy cover (Hansen definition)
glad_forest <- glad_height$gt(25L)

# Compare with MODIS forest mask
modis_forest_binary <- forest_mask$unmask(0L)
glad_forest_binary  <- glad_forest$unmask(0L)

agreement  <- modis_forest_binary$And(glad_forest_binary)
modis_only <- modis_forest_binary$And(glad_forest_binary$Not())
glad_only  <- glad_forest_binary$And(modis_forest_binary$Not())

val_stats <- ee$Image$cat(
  agreement$rename("agree"),
  modis_only$rename("modis_only"),
  glad_only$rename("glad_only")
)$reduceRegion(
  reducer   = ee$Reducer$sum(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

total <- as.integer(val_stats$agree %||% 0) +
  as.integer(val_stats$modis_only %||% 0) +
  as.integer(val_stats$glad_only %||% 0)

cat("Honduras — Forest Classification Agreement\n")
## Honduras — Forest Classification Agreement
cat("  Both agree (forest):  ", val_stats$agree, "pixels\n")
##   Both agree (forest):   110469 pixels
cat("  MODIS-only forest:    ", val_stats$modis_only, "pixels\n")
##   MODIS-only forest:     2643 pixels
cat("  Hansen-only forest:   ", val_stats$glad_only, "pixels\n")
##   Hansen-only forest:    236986 pixels
if (total > 0) {
  cat("  Agreement rate:       ",
      round(as.integer(val_stats$agree %||% 0) / total * 100, 1), "%\n")
}
##   Agreement rate:        31.6 %

# Visualize agreement
tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      agreement$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "green", min = 1L, max = 1L)
    ), group = "Both Agree: Forest"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      modis_only$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "blue", min = 1L, max = 1L)
    ), group = "MODIS-only Forest"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      glad_only$selfMask()$clip(aoi_country_ee$geometry()),
      list(palette = "red", min = 1L, max = 1L)
    ), group = "Hansen-only Forest"
  ) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("green", "blue", "red"),
    labels = c("Both Agree: Forest", "MODIS-only Forest", "Hansen-only Forest"),
    title = "Forest Classification Agreement"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.7.2 ESA WorldCover Validation

# ESA WorldCover 2020 (10m) in GEE
esa_2020 <- ee$Image("ESA/WorldCover/v100/2020")$
  clip(aoi_country_ee$geometry())

# ESA classes: 10=Tree cover, 20=Shrubland, 30=Grassland,
# 40=Cropland, 50=Built-up, 60=Bare, 70=Snow/Ice,
# 80=Water, 90=Herbaceous wetland, 95=Mangroves, 100=Moss/Lichen

# ESA forest (Tree cover)
esa_forest <- esa_2020$eq(10L)

# ESA savanna proxy (Shrubland + Grassland)
esa_savanna <- esa_2020$eq(20L)$Or(esa_2020$eq(30L))

# Compare ESA forest with MODIS forest
esa_forest_binary <- esa_forest$unmask(0L)

agree_esa  <- modis_forest_binary$And(esa_forest_binary)
modis_esa  <- modis_forest_binary$And(esa_forest_binary$Not())
esa_only   <- esa_forest_binary$And(modis_forest_binary$Not())

esa_stats <- ee$Image$cat(
  agree_esa$rename("agree"),
  modis_esa$rename("modis_only"),
  esa_only$rename("esa_only")
)$reduceRegion(
  reducer   = ee$Reducer$sum(),
  geometry  = aoi_country_ee$geometry(),
  scale     = 500L,
  maxPixels = 1e9
)$getInfo()

cat("Honduras — MODIS vs ESA WorldCover Forest Agreement\n")
cat("  Both agree (forest): ", esa_stats$agree, "pixels\n")
cat("  MODIS-only forest:   ", esa_stats$modis_only, "pixels\n")
cat("  ESA-only forest:     ", esa_stats$esa_only, "pixels\n")

---

2.8 Uncertainty & QA/QC

2.8.1 Accuracy Assessment

Sources of uncertainty:

1. Land cover classification error: ±10-15% (MCD12Q1 validation)
2. Burned area overlay error: ±2-5% (geometric registration)
3. Ecological zone delineation: ±5-10% (climate proxy vs. official zones)
4. Temporal mismatch: Land cover annual vs. monthly burns

Total uncertainty (quadrature): \[ U_{classification}^2 = U_{landcover}^2 + U_{overlay}^2 + U_{ecozone}^2 \]

Typical result: ±15-20% uncertainty in vegetation type assignment

2.8.2 QC Checklist

Pre-classification:

• Burned area and land cover from same year
• Land cover projection matches burned area
• All forest classes (1-5) included in forest mask
• All savanna classes (6-10) included in savanna mask

Classification:

• Mutually exclusive categories (no pixel in multiple categories)
• Complete coverage (all burned pixels classified)
• Climate zones applied correctly
• Organic soils restricted to validated regions

Post-classification:

• Area totals match Chapter 1 burned area
• No negative or unrealistic values
• Results comparable to previous years (if available)
• Validation with independent data completed

---

2.9 Reproducible Workflow

2.9.1 Consolidated Script

# IPCC Tier 1 Vegetation Classification Pipeline
# ============================================================\


# 1. Classify Burns
# Unclassified: burned but not forest or savanna
unclassified_mask <- forest_mask$Or(savanna_mask)$Not()
unclassified_fires <- annual_burned_2020$updateMask(unclassified_mask)

# 2. Areas Estimates
calcArea <- function(fire_img, label) {
  result <- fire_img$
    multiply(pixelArea_ha)$
    reduceRegion(
      reducer   = ee$Reducer$sum(),
      geometry  = aoi_country_ee$geometry(),
      scale     = 500L,
      maxPixels = 1e9
    )$getInfo()
  round(result[["BurnDate"]] %||% 0, 1)
}

total_burned_ha <- calcArea(annual_burned_2020, "Total")
htf_ha          <- calcArea(htf_fires_2020, "HTF")
of_ha           <- calcArea(other_fires_2020, "Other Forest")
sav_ha          <- calcArea(savanna_fires_2020, "Savanna")
unclass_ha      <- calcArea(unclassified_fires, "Unclassified")

cat("Honduras 2020 — IPCC Vegetation Classification Summary\n")
## Honduras 2020 — IPCC Vegetation Classification Summary
cat(sprintf("  %-30s %10.1f ha\n", "Total burned area:", total_burned_ha))
##   Total burned area:               747693.6 ha
cat(sprintf("  %-30s %10.1f ha\n", "Humid Tropical Forest fires:", htf_ha))
##   Humid Tropical Forest fires:     124749.6 ha
cat(sprintf("  %-30s %10.1f ha\n", "Other Forest fires:", of_ha))
##   Other Forest fires:               43816.5 ha
cat(sprintf("  %-30s %10.1f ha\n", "Savanna fires:", sav_ha))
##   Savanna fires:                   559731.0 ha
cat(sprintf("  %-30s %10.1f ha\n", "Organic soil fires:", 0))
##   Organic soil fires:                   0.0 ha
cat(sprintf("  %-30s %10.1f ha\n", "Unclassified (excluded):", unclass_ha))
##   Unclassified (excluded):          19396.4 ha
classified_ha <- htf_ha + of_ha + sav_ha
cat(sprintf("  %-30s %10.1f ha\n", "Total classified:", classified_ha))
##   Total classified:                728297.1 ha
if (total_burned_ha > 0) {cat(sprintf("  %-30s %9.1f%%\n", "Classification coverage:",
      classified_ha / total_burned_ha * 100))}
##   Classification coverage:            97.4%

# 3. Visualize
tmap::tm_shape(aoi_states_sf) + tmap::tm_borders(col = "MAP_COLORS", lwd = 1) +
  tmap::tm_shape(aoi_country_sf) + tmap::tm_borders(col = "white", lwd = 1.5) +
  tmap::tm_tiles(
    ee_tile_url(
      htf_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "darkred", min = 0L, max = 1L)
    ), group = "Humid Tropical Forest Fires"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      other_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "ff8800", min = 0L, max = 1L)
    ), group = "Other Forest Fires"
  ) +
  tmap::tm_tiles(
    ee_tile_url(
      savanna_fires_2020$clip(aoi_country_ee$geometry()),
      list(palette = "ffff00", min = 0L, max = 1L)
    ), group = "Savanna Fires"
  ) +
  tmap::tm_add_legend(
    type = "polygons",
    fill = c("darkred", "#ff8800", "#ffff00"),
    labels = c("Humid Tropical Forest Fires", "Other Forest Fires", "Savanna Fires"),
    title = "IPCC Fire Classification (2020)"
  ) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = 0.5) +
  tmap::tm_basemap("Esri.WorldImagery")

2.9.2 Output Format

Required outputs for fuel biomass calculation (Chapter 3):

# Build results data frame
output_table <- data.frame(
  country = "Honduras",
  year = 2020L,
  fire_type = c("Humid_Tropical_Forest", "Other_Forest",
                "Savanna", "Organic_Soil"),
  burned_area_ha = c(htf_ha, of_ha, sav_ha, 0),
  climate_zone = c("Tropical", "Tropical", "Tropical", NA),
  fuel_consumption_tDM_ha = c(119.6, 69.4, 7.0, 0),
  notes = c(
    "IGBP 1-5 in humid tropical zone",
    "IGBP 1-5 outside humid tropical zone",
    "IGBP 6-10 (all sub-types)",
    "Set to 0 per FAOSTAT (non-SE Asia)"
  )
)

knitr::kable(output_table,caption = "2020 Vegetation Classification Output",digits=1)

2020 Vegetation Classification Output

country	year	fire_type	burned_area_ha	climate_zone	fuel_consumption_tDM_ha	notes
Honduras	2020	Humid_Tropical_Forest	124750	Tropical	119.6	IGBP 1-5 in humid tropical zone
Honduras	2020	Other_Forest	43816	Tropical	69.4	IGBP 1-5 outside humid tropical zone
Honduras	2020	Savanna	559731	Tropical	7.0	IGBP 6-10 (all sub-types)
Honduras	2020	Organic_Soil	0	NA	0.0	Set to 0 per FAOSTAT (non-SE Asia)

---

2.10 Notes & Next Steps

1. MODIS MCD12Q1 IGBP classification provides globally consistent vegetation types
2. Three IPCC fire categories: Forest (Humid Tropical + Other), Savanna, Organic Soils
3. Climate stratification required for accurate fuel consumption assignment
4. Validation using Hansen forest cover and ESA WorldCover reduces uncertainty
5. Uncertainty ~±15-20% in vegetation classification acceptable for Tier 1

Chapter 3: Fuel Biomass Consumption will use vegetation classifications to:

1. Assign IPCC default fuel consumption values (Table 2.4) by vegetation type
2. Apply climate-dependent stratification
3. Calculate total burned biomass: \(A \times M_B\) (Activity Data × Fuel Biomass)
4. Prepare for emission factor application (Chapter 4)

Critical linkage: Each vegetation type carries forward:

• Area burned (hectares) by vegetation type
• Climate zone (for fuel consumption stratification)
• Seasonal information (early/late dry season for savannas)
• Geographic location (for validation and reporting)

---

IPCC. (2006). Chapter 3: Consistent representation of lands. In 2006 IPCC guidelines for national greenhouse gas inventories 3.1 (Agriculture, Forestry and Other Land Use, Vol. 4). Intergovernmental Panel on Climate Change. https://www.ipcc-nggip.iges.or.jp/public/2006gl/pdf/4_Volume4/V4_03_Ch3_Representation.pdf

Sulla-Menashe, D., Friedl, M. A., et al. (2018). User guide to collection 6 MODIS land cover (MCD12Q1 and MCD12C1) product. USGS: Reston, Va, USA, 1, 18.