Mapas no R, parte 2: utilizando shapes

brasilTrabalhar com mapas ou dados georreferenciados frequentemente significa fazer uso de shapes. Em linhas gerais, shapes são de arquivos de mapas formados por polígonos geocodificados (ie., suas coordenadas são latitudes e longitudes) que delineiam os contornos de unidades espaciais de interesse, como limites municipais, estaduais etc. Basicamente o que se pretende nesses casos é “colorir” os polígonos de acordo com alguma característica ou atributo relevante. Este é o tipo de mapa que Leonardo utilizou em seus posts sobre médicos no Brasil e sobre transportes na cidade de SP.

Trabalhar com mapas como polígonos significa manusear (ao menos) três arquivos conjuntamente: um que possui a extensão .shp (que traz o desenho do mapa propriamente), outro com extensão .shx (que é um indexador de informações, para facilitar buscas) e um terceiro com extensão .dbf (que é um arquivo de banco de dados, trazendo informações e atributos das unidades espaciais – p.ex.: população do município, renda per capita etc.) O que liga esses três arquivos é um ID, único para cada unidade espacial, geralmente chamado de Geocódigo. Há uma explicação mínima sobre shapefiles na Wikipédia em português e uma explicação mais completa na versão em inglês.

O objetivo desse post é replicar este mapa simples sobre PIBs municipais, disponível no site da Wikipédia. O exemplo é de Diogo Ferrari.

Precisamos então seguir os seguintes passos:

  1. Fazer download de um shape de estados brasileiros
  2. Fazer um webscrapping para obter as informações dispostas na tabela do Wiki (e “limpar” tais informações)
  3. Definir aspectos da formatação do mapa (cores, rótulos da legenda)
  4. Parear os dados da tabela e as informações do arquivo de atributos (.dbf)
  5. Plotar

Primeiramente, carregamos os pacotes necessários:

require(XML)
require(RCurl)
require(maptools)
require(RColorBrewer)

Em seguida, para guardar os arquivos do mapa que serão baixados, criamos uma pasta e a definimos como diretório de trabalho (esse passo é opcional):

dir.create("c:/mapas/")
setwd("c:/mapas/")

Agora utilizaremos a função download.file para baixar um arquivo compactado (.zip) que contém o mapa (.shp, .shx e .dbf).  Isto está disponível no site Gismaps. Em seguida aplicamos a função unzip para extrair o conteúdo compactado. Por fim, para ler o mapa e guardá-lo num objeto do R, aplicamos a função readShapePoly do pacote maptools. Lendo o arquivo .shp já estamos implicitamente fazendo referência ao .shx e ao .dbf.

# Carregando dados do Mapa estadual
download.file("https://dl.dropboxusercontent.com/u/35964572/estados_2010.zip",destfile="estados_2010.zip")

unzip("estados_2010.zip")

mapaUF = readShapePoly("estados_2010.shp")

Nesse ponto, o mapa já é “plotável”:

plot(mapaUF)

plot_zoom

Mas só conseguimos ver seus formatos… não há cores que indiquem quaisquer informações sobre cada polígono (unidade da federação). O arquivo de atributos (.dbf) não contém ainda dados interessantes. Por isso, vamos ao passo 2: fazer o download das informações sobre PIBs estaduais.

Vamos usar a função readHTMLTable, já apresentada num post anterior.

PIBEstadual = readHTMLTable(getURL("https://pt.wikipedia.org/wiki/Anexo:Lista_de_estados_do_Brasil_por_PIB", ssl.verifypeer = FALSE),which=2)

Como praticamente toda informação baixada na internet, os dados da tabela chegam “sujos”  ou desformatados no R. É preciso fazer alguns procedimentos antes de utilizá-la:

# Mantemos apenas as colunas 3 e 4, que contém o nome da UF e seu PIB
PIBEstadual = PIBEstadual[2:nrow(PIBEstadual),3:4]

#renomeando as colunas
names(PIBEstadual)=c("UF","PIB")

#retira o símbolo '.' , que estava dividindo os milhares
PIBEstadual$PIB = gsub("\\.", "",PIBEstadual$PIB)

