Tutorial em R para plotar mapa de distribuição de uma espécie de Leguminosae está disponível no GitHub

Disponível tutorial para gerar mapa publicado em artigo publicado por membros do LABOTAM/INPA no periódico Phytotaxa em 2018

Distribuição geográfica de Macrolobium aracaense.

Nesta postagem, apresento um tutorial em linguagem R para gerar o mapa publicado no artigo "A new species of Macrolobium (Fabaceae, Detarioideae) endemic on a Tepui of the Guyana Shield in Brazil" (Farroñay et al. 2018). Este trabalho foi publicado em 2018 por membros do Laboratório de Ecologia e Evolução de Plantas Amazônicas, grupo de pesquisa filiado ao Instituto Nacional de Pesquisas da Amazônia (INPA), Manaus, Amazonas, Brasil. Tanto o tutorial quanto todos os arquivos necessários para gerar o mapa, como dados de espécimes botânicos e shapefiles, estão hospedados em um repositório de minha conta no GitHub. Clique AQUI para acessar o tutorial.
O tutorial é apresentado no formato de um Jupyter Notebook, que é uma ferramenta super legal para aprendizado, especialmente na Ciência. Vale a pena checar o link anterior e os disponíveis ao fim desta postagem para que vocês aprendam um pouco mais sobre esta ferramenta. Caso haja alguma dificuldade em fazer uso dela em seu computador, também forneço na mesma pasta um script simples em R para que vocês possam executar os mesmos passos apresentados no Jupyter Notebook. De maneira geral, para poder fazer uso de um Jupyter Notebook em seu computador, é necessário ter uma distribuição Python instalada no computador (veja aqui o passo a passo para instalar e rodar o Jupyter Notebook localmente).
Para gerar o mapa, fiz uso de alguns pacotes de R que gostaria de destacar abaixo:

• ggmap, que permite utilizar mapas em formato raster oriundos de serviços de mapeamento online populares, como o Google Maps e Stamen Maps;
• alguns da coleção de pacotes tidyverse (ggplot2, dplyr, stringr, readr, purrr, magrittr, e broom);
• cowplot, utilizado para mesclar diferentes plots oriundos do ggmap/ggplot2;
• ggsn, utilizado para plotar a seta de Norte e a escala.

Espero que aproveitem e façam bom uso, enquanto ferramenta de ensino e aprendizado. Se tiverem dúvidas, encontrarem algum erro, quiserem elogiar, favor deixar uma mensagem aqui, ou então no GitHub.

Em breve, postarei mais novidades.  Até mais!

Saiba mais

• Jupyter Notebook
• Como instalar Python e Jupyter Notebook usando Anaconda
• Uma modesta introdução ao Jupyter Notebook
•  A new species of Macrolobium (Fabaceae, Detarioideae) endemic on a Tepui of the Guyana Shield in Brazil" (Farroñay et al. 2018)

Obtendo dados e plotando mapas no R - Versão 2

Continuando a reformulação de postagens antigas, reapresento uma publicada em 2016 sobre obtenção de dados e plotagem de mapas usando o R. Fiz correções no script lá apresentado, e eliminei o desnecessário (como atribuição de pasta de trabalho e caminhos relativos para ler os dados), de forma a tornar o tutorial abaixo totalmente reprodutível. Por isso, peço encarecidamente que, uma vez encontrado algum erro, favor me avisar, para que eu corrija o mais rapidamente possível e você possa fazer bom uso desta postagem.

---

PARA BAIXAR A POSTAGEM COMPLETA COM O TUTORIAL, SIGA A INSTRUÇÃO ABAIXO:

PARA REPRODUZIR A POSTAGEM EM SEU COMPUTADOR, BAIXE o arquivo .Rmd juntamente com todos os dados necessários para gerar o mapa reproduzido no fim desta postagem NESTE LINK, procure o arquivo de extensão “.Rmd” e abra-o em seu computador, preferencialmente usando o RStudio. Execute o arquivo usando a função Preview presente no RStudio, certificando-se de que você tenha os pacotes knitr, rmarkdown e todas as dependências devidamente instaladas em seu computador. Esta postagem foi construída usando o formato R Notebook. Visite os links ao fim desta postagem para obter mais informações sobre este formato e seus benefícios para a reprodutibilidade na Ciência.

