Skip to content

Code Repository for Bootstrap Assessment of Crop Area Estimates using Satellite Pixels Counting

License

Notifications You must be signed in to change notification settings

castlaboratory/croparea

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

11 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Bootstrap Assessment of Crop Area Estimates using Satellite Pixels Counting

André Leite, Cristiano Ferraz, Jacques Delincé & Raydonal Ospina

2022-04-15

Packages

library(pracma)
library(tidyverse)
library(readxl)
library(glmnet)
library(knitr)
library(kableExtra)
library(plotly)
library(paletteer)

Function to generate artificial population

gen_U_star <- function(a_matrix, size, strategy){
  stopifnot(strategy %in% c('one', 'two', 'three'))
  if (strategy == 'one') {
    U_star <- 
    expand_grid(
      c = colnames(a_matrix), 
      g = rownames(a_matrix)) %>% 
    mutate(n = c(a_matrix), 
           pi = n/sum(n), 
           Code = paste0(c, ' - ', g)) %>% 
    sample_n(size = size, replace = TRUE, weight = pi) %>% 
    select(g, c, pi)
  } else if (strategy == 'two') {
    sample_size_group <- size/nrow(a_matrix)
    U_star <-
      a_matrix %>% 
      as_tibble(rownames = 'g') %>% 
      pivot_longer(-g, names_to = 'c', values_to = 'n') %>% 
      group_by(c) %>% 
      mutate(pi = n/sum(n)) %>% 
      group_split() %>% 
      map_dfr( ~ sample_n(.x, size = sample_size_group, replace = TRUE, weight = pi))  %>% 
      select(g, c, pi)
  } else if (strategy == 'three') {
    sample_size_group <- size/ncol(a_matrix)
    U_star <-
      a_matrix %>% 
      as_tibble(rownames = 'g') %>% 
      pivot_longer(-g, names_to = 'c', values_to = 'n') %>% 
      group_by(g) %>% 
      mutate(pi = n/sum(n)) %>% 
      group_split() %>% 
      map_dfr( ~ sample_n(.x, size = sample_size_group, replace = TRUE, weight = pi)) %>% 
      select(g, c, pi)
  }
  return(U_star)
}

Function to get a sample from artificial population

# get_a_sample <- function(U_star, N_0, strategy, order){
#   a_sample <- U_star %>% 
#     sample_n(N_0) %>% 
#     count(c, g) %>% 
#     pivot_wider(names_from = 'c', values_from = n, values_fill = 0) %>% 
#     column_to_rownames(var = 'g') %>% 
#     as.matrix()  
#   if (is.null(order)) return(a_sample)
#   else return(a_sample[order, order])
# }

### Three strategies function

get_a_sample_3st <- function(U_star, N_0, strategy = 'one', order){
  stopifnot(strategy %in% c('one', 'two', 'three'))
  if (strategy == 'one') {
    a_sample <- U_star %>% 
      sample_n(N_0) %>% 
      count(c, g) %>% 
      pivot_wider(names_from = 'c', values_from = n, values_fill = 0) %>% 
      arrange(match(g, order)) %>% 
      select(crop = g, all_of(order)) %>% 
      column_to_rownames(var = 'crop') %>% 
      as.matrix() 
  } else if (strategy == 'two') {
    sample_size_group <- N_0/nrow(a_0)
    a_sample <- 
      U_star %>% 
      count(g, c) %>% 
      group_by(c) %>% 
      select(c, g, n) %>% 
      mutate(pi = n/sum(n)) %>% 
      group_split() %>% 
      map_dfr( ~ sample_n(.x, size = sample_size_group, replace = TRUE, weight = pi)) %>% 
      count(g, c) %>% 
      pivot_wider(names_from = "c", values_from = "n", values_fill = 0) %>% 
      arrange(match(g, order)) %>% 
      select(crop = g, all_of(order)) %>% 
      column_to_rownames(var = 'crop') %>% 
      as.matrix()
  } else if (strategy == 'three') {
    sample_size_group <- N_0/ncol(a_0)
    a_sample <- 
      U_star %>% 
      count(g, c) %>% 
      group_by(g) %>% 
      select(g, c, n) %>% 
      mutate(pi = n/sum(n)) %>% 
      group_split() %>% 
      map_dfr( ~ sample_n(.x, size =  sample_size_group, replace = TRUE, weight = pi)) %>% 
      count(g, c) %>% 
      pivot_wider(names_from = "c", values_from = "n", values_fill = 0) %>% 
      arrange(match(g, order)) %>% 
      select(crop = g, all_of(order)) %>% 
      column_to_rownames(var = 'crop') %>% 
      as.matrix()
  }
  return(a_sample)
}

