0havest_area.R

library(readxl)  
library(sf)
library(ggplot2)
library(dplyr)
#AMC and its shape  - 2010 and 1991
amc_mun_sf_00 <- st_read("G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/AMC_sf.gpkg")
amc_mun_sf_00 <- amc_mun_sf_00 %>% mutate(MUNICIPALITY_CODE =as.numeric(MUNICIPALITY_CODE ))


############ use data from the site - Agricultural Census 1996 (ref 1995)
#######https://sidra.ibge.gov.br/Tabela/501
#######https://sidra.ibge.gov.br/Tabela/512
### can also use package library(sidrar)   ==> limit of 50.000 for use of API

t501 <- read_xlsx("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/SIDRA/tabela501.xlsx")
t512 <- read_xlsx("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/SIDRA/tabela512.xlsx")

censo95 <- left_join(t501, t512, by = "Cód.")
names(censo95) <- iconv(names(censo95), from = "UTF-8", to = "ASCII//TRANSLIT")

censo95 <- censo95 %>%
  mutate(across(everything(), ~as.numeric(as.character(.))))
censo95 <- censo95 %>%
  mutate(
    area95_ban = rowSums(select(., contains("banana")), na.rm = TRUE),
    area95_cit = rowSums(select(., contains("laranja")), na.rm = TRUE),
    area95_coc = rowSums(select(., contains("cacau")), na.rm = TRUE),
    area95_cof = rowSums(select(., contains("cafe")), na.rm = TRUE),  # cuidado com acento
    area95_cot= rowSums(select(., contains("algodao")), na.rm = TRUE),
    area95_mze = rowSums(select(., contains("milho")), na.rm = TRUE),
    area95_rcd = rowSums(select(., contains("arroz")), na.rm = TRUE),
    area95_soy = rowSums(select(., contains("soja")), na.rm = TRUE),
    area95_suc = rowSums(select(., contains("cana")), na.rm = TRUE),
    area95_tob = rowSums(select(., contains("fumo")), na.rm = TRUE),
    area95_whe = rowSums(select(., contains("trigo")), na.rm = TRUE)
  )

censo95 <- censo95 %>% rename(MUNICIPALITY_CODE = "Cod.")
censo95 <- censo95 %>% select(MUNICIPALITY_CODE, area95_ban, area95_cit, area95_coc, area95_cof, 
                              area95_cot, area95_mze, area95_rcd, area95_soy, area95_suc, area95_tob, area95_whe)


############ SIDRA 2000
########If I use this, nee to have data from previous years and do *avarage*. Using only 1 year will give a lot of empty states 

######## plantation
# https://sidra.ibge.gov.br/Tabela/5457 => PLANTATION
# https://sidra.ibge.gov.br/tabela/3939 => HEAD OF PASTO
# minerais: https://www.ibge.gov.br/novo-portal-destaques/33322-substituicao-das-tabelas-no-sidra-para-a-pesquisa-industrial-mensal-producao-fisica-pim-pf.html
#https://analisemacro.com.br/economia/pib/analise-da-producao-industrial-com-o-r/
sidra2000 <- read_xlsx("G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/SIDRA/platacao_2000_hc_munic.xlsx") 
sidra2000 <- sidra2000  %>% select(-Total) %>% mutate(across(everything(), ~as.numeric(as.character(.)))) 
sidra2000 <- sidra2000 %>%
  rename(
    area00_ban = banana,
    area00_cit = orange,
    area00_coc = cocoa,
    area00_cof = coffee,
    area00_cot= cotton,
    area00_mze = maize,
    area00_rcd = rice,
    area00_soy = soybean,
    area00_suc = sugar,
    area00_tob = tobacco,
    area00_whe = wheat)