---

O tutorial é bem simples e 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 complexo Protium aracouchini, cuja sistemática e taxonomia vem sendo objeto de estudo em meu doutorado. A segunda espécie, P. heptaphyllum, é objeto de estudo de Gabriel Damasco, colaborador do LABOTAM e aluno de doutorado na Universidade da Califórnia, Berkeley.
Os dados foram baixados do herbário virtual do Jardim Botânico de Nova Iorque que é onde trabalha o especialista em Burseraceae, Dr. Douglas Daly, e cuja coleção de espécimes e dados vêm sendo bem cuidada há algumas décadas. Para esta postagem, filtramos apenas os dados para as duas espécies citadas acima.

---

Esta postagem foi escrita em ambiente R utilizando os pacotes R Markdown e knitr.
Dúvidas, sugestões, erros no script, favor entrar em contato comigo, via comentários, ou por email.
Esta postagem completa, incluindo os blocos de código, pode ser lida inteiramente neste link. Lá você pode copiar o arquivo bruto para o seu computador e modificá-lo à vontade.
Divirta-se!!!

Composição de mapas em R

Esta postagem é uma reformulação de uma anterior, publicada em 2017, neste mesmo blog (veja-a aqui). Aqui ela é apresentada na forma de um R notebook. Foram feitas algumas modificações a fim de torná-la totalmente reprodutível, isto é, você pode executar este arquivo do início ao fim que ele vai rodar, desde que você tenha os pacotes necessários instalados em seu computador.
Baixe o arquivo .Rmd zipado e com todos os dados e shapes necessários para gerar o mapa reproduzido abaixo neste link. Procure o arquivo de extensão ".Rmd" e abra-o em seu computador, preferencialmente usando o RStudio.
Em caso de dúvidas, sugestões, erros na execução do código, entre em contato deixando seu comentário. Até a próxima!
Divirta-se!

---

Composição de mapas em R com a função layout

Nesta postagem, quero apresentar um script em linguagem R que utilizei para gerar o mapa publicado recentemente em Rodrigues et al. 2017 (caso você possua conta no ResearchGate, pode visualizar o trabalho neste link). Nele, faço uso da função layout (pacote graphics) para compor a figura abaixo, unindo dois diferentes mapas e uma legenda através de uma única sequência de código.

O mapa principal apresenta localidades de coleta de espécimes de plantas em duas unidades de conservação do estado de Roraima. A maioria dessas coletas foram feitas durante a expedição Terra Incognita, ocorrida em dezembro de 2013, organizada pelas equipes gestoras do Parque Nacional Serra da Mocidade e Estação Ecológica de Niquiá, com financiamento do ARPA. Veja ao fim desta postagem mais detalhes sobre a expedição.

---

Partes do tutorial

• Passo 1 - Estabelecer o diretório de trabalho
• Passo 2 - Chamar os pacotes
• Passo 3 - Importar os dados e criar variáveis para os plots
• Passo 4 - Gera uma composição gráfica utilizando a função layout
• Passo 5 - Plota os mapas individualmente para serem alocados em cada painel

---

Avisos importantes!

Conhecimentos básicos de R

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.

BAIXE OS DADOS DESTA POSTAGEM!!!

Baixe os dados AQUI!!!!
Eles estão reunidos em quatro pastas dentro de um arquivo .zip. Elas contêm:

• o código completo (pasta code), disponibilizado ao fim desta postagem;
• dados de localidades de coleta (pasta data) incluídos em um arquivo .csv;
• a figura gerada nesta postagem (pasta figures), exatamente igual ao publicado em Rodrigues et al. 2017;
• os arquivos de SIG (pasta gis);
• Fora das pastas, disponibilizo também o arquivo .Rmd, que contem o texto desta postagem, e o arquivo .html, que é o output do arquivo .Rmd.

Baixe este arquivo .zip para seu computador, descomprima-o e coloque tudo dentro de uma pasta. Use esta pasta para trabalhar pois meu código faz uso do caminho que leva a cada um desses arquivos (path), por isso, NÃO RENOMEIE AS PASTAS!

Sintaxe utilizada para algumas funções

Para algumas das funções utilizadas neste script, eu adicionei o nome do pacote antes de cada uma dessas funções. Assim, fica claro a qual pacote pertence aquela determinada função. Por exemplo, ao usar a função