Function to get estimates from sample confution matrix and artificial population

get_estimates <- function(a_matrix, R_U_0){
  e_g_dado_c <- a_matrix %*% solve(diag(colSums(a_matrix)))
  e_c_dado_g <- t(a_matrix) %*% solve(diag(rowSums(a_matrix)))
  inverse_e_c_dado_g <- solve(e_c_dado_g)
  rownames(inverse_e_c_dado_g) <- rownames(e_c_dado_g)
  
  direct <- (e_g_dado_c %*% R_U_0) %>% 
    as_tibble(rownames = 'crop') %>% 
    mutate(estimator = 'direct')
  
  inverse <- (inverse_e_c_dado_g %*% R_U_0) %>%
    as_tibble(rownames = 'crop') %>% 
    mutate(estimator = 'inverse') 
  
  return(bind_rows(direct, inverse))
}

Function to get bootstrap replicas

get_bootstrap <- function(B, U_star, N_0, R_U_0, strategy, order) {
  1:B %>% 
    map(~ get_a_sample_3st(U_star, N_0, strategy, order)) %>% 
    map_dfr(~ get_estimates(.x, R_U_0)) %>% 
    mutate(replica = ceiling(row_number()/6))
}

Simulated data

jacques_tibble <- read_excel("data/Simdata.xlsx", 
                      sheet = "Sheet1", 
                      range = "B3:G8",
                      .name_repair = ~ if_else(.x == '', 'crop', .x))


jacques_matrix <- jacques_tibble %>% 
  column_to_rownames(var = "crop") %>% 
  as.matrix() 
jacques_matrix %>% 
  as.data.frame() %>% 
  rownames_to_column(var = "crop") %>% 
  kable(caption = "Jacques' matrix") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE)
Jacques’ matrix
crop Wheat Rapeseed Corn Sugarbeet Others
Wheat 0.200 0.020 0.005 0.005 0.020
Rapeseed 0.010 0.030 0.000 0.000 0.010
Corn 0.001 0.010 0.080 0.005 0.004
Sugarbeet 0.005 0.015 0.040 0.120 0.020
Others 0.100 0.020 0.010 0.030 0.240

T and V vectors

Size <- 1e6
R_vector <- jacques_matrix %>% 
  { matrix(colSums(.), ncol = 1, dimnames = list(colnames(.), 'value')) * Size } %>%
  t() 
T_vector <- jacques_matrix %>% 
  { matrix(rowSums(.), ncol = 1, dimnames = list(colnames(.), 'value')) * Size } 

Get a artificial population and associated R vector

# Population Size
Size = 1000000
# Seed
set.seed(112358)

U_0 <- gen_U_star(jacques_matrix, strategy = 'one', size = Size)
R_U_0 <- U_0 %>% 
  count(c) %>% 
  arrange(match(c, colnames(jacques_matrix))) %>% 
  column_to_rownames(var = 'c') %>% 
  as.matrix() 

T_U_0 <- U_0 %>% 
  count(g) %>% 
  arrange(match(g, colnames(jacques_matrix))) %>% 
  column_to_rownames(var = 'g') %>% 
  as.matrix() 

Reference <- 
  bind_rows(
    T_U_0 %>% as_tibble(rownames = "crop") %>% mutate(References = "T"),
    R_U_0 %>% as_tibble(rownames = "crop") %>% mutate(References = "R")
  ) %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE))

# a_0 <- U_0 %>% 
#   sample_n(N_0) %>% 
#   count(g, c) %>%   
#   pivot_wider(names_from = "c", values_from = "n",values_fill = 0) %>% 
#   arrange(match(g, colnames(jacques_matrix))) %>% 
#   select(crop = g, all_of(colnames(jacques_matrix))) %>% 
#   column_to_rownames(var = 'crop') %>% 
#   as.matrix() 

T_0 and R_0 tables

T_U_0 %>% 
  kable(caption = "T_0 vector") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) 
T_0 vector
n
Wheat 249615
Rapeseed 50086
Corn 99736
Sugarbeet 200341
Others 400222
R_U_0 %>% 
  kable(caption = "R_0 vector") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) 
R_0 vector
n
Wheat 316083
Rapeseed 94924
Corn 135024
Sugarbeet 160274
Others 293695

Strategy 1: Bivariate

a_0, r_0, t_0

## Sample
N_0 <- 1000
a_0 <- get_a_sample_3st(U_star = U_0, N_0 = 1000, order = colnames(jacques_matrix), strategy = "one")

