Mapas temáticos com R e ggplot2
Versão: 05/08/2014
Este tutorial objetiva ilustrar como uilizar o R para apresentar dados quantitativos em mapas. Mais especificamente será ilustrado como produzir mapas temáticos.
Presume-se que o usuário tenha alguma familiaridade com o R.
Especificamente será ilustrado neste pequeno tutorial como apresentar dados relativos a alguns indicadores municipais do rol de diretrizes, objetivos, metas e indicadores 2013 - 2015 do Estado do Rio de Janeiro colhidos junto ao TABNET do DATASUS.
Uma vez colhidos os dados, no formato csv, será necessário obter o mapa no qual desejamos apresentar os dados. Trabalharemos com mapas no formato shapefiles, que podem ser baixados no site do IBGE.
De posse dos dados e dos gráficos o passo seguinte é certificar-se de que os seguintes pacotes estejam instalados no R: ggplot2, rgdal e RColorBrewer.
Após exportar os dados necessários no formato csv, a etapa seguinte é importar estes dados para o R. Isto pode ser feito da seguinte forma:
Parte 1 - Importação dos dados
# Definir o diretório de trabalho
setwd("E:\\Mapas com R\\dados")
## SIFILIS CONGÊNITA
##------------------------------
# leitura dos dados
sifilis <- read.csv2('Casos de Sífilis Congênita.csv', header=FALSE, na.string=c('-', '...'))
# filtro para obter apenas as linhas de interesse
sifilis <- subset(sifilis, grepl('^33', V1))
# renomeando as colunas
names(sifilis) <- c('Municipio', 'sifilis2010', 'sifilis2011', 'sifilis2012', 'sifilis2013')
# apresentação dos registros iniciais para checar se tudo ocorreu bem.
head(sifilis)## Municipio sifilis2010 sifilis2011 sifilis2012
## 5 330010 Angra dos Reis 1 1 2
## 6 330015 Aperibé 1 0 0
## 7 330020 Araruama 1 2 2
## 8 330022 Areal 0 0 0
## 9 330023 Armação dos Búzios 0 0 0
## 10 330025 Arraial do Cabo 0 0 0
## sifilis2013
## 5 7
## 6 NA
## 7 1
## 8 NA
## 9 2
## 10 1Procedimento idêntico será feito para cada uma dos demais indicadores:
setwd("E:\\Mapas com R\\dados")
## COBERTURA DE VACINA
##------------------------------
cobvacina <- read.csv2("cobertura de vacina.csv", header=FALSE, na.string='...')
cobvacina <- subset(cobvacina, grepl('^33', V1))
names(cobvacina) <- c('Municipio', 'cobvacina2010', 'cobvacina2011', 'cobvacina2012')
## COBERTURA EQUIPES ATENÇÃO BÁSICA
##---------------------------------
cobeqab <- read.csv2("cobertura Equipes de Atenção Básica.csv", header=FALSE, na.string='...')
cobeqab <- subset(cobeqab, grepl('^33', V1))
names(cobeqab) <- c('Municipio', 'cobeqab2010', 'cobeqab2011', 'cobeqab2012', 'cobeqab2013')
## COBERTURA EQUIPES SAÚDE BUCAL
##---------------------------------
cobeqsb <- read.csv2("cobertura Equipes Saúde Bucal.csv", header=FALSE, na.string='...')
cobeqsb <- subset(cobeqsb, grepl('^33', V1))
names(cobeqsb) <- c('Municipio', 'cobeqsb2010', 'cobeqsb2011', 'cobeqsb2012', 'cobeqsb2013')
## CURA HANSENÍASE
##---------------------------------
hansen <- read.csv2("Cura Hanseníase.csv", header=FALSE, na.string=c('-', '...'))
hansen <- subset(hansen, grepl('^33', V1))
names(hansen) <- c('Municipio', 'hansen2010', 'hansen2011', 'hansen2012', 'hansen2013')
## CURA TUBERCULOSE
##---------------------------------
tuberc <- read.csv2("Cura Tuberculose.csv", header=FALSE, na.string=c('-', '...'))
tuberc <- subset(tuberc, grepl('^33', V1))
tuberc <- tuberc[, c('V1', 'V2', 'V3', 'V4', 'V5')]
names(tuberc) <- c('Municipio', 'tuberc2010', 'tuberc2011', 'tuberc2012', 'tuberc2013')
## INTERNAÇÕES POR CONDIÇÕES SENSÍVEIS ATENÇÃO BÁSICA
##-----------------------------------------------------
iconsen <- read.csv2("ISABs.csv", header=FALSE, na.string='...')
iconsen <- subset(iconsen, grepl('^33', V1))
names(iconsen) <- c('Municipio', 'iconsen2010', 'iconsen2011', 'iconsen2012', 'iconsen2013')
## MÉDIA ESCOVAÇÃO DENTAL SUPERVISIONADA
##-----------------------------------------------------
escden <- read.csv2("Media de escovação supervisionada.csv", header=FALSE, na.string='...')
escden <- subset(escden, grepl('^33', V1))
escden <- escden[, c('V1', 'V2', 'V3', 'V4', 'V5')]
names(escden) <- c('Municipio', 'escden2010', 'escden2011', 'escden2012', 'escden2013')
## NASCIDOS VIDOS COM 7 OU MAIS CONSULTAS PRÉ NATAL
##---------------------------------------------------
nascvivo <- read.csv2("Nascidos Vivos de 7 ou mais consultas.csv", header=FALSE, na.string='...')
nascvivo <- subset(nascvivo, grepl('^33', V1))
nascvivo <- nascvivo[, c('V1', 'V2', 'V3', 'V4', 'V5')]
names(nascvivo) <- c('Municipio', 'nascvivo2010', 'nascvivo2011', 'nascvivo2012', 'nascvivo2013')
## EXODONTIAS
##---------------------------------------------------
exodont <- read.csv2("Proporção de exodontia.csv", header=FALSE, na.string=c('-', '...'))
exodont <- subset(exodont, grepl('^33', V1))
names(exodont) <- c('Municipio', 'exodont2010', 'exodont2011', 'exodont2012', 'exodont2013')Parte 2 - Preparo da base de dados
Feita a importação das bases de dados, o passo seguinte é combinar todas as beses em uma base única. Isso pode ser feito da seguinte forma:
## COMBINAR TODAS AS BASES DE DADOS
##---------------------------------------------------
indic <- list(sifilis, cobvacina, cobeqab, cobeqsb, escden, hansen,iconsen, nascvivo, tuberc, exodont)
# combinar todas as bases de dados
base <- Reduce(merge, indic)
# extrair o código do ibge da variável Municipio
base <- transform(base, ibge = gsub('^(\\d{6}).*', '\\1', base$Municipio))A base de dados precisa conter informação acerca de quais municípios integram uma determinada região de saúde. Será necessário incluir esta informação na base dados. Para tanto, criaresmos uma variável denominada regiao, que indicará a qual região de saude pertence o município.
O código a seguir tem esta finalidade;
##DEFINIR AS REGIÕES DE SAÚDE
##---------------------------
# Passo 1 - criar uma lista definindo os municípios que compõem cada região de saúde:
regioes <- list(
BXL = c('330130', '330452', '330070', '330023', '330520', '330025', '330187', '330020', '330550'),
BIG = c('330260', '330010', '330380'),
CSU = c('330540', '330600', '330095', '330022', '330370', '330385', '330290', '330620', '330180', '330280', '330360'),
MPB = c('330450', '330610', '330030', '330630', '330395', '330400', '330440', '330040', '330411', '330412', '330420', '330225'),
MT1 = c('330250', '330170', '330045', '330510', '330320', '330455', '330285', '330414', '330227', '330350', '330555', '330200'),
MT2 = c('330560', '330430', '330575', '330190', '330270', '330490', '330330'),
NRO = c('330410', '330615', '330310', '330060', '330220', '330230', '330205', '330115', '330090', '330513', '330300', '330470', '330015', '330210'),
NTE = c('330475', '330500', '330100', '330480', '330415', '330093', '330140', '330240'),
SER = c('330390', '330185', '330080', '330580', '330515', '330340', '330570', '330050', '330160', '330120', '330150', '330245', '330590', '330110', '330530', '330460')
)
# Passo 2 - converter a lista em um data frame
regioes <- data.frame(regiao = names(unlist(regioes)), ibge = unname(unlist(regioes)))
regioes$regiao <- substr(regioes$regiao, 1, 3)
# Passo 3 - juntar com a base de indicadores
base <- merge(base, regioes, by='ibge')
base$ibge <- as.character(base$ibge)Agora que criamos uma base de dados completa, podemos remover o que não mais nos interessa.
## Remover objetos desnecessários
rm(indic, tuberc, sifilis, regioes, nascvivo, iconsen, hansen, exodont, escden, cobvacina, cobeqsb, cobeqab)Agora que já temos uma base de dados contendo os valores dos indicadores para alguns anos, vamos contruir uma nova base de dados contendo a média dos valores de cada indicador, nos anos para os quais possuimos informação. Esta nova base de dados pode ser construída da seguinte forma:
# Criar uma base de dados com a média dos valores das variáveis nos anos
#-----------------------------------------------------------------------
base.media <- data.frame(
ibge = base$ibge,
regiao = base$regiao,
sifilis = apply(base[, grepl('^sifilis', names(base))], 1, mean, na.rm=TRUE),
cobvacina = apply(base[, grepl('^cobvacina', names(base))], 1, mean, na.rm=TRUE),
cobeqab = apply(base[, grepl('^cobeqab', names(base))], 1, mean, na.rm=TRUE),
cobeqsb = apply(base[, grepl('^cobeqsb', names(base))], 1, mean, na.rm=TRUE),
escden = apply(base[, grepl('^escden', names(base))], 1, mean, na.rm=TRUE),
hansen = apply(base[, grepl('^hansen', names(base))], 1, mean, na.rm=TRUE),
iconsen = apply(base[, grepl('^iconsen', names(base))], 1, mean, na.rm=TRUE),
nascvivo = apply(base[, grepl('^nascvivo', names(base))], 1, mean, na.rm=TRUE),
tuberc = apply(base[, grepl('^tuberc', names(base))], 1, mean, na.rm=TRUE),
exodont = apply(base[, grepl('^exodont', names(base))], 1, mean, na.rm=TRUE)
)
head(base.media)## ibge regiao sifilis cobvacina cobeqab cobeqsb escden hansen iconsen
## 1 330010 BIG 2.7500 82.51 94.04 76.70 0.7275 95.22 22.80
## 2 330015 NRO 0.3333 118.68 100.00 100.00 1.1375 66.67 29.02
## 3 330020 BXL 1.5000 93.75 32.33 26.54 0.0000 62.68 25.30
## 4 330022 CSU 0.0000 103.28 100.00 100.00 0.8100 72.97 26.84
## 5 330023 BXL 0.5000 82.37 100.00 100.00 0.0000 100.00 20.38
## 6 330025 BXL 0.2500 91.39 78.38 67.73 0.5275 75.00 20.04
## nascvivo tuberc exodont
## 1 77.59 67.41 6.310
## 2 71.23 75.00 6.885
## 3 43.68 78.45 8.578
## 4 73.33 100.00 2.210
## 5 63.49 84.38 14.660
## 6 56.96 72.08 6.815Agora podemos plotar estes valores médios no mapa.
Parte 3 - Plotagem dos dados
Bem, agora que possuímos os dados a serem plotados no mapa, o passo seguinte é importar o mapa para o R, juntar a nossa base de dados ao mapa e plotar. Vamos ir passo a passo.
No site do IBGE, os mapas estão em um arquivo zipado. Será necessário descompactar os arquivos. A importação do mapa será feita utilizando o pacoge rgdal, conforme ilustrado no código a seguir:
setwd('E:\\Mapas com R')
library(ggplot2)
library(rgdal)## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
## Path to GDAL shared files: C:/Users/Marcos/Documents/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/Marcos/Documents/R/win-library/3.1/rgdal/projlibrary(RColorBrewer)
# Leitura do arquivo shapefile
rj <- readOGR("mapas rj", "33MUE250GC_SIR") # Municipios## OGR data source with driver: ESRI Shapefile
## Source: "mapas rj", layer: "33MUE250GC_SIR"
## with 92 features and 3 fields
## Feature type: wkbPolygon with 2 dimensionsFeita a importação dos dados, podemos visualizar a base de dados georeferenciada da seguinte forma:
head(rj@data)## ID CD_GEOCODM NM_MUNICIP
## 0 1468 3300100 ANGRA DOS REIS
## 1 1469 3300159 APERIBÉ
## 2 1470 3300209 ARARUAMA
## 3 1471 3300225 AREAL
## 4 1472 3300233 ARMAÇÃO DOS BÚZIOS
## 5 1473 3300258 ARRAIAL DO CABOPodemos visualizar diversas variáveis, mas a que vai nos interessar mais é a variável CD_GEOCODM, que contém o código do município. Olhando os valores da variável é possível ver que o código tem um dígito a mais, do que o código do muncípio contida na base da dados que acabamos de prepara. Será necessário um pequeno ajuste com vistas à compatibilizar os valores. O cádigo a seguir faz isso:
# Ajuste no código dos municípios
rj$CD_GEOCODM <- substr(rj$CD_GEOCODM, 1, 6)Feito este pequeno ajuste, é possível agora juntar nossa base de dados aos dados georeferenciados contidos no mapa. Isso pode ser feito da seguinte forma:
## Juntar ao mapa os dados a serem plotados
rj <- merge(rj, base.media, by.x='CD_GEOCODM', by.y='ibge') Para usarmos o pacote ggplot2 para plotar os dados no gráfico, será necessário extrair da base geoferenciada um data frame com as informações a serem plotadas. Para fazer esta extração utiliza-se a função fortify() do pacote ggplot2.
# Extrai um data frame com coordenadas - variáveis: long lat order hole piece group id
# gpclibPermitStatus(TRUE)
rj.f <- fortify(rj, region = "CD_GEOCODM") ## Loading required package: rgeos
## rgeos version: 0.3-6, (SVN revision 450)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUErj.f <- merge(rj.f, rj@data, by.x = "id", by.y = "CD_GEOCODM")
# inspeçao da base de dados
str(rj.f)## 'data.frame': 181143 obs. of 20 variables:
## $ id : chr "330010" "330010" "330010" "330010" ...
## $ long : num -44.4 -44.4 -44.4 -44.4 -44.3 ...
## $ lat : num -22.9 -22.9 -22.8 -22.8 -22.8 ...
## $ order : int 1 2 3 4 5 6 7 8 9 10 ...
## $ hole : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ piece : Factor w/ 106 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ group : Factor w/ 340 levels "330010.1","330010.2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ID : num 1468 1468 1468 1468 1468 ...
## $ NM_MUNICIP: Factor w/ 92 levels "ANGRA DOS REIS",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ regiao : Factor w/ 9 levels "BIG","BXL","CSU",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sifilis : num 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 ...
## $ cobvacina : num 82.5 82.5 82.5 82.5 82.5 ...
## $ cobeqab : num 94 94 94 94 94 ...
## $ cobeqsb : num 76.7 76.7 76.7 76.7 76.7 76.7 76.7 76.7 76.7 76.7 ...
## $ escden : num 0.728 0.728 0.728 0.728 0.728 ...
## $ hansen : num 95.2 95.2 95.2 95.2 95.2 ...
## $ iconsen : num 22.8 22.8 22.8 22.8 22.8 ...
## $ nascvivo : num 77.6 77.6 77.6 77.6 77.6 ...
## $ tuberc : num 67.4 67.4 67.4 67.4 67.4 ...
## $ exodont : num 6.31 6.31 6.31 6.31 6.31 6.31 6.31 6.31 6.31 6.31 ...Agora podemos partir para a apresentação dos dados no mapa. Vamos ilustrar o procedimento com a variável exodont. Para tanto vamos inicialmente categorizar a variável. Vamos criar uma nova variável exodontCateg que assumirá os valores ‘0-5%’, ‘5-10%’, ‘10-15%’, ‘15-20%’, ‘20-25%’, ‘+25’. Esta categorização pode ser feita da seguinte forma:
# Categorização da variável exodont
rj.f$exodontCateg <- cut(rj.f$exodont, breaks = c(0, 5, 10, 15, 20, 25, Inf),
labels=c('0-5%', '5-10%', '10-15%', '15-20%', '20-25%', '+25'), include.lowest=TRUE)Feita a categorizão, o mapa pode ser obtido da seguinte forma:
ggplot(rj.f, aes(long, lat, group = group, fill = exodontCateg)) +
geom_polygon(colour='black') +
coord_equal() +
ggtitle('Exodontias') +
#theme(plot.title=element_text(size=rel(1), lineheight=.9, face="bold", colour="blue")) +
labs(x = "", y = "", fill = "%") +
scale_fill_manual(values=brewer.pal(9, 'Greens')[5:9])
Támbem é possível visualizar a variável por regiões de saúde. O código a seguir ilustra como é fazer isso:
ggplot(rj.f, aes(long, lat, group = group, fill = exodontCateg)) +
geom_polygon(colour='black') +
coord_equal() +
ggtitle('Exodontias por Região de Saude') +
#theme(plot.title=element_text(size=rel(1), lineheight=.9, face="bold", colour="blue")) +
labs(x = "", y = "", fill = "%") +
facet_wrap( ~ regiao, nrow=3, scales='free') +
scale_fill_manual(values=brewer.pal(9, 'Greens')[5:9])
Leituras adicionais:
http://metodologiapolitica.com/como-gerar-mapas/#sthash.PLzgFHxS.dpbs
http://sociaisemetodos.wordpress.com/2013/09/06/mapas-e-dados-georreferenciados-no-r-parte-1/
http://sociaisemetodos.wordpress.com/2013/09/15/mapas-no-r-parte-2-utilizando-shapes/