2_Aula Prática de Mapas

Pacotes necessários

library(rnaturalearth)
library(sf)
library(terra)
library(dplyr)
library(ggsn)
library(ggmap)
library(ggplot2)
library(ggspatial)
library(geodata)
library(rgbif)
library(sp)
library(tmap)

Criando mapas com o R

ggplot

Mapas podem ser visualizados através do ggplot, da mesma maneira que fazemos para criar outros tipos de gráficos. Vamos carregar o shapefile do neotrópico, como aprendemos na aula passada. Depois, podemos usar a função geom_sf() para plotar ele com ggplot:

neot<-st_read("dados/shapes/neotropica.shp")

Reading layer `neotropica' from data source 
  `C:\Users\amabi\OneDrive\Área de Trabalho\Aula_mapas\dados\shapes\neotropica.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 49 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -151.4978 ymin: -58.49861 xmax: -26.24139 ymax: 32.71846
CRS:           NA

ggplot(neot)+
  geom_sf()

Se quisermos adicionar uma nova camada vetorial acima desta, podemos adicionar um novo geom_sf abaixo do primeiro. Vamos utilizar a função ne_countries que aprendemos como utilizar na aula passada para obter um shapefile do Brasil. Lembrem que precisamos ter o mesmo crs nas duas camadas, entao vamos atribuir o crs do objeto bra (WGS84, e o código é 4326) ao objeto neot. Então, vamos adicioná-lo ao nosso plot original com uma coloração diferente:

bra<-ne_countries(country="Brazil")
st_crs(neot) = 4326

ggplot(neot)+
  geom_sf()+
  geom_sf(data=bra, fill="darkgreen")

Para adicionar pontos de ocorrência, utilizamos o geom_sf. Vamos adicionar nossos pontos de ocorrência de Cyclodium recortados para o Brasil:

occs<-read.csv("dados/occs/cyclodium.csv")
occs<-st_as_sf(occs, coords = c("long", "lat"))
st_crs(occs)<-4326
occs_bra<-st_intersection(occs, bra)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

ggplot(neot)+
  geom_sf()+
  geom_sf(data=bra, fill="darkgreen")+
  geom_sf(data=occs_bra)

Podemos mudar alguns detalhes do nosso plot, vamos deixar o neotrópico com fundo branco e linhas tracejadas, mudar a cor do Brasil e deixá-lo sem bordas e trocar o tipo e tamanho dos pontos de ocorrência:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill=NA, linetype="dashed")+
  geom_sf(data=bra,fill="#c2a5cf", linetype=0)+
  geom_sf(data=occs_bra, colour="#40004b", pch=8, size=2)

Já que todas as nossas camadas são do tipo sf, podemos usar geom_sf para nos aproximarmos de uma região específica do mapa (“dar um zoom”). Para isso, utilizamos o coord_sf. No argumento xlim e ylim, colocamos os números em longitude e latitude dos limites onde queremos “recortar” a visualização. Vamos fazer esse recorde para nos aproximarmos da área do Brasil. Além disso, todas as espécies estão aparecendo juntas no nosso mapa. Vamos usar a função facet_wrap para criar um mapa para cada espécie:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill=NA, linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf", alpha=0.5)+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  facet_wrap(~sp)

Alternativamente, podemos fazer essa separação dos pontos de ocorrência por cores e/ou símbolos diferentes. Para fazer uma seleção por cores, vocês precisam colocar dentro de aes(colour=sp) e por tipos de símbolos aes(pch=sp). Também adicionamos um novo argumento em scale_color_hue(name=“Espécie”) para mudar o título da legenda. Isso também pode ser feito diretamente no dataframe:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+ 
  scale_color_hue(name="Espécie")

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(pch=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))