r_vector <- a_0 %>% 
  { matrix(colSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} %>%
  t() 
t_vector <- a_0 %>% 
  { matrix(rowSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} 

a_0 %>% as_tibble(rownames = 'crop') %>% 
  bind_rows(as_tibble(r_vector) %>% mutate(crop = 'R')) %>% 
  bind_cols(as_tibble(t_vector) %>% rename(T = value) %>% 
              bind_rows(tibble(`T` = 1000))) %>% 
  kable(caption = "a_0 confusion sample matrix") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) %>% 
  column_spec(c(1,7), bold = TRUE) %>% 
  row_spec(6, bold = TRUE)
a_0 confusion sample matrix
crop Wheat Rapeseed Corn Sugarbeet Others T
Wheat 201 23 6 3 19 252
Rapeseed 11 36 0 0 7 54
Corn 4 8 82 6 5 105
Sugarbeet 4 17 38 117 19 195
Others 108 31 7 29 219 394
R 328 115 133 155 269 1000

Get a artificial population for bootstrap sampling

U_star <- gen_U_star(a_0, strategy = 'one', size = Size)

Generate 1000 bootstrap replicas

#set.seed(112358)
result_one <- get_bootstrap(1000, U_star, N_0, R_U_0, strategy = 'one', order = colnames(jacques_matrix)) %>% 
  mutate(Sample = '1000')
write_rds(result_one, "data/result_one_N1m.rds")

Plot result

Palette <- palettes_d_names %>%
  filter(length == 3) %>%
  mutate(code = paste0(package, "::", palette)) %>%
  pull(code) %>%
  .[7] #3 5


result_one %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = 1, data = Reference) + 
  scale_color_paletteer_d(Palette)

result_one %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference) 

Bootstrap summary

table <- result_one %>% group_by(crop, estimator, Sample) %>% 
    summarise(Mean = mean(n), 
              Variance = var(n), 
              `Standard Deviation` = sd(n), 
              CV = `Standard Deviation`/Mean, 
              .groups = 'drop') %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  left_join(Reference %>% pivot_wider(names_from = References,  values_from = "n"), by = "crop") %>% 
  mutate(`T Bias` = `T` - Mean,
         `R Bias` = `R` - Mean) %>% 
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Mean = format(round(Mean), nsmall = 0, big.mark = ","),
         Variance = format(round(Variance), nsmall = 0, big.mark = ","),
         `Standard Deviation` = format(round(`Standard Deviation`), nsmall = 0, big.mark = ","),
         CV = format(round(CV, 4), nsmall = 2, big.mark = ","), 
         `T` = format(round(`T`), nsmall = 0, big.mark = ","),
         `T Bias` = format(round(`T Bias`), nsmall = 0, big.mark = ","),
         `R` = format(round(`R`), nsmall = 0, big.mark = ","),
         `R Bias` = format(round(`R Bias`), nsmall = 0, big.mark = ","))
table %>% 
  kable(caption = "Direct and Inverse Estimators Bootstrap Results (1000 rounds)", align = "lcrrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  column_spec(10, color = if_else(str_detect(table$`T Bias`, "-"), 'red', 'black')) %>%  
  column_spec(11, color = if_else(str_detect(table$`R Bias`, "-"), 'red', 'black')) 
Direct and Inverse Estimators Bootstrap Results (1000 rounds)
Crop Estimator Sample Mean Variance Standard Deviation CV T R T Bias R Bias
Corn direct 1000 105,769 52,343,046 7,235 0.0684 99,736 135,024 -6,033 29,255
Corn inverse 1000 106,371 155,276,956 12,461 0.1171 99,736 135,024 -6,635 28,653
Others direct 1000 405,285 169,370,330 13,014 0.0321 400,222 293,695 -5,063 -111,590
Others inverse 1000 450,186 789,857,498 28,104 0.0624 400,222 293,695 -49,964 -156,491
Rapeseed direct 1000 47,675 33,723,074 5,807 0.1218 50,086 94,924 2,411 47,249
Rapeseed inverse 1000 19,902 245,375,140 15,664 0.7871 50,086 94,924 30,184 75,022
Sugarbeet direct 1000 198,194 90,921,509 9,535 0.0481 200,341 160,274 2,147 -37,920
Sugarbeet inverse 1000 197,975 315,519,157 17,763 0.0897 200,341 160,274 2,366 -37,701
Wheat direct 1000 243,077 109,212,619 10,450 0.0430 249,615 316,083 6,538 73,006
Wheat inverse 1000 225,566 544,071,760 23,325 0.1034 249,615 316,083 24,049 90,517

Other plots

Crop: Wheat

result_one %>% filter(crop == 'Wheat', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) + 
  scale_color_paletteer_d(Palette)

result_one %>% 
  filter(crop == 'Wheat', Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) 

Crop: Rapeseed

result_one %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) + 
  scale_color_paletteer_d(Palette)

