Vamos carregar tudo que for necessário para as análises

Pacotes

library(geobr)
library(sf)
library(dplyr)
library(ggthemes)
library(tmap)    # for static and interactive maps
library(leaflet) # for interactive maps
#library(mapview) # for interactive maps
library(ggplot2) # tidyverse data visualization package
library(shiny)   # for web applications
library(colorspace)
library(tidyverse)
library(stargazer)
library(plotrix)
library(cowplot)
library(knitr)
library(rmarkdown)
library(shapefiles)
library(kableExtra)
library(janitor)

Dados

Neste caso estou carregando dados de outro

# DATA ----
#vsd_prop <- read.csv("Results/vsd_prop.txt", row.names=1, sep="")
vsd_prop_simul <- read.csv("Results/vsd_prop_simul.txt", row.names=1, sep="")
mun_rds_final<-read.table("Results/mun_caat_CF_19.08.txt",head=T)


# Aqui podemos adicionar novas variáveis como a variação da vegetação entre 2008 e 2019 por propriedade
vsd_prop_simul %>% 
  mutate(tx_desm=NV19.1-NV08.1) %>%
  #mutate(tx_desm2=tx_desm/NV19.1*100) %>% 
  mutate (propSize2 = cut(area_MF, breaks= c(-Inf, 4, Inf), labels = c("Small", "Large"))) %>%
  mutate(landtenure2 = as.factor(ifelse(landtenure=="CARpo", "Private",
                              ifelse(landtenure=="CARpr", "Private",
                              ifelse(landtenure=="SIGEF", "Private",
                              ifelse(landtenure=="SI_PL", "Simulated",
                              ifelse(landtenure == "ARU", "Rural settlements", "Quilombola"))))))) %>% 
  mutate (propSize3 = as.factor(ifelse(landtenure2 == "Private" & propSize2 == "Small", "Small", # propSize3 é a classificação considerando quilombolas e assentamentos como pequenas
                              ifelse(landtenure2 == "Private" & propSize2 == "Large", "Large",
                              ifelse(landtenure2 == "Simulated" & propSize2 == "Small", "Small",
                              ifelse(landtenure2 == "Simulated" & propSize2 == "Large", "Large",
                              ifelse(landtenure2 == "ARU" , "Small", "Small"))))))) %>% 
  mutate(app_def=area_app.1-NV_app19.1)->vsd_prop
glimpse(vsd_prop)

## Rows: 1,689,165
## Columns: 25
## $ mun_cod        <int> 2509701, 2211605, 2208874, 2923050, 2604155, 2909307, …
## $ nome_Mun       <chr> "MONTEIRO", "VILA NOVA DO PIAUÍ", "RIBEIRA DO PIAUÍ", …
## $ UF             <chr> "PB", "PI", "PI", "BA", "PE", "BA", "BA", "RN", "BA", …
## $ id_prop        <dbl> 8, 16, 17, 18, 20, 21, 22, 23, 29, 30, 38, 42, 44, 47,…
## $ area_prop.MB.1 <dbl> 2.61, 9.72, 381.51, 3.51, 8.55, 17.10, 0.99, 31.41, 1.…
## $ area_app.1     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ NV19.1         <dbl> 1.44, 8.28, 358.92, 0.00, 0.00, 2.97, 0.00, 23.85, 0.0…
## $ NV_app19.1     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ NV08.1         <dbl> 1.53, 5.58, 362.61, 0.00, 0.00, 1.44, 0.00, 24.84, 1.4…
## $ landtenure     <chr> "CARpr", "CARpr", "CARpr", "CARpr", "CARpr", "CARpr", …
## $ mf_incra       <int> 55, 70, 70, 30, 26, 65, 50, 70, 65, 35, 70, 60, 70, 70…
## $ totmf          <dbl> 0.04, 0.14, 5.37, 0.12, 0.34, 0.26, 0.02, 0.45, 0.03, …
## $ farmsizerange  <chr> "S", "S", "M", "S", "S", "S", "S", "S", "S", "S", "S",…
## $ MFmun_ha       <int> 55, 70, 70, 30, 26, 65, 50, 70, 65, 35, 70, 60, 70, 70…
## $ area_MF        <dbl> 0.05, 0.14, 5.45, 0.12, 0.33, 0.26, 0.02, 0.45, 0.03, …
## $ propSize       <chr> "S", "S", "M", "S", "S", "S", "S", "S", "S", "S", "S",…
## $ RL20_ha        <dbl> 0.522, 1.944, 76.302, 0.702, 1.710, 3.420, 0.198, 6.28…
## $ vsd_ha         <dbl> 0.918, 6.336, 282.618, -0.702, -1.710, -0.450, -0.198,…
## $ RL_Perc08      <dbl> 58.62, 57.41, 95.05, 0.00, 0.00, 8.42, 0.00, 79.08, 76…
## $ vsd_total      <dbl> 0.918, 6.336, 282.618, 0.000, 0.000, 1.530, 0.000, 17.…
## $ tx_desm        <dbl> -0.09, 2.70, -3.69, 0.00, 0.00, 1.53, 0.00, -0.99, -1.…
## $ propSize2      <fct> Small, Small, Large, Small, Small, Small, Small, Small…
## $ landtenure2    <fct> Private, Private, Private, Private, Private, Private, …
## $ propSize3      <fct> Small, Small, Large, Small, Small, Small, Small, Small…
## $ app_def        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

