Este documento demonstra como importar, processar, analisar e visualizar dados de material particulado PM2.5 em municípios brasileiros, comparando com os limites estabelecidos pela OMS e CONAMA.
# --------------------------------------------------------
# 1️⃣ Instalação (executar apenas uma vez - descomente se necessário)
# --------------------------------------------------------
# install.packages("dplyr")
# install.packages("lubridate")
# install.packages("ggplot2")
# install.packages("sf")
# install.packages("geobr")
# --------------------------------------------------------
# 2️⃣ Carregar bibliotecas
# --------------------------------------------------------
library(dplyr)
library(lubridate)
library(ggplot2)
library(sf)
library(geobr) # Para baixar malha municipal # Pacotes necessários
library(dplyr)
library(lubridate)
library(ggplot2)
library(sf)
library(geobr) # Para baixar malha municipal
#---------------------------------------------------
# 1. Carregar dados
#---------------------------------------------------
agg <- read.csv("dados_dia2.csv")
# Base OMS
oms <- agg %>%
filter(uf %in% c(13)) %>%
mutate(
ano_mes = ymd(ano_mes),
year = year(ano_mes),
month = month(ano_mes),
oms = med > 15 # Limite OMS (24h): 15 µg/m³
) %>%
filter(oms == T) %>%
group_by(cod6, month, uf, oms) %>%
summarise(n = n(), .groups = "drop")
# Base conama
conama <- agg %>%
filter(uf %in% c(13)) %>%
mutate(
ano_mes = ymd(ano_mes),
year = year(ano_mes),
month = month(ano_mes),
conama = med > 50, # Limite CONAMA (24h): 50 µg/m³
) %>%
filter(conama == T) %>%
group_by(cod6, month, uf, conama) %>%
summarise(n = n(), .groups = "drop")
#---------------------------------------------------
# 2. Classificar número de dias acima do limite OMS
#---------------------------------------------------
conama <- conama %>%
mutate(class = cut(n,
breaks = c(-Inf, 5, 10, 15, 20, 25, Inf),
labels = c("0 a 05 dias", "05 a 10 dias", "10 a 15 dias",
"15 a 20 dias", "20 a 25 dias", "acima de 25 dias"),
right = FALSE,
include.lowest = TRUE))
oms <- oms %>%
mutate(class = cut(n,
breaks = c(-Inf, 5, 10, 15, 20, 25, Inf),
labels = c("0 a 05 dias", "05 a 10 dias", "10 a 15 dias",
"15 a 20 dias", "20 a 25 dias", "acima de 25 dias"),
right = FALSE,
include.lowest = TRUE))
#---------------------------------------------------
# 3. Carregar malha municipal do Brasil (6 dígitos IBGE)
#---------------------------------------------------
mun <- read_municipality(year = 2020) %>%
mutate(code_muni6 = as.numeric(substr(code_muni, 1, 6))) %>%
filter(code_state == 13)## Using year/date 2020##
|
| | 0%
|
|=== | 4%
|
|===== | 7%
|
|======== | 11%
|
|========== | 15%
|
|============= | 19%
|
|================ | 22%
|
|================== | 26%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================= | 41%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================= | 59%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|==================================================== | 74%
|
|====================================================== | 78%
|
|========================================================= | 81%
|
|============================================================ | 85%
|
|============================================================== | 89%
|
|================================================================= | 93%
|
|=================================================================== | 96%
|
|======================================================================| 100%ufs <- read_state(year = 2020) %>%
filter(code_state == 13)## Using year/date 2020##
|
| | 0%
|
|=== | 4%
|
|===== | 7%
|
|======== | 11%
|
|========== | 15%
|
|============= | 19%
|
|================ | 22%
|
|================== | 26%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================= | 41%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================= | 59%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|==================================================== | 74%
|
|====================================================== | 78%
|
|========================================================= | 81%
|
|============================================================ | 85%
|
|============================================================== | 89%
|
|================================================================= | 93%
|
|=================================================================== | 96%
|
|======================================================================| 100%#---------------------------------------------------
# 4. Juntar com os dados
#---------------------------------------------------
mapa1 <- mun %>%
left_join(conama, by = c("code_muni6" = "cod6")) %>%
filter(!is.na(class))
mapa2 <- mun %>%
left_join(oms, by = c("code_muni6" = "cod6")) %>%
filter(!is.na(class))
#---------------------------------------------------
# 5. Gerar o mapa
#---------------------------------------------------
mapa_pm1<-ggplot() +
geom_sf(data = mapa1, aes(fill = class), color = NA) +
geom_sf(data = ufs, fill = NA, color = "black", size = .15) +
scale_fill_manual(
values = c("yellowgreen", "greenyellow", "yellow", "gold", "darkorange", "red"),
name = "Dias acima do limite"
) +
facet_wrap(~ month, ncol = 4) +
theme_bw(base_size = 8) +
theme(legend.position = "bottom") +
labs(
title = "Material particulado PM 2.5 segundo classificação CONAMA",
subtitle = "Elaborado por Observatório de Clima e Saúde - LIS/ICICT/Fiocruz",
caption = 'pm 2,5: Média de 24 horas: 50 µg/m3'
)
ggsave("mapapm1.png", plot = mapa_pm1, width = 10, height = 7, dpi = 300)
mapa_pm2<-ggplot() +
geom_sf(data = mapa2, aes(fill = class), color = NA) +
geom_sf(data = ufs, fill = NA, color = "black", size = .15) +
scale_fill_manual(
values = c("yellowgreen", "greenyellow", "yellow", "gold", "darkorange", "red"),
name = "Dias acima do limite"
) +
facet_wrap(~ month, ncol = 4) +
theme_bw(base_size = 8) +
theme(legend.position = "bottom") +
labs(
title = "Material particulado PM 2.5 segundo classificação CONAMA",
subtitle = "Elaborado por Observatório de Clima e Saúde - LIS/ICICT/Fiocruz",
caption = 'pm 2,5: média de 24 horas: 15 µg/m3'
)
ggsave("mapapm2.png", plot = mapa_pm2, width = 10, height = 7, dpi = 300)