result_one %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) 

Crop: Corn

result_one %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Corn") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) + 
  scale_color_paletteer_d(Palette)

result_one %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #ggtitle("Crop: Corn") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) 

Crop: Sugarbeet

result_one %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) + 
  scale_color_paletteer_d(Palette)

result_one %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) 

Crop: Others

result_one %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Others") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) + 
  scale_color_paletteer_d(Palette)

result_one %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Others") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) 

Strategy 2: Classification by RS

a_0, r_0, t_0

## Sample
N_0 <- 1000
a_0 <- get_a_sample_3st(U_star = U_0, N_0 = 1000, order = colnames(jacques_matrix), strategy = "two")

r_vector <- a_0 %>% 
  { matrix(colSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} %>%
  t() 
t_vector <- a_0 %>% 
  { matrix(rowSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} 

a_0 %>% as_tibble(rownames = 'crop') %>% 
  bind_rows(as_tibble(r_vector) %>% mutate(crop = 'R')) %>% 
  bind_cols(as_tibble(t_vector) %>% rename(T = value) %>% 
              bind_rows(tibble(`T` = 1000))) %>% 
  kable(caption = "a_0 confusion sample matrix") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) %>% 
  column_spec(c(1,7), bold = TRUE) %>% 
  row_spec(6, bold = TRUE)
a_0 confusion sample matrix
crop Wheat Rapeseed Corn Sugarbeet Others T
Wheat 127 44 11 4 10 196
Rapeseed 4 49 0 0 9 62
Corn 1 27 119 5 4 156
Sugarbeet 4 29 56 160 17 266
Others 64 51 14 31 160 320
R 200 200 200 200 200 1000

Get a artificial population for bootstrap sampling

U_star <- gen_U_star(a_0, strategy = 'two', size = Size)

Generate 1000 bootstrap replicas

#set.seed(112358)
result_two <- get_bootstrap(1000, U_star, N_0, R_U_0, strategy = 'two', order = colnames(jacques_matrix)) %>% 
  mutate(Sample = '1000')
write_rds(result_two, "data/result_two_N1m.rds")

Plot result

Palette <- palettes_d_names %>%
  filter(length == 3) %>%
  mutate(code = paste0(package, "::", palette)) %>%
  pull(code) %>%
  .[7] #3 5


result_two %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = 1, data = Reference) + 
  scale_color_paletteer_d(Palette)

result_two %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference) 

Bootstrap summary

table <- result_two %>% group_by(crop, estimator, Sample) %>% 
    summarise(Mean = mean(n), 
              Variance = var(n), 
              `Standard Deviation` = sd(n), 
              CV = `Standard Deviation`/Mean, 
              .groups = 'drop') %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  left_join(Reference %>% pivot_wider(names_from = References,  values_from = "n"), by = "crop") %>% 
  mutate(`T Bias` = `T` - Mean,
         `R Bias` = `R` - Mean) %>% 
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Mean = format(round(Mean), nsmall = 0, big.mark = ","),
         Variance = format(round(Variance), nsmall = 0, big.mark = ","),
         `Standard Deviation` = format(round(`Standard Deviation`), nsmall = 0, big.mark = ","),
         CV = format(round(CV, 4), nsmall = 2, big.mark = ","), 
         `T` = format(round(`T`), nsmall = 0, big.mark = ","),
         `T Bias` = format(round(`T Bias`), nsmall = 0, big.mark = ","),
         `R` = format(round(`R`), nsmall = 0, big.mark = ","),
         `R Bias` = format(round(`R Bias`), nsmall = 0, big.mark = ","))