levels(vsd_prop$propSize3) # esse propSize3 já está corrigido, considerando os assentamentos e terras quilombolas

## [1] "Large" "Small"

Análises descritivas para propriedades privadas excluindo ARU e QL (Assentamentos e Terras Quilombolas)

Quanta área está em mãos de grandes e pequenos proprietários segundo o tipo de titulação ?

vsd_prop %>% 
  filter(area_prop.MB.1>0) %>% 
  group_by(landtenure2) %>% 
  summarise_at(vars(area_prop.MB.1), funs(sum, n())) %>%
  rename(Type_of_proprerty = landtenure2, Area_of_properties=sum, Number_of_properties=n) %>% 
  rename_with( ~ tolower(gsub("_", " ", .x, fixed = TRUE))) %>% 
  adorn_totals("row") %>% 
  kable(caption = "Table 1. Share of area and number of farmlands among types of properties in the Caatinga") %>% 
  kable_paper(bootstrap_options = "striped", full_width = F)

Table 1. Share of area and number of farmlands among types of properties in the Caatinga
type of proprerty	area of properties	number of properties
Private	33206542	996182
Quilombola	396252	70
Rural settlements	3882717	1838
Simulated	43961555	668865
Total	81447066	1666955

blank_theme <- theme_minimal()+
  theme(
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  panel.border = element_blank(),
  panel.grid=element_blank(),
  axis.ticks = element_blank(),
  plot.title=element_text(size=14, face="bold")
  )
library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:plotrix':
## 
##     rescale

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

vsd_prop %>% 
  filter(area_prop.MB.1>0) %>% 
  group_by(landtenure2) %>% 
  summarise_at(vars(area_prop.MB.1), funs(sum, n())) %>% 
  ggplot(aes(x="", y=sum, fill=landtenure2))+
  geom_bar(stat = "identity", width = 1, color="white")+
  coord_polar("y", start=0)+
  blank_theme +
  theme(axis.text.x=element_blank())

Quanta Caatinga (vegetação) existe em propriedades segundo seu tamanho?

vsd_prop %>% 
  group_by(landtenure2) %>%
  filter(area_prop.MB.1>0) %>%
  summarise_at(vars(NV19.1), funs(sum, n())) %>% 
  rename(Type_of_proprerty = landtenure2, Caatinga_ha=sum, Number_of_properties=n) %>% 
  rename_with( ~ tolower(gsub("_", " ", .x, fixed = TRUE))) %>% 
  adorn_totals("row") %>% 
  kable(caption = "Table 2. Share of natural vegetation and number of farmlands among types of properties in the Caatinga") %>% 
  kable_paper(bootstrap_options = "striped", full_width = F)

