library(ggmap)
library(maps)
library(tidyverse)O R pode ser usado para fazer mapas!
Nessa aula vamos aprender algumas ferramentas básicas de representações espaciais.
Note
Essa não é minha área então considere apenas uma apresentação das ferramentas possíveis para serem exploradas independentemente!!
maps e ggmap
Podemos plotar mapas de maneira bem simples utilizando o pacote ggmap que segue a mesma estrutura do ggplot2, então todas as opções de customização que já aprendemos podem ser aplicadas
É necessário carregar (instalar previamente, caso não tenha na biblioteca) esses pacotes:
Mapa do mundo
Para gerar um mapa do mundo é BEM simples: basta chamar os dados “world” do pacote ggmap.
Tip
A função coord_quickmap() é necessária para ajustar a projeção!
ggplot()+ borders("world", fill = "grey70", colour = "black")+
coord_quickmap()
Região delimitada do mapa mundial
Para selecionar determinada região do mapa mundial, podemos usar alguns recursos:
Limitar as dimensões para a região desejada
No caso, estou colocando um ponto vermelho onde fica a cidade de Curitiba e tentando limitar apenas para mostrar o Brasil
Warning
Note que com xlim e ylim como camada não fica bom.
ggplot()+ borders("world", fill = "grey70", colour = "black")+
xlim (-100,-20)+ ylim(-45,10)+
geom_point(aes(x=-49.3044, y=-25.4809, colour = "red")) +
theme(legend.position='none')
Warning
Também continua desconfigurado com os limites inseridos na camada principal
ggplot()+ borders("world", fill = "grey70", colour = "black", xlim= c(-100,-20), ylim=c(-45,10))+
geom_point(aes(x=-49.3044, y=-25.4809, colour = "red")) +
theme(legend.position='none')
Tip
É necessário utilizar o xlim e ylim no parâmetro coord_fixed() para ficar certo.
ggplot()+ borders("world", fill = "grey70", colour = "black")+
geom_point(aes(x=-49.3044, y=-25.4809, colour = "red")) +
coord_fixed(xlim= c(-100,-20), ylim=c(-45,10))+
theme(legend.position='none')
Determinar pelo nome da região desejada, utilizando o parâmetro regions
Tip
Da mesma forma, a projeção só fica boa de verdade com coord_fixed()
ggplot()+ borders("world",regions = "Brazil", fill = "grey70", colour = "black")+
geom_point(aes(x=-49.3044, y=-25.4809, colour = "red")) +
theme(legend.position='none')
ggplot()+ borders("world",regions = "Brazil", fill = "grey70", colour = "black")+
geom_point(aes(x=-49.3044, y=-25.4809, colour = "red")) +
coord_fixed(xlim= c(-80,-30), ylim=c(-40,10))+
theme(legend.position='none')
Pacotes para dados ecológicos
Utilizando o pacote rgbif (acesso ao banco de dados da GBIF)
GBIF: Global Biodiversity Information Facility https://www.gbif.org/
- Carregar o pacote rgbif (instalar se vc não tiver na biblioteca)
- Definir a espécie que você quer dando nome científico
- Escolher dados apenas que possuam coordinadas (hasCoordinate=T)
- Pode-se limitar a região por latitude ou longitude ou nomeando o pais pela sigla internacional com 2 letras
- Transformar os dados importados em data frame usando a função as.data.frame()
Ocorrência no mundo todo
Tip
Usei o parâmetro limit para limitar para 500 observações senão demora muito pra carregar
library(rgbif)
ocorrencia<-occ_search(scientificName="Puma concolor", hasCoordinate=T, limit=500)
names(ocorrencia)[1] "meta" "hierarchy" "data" "media" "facets" localizacao <- as.data.frame(ocorrencia$data) %>% select(3:4) %>%
rename(Lat = decimalLatitude, Long=decimalLongitude)
head(localizacao) Lat Long
1 -31.25373 -61.26492
2 34.05975 -116.40297
3 -31.22785 -61.20159
4 37.58867 -122.28376
5 37.14539 -121.59656
6 40.36277 -120.62113Ocorrência no Brasil apenas
ocorrencia.br<-occ_search(scientificName="Puma concolor", country="BR", hasCoordinate=T)
localizacao.br <- as.data.frame(ocorrencia.br$data)%>% select(3:4) %>%
rename(Lat = decimalLatitude, Long=decimalLongitude)Agora a gente pode inserir os pontos de ocorrência num mapa. Podemos usar o borders() ou o que preferir!
ggplot() + borders("world", colour="grey", fill="lightgray") +
theme_bw() +
geom_point(aes(x=Long,y=Lat),data=localizacao,colour="gold",size=2,alpha=0.1) 
ggplot() + borders("world", colour="grey", fill="lightgray") +
theme_bw() +
geom_point(aes(x=Long,y=Lat),data=localizacao,colour="gold",size=2,alpha=0.1)+
coord_fixed(xlim=c(-150,-20)) 
Para ajustar exatamente as medidas do mapa com a distribuição podemos gerar um vetor para restringir as coordenadas, baseado na distribuição máxima e mínima da espécie:
distribuicao<-c(min(localizacao$Long,na.rm=T)-1,
min(localizacao$Lat,na.rm=T)-1,
max(localizacao$Long,na.rm=T)+1,
max(localizacao$Lat,na.rm=T)+1)
ggplot() + borders("world", colour="grey", fill="lightgray") +
theme_bw() +
geom_point(aes(x=Long,y=Lat),data=localizacao,colour="black",size=2,alpha=0.1)+
coord_fixed(xlim=distribuicao[c(1,3)], ylim=distribuicao[c(2,4)])
Projeções
Podemos nos deparar com dados de ocorrência no oceano, ou lugares muito improváveis. Nesse caso, o motivo pode ser diferenças nas projeções.
Projeção é uma novela e eu não entendo NADA!
Só sei que para arrumar existe um pacote próprio para lidar com isso, o rgdal
Com ele, podemos mdificar o CRS (coordinate reference system) de qualquer arquivo de dados espacializado.
O google usa WGS84 - projeção cillíndrica
library(rgdal)
class(localizacao)[1] "data.frame"#converter dataframe para dados espaciais
coordinates(localizacao)<-~Long+Lat
class(localizacao)[1] "SpatialPoints"
attr(,"package")
[1] "sp"Qual projeção usar?
WebMerc<-CRS(“+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs”)
Google<-CRS(“+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0”)
Earth<-CRS(“+proj=longlat +init=epsg:4326”) # Alternate Google Earth Projection
Geo <- CRS(“+proj=longlat +ellps=WGS84 +datum=WGS84”)# Projection from: http://pakillo.github.io/R-GIS-tutorial/
GMaps<-CRS(“+proj=longlat +init=epsg:3857 +datum=WGS84”) # Google Maps Projection
GBIF<-CRS(“+init=epsg:4326”)
summary(localizacao) # perceba que não tem projeção definidaObject of class SpatialPoints
Coordinates:
min max
Long -125.98052 -42.49123
Lat -52.03011 56.11654
Is projected: NA
proj4string : [NA]
Number of points: 500proj4string(localizacao)<-CRS("+init=epsg:4326") #assim vc define a projeção usada pelo GBIF (veja acima)
summary(localizacao)Object of class SpatialPoints
Coordinates:
min max
Long -125.98052 -42.49123
Lat -52.03011 56.11654
Is projected: FALSE
proj4string : [+proj=longlat +datum=WGS84 +no_defs]
Number of points: 500GMaps<-CRS("+proj=longlat +init=epsg:3857 +datum=WGS84")Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
+wktext +no_defsWarning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded datum WGS_1984 in Proj4 definitionloc_proj<-spTransform(localizacao,GMaps) # aqui você transforma a projeção dos dados de localizacao para dados do GMAPS
loc_proj_data <- data.frame(loc_proj)#transformar de volta em dataframe para poder colocar no mapaPacotes para dados geográficos
Alguns pacotes são bem úteis para baixar dados geográficos diretamente, e vamos conhecer alguns deles que possam ter utilidade na área de ecologia
Pacote geobr
Com vários dados oficiais para o Brasil!
library(geobr)
datasets <- list_geobr()
View(datasets)Vamos criar um shapefile com os municípios do estado do Paraná e plotar usando a função geom_sf do ggplot.
muni_PR <- read_municipality(code_muni = "PR",
year=2020, showProgress = FALSE)Using year 2020class(muni_PR)[1] "sf" "data.frame"ggplot() +
geom_sf(data=muni_PR, fill="lightgrey", color="black", show.legend = FALSE) +
theme_minimal()
Também é possivel obter os limites dos biomas com geom_br
biomas <- read_biomes(year=2019, showProgress = FALSE)Using year 2019head(biomas)Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -73.98304 ymin: -33.75115 xmax: -28.84785 ymax: 5.269581
Geodetic CRS: SIRGAS 2000
name_biome code_biome year geom
1 Amazônia 1 2019 MULTIPOLYGON (((-44.08515 -...
2 Caatinga 2 2019 MULTIPOLYGON (((-41.7408 -2...
3 Cerrado 3 2019 MULTIPOLYGON (((-43.39009 -...
4 Mata Atlântica 4 2019 MULTIPOLYGON (((-48.70814 -...
5 Pampa 5 2019 MULTIPOLYGON (((-52.82472 -...
6 Pantanal 6 2019 MULTIPOLYGON (((-57.75946 -...ggplot() +
geom_sf(data=biomas, color="black", aes(fill=name_biome)) +
theme_minimal()
Para tirar o sistema costeiro
ggplot() +
geom_sf(data= slice(biomas,-7), color="black", aes(fill=name_biome)) +
theme_minimal()
Pacote mapview
Muito legal para fazer mapas interativos! Não necessariamente para publicações, mas talvez apresentações e visualizações diversas.
library(mapview)
# mapview(biomas, zcol = "name_biome")# mapview(muni_PR, zcol = "name_muni")Customização de mapas
Agora vamos aprender a deixar os mapas com qualidade para publicação! Alguns requisitos são essenciais nessa etapa:
- Todo mapa DEVE conter uma barra de escala!
- Todo mapa deve conter também a orientação (apontando para o norte)
- Interessante prestar bastante atenção nas cores
- Não representar países que não são ilhas isoladamente!!
ggplot() +
geom_sf(data= slice(biomas,-7), aes(fill=name_biome, alpha=0.7)) +
theme_minimal()+
scale_fill_manual(values = terrain.colors(7))
Agora do jeito correto, usaremos o pacote ggsn para inserir o símbolo apontando para o norte e a escala ao mapa (funções north e scalebar respectivamente).
library(ggsn)Loading required package: gridnorthSymbols() 
Figura final
ggplot()+ borders("world", fill = "lightgray", colour = "gray")+
geom_sf(data= slice(biomas,-7), aes(fill=name_biome, alpha=0.7)) +
scale_fill_manual(values = terrain.colors(7))+
coord_sf(xlim= c(-80,-30), ylim=c(-40,10))+
theme(legend.position = c(0.8,0.9), legend.background = element_blank(), legend.text=element_text(size=10),
legend.title = element_blank(),
legend.key.size = unit(0.4, "cm"),
panel.grid.major = element_line(color = "darkgray", linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue"),
axis.title = element_blank())+
guides(fill = guide_legend(nrow = 3), alpha= "none")+
north(biomas, symbol = 3, location = "bottomright") +
scalebar(biomas, location = "bottomleft", dist = 500, dist_unit = "km",
transform = TRUE, model = "WGS84")
Mapas a partir de shapefiles
Shapefiles (.shp) podem ser encontrados em diversos banco de dados na internet. Alguns comumente utilizados no Brasil são o do IBGE, INPE, ANA e diversos sites do governo. Para dados internacionais, temos vários disponíveis gratuitamente
> confira uma tabela nesse link: https://www.oryxthejournal.org/writing-for-conservation-guide/map-with-a-message.html#tab:data-sources
Pacote terra como subtituto do finado raster
Muito conteúdo BEM completo no site deles https://rspatial.org/pkg/index.html
Vários dados disponíveis no site do INPE http://www.dpi.inpe.br/Ambdata/download.php
Outra opção (que eu achei bem mais pesada) é usar o pacote geodata Útil para dados globais de clima (wordclim), elevação, uso da terra, solo, ocorrência de espécies pelo gbif, limites administrativos, etc.
Opções para dados globais (função worldclim_global), por país (worldclim_country) ou para determinada latitude e longitude (worldclim_tile)
library(geodata)
clim_brasil <- worldclim_country(country= “BR”, var=“bio”, res=10, download=T, path= tempdir(), version = “2.1”)
Vamos baixar apenas a bio1 - Temperatura média pra ir mais rápido, ams se for baixar todas vai vir num arquivo .zip e será nescessário unzip na sua pasta de destino para utilizar.
download.file("http://www.dpi.inpe.br/amb_data/Brasil/bio1_br.asc", destfile="bio1_br.asc")library(terra)
temp <- rast('bio1_br.asc')
class(temp)[1] "SpatRaster"
attr(,"package")
[1] "terra"plot(temp)
Para plotar com ggmap precisa transformar em dataframe
Tip
Para saber as unidades e quais são as variáveis verificar em https://www.worldclim.org/
temp_df <- as.data.frame(temp, xy=T)
ggplot() +
geom_tile(data = temp_df,
aes(x = x, y = y, fill = bio1_br/10))+ #unidade do worldclim é oC*10
theme_bw()
Agora se a gente quiser só para alguma localização específica,temos que editar o arquivo. Para isso temos a função crop
ext <- c(-50,-40, -25,-20) # Para alguma localização específica
temp_x <- crop(temp, ext)
plot(temp_x)
Para o Paraná vamos usar o shape do geobr para cortar o raster no tamanho certo
estado_PR <- read_state(code_state = "PR",
year=2020,
showProgress = FALSE)Using year 2020ggplot() +
geom_sf(data=estado_PR, fill="lightgrey", color="black", show.legend = FALSE) +
theme_minimal()
temp_pr <- crop(temp, estado_PR)
temp_pr_df <- as.data.frame(temp_pr, xy=T)
ggplot() +
geom_tile(data=temp_pr_df, aes(x,y, fill=bio1_br/10))+
theme_minimal()+
coord_equal()
Tip
Atenção para customizações de escala!!! O resultado do seu mapa depende muito disso!
ggplot()+
geom_tile(data=temp_pr_df, aes(x,y, fill=bio1_br/10))+
scale_fill_gradient2(na.value="white", low= "steelblue3", mid="gray94", high="brown3", limits=c(10, 29), midpoint = 18, breaks = c(10, 15,20, 25))+
geom_sf(data=estado_PR, fill=NA, color="black", show.legend = FALSE) +
theme(axis.text = element_blank(), axis.ticks = element_blank(),
axis.title = element_blank(), panel.background = element_blank())
Alguns sites/livros úteis para aprofundar em R espacial
https://www.ecologi.st/spatial-r/
https://tmieno2.github.io/R-as-GIS-for-Economists/
https://damariszurell.github.io/EEC-MGC/a1_SpatialData.html