table %>% 
  kable(caption = "Direct and Inverse Estimators Bootstrap Results (1000 rounds)", align = "lcrrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  column_spec(4, color = if_else(str_detect(table$Mean, "-"), 'red', 'black')) %>%  
  column_spec(7, color = if_else(str_detect(table$CV, "-"), 'red', 'black')) %>%  
  column_spec(10, color = if_else(str_detect(table$`T Bias`, "-"), 'red', 'black')) %>%  
  column_spec(11, color = if_else(str_detect(table$`R Bias`, "-"), 'red', 'black')) 
Direct and Inverse Estimators Bootstrap Results (1000 rounds)
Crop Estimator Sample Mean Variance Standard Deviation CV T R T Bias R Bias
Corn direct 1000 104,645 42,450,736 6,515 0.0623 99,736 135,024 -4,909 30,379
Corn inverse 1000 77,018 113,192,305 10,639 0.1381 99,736 135,024 22,718 58,006
Others direct 1000 394,267 203,358,078 14,260 0.0362 400,222 293,695 5,955 -100,572
Others inverse 1000 567,767 913,444,669 30,223 0.0532 400,222 293,695 -167,545 -274,072
Rapeseed direct 1000 42,951 36,297,509 6,025 0.1403 50,086 94,924 7,135 51,973
Rapeseed inverse 1000 -126,446 432,983,256 20,808 -0.1646 50,086 94,924 176,532 221,370
Sugarbeet direct 1000 210,681 95,180,363 9,756 0.0463 200,341 160,274 -10,340 -50,407
Sugarbeet inverse 1000 158,880 236,696,831 15,385 0.0968 200,341 160,274 41,461 1,394
Wheat direct 1000 247,457 147,498,694 12,145 0.0491 249,615 316,083 2,158 68,626
Wheat inverse 1000 322,781 517,950,752 22,759 0.0705 249,615 316,083 -73,166 -6,698

Other plots

Crop: Wheat

result_two %>% filter(crop == 'Wheat', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) + 
  scale_color_paletteer_d(Palette)

result_two %>% 
  filter(crop == 'Wheat', Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) 

Crop: Rapeseed

result_two %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) + 
  scale_color_paletteer_d(Palette)

result_two %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) 

Crop: Corn

result_two %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Corn") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) + 
  scale_color_paletteer_d(Palette)

result_two %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #ggtitle("Crop: Corn") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) 

Crop: Sugarbeet

result_two %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) + 
  scale_color_paletteer_d(Palette)

result_two %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) 

Crop: Others

result_two %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Others") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) + 
  scale_color_paletteer_d(Palette)

result_two %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #ggtitle("Crop: Others") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) 

Strategy 3: Classification by Ground

a_0, r_0, t_0

## Sample
N_0 <- 1000
a_0 <- get_a_sample_3st(U_star = U_0, N_0 = 1000, order = colnames(jacques_matrix), strategy = "three")

r_vector <- a_0 %>% 
  { matrix(colSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} %>%
  t() 
t_vector <- a_0 %>% 
  { matrix(rowSums(.), ncol = 1, dimnames = list(colnames(.), 'value'))} 

a_0 %>% as_tibble(rownames = 'crop') %>% 
  bind_rows(as_tibble(r_vector) %>% mutate(crop = 'R')) %>% 
  bind_cols(as_tibble(t_vector) %>% rename(T = value) %>% 
              bind_rows(tibble(`T` = 1000))) %>% 
  kable(caption = "a_0 confusion sample matrix") %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE, font_size = 12) %>% 
  column_spec(c(1,7), bold = TRUE) %>% 
  row_spec(6, bold = TRUE)
a_0 confusion sample matrix
crop Wheat Rapeseed Corn Sugarbeet Others T
Wheat 163 14 3 6 14 200
Rapeseed 37 124 0 0 39 200
Corn 3 25 150 15 7 200
Sugarbeet 7 15 37 117 24 200
Others 57 15 3 12 113 200
R 267 193 193 150 197 1000

Get a artificial population for bootstrap sampling

U_star <- gen_U_star(a_0, strategy = 'three', size = Size)

Generate 1000 bootstrap replicas

#set.seed(112358)
result_three <- get_bootstrap(1000, U_star, N_0, R_U_0, strategy = 'three', order = colnames(jacques_matrix)) %>% 
  mutate(Sample = '1000')
write_rds(result_three, "data/result_three_N1m.rds")

Plot result

Palette <- palettes_d_names %>%
  filter(length == 3) %>%
  mutate(code = paste0(package, "::", palette)) %>%
  pull(code) %>%
  .[7] #3 5


result_three %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = 1, data = Reference) + 
  scale_color_paletteer_d(Palette)

result_three %>% 
  filter(Sample == "1000") %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ crop, scales = 'free', ncol = 2) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference) 

Bootstrap summary

table <- result_three %>% group_by(crop, estimator, Sample) %>% 
    summarise(Mean = mean(n), 
              Variance = var(n), 
              `Standard Deviation` = sd(n), 
              CV = `Standard Deviation`/Mean, 
              .groups = 'drop') %>% 
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  left_join(Reference %>% pivot_wider(names_from = References,  values_from = "n"), by = "crop") %>% 
  mutate(`T Bias` = `T` - Mean,
         `R Bias` = `R` - Mean) %>% 
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Mean = format(round(Mean), nsmall = 0, big.mark = ","),
         Variance = format(round(Variance), nsmall = 0, big.mark = ","),
         `Standard Deviation` = format(round(`Standard Deviation`), nsmall = 0, big.mark = ","),
         CV = format(round(CV, 4), nsmall = 2, big.mark = ","), 
         `T` = format(round(`T`), nsmall = 0, big.mark = ","),
         `T Bias` = format(round(`T Bias`), nsmall = 0, big.mark = ","),
         `R` = format(round(`R`), nsmall = 0, big.mark = ","),
         `R Bias` = format(round(`R Bias`), nsmall = 0, big.mark = ","))