Table 2. Share of natural vegetation and number of farmlands among types of properties in the Caatinga
type of proprerty	caatinga ha	number of properties
Private	19823458	996182
Quilombola	245588	70
Rural settlements	2692655	1838
Simulated	27825523	668865
Total	50587224	1666955

vsd_prop %>% 
  group_by(landtenure2) %>%
  filter(area_prop.MB.1>0) %>%
  summarise_at(vars(NV19.1), funs(sum, n())) %>% 
  ggplot(aes(x="", y=sum, fill=landtenure2))+
  geom_bar(stat = "identity", width = 1, color="white")+
  coord_polar("y", start=0)+
  blank_theme +
  theme(axis.text.x=element_blank())

Em resumo, temos terras em sua imensa maioria privadas e mais da metade não possui informação que permita ter segurança de propŕiedade (Table 1). O mesmo padrão portanto, se repete para a quantidade de vegetação (Table 2). Dá pra resumir essa tabela tods em uma só e talvez concentrar essa pŕimeira parte numa descrição da questão fundiária e vegetação da Caatinga, ainda sem focar na LPVN. Ver abaixo

vsd_prop %>% 
  group_by(landtenure2) %>%
  filter(area_prop.MB.1>0) %>%
  summarise_at(vars(area_prop.MB.1, NV19.1), funs(sum, n())) %>% 
  mutate(Perc_veg=NV19.1_sum/area_prop.MB.1_sum*100) %>% 
  select(1:4,6) %>% 
  rename(Type_of_proprerty = landtenure2,
         Area_of_properties = area_prop.MB.1_sum, 
         Vegetation_in_ha = NV19.1_sum, 
         Number_of_properties = area_prop.MB.1_n) %>% 
  rename_with( ~ tolower(gsub("_", " ", .x, fixed = TRUE))) %>% 
  adorn_totals("row") %>% 
  kable(caption = "Table 3. Share of natural vegetation and number of farmlands among types of properties in the Caatinga") %>% 
  kable_paper(bootstrap_options = "striped", full_width = F)

Table 3. Share of natural vegetation and number of farmlands among types of properties in the Caatinga
type of proprerty	area of properties	vegetation in ha	number of properties	perc veg
Private	33206542	19823458	996182	59.7
Quilombola	396252	245588	70	62.0
Rural settlements	3882717	2692655	1838	69.3
Simulated	43961555	27825523	668865	63.3
Total	81447066	50587224	1666955	254.3

# atenção: desconsiderar o total de perc_veg que e na verdade 62.1%

Agora vamos ao deficit legal de APP e RL

Nestas análises vamos destrinchar por categoria de tamanho de prorpiedade. Dessa vez, os simulados vão ficar de fora, mas acredito que podemos pensar que o padrão encontrado para o resto da Caaringa, deve se repetir para eles.

Qual o deficit por propriedade da Caatinga?

Note que aqui vamos tomando somente aquelas propriedade com “vsd_total” negativo.

vsd_prop %>% 
  group_by(UF,landtenure2) %>% 
  filter(area_prop.MB.1>0) %>%
  filter(vsd_total<0) %>% # resumo aqui àquelas propriedades que apresentam déficit
  summarise_at(vars(vsd_total), sum) %>% 
  ggplot(aes(landtenure2, vsd_total))+
  geom_col()+facet_wrap(~UF)

# Com simulados

vsd_prop %>% 
  #group_by(UF,propSize3) %>% 
  filter(landtenure2 == "Simulated" ) %>% 
  filter(area_prop.MB.1>0) %>% # não muda nada porque as simuladas não possuem déficit algum, porque provavelmente não podemos calcular
  filter(vsd_total<0) %>% # resumo aqui àquelas propriedades que apresentam déficit
  summarise_at(vars(vsd_total), sum)

