VM0048: Jurisidictional Deforestation Risk Allocation

Author

Seamus Murphy

Published

January 4, 2025

Abstract

A workflow for deriving jurisidictional risk-allocated deforestation mapping compliant with Verra’s VMD0055 (V1.1) module and the VM0048 (V1.0) consolidated methodology. This includes baseline estimates of forest change and jurisdictional allocated deforestation.

Keywords

REDD+, VCS, Verra, Carbon verification, Jurisdictional

1. Introduction

This assessment develops a spatially explicit deforestation risk map and allocation framework for Bong County, Liberia, following Verra’s VM0048 methodology requirements. The analysis integrates the baseline emissions estimates developed in the previous workflow with spatial risk modeling to allocate jurisdictional deforestation across the landscape based on empirically-derived risk factors.

This analysis maintains consistency with the baseline emissions assessment by focusing exclusively on Bong County (8,772 km²), located in west-central Liberia. The county encompasses diverse forest-agricultural landscapes within the Upper Guinea Forest biodiversity hotspot, representing typical West African deforestation pressures including:

  • Agricultural expansion (primarily rice and cassava)
  • Small-scale logging and charcoal production
  • Road development and settlement expansion
  • Mining activities (iron ore and artisan mining)

2. Method

Four-step process:

  1. Covariate Development: Process infrastructure, demographic, and biophysical risk data layers
  2. Statistical Evaluation: Check covariate magnitude to historical deforestation
  3. Model Development: Empirically weighted risk model generation
  4. Deforestation Allocation: Spatial distribution of baseline deforestation rates based on risk surfaces
# Import & reproject to OSM grid systems
aoi_country = geodata::gadm(country="LBR", level=0, path="./assets/AOI") |> sf::st_as_sf() |> sf::st_transform("EPSG:3857")
aoi_states  = geodata::gadm(country="LBR", level=1, path="./assets/AOI") |> sf::st_as_sf() |>
  dplyr::filter(NAME_1 == "Bong") |> sf::st_transform("EPSG:3857")
crs_master = sf::st_crs(aoi_country)
  
# Total extent
area_km2 = as.numeric(sf::st_area(aoi_states)) / 1e6
cat("Study Area: Bong County, Liberia\n")
NA Study Area: Bong County, Liberia
cat("Study Area:", round(area_km2, 0), "km²\n")
NA Study Area: 8565 km²
# Visualize AOI
tmap::tmap_mode("view")
tmap::tm_shape(aoi_country) + tmap::tm_borders(col="purple", lwd=2) +
  tmap::tm_shape(aoi_states) + tmap::tm_borders(col="red", lwd=1) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_title("AOI Jurisdictional Boundaries", size=.8) + 
  tmap::tm_basemap("OpenTopoMap")

Area check

In Liberia, the official definition of forest land is provided by the Forestry Development Authority (Government of Liberia 2019), including areas of land that meet the following criteria:

  • Canopy cover of minimum 30%;
  • Canopy height of minimum 5m or the capacity to reach it;
  • Covering a minimum of 1 hectare of land.
forest_2024=terra::rast("./assets/LULC/outputs/forest_2024.tif")
forest_2024_poly=terra::as.polygons(forest_2024) |> sf::st_as_sf()
forest_2024_poly$area_ha = round(as.numeric(sf::st_area(forest_2024_poly) * 0.0001, 4))
forest_2024_poly |> sf::st_drop_geometry() |> janitor::adorn_totals()
label area_ha
0 13648
1 509738
Total 523386 
slivers = forest_2024_poly |> dplyr::filter(as.numeric(area_ha) < 1)  
slivers
label geometry area_ha 
# no artefacts found??

Spatial Covariates

The following spatial covariates were processed as potential drivers of deforestation risk. Covariates were merged between sociodemographic and geographic datasets surrounding the project area and national level datasets beyond the project area in order to enable jurisdictions analysis.