readOGR

eu adicionei o nome do pacote antes da função, ficando a sintaxe assim:

rgdal::readOGR

---

Passo a passo

Passo 1 - Estabelecer o diretório de trabalho

Geralmente estabeleço inicialmente o diretório de trabalho usando a função

setwd().

e acrescentando o local onde está minha pasta de trabalho. Alguns não recomendam esta ação por não fornecer portabilidade ao script. Porém, eu gosto de usá-la e aqui eu coloco. Gosto pois posso executar o script por inteiro de uma vez só, ao selecionar todo o conteúdo e executá-lo.
Porém, você deve se atentar a SEMPRE colocar o caminho do seu diretório de trabalho dentro dessa função. Ressalto que isso deve alterado por você para o seu diretório de trabalho!
Caso você se enrole demais com essa parte do script, sugiro que apague esta linha e, sempre ao iniciar a sessão, use os botões para selecionar o diretório de trabalho.
Pois bem, também possuo o costume de, ao iniciar a sessão, limpar completamente a área de trabalho.
Por isso, FIQUE ATENTO caso você tenha algum objeto de importância para algum trabalho em sua área de trabalho. Pois após executar o comando abaixo,

rm(list=ls())

qualquer objeto deixado na área de trabalho sumirá completamente…

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

# esse e o meu diretorio de trabalho
# ATUALIZEM ESSE LOCAL, colocando o SEU DIRETORIO DE TRABALHO
setwd('~/Documents/DOC/PROJETO_DOC/LABOTAM/posts_blog/post_layout/')

Passo 2 - Chamar os pacotes

Para gerar este mapa, fiz uso dos pacotes abaixo:

• broom;
• dplyr;
• GISTools;
• maps;
• rgdal.

É imprescindível que todos estejam instalados, atualizados e com as dependências instaladas.

#chama os pacotes

# primeiro estocamos os nomes dos pacotes em um objeto ...
libs.maps <- c('broom','dplyr', 'GISTools', 'maps', 'rgdal')
# ... para em seguidar chamá-los para a sessão utilizando o sapply
sapply(libs.maps, library, character.only=T,logical.return=T)

Passo 3 - Importar os dados e criar variáveis para os plots

# ----------------------------------------------------------------------------------------
#le os dados
# ----------------------------------------------------------------------------------------
#dados dos pontos de coleta dos vouchers
new.reg.voucher <- read.table('data/rodrigues_etal_2017_coordenadas_novos_registros.csv',header=T,as.is=T,sep='\t')

#kml do PARNA S Mocidade e ESEC Niquia
mocidade <- rgdal::readOGR(dsn='gis/parna_serra_da_mocidade.kml',layer='sql_statement')
niquia <- rgdal::readOGR(dsn='gis/esec_niquia.kml',layer='sql_statement')

#shape da Am Sul e Central
area.mapa <- rgdal::readOGR(dsn = 'gis/SAm_CAm_shape.shp')

#shape de RR
rr <- rgdal::readOGR(dsn = 'gis/roraima.shp')

#shape de rios em RR
rios.main <- rgdal::readOGR('gis','rios',encoding='latin1')
rios.sec <- rgdal::readOGR('gis','rios2',encoding='latin1')

# ----------------------------------------------------------------------------------------
# cria variaveis
# ----------------------------------------------------------------------------------------

#mapa de referencia
#amplitude de long
x1 <- c(-70,-50)
y1 <- c(-10,5)

#vetor de tamanho para o cex
cex.tam <- seq(0.8,1.6,by=0.1)

#vetor com tipos de pontos
pontos <- c(17:25)


#amplitude de lat e long - vetores
lat_long <- dplyr::full_join(tidy(mocidade), tidy(niquia))
y.geral <- range(lat_long$lat) + c(-0.01,0.01)
x.geral <- range(lat_long$long) + c(-0.01,0.01)

Passo 4 - Gera uma composição gráfica utilizando a função layout

Chegamos ao que interessa. A função layout permite compor complexos gráficos em R, pois habilita o usuário a criar uma matriz em que será especificado onde cada plot será alocado na imagem. Dê uma olhada no help file da função layout com o comando abaixo:

?layout

Veja o código abaixo:

#layout do mapa
mapa.layout <- layout(
    mat = matrix(c(1,1,2,2,2,2,2,2,3,3,2,2,2,2,2,2,3,3,2,2,2,2,2,2),nrow=3,byrow=T),
    widths=rep(2/8,8),
    heights=rep(1/3,3),
    respect=T)

O primeiro argumento de layout é mat, que é a matriz. É nesta matriz que se especifica onde cada plot será alocado. Plot número 1 será alocado no painel 1; plot número 2, painel 2; plot número 3, painel 3; e assim por diante.
Para poder visualizar como ficou seu esquema de composição de plots, a partir do objeto mapa.layout, use o comando abaixo:

layout.show(mapa.layout)

Você obterá como resultado o gráfico abaixo:

Passo 5 - Plota os mapas individualmente para serem alocados em cada painel

Agora vamos plotar cada painel. Neste caso, eu escolhi o painel 1 para situar o leitor quanto à localização de Roraima e das unidades de conservação visitadas durante a expedição Terra Incognita. O painel 2 foi o escolhido para evidenciar com mais detalhes as localidades de coleta das plantas, enquanto o painel 3 foi utilizado para alocar a legenda. Veja abaixo o restante do script.

# ----------------------------------------------------------------------------------------
#plot 1
# ----------------------------------------------------------------------------------------


par(mar=c(0,2,2,2))
plot(area.mapa,xlim=x1,ylim=y1,lwd=0.2)
plot(rr,add=T,col='gray90')
#plota um retangulo na area do mapa da direita
rect(xleft = min(x.geral), ybottom = min(y.geral), xright = 
         max(x.geral), ytop = max(y.geral), density = 25)
#coloca escala do mapa
par(cex=1, las=1)
#library(GISTools)
# inclui o NORTE
GISTools::north.arrow(max(x1)-5,max(y1)-10,len=1,cex.lab=0.9,lab='N', col = 'black')
#coloca um texto para avisar quem é Brasil e Guiana
text(-58,-10,labels='BRASIL',pos=4,cex=cex.tam[1],font=2)
box()

# ----------------------------------------------------------------------------------------
#plot 2
# ----------------------------------------------------------------------------------------


par(mar=c(2,2,2,3))
#plota o mapa das uc`s com dados específicos
plot(rr,col='gray90',xlim=x.geral,ylim=y.geral,lwd=1,lty=1) #default lwd=1,lty=1

plot(mocidade,add=T,lty=0,col='gray80')
plot(niquia,add=T,lty=0,col='gray80')

#plota os rios
plot(rios.main,add=T,col='gray60')
plot(rios.sec,add=T,col='gray60')

#plota as areas dos mapas
plot(mocidade,add=T,lwd=1.4,lty=2)
plot(niquia,add=T,lwd=1.4,lty=3)

#plota os pontos de coleta
points(new.reg.voucher$long,new.reg.voucher$lat,pch=pontos[3],col='black',cex=1.4)

#coloca nome dos rios
text(x=-61.42,y=0.8,labels='Rio Branco',pos=1,srt=60,font=3)
text(x=-61.65,y=1.3,labels='Rio Água Boa do Univini',pos=1,srt=70,font=3)
text(x=-61.97,y=0.77,labels='Rio Catrimani',pos=3,srt=-50,font=3)
text(x=-61.9,y=1.2,labels='Rio Capivara',pos=3,srt=-60,font=3)
#coloca nome das UC's
text(x=-61.80,y=1.5,labels='PARNA\nS. Mocidade',pos=1,font=2)
text(x=-61.50,y=1.3,labels='ESEC\nNiquiá',pos=1,font=2)

# coloca um texto para informar área do Amazonas e Roraima
text(-62.845,sum(range(y.geral))/2,labels='Amazonas',pos=4,cex=1.4,font=2)
text((sum(range(x.geral))/2)-0.1,max(y.geral)-0.02,labels='Roraima',pos=4,cex=1.4,font=2)

#coloca escala do mapa
par(cex=1, las=1)
maps::map.scale(max(x.geral) - 0.15, min(y.geral)+0.02, relwidth = 0.10, ratio = F, cex = 1, metric = T, col = 'black')
#coloca o Norte
GISTools::north.arrow(max(x.geral),min(y.geral)+0.07,len=0.02,lab='N',cex.lab=1.2,lwd=1.5,col='black')
#coloca eixos das coordenadas
maps::map.axes()
axis(side=4,las=1)
axis(side=3,las=1)