##   vsd_total
## 1         0

Em forma de tabela, temos esses valores de déficit de RL

vsd_prop %>% 
  group_by(landtenure2, propSize3) %>%
  filter(area_prop.MB.1>0) %>%
  filter(vsd_total<0) %>%
  summarise_at(vars(vsd_total), funs(sum,mean, n())) %>%
  unite(land_prop, landtenure2:propSize3, remove = FALSE) %>%
  #group_by(UF,land_prop) %>%
  #summarise_at(vars(vsd_total), sum) %>%
  #spread(land_prop, vsd_total) %>%
  #mutate_if(is.numeric, funs(replace_na(., 0))) %>%
  adorn_totals("row") %>%
  kable(caption ="Table 5. Deficit of natural vegetation between large and smallholders in legal reserves (LR) among states in the Caatinga") %>% 
  kable_paper(bootstrap_options = "striped", full_width = F)

Table 5. Deficit of natural vegetation between large and smallholders in legal reserves (LR) among states in the Caatinga
land_prop	landtenure2	propSize3	sum	mean	n
Private_Large	Private	Large	-139644	-46.77	2986
Private_Small	Private	Small	-143501	-1.02	141144
Quilombola_Small	Quilombola	Small	-373	-62.24	6
Rural settlements_Small	Rural settlements	Small	-19807	-109.43	181
Total			-303325	-219.46	144317

Portanto, o puco deficit de reserva legal da Caatinga está quase igualmente distribuído entre pequenas e grandes propriedades. Claramente, vejo duas oportunidades aqui. As grande propriedades podem focar em restauração com objetivos de biodiversidade. Jà os pequenos podem (têm esse direito) de fazer uso da RL na forma de SAFs

Deficit em APP segundo tipo e tamanho de propriedade

#Area total de APP
vsd_prop %>% 
  summarise_at(vars(area_app.1), funs(sum))

##   area_app.1
## 1     208370

vsd_prop %>% 
  group_by(landtenure2, propSize3) %>%
  filter(area_prop.MB.1>0) %>%
  summarise_at(vars(app_def), funs(sum,mean, n())) %>%
  unite(land_prop, landtenure2:propSize3, remove = FALSE) %>%
  #group_by(UF,land_prop) %>%
  #summarise_at(vars(vsd_total), sum) %>%
  #spread(land_prop, vsd_total) %>%
  #mutate_if(is.numeric, funs(replace_na(., 0))) %>%
  adorn_totals("row") %>%
  kable(caption ="Table 6. Deficit of natural vegetation between large and smallholders in APP") %>% 
  kable_paper(bootstrap_options = "striped", full_width = F)

Table 6. Deficit of natural vegetation between large and smallholders in APP
land_prop	landtenure2	propSize3	sum	mean	n
Private_Large	Private	Large	11959	0.567	21104
Private_Small	Private	Small	18433	0.019	975078
Quilombola_Small	Quilombola	Small	1090	15.566	70
Rural settlements_Small	Rural settlements	Small	7044	3.833	1838
Simulated_Large	Simulated	Large	56642	1.649	34349
Simulated_Small	Simulated	Small	24137	0.038	634516
Total			119306	21.671	1666955

Agora podemos, inclusive espacializar esse déficit na Caatinga

agora vem os mapas teste

Porcentagem de Caatinga média mas RLs por propriedade GRANDE por município

vsd_prop_map %>% 
  filter(propSize3 == "Large") %>% 
  group_by(code_muni) %>% 
  summarise_at(vars(RL_Perc08), funs(mean)) %>% 
  inner_join(semiarid, by = "code_muni")->RL_perc
  ggplot()+  geom_sf(data=RL_perc$geom, aes(fill=RL_perc$RL_Perc08), lwd = 0)+
  scale_fill_gradient(low = "red", high = "blue", na.value = NA)+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Mean % of forest cover in LR of Large proporties", size=8, fill = "Forest cover (%)") +
  theme_map()+
  theme(legend.position = "right")->fc_l
  fc_l