table %>% 
  kable(caption = "Direct and Inverse Estimators Bootstrap Results (1000 rounds)", align = "lcrrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  column_spec(5, color = if_else(str_detect(table$Mean, "-"), 'red', 'black')) %>%  
  column_spec(7, color = if_else(str_detect(table$CV, "-"), 'red', 'black')) %>%  
  column_spec(10, color = if_else(str_detect(table$`T Bias`, "-"), 'red', 'black')) %>%  
  column_spec(11, color = if_else(str_detect(table$`R Bias`, "-"), 'red', 'black')) 
Direct and Inverse Estimators Bootstrap Results (1000 rounds)
Crop Estimator Sample Mean Variance Standard Deviation CV T R T Bias R Bias
Corn direct 1000 147,040 27,286,155 5,224 0.0355 99,736 135,024 -47,304 -12,016
Corn inverse 1000 116,631 163,188,963 12,775 0.1095 99,736 135,024 -16,895 18,393
Others direct 1000 258,454 50,499,647 7,106 0.0275 400,222 293,695 141,768 35,241
Others inverse 1000 432,320 1,299,969,972 36,055 0.0834 400,222 293,695 -32,098 -138,625
Rapeseed direct 1000 162,789 53,864,162 7,339 0.0451 50,086 94,924 -112,703 -67,865
Rapeseed inverse 1000 29,006 340,608,588 18,456 0.6363 50,086 94,924 21,080 65,918
Sugarbeet direct 1000 202,142 45,905,873 6,775 0.0335 200,341 160,274 -1,801 -41,868
Sugarbeet inverse 1000 202,983 427,215,101 20,669 0.1018 200,341 160,274 -2,642 -42,709
Wheat direct 1000 229,574 48,887,796 6,992 0.0305 249,615 316,083 20,041 86,509
Wheat inverse 1000 219,060 830,305,239 28,815 0.1315 249,615 316,083 30,555 97,023

Other plots

Crop: Wheat

result_three %>% filter(crop == 'Wheat', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) + 
  scale_color_paletteer_d(Palette)

result_three %>% 
  filter(crop == 'Wheat', Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Wheat") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Wheat')) 

Crop: Rapeseed

result_three %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) + 
  scale_color_paletteer_d(Palette)

result_three %>% filter(crop == 'Rapeseed', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Rapeseed") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Rapeseed')) 

Crop: Corn

result_three %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Corn") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) + 
  scale_color_paletteer_d(Palette)

result_three %>% filter(crop == 'Corn', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #ggtitle("Crop: Corn") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Corn')) 

Crop: Sugarbeet

result_three %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) + 
  scale_color_paletteer_d(Palette)

result_three %>% filter(crop == 'Sugarbeet', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_brewer(palette = "Greys", direction = -1) +
  #ggtitle("Crop: Sugarbeet") +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Sugarbeet')) 

Crop: Others

result_three %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  #ggtitle("Crop: Others") +
  geom_vline(aes(xintercept = n/1000, color = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) + 
  scale_color_paletteer_d(Palette)

result_three %>% filter(crop == 'Others', Sample == '1000') %>% ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  theme(legend.position="bottom") + 
  labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', fill = 'Estimator') + 
  #ggtitle("Crop: Others") +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), 
             size = 1, data = Reference %>% filter(crop == 'Others')) 

Comparing Strategies

Crop: Wheat

Crop =  colnames(jacques_matrix)[1]

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = .75, data = Reference %>% filter(crop == Crop)) +
  scale_color_paletteer_d(Palette) 

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference %>% filter(crop == Crop)) 

Crop: Rapeseed

Crop =  colnames(jacques_matrix)[2]

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = .75, data = Reference %>% filter(crop == Crop)) +
  scale_color_paletteer_d(Palette) 

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference %>% filter(crop == Crop)) 

Crop: Corn

Crop =  colnames(jacques_matrix)[3]

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = .75, data = Reference %>% filter(crop == Crop)) +
  scale_color_paletteer_d(Palette) 

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference %>% filter(crop == Crop)) 

Crop: Sugarbeet

Crop =  colnames(jacques_matrix)[4]

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = .75, data = Reference %>% filter(crop == Crop)) +
  scale_color_paletteer_d(Palette) 

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference %>% filter(crop == Crop)) 

