In albhasan/sits.prodes: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

library(dplyr)
library(ggplot2)
library(kableExtra)
library(knitr)
library(magrittr)
library(purrr)
library(tidyr)
library(sits.prodes)
library(xtable)

utils::data("prodes_mapbiomas", package = "sits.prodes")
utils::data("prodes_labels", package = "sits.prodes")
utils::data("mapbiomas_labels", package = "sits.prodes")

Mapbiomas

Their Landsat images' maximum cloud cover is 30%.
Their classification algorithm is random forest.
Their maps accuracy is unavailable (under construction).
Their maps have 6 (plus 'No information') classes (out of 27 for the whole Mapbiomas).
They use 10 variables for classification:
- GVS Green Vegetation Shade (mean and sd).
- NFVI Normalized Fraction Difference Index (mean and sd).
- NPV Non-Photosynthetic Vegetation (mean and sd).
- Shade and Soil fractions (means and sds).
During training, they sampled TerraClass, GlobeLand30, PRODES, Global Forest Cover, ESA CCI (among others) using a majority rule.
They used postprocessing rules for improving their results.

Label mappings

PRODES - Portuguese to English

prodes_labels %>% 
    dplyr::select(-id_pd) %>% 
    knitr::kable() %>% 
    kableExtra::kable_styling(full_width = F)

Mapbiomas to PRODES

mapbiomas_labels %>% 
    dplyr::filter(id_mb %in% c(3, 13, 21, 25, 26, 27)) %>% 
    dplyr::select(-id_mb) %>%
    knitr::kable() %>% 
    kableExtra::kable_styling(full_width = F)

Results

PRODES is the reference map (it is represented in the columns)
Data given in numer of pixels (except for user and producer accuracies)

# sort by scene-year
prodes_mapbiomas <- prodes_mapbiomas %>% 
  dplyr::arrange(tile, year_pd)

# extract the confusion matrices
prodes_mapbiomas$cm <- lapply(prodes_mapbiomas$conmat, function(x){
    x %>% .$table %>% as.matrix() %>% return()    
})

# compute user-producer accuracies
prodes_mapbiomas$up_acc <- add_upacc(prodes_mapbiomas$cm)

# get data for plots
acc_tb <- prodes_mapbiomas %>% 
    dplyr::mutate(up_acc2 = purrr::map(.$conmat, function(x){
        res <- sits.prodes::asses_accuracy_simple(x$table)
        c(res$user, res$producer) %>%
            tibble::enframe(name = NULL) %>%
            dplyr::bind_cols(cname =  c(paste0("ua_", names(res$user)), 
                                        paste0("pa_", names(res$producer)))) %>%
            tidyr::spread("cname", "value") %>%
            return()
    })) %>% 
    dplyr::select(tile, year_pd, up_acc2) %>% 
    tidyr::unnest() %>% 
    tidyr::gather(tidyselect::starts_with("pa_"), tidyselect::starts_with("ua_"), 
                  key = "variable", value = "value") %>% 
    dplyr::mutate(acc_type = ifelse(stringr::str_sub(variable, 1, 3) == "ua_", 
                                    "user", "producer"),
                  variable = stringr::str_sub(variable, 4)) 
acc_pr <- acc_tb %>% dplyr::filter(acc_type == "producer") %>% 
    dplyr::rename(prod_acc = value) %>% 
    dplyr::select(-acc_type)
acc_ur <- acc_tb %>% dplyr::filter(acc_type == "user") %>% 
    dplyr::rename(user_acc = value) %>% 
    dplyr::select(-acc_type)
pa_acc_tb <- acc_pr %>% 
    dplyr::inner_join(acc_ur, by = c("tile", "year_pd", "variable"))
prodes_years <- pa_acc_tb %>% dplyr::pull(year_pd) %>% unique()
pa_acc_tb %>% 
  ggplot2::ggplot(aes(x = prod_acc, y = user_acc, 
                      shape = variable, color = variable)) +
  ggplot2::geom_point(size = 3) +
  ggplot2::coord_fixed() +
  ggplot2::xlim(0, 1) +
  ggplot2::ylim(0, 1) +
  ggplot2::facet_wrap(tile ~ year_pd, ncol = length(prodes_years)) + 
  ggplot2::guides(shape = guide_legend(override.aes = list(size = 4))) +
  ggplot2::labs(title = "PRODES - Mapbiomas", 
                x = "Producer accuracy", 
                y = "User accuracy") +
  ggplot2::theme(axis.text.x = element_text(angle = 90, hjust = 1))

for (i in 1:nrow(prodes_mapbiomas)) {
    cap <- sprintf("Confusion matrix %s %s", unlist(prodes_mapbiomas[i, "tile"]), 
                   unlist(prodes_mapbiomas[i, "year_pd"]))
    cat(paste0("### ", cap), "\n")
    p_conmat <- prodes_mapbiomas$up_acc[[i]] %>% 
        apply(2, function(x){x/1}) %>% 
        tibble::as_tibble(rownames = NA, caption = cap) %>% 
        knitr::kable(digits = 2, row.names = TRUE) %>% 
        kableExtra::kable_styling(full_width = F)
    print(p_conmat)
    cat("\n\n")
}