Porcentagem de Caatinga média por propriedade PEQUENA por município

vsd_prop_map %>% 
  filter(propSize3 == "Small") %>% 
  group_by(code_muni) %>% 
  summarise_at(vars(RL_Perc08), funs(mean)) %>% 
  inner_join(semiarid, by = "code_muni")->RL_perc
  ggplot()+  geom_sf(data=RL_perc$geom, aes(fill=RL_perc$RL_Perc08), lwd = 0)+
  scale_fill_gradient(low = "red", high = "blue", na.value = NA)+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Mean % of forest cover in LR of Small properties", size=8, fill = "Forest cover (%)") +
  theme_map()+
  theme(legend.position = "right")->fc_s
  fc_s

Qual a proporção de Caatinga por cataegoria de propriegade (G e P) por município

Aqui vemos que muitos municípios estão com a maior parte da Caatinga remanescente em grande propriedades

options(scipen=999)
vsd_prop_map %>% 
  group_by(code_muni) %>% 
  spread(propSize2, NV19.1) %>% 
  replace(is.na(.), 0) %>% 
  summarize_at(vars("Small","Large"), sum) %>% 
  mutate(LS_ratio = (Large/(Small+Large)*100)) %>% 
  inner_join(semiarid [,1:5], by = "code_muni")->LS_ratio
  ggplot()+  geom_sf(data=LS_ratio$geom, aes(fill=LS_ratio$LS_ratio), lwd = 0)+
  scale_fill_gradient(low = "orange", high = "red", na.value = NA)+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Share of forest in large properties", size=8, fill = "Forest in large\nproperties (%)") +
  theme_map()+
  theme(legend.position = "right")->share_fc
  share_fc

Agora temos o déficit total de RL em hectares (ha)

library(viridis)

## Loading required package: viridisLite

## 
## Attaching package: 'viridis'

## The following object is masked from 'package:scales':
## 
##     viridis_pal

options(scipen=999)
vsd_prop_map %>% 
  group_by(code_muni) %>% 
  filter(vsd_total<0) %>% 
  #spread(propSize2, vsd_total) %>% 
  #replace(is.na(.), 0) %>% 
  summarize_at(vars("vsd_total"), sum) %>% 
  #mutate(vsd_LR = (Large/(Small+Large)*100)) %>% 
  inner_join(semiarid [,1:5], by = "code_muni")->vsd_spatial
  ggplot()+  geom_sf(data=vsd_spatial$geom, aes(fill=vsd_spatial$vsd_total), lwd = 0)+
  scale_fill_gradient(low = "red", high = "orange", na.value = NA)+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Total deficit per municipality", size=8, fill = "Legal deficit (ha)") +
  theme_map()+
  theme(legend.position = "right")->deficit_tot
  deficit_tot

Agora os 10% de propriedades com mais deficit

options(scipen=999)
vsd_prop_map %>% 
  group_by(code_muni) %>% 
  filter(vsd_total<0) %>% 
  spread(propSize3, vsd_total) %>% 
  replace(is.na(.), 0) %>% 
  summarize_at(vars("Small","Large"), sum) %>% 
  mutate(vsd_LR = (Large/(Small+Large)*100)) %>% 
  inner_join(semiarid [,1:5], by = "code_muni")->vsd_spatial2
  ggplot()+  geom_sf(data=vsd_spatial2$geom, aes(fill=vsd_spatial2$vsd_LR), lwd = 0)+
  scale_fill_gradient(low = "orange", high = "red", na.value = NA)+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Share of deficit in large properties", size=8, fill = "Deficit in\nlarge properties (%)") +
  theme_map()+
  theme(legend.position = "right")->share_large
  share_large

