Pisf Callaway Sant’Anna
Thomas Daher
15/03/2023
Install packages
Loading shapefile
municipalities.poly <- read_sf(dsn = "d:/Users/thomas.daher/Desktop/PISF/Municipios_beneficiados", layer = "Atlas_%C3%81gua_2021_-_PISF_-_Munic%C3%ADpios_Beneficiados")
channel.poly <- read_sf(dsn = "d:/Users/thomas.daher/Desktop/PISF/Canais_PISF", layer = "Atlas_%C3%81gua_2021_-_PISF_-_Canais")
adutoras.poly <- read_sf(dsn = "d:/Users/thomas.daher/Desktop/PISF/Adutoras_PISF", layer = "Atlas_%C3%81gua_2021_-_PISF_-_Adutoras")
bacia.poly <- read_sf(dsn = "d:/Users/thomas.daher/Desktop/PISF/Bacias_semi arido_nordeste1", layer = "GEOFT_BHO_ANO_TDR")
municipalities.df <- dplyr::select(as.data.frame(municipalities.poly), -geometry, -ben_pisf, -FID)extract data
### write.csv(municipalities.df, file = "d:/Users/thomas.daher/Desktop/PISF/Municipios_beneficiados/municipios_df.csv")Plot municipalities and channel with ANA definition
### sem audtoras
ggplot() +
geom_sf(data = municipalities.poly, aes(fill = municipalities.poly$ben_pisf), show.legend = TRUE, size = 0.1) +
geom_sf(data = channel.poly, colour = "blue", size = 2)## Warning: Use of `municipalities.poly$ben_pisf` is discouraged.
## i Use `ben_pisf` instead.
### com adutoras
ggplot() +
geom_sf(data = municipalities.poly, aes(fill = municipalities.poly$ben_pisf), show.legend = TRUE, size = 0.1) +
geom_sf(data = channel.poly, colour = "blue", size = 2) +
geom_sf(data = adutoras.poly, colour = "green", size = 0.5)## Warning: Use of `municipalities.poly$ben_pisf` is discouraged.
## i Use `ben_pisf` instead.
## Grupos de controle e tratamento
Queen weight matriz
list.queen<-poly2nb(municipalities.poly, queen=TRUE)
W<-nb2listw(list.queen, style="W", zero.policy=TRUE)Loading data
### grupo
grupo <- read_csv("d:/Users/thomas.daher/Desktop/PISF/Grupos - Grupos2.csv")## Rows: 762 Columns: 7
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (4): nm_municipio, uf, grupo, eixo
## dbl (3): cod_ibge, cod_grupo, cod_eixo
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.### pib_municipal
pib <- read.csv2("d:/Users/thomas.daher/Desktop/PISF/pesquisas_anuais/pib_municipal_serie_historica.csv")
colnames(pib)[4] <- "pib_municipal"
colnames(pib)[1] <- "cod_ibge"
pib <- na.omit(pib)
### PAM
pam <- read.csv2("d:/Users/thomas.daher/Desktop/PISF/pesquisas_anuais/Pam_serie_temporal.csv")
colnames(pam)[4] <- "area_plantada"
colnames(pam)[1] <- "cod_ibge"
pam <- na.omit(pam)important data cleaning
cod_ibge_geom <- municipalities.poly %>% select(cod_ibge,geometry)
grupo_geom <- merge(grupo, cod_ibge_geom, by.x = "cod_ibge",
all.x = TRUE)merging data
### pib
pib_grupo <- merge(pib, grupo, by.x = "cod_ibge",
all.x = TRUE)
pib_grupo <- pib_grupo %>% filter(!is.na(grupo)) %>% select(-Município)
### pam
pam <- pam %>% filter(cod_ibge != 1)
pam_grupo <- merge(pam, grupo, by.x = "cod_ibge", all.x = TRUE)
pam_grupo <- pam_grupo %>% filter(!is.na(grupo)) %>% select(-Brasil.e.Município)
pam_grupo$area_plantada <- as.numeric(pam_grupo$area_plantada)## Warning: NAs introduzidos por coerçãopam_grupo <- pam_grupo %>% filter(!is.na(area_plantada))
summary(pam_grupo)## cod_ibge Ano area_plantada nm_municipio
## Length:15132 Min. :2002 Min. : 1 Length:15132
## Class :character 1st Qu.:2006 1st Qu.: 833 Class :character
## Mode :character Median :2011 Median : 2303 Mode :character
## Mean :2011 Mean : 4704
## 3rd Qu.:2016 3rd Qu.: 5877
## Max. :2021 Max. :55340
##
## uf grupo cod_grupo eixo
## Length:15132 Length:15132 Min. :0.000 Length:15132
## Class :character Class :character 1st Qu.:2.000 Class :character
## Mode :character Mode :character Median :3.000 Mode :character
## Mean :2.858
## 3rd Qu.:4.000
## Max. :4.000
##
## cod_eixo
## Min. :0.000
## 1st Qu.:1.000
## Median :2.000
## Mean :2.486
## 3rd Qu.:4.000
## Max. :5.000
## NA's :14474Plot municipalities and channel with research definition
ggplot() +
geom_sf(data = grupo_geom, aes(fill = grupo_geom$grupo, geometry = geometry), show.legend = TRUE, size = 0.1) +
geom_sf(data = channel.poly, colour = "blue", size = 2)## Warning: Use of `grupo_geom$grupo` is discouraged.
## i Use `grupo` instead.
Descriptive statistics
pib mean by group
pib_grupo$pib_municipal <- as.numeric(pib_grupo$pib_municipal)
mean_pib <- pib_grupo %>% group_by(grupo, Ano) %>% summarise(pib_mean = mean(pib_municipal))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(mean_pib, mapping = aes(x = Ano, y = pib_mean, group = grupo)) + geom_line(aes(colour = grupo)) + ggtitle("média do pib municipal") + geom_vline(xintercept = 2008, linetype="dotted",
color = "blue", size=1) +
geom_vline(xintercept = 2012, linetype="dotted",
color = "red", size=1) +
geom_vline(xintercept = 2014, linetype="dotted",
color = "blue", size=1)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
### scatter plot of cities mean
mean_pib <- pib_grupo %>% group_by(grupo, nm_municipio) %>% summarise(pib_mean = mean(pib_municipal))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(mean_pib, mapping = aes(x = grupo, y = pib_mean, group = grupo)) + geom_point(aes(colour = grupo)) + ggtitle("média do pib municipal")
pib sum by group
pib_grupo$pib_municipal <- as.numeric(pib_grupo$pib_municipal)
sum_pib <- pib_grupo %>% group_by(grupo, Ano) %>% summarise(pib_sum = log(sum(pib_municipal)))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(sum_pib, mapping = aes(x = Ano, y = pib_sum, gorup = grupo)) + geom_line(aes(colour = grupo)) + ggtitle("log da soma do pib municipal") + geom_vline(xintercept = 2008, linetype="dotted",
color = "blue", size=1) +
geom_vline(xintercept = 2012, linetype="dotted",
color = "red", size=1) +
geom_vline(xintercept = 2014, linetype="dotted",
color = "blue", size=1)
pam mean by group
mean_area <- pam_grupo %>% group_by(grupo, Ano) %>% summarise(pam_mean = mean(area_plantada))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(mean_area, mapping = aes(x = Ano, y = pam_mean, gorup = grupo)) + geom_line(aes(colour = grupo)) + ggtitle("média da área plantada") + geom_vline(xintercept = 2008, linetype="dotted",
color = "blue", size=1) +
geom_vline(xintercept = 2012, linetype="dotted",
color = "red", size=1) +
geom_vline(xintercept = 2014, linetype="dotted",
color = "blue", size=1)
### scatter plot of cities mean
mean_area <- pam_grupo %>% group_by(grupo, nm_municipio) %>% summarise(pam_mean = mean(area_plantada))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(mean_area, mapping = aes(x = grupo, y = pam_mean, group = grupo)) + geom_point(aes(colour = grupo)) + ggtitle("média da área plantada")
pam log sum by group
sum_area <- pam_grupo %>% group_by(grupo, Ano) %>% summarise(pam_sum = log(sum(area_plantada)))## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.ggplot(sum_area, mapping = aes(x = Ano, y = pam_sum, gorup = grupo)) + geom_line(aes(colour = grupo)) + ggtitle("log da soma da área plantada") + geom_vline(xintercept = 2008, linetype="dotted",
color = "blue", size=1) +
geom_vline(xintercept = 2012, linetype="dotted",
color = "red", size=1) +
geom_vline(xintercept = 2014, linetype="dotted",
color = "blue", size=1)