Obtendo dados e plotando mapas no R

• Show All
• Hide All

O R é uma linguagem ou software que tem se tornado importante ferramenta para análise de dados espaciais e produção de mapas. Como exemplo, sugiro visitar o sítio web R-bloggers e fazer uma pequena busca utilizando as palavras “plot maps”. Vejam quantos resultados aparecerão. Para mais usos e maneiras de produzir mapas no R, veja também o sítio web Spatial data in R: Using R as a GIS.
Nesta postagem, quero demonstrar o quanto podemos produzir mapas informativos de uma maneira bem simples, com poucas linhas de comando. Como exemplo, vejam que com duas linhas apenas, podemos plotar o mapa mundi.

library(maps)
map()

O tutorial visa a produção de um mapa de distribuição geográfica utilizando dados que abrangem a distribuição espacial de duas espécies de Burseraceae, Protium aracouchini (Aubl.) Marchand e P. heptaphyllum (Aubl.) Marchand. A primeira planta faz parte de um complexo de espécies informalmente denominado completo Protium aracouchini que é o modelo de estudo de minha tese de doutorado. A segunda espécie é o modelo de estudo de um amigo, Gabriel Damasco, colaborador do LABOTAM e sob orientação do Dr. Paul Fine, também meu orientador.

Os dados foram obtidos através da base de dado do herbário virtual do Jardim Botânico de Nova Iorque.
Os tópicos a serem explorados neste tutorial são:

1. Importar os dados para o R;
2. Verificar e limpar os dados;
3. Preparação de variáveis a serem utilizadas no mapa;
4. Plotar o mapa.

Aviso importante!

Este tutorial necessita de conhecimentos básicos da linguagem por parte do usuário. Desta forma, caso você queira uma boa fonte para se iniciar na linguagem, recomendo visitar o sítio web da disciplina Preparação de dados para Análise Estatística do Programa de Pós-graduação em Ciências Biológicas (Botânica) do INPA. Lá você encontrará rica oferta de materiais para se familiarizar com o R e ser capaz de seguir este tutorial.

Passo 1 - Instalar o R

O passo essencial deste tutorial é ter a última versão do R instalado em seu computador. Para obtê-lo, visite a página do projeto R-project.
Caso haja dúvidas sobre o processo de instalação, favor visitar este link.

Passo 2 - Chamar os pacotes

Faremos uso dos pacotes maps (para plotar o mapa), RColorBrewer (paleta de cores), magrittr (pipeline) e dplyr (manipulação de ‘dataframes’) . É imprescindível que todos estejam atualizados e com as dependências instaladas.
Iniciando a sessão, execute os comandos abaixo para limpar a área de trabalho E chamar os pacotes para dentro da sessão.

rm(list=ls()) #limpa COMPLETAMENTE o ambiente de trabalho

#chama os pacotes
libs <- c('maps','RColorBrewer','magrittr','dplyr')
sapply(libs,library,character.only=T,logical.return=T)

Passo 3 - Importar os dados para a sessão

Os dados estão reunidos em dois arquivos ‘.csv’ separados por tabulação e codificados em UTF-8. Baixe-os para seu computador:

• Protium aracouchini
• Protium heptaphyllum

Cabe ressaltar que vários formatos podem ser utilizados para importar dados ao R. Utilizo o ‘CSV’ porque é mais simples.
Seguindo o tutorial, geralmente estabeleço inicialmente o diretório de trabalho usando a função setwd(). Desta forma, posso executar o script por inteiro, sem me preocupar em escolher o diretório utilizando botões, enfim, é mais prático. Ressalto que isso deve alterado por você para o seu diretório de trabalho!

#estabelece o diretorio de trabalho
setwd('/Users/ricoperdiz/Documents/DOC/PROJETO_DOC/LABOTAM/data')

Chamamos o arquivo para dentro da sessão e utilzamos o operador %>% para passar o resultado da importação para a próxima ação, que é a de ordenar esses dados segundo o epíteto específico, nome do coletor e número de coleta. Ao resultado dessa ação é atribuído um objeto, ‘cpa’. Veja como estão as primeiras linhas desse objeto após a ordenação.

#importa os dados - importante que haja uma coluna com dados de lat e outra de long
#protium aracouchini
read.table("nybg_paracouchini.csv", header = T, as.is = T, sep = '\t', dec = '.') %>%
    select(recordedBy, recordNumber, decimalLatitude, decimalLongitude, identifiedBy, specificEpithet) %>%
    arrange(recordedBy, recordNumber) -> aracouch
#protium heptaphyllum
read.table("nybg_pheptaphyllum.csv", header = T, as.is = T, sep = '\t', dec = '.') %>%
    select(recordedBy, recordNumber, decimalLatitude, decimalLongitude, identifiedBy, specificEpithet) %>%
    arrange(recordedBy, recordNumber) -> heptaphy

Passo 4 - Verifica os dados

Em nosso caso, devemos verificar: * existência valores vazios; em caso positivo, devemos eliminá-los; * confiabilidade dos valores de latitude e longitude, às vezes, por diversos fatores, há troca de sinais (negativos e positivos) ocasionando equívocos quanto à ocorrência exata da amostra. Se for percebido algo assim, é bom checar os dados e buscar corrigí-los.

############################################################
############################################################
### Limpeza e verificação de dados
#une os dados e passa para a proxima acao
rbind(aracouch,heptaphy) %>%
    #elimina os registros vazios de coletor
    filter(recordedBy != '') %>%
    #elimina os registros sem lat ou long
    filter(decimalLatitude != '' | decimalLongitude != '') %>%
    #filtra apenas os especimes identificados pelo especialista
    filter(identifiedBy == 'D. C. Daly') -> dados