Crop: Others

Crop =  colnames(jacques_matrix)[5]

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) + 
  geom_vline(aes(xintercept = n/1000, color = References), size = .75, data = Reference %>% filter(crop == Crop)) +
  scale_color_paletteer_d(Palette) 

result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground")) %>% 
  filter(crop == Crop) %>% 
  filter(Sample == '1000') %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_wrap(~ Strategy, scales = "free_y", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference %>% filter(crop == Crop)) 

Previous results on sample size impact

result <- read_rds("data/result1k10k100k.rds")

Crop: Wheat

Crop =  colnames(jacques_matrix)[1]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) +
  geom_vline(aes(xintercept = n/1000, color = References, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop)) +
  scale_color_paletteer_d(Palette)

result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop))

Crop: Rapeseed

Crop =  colnames(jacques_matrix)[2]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) +
  geom_vline(aes(xintercept = n/1000, color = References, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop)) +
  scale_color_paletteer_d(Palette)

result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop))

Crop: Corn

Crop =  colnames(jacques_matrix)[3]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) +
  geom_vline(aes(xintercept = n/1000, color = References, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop)) +
  scale_color_paletteer_d(Palette)

Crop =  colnames(jacques_matrix)[3]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop))

Crop: Sugarbeet

Crop =  colnames(jacques_matrix)[4]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) +
  geom_vline(aes(xintercept = n/1000, color = References, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop)) +
  scale_color_paletteer_d(Palette)

result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop))

Crop: Others

Crop =  colnames(jacques_matrix)[5]
result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_paletteer_d(`"ggthemes::wsj_red_green"`) +
  geom_vline(aes(xintercept = n/1000, color = References, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop)) +
  scale_color_paletteer_d(Palette)

result %>%
  filter(crop == Crop) %>%
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') +
  facet_wrap(~ Sample, scales = "fixed", ncol = 1) +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density',
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = 1, data = Reference %>%
               filter(crop == Crop))

Outros Resultados

result <- result_one %>% 
  mutate(Strategy = "Strategy 1: Bivariate") %>% 
  bind_rows(result_two %>% mutate(Strategy = "Strategy 2: Classification by RS")) %>% 
  bind_rows(result_three %>% mutate(Strategy = "Strategy 3: Classification by Ground"))



result %>% filter(Strategy == "Strategy 1: Bivariate") %>%  
  group_by(crop, estimator, Strategy) %>% 
  summarise(Estimate = mean(n/1000),
            `Std.Dev.` = sd(n/1000),
            `CV(%)` = (`Std.Dev.`/Estimate)*100, .groups = 'drop') %>% 
  left_join(T_U_0 %>% as_tibble(rownames = "crop") %>% rename(T = n),
            by = 'crop') %>% 
  mutate(Bias = Estimate - T/1000) %>% 
  mutate(Estimate = round(Estimate, 1),
         `Std.Dev.` = round(`Std.Dev.`, 2),
         `CV(%)` = round(`CV(%)`, 1),
         Bias = round(Bias, 1)) %>% 
  select(-Strategy, -T) %>% 
  arrange(match(crop, colnames(jacques_matrix))) %>%
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Estimator = recode(Estimator, 
                              direct = 'Direct', 
                              inverse = 'Inverse')) %>% 
  kable(caption = "Strategy 1: Bivariate (x 1,000 hectares)", align = "lcrrrr") %>%
kable_classic_2() %>% 
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
              full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  row_spec(0, color = "black", background = 'white')
Strategy 1: Bivariate (x 1,000 hectares)
Crop Estimator Estimate Std.Dev. CV(%) Bias
Wheat Direct 243.1 10.45 4.3 -6.5
Wheat Inverse 225.6 23.33 10.3 -24.0
Rapeseed Direct 47.7 5.81 12.2 -2.4
Rapeseed Inverse 19.9 15.66 78.7 -30.2
Corn Direct 105.8 7.23 6.8 6.0
Corn Inverse 106.4 12.46 11.7 6.6
Sugarbeet Direct 198.2 9.54 4.8 -2.1
Sugarbeet Inverse 198.0 17.76 9.0 -2.4
Others Direct 405.3 13.01 3.2 5.1
Others Inverse 450.2 28.10 6.2 50.0
result %>% filter(Strategy == "Strategy 2: Classification by RS") %>%  
  group_by(crop, estimator, Strategy) %>% 
  summarise(Estimate = mean(n/1000),
            `Std.Dev.` = sd(n/1000),
            `CV(%)` = (`Std.Dev.`/Estimate)*100, .groups = 'drop') %>% 
  left_join(T_U_0 %>% as_tibble(rownames = "crop") %>% rename(T = n),
            by = 'crop') %>% 
  mutate(Bias = Estimate - T/1000) %>% 
  mutate(Estimate = round(Estimate, 1),
         `Std.Dev.` = round(`Std.Dev.`, 2),
         `CV(%)` = round(`CV(%)`, 1),
         Bias = round(Bias, 1)) %>% 
  select(-Strategy, -T) %>% 
  arrange(match(crop, colnames(jacques_matrix))) %>%
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Estimator = recode(Estimator, 
                              direct = 'Direct', 
                              inverse = 'Inverse')) %>% 
  kable(caption = "Strategy 2: Classification by RS (x 1,000 hectares)", align = "lcrrrr") %>%
