knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
suppressMessages(library(kableExtra))
suppressMessages(library(knitr))
suppressMessages(library(magrittr))
suppressMessages(library(purrr))
suppressMessages(library(tidyr))
suppressMessages(library(sf))
suppressMessages(library(sits.prodes))
suppressMessages(library(xtable))

utils::data("expert_validation", package = "sits.prodes")

# Get the accuracies out of a confusionMatrix object. 
# 
# @param  x A confusionMatrix object (see caret::confusionMatrix).
# @return A list of numeric. The overall, user, and producer accuracies.
.get_accuracies <- function(x){
    if (length(x) == 1 && is.na(x)) 
        return(tibble::tibble(NA))
    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()
}


# Check if the confusion matrix's values are greater than zero. 
# 
# @param  x         A confusionMatrix object (see caret::confusionMatrix).
# @param  label_vec A chatacter. Names of the labels to consider.
# @return           A length-one logical.
.is_mt_complete <- function(x, label_vec){
    if (length(x) == 1 && is.na(x)) 
        return(FALSE)
    mt <- x$table %>% as.matrix()
    if (!all(label_vec %in% colnames(mt)) && all(label_vec %in% rownames(mt)))
        return(FALSE)
    mt_lab <- mt[label_vec, label_vec]
    if (any(is.na(mt_lab)))
        return(FALSE)
    if (any(apply(mt_lab, 2, sum) == 0))
        return(FALSE)
    return(TRUE)
}

PRODES mapping to English

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

Results

# get data for plots
acc_tb <- expert_validation %>%
    dplyr::mutate(
        is_mt_complete = purrr::map_lgl(.$confusion_matrix, .is_mt_complete, 
                                    label_vec = c("deforestation", "forest")),
        up_acc2 = purrr::map(.$confusion_matrix, .get_accuracies)
    ) %>%
    dplyr::filter(is_mt_complete == TRUE) %>%
    tidyr::drop_na(confusion_matrix) %>% 
    dplyr::select(experiment, algorithm, smooth, scene, pyear, up_acc2) %>%
    tidyr::unnest() %>%
    dplyr::select(-tidyselect::ends_with(match = "no forest"), 
                  -tidyselect::ends_with(match = "water")) %>%
    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("experiment", "algorithm", "smooth", 
                                     "scene", "pyear", "variable"))
f_plot <- function(data_tb, my_algorithm, my_smooth){
  prodes_years <- data_tb %>% 
    dplyr::pull(pyear) %>% 
    unique() %>% 
    sort() 
  data_tb <- data_tb %>%
    dplyr::filter(smooth == !!my_smooth, 
                  algorithm == !!my_algorithm)
  if (nrow(data_tb) < 1 || ncol(data_tb) < 1) {
    warning(sprintf("No data found for  %s %s", my_algorithm, my_smooth))
    return()
  }
  plot(
    data_tb %>% 
      ggplot2::ggplot(aes(x = prod_acc, y = user_acc,
                          shape = variable, color = experiment)) +
      ggplot2::geom_point(size = 2) +
      ggplot2::coord_fixed() +
      ggplot2::xlim(0, 1) +
      ggplot2::ylim(0, 1) +
      ggplot2::guides(shape = guide_legend(override.aes = list(size = 3))) +
      ggplot2::facet_wrap(scene ~ pyear, ncol = length(prodes_years)) +
      ggplot2::labs(title = sprintf("Expert validation %s %s", 
                                    my_algorithm, my_smooth),
                    x = "Producer accuracy",
                    y = "User accuracy")
  )
}

for (my_algorithm in unique(pa_acc_tb$algorithm)) {
    cat("\n\n")
    cat(paste0("#### ", my_algorithm), "\n")
    for (my_smooth in unique(pa_acc_tb$smooth)) {
        print(paste(my_algorithm, my_smooth, sep = "-"))
      #-----
      #### results_vote 
#[1] "results_vote-smooth_3x3_n10"
#[1] "results_vote-smooth_5x5_n10"
#[1] "results_vote-smooth_7x7_n10"
#[1] "results_vote-smooth_9x9_n10"
#[1] "results_vote-smooth_11x11_n10"
#[1] "results_vote-"
        suppressWarnings(f_plot(data_tb = pa_acc_tb, my_algorithm, my_smooth))
    }
}
for (i in 1:nrow(expert_validation)) {
    cap <- sprintf("Confusion matrix %s %s", unlist(expert_validation[i, "scene"]), 
                   unlist(expert_validation[i, "pyear"]))
    cat(paste0("### ", cap), "\n")
    if (!(length(expert_validation$confusion_matrix[[i]]) == 1 &&
          is.na(expert_validation$confusion_matrix[[i]]))) {
        p_conmat <- expert_validation$confusion_matrix[[i]]$table %>%
            apply(2, function(x){x/1}) %>% 
            add_upacc() %>%
            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")
}  


albhasan/sits.prodes documentation built on Dec. 29, 2019, 5:03 a.m.