################# comparing with mean 95-99 : NOT MUCH DIFFERENCE
sidra_95_99 <- read.csv2("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/SIDRA/tabela5457_95_99.csv") 
names(sidra_95_99) <- iconv(names(sidra_95_99), from = "UTF-8", to = "ASCII//TRANSLIT")
sidra_95_99 <- sidra_95_99 %>% mutate(across(everything(), ~as.numeric(as.character(.))))
sidra_95_99 <- sidra_95_99 %>%
  mutate(
    area_m_ban = rowSums(select(., contains("banana")), na.rm = TRUE),
    area_m_cit = rowSums(select(., contains("laranja")), na.rm = TRUE),
    area_m_coc = rowSums(select(., contains("cacau")), na.rm = TRUE),
    area_m_cof = rowSums(select(., contains("cafe")), na.rm = TRUE),  # cuidado com acento
    area_m_cot= rowSums(select(., contains("algodao")), na.rm = TRUE),
    area_m_mze = rowSums(select(., contains("milho")), na.rm = TRUE),
    area_m_rcd = rowSums(select(., contains("arroz")), na.rm = TRUE),
    area_m_soy = rowSums(select(., contains("soja")), na.rm = TRUE),
    area_m_suc = rowSums(select(., contains("cana")), na.rm = TRUE),
    area_m_tob = rowSums(select(., contains("fumo")), na.rm = TRUE),
    area_m_whe = rowSums(select(., contains("trigo")), na.rm = TRUE)
  )

sidra_95_99 <- sidra_95_99 %>% rename(MUNICIPALITY_CODE = "Cod.")
sidra_95_99 <- sidra_95_99 %>% select(MUNICIPALITY_CODE, area_m_ban, area_m_cit, area_m_coc, area_m_cof, 
                                      area_m_cot, area_m_mze, area_m_rcd, area_m_soy, area_m_suc, area_m_tob, area_m_whe)

all_area <- left_join(sidra_95_99, sidra2000, by = "MUNICIPALITY_CODE")
all_area <- all_area %>% filter(!is.na(MUNICIPALITY_CODE))  %>%  filter(!(MUNICIPALITY_CODE==0))
all_area <- left_join(all_area, censo95, by = "MUNICIPALITY_CODE")
all_area <- all_area %>% filter(!is.na(MUNICIPALITY_CODE))  %>%  filter(!(MUNICIPALITY_CODE==0))

all_area <- left_join(all_area, amc_mun_sf_00, by = "MUNICIPALITY_CODE")

## gENERATING VALUE BY AMC
all_area <- all_area %>% 
  select(-MUNICIPALITY_CODE) %>%
  mutate(across(.cols = !c(id_amc, geom), .fns = ~ as.numeric(.)))%>%
  group_by(id_amc, geom) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop")


FAO2000 <- read.csv("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/FAO_GAEZ/muni_with_yield_actual2000_values_amc.csv") 
FAO2000 <- FAO2000 %>% rename(id_amc = code_amc) %>% mutate(id_amc = as.numeric(id_amc)) 
all_area <- left_join(all_area, FAO2000, by = "id_amc")

cor(all_area$soy_2000_har, all_area$area95_soy, use = "complete.obs") #FAO 2000 and CENSUS 95    0.471895
cor(all_area$area00_soy, all_area$area_m_soy, use = "complete.obs") #PMM 2000 and mean 95-99      0.9847395
cor(all_area$soy_2000_har, all_area$area_m_soy, use = "complete.obs")  #FAO 2000 andmean 95-99   0.4628705
cor(all_area$area95_soy, all_area$area_m_soy, use = "complete.obs")  #FAO 2000 and mean 95-99    0.9847926


cor(all_area$soy80, all_area$soy, use = "complete.obs") #PMM 1980 and CENSO 95                0.5500856
cor(all_area$soybean, all_area$soy, use = "complete.obs") #PMM 2000 and CENSO 95              0.9635124
cor(all_area$soy_2000_har, all_area$soy, use = "complete.obs")  #FAO 2000 and CENSO 95        0.471895
cor(all_area$soy_2000_har, all_area$soybean, use = "complete.obs")  #FAO 2000 and PMM 2000    0.4336284

all_area_map <- st_as_sf(all_area) 

ggplot() +
  geom_sf(data = all_area_map, aes(fill =  as.numeric(area95_soy)), color = NA) +
  scale_fill_viridis_c(option = "magma", na.value = "white") +  # Enhancing visualization
  theme_void() +
  labs(fill = "Hac by municipality - Censo 1995") 

