Mapas raciais de pontos no R com o pacote sf

Da obtenção dos dados até os mapas

Mapas raciais de pontos no R com o pacote sf

Da obtenção dos dados até os mapas

O trabalho com mapas raciais feito por Hugo Nicolau e divulgado em seu blog me chamou a atenção para a metodologia conhecida como mapas de pontos de contagem, potente para a análise da densidade de fenômenos.

Ela é usada em uma porção de áreas como epidemiologia, demografia, logística, biologia e tem grande espaço para aplicações no campo das políticas públicas, pois é potente para abordar uma série de problemas, entre eles associações de desigualdade.

Recentemente busquei formas de reproduzir o trabalho do Hugo no R e, fortunamente, encontrei que é na verdade um procedimento relativamente simples. Neste post vou compartilhar como construir mapas de ponto de contagens no R em base tidy, criando mapas da cidade de Peruíbe, no litoral sul de São Paulo.

Para isto utilizarei os pacotes sf para manipular os dados espaciais e o ggplot2 para gerar os mapas. Além disso, farei uso também de funções criadas por mim, agrupadas por hora no pacote seda, para obter os dados do IBGE.

Mapas de pontos de contagem

Criar um mapa de pontos é simples quando você já tem o layer de pontos. O que esta metodologia propõe é recriar representações espaciais dadas por polígonos em representações espaciais dadas por pontos, baseadas nas informações associadas aos polígonos.

A ideia é simples: cada ponto representará uma quantidade n de observações (ex.: um ponto para cada indivíduo, para cada 10 etc) e eles serão distribuídos aleatoriamente dentro de polígonos espaciais aos quais estão associados.

Exemplo. Quando faz o Censo o IBGE registra dados demográficos e socioeconômicos de toda a população brasileira e cada pessoa é registrada no setor censitário de sua residência. O setores censitários são pequenas porções territoriais contíguas utlizadas para fins cadastrais e constituem a unidade básica de desagregação espacial do Censo.

Com estes dados construir um mapa de pontos significa distribuir aleatoriamente pontos representando pessoas em seus respectivos setores censitários.

Passo 1: obter os dados

A ideia é criar um mapa que permita analisar a associação espacial entre raça e renda dentro dos limites da cidade de Peruíbe.

O procedimento envolve especificamente variáveis de raça e renda dos resultados do universo do Censo 2010. Estes dados estão disponíveis no ftp do IBGE e no seu geoftp.

Os dados demográficos fazem parte dos resultados do universo e os dados espaciais, isto é, os polígonos dos setores censitários, integram o conjunto de dados sobre malhas territoriais.

Obter dados do IBGE, embora fácil, é uma tarefa dispendiosa. São muitos arquivos, é muito pesado para guardar tudo (além de desnecessário já que o IBGE já faz o serviço de armazená-los) e quem trabalha com estes dados vira e mexe precisa novamente baixá-los.

Para facilitar esta tarefa eu criei algumas funções que facilitam a obtenção dos dados do IBGE diretamente no R, particularmente os dos resultados do universo e os dados espaciais. Elas fazem parte de um pacote chamado seda. Para construir os mapas não é necessário instalá-lo, basta ter os dados e você pode obtê-los nos links acima. Mas se você preferir seguir pelo caminho mais fácil:

# para instalar o pacote seda
devtools::install_github("bruno-pinheiro/seda")

Obs.: o pacote seda foi testado apenas em linux Debian 9

# carregar pacotes
library(seda)
library(sf)
library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)

Quero a malha de setores censitários e as seguintes variáveis populacionais:

  • número de habitantes nos setores censitários por raça
  • número de habitantes nos setores censitários por classe social

Como a análise envolve renda, trabalharei apenas com dados a respeito de pessoas com 10 anos ou mais, já que esta é a faixa etária a partir do qual o IBGE passa a coletar informações de rendimento.

Usando o pacote seda vou baixar e importar o shapefile de setores censitários do estado de São Paulo. Esta etapa demorará um pouco, pois os dados são pesados. Dependendo da sua conexão, dá tempo de pegar (ou tomar) um café.

# malha de setores
url <- "ftp://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2010/setores_censitarios_shp/sp/sp_setores_censitarios.zip"
get_ibge(url)

# dados tabulares
url <- "ftp://ftp.ibge.gov.br/Censos/Censo_Demografico_2010/Resultados_do_Universo/Agregados_por_Setores_Censitarios/SP_Exceto_a_Capital_20190207.zip"
get_ibge(url)

rm(url)

