Cantinho do R

Mapinhas no R

Yes, ele faz isso sim!

Por Marcos Vital, do LEQ-UFAL

Este material foi construído com a ajuda de muitas pessoas que acreditam no LEQ e em Ciência Livre. Muito obrigado!

Para mais material, visite o Cantinho do R


Eu falo (e escrevo!) com bastante frequência que uma das coisas bacanas do R é a flexibilidade, esta capacidade de fazer um monte de análises diferentes. E este script é um pouco sobre isso! Você com certeza já espera que o R seja capaz de fazer tudo quanto é análise estatística, né? Mas esperava que ele fosse capaz de criar mapas como um SIG? Pois é, ele faz isso sim! (se você já sabia, seja educado, e faça uma cara de supresa, tipo quando você tem uma festa surpresa de aniversário mas já sacou e finge supresa mesmo assim pra galera não ficar chateada…)

(As coisas que a gente aprende totalmente por acaso: fui “linkar” a página da wikipedia que explica o que é um SIG ali em cima, e abri pra dar uma olhada (sabe como é, né? Wikipedia às vezes tem umas besteiras, é sempre bom conferir). Aí eu descrubro que um certo Dr. John Snow desenvolveu um trabalho no qual mapeou surtos de cólera em Londres, e que seu trabalho é considerado um dos precursores dos SIGs. Enfim… Então John Snow sabia de alguma coisa, hein? Rá!)

Mas vamos deixar de papo, e bora trabalhar com alguns exemplos, ok? Este script é basicão, quero mostrar primeiro algumas coisas beeeeem simples. Mas aí, mais pra frente, a gente vai aos poucos criando novos tutoriais e avançando mais um poquin com essa coisa de SIG no R, pois as possibilidades são grandonas.

 


Do que vamos precisar

Bora instalar um monte de pacotes? Pois é, tem vários recursos bacanudos de SIG no R, mas eles estão em vários pacotes diferentes por aí afora. Então, para já ter uma boa parte deles, comece instalando:

  • dismo
  • raster
  • sp
  • maptools
  • rgeos

Alguns deles estão ligados por dependências: tipo, o dismo precisa do raster e do sp, então se você instalá-lo primeiro os outros dois devem vir automaticamente. Mas, na dúvida, confira se você tem todos, beleza?

Ei, já que falei nos pacotes, uma dica legal: confira o tutorial do dismo sobre distribuição potencial de espécies no R! É bem legal, e eu aprendi o basicão de SIG aqui no R por lá, pra depois ir caçar outras fontes.

Quando estiver com tudo pronto, comece pelo começo e carregue a pacotaiada:

library(dismo)
## Loading required package: raster
## Loading required package: sp
library(maptools)
## Checking rgeos availability: TRUE

Viu só o que eu fiz? Como eu sei que rolam ali algumas dependências, carreguei os dois pois sabia que eles iriam carregar os outros sozinhos. Espertinho, né?

Mas vamos que vamos!

 


Meu primeiro mapinha no R!

Booooooooora, galera. Para este tutorial, quero fazer tudinho que eu puder usando apenas dados da memória do R. Então vamos carregar um conjuntinho esperto de dados do maptools que vai nos trazer prontinho de uma vez um shapefile do mundo todo aqui pra gente. Lá vai:

data(wrld_simpl)

Facin, né? Então bora plotar. Uma coisa bacana do R é que ele vai tratar os mapas como gráficos, então os comandos que você já conhece vão funcionar por aqui. Isso é ótimo, pois torna tudo bem mais simples. Ah sim, e o mesmo vale para os mapas que você pode fazer pelo ggplot2! Bom demais, né não?

Então vamos lá, vai ser facin, facin:

plot(wrld_simpl)

Tcharam!!!!!!!

Falei que era fácil!

Mas agora bora brincar um pouquinho. Veja como fica isso aqui, ó:

windows(width=6, height=6, rescale="fixed")

plot(wrld_simpl, xlim=c(-90, -30), ylim=c(-60, 20),
axes=T, col="gray")