#ha registros duplicados
#busca-se entao apenas os dados unicos
#limpa os dados de coletor e numero para poder criar um identificador
#limpeza consiste em eliminar '.', espacos vazios, apostrofe e '_' duplos
gsub('\\.', '_', dados$recordedBy) %>%
    gsub(' ', '_', .) %>%
    gsub("'", '_', .) %>%
    gsub('__', '_', .) -> dados$recordedBy
dados$recordedBy
#faz-se o mesmo para os numeros de coleta
gsub('/','_', dados$recordNumber) %>%
    gsub(' ', '_', .) %>%
    gsub('\\.', '_', .) -> dados$recordNumber
dados$recordNumber
#cria o identificador de coleta e especie
dados %>%
    mutate( ID = paste(recordedBy,recordNumber, sep = '_'), Species = paste('Protium', specificEpithet, sep = ' ')) -> dados
dados$ID <- gsub('__','_',dados$ID)
#quem sao os dados unicos
unicos <- unique(dados$ID)
unicos
#agora filtra os dados unicos no dataframe, eliminando os duplicados
#faz-se uso da funcao match para obter esse resultado
match(unicos, dados$ID) %>%
    dados[.,] -> prot
#verifica a cobertura de lat e long para ver se estao dentro
# dos limites da America do Sul
range(prot$decimalLatitude)
range(prot$decimalLongitude) #aqui tem algo estranho
#percebe-se aqui que ha valores que caem fora da Am Sul
#limite e pouco mais de -80
head(sort(prot$decimalLongitude))
#devemos eliminar
prot %>%
    filter(decimalLongitude > -80) -> prot

Passo 5 - Prepara variáveis para plotar o mapa

#vetor para lat e long
lat <- prot$decimalLatitude
long <- prot$decimalLongitude

#vetor com epitetos especificos
spp <- unique(prot$Species)

#vetor de cores para utilizar no mapa
#para cada especie, uma cor
prot$cores.map <- ifelse(prot$Species == spp[1], 'red','black')
cores.map <- prot$cores.map
#cria um vetor de tamanho para cada especie
#como uma coluna de prot
#P heptaphyllum possui uma distribuicao mais ampla
#por isso atribuo um tamanho menor pra ela
prot$cex.p <- ifelse(prot$Species == spp[1], 1, 0.8)

#vetor para lat e long
y1 <- range(lat) + c(-1,1)
x1 <- range(long) + c(-1,1)

#vetor para os simbolos pch
prot$pontos <- ifelse(prot$Species == spp[1], 21, 24)

#vetor para paises
paises <- c('Brazil','Argentina','Peru','Paraguay','Ecuador','Chile','Uruguay','French Guiana','Suriname','Venezuela','Colombia','Guyana','Bolivia','Panama','Costa Rica')

Passo 6 - Plota o mapa

#plota o mapa
map(regions = paises, fill = F, xlim = x1, ylim = y1) 

#plota os pontos de ocorrencia de cada especie
points(long,lat, pch = prot$pontos, col = cores.map, bg = cores.map, cex = prot$cex.p)

#coloca eixos das coordenadas
map.axes()
axis(side=4,las=1)
axis(side=3,las=1)

#coloca escala do mapa
par(cex=1, las=1)
map.scale(max(long) - 12, ratio = F, cex = 1, metric = T)

#plota a linha do equador
abline(h=0,lwd=0.5,lty="dotted")
#nomeia a linha do equador
text(x = max(long) - 1, y = 1,labels="Equador", cex=1, adj=c(1,0.5))
#plota uma legenda
legend(max(long) - 16 , max(lat),legend = spp, pch = unique(prot$pontos), pt.bg = unique(cores.map), cex = 0.8, x.intersp = 0.8, text.font = 3)

Pode-se salvar também este mapa em um pdf, utilizando a função pdf(). Veja abaixo.

### Cria um pdf e Plota o mapa 

#inicia o pdf
pdf('meu_mapa_complexo_protium_aracouchini.pdf',height=8,width=10)
map(regions = paises, fill = F, xlim = x1, ylim = y1) 
points(long,lat, pch = prot$pontos, col = cores.map, bg = cores.map, cex = prot$cex.p)
map.axes()
axis(side=4,las=1)
axis(side=3,las=1)
par(cex=1, las=1)
map.scale(max(long) - 12, ratio = F, cex = 1, metric = T)
abline(h=0,lwd=0.5,lty="dotted")
text(x = max(long) - 1, y = 1,labels="Equador", cex=1, adj=c(1,0.5))
legend(max(long) - 16 , max(lat),legend = spp, pch = unique(prot$pontos), pt.bg = unique(cores.map), cex = 0.8, x.intersp = 0.8, text.font = 3)
dev.off()

---

Conclusão

Simples demais… Brincadeira! Com o tempo fica fácil.
Todavia, há muito outras possibilidades, como, por exemplo, pode-se importar shapefiles, plotar rios, colocar cores etc. Em próximo post, pretendo demonstrar como importar shapefiles para dentro da sessão e utilizá-los para incrementar o mapa.
Dúvidas, sugestões, favor entrar em contato via comentários.
O script completo deste tutorial pode ser encontrado aqui.
Esta postagem foi escrita em ambiente R utilizando os pacotes R Markdown e knitr. Você pode reproduzir a postagem baixando o arquivo ‘Rmd’ em seu computador (baixe-o aqui, em arquivo zipado; faça a descompressão antes de fazer uso do mesmo) e processando-a com o pacote ‘knitr’. Por questões de portabilidade entre R e a plataforma Blogger, este arquivo ‘Rmd’ não contem as imagens desta postagem.
Divirta-se!!!

---

Para saber mais

• Making maps with R
• Maps in R: Introduction