options(scipen=999)
vsd_prop_map %>% 
  group_by(code_muni) %>% 
  filter(vsd_total<0) %>% 
  spread(propSize3, vsd_total) %>% 
  replace(is.na(.), 0) %>% 
  summarize_at(vars("Small","Large"), sum) %>% 
  mutate(vsd_LR = (Large/(Small+Large)*100)) %>% 
  inner_join(semiarid [,1:5], by = "code_muni")->vsd_spatial2
  ggplot()+  geom_sf(data=vsd_spatial2$geom, aes(fill=100-vsd_spatial2$vsd_LR), lwd = 0)+
  scale_fill_gradient(low = "orange", high = "red", na.value = "grey")+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Share of deficit in small properties", size=8, fill = "Deficit in\nsmall properties (%)") +
  theme_map()+
  theme(legend.position = "right")->share_small
  share_small

plot_grid(share_fc, share_large, share_small, deficit_tot)

Agora podemos resumir e quantificar onde estão os 1% de propriedades com maior déficit

# Tabela resumida
vsd_prop %>% 
  filter(vsd_total<(-100)) %>% 
  filter(area_prop.MB.1>0) %>% 
  #slice_max(-vsd_total, prop = 0.01) %>% 
  group_by(landtenure2, propSize3) %>% 
  summarise_at(vars(vsd_total), funs(sum, n())) %>% 
  rename(area_ha=sum, N_prop=n) %>% 
  mutate(per_def=area_ha/-132805)

## # A tibble: 3 x 5
## # Groups:   landtenure2 [3]
##   landtenure2       propSize3 area_ha N_prop per_def
##   <fct>             <fct>       <dbl>  <int>   <dbl>
## 1 Private           Large     -55247.    267 0.416  
## 2 Quilombola        Small       -281.      2 0.00212
## 3 Rural settlements Small     -15114.     57 0.114

# Mapa
vsd_prop_map %>% 
  filter(vsd_total<0) %>% 
  filter(area_prop.MB.1>0) %>% 
  slice_max(-vsd_total, prop = 0.01) %>% 
  inner_join(semiarid [,1:5], by = "code_muni") %>% 
  arrange(desc(vsd_total))->top_1_map
ggplot()+  geom_sf(data=top_1_map$geom, aes(fill=-top_1_map$vsd_total), lwd = 0)+
  scale_fill_gradientn(colours = colorspace::diverge_hcl(2))+
  geom_sf(data=uf_caat, fill="transparent")+
  geom_sf(data=caatinga, fill="transparent")+
  coord_sf(expand = FALSE)+
  labs(subtitle="Share of deficit in top 10% properties with highest deficit", size=8, fill = "Deficit of vegetation (ha)") +
  theme_map()+
  theme(legend.position = "right")->map_top_1
  map_top_1

RL_APP

Felipe Melo & Adriana Pallegrini

07/10/2020

Vamos carregar tudo que for necessário para as análises

Pacotes

Dados

Neste caso estou carregando dados de outro

Análises descritivas para propriedades privadas excluindo ARU e QL (Assentamentos e Terras Quilombolas)

Quanta área está em mãos de grandes e pequenos proprietários segundo o tipo de titulação ?

Quanta Caatinga (vegetação) existe em propriedades segundo seu tamanho?

Agora vamos ao deficit legal de APP e RL

Nestas análises vamos destrinchar por categoria de tamanho de prorpiedade. Dessa vez, os simulados vão ficar de fora, mas acredito que podemos pensar que o padrão encontrado para o resto da Caaringa, deve se repetir para eles.

Qual o deficit por propriedade da Caatinga?

Em forma de tabela, temos esses valores de déficit de RL

Deficit em APP segundo tipo e tamanho de propriedade

Agora podemos, inclusive espacializar esse déficit na Caatinga

agora vem os mapas teste

Porcentagem de Caatinga média mas RLs por propriedade GRANDE por município

Porcentagem de Caatinga média por propriedade PEQUENA por município

Qual a proporção de Caatinga por cataegoria de propriegade (G e P) por município

Agora temos o déficit total de RL em hectares (ha)

Agora os 10% de propriedades com mais deficit

Agora podemos resumir e quantificar onde estão os 1% de propriedades com maior déficit