# Forest edge
forest_2024 = terra::rast("./assets/LULC/outputs/forest_2024.tif")
forest <- terra::classify(forest_2024, cbind(c(1,0), c(1,NA)))
forest_edge = terra::boundaries(forest, inner=FALSE)
distance_forest_edge = terra::distance(forest_edge)
distance_forest_edge = distance_forest_edge / 1000
names(distance_forest_edge) = "edge_distance_km"
writeRaster(distance_forest_edge, "./assets/RISK/distance_to_edge.tif", overwrite=T)

# Roads
osmdata::available_features()  
osmdata::available_tags("route") # "highway", "public_transport", "side_road", "tracks"
osm_bbox = osmdata::getbb("Liberia")
roads_osm = osm_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "route") |>
  osmdata::osmdata_sf()
template <- forest
roads_rast <- rasterize(vect(roads_osm), template, field=1, background=0)
roads_rast <- subst(roads_rast, from=0, to=NA)
distance_to_roads <- distance(roads_rast)
writeRaster(distance_to_roads, "./assets/RISK/distance_to_roads.tif", overwrite=T)

# Hospitals
buildings_hospital = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "amenity", value = "hospital") |>
  osmdata::osmdata_sf() # |> st_transform(crs_master)

# Housing
buildings_residential = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "building", value = "residential") |>
  osmdata::osmdata_sf() # |> st_transform(crs_master)

# Religion
buildings_workship = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "building", value = "religious") |>
  osmdata::osmdata_sf() 

# Towns 
townnames = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "place", value = "town") |>
  osmdata::osmdata_sf()

# Villages
villages = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "place", value = "village") |>
  osmdata::osmdata_sf()

# Places
places_merged = townnames$osm_points |> bind_rows(villages$osm_points) |>
  group_by(across(-geometry)) |> sf::st_cast("POINT") |>
  summarise(geometry = st_union(geometry), .groups = "drop")
sf::st_write(places_merged, "./assets/RISK/distance_to_places.shp", delete_dsn=T) 

places = buildings_hospital$osm_points |>
  bind_rows(buildings_residential$osm_points, buildings_workship$osm_points, places_merged$osm_points) |>
  group_by(across(-geometry)) |> sf::st_cast("POINT") |>
  summarise(geometry = st_union(geometry), .groups = "drop")
sf::st_write(places, "./assets/RISK/distance_to_places.shp", delete_dsn=T) 

waterways_rast<-rasterize(vect(waterways_rast), template, field=1, background=0)
waterways_rast<-subst(waterways_rast, from=0, to=NA)
distance_to_waterways<-distance(waterways_rast)
writeRaster(distance_to_waterways, "./assets/RISK/distance_to_waterways.tif", overwrite=T)

# Normalize covariates 
normalize <- function(x){
  (x-global(x,"min",na.rm=T))/
    (global(x,"max",na.rm=T) -
       global(x, "min",na.rm=T))}

dem   = raster::subset(STACK, "DEM")
slope_tangent = raster::terrain(
  dem, opt="slope",unit="tangent",neighbors=8,filename="./assets/DEM/slope_tangent.tif")
slope_tangent = terra::rast("./assets/DEM/slope_tangent.tif")
slope_percent = slope_tangent * 100
slope_percent = terra::clamp(slope_percent, 0, 100) 
slope_percent = raster::raster(slope_percent)
slope   = raster::raster(slope)

slope               = normalize(slope)
distance_to_forest  = normalize_function(distance_to_forest)
distance_to_roads   = normalize_function(distance_to_roads)
distance_to_places  = normalize_function(distance_to_places)
distance_to_water   = normalize_function(distance_to_waterways)

# invert all risk-producing covariates
distance_to_edge_inv   <- 1 - distance_to_forest
distance_to_roads_inv  <- 1 - distance_to_roads
distance_to_places_inv <- 1 - distance_to_places
distance_to_water_inv  <- 1 - distance_to_waterways