#substitui o nome pela sigla das UFs
PIBEstadual$UF = c('SP', 'RJ', 'MG', 'RS', 'PR', 'BA', 'SC', 'DF', 'GO', 'PE', 'ES', 'CE', 'PA', 'AM', 'MT', 'MA', 'MS', 'RN', 'PB', 'AL', 'SE', 'RO', 'PI', 'TO', 'AC', 'AP', 'RR')

# informa que os dados contidos na coluna 'PIB' são números
PIBEstadual$PIB = as.numeric(PIBEstadual$PIB)

# Divide o valor do PIB por 10.000, para que seja expresso em bilhões de reais
# (já está em R$ 1.000)
PIBEstadual$PIB = PIBEstadual$PIB/10^6

# Transforma os dados do PIB em uma variável categórica.
PIBEstadual$PIB_cat = cut(PIBEstadual$PIB, breaks=c(0,30,60,90,120,500,5000),
 labels=c('até 30','+ 30', '+ 60', '+ 90', '+ 120', '+ 500'))

Passamos ao passo 3: definir a formatação.

# Selecionamos algumas cores de uma paleta de cores do pacote RColorBrewer
paletaDeCores = brewer.pal(9, 'OrRd')
paletaDeCores = paletaDeCores[-c(3,6,8)]

# Agora fazemos um pareamento entre as faixas da variável sobre PIB (categórica) e as cores:
coresDasCategorias = data.frame(PIB_cat=levels(PIBEstadual$PIB_cat), Cores=paletaDeCores)
PIBEstadual = merge(PIBEstadual, coresDasCategorias)

Agora o passo 4: fazer o pareamento entre os dados da tabela e o mapa:

# Primeiramente, guardamos os dados do .dbf num objeto do R.
# Ele é um atributo do objeto mapaUF
mapaData = attr(mapaUF, 'data')

# Guardamos o número das linhas numa nova variável
# Esse passo é necessário pois o pareamento entre esses dados e a tabela do PIB
# muda a ordem os casos, o que prejudica, depois, a construção do mapa
mapaData$Index = row.names(mapaData)

# Mudando o nome da variável que indica a sigla dos estados
names(mapaData)[3] = "UF"

# Fundimos então as duas tabelas:
mapaData = merge(mapaData, PIBEstadual, by="UF")

# Reordenamos os dados do mapa
mapaData = mapaData[order(as.numeric(mapaData$Index)),]

# E guardamos essas informações novamente como sendo um atributo do arquivo de mapa.
attr(mapaUF, 'data') = mapaData

Agora, o passo final, plotar:


# Configurando tela (reduzindo as margens da figura)
parDefault = par(no.readonly = T)
layout(matrix(c(1,2),nrow=2),widths= c(1,1), heights=c(4,1))
par (mar=c(0,0,0,0))

# Plotando mapa
plot(mapaUF, col=as.character(mapaData$Cores))
plot(1,1,pch=NA, axes=F)
legend(x='center', legend=rev(levels(mapaData$PIB_cat)),
 box.lty=0, fill=rev(paletaDeCores),cex=.8, ncol=2,
 title='Mapa dos estados brasileiros segundo o PIB em\n2010. Em bilhões de reais:')

Pronto:

mapa.pib.estados

Anúncios

Mapas e dados georreferenciados no R – parte 1

Existem diversos softwares especializados para trabalhar com dados geográficos, que possuem funções das mais diversas, dedicadas desde desenhar e/ou corrigir polígonos, passando por gerenciamento de bancos dados até executar complexos modelos de estatística espacial. No entanto, pouco a pouco, a comunidade de usuários do R vem implementando pacotes, funções e interfaces para lidar com sistemas de informação geográficas (GIS). Talvez ainda hoje não possibilitam todos os usos que um software como o ArcGIS. Mas conhecendo um pouco da grande velocidade e gana com que se desenvolvem pacotes para o R, é bem provável daqui a pouco isso seja revertido (se é que é verdade!).

Neste post ilustro uma forma de lidar com dados georreferenciados (latitude/longitude) e plotar mapas. O pacote a ser usado é o ggmap, que possui uma função chamada geocode. Seu uso é bem simples: basta informar um endereço postal e recebemos de volta sua latitude e longitude. Detalhe: não pode haver acentos nas vogais, nem “ç”. Os caracteres aceitos são aqueles do teclado americano. :-/