# ----------------------------------------------------------------------------------------
#plot 3 - LEGENDA
# ----------------------------------------------------------------------------------------


par(mar = c(1,1,0,1))
plot.new()
par(cex = 1.2)
legend(x = 'center', ncol = 1, legend = 
           c('Pontos de coleta dos\nvouchers','Parque Nacional Serra da\nMocidade','Estação Ecológica de\nNiquiá','Fronteira entre estados','Cursos de água'),
       pch = c(pontos[3], NA, NA, NA, NA), col = c('black','black','black','black','gray60'), title = 'LEGENDA',
       cex = cex.tam[3], lty = c(0,2,3,1,1), pt.lwd = 2, text.width = 0.8, lwd = 3, y.intersp = 1.8)

Por fim…

Para terminar, lembre-se que você pode salvar este mapa em um pdf, utilizando a função pdf(). Dê uma olhada no post Obtendo dados e plotando mapas no R para aprender (ou se lembrar) como fazer. O script completo que gera este mapa já contem a função pdf() dentro do código. Você pode baixá-lo seguindo os links que estão no início e fim desta postagem.

---

Informações da sessão

---

Conclusão

Por favor, me avise caso algum erro seja encontrado. Dúvidas, sugestões, favor entrar em contato via e-mail!

Informações sobre esta postagem

Esta postagem foi escrita em ambiente R utilizando os pacotes R Markdown e knitr.
Novamente, caso tenha interesse nos produtos desta postagem, que inclui o script completo deste tutorial, o arquivo ‘Rmd’ que gerou esta postagem em html, a figura aqui produzida, dados de localidade, e arquivos de SIG podem ser baixados para seu computador seguindo este link. Os dados estão zipados, portanto, faça a descompressão antes de fazer uso do mesmo.
Você pode reproduzir esta postagem utilizado o arquivo ‘Rmd’ em seu computador 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!!!

---

Sobre as unidades de conservação aqui citadas

O Parque Nacional (PARNA) Serra da Mocidade e a Estação Ecológica (ESEC) de Niquiá são unidades de conservação federal situadas no município de Caracaraí, no estado de Roraima, Brasil. O PARNA Serra da Mocidade possui uma área de 350.960 hectares e caracteriza-se pela ocorrência de uma serra isolada, com aproximadamente 1.800 metros de altitude, localizada no limite da unidade conservação com a área indígena Yanomami. Por outro lado, a ESEC Niquiá possui 284.787,42 hectares, apresentando na sua maior parte áreas planas, com duas elevações ente 200-500 m de altitude.

Sobre a expedição

No período de 04 a 18 de dezembro o Parque Nacional Serra da Mocidade e a Estação Ecológica Niquiá realizaram a Expedição de Pesquisa Científica Terra incognita, nomeada em alusão à Hamilton Rice, explorador americano que assim chamou a região oeste de Roraima em seu livro Exploração à Guiana Brasileira. Financiada com recursos do Programa ARPA (Áreas Protegidas da Amazônia), a expedição reuniu várias instituições (Instituto Chico Mendes de Conservação da Biodiversidade/ICMBio, Universidade Federal de Roraima/UFRR, Instituto de Amparo à Ciência e Tecnologia do Estado de Roraima/IACTI, Universidade Federal do Amazonas/UFAM, Centro de Estudos da Biodiversidade Amazônica/CENBAM/PPBio, Secretaria do Patrimônio da União/SPU e Instituto Federal de Roraima/IFRR) e pesquisadores de diversas áreas do conhecimento, que se propuseram a levantar informações sobre uma região pouco conhecida pela ciência, buscando contribuir com a elaboração dos planos de manejo das unidades, ora em execução.
Vejam este link para ter acesso a um pequeno documentário sobre algumas das descobertas dessa viagem.

---

Saiba mais

• R plot layout;
• Multipanel plotting in R (with base graphics).

Para conferir o trabalho onde esse mapa foi publicado siga o link abaixo:

• Rodrigues et al. 2017. Novas ocorrências de angiospermas para o estado de Roraima, Brasil. Rodriguésia - Link no ResearchGate.