all_area <- st_as_sf(all_area) 
all_area <- st_drop_geometry(all_area)
write.csv(all_area, "G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/SIDRA/all_area.csv", row.names = TRUE)
all_area <- read.csv("G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/SIDRA/all_area.csv") 




################ took out SIDRA 1980 as municiplaity / amc shoudl be different
AMC <- read.csv("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/br_ibge_amc_municipio_de_para.csv")

AMC_1980 <- AMC %>%
  filter(ano_de == 1980)%>%
  filter(ano_para == 2010) %>% 
  select(id_municipio, id_amc)  %>%
  rename(MUNICIPALITY_CODE = id_municipio) %>%
  mutate(MUNICIPALITY_CODE = as.numeric(MUNICIPALITY_CODE))

sidra1980 <- read_xlsx("G:/Meu Drive/MASTER 23 25/Tese/0_these/0raw_data/SIDRA/tabela5457_COLHIDA.xlsx") 
sidra1980 <- sidra1980  %>% mutate(across(everything(), ~as.numeric(as.character(.))))  %>% 
  mutate(
    area80_ban = rowSums(select(., contains("banana")), na.rm = TRUE),
    area80_cit = rowSums(select(., contains("laranja")), na.rm = TRUE),
    area80_coc = rowSums(select(., contains("cacau")), na.rm = TRUE),
    area80_cof = rowSums(select(., contains("cafe")), na.rm = TRUE),  # cuidado com acento
    area80_cot = rowSums(select(., contains("algodao")), na.rm = TRUE),
    area80_mze = rowSums(select(., contains("milho")), na.rm = TRUE),
    area80_rcd = rowSums(select(., contains("arroz")), na.rm = TRUE),
    area80_soy = rowSums(select(., contains("soja")), na.rm = TRUE),
    area80_suc = rowSums(select(., contains("cana")), na.rm = TRUE),
    area80_tob = rowSums(select(., contains("fumo")), na.rm = TRUE),
    area80_whe = rowSums(select(., contains("trigo")), na.rm = TRUE)
  )

sidra1980 <- sidra1980 %>% rename(MUNICIPALITY_CODE = `Cód.`) %>%
  select(MUNICIPALITY_CODE, area80_ban, area80_cit, area80_coc, area80_cof, area80_cot, area80_mze, 
         area80_rcd, area80_soy, area80_suc, area80_tob, area80_whe)
sidra1980 <- left_join(sidra1980, AMC_1980, by = "MUNICIPALITY_CODE")
sidra1980_g <- sidra1980 %>% 
  select(-MUNICIPALITY_CODE) %>%
  mutate(across(.cols = !c(id_amc), .fns = ~ as.numeric(.)))%>%
  group_by(id_amc) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop")

sidra1980_g <- left_join(AMC_1980, sidra1980_g, by = "id_amc")  #AMC 80-10 WITH MUNICIPALITY CODE 2010

AMC_1980_91 <- AMC %>%
  filter(ano_de == 1991)%>%
  filter(ano_para == 2010) %>% 
  select(id_municipio, id_amc)  %>%
  rename(MUNICIPALITY_CODE = id_municipio) %>%
  mutate(MUNICIPALITY_CODE = as.numeric(MUNICIPALITY_CODE)) #AMC 91-10, MUNIC_CODE 2010
sidra1980_g <- sidra1980_g %>% rename(id_amc_80 = id_amc)

sidra1980_g <- left_join(AMC_1980_91, sidra1980_g, by = "MUNICIPALITY_CODE")
sidra1980_g <- sidra1980_g %>% select(-id_amc_80, -MUNICIPALITY_CODE)

all_area2 <- left_join(all_area, sidra1980_g, by = "id_amc") 

#CHECKING IF SAME AMC HAVE DIFF VALUES
K<-all_area2 %>%
  group_by(id_amc) %>%
  summarise(across(everything(), ~ n_distinct(.))) %>%
  filter(if_any(-id_amc, ~ . > 1))    # ONLY 2 AMC ARE REPEATED WITH DIF VALUES, IN 1980. IGNORE THATA

all_area2 <- all_area2 %>%
  distinct(id_amc, .keep_all = TRUE)

write.csv(all_area2, "G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/SIDRA/all_area2.csv", row.names = TRUE)