writeRaster(distance_to_edge_inv, "./assets/RISK/distance_to_edge_inverted.tif", overwrite=T)
writeRaster(distance_to_roads_inv, "./assets/RISK/distance_to_roads_inverted.tif", overwrite=T)
writeRaster(distance_to_places_inv, "./assets/RISK/distance_to_places_inverted.tif", overwrite=T)
writeRaster(distance_to_water_inv, "./assets/RISK/distance_to_water_inverted.tif", overwrite=T)
writeRaster(slope_norm, "./assets/RISK/slope_norm.tif", overwrite=T)

Deforestation risk

Two methods were explored for weighting variables and creating a generalized deforestation risk index. We could consider developing a spatial risk model using the spatstat package or logistic regression, as has been cited in recent Verra guides. In addition, some of the heavy lifting with input formatting and data wrangling has already been completed.

However, spatial modelling has tended to produce challenges when fitting such large country-wide covariates. Moreover, these kinds of spatialy driven models tend to require longer training procedures due to their intercept-based spatial kernels and slower resampling patterns.

Alternatively, we have drafted a tentative risk indexing approach based on a weighted sum of subjectively scored covariate effects. While each variable would still need a carefully reasoned score, this option offers a more streamlined method that is easier to adjust. We applied this risk index to inform a risk weighted allocation of the 10-year deforestation rate, first by multiplying the fraction of pixel risk by zonal forest loss, and second by factoring out annual zonal loss by multiplying by pixel risk values, as shown below

\[ \mathrm{AllocatedLoss}_{\mathrm{pixel}} = \left( \frac{\mathrm{risk}_{\mathrm{pixel}}}{\sum \mathrm{risk}_{\mathrm{zone}}} \right) \times \mathrm{annual\_loss\_10yr}_{\mathrm{zone}} \]

\[ \mathrm{allocated\_loss}_{\mathrm{pixel}} = \mathrm{risk}_{\mathrm{pixel}} \times \left( \frac{\mathrm{annual\_loss\_10yr}_{\mathrm{zone}}}{\sum \mathrm{risk}_{\mathrm{zone}}} \right) \]

Both formulas describe the same operation in different orders of multiplication: each pixel in a given zone Z receives a share of annual_loss_10yrZ based on its proportional risk (the pixel’s risk relative to the sum of all pixel risks in that zone). This ensures that higher-risk pixels are allocated more deforestation, in line with the Verra guidance for an allocated deforestation risk map.

We intend to present both of these approaches for broader review and discussion in our upcoming meeting.

risk_index <- (0.2 * distance_to_edge) +
              (0.2 * distance_to_roads) +
              (0.2 * distance_to_places) +
              (0.1 * distance_to_urban) +
              (0.1 * distance_to_water) +
              (0.1 * slope)

# Re-normalize the index to between 0 and 1
rmin <- global(risk_index, "min", na.rm=TRUE)[1]
rmax <- global(risk_index, "max", na.rm=TRUE)[1]
risk_index_norm <- (risk_index - rmin) / (rmax - rmin)
writeRaster(risk_index_norm, "./assets/RISK/deforestation_risk_index.tif",overwrite=T)

# Returns a data.frame with columns: ID, risk_index_norm_sum
risk_sum <- extract(risk_index_norm, zones_sv, fun = sum, na.rm=TRUE)
colnames(risk_sum) <- c("zone_id","sum_risk")
zones_stats <- merge(zones_sf, risk_sum, by="zone_id", all.x=TRUE)
zones_rast <- rasterize(vect(zones_stats), risk_index_norm, field="zone_id")
zones_stats$loss_factor <- zones_stats$annual_loss_10yr / zones_stats$sum_risk
loss_factor_rast <- rasterize(vect(zones_stats), risk_index_norm, field="loss_factor")
allocated_loss <- risk_index_norm * loss_factor_rast
writeRaster(allocated_loss, "./assets/RISK/allocated_deforestation.tif", overwrite=T)

# Visualize
terra::plot(risk_index_norm, main="Deforestation risk map")
terra::plot(allocated_loss, main="Allocated deforestation map")

Runtime log
devtools::session_info()

References

Government of Liberia. 2019. “Liberia’s Forest Reference Emission Level Submission to the UNFCCC.”