Os arquivos foram baixados na pasta temporária e preciso saber seus nomes para importar no R.

Vou então buscar o nome do shapefile…

#listar arquivo .shp
list_tmp("shp")
## [1] "/tmp/RtmpxdkmxU/35SEE250GC_SIR.shp"

e listar os arquivos em formato .csv que foram baixados.

# listar arquivos .csv
list_tmp("csv")
##  [1] "/tmp/RtmpxdkmxU/Basico_SP2.csv"          
##  [2] "/tmp/RtmpxdkmxU/Domicilio01_SP2.csv"     
##  [3] "/tmp/RtmpxdkmxU/Domicilio02_SP2.csv"     
##  [4] "/tmp/RtmpxdkmxU/DomicilioRenda_SP2.csv"  
##  [5] "/tmp/RtmpxdkmxU/Entorno01_SP2.csv"       
##  [6] "/tmp/RtmpxdkmxU/Entorno02_SP2.csv"       
##  [7] "/tmp/RtmpxdkmxU/Entorno03_SP2.csv"       
##  [8] "/tmp/RtmpxdkmxU/Entorno04_SP2.csv"       
##  [9] "/tmp/RtmpxdkmxU/Entorno05_SP2.csv"       
## [10] "/tmp/RtmpxdkmxU/Pessoa01_SP.csv"         
## [11] "/tmp/RtmpxdkmxU/Pessoa02_SP.csv"         
## [12] "/tmp/RtmpxdkmxU/Pessoa03_SP.csv"         
## [13] "/tmp/RtmpxdkmxU/Pessoa04_SP.csv"         
## [14] "/tmp/RtmpxdkmxU/Pessoa05_SP.csv"         
## [15] "/tmp/RtmpxdkmxU/Pessoa06_SP.csv"         
## [16] "/tmp/RtmpxdkmxU/Pessoa07_SP.csv"         
## [17] "/tmp/RtmpxdkmxU/Pessoa08_SP.csv"         
## [18] "/tmp/RtmpxdkmxU/Pessoa09_SP.csv"         
## [19] "/tmp/RtmpxdkmxU/Pessoa10_SP.csv"         
## [20] "/tmp/RtmpxdkmxU/Pessoa11_SP.csv"         
## [21] "/tmp/RtmpxdkmxU/Pessoa12_SP.csv"         
## [22] "/tmp/RtmpxdkmxU/Pessoa13_SP.csv"         
## [23] "/tmp/RtmpxdkmxU/PessoaRenda_SP2.csv"     
## [24] "/tmp/RtmpxdkmxU/responsavel01_sp2.csv"   
## [25] "/tmp/RtmpxdkmxU/responsavel02_sp2.csv"   
## [26] "/tmp/RtmpxdkmxU/ResponsavelRenda_SP2.csv"

Pronto. O shapefile está nomeado como “35SEE250GC_SIR.shp” e as variáveis tabulares de interesse estão nos seguintes arquivos (sim, é preciso ler o dicionário de variáveis):

  • Pessoa03_SP.csv
    • raça: V002:V006
  • Pessoa11_SP.csv e Pessoa12_SP.csv:
    • sexo: V001
  • PessoaRenda_SP.csv:
    • renda: V001 a V009 (excluindo pessoas sem rendimento)

Então já posso importar os dados:

# importar dados em lista
camadas <- list(
    # raca
  raca = read_ibge("Pessoa03_SP.csv") %>% select(Cod_setor, V002:V006),
  # renda
  renda = read_ibge("PessoaRenda_SP2.csv") %>% select(Cod_setor, V001:V009),
  # setores censitários
  setores = read_geoibge("35SEE250GC_SIR")
  )
