Gráficos Miojo: mapas temáticos, the tidy way

Resposta

Para quem não tem paciência de ler o post, a resposta está aqui, em 4 linhas de código! Carregue o tidyverse e o brazilmaps e rode o código abaixo.

library(tidyverse)
library(brazilmaps)

get_brmap("City", geo.filter = list(State = 33)) %>% 
  left_join(pop2017, c("City" = "mun")) %>% 
  ggplot() +
  geom_sf(aes(fill = pop2017/1e6))

No caso, estamos fazendo um mapa do Rio de Janeiro (Unidade Federativa de código 33) e colorindo o mapa dos municípios com a população das cidades em 2017 (pop2017, dividida por um milhão). A cidade em azul claro é a capital.

OK, mas como fizemos isso, e como podemos criar mapas mais bonitos? É o que veremos a seguir!

Objetivo

Nosso objetivo é fazer mapas temáticos, que são mapas que levam alguma característica das regiões a cores dessas regiões. O nome específico desse gráfico é mapa coroplético (choropleth map). Vamos fazer dois deles: um de unidades federativas e um de municípios.

Para fazer um mapa desse tipo, precisamos de dois insumos principais: i) os dados que vão compor as cores de cada região e ii) os polígonos para desenhar o mapa.

Geralmente esses dois insumos vêm de bases de dados diferentes, então é importante que exista uma chave para ligá-las através de um join. Se você obteve seus dados de uma fonte oficial, provavelmente você terá o código do município ou da unidade federativa. Se não, provavelmente você terá de cruzar pelos nomes, o que pode ser uma fonte de dores de cabeça.

Para reproduzir as análises que seguem, carregue os pacotes tidyverse e brazilmaps:

library(tidyverse)
library(brazilmaps)

Polígonos

A melhor fonte oficial para obter os polígonos do Brasil é o Instituto Brasileiro de Geografia e Estatística (IBGE). Lá, você vai encontrar malhas territoriais em um formato chamado shapefile. Esse arquivo, que na verdade é um conjunto de arquivos com dados e metadados, contém as informações necessárias para desenhar os polígonos no seu mapa.

Por muito tempo, a leitura desse tipo de arquivo no R foi uma dor de cabeça. Até 2015 pelo menos, eu utilizava o sp e o maptools, e sempre sofria um pouco para carregar a base e fazer mapas.

Agora, com o pacote sf, isso ficou bem mais fácil. Atendendo aos princípios tidy, o sf armazena mapas em data.frames, onde cada cada linha é um território, cada coluna representa uma característica desse território, e a lista-coluna geometry contém as informações dos polígonos que representam esse território.

Como objetos de classe sf são data.frames, estamos em casa. Todas operações dos amados pacotes dplyr e ggplot2 já estão preparadas para lidar com esses objetos, logo a tarefa de trabalhar com mapas se resume à tarefa de transformação e visualização de dados, do jeito que gostamos.

E como faz para ler esses dados no R? Com o pacote brazilmaps, basta utilizar a função get_brmap(). Essa função recebe parâmetros do tipo de região que se deseja baixar, e também algum filtro que se deseja aplicar. Por exemplo, se você quiser carregar a base no nível Região, utilize

library(brazilmaps)
mapa <- get_brmap("Region")
glimpse(mapa)

Note que a base de dados contém apenas cinco linhas, uma para cada região, e a coluna especial geometry. Se você deseja obter todos os municípios do estado de São Paulo, utilize

mapa <- get_brmap("State", geo.filter = list(State = 35))
glimpse(mapa)

O pacote brazilmaps permite que você obtenha as malhas em seis níveis diferentes: Brasil, Regiões, Estados, Municípios, Mesorregiões e Microrregiões.

Dados

Para preencher os dados, vamos utilizar a base pnud_muni do pacote abjData. Esse pacote foi desenvolvido pela Associação Brasileira de Jurimetria (ABJ) para guardar e organizar dados úteis para a Jurimetria. O pacote não está disponível no CRAN, então para utilizá-lo você precisará instalar o pacote do GitHub com

devtools::install_github("abjur/abjData")

A base de dados pnud_muni, em particular, contém as informações utilizadas para calcular o Índice de Desenvolvimento Humano (IDH) a nível municipal, com base nos dados dos Censos de 1991, 2000 e 2010. É uma base muito rica e pode ser utilizada para realizar diversas análises interessantes.

No caso, utilizaremos apenas as colunas identificadoras de município, unidade federativa, região e IDH, e apenas o ano de 2010. Faremos isso rodando:

pnud_minima <- abjData::pnud_muni %>% 
  filter(ano == 2010) %>% 
  select(cod_uf = uf, 
         cod_muni = codmun7, 
         starts_with("idh"))

glimpse(pnud_minima)

Essa base de dados contém somente 5565 municípios, que eram os que existiam na época que os dados foram levantados. Atualmente, temos 5570 municípios no Brasil.

Juntando as bases

