16 Deforestation Emissions
16.1 Land Cover Classification
This target extracts the annual deforested extents from land-cover classification (LCC) maps which were generated in Google Earth Engine from a combination of Landsat 8 and 9 (L8 & L9), sentinel 2 (S2) and Harmonized Landsat-Sentinel (HLS) data. For each year, either L8, L9, S2 or HLS were used to build annual composites; data were not mixed within years. Annual median composites were generated for years 2015-2024 from cloud masked images. For 2016,it was necessary to split the region into two areas to ensure that a known region of deforestation was correctly captured before green-up. HLS data were downloaded and pre-processed for years 2023 and 2024 using the open data cube python package ecosystem, saved and then ingested into Google Earth Engine; processing workflow is available from here. A random forest model was trained for each year individually using ~100 labelled training polygons which were digitised by an analyst through comparison with both the annual composite and high resolution Planet labs data. labelled points were assigned a category of either: forest, non-forest, water, or cloud. Classification accuracy for all years was greater than 90% when evaluated with a 70-30% train-test split, which was derived from random sampling.
16.2 Deforestation Detection
Areas of deforestation were extracted from these land classification maps using a cumulative differencing approach where intersecting regions of non-forest in consecutive years were considered stable and not deforested. Novel regions of non-forest land cover, that were not present in the preceding year, were classified as deforested. A sieve filter was used to remove areas of potential deforestation with an area < 2.5 ha. This approach was used as it was considered that, at the spatial grain of the satellite imagery that was used (10-30 m), it was not possible to differentiate between deforestation and forest degradation for areas < 2.5 ha. However, any degradation occurring in these regions is also captured in the biomass modelling, presented in Chapter 7, and will therefore be reflected in the estimated ongoing forest growth rate.
16.3 Manual Digitisation of Deforestation
In one isolated case, in 2023, the random forest LCC was not able to appropriately identify a known region of deforestation that occurred as a result of flood inundation. This omission was most likely due to the unique spectral signature of the event (stagnant flood waters) being absent from the training data. The deforestation was identified during routine monitoring of high resolution Planet Labs imagery - a field inspection was carried out to confirm the deforestation. The deforestation was manually digitised from high resolution imagery and merged with the LCC-derived deforestation dataset.
Figure 16.1 presents a map of annual deforestation area across the project area. Note that much of the illegal deforestation has occurred along main rivers and access routes in the central region of the project. No deforestation with an area >= 2.5 ha has been detected in the project area during 2021 and 2022.
16.4 Deforestation Emissions
The net greenhouse gas emissions due to deforestation within the project area were calculated in line with Equation 35 of the VM0010 methods as shown in Equation 16.1
\[ \Delta C_{DIST \space t \space PRJ} = C_{AB \space i \space t \space PRJ} \times \frac{44}{12} \tag{16.1}\]
Where \(\Delta C_{DIST \space t \space PRJ}\) is the net greenhouse gas emissions resulting from disturbance (including illegal logging, fire and natural disturbance) each year, \(C_{AB \space i \space t \space PRJ}\) is total the carbon stock in aboveground biomass in the deforested area for event \(i\) in the year \(t\), and \(\frac{44}{12}\) is the ratio of molecular weights of carbon dioxide and carbon.
Unlike Equation 35 in the VCS methodology, it is not necessary to split this calculation further by strata as our biomass estimates are derived directly from the biomass maps (Chapter 7). Our estimates of the Above Ground carbon emissions are based on the carbon stock in the year of the event, rather than the baseline year, and therefore our approach captures growth since the start of the project, in addition to the baseline stocks.
The annual calculations for net greenhouse gas emissions are presented in the following table
16.5 Earth engine app
To view these data on an interactive google earth engine app follow this link
16.6 Code
Deforested Area Extraction Code R/deforest_extract.R
#' helper for getting site specific LULC maps.
#' to be used for years 2016-2022
#' @param yr integer - year of landclass map.
#' @return character - url to the landclass map.
kuamut_kdefor_url_mr1 <- function(yr) {
sprintf("https://storage.googleapis.com/rgee_store/Kuamut-Defor-L8-L9-S2/Kuamut_LC_%s.tif", yr)
}
#' helper for getting site specific LULC maps.
#' to be used for years 2023-2024
#' @param yr integer - year of landclass map.
#' @return character - url to the landclass map.
kuamut_kdefor_url_mr2 <- function(yr) {
sprintf("https://storage.googleapis.com/rgee_store/Kuamut-HSL-Defor/Kuamut_HSL_LC_%s.tif", yr)
}
#' Function to extract deforestation areas from Land Use Land Cover (LULC) maps.
#'
#' This function extracts deforestation areas from the LULC maps. These
#' maps were generated using a Random Forest classifier trained on a combination
#' of Landsat 8, Sentinel 2 and Harmonized Sentinel-Landsat (HLS) annual
#' composites.
#'
#' @param rds_path Path to the roads detail file.
#' @param yr_range Range of years to extract deforestation from.
#' @param area.thresh Area of contiguous non-forest since previous year
#' considered to be deforestation.
#' @param sub_fold Subfolder to save the output file.
#' @param outfile_version Version of the output file.
#'
#' @return A vector spatial file of deforested areas by year.
deforest_extract <- function(
rds_path,
yr_range = 2016:2024,
area.thresh = 2.5,
sub_fold = "data_out/RoutputProducts",
outfile_version = "v0_5") {
output_dir <- file.path(sub_fold, "Deforest_Vector")
terra_filter_patches <- function(r, .thresh) {
y <- terra::patches(r, directions = 8L)
# remove patches smaller than x ha
rz <- terra::zonal(
terra::cellSize(y, unit = "ha"),
y, sum,
as.raster = TRUE
)
terra::ifel(rz < .thresh, NA, y)
}
get_non_forest <- function(yr,
kuam_rds_path = rds_path) {
if (yr == 2016) {
kdf <-
ezwarp(
list(kuamut_kdefor_url_mr1(2015), kuamut_kdefor_url_mr1(2016)),
kuamut_kdefor_url_mr1(2016),
resample = "near"
) |>
project("EPSG:32650")
} else if (yr %in% 2017:2022) {
kdf <-
ezwarp(kuamut_kdefor_url_mr1(yr),
kuamut_kdefor_url_mr1(2016),
resample = "near"
) |>
project("EPSG:32650")
} else if (yr %in% 2023:2024) {
kdf <-
ezwarp(kuamut_kdefor_url_mr2(yr),
kuamut_kdefor_url_mr1(2016),
resample = "near"
) |>
project("EPSG:32650")
}
if (yr == 2016) {
k.rds <- vect(kuam_rds_path)
rds <- rasterize(k.rds, kdf, touches = TRUE)
kdf[rds == 1] <- 2
}
kdf[kdf != 2] <- NA
return(kdf)
}
insist_get_non_forest <- purrr::insistently(get_non_forest,
rate = purrr::rate_backoff(max_times = 6),
quiet = FALSE
)
get_non_forest_all <- function(yr_range) {
it.pr <- cli::cli_progress_along(yr_range,
name = "non-forest extract"
)
x <- purrr::map(
it.pr,
~ insist_get_non_forest(yr_range[.x])
) |>
terra::rast()
names(x) <- paste0("YEAR_", yr_range)
return(x)
}
calc_deforest_rast <- function(nf_rast, .thresh = 1, years) {
make_yr_names <- function(yr) {
paste(yr, yr + 1, sep = "_")
}
yr_names <- years[seq_along(years) - 1] |>
purrr::map_chr(~ make_yr_names(.x))
accum_deforest <- nf_rast[[1]]
deforested <- terra::rast(nf_rast[[1]], vals = NA, nlyrs = length(yr_names))
for (b in 2:dim(nf_rast)[3]) {
deforested[[b - 1]][is.na(accum_deforest) & !is.na(nf_rast[[b]])] <- 1
deforested[[b - 1]] <- terra_filter_patches(deforested[[b - 1]],
.thresh = .thresh
)
all_deforest <- sum(deforested, na.rm = TRUE)
cum_nf <- sum(nf_rast[[1:b]], na.rm = TRUE)
accum_deforest <- sum(c(accum_deforest, cum_nf), na.rm = TRUE)
# Run Garbage collection at the end of each run.
gc()
}
names(deforested) <- yr_names
stack_to_sf <- function(deforested_raster) {
polys <- terra::as.polygons(deforested_raster) |>
sf::st_as_sf()
if (nrow(polys) == 0) {
return(NULL)
}
polys |>
dplyr::rename(PATCH = 1) |>
mutate(PERIOD = names(deforested_raster))
}
cdf_v <- seq_along((deforested)[3]) |>
purrr::map_dfr(~ stack_to_sf(deforested[[.x]])) |>
group_by(PERIOD) |>
dplyr::summarise()
.areas <- units::set_units(sf::st_area(cdf_v), "ha")
cdf_v |>
mutate(DEFOREST_AREA = .areas)
}
fwrap <- function(nf_ras, thresh, .v) {
x <- nf_ras |>
calc_deforest_rast(.thresh = thresh, years = yr_range)
out.file <- set_path(sprintf("kuamut_deforest_%s.fgb", .v), output_dir)
sf::write_sf(x, out.file,
delete_dsn = TRUE
)
return(out.file)
}
nf <- get_non_forest_all(yr_range)
fwrap(nf, area.thresh, .v = outfile_version)
}
#' function to merge auto and manual deforestation data
#' @param auto_def a file path for the spatial vector file comprising
#' annual deforestation areas, detected from the random forest LULC map. The
#' output of `deforest_extract`.
#' @param ... file path(s) for the spatial vector files comprising manually
#' digitised deforestation areas. must contain a column for the year of the
#' event named in `year_col`.
#' @param year_col a character string for the column name containing the year
#' of the manual deforestation event. either length 1 if all manual files are
#' have the same naming (preferred) or the same length as the number of manual
#' files.
#' @param sub_fold a character string for the subfolder to save the output file.
#' @param outfile_version a character string for the version of the output file.
#' @return a file path to the output file.
#' @details
#' This function is required where the LULC classification doesn't adequately
#' capture specfic deforestation events. Such cases include localised flooding.
#'
rbind_auto_manual_deforest <- function(
auto_def, ...,
year_col = "year",
sub_fold = "data_out/RoutputProducts",
outfile_version = "v0_5") {
output_dir <- file.path(sub_fold, "Deforest_Vector")
ad <- read_sf(auto_def) |>
mutate(DETECT_TYPE = "auto")
mdlist <- rlang::dots_list(...)
purrr::map(mdlist, assert_file_exists)
if (!length(year_col) %in% c(1, length(mdlist))) {
cli::cli_abort(
c("x" = "The length of `year_col` must be 1 or equal to the number
of manual deforestation files.")
)
}
read_man_def <- function(x, yc) {
yc <- ensym(yc)
man_def <- read_sf(x)
man_def <- man_def |>
mutate(
PERIOD = paste(!!yc, !!yc + 1, sep = "_"),
DEFOREST_AREA = round(
as.numeric(
units::set_units(sf::st_area(man_def), "ha")
), 1
),
DETECT_TYPE = "manual"
) |>
dplyr::select(PERIOD, DEFOREST_AREA, DETECT_TYPE)
}
md <- purrr::map2(mdlist, year_col, read_man_def) |>
purrr::list_rbind()
auto_man_defor <- rbind(ad, md) |>
sf::st_cast("MULTIPOLYGON")
out_file <- set_path(
sprintf("kuamut_deforest_combined_%s.fgb", outfile_version), output_dir
)
sf::write_sf(auto_man_defor, out_file,
delete_dsn = TRUE
)
return(out_file)
}Deforestation Emission Calculations Code R/deforest_emissions.R
#' Calculate the total C lost each year due to deforestation.
#'
#' @param acd_rasters list of ACD raster paths
#' @param deforest_vect spatial vector defining the extent of deforestation
#' each year
#'
#' @return A list with summary tibble and exported table.
deforest_emissions <- function(acd_rasters, deforest_vect) {
deforest_vect <- deforest_vect |>
sf::read_sf()
# get raster years
yr_names <- purrr::imap_dbl(acd_rasters, ~ readr::parse_number(.y)) |>
unname()
# make temporal raster stack
acd_stack <- acd_rasters |>
purrr::map(function(x) {
rast(x)[[1]]
}) |>
rast()
# function to calculate carbon emissions due to deforestation.
defor_loss <- function(s.yr) {
# select the corresponding raster
s.band <- acd_stack[[paste0("ACD_map_", s.yr)]]
# select the correct year from def. area
vect_name <- paste(s.yr, s.yr + 1, sep = "_")
yr_area <- deforest_vect[deforest_vect$PERIOD == vect_name, ]
# get deforest area
ext_area <- sf::st_area(yr_area) |>
units::set_units("ha") |>
as.numeric()
# get raster resolution for total C conversion
.res <- terra::res(s.band)[1]
# mask the raster with the def. area
defor_only <- terra::mask(s.band, terra::vect(yr_area))
# calculate total C lost. This is fine - exactextract not needed because the
# polygon is derived from the raster anyway so makes no difference
totC <- sum(defor_only[], na.rm = TRUE) * ((.res * .res) / 10000)
# convert deforest_emission_calcs to the required units for annual calcs.
CO2_frac <- totC * (44 / 12)
# build df.
tibble(
Period = s.yr,
`Deforest Area` = ext_area,
`Total C lost` = totC,
`C_DIST_IL_i_t_PRJ` = CO2_frac
) |>
mutate(across(where(is.numeric), \(x) round(x, digits = 1)))
}
defor_emit_df <- purrr::map_df(yr_names, defor_loss)
export_tab_n_csv(
defor_emit_df, "Deforest_Emissions",
"data_out/Routputs/Tables"
)
return(defor_emit_df)
}