Tiago Mendonça - Mapas com R

Atualização: como o pacote brazilmaps não é compatível com as verões mais recentes do R, atualizei os códigos com o pacote geobr

Uma excelente referência para o assunto tratado aqui é o livro Geocomputation with R.

Utilizaremos arquivos vetorizados para trabalhar com mapas. Um arquivo vetorizado pode ser definido de acordo com geometrias que indicam como o mapa é formado. As geometrias utilizadas podem ser dadas pelas seguintes categorias:

POINT: um único ponto,
LINESTRING: uma sequência de pontos conectados por uma linha reta e
POLYGON: um polígono fechado (pode conter buracos internamente)

Cada geometria pode ser composta por uma série de elementos pertencentes a uma mesma geometria. Nesse caso, teríamos mulitpoint, multilinestring e multipolygon. Além disso, podemos ter um arquivo com uma combinação de diferentes geometrias formando uma geometrycollection.

vamos carregar um mapa com o pacote brazilmaps e verificar a classe do objeto com o comando class().

library(ggplot2)
library(dplyr)
library(viridis)
# library(brazilmaps) 
library(geobr)
library(sf)
library(maptools)
library(leaflet)

theme_set(theme_bw())

# mapa <- brazilmaps::get_brmap("State")

mapa <- read_state(showProgress = FALSE)

class(mapa)

[1] "sf"         "data.table" "data.frame"

Note que o objeto mapa pertence às classes sf (simple feature) e data.frame. A classe sf faz referência a propriedades do mapa e o data frame apresenta informações associadas a cada uma das geometrias do mapa.

A seguir vamos verificar algumas propriedades importantes desse objeto. Incialmente temos as características do mapa (tipo de geometria, epsg e proj4string) e, posteriormente, uma tabela com os dados associados a cada geometria.

mapa

Simple feature collection with 27 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -66.81025 ymin: -13.6937 xmax: -59.77435 ymax: -7.969294
Geodetic CRS:  SIRGAS 2000
First 10 features:
   code_state abbrev_state name_state code_region name_region
1          11           RO   Rondônia           1       Norte
2          12           AC       Acre           1       Norte
3          13           AM   Amazonas           1       Norte
4          14           RR    Roraima           1       Norte
5          15           PA       Pará           1       Norte
6          16           AP      Amapá           1       Norte
7          17           TO  Tocantins           1       Norte
8          21           MA   Maranhão           2    Nordeste
9          22           PI      Piauí           2    Nordeste
10         23           CE      Ceará           2    Nordeste
                             geom
1  MULTIPOLYGON (((-63.32721 -...
2  MULTIPOLYGON (((-73.18253 -...
3  MULTIPOLYGON (((-67.32609 2...
4  MULTIPOLYGON (((-60.20051 5...
5  MULTIPOLYGON (((-54.95431 2...
6  MULTIPOLYGON (((-51.1797 4....
7  MULTIPOLYGON (((-48.35878 -...
8  MULTIPOLYGON (((-45.84073 -...
9  MULTIPOLYGON (((-41.74605 -...
10 MULTIPOLYGON (((-41.16703 -...

Portanto, para fazer um gráfico de acordo com alguma medida, podemos associá-la a tabela e, assim, estará diretamente associada a geometria.

Primeiro, vamos fazer um gráfico apenas com as geometrias.

ggplot(mapa)+ 
  geom_sf()

Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação ( fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.

acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51, 50, 31, 15, 
                                   25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), 
                         com_rede = c(0.273, 0.412, 0.313, 0.177, 0.513, 0.696, 1.000, 0.974, 0.280, 0.065, 
                                      0.191, 0.449, 0.916, 0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 
                                      0.405, 0.096, 0.400, 0.352, 0.998, 0.347, 0.129))

mapa %>% 
  left_join(acesso_san, by = "code_state") %>% 
  mutate(com_rede = 100*com_rede) %>%  
  ggplot(aes(fill = com_rede), color = "black") +
    geom_sf() + 
    scale_fill_viridis(name = "Municípios com rede de esgoto (%)", direction = -1)

Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF). Primeiro vamos gerar um objeto com as coordenadas dos centroides e depois adicionar os pontos ao mapa.

coord_pontos <- mapa %>% 
                  left_join(acesso_san, by = "code_state") %>% 
                  mutate(com_rede = 100*com_rede) %>% 
                  st_centroid()

ggplot(mapa)+ 
  geom_sf() + 
  geom_sf(data = coord_pontos, aes(size = com_rede), col = "blue", alpha = .65,
          show.legend = "point") + 
  scale_size_continuous(name = "Municípios com rede de esgoto (%)")

Leaflet

Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet. Com base nesse formato, retomaremos o exemplo utilizado anteriormente da porcentagem de municípios com rede de esgoto de acordo com a UF e adicionaremos círculos com os raios proporcionais a essa porcentagem nos centroides das respectivas UF. Utilizaremos o comando st_coordinates() para obter um data frame com as coordenadas desses centroides.

data.frame(st_coordinates(coord_pontos), 
           com_rede = coord_pontos$com_rede, 
           UF = coord_pontos$name_state) %>% 
  leaflet() %>% 
    addTiles() %>%
    addCircleMarkers(~ X, ~ Y,
                     label = ~ as.character(paste0(UF, ": ", com_rede, "%")),
                     labelOptions = labelOptions(textsize = "13px"),
                     radius = ~ com_rede/10,
                     fillOpacity = 0.5)

É possível criar, utilizando poucas linhas, um mapa apenas com a latitude e longitude dos pontos. A seguir faremos um gráfico marcando três estações do metrô em São Paulo.

estacoes <- data.frame(estacao = c("Saúde", "Santa Cruz", "Paraíso"),
                       lat = c(-23.6185, -23.5989, -23.5765),
                       long = c(-46.6393, -46.6366, -46.6408))

leaflet(estacoes) %>% 
  addTiles() %>%
  addMarkers(~long, ~lat, label = ~as.character(estacao))

É possível calcular a distância entre os pontos com a função st_distance. Para uma introdução sobre projeções e sistemas de coordenadas, veja esse vídeo

coordinates(estacoes) <- c('long', 'lat')
proj4string(estacoes) <- CRS("+init=epsg:4326")

distance <- estacoes %>% 
              st_as_sf() %>% 
              st_distance 

dimnames(distance) <- list(estacoes$estacao, estacoes$estacao)
distance

Units: [m]
              Saúde Santa Cruz  Paraíso
Saúde         0.000   2196.718 4672.695
Santa Cruz 2196.718      0.000 2527.275
Paraíso    4672.695   2527.275    0.000

Caso tenha alguma crítica, sugestão ou comentário, me envie uma mensagem.