Vamos obter a latitude e a longitude do MASP, Museu de Arte de São Paulo:

require(ggmap)
geocode("Av. Paulista, 1578 - Bela Vista, Sao Paulo, Brasil")

Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Av.+Paulista,+1578+-+Bela+Vista,+Sao+Paulo,+Brasil&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
lon lat
1 -46.65726 -23.56049

A função geocode acessa o API do Google Maps e faz uma pesquisa pelo endereço informado. Obviamente, nem todos os endereços são unívocos. Qualquer um que já utilizou o Google Maps sabe que esse aplicativo retorna um rol de alternativas em caso de múltiplas possibilidades. O resultado acima é a alternativa mais provável (e, na realidade, é a correta). Mas se for o caso, é possível acessar, dentro do objeto que é resultado dessa pesquisa, as demais alternativas. É possível também procurar vários endereços de uma vez só. Basta que cada endereço esteja gravado dentro de um elemento de um vetor de caracteres (string ou character). O uso do API, no entanto, tem uma limitação: só é possível fazer 2500 consultas por dia. Limites são muito comuns em APIs de acesso gratuito. Há planos pagos que permitem elevar aquele número.

O que fazer com a latitude e a longitude? É possível plotá-las num mapa!

No exemplo abaixo, vou até o site da prefeitura de SP, baixo o código fonte da página que contém uma lista de endereços de 13 unidades do SESC na região metropolitana (existem mais unidades… mas para o exemplo, essas bastam).

# Baixando, na página da prefeitura de SP, os endereços dos SESCs
pag_pref = readLines("http://www9.prefeitura.sp.gov.br/secretarias/smpp/sites/saopaulomaisjovem/cultura/index.php?p=25")[71]

Reparem que já extraio diretamente a linha de número 71. É exatamente nela que estão todos os endereços. Então baixar o resto era desnecessário. Mas esse conteúdo não está formatado:

pag

É preciso fazer uma limpeza do código HTML e extrair apenas o conteúdo desejado:

# Separando as linhas e removendo conteúdos desnecessários
pag_pref = unlist(strsplit(pag_pref,"<br />|/ Tel"))
pag_pref =gsub("<strong>|</strong>","",pag_pref)

# Mantém apenas as linhas que contêm a expressão "Endereço"
pag_pref = pag_pref[grep("Endereço",pag_pref)]

# Remove a expressão "Endereço"
pag_pref = gsub("Endereço: ","",pag_pref)

# Retira todos os caracteres especiais
pag_pref = gsub("[[:punct:]]", ",", pag_pref)

# Remove conteúdo desnecessário da linha 1
pag_pref = gsub("esquina com a Rua 24 de Maio","Sao Paulo, SP",pag_pref)

# Adiciona a cidade à linha 8
pag_pref[8] = paste(pag_pref[8],", Sao Paulo, SP")

# Adiciona o país a todas as linhas
pag_pref = paste(pag_pref,", Brasil")

# Remove todos os acentos
pag_pref = iconv(pag_pref, to="ASCII//TRANSLIT")

Pronto, esse é o resultado, uma lista de endereços guardada num vetor de caracteres:

pag

Agora é só aplicar a função geocode para obter todas as latitudes e longitudes:

latlon = geocode(pag_pref)

Falta o mapa.

A função get_map do pacote ggmap acessa um repositório público de mapas, o stamen, copia a imagem (geocodificada) da região desejada e a retorna como um objeto do R. Há várias opções de formatação, cores etc. E há também outros repositórios de mapas (inclusive o próprio Google Maps).

# Baixa o mapa de SP (centrado na Sé - isto pode ser alterado)
sp.map =get_map(location="Sao Paulo", zoom = 11,
 source = "stamen", maptype = "toner", color = "color")

# Transforma o arquivo de mapa em um grafico (ggplot2)
sp.map.2012 <- ggmap(sp.map, base_layer =
 ggplot(aes(x = lon, y = lat), data = latlon), extent = "device")

#Plota o resultado (os pontos dos endereços)
sp.map.2012 + geom_point(size = I(4),colour="red", alpha = 2/3)

Pronto, o mapa final é esse aqui:

plot_zoom

Mas o assunto sobre dados georreferenciados e mapas está longe de ser esgotado aqui…