Para ver outros mapas produzidos da mesma maneira:

• Lavor et al. 2016. Rediscovery of Pilosocereus oligolepsis (Cactaceae) in the State of Roraima, Brazil. Cactus and Succulent Journal 88(3):137-143 - Link no ResearchGate;
• Barbosa et al. 2017. Decomposition rates of coarse woody debris in undisturbed Amazonian seasonally flooded and unflooded forests in the Rio Negro-Rio Branco Basin in Roraima, Brazil. Forest Ecology and Management 397(1) - Link no ResearchGate.

Referências dos pacotes de R

• Original S code by Richard A. Becker, Allan R. Wilks. R version by Ray Brownrigg. Enhancements by Thomas P Minka and Alex Deckmyn. (2017). maps: Draw Geographical Maps. R package version 3.2.0. https://CRAN.R-project.org/package=maps
• Roger Bivand, Tim Keitt and Barry Rowlingson (2017). rgdal: Bindings for the Geospatial Data Abstraction Library. R package version 1.2-7. https://CRAN.R-project.org/package=rgdal
• Chris Brunsdon and Hongyan Chen (2014). GISTools: Some further GIS capabilities for R. R package version 0.7-4. https://CRAN.R-project.org/package=GISTools
• David Robinson (2017). broom: Convert Statistical Analysis Objects into Tidy Data Frames. R package version 0.4.2. https://CRAN.R-project.org/package=broom
• Hadley Wickham, Romain Francois, Lionel Henry and Kirill Müller (2017). dplyr: A Grammar of Data Manipulation. R package version 0.7.0. https://CRAN.R-project.org/package=dplyr

Importando shapefiles no R

• Show All
• Hide All

Shapefile é um formato de dados vetoriais geoespaciais muito utilizado em ambientes SIG. O ambiente R permite que se trabalhe com esse tipo de dados. Nesta postagem, quero demonstrar como podemos importar esse tipo de arquivo ao R em duas diferentes maneiras, a fim de plotar mapas utilizando esses dados.

Os tópicos a serem explorados neste tutorial são:

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

---

O que é necessário ter instalado em seu computador

• R instalado e atualizado;
• Pacotes ‘maptools’ e ‘rgdal’ atualizados.

---

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.

---

O passo a passo abaixo disponibilizado tenta apresentar os dados de maneira simples. Mas vejam que podemos cumprir a idéia central desta postagem com apenas três linhas.

library(maptools)
paises_shape <- readShapePoly('countries.shp')
plot(paises_shape)

---

Passo 01 - Estabelecer o diretório de trabalho

Tenho por costume o hábito de estabelecer o diretório de trabalho no início de todo script. Lembre-se de SEMPRE COLOCAR OS ARQUIVOS NECESSÁRIOS para rodar o tutorial DENTRO desta PASTA DE TRABALHO. ADAPTE o TEXTO ABAIXO para a SUA PASTA DE TRABALHO!!!

rm(list = ls()) #limpa a area de trabalho
#estabelece o diretorio de trabalho
setwd('/Users/ricoperdiz/Documents/DOC/PROJETO_DOC/LABOTAM/posts_blog/post_shapefile')

Passo 02 - Baixar shapefiles da web

Existem vários sítios web que disponibilizam shapefiles gratuitamente. Aqui vamos utilizar os dados disponíveis em DIVA-GIS. Clique no link abaixo para obter um arquivo ‘.zip’ contendo os shapefiles a serem utilizados neste tutorial.

• Dados

Após baixá-lo, deve-se descomprimir os arquivos e colocá-los em sua pasta de trabalho.

ATENÇÃO!!! ACESSE a pasta countries_shp e retire os quatros arquivos, de extensão ‘.shp’, ‘.dbf’, ‘.shx’ e ‘.prj’. PEGUE-OS e COLOQUE-OS na pasta de trabalho. APENAS os arquivos, NÃO INCLUA a pasta countries_shp.

Passo 03 - Importar os dados para a sessão

Os dados estão dispostos em um arquivo ‘.csv’ e estão separados por tabulação e codificados em UTF-8. Baixe-o para seu computador:

• Dados de espécimes de herbário de Protium heptaphyllum e Protium aracouchini

Cabe ressaltar que vários formatos podem ser utilizados para importar dados ao R. Utilizo o ‘CSV’ porque é mais simples.