############### maps
library(sf)

all_area2 <- read.csv("G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/SIDRA/all_area2.csv")
#amc
AMC_sf <- st_read("G:/Meu Drive/MASTER 23 25/Tese/0_these/1data/AMC_sf.gpkg")
#state limit
states_sf <- geobr::read_state(year = 2010, showProgress = TRUE) 


all_area2 <- left_join(AMC_sf, all_area2, by = "id_amc")
all_area2 <- st_as_sf(all_area2) 


library(ggplot2)
library(patchwork)  # para combinar os mapas
library(scales)     # para formatação se necessário

# Lista de sufixos e nomes legíveis das commodities
commodity_vars <- c("ban", "cit", "coc", "cof", "cot", "mze", "rcd", "soy", "suc", "tob", "whe")
commodity_names <- c(
  ban = "Banana", cit = "Orange", coc = "Cocoa", cof = "Coffee",
  cot = "Cotton", mze = "Maize", rcd = "Rice", soy = "Soybean",
  suc = "Sugarcane", tob = "Tobacco", whe = "Wheat"
)

# Função para criar 1 mapa
plot_map <- function(var, name) {
  all_area2 <- all_area2 %>%
    mutate(tmp_area = na_if(.data[[paste0("area_m_", var)]], 0))  # trata zeros como NA
  
  ggplot() +
    geom_sf(data = all_area2, aes(fill = tmp_area / 1000), color = NA) +
    geom_sf(data = states_sf, fill = NA, color = "lightgray", size = 0.5) +
    scale_fill_gradient(
      low = "beige", high = "darkgreen",
      name = "ha",
      na.value = "white"  
    ) +
    labs(title = name) +
    theme_void() +
    theme(
      legend.position = "bottom",
      legend.key.width = unit(0.8, "cm"),
      plot.title = element_text(hjust = 0.5, size = 12)
    )
}


# Criar todos os mapas
plots <- mapply(plot_map, var = commodity_vars, name = commodity_names, SIMPLIFY = FALSE)

# Organizar: 6 mapas em cima, 5 embaixo
top_row <- wrap_plots(plots[1:6], ncol = 6)
bottom_row <- wrap_plots(plots[7:11], ncol = 5)

# Combinar com título geral
final_plot <- (top_row / bottom_row) +
  plot_annotation(
    #title = "Harvest area at the municipality level, average across 1995–1999, in 1000 hac",
    theme = theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))
  )

# Exibir
#final_plot
ggsave("G:/Meu Drive/MASTER 23 25/Tese/0_these/3output/descriptive/harvest_area_maps_beige_green1.png", final_plot, width = 16, height = 8, dpi = 500)


plot_map <- function(var, name) {
  all_area2 <- all_area2 %>%
    mutate(tmp_area = na_if(.data[[paste0("area_m_", var)]], 0))  # trata zeros como NA
  ggplot() +
    geom_sf(data = all_area2, aes(fill = (.data[[paste0("area_m_", var)]]/1000)), color = NA) +
    geom_sf(data = states_sf, fill = NA, color = "lightgray", size = 0.05) +
    scale_fill_distiller(
      palette = "Greens", direction = 1,
      na.value = "white",
      name = "ha"
    ) +
    labs(title = name) +
    theme_void() +
    theme(
      legend.position = "bottom",
      legend.key.width = unit(0.5, "cm"),
      plot.title = element_text(hjust = 0.5, size = 12)
    )
}


# Criar todos os mapas
plots <- mapply(plot_map, var = commodity_vars, name = commodity_names, SIMPLIFY = FALSE)

# Organizar: 6 mapas em cima, 5 embaixo
top_row <- wrap_plots(plots[1:6], ncol = 6)
bottom_row <- wrap_plots(plots[7:11], ncol = 5)

# Combinar com título geral
final_plot <- (top_row / bottom_row) +
  plot_annotation(
    #title = "Harvest area at the municipality level, average across 1995–1999, in 1000 hac",
    theme = theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))
  )

ggsave("G:/Meu Drive/MASTER 23 25/Tese/0_these/3output/descriptive/harvest_area_maps_greens.png", final_plot, width = 16, height = 8, dpi = 500)