kable_classic_2() %>% 
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
              full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  row_spec(0, color = "black", background = 'white')
Strategy 2: Classification by RS (x 1,000 hectares)
Crop Estimator Estimate Std.Dev. CV(%) Bias
Wheat Direct 247.5 12.14 4.9 -2.2
Wheat Inverse 322.8 22.76 7.1 73.2
Rapeseed Direct 43.0 6.02 14.0 -7.1
Rapeseed Inverse -126.4 20.81 -16.5 -176.5
Corn Direct 104.6 6.52 6.2 4.9
Corn Inverse 77.0 10.64 13.8 -22.7
Sugarbeet Direct 210.7 9.76 4.6 10.3
Sugarbeet Inverse 158.9 15.38 9.7 -41.5
Others Direct 394.3 14.26 3.6 -6.0
Others Inverse 567.8 30.22 5.3 167.5
result %>% filter(Strategy == "Strategy 3: Classification by Ground") %>%  
  group_by(crop, estimator, Strategy) %>% 
  summarise(Estimate = mean(n/1000),
            `Std.Dev.` = sd(n/1000),
            `CV(%)` = (`Std.Dev.`/Estimate)*100, .groups = 'drop') %>% 
  left_join(T_U_0 %>% as_tibble(rownames = "crop") %>% rename(T = n),
            by = 'crop') %>% 
  mutate(Bias = Estimate - T/1000) %>% 
  mutate(Estimate = round(Estimate, 1),
         `Std.Dev.` = round(`Std.Dev.`, 2),
         `CV(%)` = round(`CV(%)`, 1),
         Bias = round(Bias, 1)) %>% 
  select(-Strategy, -T) %>% 
  arrange(match(crop, colnames(jacques_matrix))) %>%
  rename(Crop = crop, Estimator = estimator) %>% 
  mutate(Estimator = recode(Estimator, 
                              direct = 'Direct', 
                              inverse = 'Inverse')) %>% 
  kable(caption = "Strategy 3: Classification by Ground (x 1,000 hectares)", align = "lcrrrr") %>%
kable_classic_2() %>% 
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
              full_width = FALSE, font_size = 12) %>% 
  column_spec(1, bold = TRUE) %>% 
  row_spec(0, color = "black", background = 'white')
Strategy 3: Classification by Ground (x 1,000 hectares)
Crop Estimator Estimate Std.Dev. CV(%) Bias
Wheat Direct 229.6 6.99 3.0 -20.0
Wheat Inverse 219.1 28.82 13.2 -30.6
Rapeseed Direct 162.8 7.34 4.5 112.7
Rapeseed Inverse 29.0 18.46 63.6 -21.1
Corn Direct 147.0 5.22 3.6 47.3
Corn Inverse 116.6 12.77 11.0 16.9
Sugarbeet Direct 202.1 6.78 3.4 1.8
Sugarbeet Inverse 203.0 20.67 10.2 2.6
Others Direct 258.5 7.11 2.7 -141.8
Others Inverse 432.3 36.06 8.3 32.1
result %>%  
  mutate(crop = factor(crop, levels = colnames(jacques_matrix), ordered = TRUE)) %>% 
  ggplot() +
  theme_bw() +
  geom_density(aes(x = n/1000, fill = estimator), alpha = .75, position = 'identity') + 
  facet_grid(crop ~ Strategy, scales = "free") +
  theme(legend.position="bottom") + labs(x = 'Estimates (x 1000 pixels)', y  = 'Density', 
                                         fill = 'Estimator') +
  scale_fill_brewer(palette = "Greys", direction = -1) +
  geom_vline(aes(xintercept = n/1000, linetype = References), size = .75, data = Reference) 

Stay Tuned

Please visit the CastLAB page for latest updates and news. Comments, bug reports and pull requests are always welcome.

About

Code Repository for Bootstrap Assessment of Crop Area Estimates using Satellite Pixels Counting

Resources

License

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published