library(rnaturalearth)
library(sf)
library(terra)
library(dplyr)
library(ggsn)
library(ggmap)
library(ggplot2)
library(ggspatial)
library(geodata)
library(rgbif)
library(sp)
library(tmap)Aula Prática de Mapas
Pacotes necessários
Abrindo dados espaciais
Dados Próprios
Vamos primeiro aprender a abrir dados que vocês possam ter e queiram utilizar para construir os mapas:
Vetores
Vocês podem abrir dados de vetores através do pacote sf ou terra e plotar com o r base para conseguir observar os dados. O sf é um pacote mais específico para dados vetoriais, enquanto o terra pode ser utilizado para diferentes tipos de dados espaciais.
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: NAterra::crs(neot)[1] ""plot(neot)Warning: plotting the first 9 out of 11 attributes; use max.plot = 11 to plot
all
Um segundo exemplo com os biomas do Brasil:
biomas<-st_read("dados/shapes/lm_bioma_250.shp")Reading layer `lm_bioma_250' from data source
`C:\Users\amabi\OneDrive\Área de Trabalho\Aula_mapas\dados\shapes\lm_bioma_250.shp'
using driver `ESRI Shapefile'
Simple feature collection with 6 features and 0 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -73.98318 ymin: -33.75118 xmax: -28.84777 ymax: 5.269581
Geodetic CRS: SIRGAS 2000terra::crs(biomas)[1] "GEOGCRS[\"SIRGAS 2000\",\n DATUM[\"Sistema de Referencia Geocentrico para las AmericaS 2000\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore.\"],\n BBOX[-59.87,-122.19,32.72,-25.28]],\n ID[\"EPSG\",4674]]"plot(biomas)
Notem que com o sf abre nossos dados como um dataframe e a última coluna mostra a geometria do vetor, possibilitando que ele seja visualizado. Notem também, que por conta desse formato, quando plotamos o vetor carregado com o st_read, todos os atributos (colunas) geram um mapa separado. Para plotar apenas um dos atributos, podemos usar o plot(neot$nomedacoluna) para mostrar apenas um deles. No caso dos biomas, apenas as informações de geometria estão presentes no arquivo original e isso não acontece.
Também podemos abrir dados vetorias com o R, através da função vect:
neot_terra<-vect("dados/shapes/neotropica.shp")
terra::crs(neot_terra)[1] ""plot(neot_terra)
E com os biomas do Brasil
setwd("C:/Users/amabi/OneDrive/Área de Trabalho/Aula_mapas/")
biomas_terra<-vect("dados/shapes/lm_bioma_250.shp")
terra::crs(biomas_terra)[1] "GEOGCRS[\"SIRGAS 2000\",\n DATUM[\"Sistema de Referencia Geocentrico para las AmericaS 2000\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore.\"],\n BBOX[-59.87,-122.19,32.72,-25.28]],\n ID[\"EPSG\",4674]]"plot(biomas_terra)
Já o terra abre o vetor na extensão SpatVector, que tem características diferentes. Ao plotar, podemos ver que ele não identifica as colunas como atributos diferentes e plota apenas as linhas do vetor. Por fim, nós usamos o crs do pacote terra para conferir o sistema de coordenadas do shapefile. Nesse caso, o shape neotropica não tem nenhum csr associado e precisaremos lidar com isso no futuro. Já o shape dos biomas foi construído no sistema Sirgas 2000.
Rasters
Para abrir um raster da temperatura média anual do ar, vamos utilizar a função rast do pacote terra:
setwd("C:/Users/amabi/OneDrive/Área de Trabalho/Aula_mapas/")
bio1<-rast("dados/rasters/CHELSA_bio1_1981-2010_V.2.1.tif")
terra::crs(bio1)[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"plot(bio1)
A função rast do terra abre o raster como um SpatRaster. O CSR desse arquivo é o WGS84.
Pontos de Ocorrência
Pontos de ocorrência são geralmente arquivos de texto que possuem, ao menos, o nome da espécie/localidade específica (parcelas, por exemplo), e a longitude e latitude para cada uma das observações (linhas).
occs<-read.csv("dados/occs/cyclodium.csv")
head(occs) sp long lat
1 Cyclodium akawaiorum -59.83 5.08
2 Cyclodium akawaiorum -59.83 5.08
3 Cyclodium akawaiorum -59.84 5.10
4 Cyclodium akawaiorum -59.88 5.07
5 Cyclodium akawaiorum -59.88 5.07
6 Cyclodium akawaiorum -59.88 5.05Como este não é um dado espacial, o pacote sf possui uma função chamada st_as_sf, que transforma alguns tipos de dados em sf, fazendo com que as informações de longitude e latitude se tornem geometrias (chequem a última coluna do df):
occs<-st_as_sf(occs, coords = c("long", "lat"), crs = 4326)
head(occs)Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -59.88 ymin: 5.05 xmax: -59.83 ymax: 5.1
Geodetic CRS: WGS 84
sp geometry
1 Cyclodium akawaiorum POINT (-59.83 5.08)
2 Cyclodium akawaiorum POINT (-59.83 5.08)
3 Cyclodium akawaiorum POINT (-59.84 5.1)
4 Cyclodium akawaiorum POINT (-59.88 5.07)
5 Cyclodium akawaiorum POINT (-59.88 5.07)
6 Cyclodium akawaiorum POINT (-59.88 5.05)Dados Externos
Também é possível obter alguns dados a partir do R, caso você não tenha as camadas espaciais.
Vetores
Uma ótima opção alternativa para obter dados de vetores é o pacote rnaturalearth. Ao carregar o pacote com library e escrevendo rnaturalearth:: e apertando tab, vocês podem conferir que tipos de dados estão disponíveis para serem obtidos. Aqui, como exemplo, vamos baixar uma camada de países do mundo e estados brasileiros:
Países, com ne_countries:
paises<-ne_countries(scale="medium", returnclass = "sf")
plot(paises$geometry)
Estados brasileiros, com ne_states
estados<-ne_states(country="Brazil", returnclass = "sf")
plot(estados$geometry)
Rasters
Há pouco tempo, o pacote raster era uma das possibilidades para obtermos dados externos. É possível que vocês encontrem muito material na internet oferecendo essa opção. Porém, o pacote foi descontinuado e muitas funções já não estão funcionando mais. Uma ótima alternativa é o pacote geodata. Aqui, vocês podem ir no help do pacote e olhar quais tipos de dados estão disponíveis para baixar. Vamos baixar uma camada de solo do SoilGrids através desse pacote:
soil<-soil_world(var="ocd", depth=30, path="C:/Users/amabi/OneDrive/Área de Trabalho/Aula_mapas/dados/rasters")
plot(soil)Warning: TIFFFillStrip:Read error at scanline 8094; got 0 bytes, expected 10665
(GDAL error 1)Warning: TIFFReadEncodedStrip() failed. (GDAL error 1)Warning: C:/Users/amabi/OneDrive/Área de
Trabalho/Aula_mapas/dados/rasters/soil_world/ocd_15-30cm_mean_30s.tif, band 1:
IReadBlock failed at X offset 0, Y offset 8133: TIFFReadEncodedStrip() failed.
(GDAL error 1)
Cada tipo de dado vai ter argumentos diferentes, então é necessário sempre conferir o help da função para saber o que é necessário incluir. Nesse caso, falamos qual a informação do solo que queremos (var=“ocd” - densidade de carbono do solo), a profundidade de solo (depth=30) e onde queremos salvar esses dados no nosso computador (path=..)
Pontos de Ocorrência
Por fim, também é possível obter dados de ocorrência para várias espécies a partir de diferentes pacotes no R. Aqui, vou mostrar a opção do pacote rgbif, que faz download das ocorrências da base de dados do GBIF.
araucaria<-occ_search(scientificName = "Araucaria angustifolia", hasCoordinate = T)
head(araucaria$data)# A tibble: 6 × 105
key scientificName decimalLatitude decimalLongitude issues datasetKey
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 4512106571 Araucaria angus… -29.2 -50.1 cdc,c… 50c9509d-…
2 4512231535 Araucaria angus… -22.4 -43.2 cdc,c… 50c9509d-…
3 4510121216 Araucaria angus… -26.5 -52.6 cdc,c… 50c9509d-…
4 4599988303 Araucaria angus… -22.7 -45.6 cdc,c… 50c9509d-…
5 4936471725 Araucaria angus… -22.7 -45.6 cdc,c… 50c9509d-…
6 4510379393 Araucaria angus… -29.4 -50.9 cdc,c… 50c9509d-…
# ℹ 99 more variables: publishingOrgKey <chr>, installationKey <chr>,
# hostingOrganizationKey <chr>, publishingCountry <chr>, protocol <chr>,
# lastCrawled <chr>, lastParsed <chr>, crawlId <int>, basisOfRecord <chr>,
# occurrenceStatus <chr>, taxonKey <int>, kingdomKey <int>, phylumKey <int>,
# classKey <int>, orderKey <int>, familyKey <int>, genusKey <int>,
# speciesKey <int>, acceptedTaxonKey <int>, acceptedScientificName <chr>,
# kingdom <chr>, phylum <chr>, order <chr>, family <chr>, genus <chr>, …occs_gbif <-as.data.frame(araucaria$data)
occs_gbif<-occs_gbif %>%
dplyr::select(scientificName, decimalLongitude, decimalLatitude) %>%
dplyr::rename(sp=scientificName, lat=decimalLatitude, long=decimalLongitude)
head(occs_gbif) sp long lat
1 Araucaria angustifolia (Bertol.) Kuntze -50.06277 -29.16073
2 Araucaria angustifolia (Bertol.) Kuntze -43.18344 -22.39742
3 Araucaria angustifolia (Bertol.) Kuntze -52.56508 -26.46005
4 Araucaria angustifolia (Bertol.) Kuntze -45.57525 -22.70443
5 Araucaria angustifolia (Bertol.) Kuntze -45.57405 -22.70446
6 Araucaria angustifolia (Bertol.) Kuntze -50.87736 -29.38844O occ_search retorna uma lista de listas, com diferentes informações. O que nos interessa são os dados que estão arquivados no $data e as colunas de espécie, longitude e latitude. Então já arrumamos os dados.
Manipulando Dados Espaciais
Muitas vezes, os dados espaciais que obtemos precisam ser manipulados para o que estamos mais interessados em demonstrar. Nessa seção, vamos aprender como fazer algumas dessas etapas.
Recortar vetores
Quando temos um arquivo sf, já que ele é um dataframe, podemos simplesmente utilizar as funções do dplyr para filtrar os dados que queremos:
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: NAbra<-dplyr::filter(neot, CNTRY_NAME=="Brazil")
plot(bra$geometry)
Com o SpatVector carregado pelo terra, podemos fazer uma intersecção de um shapefile com outro (neste caso, ambos precisam ser um SpatVector) para fazer essa seleção. Então vamos pegar o objeto bra (sf) que criamos anteriormente e criar outro objeto, transformando-o em spatVect com a função vect:
neot_terra<-vect("dados/shapes/neotropica.shp")
bra_vect<-vect(bra)
bra_terra<-terra::intersect(neot_terra, bra_vect)
plot(bra_terra)
Agregar vetores (remover as subdivisões)
Algumas vezes, podemos obter shapefiles que possuem algumas subdivisões (por exemplo, o Brasil com as subdivisões dos estados), mas queremos apenas um contorno específico. Podemos realizar isso com os dados do sf, através da manipulação do dataframe, juntamente à função st_union:
plot(neot$geometry)
neot_ag<-neot %>%
summarise(geometry=st_union(geometry))
plot(neot_ag)
Também existe uma função específica para isso no pacote terra, com a função aggregate:
plot(neot_terra)
neot_vect_ag<-terra::aggregate(neot_terra)
plot(neot_vect_ag)
Recortar Rasters (crop e mask)
Para recortar um raster, precisamos de um vetor que delimite a área de interesse. Esse recorte pode ser feito de duas maneiras. Vamos cortar o nosso bio1 com o shapefile do Brasil utilizando a função crop, do terra:
bio1_bra<-crop(bio1, bra_terra)
plot(bio1_bra)
Notem que o crop utiliza os valores de extensão (xmin, xmax, ymin, ymax) para fazer o recorte. Por isso, os dados de outros países da América do Sul e também do mar também aparecem aqui.
Para realizar um recorte através da geometria do vetor do Brasil, utilizamos a função mask, do terra:
bio1_bra_mask<-mask(bio1_bra, bra_terra)Warning: [mask] CRS do not matchplot(bio1_bra_mask)
Mudando a resolução de um raster
Como comentei, rasters são uma imagem composta por uma matriz de pixels. O tamanho destes pixels pode ser alterado dependendo da sua necessidade. Pixels maiores possuem baixa resolução e pixels menores possuem alta resolução. Nosso bio1 é um raster de 30 segundos, isso quer dizer que cada pixel tem, aproximadamente, 1km. Vamos utilizar um raster com resolução de 10km para mudar a nossa resolução:
rast<-rast("dados/rasters/wc2.1_10m_bio_1.tif")
plot(rast)
Esse raster em 10km contém informações de todo o mundo, então vamos utilizar crop e mask para deixá-lo igual ao bio1 de 1km:
rast_bra<-crop(rast, bra_terra)
rast_bra<-mask(rast_bra, bra_terra)Warning: [mask] CRS do not matchplot(rast_bra)
E agora podemos utilizar a função resample, do pacote terra, para mudar o tamanho das células de bio1 (1km) para 10km:
raster_10<-resample(bio1_bra_mask, rast_bra, method="bilinear")
plot(raster_10)
Aqui podemos ver uma mudança sutil no tamanho dos pixels, mas vamos criar uma outra opção. Vocês podem criar um raster com o tamanho desejado para fazer a mudança da resolução do seu raster original. Vamos criar um raster que tenha apenas 50 linhas e 50 colunas, com a mesma extensão do nosso raster anterior.
raster<-rast(ncol=50, nrow=50, ext=ext(bio1_bra_mask))
plot(raster)
Notem que ao tentar plotar esse objeto SpatRaster, nada aparecerá. Isso é porque não demos nenhuma informação de dados para cada pixel, só definimos o tamanho dos pixels. Agora, podemos usar ele para mudar a resolução do nosso raster original:
raster_res<-resample(bio1_bra_mask, raster, method="bilinear")
plot(raster_res)
Neste caso, conseguimos ver bem mais o efeito do tamanho dos pixels.
Recortar ocorrências
Por fim, às vezes se faz necessário retirar alguns pontos de ocorrência que não estão dentro da nossa área de interesse. Vamos dar uma olhada nos nossos pontos de ocorrência de araucárias que baixamos anteriormente do GBIF:
occs_gbif_sp<-st_as_sf(occs_gbif, coords= c("long", "lat"))
plot(occs_gbif_sp$geometry)
O quarto não permite o plot simultâneo do shapefile do mundo, mas conseguimos ver uma concentração de pontos no neotrópico e alguns outros espalhados pelo mundo. Vamos recortar nossos pontos apenas para o Brasil, através do st_intersection:
occs_recortada<-sf::st_intersection(occs_gbif_sp, bra)Warning: attribute variables are assumed to be spatially constant throughout
all geometriesplot(occs_recortada$geometry)
Pronto, agora temos apenas os pontos de ocorrência de Araucárias para o Brasil. O mesmo pode ser feito com os dados no terra, mas precisamos transformar nossos dados de ocorrência para spatVector, com a função vect. Depois utilizamos intersect do terra para realizar a remoção dos pontos fora do Brasil:
occs_recortar<-vect(occs_gbif_sp)
occs_recortada_terra<-terra::intersect(occs_recortar, bra_terra)
plot(occs_recortada_terra)
Ótimo! Agora vocês sabem como carregar dados próprios, obter dados externos pelo R e o básico de manipulação os dados espaciais. Na próxima aula, vamos explorar como podemos plotar os mapas de forma mais elegante!