camadas
## $raca
## # A tibble: 47,733 x 6
##    Cod_setor V002  V003  V004  V005  V006 
##        <dbl> <chr> <chr> <chr> <chr> <chr>
##  1   3.50e14 249   3     28    24    0    
##  2   3.50e14 550   13    86    136   0    
##  3   3.50e14 432   3     65    24    0    
##  4   3.50e14 613   8     63    111   0    
##  5   3.50e14 744   4     31    33    1    
##  6   3.50e14 740   4     52    71    0    
##  7   3.50e14 438   13    15    96    0    
##  8   3.50e14 656   17    41    112   1    
##  9   3.50e14 424   20    7     113   0    
## 10   3.50e14 465   11    2     256   1    
## # … with 47,723 more rows
## 
## $renda
## # A tibble: 47,733 x 10
##    Cod_setor V001  V002  V003  V004  V005  V006  V007  V008  V009 
##        <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
##  1   3.50e14 0     43    63    33    47    34    4     3     4    
##  2   3.50e14 10    135   168   78    83    56    10    9     2    
##  3   3.50e14 0     57    94    45    79    61    16    22    3    
##  4   3.50e14 8     122   125   78    83    88    16    10    12   
##  5   3.50e14 5     133   179   74    98    55    5     8     7    
##  6   3.50e14 11    147   175   90    84    65    8     9     3    
##  7   3.50e14 15    159   118   46    24    22    1     2     2    
##  8   3.50e14 26    152   172   90    87    42    4     2     1    
##  9   3.50e14 10    109   156   41    49    13    2     3     0    
## 10   3.50e14 14    160   217   48    20    15    2     1     0    
## # … with 47,723 more rows
## 
## $setores
## Simple feature collection with 68296 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -335491.8 ymin: 7196857 xmax: 586149.1 ymax: 7803609
## epsg (SRID):    31983
## proj4string:    +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 68,296 x 14
##       ID CD_GEOCODI TIPO  CD_GEOCODB NM_BAIRRO CD_GEOCODS NM_SUBDIST
##    <dbl> <chr>      <chr> <chr>      <chr>     <chr>      <chr>     
##  1 67103 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  2 67104 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  3 67105 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  4 67106 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  5 67107 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  6 67108 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  7 67109 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  8 67110 350280405… URBA… <NA>       <NA>      350280405… <NA>      
##  9 67111 350280405… URBA… <NA>       <NA>      350280405… <NA>      
## 10 67112 350280405… URBA… <NA>       <NA>      350280405… <NA>      
## # … with 68,286 more rows, and 7 more variables: CD_GEOCODD <chr>,
## #   NM_DISTRIT <chr>, CD_GEOCODM <chr>, NM_MUNICIP <chr>, NM_MICRO <chr>,
## #   NM_MESO <chr>, geometry <MULTIPOLYGON [m]>

Como resultado cada uma das três bases é um elemento dentro da lista camadas.

O shapefile contem os polígonos dos setores censitários e as duas bases tabulates contém informações sobre raça e renda dos habitantes dos 68296 setores censitários do estado de São Paulo, embora eu queira apenas os de Peruíbe.

O próximo passo é preparar os dados.

Passo 2: preparar os dados

Farei o filtro e a limpeza das observações e variáveis que não serão utilizadas, além dos tratamentos necessários. Ao fim quero ter as seguintes variáveis na base de setores censitários:

  • raca: brancos, pretos, amarelos, pardos, indigenas
  • renda: classeA, classeB, classeC, classeD, classeE

Com relação à renda, vale destacar que as faixas de renda são baseadas em salários mínimos e que em 2010 o salário mínimo era de R$ 510. Especificamente, as faixas de renda serão:

  • classeA: > 10 salários mínimos
  • classeB: 5 a 10 salários mínimos
  • classeC: 2 a 5 salários mínimos
  • classeD: 1 a 2 salários mínimos
  • classeE: até 1 salário mínimo
# filtrar setores censitários de Peruíbe
camadas$setores <- camadas$setores %>% 
  filter(NM_MUNICIP == "PERUÍBE") %>% 
  select(CD_GEOCODI)

# manipular variáveis de raça e renda
camadas$raca <- camadas$raca %>%  
  filter(Cod_setor %in% camadas$setores$CD_GEOCODI) %>% 
  mutate(Cod_setor = as.character(Cod_setor)) %>% 
  mutate_at(vars(-Cod_setor), as.numeric) %>% 
  rename(CD_GEOCODI = Cod_setor, brancos = V002, pretos = V003,
         amarelos = V004, pardos = V005, indigenas = V006)

camadas$renda <- camadas$renda %>%  
      filter(Cod_setor %in% camadas$setores$CD_GEOCODI) %>% 
      mutate(Cod_setor = as.character(Cod_setor)) %>% 
      mutate_at(vars(-Cod_setor), as.numeric) %>% 
  mutate(classeE = V001 + V002,
         classeC = V004 + V005,
         classeA = V007 + V008 + V009) %>% 
  select(CD_GEOCODI = Cod_setor, classeA, classeB = V006,
         classeC, classeD = V003, classeE)

# mesclar dados de raça e renda na malha de setores
camadas$setores <- camadas$setores %>% 
  select(CD_GEOCODI) %>% 
  left_join(camadas$raca, by = "CD_GEOCODI") %>% 
  left_join(camadas$renda, by = "CD_GEOCODI") %>%
  arrange(CD_GEOCODI)