E aí, gostou? Então vamos conversar a respeito do que fizemos:

  • a função windows foi usada antes pra abrir uma janelinha gráfica com as dimensões que eu pedi, e com o argumento rescale usado para que a janela não mude o tamanho do gráfico. Você pode brincar com as dimensões da janela para ver o efeito. É bem útil, às vezes.
  • já no “gráfico que é mapa que é gráfico”, eu usei os bons e velhos conhecidos xlim e ylim para determinar os limites do meu mapa. Veja que legal: eles são normalmente usados para colocar os limites dos eixos, mas como aqui os eixos são coordenadas, eu usei valores de latitude e longitude para “cortar” o mapa no pedaço que eu queria. Bacaninha, né?
  • finalmente, axis e col são pura estética: colocar os valores de coordenadas nos eixos e tacar uma corzinha no mapa.

Maneiro, que bom que você gostou. Agora vamos seguir!

 


Eu quero uns pontinhos de ocorrência!

Yeah, vamos lá adicionar a informação sobre a distribuição de uma espécie, vamos? ra isso vamos usar um conjuntinho de dados de exemplo bem manjadão, que aparece por aí em tutoriais do Maxent e do Diva-GIS, e que está dentro da pasta do pacote dismo no seu computador: o arquivo bradypus.csv. Se você não quer fuçar a pasta do pacote, pode pegar ele aqui.

De posso do arquivo (que tem dados de distribuição de bichos preguiça fofinhos), você pode carregá-lo e conferir rapidão o que tem dentro:

dados<-read.csv("bradypus.csv")
summary(dados)
##                 species       dd.long           dd.lat       
##  bradypus_variegatus:116   Min.   :-85.93   Min.   :-23.450  
##                            1st Qu.:-79.20   1st Qu.: -2.896  
##                            Median :-75.08   Median :  4.817  
##                            Mean   :-71.79   Mean   :  2.449  
##                            3rd Qu.:-65.51   3rd Qu.:  9.150  
##                            Max.   :-40.07   Max.   : 13.950

Tá aí. É uma coluna com o nome da espécie (que aqui é meio inútil, uma vez que temos uma espécie só; mas quem já usou o Maxent, por exemplo, sabe que esta coluna é necessária por lá), a longitude e a latitude.

Então, só pra simplificar, eu vou fazer uma malandragem aqui, e jogar fora o nome da espécie. Vou criar um objeto novo, só com as colunas 2 e 3:

bradypus<-dados[,2:3]
summary(bradypus)
##     dd.long           dd.lat       
##  Min.   :-85.93   Min.   :-23.450  
##  1st Qu.:-79.20   1st Qu.: -2.896  
##  Median :-75.08   Median :  4.817  
##  Mean   :-71.79   Mean   :  2.449  
##  3rd Qu.:-65.51   3rd Qu.:  9.150  
##  Max.   :-40.07   Max.   : 13.950

Agora vamos tacar tudo junto no nosso mapa, pode ser?

Se você fechou a janela com o mapinha (coitado, ele tava tão bonitinho, pra que fechar?!), ent]ao vai ter que fazer de novo. Na dúvida, vou colocar o comando todinho aqui, com o mapa outra vez:

windows(width=6, height=6, rescale="fixed")

plot(wrld_simpl, xlim=c(-90, -30), ylim=c(-60, 20),
axes=T, col="gray")

points(bradypus, pch=16, col="red", cex=0.5)

Os comandos aqui não devem ser novidade, mas bora explicar mesmo assim. Usamos a função points, que serve para se plotar pontos em um gráfico qualquer. Aí usamos:

  • pch pro tipo de símbolo;
  • col pra corzinha vermelha bonitinha;
  • cex (com c! Você pensou besteira, né?) pra definir o tamanho da bolinha.

Bacana? Pois vamos avançar mais um pouquinho, pode ser?

 


Diversão com os pontinhos!

Maneiro, fizemos nossos mapas, incluímos os registros de ocorrência e tal e coisa. Mas e aí, o que mais podemos fazer? Uma das coisas bacanas no R quanod usado como SIG é que ele continua sendo o R de sempre, hehe. Então podemos explorar coisas dos nossos dados como sempre, ao mesmo tempo que explorarmos isso do ponto de vista de um mapa, por exemplo.