Vamos primeiro calcular o IDH médio por unidade federativa, depois juntar com a base de estados. Rigorosamente falando, nós não deveríamos fazer isso, pois a média dos IDHs municipais de uma unidade federativa não é o IDH da unidade federativa. Mas vida que segue:

pnud_minima_estado <- pnud_minima %>% 
  group_by(cod_uf) %>% 
  summarise(idh = mean(idhm))

glimpse(pnud_minima_estado)

Agora vamos juntar com o mapa de unidades federativas:

uf_map <- get_brmap("State") %>% 
  inner_join(pnud_minima_estado, c("State" = "cod_uf"))

glimpse(uf_map)

Parece que funcionou! Agora, temos uma coluna idh. Agora vamos fazer um mapa. Para isso, usaremos uma operação geométrica especial do ggplot2: o geom_sf(). Essa função permite desenhar o mapa diretamente, sem a necessidade de realizar transformações nos dados nem mapear longitude e latitude.

uf_map %>% 
  ggplot() +
  geom_sf(aes(fill = idh))

O gráfico mostra claramente o aspecto regional do IDH nos estados do Brasil.

Para fazer o gráfico com os municípios, a lógica é a mesma: i) fazer um join das bases com dplyr e ii) usar o ggplot2 para desenhar o mapa. O único problema é que o gráfico demora para rodar¹.

muni_map <- get_brmap("City") %>% 
  left_join(pnud_minima, c("City" = "cod_muni"))

muni_map %>% 
  ggplot() +
  geom_sf(aes(fill = idhm))

Podemos também quebrar pelas regiões com um simples facet_wrap():

muni_map %>% 
  ggplot() +
  geom_sf(aes(fill = idhm)) +
  facet_wrap(~Region)

Parece que estamos indo bem. Alguns problemas que podemos notar:

1. As linhas estão muito largas.
2. As cores ajudam a identificar tendências regionais, mas parece que a imagem está embaçada.
3. Os mapas não estão centralizados em cada facet. E não adianta usar o parâmetro scales="free", pois infelizmente o sistema de coordenadas do sf não suporta essa funcionalidade.

Resolveremos esses problemas na próxima seção.

Masterizando os mapas

muni_map %>% 
  ggplot() +
  geom_sf(aes(fill = idhm), 
          # ajusta tamanho das linhas
          colour = "black", size = 0.1) +
  # muda escala de cores
  scale_fill_viridis_c(option = 2) +
  # tira sistema cartesiano
  theme(panel.grid = element_line(colour = "transparent"),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

muni_map %>% 
  ggplot() +
  geom_sf(aes(fill = idhm), 
          # ajusta tamanho das linhas
          colour = "black", size = 0.1) +
  geom_sf(data = get_brmap("State"),
          fill = "transparent",
          colour = "black", size = 0.5) +
  # muda escala de cores
  scale_fill_viridis_c(option = 2) +
  # tira sistema cartesiano
  theme(panel.grid = element_line(colour = "transparent"),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

muni_map %>% 
  mutate(idhm = cut(idhm, c(0, 0.5, 0.6, 0.75, 1.0))) %>% 
  ggplot() +
  geom_sf(aes(fill = idhm), 
          # ajusta tamanho das linhas
          colour = "black", size = 0.1) +
  geom_sf(data = get_brmap("State"),
          fill = "transparent",
          colour = "black", size = 0.5) +
  # muda escala de cores
  scale_fill_viridis_d(option = 2, begin = 0.2, end = 0.8) +
  # tira sistema cartesiano
  theme(panel.grid = element_line(colour = "transparent"),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

grafico_regiao <- function(regiao) {
  # titulo do mapa
  titulo <- switch(regiao, "1" = "Norte", "2" = "Nordeste",
                   "3" = "Sudeste", "4" = "Sul", "5" = "Centro-Oeste")
  muni_map %>% 
    filter(Region == regiao) %>% 
    mutate(idhm = cut(idhm, c(0, 0.5, 0.6, 0.75, 1.0))) %>% 
    ggplot() +
    geom_sf(aes(fill = idhm), 
            # ajusta tamanho das linhas
            colour = "black", size = 0.1) +
    geom_sf(data = filter(get_brmap("State"), Region == regiao),
            fill = "transparent",
            colour = "black", size = 0.5) +
    # muda escala de cores
    scale_fill_viridis_d(option = 2, begin = 0.2, end = 0.8) +
    guides(fill = FALSE) +
    # tira sistema cartesiano
    theme(panel.grid = element_line(colour = "transparent"),
          panel.background = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank()) +
    ggtitle(titulo)
}
mapas <- map(1:5, grafico_regiao)

library(patchwork)
reduce(mapas, magrittr::add)

Wrap-up

1. Use brazilmaps para carregar os mapas.
2. Faça join com a sua base de interesse.
3. Use geom_sf() do ggplot2 para plotar os mapas.
4. Customize seus mapas com as opções do ggplot2 e com o pacote patchwork!

Se quiserem ver uma aplicação legal disso, dêem uma olhada nessa análise da Bruna Wundervald elaborado durante a Datathon do RLadies-SP.

É isso pessoal. Happy coding ;)