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

and 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.

Files

README.md

Latest commit

History

README.md

File metadata and controls

Bootstrap Assessment of Crop Area Estimates using Satellite Pixels Counting

Packages

Function to generate artificial population

Function to get a sample from artificial population

Function to get estimates from sample confution matrix and artificial population

Function to get bootstrap replicas

Simulated data

and vectors

Get a artificial population and associated R vector

T_0 and R_0 tables

Strategy 1: Bivariate

a_0, r_0, t_0

Get a artificial population for bootstrap sampling

Generate 1000 bootstrap replicas

Plot result

Bootstrap summary

Other plots

Crop: Wheat

Crop: Rapeseed

Crop: Corn

Crop: Sugarbeet

Crop: Others

Strategy 2: Classification by RS

a_0, r_0, t_0

Get a artificial population for bootstrap sampling

Generate 1000 bootstrap replicas

Plot result

Bootstrap summary

Other plots

Crop: Wheat

Crop: Rapeseed

Crop: Corn

Crop: Sugarbeet

Crop: Others

Strategy 3: Classification by Ground

a_0, r_0, t_0

Get a artificial population for bootstrap sampling

Generate 1000 bootstrap replicas

Plot result

Bootstrap summary

Other plots

Crop: Wheat

Crop: Rapeseed

Crop: Corn

Crop: Sugarbeet

Crop: Others

Comparing Strategies

Crop: Wheat

Crop: Rapeseed

Crop: Corn

Crop: Sugarbeet

Crop: Others

Previous results on sample size impact

Crop: Wheat

Crop: Rapeseed

Crop: Corn

Crop: Sugarbeet

Crop: Others

Outros Resultados

Stay Tuned