Warning: The shape palette can deal with a maximum of 6 discrete values because more
than 6 becomes difficult to discriminate
ℹ you have requested 7 values. Consider specifying shapes manually if you need
  that many have them.

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_sf()`).

Também podemos mexer no tema do plot. Aqui, temos o plot padrão do ggplot, mas vocês podem brincar com o argumento theme para conferir outros temas pré-determinados. Vou dar dois exemplos abaixo:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+ 
  scale_color_hue(name="Espécie")+
  theme_bw()

O theme_bw() deixa o fundo branco e as linhas das coordenadas em cinza, além de incluir uma delimitação de uma caixa em volta do mapa. Outra possibilidade:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+ 
  scale_color_hue(name="Espécie")+
  theme_ps()

O theme_ps() do pacote ggmap remove as linhas de grade e coordenadas e também a legenda. É uma opção para salvar o mapa separadamente para apresentações.

Também podemos incluir um título geral ao nosso gráfico através do argumento title com a função labs. Além disso, vamos incluir um fundo azul para simular as regiões de mar, dentro da função theme, utilizando panel_background:

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  theme_bw()+ 
  scale_color_hue(name="Espécie")+
  labs(title="Distribuição das espécies de Cyclodium no Brasil")+
  theme(panel.background = element_rect(fill="lightblue"))

Notem que as linhas brancas ficam bem visíveis com esse fundo, vamos corrigir isso trocando a cor por cinza (com o argumento panel_grid) e diminuir sua espessura (com o argumento linewidth):

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  theme_bw()+ 
  scale_color_hue(name="Espécie")+
  labs(title="Distribuição das espécies de Cyclodium no Brasil")+
  theme(panel.background = element_rect(fill="lightblue"), panel.grid=element_line(colour="grey", linewidth = 0.2))

Por fim, vamos mexer na legenda. Vamos tirar esse fundo azul que foi incluído ao pintarmos o fundo do painel (com o argumento legend.key) e colocar os nomes em itálico, já que são nomes científicos (com o argumento legend.text):

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  theme_bw()+
  labs(title="Distribuição das espécies de Cyclodium no Brasil")+
  scale_color_hue(name="Espécie")+
  theme(panel.background = element_rect(fill="lightblue"),  panel.grid=element_line(colour="grey", linewidth = 0.2),legend.key=element_rect(fill="white"),legend.text = element_text(face="italic"))

Por fim, vamos adicionar dois elementos muito importantes em um mapa: a indicação do norte e a escala. Para isso vamos usar duas funções do pacote ggspatial (annotation_north_arrow e annotation_scale):

ggplot()+
  geom_sf(data=neot, colour="darkgrey",fill="lightgrey", linetype="dashed")+
  geom_sf(data=bra, fill="#c2a5cf")+
  geom_sf(data=occs_bra, aes(colour=sp))+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  theme_bw()+
  labs(title="Distribuição das espécies de Cyclodium no Brasil")+
  scale_color_hue(name="Espécie")+
  theme(panel.background = element_rect(fill="lightblue"),  panel.grid=element_line(colour="grey", linewidth = 0.2),legend.key=element_rect(fill="white"),legend.text = element_text(face="italic"))+
  annotation_north_arrow(style=north_arrow_nautical, location="tr", height=unit(2, "cm"), width = unit(2, "cm"))+annotation_scale(location="br", bar_cols= c("white", "black"), line_width=0.5)

Scale on map varies by more than 10%, scale bar may be inaccurate

Ao consultar o help de cada uma das funções, podemos ver mais detalhes do que podemos estilizar do símbolo de norte e da escala. Para o símbolo de norte, em style, escolhemos que tipo de símbolo queremos (há mais 3 possibilidades que podem ser escolhidas). Em location=“tr”, dissemos que queremos nosso símbolo no top-right (as outras opções são tl, bl, e br, respectivamente top-left, bottom-left e bottom-right). E em height e width, nós ajustamos o tamanho do símbolo. Para escolher a localização da escala, utilizamos da mesma forma o argumento location, nesse caso, escolhemos o bottom-right. Também indicamos a cor da caixa com bar_cols e deixamos as bordas caixa mais finas com line_width. Explorem as possibilidades!

Agora, vamos ver como podemos plotar um raster. Nesse caso, nosso raster foi carregado com o pacote terra, já que o pacote sf funciona apenas para dados vetoriais. Dessa forma, o argumento que utilizamos para abrir um raster é o geom_raster. Precisamos recortar novamente nosso raster para os limites do Brasil com crop e raster. Para isso, vamos carregar um shapefile do Brasil como spatVector (com ne_countries, returnclass=“sv” - spatVector, por padrão ele baixa um arquivo sf, então precisamos mudar aqui). Também precisamos transformar esse raster em um dataframe para plotar os dados, com o as.data.frame do pacote terra, indicando que é um dado espacial com xy=T. Nesse caso, ao plotar o geom_raster, precisamos indicar no aes qual é o valor de longitude (x), latitude (y), e o fill (valor do raster) para que possamos ver as informações de temperatura média anual do ar:

bio1<-rast("dados/rasters/wc2.1_10m_bio_1.tif")
plot(bio1)

bra_terra<-ne_countries(country="Brazil", returnclass = "sv")
bio1_bra<-terra::crop(bio1, bra_terra)
bio1_bra<-terra::mask(bio1_bra, bra_terra)
bio1_bra_df<-terra::as.data.frame(bio1_bra, xy=T)
  
ggplot()+
  geom_sf(data=neot)+
  geom_raster(data=bio1_bra_df, aes(x=x, y=y, fill=wc2.1_10m_bio_1))

Por padrão, as cores mostradas no geom_raster é esse gradiente de azul. Podemos mudar isso ao utilizar a função scale_fill_gradient2(). Vamos usar essa opção, pois ela permite que escolhamos a cor de valores baixos (no nosso caso, temperaturas mais frias), do valor mediano (temperaturas intermediárias) e dos valores altos (temperaturas altas). Vamos ver como ficaria e mudar nosso tema para theme_bw:

ggplot()+
  geom_sf(data=neot)+
  geom_raster(data=bio1_bra_df, aes(x=x, y=y, fill=wc2.1_10m_bio_1))+
  scale_fill_gradient2(na.value = "white", low="steelblue3", mid="gray94", high="brown3", limits=c(10,30), midpoint=20, breaks=c(10, 20, 30))+
  theme_bw()

Agora, vamos aproximar a visualização para o Brasil (com coord_sf) e mudar o nome do título da legenda (ao adicionar o argumento name no scale_fill_gradient2), e os nomes do eixo x (xlab) e y (ylab):

ggplot()+
  geom_sf(data=neot)+
  geom_raster(data=bio1_bra_df, aes(x=x, y=y, fill=wc2.1_10m_bio_1))+
  scale_fill_gradient2(na.value = "white", low="steelblue3", mid="gray94", high="brown3", limits=c(10,30), midpoint=20, breaks=c(10, 20, 30),name="Temperatura média anual")+
  theme_bw()+
  coord_sf(xlim=c(-80,-35), ylim=c(-40,10))+
  xlab("Longitude")+
  ylab("Latitude")

Por fim, para salvar cada um dos mapas criados com ggplot, utilizamos o ggsave. Podemos criar um objeto para cada plot final e salvar o objeto com plot=nome_do_objeto, ou com plot=last_plot(), para salvar o último plot carregado. Vamos salvar em png (ao incluir .png no nome do arquivo), mudar o tamanho para 10 cm de largura e 6 de altura, com um dpi de 300 (geralmente o mínimo que revistas pedem):

ggsave("mapa_ggplot.png", plot=last_plot(), width=10, height=6, dpi=300)

tmap

Agora, vamos ver como fazer mapas com o pacote tmap, que é uma alternativa ao ggplot e específico para plotar mapas. Vamos carregar os dados de estados do Brasil para entender como as funções funcionam. Além disso, precisamos usar a função st_make_valid, pois as funções do tmap não conseguem ler esse sf. Para fazer um plot simples com o tmap, precisamos primeiro indicar que estamos incluindo uma camada com tm_shape e depois dizer que tipo de plotagem queremos com tm_borders, tm_fill ou tm_polygons. Vamos ver a diferença de cada um deles. O tm_borders mostra as bordas do shapefile com todas as subdivisões.

estados<-ne_states(country="Brazil")
estados<-st_make_valid(estados)

tm_shape(estados)+
  tm_borders()

O tm_fill preenche toda a extensão do shapefile:

tm_shape(estados)+
  tm_fill()

O tm_polygons identifica cada um dos polígonos, mantendo as bordas e preenchendo eles. É uma mistura do borders com o fill:

tm_shape(estados)+
  tm_polygons()

Para adicionar mais de um shapefile ao mesmo tempo, precisamos repetir o tm_shape e o tipo de polígono que queremos, um abaixo do outro. Lembrem que aqui a ordem importa, então, se queremos que o shapefile do neotrópico fique abaixo do shapefile dos estados do Brasil, precisamos adicioná-lo primeiro:

tm_shape(neot)+
  tm_borders()+
tm_shape(estados)+
  tm_polygons()

Podemos customizar as cores e bordas de cada uma das camadas, assim como fizemos no ggplot. Aqui, vamos fazer essas modificações dentro do tm_borders (para customizar o neotrópico) e tm_polygons (para customizar os estados do Brasil):

tm_shape(neot)+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn")

Aqui, nós mudamos a cor (col) e tipo de linha (lty) da camada do neotrópico em tm_borders e indicamos a coluna de estados (“name”) para serem pintadas através da paleta especificada (em palette).

Para dar um zoom na área do Brasil, podemos usar o argumento bbox na primeira camada de shapefile. Vamos também mudar o título da nossa legenda de name para estados, incluindo o argumento title no tm_polygons do estados:

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")

Dica: Existem várias paletas possíveis para serem utilizadas. Se vocês usarem a função tmaptools::palette_explorer(), uma janela abrirá com todas as opções que podem ser utilizadas.

Podemos adicionar os nossos pontos de ocorrência, adicionando um novo tm_shape e em seguida tm_dots. Dentro de tm_dots, vamos pedir para ele mudar o tipo (shape), tamanho (size) e cor (col) dos nossos pontos:

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape=21, size=0.5, col="white")

Podemos também mudar a cor e símbolo pelo tipo de espécie, como fizemos com o ggplot, mudando dentro de tm_dots os argumentos col e shape = sp.

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape="sp", size=0.5, col="sp")

Warning: Number of levels (unique values) is 7, which is larger than number of
symbol shapes (5).

Além disso, também podemos criar facetas para separar um mapa por espécie como fizemos com facet_wrap no ggplot

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape=21, size=0.5, col="white")+
  tm_facets(by="sp")

O tmap possui argumentos próprios para adicionar o símbolo de norte (tm_compass) e escala (tm_scale_bar):

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape=21, size=0.5, col="white")+
tm_compass()+
tm_scale_bar()

Scale bar set for latitude km and will be different at the top and bottom of the map.

Explorem o help de cada uma dessas funções para conferir os tipos de customizações que podem ser feitos. Vamos mexer no tipo de símbolo (type), tamanho (size) e posição (position) do símbolo de norte e o diminuir o número da escala para 1000 (breaks):

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape=21, size=0.5, col="white")+
tm_compass(type="rose", size=2, position=c(0.9, 0.85))+
  tm_scale_bar(breaks=1000)

Warning: First scale_bar breaks value should be 0.

Scale bar set for latitude km and will be different at the top and bottom of the map.

Também podemos ver esse mapa de forma interativa, ao inserir tmap_mode(“view”) antes de plotar nosso mapa. Lembrem-se de utilizar tmap_mode(“plot”) novamente para sair do modo interativo:

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(neot, bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_borders(col="grey", lty="dashed")+
tm_shape(estados)+
  tm_polygons("name",palette="RdYlGn", title="Estados")+
tm_shape(occs_bra)+
  tm_dots(shape=21, size=0.5, col="white")+
tm_compass(type="rose", size=2, position=c(0.9, 0.85))+
  tm_scale_bar(breaks=1000)

Compass not supported in view mode.

Warning: In view mode, scale bar breaks are ignored.

Por fim, vamos ver como podemos adicionar camadas raster ao nosso mapa tmap. Vamos carregar uma base de dados de elevação do mundo:

data("land")

Notem que nessa base de dados, temos 4 itens dentro de uma lista. Os dados de elevação estão dentro da camada “elevation”. Vamos ter que especificar isso ao carregar o raster:

tm_shape(land)+
  tm_raster("elevation")

Variable(s) "elevation" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Esse é o padrão de cores que o tmap tem ao plotar informações de raster. Podemos trocar isso ao incluir uma outra paleta e pedir para a legenda não aparecer:

tm_shape(land)+
  tm_raster("elevation", palette=terrain.colors(10), legend.show = F)

Vamos dar um zoom na nossa área de interesse e adicionar mais camadas de shapefiles acima deste raster. Além disso, vamos incluir a escala e o símbolo de norte:

tm_shape(land,bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_raster("elevation", palette=terrain.colors(10), legend.show = F)+
tm_shape(neot)+
  tm_borders(col="grey")+
  tm_shape(estados)+
  tm_borders(col="black")+
tm_compass(type="rose", size=2, position=c(0.9, 0.85))+
  tm_scale_bar(breaks=1000)

Compass not supported in view mode.

Warning: In view mode, scale bar breaks are ignored.

Por fim, para salvar nosso mapa, utilizamos o função tmap_save. Nesse caso, vocês precisam criar um objeto com o seu mapa pronto para pedir para salvá-lo. Os mesmos argumentos do ggsave podem ser utilizados aqui.

mapa<-tm_shape(land,bbox = c(xmin = -95, xmax = -34, ymin = -33, ymax = 10))+
  tm_raster("elevation", palette=terrain.colors(10), legend.show = F)+
tm_shape(neot)+
  tm_borders(col="grey")+
  tm_shape(estados)+
  tm_borders(col="black")+
tm_compass(type="rose", size=2, position=c(0.9, 0.85))+
  tm_scale_bar(breaks=1000)

tmap_save(mapa, filename="mapa_tmap_altitude.png", width = 10, height= 6, units = "cm", dpi=300)

Warning: First scale_bar breaks value should be 0.

Scale bar set for latitude km and will be different at the top and bottom of the map.

Map saved to C:\Users\amabi\OneDrive\Área de Trabalho\Aula_mapas\mapa_tmap_altitude.png

Resolution: 1181.102 by 708.6614 pixels

Size: 3.937008 by 2.362205 inches (300 dpi)

Lembrem que essa é uma aula bem introdutória. Apesar de ter bastante conteúdo e muitos detalhes, ainda existem inúmeras customizações que podem ser feitos utilizando o ggplot e o tmap. Além disso, também existem outros pacotes especializados em mapas, que não abordamos aqui! Espero que tenham ficado inspirados a construir mapas no R!