# checar dados
camadas$setores
## Simple feature collection with 130 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 284871.1 ymin: 7294947 xmax: 307339.1 ymax: 7330827
## epsg (SRID):    31983
## proj4string:    +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 130 x 12
##    CD_GEOCODI brancos pretos amarelos pardos indigenas classeA classeB
##    <chr>        <dbl>  <dbl>    <dbl>  <dbl>     <dbl>   <dbl>   <dbl>
##  1 353760205…     424      5       31     33         0      24      51
##  2 353760205…     154      5        2     41         0       4      23
##  3 353760205…     125      1       33      7         0       4       8
##  4 353760205…     333     13       17     82         0       6      27
##  5 353760205…     267      9       27     86         0       9      21
##  6 353760205…     268     26       15    110         0       1       8
##  7 353760205…     257     13       10     31         0       6      17
##  8 353760205…     338     12       17     81         0      10      35
##  9 353760205…     448     10       23     87         0      12      47
## 10 353760205…     434      4        8     41         0      10      43
## # … with 120 more rows, and 4 more variables: classeC <dbl>,
## #   classeD <dbl>, classeE <dbl>, geometry <MULTIPOLYGON [m]>

Aí está a variável de identificação do setor censitário reunida às 10 variáveis que indicam o número de pessoas de cada grupo populacional em cada setor censitário.

Passo 3: Processamento espacial

Entendendo o procedimento

Uma forma bastante simples de realizar a atribuição espacial aleatória de pontos no R é com a função st_sample() do pacote sf. O que a função faz é pegar um argumento numérico associado a um polígono, que especifique o número de pontos a serem criados, e atribuir aleatoriamente o número determinado de pontos dentro deste polígono.

Este número determinado de pontos \(n\) é dado a partir de uma razão \(b\) do número de habitantes dos grupos populacionais nos setores censitários \(k_{mi}\), em que \(k\) é o número de habitantes, \(m\) é o grupo populacional e \(i\) é o setor censitário. Ou seja, é simples como:

\[n_{mi} = \frac{k_{mi}}{b}\]

São mais de 60 mil habitantes, o que geraria uma sobreposição desnecessária no mapa. Por isso vamos reduzir esta densidade em uma razão \(b = 5\), isto é, um ponto para cada 5 pessoas.

Atribuição aleatória dos pontos

Para realizar a atribuição aleatória dos pontos é preciso tratar os valores ausentes, pois a função st_sample() precisa do número de pontos a serem atribuídos para cada setor. Se não há o número de pessoas indíginas em um setor, por exemplo, a função falhará.

Além de distribuir aleatoriamente os pontos nos setores censitários, eu quero:

  • realizar a itersecção entre os pontos e os polígonos, para identificar em qual setor cada ponto está
  • criar uma variável indicando o grupo populacional que os pontos representam

E por fim reagruparei os pontos representando a distribuição espacial da população por raça e aqueles que representam a distribuição por renda.

É necessário repetir este procedimento para cada um dos grupos, isto é, 10 vezes. Aplicando a função lapply() o processo é feito em looping, os objetos são estocando numa lista e evita-se repetição de código.

# separar dados missing
nas <- camadas$setores[is.na(camadas$setores$brancos), ]
camadas$setores <- camadas$setores[!is.na(camadas$setores$brancos), ]

# separar nomes das variáveis
v <- names(as_tibble(camadas$setores)[, 2:11])
grupos <- c("Brancos", "Pretos", "Amarelos", "Pardos", "Indígenas",
            "Classe A", "Classe B", "Classe C", "Classe D", "Classe E")

# realizar a atribuição dos pontos para cada grupo populacional
pontos <- lapply(1:10, function(i) {
    st_sample(camadas$setores, round(camadas$setores[, v[i]][[v[i]]]) / 5) %>% 
      st_sf() %>% 
      st_join(camadas$setores %>% select(CD_GEOCODI), join = st_intersects) %>%
      mutate(categoria = grupos[i])
})

# unir grupos de faixa e renda
# e reunir setores censitários
camadas <- list(
  raca = do.call(rbind, pontos[(1:5)]),
  renda = do.call(rbind, pontos[6:10]),
  setores = rbind(camadas$setores, nas)
  )

rm(nas, pontos, v)