Imagine, por exemplo, que para cada pontinho de ocorrência dos bichinhos preguiça fofinhos você tenha também outras informações, como alguma estimativa de abundância. Então para cada ponto, além do registro em si você tem um valor do número de indivíduos no local de coleta. Será que podemos incluir isso no nosso mapa? Mas é claro que sim!

Como os dados originais não possuem esta informação, para este exemplo eu vou simplesmente inventar dados de abundância para cada ponto. Tipo, vou inventar mesmo: vamos pedir para o R gerar um conjunto de números inteiros ao acaso, e aí vamos fingir que eles são dados de abundância para continuarmos com o nosso exemplo. Bora, é fácil:

set.seed(1)

abund<-rpois(nrow(dados), 5)

Explicando explicadinho:

  • a função set.seed foi usada pelo seguinte: quando mandamos o R gerar números ao acaso, é claro que a cada vez que mandarmos ele gera números diferentes, certo? Só que aí, para este exemplo, cada um de vocês replicando o script teria um resultado diferente do meu! Como quero vocês olhando pro mesmo mapinha que eu, o set.seed fixa a próxima geração de números dentro de um padrão, e aí todo mundo vai ter um conjuntinho igual.
  • a função rpois gera números ao acaso seguindo a distribuição de Poisson. Eu escolhi esta distribuição porque ela é mesmo a cara de dados de abundância. Dentro da função eu usei a função nrows para o R contar o número de linhas de dados e criar os dados de abundância no tamanho certo; e o 5 ali é o lambda, o parâmetro da distribuição de poisson.
  • ufa, acho que é isso, bora prosseguir.

Vamos fazer o mapa de novo outra vez, mas agora usando pontinhos com tamanho proporcional à abundância de invidíduos em cada localidade:

windows(width=6, height=6, rescale="fixed")

plot(wrld_simpl, xlim=c(-90, -30), ylim=c(-60, 20),
axes=T, col="gray")

points(bradypus, pch=16, col="red", cex=abund/5)

Yeaaaahhhhhh! Deu certo! :D

O que fizemos aqui:

  • Usamos o argumento cex de novo, mas desta vez usando a variável com valores de abundância. Aqui eu dividi o valor por 5, e isto foi feito na base da tentativa e erro mesmo. Tipo, se for apenas o valor puro de abund, as bolinhas ficam grandes demais. Se dividir por 10, ficam pequenas demais. Aí o esquema é ir tentando (teste aí em casa e veja!).

Mas assim, acho que ainda podemos dar aquela melhoradinha básica. Bora?

É que com os pontos com tamanhos variados vem um problema: alguns deles se sobrepõe, e a gente pode estar deixando de enxergar os registros por conta da sobreposição de bolinhas grandes, que podem estar engolindo as menores. Uma solução simples seria tacar um pontinho preto dentro de cada círculo vermelho, de forma a funcionar como um “núcleo” do registro de ocorrência. Aí temos dois elementos: o ponto preto central indica a localidade, e o tamanho da bola vermelha em torno indica a abundância. É super simples, olha:

(sim, fazendo tudo de novo, aiaiai)

windows(width=6, height=6, rescale="fixed")

plot(wrld_simpl, xlim=c(-90, -30), ylim=c(-60, 20),
axes=T, col="gray")

points(bradypus, pch=16, col="red", cex=abund/5)

points(bradypus, pch=16, cex=0.1)

  • Aqui eu apenas dei outro comando points, e ajustei o cex pra um tamanho miúdo. ;)

 


Eu quero mais diversão com SIG no R!

Eu também! Mas agora vamos parar um pouquinho. O R tem mais um monte de coisas pra nos mostrar neste meio, e em algum momento eu volto aqui e mostro para vocês, certinho?

Um abração, e até o próximo script!

Prof Marcos

 


Este documento faz parte do material que disponibilizo no meu blog, o “Cantinho do R”, e foi foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais sobre eles, acesse:

http://www.r-project.org/

http://www.rstudio.com

http://rmarkdown.rstudio.com/

 

Quer comentar ou perguntar? Faça isso no blog, aqui