---

Passo 04 A

USANDO O PACOTE rgdal

#chama o pacote
library('rgdal')

#importa os dados para dentro da sessao
paises_shape <- readOGR(dsn = getwd(), layer = 'countries')

#cria um vetor com nomes de paises que queremos plotar no mapa
paises <- c('Brazil','Argentina','Peru','Paraguay','Ecuador','Chile','Uruguay','French Guiana','Suriname','Venezuela','Colombia','Guyana','Bolivia','Panama','Costa Rica')

#filtra os dados para uma lista de paises
dados <- paises_shape[paises_shape$NAME%in%paises,]

#apaga os dados desnecessarios do shape original
rm(paises_shape)

# DADOS DE OCORRENCIA
# chama os dados de ocorrencia de protium heptaphyllum
# obtido na base de dados do herbario virtual do NYBG
protium <- read.table("dados_phepta_aracouc.csv", header = T, as.is = T, sep = '\t', dec = '.')

############################################################
############################################################
### Vetores

#vetores de lat e long
lat <- protium$decimalLatitude
long <- protium$decimalLongitude

#vetor de amplitude
ylat <- range(lat) + c(-1,1)
xlong <- range(long) + c(-1,1)

#vetor para nomes das especies
spp <- unique(protium$Species)

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

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

############################################################
############################################################
### Plota os dados

plot(dados, xlim = xlong, ylim = ylat)
points(long, lat, pch = protium$pontos , col = protium$cores.map, bg = protium$cores.map, cex = protium$cex.p)
legend(max(long) - 16 , max(lat),legend = spp, pch = unique(protium$pontos), pt.bg = unique(cores.map), cex = 1, x.intersp = 0.8, text.font = 3)

Passo 04 B

USANDO O PACOTE maptools

#chama o pacote
library('maptools')

#le o shapefile e cria um objeto com esses dados
paises_shape <- readShapePoly('countries.shp')

#cria um vetor com nomes de paises que queremos plotar no mapa
paises <- c('Brazil','Argentina','Peru','Paraguay','Ecuador','Chile','Uruguay','French Guiana','Suriname','Venezuela','Colombia','Guyana','Bolivia','Panama','Costa Rica','Belize','El Salvador','Guatemala','Honduras','Nicaragua','Mexico')

#filtra os dados para uma lista de paises
dados <- paises_shape[paises_shape$NAME%in%paises,]

#apaga os dados desnecessarios do shape original
rm(paises_shape)

# DADOS DE OCORRENCIA
# chama os dados de ocorrencia de protium heptaphyllum
# obtido na base de dados do herbario virtual do NYBG
protium <- read.table("dados_phepta_aracouc.csv", header = T, as.is = T, sep = '\t', dec = '.')

############################################################
############################################################
### Vetores

#vetores de lat e long
lat <- protium$decimalLatitude
long <- protium$decimalLongitude

#vetor de amplitude
ylat <- range(lat) + c(-1,1)
xlong <- range(long) + c(-1,1)

#vetor para nomes das especies
spp <- unique(protium$Species)

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

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

############################################################
############################################################
### Plota os dados

plot(dados, xlim = xlong, ylim = ylat)
points(long, lat, pch = protium$pontos , col = protium$cores.map, bg = protium$cores.map, cex = protium$cex.p)
legend(max(long) - 16 , max(lat),legend = spp, pch = unique(protium$pontos), pt.bg = unique(cores.map), cex = 1, x.intersp = 0.8, text.font = 3)

---

Conclusão

Os comandos aqui apresentados são muito simples, há muito a se explorar em novas postagens, tanto na utilidade das funções aqui apresentadas dos dois pacotes citados (rgdal e maptools) quanto na variedade de shapefiles a serem utilizados. Deixemos para próximas postagens.

Dúvidas, sugestões, elogios, críticas, correções, favor entrar em contato via comentários. Sua opinião é sempre bem vinda.

Caso você encontre algum erro ou tenha tido problemas em rodar os comandos, ficarei muito feliz de poder consertar ou ajudar!

Tenha um bom divertimento e aprendizado!

---

Para saber mais:

• Things I Forget: Reading a Shapefile in R with readOGR
• Read and write ESRI Shapefiles with R
• Shapefiles and projections in R