camadas
## $raca
## Simple feature collection with 12086 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 285443.3 ymin: 7295754 xmax: 307151.4 ymax: 7330695
## epsg (SRID):    31983
## proj4string:    +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## First 10 features:
##         CD_GEOCODI categoria                 geometry
## 1  353760205000001   Brancos POINT (297360.1 7308998)
## 2  353760205000001   Brancos POINT (297393.3 7308714)
## 3  353760205000001   Brancos POINT (297253.3 7308865)
## 4  353760205000001   Brancos POINT (297167.7 7308978)
## 5  353760205000001   Brancos POINT (297314.5 7308839)
## 6  353760205000001   Brancos   POINT (297317 7309213)
## 7  353760205000001   Brancos POINT (297354.5 7308958)
## 8  353760205000001   Brancos POINT (297652.6 7309018)
## 9  353760205000001   Brancos POINT (296996.9 7308825)
## 10 353760205000001   Brancos POINT (297408.1 7309041)
## 
## $renda
## Simple feature collection with 6181 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 285515.7 ymin: 7296513 xmax: 306920.8 ymax: 7330286
## epsg (SRID):    31983
## proj4string:    +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## First 10 features:
##         CD_GEOCODI categoria                 geometry
## 1  353760205000001  Classe A POINT (297156.5 7308927)
## 2  353760205000001  Classe A POINT (297319.6 7309036)
## 3  353760205000005  Classe A POINT (296755.2 7308705)
## 4  353760205000007  Classe A POINT (296853.4 7309423)
## 5  353760205000008  Classe A POINT (297394.4 7309465)
## 6  353760205000008  Classe A POINT (297426.5 7309398)
## 7  353760205000009  Classe A POINT (297282.1 7310019)
## 8  353760205000009  Classe A POINT (297373.4 7309870)
## 9  353760205000009  Classe A   POINT (297687 7309750)
## 10 353760205000009  Classe A POINT (297624.6 7309765)
## 
## $setores
## Simple feature collection with 130 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 284871.1 ymin: 7294947 xmax: 307339.1 ymax: 7330827
## epsg (SRID):    31983
## proj4string:    +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 130 x 12
##    CD_GEOCODI brancos pretos amarelos pardos indigenas classeA classeB
##    <chr>        <dbl>  <dbl>    <dbl>  <dbl>     <dbl>   <dbl>   <dbl>
##  1 353760205…     424      5       31     33         0      24      51
##  2 353760205…     154      5        2     41         0       4      23
##  3 353760205…     125      1       33      7         0       4       8
##  4 353760205…     333     13       17     82         0       6      27
##  5 353760205…     267      9       27     86         0       9      21
##  6 353760205…     268     26       15    110         0       1       8
##  7 353760205…     257     13       10     31         0       6      17
##  8 353760205…     338     12       17     81         0      10      35
##  9 353760205…     448     10       23     87         0      12      47
## 10 353760205…     434      4        8     41         0      10      43
## # … with 120 more rows, and 4 more variables: classeC <dbl>,
## #   classeD <dbl>, classeE <dbl>, geometry <MULTIPOLYGON [m]>

Os dados agora estão prontos para ser visualizados.

Passo 4: Criar os mapas

Mapas de raças

pal <- list(
  raca = c("yellow2", "dodgerblue", "red", "orange4", "black"),
  renda = rev(c("#ffffb2","#fecc5c","#fd8d3c","#f03b20","#bd0026"))
  )

leg <- c("Raça", "Classe social")

p <- lapply(1:2, function(i) {
  ggplot() +
    geom_sf(data = st_geometry(camadas$setores),
            col = "grey30", alpha = 0, lwd = .1) +
    geom_sf(data = camadas[[i]], aes(colour = categoria),
            size = .3, show.legend = "point") +
    scale_colour_manual(values = pal[[i]], leg[i]) +
    guides(color = guide_legend(override.aes = list(size=2, alpha = 1))) +
    ggtitle(paste("Distribuição espacial dos habitantes por", leg[i])) +
    light_map_theme()
})

do.call(grid.arrange, c(p, ncol = 2))

Mapas dos grupos populacionais

p <- lapply(1:2, function(i)
    ggplot() +
      geom_sf(data = camadas$setores,
              color = "grey30", lwd = .06, alpha = 0) +
      geom_sf(data = camadas[[i]], aes(colour = categoria),
              size = 1.2, alpha = 0.08) +
      scale_colour_manual(values = pal[[i]], guide = FALSE) +
      ggtitle(leg[i]) +
      facet_wrap(~categoria, ncol = 5) +
      light_map_theme()
    )

do.call(gridExtra::grid.arrange, c(p, nrow = 2))

Depois de ver esses mapas, o que você pensa sobre a relação entre renda, raça e local de moradia na cidade de Peruíbe?


Veja também