23 junho 2017

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

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:
É 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.

Referências dos pacotes de R

Publicado por ::
Tags :: , , , , , , , , , , , , , ,

Um comentário:

  1. gdvaleoutubro 12, 2017 1:36 PM

    Demais galera!! Muito legal esse tutorial!! Saci cade vc meu filho?

    ResponderExcluir

Assinar: Postar comentários (Atom)