Capítulo 15 Dados geoespaciais
Pré-requisitos do capítulo
Pacotes e dados que serão utilizados neste capítulo.
## Pacotes
library(ecodados)
library(here)
library(tidyverse)
library(sf)
library(raster)
library(rgdal)
library(spData)
library(rnaturalearth)
library(geobr)
library(ggplot2)
library(ggspatial)
library(tmap)
library(tmaptools)
library(grid)
library(mapview)
library(leaflet)
library(viridis)
library(knitr)
library(sidrar)
library(landscapetools)
## Dados
# world <- world
# volcano <- volcano
# geo_anfibios_locais <- ecodados::geo_anfibios_locais
# geo_anfibios_especies <- ecodados::geo_anfibios_especies
# geo_vetor_nascentes <- ecodados::geo_vetor_nascentes
# geo_vetor_hidrografia <- ecodados::geo_vetor_hidrografia
# geo_vetor_cobertura <- ecodados::geo_vetor_cobertura
# geo_vetor_rio_claro <- ecodados::geo_vetor_rio_claro
# geo_vetor_brasil <- ecodados::geo_vetor_brasil
# geo_vetor_brasil_anos <- ecodados::geo_vetor_brasil_anos
# geo_vetor_am_sul <- ecodados::geo_vetor_am_sul
# geo_vetor_biomas <- ecodados::geo_vetor_biomas
# geo_vetor_mata_atlantica <- ecodados::geo_vetor_mata_atlantica
# geo_raster_srtm <- ecodados::geo_raster_srtm
# geo_raster_bioclim <- ecodados::geo_raster_bioclim
# geo_raster_globcover_mata_atlantica <- ecodados::geo_raster_globcover_mata_atlantica15.1 Contextualização
Nesta seção, vamos fazer uma breve introdução aos principais conceitos sobre a manipulação e visualização de dados geoespaciais no R. Iremos abordar temas de forma teórica e prática, utilizando a linguagem R, focando em: i) formatos de dados vetoriais e dados raster, ii) Sistemas de Referências de Coordenadas e unidades (geográficas e projetadas), iii) fontes de dados, iv) importar e exportar dados, v) descrição de objetos geoespaciais e vi) principais operações (atributos, espaciais e geométricas). Num segundo momento, criaremos mapas com seus principais elementos como mapas principal e secundário, título, legenda, barra de escala, indicador de orientação (Norte), gride de coordenadas, descrição do Sistema de Referência de Coordenadas e informações de origem dos dados. Por fim, apresentaremos exemplos de aplicações de análises geoespaciais para dados ecológicos, focadas em: i) agregar informações sobre a biodiversidade, ii) preparar dados para compor variáveis preditoras, e iii) como fazer predições espaciais de distribuição de uma espécie e riqueza de espécies.
Esse capítulo segue parte da estrutura organizada por Lovelace et al. (2019), principalmente os Capítulos 2 a 8, sendo adaptado para atender aos principais requisitos que julgamos necessários a estudos ecológicos. Entretanto, não foi possível cobrir todos os assuntos sobre o uso de dados geoespaciais no R, sendo um tema muito extenso que requer a leitura de livros especializados na área como: i) Mas et al. (2019) Análise espacial com R, ii) Wegmann, Leutner & Dech (2016) Remote Sensing and GIS for Ecologists: Using Open Source Software, iii) Wegmann, Schwalb-Willmann & Dech (2020) An Introduction to Spatial Data Analysis Remote Sensing and GIS with Open Source Software, e iv) Fletcher & Fortin (2018) Spatial ecology and conservation modeling: Applications with R. Outros livros sobre a análise geoespacial no R podem ser consultados no Capítulo 11 - Geospatial do Big Book of R.
15.2 Vetor
Dados vetoriais são usados para mapear fenômenos ou objetos espacialmente explícitos que possuem localização ou dimensões bem definidas, representado a partir de formas geométricas (como pontos, linhas e polígonos) e possuem a possibilidade de ter associado a eles informações tabulares. A tabela de atributos é uma tabela que inclui dados geoespaciais e dados alfanuméricos. Os dados geoespaciais são representados por feições geolocalizadas espacialmente (ponto, linha ou polígono), e os dados alfanuméricos (tabela de dados). Dessa forma, a tabela de atributos reúne informações sobre cada feição e pode ser utilizada para realizar filtros ou agregações dos dados de cada feição (Figura 15.1).
![Representação das geometrias de ponto, linha e polígono e atributos. Adaptado de: Olaya [-@olaya2020].](img/cap15_fig01.png)
Figura 15.1: Representação das geometrias de ponto, linha e polígono e atributos. Adaptado de: Olaya (2020).
15.2.1 sf: principal pacote no R para dados vetoriais
Atualmente o principal pacote para trabalhar com dados vetoriais no R é o sf, que implementou o Simple Feature (E. Pebesma 2018). Entretanto, outro pacote pode ser tão versátil quanto o sf, no caso o terra, com algumas mudanças na sintaxe que não abordaremos nesse livro por questões de redução de espaço.
Os tipos de geometrias apresentadas são representados por diferentes classes: POINT, LINESTRING e POLYGON para apenas uma feição de cada tipo de geometria; MULTIPOINT, MULTILINESTRING e MULTIPOLYGON para várias feições de cada tipo de geometria e; GEOMETRYCOLLECTION para várias feições e tipos de geometrias e classes.
Ao olharmos as informações de um objeto da classe sf, podemos notar diversas informações que descrevem o mesmo, numa espécie de cabeçalho:
- resumo do vetor: indica o número de feições (linhas) e campos (colunas)
- tipo da geometria: umas das sete classes (ou mais outras) listadas anteriormente
- dimensão: número de dimensões, geralmente duas (XY)
- bbox (bordas): coordenadas mínimas e máximas da longitude e latitude
-
informação do CRS:
epsgouproj4stringindicando o CRS (Coordinate Reference System) -
tibble: tabela de atributos, com destaque para a coluna
geomougeometryque representa cada feição ou geometria
## Dados vetoriais de polígonos do mundo
data(world)
world
#> Simple feature collection with 177 features and 10 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
#> Geodetic CRS: WGS 84
#> # A tibble: 177 × 11
#> iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp gdpPercap geom
#> * <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <MULTIPOLYGON [°]>
#> 1 FJ Fiji Oceania Oceania Melanesia Sovereign coun… 1.93e4 8.86e5 70.0 8222. (((-180 -16.55522, -179.…
#> 2 TZ Tanzania Africa Africa Eastern Africa Sovereign coun… 9.33e5 5.22e7 64.2 2402. (((33.90371 -0.95, 31.86…
#> 3 EH Western Sahara Africa Africa Northern Africa Indeterminate 9.63e4 NA NA NA (((-8.66559 27.65643, -8…
#> 4 CA Canada North America Americas Northern America Sovereign coun… 1.00e7 3.55e7 82.0 43079. (((-132.71 54.04001, -13…
#> 5 US United States North America Americas Northern America Country 9.51e6 3.19e8 78.8 51922. (((-171.7317 63.78252, -…
#> 6 KZ Kazakhstan Asia Asia Central Asia Sovereign coun… 2.73e6 1.73e7 71.6 23587. (((87.35997 49.21498, 86…
#> 7 UZ Uzbekistan Asia Asia Central Asia Sovereign coun… 4.61e5 3.08e7 71.0 5371. (((55.96819 41.30864, 57…
#> 8 PG Papua New Guinea Oceania Oceania Melanesia Sovereign coun… 4.65e5 7.76e6 65.2 3709. (((141.0002 -2.600151, 1…
#> 9 ID Indonesia Asia Asia South-Eastern Asia Sovereign coun… 1.82e6 2.55e8 68.9 10003. (((104.37 -1.084843, 104…
#> 10 AR Argentina South America Americas South America Sovereign coun… 2.78e6 4.30e7 76.3 18798. (((-68.63401 -52.63637, …
#> # … with 167 more rowsPodemos fazer um mapa simples utilizando a função plot() desse objeto. Para facilitar, escolheremos apenas a primeira coluna [1] (Figura 15.2). Caso não escolhermos apenas uma coluna, um mapa para cada coluna será plotado.
📝 Importante
Faremos mapas mais elaborados na seção de visualização de dados geoespaciais deste capítulo.
Figura 15.2: Mapa vetorial do mundo.
15.3 Raster
Os dados no formato raster consistem em uma matriz (com linhas e colunas) em que os elementos representam células, geralmente igualmente espaçadas (pixels; Figura 15.3). As células dos dados raster possuem duas informações: i) identificação das células (IDs das células) para especificar sua posição na matriz (Figura 15.3 A) e; ii) valores das células (Figura 15.3 B), que geralmente são coloridos para facilitar a interpretação da variação dos valores no espaço (Figura 15.3 C). Além disso, valores ausentes ou não amostrados são representados por NA, ou seja, not available (Figura 15.3 B e C).
![Raster: (A) IDs das células, (B) valores das células, (C) células coloridas. Adaptado de: Lovelace et al. [-@lovelace2019].](livro_files/figure-html/fig-raster-1.png)
Figura 15.3: Raster: (A) IDs das células, (B) valores das células, (C) células coloridas. Adaptado de: Lovelace et al. (2019).
Podemos ainda fazer uma comparação com as representações de dados vetoriais vistos na Figura 15.1, mas agora no formato raster (Figura 15.4).
![Representação das geometrias de ponto, linha e polígono no formato raster. Adaptado de: Olaya [-@olaya2020].](img/cap15_fig02.png)
Figura 15.4: Representação das geometrias de ponto, linha e polígono no formato raster. Adaptado de: Olaya (2020).
15.3.1 raster: principal pacote no R para dados raster
Atualmente, o principal pacote para trabalhar com dados raster é o raster, apesar de existir outros dois: terra e stars, com algumas mudanças na sintaxe que não abordaremos neste livro.
O pacote raster fornece uma ampla gama de funções para criar, importar, exportar, manipular e processar dados raster no R. O objeto raster criado à partir do pacote raster pode assumir três classes: RasterLayer, RasterStack e RasterBrick.
A classe RasterLayer representa apenas uma camada raster. Para criar ou importar um raster no R podemos utilizar a função raster::raster(). Observando essa classe, podemos notar as seguintes informações:
- class: classe raster do objeto raster
- dimensions: número de linhas, colunas e células
- resolution: largura e altura da célula
- extent: coordenadas mínimas e máximas da longitude e latitude
- crs: Sistema de Referência de Coordenadas (CRS)
- source: fonte dos dados (memória ou disco)
- names: nome das camadas
- values: valores máximos e mínimos das células
Vamos utilizar os dados volcano, que possui informações topográficas (elevação) do vulcão Maunga Whau de Auckland na Nova Zelândia.
## Dados de altitude de um vulcão
volcano[1:5, 1:5]
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 100 100 101 101 101
#> [2,] 101 101 102 102 102
#> [3,] 102 102 103 103 103
#> [4,] 103 103 104 104 104
#> [5,] 104 104 105 105 105Vamos transformar essa matriz de dados em um raster com a função raster::raster().
## Rasterlayer
raster_layer <- raster::raster(volcano)
raster_layer
#> class : RasterLayer
#> dimensions : 87, 61, 5307 (nrow, ncol, ncell)
#> resolution : 0.01639344, 0.01149425 (x, y)
#> extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#> crs : NA
#> source : memory
#> names : layer
#> values : 94, 195 (min, max)Um mapa simples do objeto raster pode ser obtido utilizando a função plot(), do próprio pacote raster (Figura 15.5).
Figura 15.5: Mapa simples de um RasterLayer.
Além da classe RasterLayer, há mais duas classes que trabalham com múltiplas camadas: RasterBrick e RasterStack. Elas diferem em relação ao formato dos arquivos suportados, tipo de representação interna e velocidade de processamento.
A classe RasterBrick geralmente corresponde à importação de um único arquivo de imagem de satélite multiespectral (multicamadas) ou a um único objeto com várias camadas na memória. A função raster::brick() cria um objeto RasterBrick.
## Raster layers
raster_layer1 <- raster_layer
raster_layer2 <- raster_layer * raster_layer
raster_layer3 <- sqrt(raster_layer)
raster_layer4 <- log10(raster_layer)
## Raster brick
raster_brick <- raster::brick(raster_layer1, raster_layer2,
raster_layer3, raster_layer4)
raster_brick
#> class : RasterBrick
#> dimensions : 87, 61, 5307, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 0.01639344, 0.01149425 (x, y)
#> extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#> crs : NA
#> source : memory
#> names : layer.1, layer.2, layer.3, layer.4
#> min values : 94.000000, 8836.000000, 9.695360, 1.973128
#> max values : 195.000000, 38025.000000, 13.964240, 2.290035Ao utilizarmos a função plot() do pacote raster, podemos visualizar os raster contidos no objeto RasterBrick (Figura 15.6).
Figura 15.6: Mapas simples de um raster RasterBrick.
Já a classe RasterStack permite conectar vários objetos raster armazenados em arquivos diferentes ou vários objetos no ambiente do R. Um RasterStack é uma lista de objetos RasterLayer com a mesma extensão, resolução e CRS. Uma maneira de criá-lo é com a junção de vários objetos geoespaciais já existentes no ambiente do R ou listar vários arquivos raster em um diretório armazenado no disco. A função raster::stack() cria um objeto RasterStack.
Outra diferença é que o tempo de processamento, para objetos RasterBrick geralmente é menor do que para objetos RasterStack. A decisão sobre qual classe Raster deve ser usada depende principalmente do caráter dos dados de entrada.
## Raster layers
raster_layer1 <- raster_layer
raster_layer2 <- raster_layer * raster_layer
raster_layer3 <- sqrt(raster_layer)
raster_layer4 <- log10(raster_layer)
## Raster stack
raster_stack <- raster::stack(raster_layer1, raster_layer2,
raster_layer3, raster_layer4)
raster_stack
#> class : RasterStack
#> dimensions : 87, 61, 5307, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 0.01639344, 0.01149425 (x, y)
#> extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#> crs : NA
#> names : layer.1, layer.2, layer.3, layer.4
#> min values : 94.000000, 8836.000000, 9.695360, 1.973128
#> max values : 195.000000, 38025.000000, 13.964240, 2.290035Da mesma forma, ao utilizar a função plot() do pacote raster, podemos visualizar os raster contidos no objeto RasterStack (Figura 15.7).
Figura 15.7: Mapas simples de um raster RasterStack.
15.4 Sistema de Referência de Coordenadas e Unidades
Os dados geoespaciais (vetor e raster) possuem ainda um outro componente fundamental que é o Sistema de Referência de Coordenadas, ou do inglês Coordinate Reference System (CRS). Esse componente define a referência espacial dos elementos geoespaciais (vetor e raster) na superfície da Terra. Ele é composto por dois principais conceitos: primeiro, que tipos de unidades estão sendo utilizadas para a representação geográfica, podendo assumir dois tipos - ângulos ou metros, que definem o Sistema de Coordenadas Geográficas e o Sistema de Coordenadas Projetadas, respectivamente. O segundo componente é o datum, que é a relação do sistema de coordenadas (geográfica ou projetada) com a superfície da Terra. Esse último componente faz parte de uma área da Cartografia denominada Geodésia que estuda a forma e dimensões da Terra, campo gravitacional e a localização de pontos fixos e sistemas de coordenadas. O livro de Lapaine & Usery (2017) é um excelente material para se aprofundar nesse assunto.
15.4.1 Sistema de Coordenadas Geográficas
O Sistema de Coordenadas Geográficas utiliza ângulos (graus) para representar feições na superfície da Terra através de dois valores: longitude e latitude. A longitude representa o eixo Leste-Oeste e a latitude o eixo Norte-Sul. Nesse sistema, a superfície da Terra é representada geralmente por uma superfície elipsoidal, pois a Terra é ligeiramente achatada nos polos devido ao movimento de rotação.
15.4.2 Sistema de Coordenadas Projetadas
O Sistema de Coordenadas Projetadas utiliza um Sistema Cartesiano de Coordenadas em uma superfície plana. Dessa forma, a partir de uma origem traçam-se eixos X e Y e uma unidade linear é utilizada, como o metro. Todos as projeções feitas de sistemas geoespaciais convertem uma superfície tridimensional em uma superfície plana bidimensional. Sendo assim, essa conversão traz consigo algum tipo de distorção em relação à porção real, podendo ser distorções em: i) formas locais, ii) áreas, iii) distâncias, iv) flexão ou curvatura, v) assimetria ou vi) lacunas de continuidade. Dessa forma, um sistema de coordenadas projetadas pode preservar somente uma ou duas dessas propriedades.
Existem três grandes grupos de projeções: i) cilíndricos, ii) cônicos e iii) planares. Na projeção cilíndrica, a superfície da Terra é mapeada em um cilindro, criada tocando a superfície da Terra ao longo de uma ou duas linhas de tangência, sendo utilizada com mais frequência para mapear todo o globo tendo como exemplo mais conhecido a Projeção Universal Transversa de Mercator (UTM). Na projeção cônica, a superfície da Terra é projetada em um cone ao longo de uma linha ou duas linhas de tangência, de modo que as distorções são minimizadas ao longo das linhas e aumentam com a distância das mesmas sendo, portanto, mais adequada para mapear áreas de latitudes médias, tendo como exemplo mais conhecido a Projeção Cônica Equivalente de Albers e a Projeção Cônica Conforme de Lambert. E na projeção plana, também denominada Projeção Azimutal, o mapeamento toca o globo em um ponto ou ao longo de uma linha de tangência, sendo normalmente utilizado no mapeamento de regiões polares, sendo a mais comum a Projeção Azimutal Equidistante, a mesma utilizada na bandeira da ONU.
15.4.3 Datum
Como dito anteriormente, o datum é a relação do sistema de coordenadas com a superfície da Terra. Ele representa o ponto de intersecção do elipsoide de referência com a superfície da Terra (geoide, a forma verdadeira da Terra), compensando as diferenças do campo gravitacional da Terra. Existem dois tipos de datum: i) local e ii) geocêntrico. Em um datum local, como o SAD69 - South American Datum 1969, o elipsoide de referência é deslocado para se alinhar com a superfície em um determinado local, por exemplo, na América do Sul. Já em um datum geocêntrico, como WGS84 - World Geodetic System 1984, o centro do elipsoide é o centro de gravidade da Terra e a precisão das projeções não é otimizada para um local específico do globo.
No Brasil, desde 2015, o Instituto Brasileiro de Geografia e Estatística (IBGE) ajudou a desenvolver e reafirmou o uso do datum SIRGAS2000 - Sistema de Referencia Geocéntrico para las Américas 2000 para todos os mapeamentos realizados no Brasil, um esforço conjunto para adotar o mesmo datum em toda a América. Mais sobre esse datum pode ser lido aqui: SIRGAS2000.
15.4.4 Sistema de Referência de Coordenadas (CRS) no R
No R, há duas formas principais de representar um Sistema de Referência de Coordenadas: i) código epsg e ii) proj4string. O código EPSG (European Petroleum Survey Group) é uma sequência de números curta, referindo-se apenas a um CRS. O site epsg.io permite consultar diversas informações como procurar por um código, representação de mapas e fazer transformações de CRS.
Já o proj4string permite mais flexibilidade para especificar diferentes parâmetros, como o tipo de projeção, datum e elipsoide. Dessa forma, é possível especificar muitas projeções, ou mesmo modificar as projeções existentes, tornando a representação proj4string mais complexa e flexível.
Além disso, ainda é possível consultar uma extensa lista de CRSs no site spatialreference.org, que fornece descrições em diversos formatos, baseados em GDAL e Proj.4. Essa abordagem permite consultar uma URL que pode produzir uma referência espacial em um formato que seu software SIG ou o R pode utilizar como referência.
Os pacotes (geo)espaciais no R suportam uma ampla variedade de CRSs e usam a biblioteca PROJ. A função rgdal::make_EPSG() retorna um data frame das projeções disponíveis, com informações dos códigos epsg e proj4string numa mesma tabela, facilitando a busca e uso de CRSs (Tabela 15.1).
## Listagem dos Sistemas de Referências de Coordenadas no R
crs_data <- rgdal::make_EPSG()
head(crs_data)| code | note | prj4 | prj_method |
|---|---|---|---|
| 2000 | Anguilla 1957 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
| 2001 | Antigua 1943 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
| 2002 | Dominica 1945 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
| 2003 | Grenada 1953 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
| 2004 | Montserrat 1958 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
| 2005 | St. Kitts 1955 / British West Indies Grid | +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs +type=crs | Transverse Mercator |
15.5 Principais fontes de dados geoespaciais
Existem diversas fontes de dados geoespaciais em diferentes bases de dados disponíveis gratuitamente. Geralmente essas bases de dados são disponibilizadas separadamente em apenas dados vetoriais e dados raster. Para dados vetoriais, grande parte dos dados disponibilizados são utilizados em mapas como limites políticos, limites de biomas ou distribuição de espécies para polígonos; estradas e rios para dados lineares, ou ainda pontos de ocorrência de espécies ou comunidades, ou medidas tomadas em campo sobre condições naturais como clima ou relevo, como pontos. Entretanto, é sempre recomendado o uso de bases oficiais, principalmente em relação a dados vetoriais de limites políticos. Para tanto, é fundamental buscar as bases oficiais de cada país, entretanto, há bases que podem ser utilizadas globalmente, como veremos.
Sobre as bases de dados raster, há uma infinidade de dados para diferentes objetivos, mas grande parte deles são relativos a condições ambientais, representando uma variável de interesse de forma contínua no espaço, como temperatura, precipitação, elevação, etc.
Há uma compilação de dados geoespaciais vetoriais e raster feita por Marcus Vinícius Alves de Carvalho e Angelica Carvalho Di Maio, chamada GeoLISTA. Entretanto, como as bases de dados tendem a ser muito dinâmicas, é possível que muitas bases tenham surgido e desaparecido desde a listagem realizada.
Além das bases de dados, há pacotes específicos no R que fazem o download de dados vetoriais e rasters, facilitando a aquisição e reprodutibilidade. Para conferir uma listagem completa de pacotes para diversas análises espaciais, veja CRAN Task View: Analysis of Spatial Data.
15.5.1 Vetor
Dentre as bases vetoriais, destacamos as seguintes na Tabela 15.2.
| Bases de dados | Descrição |
|---|---|
| IBGE | Limites territoriais e censitários do Brasil |
| FBDS | Uso da terra, APP e hidrografia - Mata Atlântica e Cerrado |
| GeoBank | Dados geológicos do Brasil |
| Pastagem.org | Dados de pastagens e gado para o Brasil |
| CanaSat | Dados de cana-de-açúcar para o Brasil |
| CSR Maps | Diversos dados vetoriais e raster para o Brasil |
| Ecoregions | Dados de biorregiões e biomas do mundo |
| UN Biodiversity Lab | Diversas bases de dados para o mundo |
| Biodiversity Hotspots | Dados dos limites dos Hotspots de Biodiversidade |
| IUCN Red List of Threatened Species | Dados dos limites das distribuições das espécies para o mundo |
| Map of Life (MOL) | Dados da distribuição de espécies e outros dados para o mundo |
| Key Biodiversity Areas | Dados dos limites das Key Biodiversity Areas |
| HydroSHEDS | Informações hidrológicas do mundo |
| Global Roads Inventory Project (GRIP) | Dados de estradas do mundo todo |
| Database of Global Administrative Areas (GADM) | Limites de áreas administrativas do mundo |
| Natural Earth | Diversos limites para o mundo |
| Protected Planet | Limites de áreas protegidas para o mundo |
| Global Biological Information Facility (GBIF) | Dados de ocorrências de espécies para o mundo |
| Species Link | Dados de ocorrências de espécies para o Brasil |
| Global Invasive Species Information Network (GISIN) | Dados de ocorrências de espécies invasoras para o Mundo |
15.5.2 Raster
Dentre as bases raster, destacamos as seguintes na Tabela 15.3.
| Bases de dados | Descrição |
|---|---|
| MapBiomas | Uso e cobertura da terra para o Brasil, Panamazonia Legal, Chaco e Mata Atlântica de 1985 a 2020 |
| Bahlu | Distribuições históricas de terras agrícolas e pastagens para todo o Brasil de 1940 a 2012 |
| USGS | Dados de diversos satélites livres para o mundo |
| SRTM | Dados de elevação para o mundo |
| Geoservice Maps | Dados de elevação e florestas para o mundo |
| Global Forest Watch | Dados de florestas para o mundo |
| GlobCover | Dados de uso e cobertura da terra para todo o planeta |
| Landcover | Dados de uso e cobertura da terra para todo o planeta |
| Global Human Footprint | Dados de pegada ecológica para o mundo |
| GHSL - Global Human Settlement Layer | Dados e ferramentas abertos e gratuitos para avaliar a presença humana no planeta |
| Land-Use Harmonization (LUH2) | Dados atuais e previsões de uso da terra |
| ESA Climate Change Initiative | Arquivos globais de observação da Terra nos últimos 30 anos da Agência Espacial Europeia (ESA) |
| WorldClim | Dados climáticos para o mundo |
| CHELSA | Dados climáticos para o mundo |
| EarthEnv | Dados de cobertura da terra, nuvens, relevo e hidrografia |
| SoilGrids | Dados de solo para o mundo |
| Global Wetlands | Dados de áreas úmidas para o mundo |
| Global Surface Water Explorer | Dados de águas superficiais para o mundo |
| MARSPEC | Dados de condições do oceano para o mundo |
| Bio-ORACLE | Dados de condições do oceano para o mundo |
15.5.3 Pacotes do R
Dentre os pacotes no R para download de dados geoespaciais, destacamos os seguintes na Tabela 15.4.
| Pacotes | Descrição |
|---|---|
| geobr | Carrega Shapefiles de Conjuntos de Dados Espaciais Oficiais do Brasil |
| rnaturalearth | Dados do mapa mundial da Natural Earth |
| geodata | Diversas bases de dados para o mundo |
| rworldmap | Visualização de dados globais |
| spData | Conjuntos de dados para análise espacial |
| OpenStreetMap | Acesso para abrir imagens raster de mapas de ruas |
| osmdata | Baixa e importa dados do OpenStreetMap |
| geonames | Interface para o serviço da Web de consulta espacial ‘Geonames’ |
| rgbif | Interface para o Global ‘Biodiversity’ Information Facility API |
| maptools | Ferramentas para lidar com objetos geoespaciais |
| marmap | Importar, traçar e analisar dados batimétricos e topográficos |
| oce | Fonte e processamento de dados oceanográficos |
| envirem | Geração de variáveis ENVIREM |
| sdmpredictors | Conjuntos de dados preditor de modelagem de distribuição de espécies |
| metScanR | Encontra, mapeia e coleta dados e metadados ambientais |
| ClimDown | Biblioteca de redução de escala do clima para a produção diária do modelo climático |
| rWBclimate | Acessa dados climáticos do Banco Mundial |
| rnoaa | Dados meteorológicos ‘NOAA’ de R |
| RNCEP | Obtém, organiza e visualiza dados meteorológicos NCEP |
| smapr | Aquisição e processamento de dados ativos-passivos (SMAP) de umidade do solo da NASA |
15.6 Importar e exportar dados geoespaciais
Agora que sabemos o que são dados geoespaciais e em quais bases de dados podemos buscar e baixar esses dados, veremos seus principais formatos e como importá-los e exportá-los do R.
15.6.1 Principais formatos de arquivos geoespaciais
Há diversos formatos de arquivos geoespaciais, alguns específicos para dados vetoriais e raster, e outros no formato de banco de dados geoespaciais, como PostGIS, que podem armazenar ambos os formatos.
Entretanto, todos os formatos para serem importados para o R usam o GDAL (Geospatial Data Abstraction Library), uma interface unificada para leitura e escrita de diversos formatos de arquivos geoespaciais, sendo utilizado também por uma série de softwares de GIS como QGIS, GRASS GIS e ArcGIS.
Dentre esses formatos, destacamos os seguintes na Tabela 15.5.
| Nome | extensão | Descrição | Tipo | Modelo |
|---|---|---|---|---|
| ESRI Shapefile | .shp (arquivo principal) | Formato popular que consiste em pelo menos quatro arquivos: .shp (feição), .dbf (tabela de atributos), .shx (ligação entre .shp e .dbf) e .prj (projeção) | Vetor | Parcialmente aberto |
| GeoJSON | .geojson | Estende o formato de troca JSON incluindo um subconjunto da representação de recurso simples | Vetor | Aberto |
| KML | .kml | Formato baseado em XML para visualização espacial, desenvolvido para uso com o Google Earth. O arquivo KML compactado forma o formato KMZ | Vetor | Aberto |
| GPX | .gpx | Esquema XML criado para troca de dados de GPS | Vetor | Aberto |
| GeoTIFF | .tif/.tiff | Formato raster popular. Um arquivo TIFF contendo metadados espaciais adicionais. | Raster | Aberto |
| Arc ASCII | .asc | Formato de texto em que as primeiras seis linhas representam o cabeçalho raster, seguido pelos valores das células raster organizadas em linhas e colunas | Raster | Aberto |
| NetCDF | .nc | NetCDF (Network Common Data Form) é um conjunto de bibliotecas de software e formatos de dados independentes que suportam a criação, acesso e compartilhamento de dados científicos orientados a arrays | Raster | Aberto |
| BIL | .bil/.hdr | BIL (Banda intercalada por linha) são métodos comuns de organização para imagens multibanda, geralmente acompanhados por um arquivo .hdr, descrevendo atributos específicos da imagem | Raster | Aberto |
| R-raster | .gri/ .grd | Formato raster nativo do raster do pacote R | Raster | Aberto |
| SQLite/SpatiaLite | .sqlite | Banco de dados relacional autônomo | Vetor e raster | Aberto |
| ESRI FileGDB | .gdb | objetos geoespaciais e não espaciais criados pelo ArcGIS. Permite: várias classes de recursos; topologia | Vetor e raster | Proprietário |
| GeoPackage | .gpkg | Contêiner de banco de dados leve baseado em SQLite permitindo uma troca fácil e independente de plataforma de geodados | Vetor e raster | Aberto |
O formato mais comum para arquivos vetoriais é o ESRI Shapefile; para arquivos raster é o GeoTIFF; e para dados climáticos em múltiplas camadas, geralmente há a disponibilização de dados no formato NetCDF. Entretanto, recentemente tivemos o surgimento do GeoPackage, que possui diversas vantagens em relação aos formatos anteriores, podendo armazenar em apenas um arquivo, dados no formato vetorial, raster e também dados não-espaciais (e.g., tabelas), além de possuir uma grande integração com diversos softwares e bancos de dados.
15.6.2 Importar dados
As principais funções para importar dados no R são: i) para vetores a função sf::st_read(), e ii) para raster a função raster::raster() e suas variações raster::brick() e raster::stack() para múltiplas camadas. Essas funções atribuem objetos ao seu espaço de trabalho, armazenando-os na memória RAM disponível em seu hardware, sendo essa a maior limitação para trabalhar com dados geoespaciais no R. Por exemplo, se um arquivo raster possui mais de 8 Gb de tamanho, e seu computador possui exatamente 8 Gb de RAM, é muito provável que ele não seja importado ou mesmo criado como um objeto dentro do ambiente R. Existem soluções para esses problemas, mas não as abordaremos neste capítulo.
Vetor
Como vimos, os arquivos vetoriais são disponibilizados em diversos formatos. Para sabermos se um determinado formato pode ser importado ou exportado utilizando o pacote sf, podemos utilizar a função sf::st_drivers(). Uma amostra desses formatos é apresentado na Tabela 15.6.
## Formatos vetoriais importados e exportados pelo pacote sf
head(sf::st_drivers())| name | long_name | write | copy | is_raster | is_vector | vsi |
|---|---|---|---|---|---|---|
| ESRIC | Esri Compact Cache | FALSE | FALSE | TRUE | TRUE | TRUE |
| FITS | Flexible Image Transport System | TRUE | FALSE | TRUE | TRUE | FALSE |
| PCIDSK | PCIDSK Database File | TRUE | FALSE | TRUE | TRUE | TRUE |
| netCDF | Network Common Data Format | TRUE | TRUE | TRUE | TRUE | TRUE |
| PDS4 | NASA Planetary Data System 4 | TRUE | TRUE | TRUE | TRUE | TRUE |
| VICAR | MIPL VICAR file | TRUE | TRUE | TRUE | TRUE | TRUE |
Importar dados vetoriais existentes
Para importar vetores existentes para o R, utilizaremos a função sf::st_read(). A estrutura é semelhante para todos os formatos descritos na Tabela 15.6, de modo que sempre preencheremos o argumento dsn (data source name) com o nome do arquivo a ser importado. Entretanto, para banco de dados, como GeoPackage, pode ser necessário especificar a camada que se tem interesse com um segundo argumento chamado layer, com o nome da camada.
Para quase todas as operações vetoriais nesse capítulo, usaremos os dados disponíveis para o município de Rio Claro/SP. Primeiramente, baixaremos esses dados da FBDS (Fundação Brasileira para o Desenvolvimento Sustentável), através desse repositório de dados. Em 2013, a FBDS deu início ao Projeto de Mapeamento em Alta Resolução dos Biomas Brasileiros, mapeando a cobertura da terra, hidrografia (nascentes, rios e lagos) e Áreas de Preservação Permanente (APPs). O mapeamento foi concluído para os municípios dos Biomas Mata Atlântica e Cerrado, e mais recentemente para os outros biomas. Para fazer o download dos arquivos de interesse, utilizaremos o R, através da função download.file().
Primeiramente, criaremos um diretório com a função dir.create(), usando a função here::here() para indicar o repositório (ver o Capítulo 5).
## Criar diretório
dir.create(here::here("dados"))
dir.create(here::here("dados", "vetor"))Em seguida, vamos fazer o download de pontos de nascentes, linhas de hidrografia e polígonos de cobertura da terra para o município de Rio Claro/SP.
## Aumentar o tempo de download
options(timeout = 1e3)
## Download
for(i in c(".dbf", ".prj", ".shp", ".shx")){
# Pontos de nascentes
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/HIDROGRAFIA/SP_3543907_NASCENTES", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_NASCENTES", i)), mode = "wb")
# Linhas de hidrografia
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/HIDROGRAFIA/SP_3543907_RIOS_SIMPLES", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_RIOS_SIMPLES", i)), mode = "wb")
# Polígonos de cobertura da terra
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/USO/SP_3543907_USO", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_USO", i)), mode = "wb")
}Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_nascentes
ecodados::geo_vetor_hidrografia
ecodados::geo_vetor_coberturaAgora podemos importar esses dados para o R. Primeiro vamos importar as nascentes (Figura 15.8).
## Importar nascentes
geo_vetor_nascentes <- sf::st_read(
here::here("dados", "vetor", "SP_3543907_NASCENTES.shp"), quiet = TRUE)
## Plot
plot(geo_vetor_nascentes[1], pch = 20, col = "blue", main = NA,
axes = TRUE, graticule = TRUE)
Figura 15.8: Mapa de nascentes de Rio Claro/SP.
Agora vamos importar a hidrografia (Figura 15.9).
## Importar hidrografia
geo_vetor_hidrografia <- sf::st_read(
here::here("dados", "vetor", "SP_3543907_RIOS_SIMPLES.shp"), quiet = TRUE)
## Plot
plot(geo_vetor_hidrografia[1], col = "steelblue", main = NA, axes = TRUE, graticule = TRUE)
Figura 15.9: Mapa da hidrografia de Rio Claro/SP.
E por fim, vamos importar a cobertura da terra (Figura 15.10).
## Importar cobertura da terra
geo_vetor_cobertura <- sf::st_read(
here::here("dados", "vetor", "SP_3543907_USO.shp"), quiet = TRUE)
## Plot
plot(geo_vetor_cobertura[5],
col = c("blue", "orange", "gray30", "forestgreen", "green"),
main = NA, axes = TRUE, graticule = TRUE)
legend(x = .1, y = .3, pch = 15, cex = .7, pt.cex = 2.5,
legend = (geo_vetor_cobertura$CLASSE_USO),
col = c("blue", "orange", "gray30", "forestgreen", "green"))
Figura 15.10: Mapa de cobertura da terra de Rio Claro/SP.
Importar utilizando pacotes
Além de dados existentes, podemos importar dados vetoriais de pacotes, como listado anteriormente na Tabela 15.4. Para o Brasil, o pacote mais interessante trata-se do geobr, do Instituto de Pesquisa Econômica Aplicada (IPEA), que possui dados oficiais do Instituto Brasileiro de Geografia e Estatística (IBGE).
É possível listar todos os dados disponíveis no pacote através da função geobr::list_geobr(). Na Tabela 15.7 é possível ver alguns desses dados.
## Listar todos os dados do geobr
geobr::list_geobr()| function | geography | years | source |
|---|---|---|---|
read_country
|
Country | 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 | IBGE |
read_region
|
Region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 | IBGE |
read_state
|
States | 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 | IBGE |
read_meso_region
|
Meso region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 | IBGE |
read_micro_region
|
Micro region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 | IBGE |
read_intermediate_region
|
Intermediate region | 2017, 2019, 2020 | IBGE |
Como exemplo, vamos fazer o download o limite do município de Rio Claro/SP, utilizando o código do município (3543907) (Figura 15.11).
📝 Importante
Para saber todos os códigos dos municípios do Brasil, recomendamos a verificação no site do IBGE.
## Polígono do limite do município de Rio Claro
geo_vetor_rio_claro <- geobr::read_municipality(code_muni = 3543907,
year = 2020, showProgress = FALSE)Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_rio_claro
## Plot
plot(geo_vetor_rio_claro[1], col = "gray", main = NA, axes = TRUE, graticule = TRUE)
Figura 15.11: Limite do município de Rio Claro/SP.
Já para o mundo, o pacote mais interessante trata-se do rnaturalearth, que faz o download de dados do Natural Earth. Vamos fazer o download do limite do Brasil (Figura 15.12).
## Polígono do limite do Brasil
geo_vetor_brasil <- rnaturalearth::ne_countries(scale = "large",
country = "Brazil", returnclass = "sf")Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_brasil
## Plot
plot(geo_vetor_brasil[1], col = "gray", main = NA, axes = TRUE, graticule = TRUE)
Figura 15.12: Limite do Brasil.
Criar um objeto espacial de uma tabela de coordenadas
É muito comum em coletas de campo ou bases de dados, ter coordenadas de locais de estudo ou de ocorrências de espécies organizadas em tabelas. Essas tabelas devem possuir duas colunas: longitude e latitude, ou X e Y para dados UTM, por exemplo. Ao importá-las para o R, o formato que assumem pode ser de uma das classes: matrix, data frame ou tibble, ou seja, ainda não são da classe vetorial sf. Nesta seção iremos ver como fazer essa conversão.
Para tanto, vamos usar os dados de comunidades de anfíbios da Mata Atlântica (Atlantic Amphibians, Vancine et al. (2018)). Faremos o download diretamente do site da fonte dos dados. Antes vamos criar um diretório.
## Criar diretório
dir.create(here::here("dados", "tabelas"))Em seguida, vamos fazer o download de um arquivo .zip e vamos extrair usando a função unzip() nesse mesmo diretório.
## Download
download.file(url = "https://esajournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2Fecy.2392&file=ecy2392-sup-0001-DataS1.zip",
destfile = here::here("dados", "tabelas", "atlantic_amphibians.zip"), mode = "wb")
## Unzip
unzip(zipfile = here::here("dados", "tabelas", "atlantic_amphibians.zip"),
exdir = here::here("dados", "tabelas"))Agora podemos importar a tabela de dados com a função readr::read_csv().
## Importar tabela de locais
geo_anfibios_locais <- readr::read_csv(
here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_sites.csv"),
locale = readr::locale(encoding = "latin1")
)
geo_anfibios_locais
#> # A tibble: 1,163 × 25
#> id reference_number species_number record sampled_habitat active_methods passive_methods complementary_meth… period month_start year_start
#> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 amp1001 1001 19 ab fo,ll as pt <NA> mo,da… 9 2000
#> 2 amp1002 1002 16 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 3 amp1003 1002 14 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 4 amp1004 1002 13 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 5 amp1005 1003 30 co fo,ll,br as <NA> <NA> mo,da… 7 1988
#> 6 amp1006 1004 42 co tp,pp,la,ll,is <NA> <NA> <NA> <NA> NA NA
#> 7 amp1007 1005 23 co sp as <NA> <NA> <NA> 4 2007
#> 8 amp1008 1005 19 co sp,la,sw as,sb,tr <NA> <NA> tw,ni 4 2007
#> 9 amp1009 1005 13 ab fo <NA> pt <NA> mo,da… 4 2007
#> 10 amp1010 1006 1 ab fo <NA> pt <NA> mo,da… 5 2011
#> # … with 1,153 more rows, and 14 more variables: month_finish <dbl>, year_finish <dbl>, effort_months <dbl>, country <chr>, state <chr>,
#> # state_abbreviation <chr>, municipality <chr>, site <chr>, latitude <dbl>, longitude <dbl>, coordinate_precision <chr>, altitude <dbl>,
#> # temperature <dbl>, precipitation <dbl>Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_anfibios_locaisPor fim, podemos facilmente criar um objeto espacial do tipo MULTIPOINT utilizando a função sf::st_as_sf(). Podemos ver essas coordenadas plotadas no mapa simples da Figura 15.13 (Vancine et al. 2018).
É necessário antes se ater ao argumento coords que deve indicar as colunas de longitude e latitude, nessa ordem; e também ao argumento crs para indicar o CRS correspondente dessas coordenadas, que aqui sabemos se tratar de coordenadas geográficas e datum WGS84. Então podemos facilmente utilizar o código EPSG 4326. Entretanto, se as coordenadas estiverem em metros, por exemplo, teremos de nos ater a qual CRS as mesmas foram coletadas, ou seja, se forem coordenadas de GPS, é preciso saber como o GPS estava configurado (projeção e datum).
## Converter dados tabulares para sf
geo_anfibios_locais_vetor <- geo_anfibios_locais %>%
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
geo_anfibios_locais_vetor
#> Simple feature collection with 1163 features and 23 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -56.74194 ymin: -33.51083 xmax: -34.79667 ymax: -3.51525
#> Geodetic CRS: WGS 84
#> # A tibble: 1,163 × 24
#> id reference_number species_number record sampled_habitat active_methods passive_methods complementary_meth… period month_start year_start
#> * <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 amp1001 1001 19 ab fo,ll as pt <NA> mo,da… 9 2000
#> 2 amp1002 1002 16 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 3 amp1003 1002 14 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 4 amp1004 1002 13 co fo,la,ll as pt <NA> mo,da… 12 2007
#> 5 amp1005 1003 30 co fo,ll,br as <NA> <NA> mo,da… 7 1988
#> 6 amp1006 1004 42 co tp,pp,la,ll,is <NA> <NA> <NA> <NA> NA NA
#> 7 amp1007 1005 23 co sp as <NA> <NA> <NA> 4 2007
#> 8 amp1008 1005 19 co sp,la,sw as,sb,tr <NA> <NA> tw,ni 4 2007
#> 9 amp1009 1005 13 ab fo <NA> pt <NA> mo,da… 4 2007
#> 10 amp1010 1006 1 ab fo <NA> pt <NA> mo,da… 5 2011
#> # … with 1,153 more rows, and 13 more variables: month_finish <dbl>, year_finish <dbl>, effort_months <dbl>, country <chr>, state <chr>,
#> # state_abbreviation <chr>, municipality <chr>, site <chr>, coordinate_precision <chr>, altitude <dbl>, temperature <dbl>, precipitation <dbl>,
#> # geometry <POINT [°]>
## Plot
plot(geo_anfibios_locais_vetor[1], pch = 20, col = "black",
main = NA, axes = TRUE, graticule = TRUE)
Figura 15.13: Coordenadas das comunidades do Atlantic Amphibians.
Converter dados espaciais sp para sf
O pacote sf é mais recente e mais fácil de manipular objetos vetoriais no R. Seu predecessor, o pacote sp possui uma classe própria e homônima. Entretanto, muitos pacotes de análises espaciais ainda utilizam essa classe em suas funções, apesar dessa migração ter ocorrido rapidamente recentemente. Dessa forma, a conversão entre essas classes pode ser necessária em alguns momentos.
Abaixo, veremos como podemos fazer essa conversão facilmente. Primeiramente, vamos importar dados sp.
## Polígonos países sp
co110_sp <- rnaturalearth::countries110
class(co110_sp)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"Agora, podemos converter facilmente com a função sf::st_as_sf().
Podemos facilmente converter esse objeto novamente para a classe sp com a função sf::as_Spatial().
## Polígonos países sp
co110_sp <- sf::as_Spatial(co110_sf)
class(co110_sp)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"Raster
Para importar dados raster no R, utilizaremos a função raster::raster(), raster::brick() ou raster::stack(). Para apenas uma camada raster, usaremos a função raster::raster(), com o argumento x sendo o nome do arquivo. Já para mais camadas, usaremos raster::brick() para um arquivo que possua múltiplas camadas, ou ainda a função raster::stack() para vários arquivos em diferentes camadas também no argumento x, sendo necessário listar os arquivos no diretório, geralmente utilizando a função dir() ou list.files(). Entretanto, para especificar uma camada, podemos utilizar o argumento band ou layer e o nome dessa camada.
Raster Layer
Primeiramente, vamos criar um diretório para os dados raster que fazeremos o download.
## Criar diretório
dir.create(here::here("dados", "raster"))Em seguida, vamos fazer o download de dados de elevação, na verdade dados de Modelo Digital de Elevação (Digital Elevation Model - DEM), localizados também para o município de Rio Claro. Utilizaremos os dados do Shuttle Radar Topography Mission - SRTM. Para saber mais sobre esses dados, recomendamos a leitura do artigo de Farr et al. (2007).
## Aumentar o tempo de download
options(timeout = 1e3)
## Download
download.file(url = "https://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_27_17.zip",
destfile = here::here("dados", "raster", "srtm_27_17.zip"), mode = "wb")
## Unzip
unzip(zipfile = here::here("dados", "raster", "srtm_27_17.zip"),
exdir = here::here("dados", "raster"))Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_raster_srtmAgora podemos importar essa camada para o R, e visualizá-la em relação ao limite do município de Rio Claro/SP (Figura 15.14).
## Importar raster de altitude
geo_raster_srtm <- raster::raster(here::here("dados", "raster", "srtm_27_17.tif"))
geo_raster_srtm
#> class : RasterLayer
#> dimensions : 6000, 6000, 3.6e+07 (nrow, ncol, ncell)
#> resolution : 0.0008333333, 0.0008333333 (x, y)
#> extent : -50, -45, -25, -20 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : srtm_27_17.tif
#> names : srtm_27_17
#> values : -32768, 32767 (min, max)
## Plot
plot(geo_raster_srtm, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.14: Camada raster do DEM em relação ao limite do município de Rio Claro/SP.
Raster Stack
Além dos dados de elevação, dados de temperatura e precipitação podem ser obtidos do WorldClim. Para saber mais sobre esses dados, recomendamos a leitura do artigo Fick & Hijmans (2017).
## Aumentar o tempo de download
options(timeout = 1e3)
## Download
download.file(url = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_bio.zip",
destfile = here::here("dados", "raster", "wc2.0_10m_bio.zip"), mode = "wb")
## Unzip
unzip(zipfile = here::here("dados", "raster", "wc2.0_10m_bio.zip"),
exdir = here::here("dados", "raster"))Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_raster_bioclimPara importar essa série de camadas, primeiramente listaremos os arquivos e depois importaremos no formato RasterStack (Figura 15.15).
## Listar arquivos
arquivos_raster <- dir(path = here::here("dados", "raster"), pattern = "wc") %>%
grep(".tif", ., value = TRUE)
arquivos_raster
#> [1] "wc2.1_10m_bio_1.tif" "wc2.1_10m_bio_10.tif" "wc2.1_10m_bio_11.tif" "wc2.1_10m_bio_12.tif" "wc2.1_10m_bio_13.tif" "wc2.1_10m_bio_14.tif"
#> [7] "wc2.1_10m_bio_15.tif" "wc2.1_10m_bio_16.tif" "wc2.1_10m_bio_17.tif" "wc2.1_10m_bio_18.tif" "wc2.1_10m_bio_19.tif" "wc2.1_10m_bio_2.tif"
#> [13] "wc2.1_10m_bio_3.tif" "wc2.1_10m_bio_4.tif" "wc2.1_10m_bio_5.tif" "wc2.1_10m_bio_6.tif" "wc2.1_10m_bio_7.tif" "wc2.1_10m_bio_8.tif"
#> [19] "wc2.1_10m_bio_9.tif"
## Importar vários rasters como stack
geo_raster_bioclim <- raster::stack(here::here("dados", "raster", arquivos_raster))
geo_raster_bioclim
#> class : RasterStack
#> dimensions : 1080, 2160, 2332800, 19 (nrow, ncol, ncell, nlayers)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> names : wc2.1_10m_bio_1, wc2.1_10m_bio_10, wc2.1_10m_bio_11, wc2.1_10m_bio_12, wc2.1_10m_bio_13, wc2.1_10m_bio_14, wc2.1_10m_bio_15, wc2.1_10m_bio_16, wc2.1_10m_bio_17, wc2.1_10m_bio_18, wc2.1_10m_bio_19, wc2.1_10m_bio_2, wc2.1_10m_bio_3, wc2.1_10m_bio_4, wc2.1_10m_bio_5, ...
#> min values : -54.724354, -37.781418, -66.311249, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 9.131122, 0.000000, -29.686001, ...
#> max values : 30.98764, 38.21617, 29.15299, 11191.00000, 2381.00000, 484.00000, 229.00169, 5284.00000, 1507.00000, 5282.00000, 4467.00000, 21.14754, 100.00000, 2363.84595, 48.08275, ...
Figura 15.15: Camadas rasters do WorldClim (BIO01 e BIO12) para o mundo.
15.6.3 Exportar dados
Saber a melhor forma de exportar dados geoespaciais de objetos recém-criados no R é fundamental, principalmente porque essa ação dependerá do tipo de dado (vetor ou raster), classe do objeto (por exemplo, MULTIPOINT ou RasterLayer) e tipo e quantidade de informações armazenadas (por exemplo, tamanho do objeto, intervalo de valores, etc.).
Vetor
Para dados vetoriais, a principal função utilizada é a sf::st_write(). Essa função permite gravar objetos sf em vários formatos de arquivos vetoriais, como .shp, .gpkg ou .geojson. O formato a ser exportado vai influenciar na velocidade do processo de gravação.
Os argumentos dessa função será o obj que é o objeto sf criado no ambiente R, e o dsn (data source name), ou seja, o nome que o arquivo terá ao ser exportado do R, de modo que o complemento .shp no nome de saída, por exemplo, definirá que o arquivo terá a extensão ESRI Shapefile. Entretanto, essa extensão pode ser definida também utilizando o argumento driver, com as possibilidades listadas nesse site.
## Exportar o polígono de Rio Claro na extensão ESRI Shapefile
sf::st_write(obj = geo_vetor_rio_claro,
dsn = here::here("dados", "vetor", "geo_vetor_rio_claro.shp"))Ou podemos ainda exportar o objeto vetorial na extensão GeoPackage. Entretanto, aqui é interessante acrescentar um argumento chamado layer para definir o nome das camadas a serem exportadas no mesmo arquivo GeoPackage, por exemplo.
## Exportar o polígono de Rio Claro na extensão Geopackage
sf::st_write(obj = geo_vetor_rio_claro,
dsn = here::here("dados", "vetor", "vetores.gpkg"),
layer = "rio_claro")Ainda sobre o formato GeoPackage, há algo muito interessante que podemos fazer: podemos acrescentar outros arquivos vetoriais ao mesmo arquivo já criado. Como exemplo, exportaremos o limite do Brasil para o mesmo arquivo.
## Exportar o polígono do Brasil na extensão Geopackage
sf::st_write(obj = geo_vetor_brasil,
dsn = here::here("dados", "vetor", "vetores.gpkg"),
layer = "brasil")Raster
Para exportar dados raster utilizamos geralmente a função raster::writeRaster(). Exportar dados raster é um pouco mais complexo que exportar dados vetoriais. Teremos de definir se exportaremos arquivos em uma ou várias camadas, quantidade de informações por pixel, e ainda diferentes extensões de saída.
📝 Importante
Arquivos raster escritos em discos geralmente ocupam bastante espaço, e dessa forma, há parâmetros específicos para certos tipos de dados, que detalharemos a seguir para contornar esse problema e comprimir os arquivos.
Na função raster::writeRaster(), o argumento x diz respeito ao objeto raster no ambiente R. O argumento filename é nome do arquivo que será exportado do R, podendo ou não possuir a extensão que se pretende que o arquivo tenha. O argumento format é o formato do arquivo, sendo as principais possibilidades resumidas na Tabela 15.8, e para saber das possibilidades suportadas, use a função raster::writeFormats(). O argumento bylayer diz se múltiplas camadas serão exportadas em arquivos diferentes ou em apenas um arquivo.
| Tipo de arquivo | Nome longo | Extensão | Suporte a múltiplas camadas |
|---|---|---|---|
| raster | Formato pacote raster | .grd | Sim |
| ascii | ESRI Ascii | .asc | Não |
| SAGA | SAGA GIS | .sdat | Não |
| IDRISI | IDRISI | .rst | Não |
| CDF | netCDF (requer ncdf4) | .nc | Sim |
| GTiff | GeoTiff (requer rgdal) | .tif | Sim |
| ENVI | ENVI .hdr | .envi | Sim |
| EHdr | ESRI .hdr | .bil | Sim |
| HFA | Erdas imagem (.img) | .img | Sim |
Dentre os argumentos adicionais, temos ainda o datatype, que faz referência a um dos nove tipos de formato de dados detalhados na Tabela 15.9, sendo que o tipo de dado determina a representação em bits (quantidade de informação) na célula do objeto raster exportado e depende da faixa de valores do objeto raster em cada pixel. Quanto mais valores um tipo de dado puder representar, maior será o arquivo exportado no disco. Dessa forma, é interessante utilizar um tipo de dado que diminua o tamanho do arquivo a ser exportado, dependendo do tipo de dado em cada pixel. Para a função raster::writeRaster(), o default é FLT4S, o que pode ocupar mais espaço em disco do que o necessário.
| Tipo de dado | Valor mínimo | Valor máximo |
|---|---|---|
| LOG1S | FALSE (0) | TRUE (1) |
| INT1S | -127 | 127 |
| INT1U | 0 | 255 |
| INT2S | -32.767 | 32.767 |
| INT2U | 0 | 65534 |
| INT4S | -2.147.483.647 | 2.147.483.647 |
| INT4U | 0 | 42.94.967.296 |
| FLT4S | -3,4e+38 | 3,4e+38 |
| FLT8S | -1,7e+308 | 1,7e+308 |
Outros argumentos de suporte são: overwrite para sobrescrever um arquivo que já exista, progress para mostrar uma barra de progresso da exportação como “text” ou “window”, e options que permite opções do GDAL. Para esse último, quando exportar especificamente na extensão GeoTIFF, podemos utilizar options = c("COMPRESS=DEFLATE", "TFW=YES") para que haja compressão do arquivo, diminuindo consideravelmente seu tamanho (cerca de um terço), aliado à criação de um arquivo auxiliar .tfw, para ser carregado em softwares específicos de SIG, como o ArcGIS.
Para exportar apenas uma camada RasterLayer, podemos utilizar a função raster::writeRaster() em um formato mais simples.
## Criar diretório
dir.create(here::here("dados", "raster", "exportados"))
## Exportar raster layer
raster::writeRaster(geo_raster_srtm,
filename = here::here("dados", "raster", "exportados", "elevation"),
format = "GTiff",
datatype = "INT2S",
options = c("COMPRESS=DEFLATE", "TFW=YES"),
progress = "text",
overwrite = TRUE)Para mais de uma camada RasterBrick ou RasterStack, podemos utilizar a função raster::writeRaster() com o bylayer = TRUE.
## Exportar raster stack
raster::writeRaster(x = geo_raster_bioclim,
filename = here::here("dados", "raster", "exportados", names(geo_raster_bioclim)),
bylayer = TRUE,
format = "GTiff",
datatype = "INT2S",
options = c("COMPRESS=DEFLATE", "TFW=YES"),
progress = "text",
overwrite = TRUE)15.7 Descrição de objetos geoespaciais
Muitas vezes precisaremos verificar as informações dos objetos geoespaciais importados para o R. Apesar de chamar o objeto trazer grande parte das informações que precisamos consultar, existem funções específicas que nos auxiliam nesse processo de descrição dos objetos.
15.7.1 Vetor
Podemos acessar as informações geoespaciais e a tabela de atributos de um objeto importado como vetor simplesmente chamando o nome do objeto no R.
## Município de Rio Claro
geo_vetor_rio_claro
#> Simple feature collection with 1 feature and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -47.76521 ymin: -22.55203 xmax: -47.46188 ymax: -22.24368
#> Geodetic CRS: SIRGAS 2000
#> code_muni name_muni code_state abbrev_state name_state code_region name_region geom
#> 493 3543907 Rio Claro 35 SP São Paulo 3 Sudeste MULTIPOLYGON (((-47.46875 -...Mas também podemos acessar informações geoespaciais com funções específicas, como tipo de geometria, limites geoespaciais do vetor (extensão), sistema de referência de coordenadas (CRS), e a tabela de atributos.
## Tipo de geometria
sf::st_geometry_type(geo_vetor_rio_claro)
#> [1] MULTIPOLYGON
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON GEOMETRYCOLLECTION CIRCULARSTRING ... TRIANGLE
## Extensão
sf::st_bbox(geo_vetor_rio_claro)
#> xmin ymin xmax ymax
#> -47.76521 -22.55203 -47.46188 -22.24368
## CRS
sf::st_crs(geo_vetor_rio_claro)
#> Coordinate Reference System:
#> User input: SIRGAS 2000
#> wkt:
#> GEOGCRS["SIRGAS 2000",
#> DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
#> ELLIPSOID["GRS 1980",6378137,298.257222101,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
#> BBOX[-59.87,-122.19,32.72,-25.28]],
#> ID["EPSG",4674]]
## Acessar a tabela de atributos
geo_vetor_rio_claro_tab <- sf::st_drop_geometry(geo_vetor_rio_claro)
geo_vetor_rio_claro_tab
#> code_muni name_muni code_state abbrev_state name_state code_region name_region
#> 493 3543907 Rio Claro 35 SP São Paulo 3 Sudeste15.7.2 Raster
Da mesma forma, podemos acessar as informações objetos raster chamando o nome do objeto.
## Raster layer
geo_raster_srtm
#> class : RasterLayer
#> dimensions : 6000, 6000, 3.6e+07 (nrow, ncol, ncell)
#> resolution : 0.0008333333, 0.0008333333 (x, y)
#> extent : -50, -45, -25, -20 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : srtm_27_17.tif
#> names : srtm_27_17
#> values : -32768, 32767 (min, max)Além disso, podemos selecionar informações desse objeto com funções específicas, tanto para RasterLayer, quanto para RasterBrick ou RasterStack como: classe, dimensões (número de linhas, colunas e camadas), número de camadas, número de linhas, número de colunas, número de células, resolução (largura e altura do tamanho do pixel), extensão (limites geoespaciais), sistema de referência de coordenadas (CRS), nome das camadas e extrair os valores de todos os pixels.
## Classe
class(geo_raster_srtm)
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"
## Dimensões
dim(geo_raster_srtm)
#> [1] 6000 6000 1
## Número de camadas
nlayers(geo_raster_srtm)
#> [1] 1
## Número de linhas
nrow(geo_raster_srtm)
#> [1] 6000
## Número de colunas
ncol(geo_raster_srtm)
#> [1] 6000
## Número de células
ncell(geo_raster_srtm)
#> [1] 3.6e+07
## Resolução
res(geo_raster_srtm)
#> [1] 0.0008333333 0.0008333333
## Extensão
extent(geo_raster_srtm)
#> class : Extent
#> xmin : -50
#> xmax : -45
#> ymin : -25
#> ymax : -20
## Projeção ou CRS
projection(geo_raster_srtm)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
## Nomes
names(geo_raster_srtm)
#> [1] "srtm_27_17"
## Valores
getValues(geo_raster_srtm) %>% head
#> [1] 382 379 379 379 379 383
values(geo_raster_srtm) %>% head
#> [1] 382 379 379 379 379 383
geo_raster_srtm[] %>% head
#> [1] 382 379 379 379 379 38315.8 Reprojeção de dados geoespaciais
Em algumas situações é necessário alterar o CRS de um objeto espacial para um novo CRS. A reprojeção é justamente a transformação de coordenadas de um CRS para outro: geoespaciais (‘lon/lat’, com unidades em graus de longitude e latitude) e projetados (normalmente com unidades de metros a partir de um datum).
Geralmente precisaremos fazer essa operação para transformar camadas vetoriais ou rasters para o mesmo CRS, de modo que possam ser exibidas conjuntamente, ou ainda que as camadas possuem CRS projetado para realizar alguma operação espacial entre camadas, ou quando precisamos calcular áreas, formatos ou distâncias, como métricas de paisagem, por exemplo. Existe uma infinidade de projeções e um excelente material de consulta é o livro de Lapaine & Usery (2017).
Podemos verificar o CRS de uma camada através da função sf::st_crs() ou raster::projection() e raster::crs(), ou ainda, saber se a mesma possui um CRS geográfico ou não, com a função sf::st_is_longlat().
Já para reprojetar um objeto sf usamos a função sf::st_transform() e para um objeto raster usamos a função raster::projectRaster().
## Projeção de vetores
sf::st_crs(geo_vetor_rio_claro)
#> Coordinate Reference System:
#> User input: SIRGAS 2000
#> wkt:
#> GEOGCRS["SIRGAS 2000",
#> DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
#> ELLIPSOID["GRS 1980",6378137,298.257222101,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
#> BBOX[-59.87,-122.19,32.72,-25.28]],
#> ID["EPSG",4674]]
## Projeção de raster
raster::projection(geo_raster_srtm)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
raster::crs(geo_raster_srtm)
#> Coordinate Reference System:
#> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
#> WKT2 2019 representation:
#> GEOGCRS["WGS 84 (with axis order normalized for visualization)",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> REMARK["Axis order reversed compared to EPSG:4326"]]
## Verificar se o CRS é geográfico
sf::st_is_longlat(geo_vetor_rio_claro)
#> [1] TRUEAs funções sf::st_transform() e raster::projectRaster() possuem dois argumentos importantes: x que é o objeto a ser reprojetado e o crs que é o CRS alvo. O argumento crs pode ser especificado de quatro maneiras: i) código EPSG (por exemplo, 4326), ii) string PROJ4 (por exemplo, + proj = longlat + datum = WGS84 + no_defs), iii) string WKT, ou iv) objeto crs de outra camada, conforme retornado por sf::st_crs() ou raster::crs(). Essas informações de EPSG, PROJ4 e WKT podem ser acessadas nas bases: epsg.io e spatialreference.org.
Dentre os possíveis CRSs a serem utilizados, alguns são mais comuns para CRSs geoespaciais e projetados. Para CRSs geoespaciais, o mais comum para o mundo é o World Geodetic System 1984 (WGS84), ou seja, geográfico com datum WGS84. Para o Brasil, o CRS adotado é o Sistema de Referencia Geocéntrico para las Américas 2000 (SIRGAS 2000), ou seja, geográfico com datum SIRGAS2000.
Para CRSs projetados, essa escolha vai depender da extensão e localização da área de interesse no globo terrestre. Aqui destacaremos os principais, para três escalas: global, regional e local. Para a escala global, geralmente usa-se umas dessas projeções, dependendo do objetivo: i) Projeção de Mollweide, ii) Projeção de Winkel Tripel, iii) Projeção de Eckert IV, iv) Projeção Azimutal de Lambert. Para a escala regional, como um hemisfério, geralmente usa-se a Projeção Cônica de Albers. Por fim, para a escala local, usa-se geralmente a Projeção Universal Transverse Mercator (UTM), um conjunto de CRSs que divide a Terra em 60 cunhas longitudinais e 20 segmentos latitudinais, como pode ser visto neste link.
Os principais CRSs são descritos na Tabela 15.10.
| CRS | Tipo de CRS | Descrição | epsg.io | spatialreference.org |
|---|---|---|---|---|
| World Geodetic System 1984 (WGS84) | Geográfico | CRS geográfico mais comum para o mundo | EPSG:4326 | EPSG:4326 |
| Sistema de Referencia Geocéntrico para las Américas 2000 (SIRGAS 2000) | Geográfico | CRS geográfico oficial para o Brasil | EPSG:4674 | EPSG:4674 |
| Projeção de Mollweide | Projetado | CRS projetado que preserva as relações de área | ESRI:54009 | SR-ORG:7099 |
| Projeção de Winkel Tripel | Projetado | CRS projetado com mínimo de distorção para área, direção e distância | NA | SR-ORG:7291 |
| Projeção de Eckert IV | Projetado | CRS projetado que preserva a área e com meridianos elípticos | EPSG:54012 | ESRI:54012 |
| Projeção Azimutal de Lambert | Projetado | CRS projetado que preserva os tamanhos relativos e senso de direção a partir do centro | NA | NA |
| Projeção Cônica de Albers | Projetado | CRS projetado para escala regional, mantendo a área constante em toda sua superfície | NA | SR-ORG:7823 |
| Projeção Universal Transverse Mercator (UTM) | Projetado | CRS projetado para escala local, distorcendo áreas e distâncias com gravidade crescente com a distância do centro da zona UTM | EPSG:31983 | EPSG:31983 |
15.8.1 Vetor
Como dissemos, para reprojetar um vetor, utilizamos a função sf::st_transform(), observando os argumentos x que é a camada a ser reprojetada, e o crs que é o CRS alvo.
Vamos reprojetar o limite do município de Rio Claro/SP do CRS SIRGAS2000/geográfico para o CRS projetado SIRGAS2000/UTM23S, com os efeitos da transformação podendo ser notados na Figura 15.16.
## Converter CRS
geo_vetor_rio_claro_sirgas2000_utm23s <- sf::st_transform(x = geo_vetor_rio_claro,
crs = 31983)
Figura 15.16: Limites do município de Rio Claro/SP com CRS SIRGAS2000/geográfico e com CRS SIRGAS2000/UTM23S.
Podemos ainda utilizar o formato proj4string no argumento crs para fazer a transformação. Vamos primeiramente plotar o mundo em WGS84/Geográfico (Figura 15.17).
## Plot
plot(co110_sf[1], col = "gray", main = "WGS84/Geográfio", graticule = TRUE)
Figura 15.17: Limite dos países do mundo com CRS geográfico e datum WGS84.
Agora, reprojetaremos utilizando a Projeção de Mollweide (Figura 15.18).
## Projeção de Mollweide
co110_sf_moll <- sf::st_transform(x = co110_sf, crs = "+proj=moll")
## Plot
plot(co110_sf_moll[1], col = "gray", main = "Projeção de Mollweide", graticule = TRUE)
Figura 15.18: Limite dos países do mundo com CRS Projeção de Mollweide.
Ou ainda podemos utilizar a Projeção Azimutal de Lambert com alguns parâmetros ajustados para centralizar a projeção no Brasil (Figura 15.19).
## Projeção Azimutal de Lambert
co110_sf_laea <- sf::st_transform(x = co110_sf,
crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=-50 +lat_0=0")
## Plot
plot(co110_sf_laea[1], col = "gray", main = "Projeção Azimutal de Lambert", graticule = TRUE)
Figura 15.19: Limite dos países do mundo com CRS Projeção Azimutal de Lambert centrado no Brasil.
15.8.2 Raster
A reprojeção de objetos raster não é uma tarefa tão simples quanto a reprojeção de vetores. Em vetores, a reprojeção altera as coordenadas de cada vértice. Entretanto, como rasters são compostos de células retangulares do mesmo tamanho,a reprojeção do raster envolve a criação de um novo objeto raster, com duas operações espaciais separadas: i) reprojeção vetorial dos centroides celulares para outro CRS (i.e., muda a posição e tamanho do pixel) e, ii) cálculo de novos valores do pixel por meio de reamostragem (i.e., muda o valor do pixel).
A função raster::projectRaster() possui alguns parâmetros que necessitam de algumas especificações. O argumento from é o objeto raster de entrada que sofre a reprojeção. O argumento to é um objeto raster do qual todas as propriedades dos CRSs, como extensão e resolução serão associadas ao objeto raster indicado no argumento from. O argumento res permite ajustar a resolução do pixel de saída do objeto raster reprojetado.
O argumento crs aceita apenas as definições de proj4string extensas de um CRS em vez de códigos EPSG concisos. Contudo, é possível usar um código EPSG em uma definição de proj4string com +init=epsg:EPSG. Por exemplo, pode-se usar a definição +init=epsg:4326 para definir CRS para WGS84 (código EPSG de 4326). A biblioteca PROJ adiciona automaticamente o resto dos parâmetros e os converte em +init=epsg:4326 +proj=longlat +datum=WGS84 + no_defs + ellps=WGS84 + towgs84=0,0,0.
O argumento method permite escolher entre os métodos ngb (vizinho mais próximo) ou biliniar (interpolação bilinear), sendo o primeiro mais indicado para reprojeção de rasters categóricos, pois os valores estimados devem ser iguais aos do raster original. O método ngb define cada novo valor de célula para o valor da célula mais próxima (centro) do raster de entrada. Já o método biliniar é indicado para raster contínuos e calcula o valor da célula de saída com base nas quatro células mais próximas no raster original, sendo a média ponderada da distância dos valores dessas quatro células. Existem outras formas de interpolação, mas não as abordaremos aqui.
Aqui, vamos reprojetar os dados de elevação para Rio Claro/SP. Para que esse processo seja mais rápido, iremos antes ajustar a extensão do raster para o limite do município usando a função raster::crop() (Figura 15.20). Essa função é melhor explicada na seção de cortes e máscaras, mais adiante.
## Ajuste do limite
geo_raster_srtm_rio_claro <- raster::crop(x = geo_raster_srtm,
y = geo_vetor_rio_claro)
geo_raster_srtm_rio_claro
#> class : RasterLayer
#> dimensions : 370, 364, 134680 (nrow, ncol, ncell)
#> resolution : 0.0008333333, 0.0008333333 (x, y)
#> extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : srtm_27_17
#> values : 491, 985 (min, max)
## Plot
plot(geo_raster_srtm_rio_claro, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.20: Ajuste da extensão do raster de elevação para o município de Rio Claro/SP.
Primeiramente, vamos reprojetar indicando uma projeção e sem especificar o tamanho da célula. Note que o tamanho da célula vai se ajustar para valores diferentes, sendo portanto, pixels retangulares e não quadrados.
## Reprojeção
geo_raster_srtm_rio_claro_sirgas2000_utm23s <- raster::projectRaster(
from = geo_raster_srtm_rio_claro,
crs = "+init=epsg:31983",
method = "bilinear")
geo_raster_srtm_rio_claro_sirgas2000_utm23s
#> class : RasterLayer
#> dimensions : 386, 381, 147066 (nrow, ncol, ncell)
#> resolution : 85.8, 92.3 (x, y)
#> extent : 214575.4, 247265.2, 7503009, 7538637 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
#> source : memory
#> names : srtm_27_17
#> values : 491.6033, 980.4151 (min, max)Agora vamos reprojetar especificando o tamanho da célula (Figura 15.21). Dessa forma, todas as células terão o mesmo, i.e., quadrados de 90 metros.
## Reprojeção
geo_raster_srtm_rio_claro_sirgas2000_utm23s <- raster::projectRaster(
from = geo_raster_srtm_rio_claro,
crs = "+init=epsg:31983",
method = "bilinear",
res = 90)
geo_raster_srtm_rio_claro_sirgas2000_utm23s
#> class : RasterLayer
#> dimensions : 396, 364, 144144 (nrow, ncol, ncell)
#> resolution : 90, 90 (x, y)
#> extent : 214554.4, 247314.4, 7502985, 7538625 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
#> source : memory
#> names : srtm_27_17
#> values : 493.2395, 986.686 (min, max)
## Plot
plot(geo_raster_srtm_rio_claro_sirgas2000_utm23s,
col = viridis::viridis(10))
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom,
col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.21: Reprojeção do raster de elevação para SIRGAS2000/UTM23S especificado por um objeto e informando o tamanho da célula.
Vamos também reprojetar uma camada mundial da média de temperatura anual (BIO01), indicando o tamanho da célula para 25.000 m (Figura 15.22).
## Reprojeção
geo_raster_bioclim_moll <- raster::projectRaster(
from = geo_raster_bioclim[[1]],
crs = "+proj=moll",
res = 25000,
method = "bilinear")
geo_raster_bioclim_moll
#> class : RasterLayer
#> dimensions : 732, 1453, 1063596 (nrow, ncol, ncell)
#> resolution : 25000, 25000 (x, y)
#> extent : -18159905, 18165095, -9154952, 9145048 (xmin, xmax, ymin, ymax)
#> crs : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> source : memory
#> names : wc2.1_10m_bio_1
#> values : -54.66752, 30.71805 (min, max)
## Plot
plot(geo_raster_bioclim_moll, col = viridis::viridis(10))
plot(co110_sf_moll[1], col = NA, add = TRUE)
Figura 15.22: Reprojeção do raster de média de temperatura anual (BIO01) para Projeção de Mollweide informando o tamanho da célula.
15.9 Principais operações com dados geoespaciais
Nesta seção veremos as principais funções para realizar operações com dados geoespaciais. Essas operações são separadas conforme Lovelace et al. (2019) em: Operações de atributos, Operações espaciais, e Operações geométricas.
15.9.1 Operações de atributos
São modificação de objetos geoespaciais baseado em informações não espaciais, como a tabela de atributos ou valores das células e nome das camadas dos rasters.
Vetor
As principais operações de atributos vetoriais são com respeito à tabela de atributos, sendo as principais: i) filtro, ii) junção, iii) agregação e iv) manipulação da tabela de atributos. A lista de possíveis operações é longa, dessa forma, apresentaremos algumas operações utilizando as principais funções e listamos as demais funções e suas operações, que dependerão de objetivos específicos.
Quase todas as operações serão as mesmas realizadas pelo pacote dplyr em uma tabela de dados (ver o Capítulo 5), sendo algumas operações específicas para alterar apenas campos da tabela de atributos e outras que refletem operações nas feições, ou seja, alterarão através da tabela de atributos as características das feições. Essas funções e suas operações são descritas com detalhes na Tabela 15.11.
| Funções | Onde atua | Descrição |
|---|---|---|
filter()
|
Feições | Selecionar feições por valores |
slice()
|
Feições | Selecionar feições pela posição na tabela de atributos |
n_sample()
|
Feições | Amostrar feições na tabela de atributos |
group_by()
|
Feições | Agrupar feições por valores da tabela de atributos |
summarise()
|
Feições | Operações com valores das feições na tabela de atributos, que acabam por dissolver as feições |
select()
|
Atributos | Selecionar colunas da tabela de atributos |
pull()
|
Atributos | Selecionar uma coluna da tabela de atributos como vetor |
rename()
|
Atributos | Renomear uma coluna da tabela de atributos |
mutate()
|
Atributos | Criar uma coluna ou alterar os valores da tabela de atributos |
*_join()
|
Atributos | Diversas funções para juntar dados de outras tabelas de dados à tabela de atributos |
Para exemplificar as operações de atributos, vamos utilizar os dados de nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Filtro
Vamos iniciar as operações fazendo o filtro de feições pela tabela de atributos, que permite selecionar feições pelos seus valores atribuídos, utilizando a função dplyr::filter(). Aqui vamos selecionar as feições de floresta do mapa de cobertura da terra para Rio Claro/SP (Figura 15.23).
## Filtro
geo_vetor_cobertura_floresta <- geo_vetor_cobertura %>%
dplyr::filter(CLASSE_USO == "formação florestal")
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_cobertura_floresta$geometry, col = "forestgreen", add = TRUE)
Figura 15.23: Filtro da classe floresta para o mapeamento de cobertura da terra para o município de Rio Claro/SP.
Junção
Uma das funções mais úteis de operações de atributos é a junção, referida em inglês como join, realizada através das funções dplyr::*_join() (ver detalhes do Capítulo 5). Nela, usamos uma coluna identificadora para atribuir dados de outra tabela de dados. Como exemplo, vamos criar uma tabela de dados com novos nomes das classes de cobertura da terra e atribuir esses novos nomes à tabela de atributos do objeto vetorial. É fundamental destacar que para que essa função funcione, precisamos de uma coluna identificadora dos valores para que a junção seja possível.
## Dados
dados_classes <- tibble::tibble(
CLASSE_USO = geo_vetor_cobertura$CLASSE_USO,
classe = c("agua", "antropico", "edificado", "floresta", "silvicultura"))
dados_classes
#> # A tibble: 5 × 2
#> CLASSE_USO classe
#> <chr> <chr>
#> 1 água agua
#> 2 área antropizada antropico
#> 3 área edificada edificado
#> 4 formação florestal floresta
#> 5 silvicultura silvicultura
## Junção
geo_vetor_cobertura_classes <- dplyr::left_join(
x = geo_vetor_cobertura,
y = dados_classes,
by = "CLASSE_USO") %>%
sf::st_drop_geometry()
geo_vetor_cobertura_classes
#> GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA classe
#> 1 3543907 RIO CLARO SP 35 água 357.027 agua
#> 2 3543907 RIO CLARO SP 35 área antropizada 37297.800 antropico
#> 3 3543907 RIO CLARO SP 35 área edificada 5078.330 edificado
#> 4 3543907 RIO CLARO SP 35 formação florestal 7017.990 floresta
#> 5 3543907 RIO CLARO SP 35 silvicultura 138.173 silviculturaAgregação
Outra função bastante útil é a agregação de atributos. Apesar de existir uma função que realiza a união de feições que veremos na próxima seção, o uso conjunto das funções dplyr::group_by() e dplyr::summarise() realizam uma tarefa semelhante. Aqui vamos agregar as nascentes para Rio Claro/SP, i.e., juntar cada ponto que estava numa linha da tabela de atributos de modo que todos fiquem numa mesma linha, com o valor da quantidade de nascentes (Figura 15.24).
## Agregar
geo_vetor_nascentes_n <- geo_vetor_nascentes %>%
dplyr::group_by(MUNICIPIO, HIDRO) %>%
dplyr::summarise(n = n())
geo_vetor_nascentes_n
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: MULTIPOINT
#> Dimension: XY
#> Bounding box: xmin: 217622.9 ymin: 7504132 xmax: 246367.4 ymax: 7537855
#> Projected CRS: SIRGAS 2000 / UTM zone 23S
#> # A tibble: 1 × 4
#> # Groups: MUNICIPIO [1]
#> MUNICIPIO HIDRO n geometry
#> <chr> <chr> <int> <MULTIPOINT [m]>
#> 1 RIO CLARO nascente 1220 ((217622.9 7528315), (217836.5 7528103), (217988.9 7528203), (218288.9 7528237), (2183...
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom,
col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_nascentes_n$geometry, pch = 20,
col = "blue", add = TRUE)
Figura 15.24: Agregação e contagem das nascentes para o município de Rio Claro/SP.
Manipulação da tabela de atributos
Por fim, é muito comum em análises de softwares SIG a criação ou atualização dos valores de colunas na tabela de atributos. Aqui, podemos utilizar a função dplyr::mutate() para criar essas novas colunas, assim como atualizar os valores de colunas existentes. Em nosso exemplo, faremos uma composição das colunas CLASSE_USO e AREA_HA numa nova coluna chamada classe_area.
## Criar coluna
geo_vetor_cobertura_cob_col_area <- geo_vetor_cobertura %>%
dplyr::mutate(classe_area = paste0(CLASSE_USO, " (", AREA_HA, " ha)")) %>%
sf::st_drop_geometry()
geo_vetor_cobertura_cob_col_area
#> GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA classe_area
#> 1 3543907 RIO CLARO SP 35 água 357.027 água (357.027 ha)
#> 2 3543907 RIO CLARO SP 35 área antropizada 37297.800 área antropizada (37297.8 ha)
#> 3 3543907 RIO CLARO SP 35 área edificada 5078.330 área edificada (5078.33 ha)
#> 4 3543907 RIO CLARO SP 35 formação florestal 7017.990 formação florestal (7017.99 ha)
#> 5 3543907 RIO CLARO SP 35 silvicultura 138.173 silvicultura (138.173 ha)Duas funções são bastante interessantes de serem integradas junto à manipulação de tabelas de atributos. Elas calculam propriedades geométricas numéricas dos vetores de linhas (comprimento) e polígonos (área): sf::st_length() e sf::st_area(). Essas funções calculam essas propriedades em metros para comprimento e em metros quadrados para área, independentemente do CRS. Para tanto, vamos utilizar as linhas de hidrografia e os polígonos de cobertura da terra para Rio Claro/SP, e atribuir esses valores à tabela de atributos de ambos os objetos geoespaciais, utilizando em conjunto a função dplyr::mutate().
## Comprimento das linhas
geo_vetor_hidrografia_comp <- geo_vetor_hidrografia %>%
dplyr::mutate(comprimento = sf::st_length(.))
geo_vetor_hidrografia_comp
#> Simple feature collection with 1 feature and 7 fields
#> Geometry type: MULTILINESTRING
#> Dimension: XY
#> Bounding box: xmin: 215155.3 ymin: 7504132 xmax: 246367.4 ymax: 7537978
#> Projected CRS: SIRGAS 2000 / UTM zone 23S
#> GEOCODIGO MUNICIPIO UF CD_UF HIDRO COMP_KM geometry comprimento
#> 1 3543907 RIO CLARO SP 35 curso d'água (0 - 10m) 1142.98 MULTILINESTRING ((231815.7 ... 1142981 [m]
## Área dos polígonos
geo_vetor_cobertura_area <- geo_vetor_cobertura %>%
dplyr::mutate(area_m2 = sf::st_area(.))
geo_vetor_cobertura_area
#> Simple feature collection with 5 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 215151.7 ymin: 7503723 xmax: 246582.4 ymax: 7537978
#> Projected CRS: SIRGAS 2000 / UTM zone 23S
#> GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA geometry area_m2
#> 1 3543907 RIO CLARO SP 35 água 357.027 MULTIPOLYGON (((235487.6 75... 3570267 [m^2]
#> 2 3543907 RIO CLARO SP 35 área antropizada 37297.800 MULTIPOLYGON (((232275 7504... 372978415 [m^2]
#> 3 3543907 RIO CLARO SP 35 área edificada 5078.330 MULTIPOLYGON (((233123.6 75... 50783283 [m^2]
#> 4 3543907 RIO CLARO SP 35 formação florestal 7017.990 MULTIPOLYGON (((232355 7504... 70179895 [m^2]
#> 5 3543907 RIO CLARO SP 35 silvicultura 138.173 MULTIPOLYGON (((243052.1 75... 1381726 [m^2]Raster
Devido à estrutura espacial do raster ser formada por uma ou mais superfícies contínuas, as manipulações como subconjunto e outras operações em objetos raster funcionam de uma maneira diferente do que em objetos vetoriais. Veremos aqui as três principais: i) subconjunto de células usando o operador [] ou subconjunto de camadas RasterStack ou RasterBrick utilizando os operadores [[]] e $, ii) renomear nomes das camadas, e iii) resumir informações de todos os pixels.
Subconjunto
Podemos fazer um subconjunto de células utilizando dentro dos operadores [] valores para indicar a posição da linha e da coluna de um raster, ou ainda a posição de uma célula utilizando apenas um número. Essas operações resultarão em valores diferentes para RasterLayer e RasterBrick ou RasterStack.
## Raster - linha 1 e columna 1
geo_raster_srtm[1, 1]
#>
#> 382
## Raster - célula 1
geo_raster_srtm[1]
#>
#> 382
## Stack - linha 1 e columna 1
geo_raster_bioclim[1, 1]
#> wc2.1_10m_bio_1 wc2.1_10m_bio_10 wc2.1_10m_bio_11 wc2.1_10m_bio_12 wc2.1_10m_bio_13 wc2.1_10m_bio_14 wc2.1_10m_bio_15 wc2.1_10m_bio_16
#> [1,] NA NA NA NA NA NA NA NA
#> wc2.1_10m_bio_17 wc2.1_10m_bio_18 wc2.1_10m_bio_19 wc2.1_10m_bio_2 wc2.1_10m_bio_3 wc2.1_10m_bio_4 wc2.1_10m_bio_5 wc2.1_10m_bio_6
#> [1,] NA NA NA NA NA NA NA NA
#> wc2.1_10m_bio_7 wc2.1_10m_bio_8 wc2.1_10m_bio_9
#> [1,] NA NA NA
## Stack - célula 1
geo_raster_bioclim[1]
#> wc2.1_10m_bio_1 wc2.1_10m_bio_10 wc2.1_10m_bio_11 wc2.1_10m_bio_12 wc2.1_10m_bio_13 wc2.1_10m_bio_14 wc2.1_10m_bio_15 wc2.1_10m_bio_16
#> [1,] NA NA NA NA NA NA NA NA
#> wc2.1_10m_bio_17 wc2.1_10m_bio_18 wc2.1_10m_bio_19 wc2.1_10m_bio_2 wc2.1_10m_bio_3 wc2.1_10m_bio_4 wc2.1_10m_bio_5 wc2.1_10m_bio_6
#> [1,] NA NA NA NA NA NA NA NA
#> wc2.1_10m_bio_7 wc2.1_10m_bio_8 wc2.1_10m_bio_9
#> [1,] NA NA NAPara selecionar uma camada de um RasterBrick ou RasterStack, podemos utilizar as funções raster::subset() ou raster::raster() com o argumento layer indicando a ordem ou o nome da camada, além dos operadores [[]] e $ (Figura 15.25).
## Seleção de camada num objeto stack utilizando a função subset
geo_raster_bioclim_bio01 <- raster::subset(geo_raster_bioclim, "wc2.1_10m_bio_1")
geo_raster_bioclim_bio01
#> class : RasterLayer
#> dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : wc2.1_10m_bio_1.tif
#> names : wc2.1_10m_bio_1
#> values : -54.72435, 30.98764 (min, max)
## Seleção de camada num objeto stack utilizando a função raster
geo_raster_bioclim_bio01 <- raster::raster(geo_raster_bioclim, layer = 1)
geo_raster_bioclim_bio01
#> class : RasterLayer
#> dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : wc2.1_10m_bio_1.tif
#> names : wc2.1_10m_bio_1
#> values : -54.72435, 30.98764 (min, max)
## Seleção de camada num objeto stack utilizando os operadores [[]] e o nome
geo_raster_bioclim_bio01 <- geo_raster_bioclim[["wc2.1_10m_bio_1"]]
geo_raster_bioclim_bio01
#> class : RasterLayer
#> dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : wc2.1_10m_bio_1.tif
#> names : wc2.1_10m_bio_1
#> values : -54.72435, 30.98764 (min, max)
## Seleção de camada num objeto stack utilizando os operadores [[]] e a posicao
geo_raster_bioclim_bio01 <- geo_raster_bioclim[[1]]
geo_raster_bioclim_bio01
#> class : RasterLayer
#> dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : wc2.1_10m_bio_1.tif
#> names : wc2.1_10m_bio_1
#> values : -54.72435, 30.98764 (min, max)
## Seleção de camada num objeto stack utilizando o operador $
geo_raster_bioclim_bio01 <- geo_raster_bioclim$wc2.1_10m_bio_1
geo_raster_bioclim_bio01
#> class : RasterLayer
#> dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : wc2.1_10m_bio_1.tif
#> names : wc2.1_10m_bio_1
#> values : -54.72435, 30.98764 (min, max)
Figura 15.25: Camada BIO01 selecionada pelas operações de subconjunto.
Renomear
Podemos ainda renomear camadas dos raster RasterLayer utilizando a função names().
## Raster - nomes
names(geo_raster_srtm_rio_claro)
#> [1] "srtm_27_17"
## Raster - renomear
names(geo_raster_srtm_rio_claro) <- "elevacao"
## Raster - nomes
names(geo_raster_srtm_rio_claro)
#> [1] "elevacao"E essa operação também funciona para RasterBrick e RasterStack.
## Stack - nomes
names(geo_raster_bioclim)
#> [1] "wc2.1_10m_bio_1" "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12" "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15"
#> [8] "wc2.1_10m_bio_16" "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19" "wc2.1_10m_bio_2" "wc2.1_10m_bio_3" "wc2.1_10m_bio_4"
#> [15] "wc2.1_10m_bio_5" "wc2.1_10m_bio_6" "wc2.1_10m_bio_7" "wc2.1_10m_bio_8" "wc2.1_10m_bio_9"
## Stack - renomear
names(geo_raster_bioclim) <- c("bio01", paste0("bio", 10:19), paste0("bio0", 2:9))
## Stack - nomes
names(geo_raster_bioclim)
#> [1] "bio01" "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16" "bio17" "bio18" "bio19" "bio02" "bio03" "bio04" "bio05" "bio06" "bio07"
#> [18] "bio08" "bio09"Resumir
Muitas vezes queremos fazer cálculos para todos as células de um raster. Podemos resumir informações de todos os pixels fazendo cálculos simples com todos os pixels de cada camada com a função raster::cellStats(), sendo x o argumento do objeto raster e stat o nome da função resumo, como “mean” ou “sum”.
## Raster - média de todas as células de altitude
raster::cellStats(x = geo_raster_srtm_rio_claro, stat = mean)
#> [1] 625.8273
## Stack - média de todas as células de cada camada bioclimática
raster::cellStats(x = geo_raster_bioclim, stat = mean)
#> bio01 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio02
#> -4.0378283 7.2035545 -13.8963286 550.0569022 93.4633916 15.3689993 74.7084151 241.6525005 55.4149542 156.4237816 108.8950626 9.9432120
#> bio03 bio04 bio05 bio06 bio07 bio08 bio09
#> 34.5221528 880.1215546 13.9386423 -19.7938943 33.7325366 -0.9226276 -5.3774489Ou ainda, podemos analisar a frequência com que cada valor dos pixels ocorre, utilizando a função raster::freq().
## Raster - frequência das células
raster::freq(x = geo_raster_srtm_rio_claro) %>% head()
#> value count
#> [1,] 491 1
#> [2,] 492 4
#> [3,] 493 9
#> [4,] 494 19
#> [5,] 495 32
#> [6,] 496 44
## Stack - frequência das células
raster::freq(x = geo_raster_bioclim[[1]]) %>% head()
#> value count
#> [1,] -55 319
#> [2,] -54 4529
#> [3,] -53 5778
#> [4,] -52 6128
#> [5,] -51 6090
#> [6,] -50 789215.9.2 Operações espaciais
As operações espaciais são modificações de objetos geoespaciais baseado em informações espaciais, como localização e formato. Seria quase impossível abordar todas as operações realizáveis nesse capítulo, então demonstraremos as principais para dados vetoriais e raster.
Vetor
As principais operações espaciais para dados vetoriais são: i) filtro espacial, ii) junção espacial, iii) agregação espacial e iv) distância espacial. Apresentaremos essas operações utilizando as principais funções utilizando os dados de nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Filtro espacial
Filtros espaciais são operações que realizam seleção de feições espaciais entre dois objetos geoespaciais (x e y). Existe uma grande quantidade de funções para realizar filtros espaciais no R, como podemos ver na Tabela 15.12. Essas funções verificam se cada feição em x mantém sua relação em y. Ao especificar o parâmetro sparse = FALSE, as funções retornam uma matriz lógica (composta por TRUE e FALSE).
| Função | Descrição | Função inversa |
|---|---|---|
sf::st_contains()
|
Nenhum dos pontos de x está fora de y |
st_within
|
sf::st_contains_properly()
|
x contém y, e y não tem pontos em comum com a fronteira de x | NA |
sf::st_covers()
|
Nenhum ponto de y se encontra no exterior de x |
st_covered_by
|
sf::st_covered_by()
|
Inverso de sf::st_covers()
|
NA |
sf::st_crosses()
|
x e y têm alguns, mas não todos os pontos internos em comum | NA |
sf::st_disjoint()
|
x e y não têm pontos em comum |
st_intersects
|
sf::st_equals()
|
x e y são geometricamente iguais; o número de pedido dos nós pode ser diferente | NA |
sf::st_equals_exact()
|
x e y são geometricamente iguais e têm ordem de nó idêntica | NA |
sf::st_intersects()
|
x e y não são separados |
st_disjoint
|
sf::st_is_within_distance()
|
x está mais perto de y do que uma determinada distância | NA |
sf::st_within()
|
Nenhum dos pontos de y está fora de x |
st_contains
|
sf::st_touches()
|
x e y têm pelo menos um ponto limite em comum, mas nenhum ponto interno | NA |
sf::st_overlaps()
|
x e y têm alguns pontos em comum; a dimensão destes é idêntica à de x e y | NA |
sf::st_relate()
|
Dado um padrão, retorna se x e y aderem a este padrão | NA |
Em nosso exemplo, utilizaremos a função sf::intersects() para filtrar as nascentes dentro de floresta para Rio Claro/SP. Essa função vai retornar a resposta binária se as nascentes estão (1) ou não (empty) dentro dos polígonos de floresta.
## Filtro espacial
sf::st_intersects(x = geo_vetor_nascentes, y = geo_vetor_cobertura_floresta)
#> Sparse geometry binary predicate list of length 1220, where the predicate was `intersects'
#> first 10 elements:
#> 1: 1
#> 2: 1
#> 3: (empty)
#> 4: 1
#> 5: (empty)
#> 6: (empty)
#> 7: (empty)
#> 8: (empty)
#> 9: 1
#> 10: (empty)Podemos usar essa mesma função em conjunto com a função dplyr::filter() para filtrar as nascentes dentro de florestas, mas agora com o argumento sparse = FALSE para valores lógicos funcionarem com o filtro.
## Filtro espacial - interno
geo_vetor_nascentes_floresta_int <- geo_vetor_nascentes %>%
dplyr::filter(sf::st_intersects(x = ., y = geo_vetor_cobertura_floresta, sparse = FALSE))Ou ainda podemos utilizar o operador [] para realizar esse filtro, como podemos notar na Figura 15.26.
## Filtro espacial com [] - interno
geo_vetor_nascentes_floresta_int <- geo_vetor_nascentes[geo_vetor_cobertura_floresta, ]
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_cobertura_floresta$geometry, col = "forestgreen", add = TRUE)
plot(geo_vetor_nascentes_floresta_int$geometry, col = "blue", pch = 20, cex = 1, add = TRUE)
Figura 15.26: Nascentes dentro de florestas no município de Rio Claro/SP.
Entretanto, muitas vezes queremos fazer o filtro de feições que estão fora de feições de outro objeto espacial. Para isso, podemos usar a função sf::st_disjoint() ou ainda utilizando o operador [], mas com o argumento op, nesse caso utilizando a mesma função sf::st_disjoint() como operação (Figura 15.27). Atentar o segundo argumento vazio nessa operação.
## Filtro espacial - externo
geo_vetor_nascentes_floresta_ext <- geo_vetor_nascentes %>%
dplyr::filter(sf::st_disjoint(x = ., y = geo_vetor_cobertura_floresta, sparse = FALSE))
## Filtro espacial com [] - externo
geo_vetor_nascentes_floresta_ext <- geo_vetor_nascentes[geo_vetor_cobertura_floresta, , op = st_disjoint]
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_cobertura_floresta$geometry, col = "forestgreen", add = TRUE)
plot(geo_vetor_nascentes_floresta_ext$geometry, col = "steelblue", pch = 20, cex = 1, add = TRUE)
Figura 15.27: Nascentes fora de florestas no município de Rio Claro/SP.
Junção espacial
Outra operação muito usada dentro de análises espaciais é a junção espacial ou do inglês spatial join. A ideia base é muito semelhante com a junção baseada em atributos, mas aqui atribuiremos o valor da tabela de atributos das feições de um objeto espacial y às feições que fazem intersecção com um objeto espacial x, de modo que esses valores sejam armazenados na tabela de atributos do primeiro objeto espacial.
Para exemplificar, vamos atribuir os valores dos polígonos de cobertura da terra aos pontos de nascentes para Rio Claro/SP, fazendo um agrupamento pela tabela de atributos para permitir criar o mapa da Figura 15.28.
## Junção espacial
geo_vetor_nascentes_cob_jun <- geo_vetor_nascentes %>%
sf::st_join(x = ., y = geo_vetor_cobertura) %>%
dplyr::group_by(CLASSE_USO) %>%
dplyr::summarise(n = n())
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_nascentes_cob_jun[1], col = c("blue", "orange", "gray30", "forestgreen", "green"),
pch = 20, add = TRUE)
legend(x = 209000, y = 7520000, pch = 15, cex = .7, pt.cex = 2.5,
legend = (geo_vetor_nascentes_cob_jun$CLASSE_USO),
col = c("blue", "orange", "gray30", "forestgreen", "green"))
Figura 15.28: Junção espacial da cobertura da terra para as nascentes no município de Rio Claro/SP.
Agregação espacial
Muitas vezes queremos contabilizar quantas feições ou agregar valores de feições para polígonos. Podemos realizar essa operação usando as funções dplyr::group_by() e dplyr::summarise, ou utilizar a função aggregate(). Nesse exemplo, vamos contabilizar quantas nascentes há por polígono de cobertura da terra para o município de Rio Claro/SP (Figura 15.29).
## Agregação espacial
geo_vetor_cobertura_nas_agre <- geo_vetor_nascentes %>%
aggregate(x = ., by = geo_vetor_cobertura, FUN = length)
## Plot
plot(geo_vetor_cobertura_nas_agre[1], axes = TRUE, graticule = TRUE, main = NA)
Figura 15.29: Agregação espacial contabilizando o número de nascentes para cada classe de cobertura da terra no município de Rio Claro/SP.
Distância espacial
A distância espacial é a distância calculada em duas dimensões (2D) entre um objeto espacial x e y baseado no CRS e para cada feição dos objetos geoespaciais. Para realizar esse cálculo, utilizamos a função sf::st_distance(). Em nosso exemplo, vamos calcular a distância das nascentes até a floresta mais próxima, e adicionar essa informação para cada ponto na tabela de atributos com a função dplyr::mutate(), para o município de Rio Claro/SP (Figura 15.30).
## Distância espacial
geo_vetor_nascentes_dist_flo <- geo_vetor_nascentes %>%
dplyr::mutate(dist_flo = sf::st_distance(geo_vetor_nascentes, geo_vetor_cobertura_floresta))
## Plot
plot(geo_vetor_nascentes_dist_flo[7], pch = 20, axes = TRUE, graticule = TRUE, main = NA)
Figura 15.30: Distância espacial das nascentes até o fragmento de floresta mais próxima no município de Rio Claro/SP.
Raster
As principais operações espaciais para dados raster podem ser classificas, segundo Lovelace et al. (2019), em: i) operações locais (por célula), ii) operações focais (por bloco de múltiplas células regulares - e.g. 3x3), iii) operações zonais (por bloco de múltiplas células irregulares) e iv) operações globais (por um ou vários rasters inteiros). Cada uma delas é aplicada para objetivos e escalas espaciais específicas. Para os exemplos desta seção, utilizaremos o dado raster de elevação para o município de Rio Claro/SP.
Operações locais
As operações locais contemplam todas as operações realizadas célula a célula em uma ou várias camadas de um objeto raster. A álgebra de raster é uma das mais comuns, simples e poderosas operações no R envolvendo rasters. Com ela podemos fazer operações simples através de operadores aritméticos (soma, subtração, multiplicação, divisão ou potenciação) entre dois ou mais objetos raster, ou utilizar funções para alterar todos os valores dos pixels como, por exemplo, as funções log10() ou sqrt(), ou ainda a função raster::scale() para padronizar ou centralizar os valores dos rasters (Figura 15.31).
## Soma
geo_raster_srtm_rio_claro2 <- geo_raster_srtm_rio_claro + geo_raster_srtm_rio_claro
## Log10
geo_raster_srtm_rio_claro_log10 <- log10(geo_raster_srtm_rio_claro)
## Plot
old_par <- par(mfrow = c(1, 2))
plot(geo_raster_srtm_rio_claro2, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
plot(geo_raster_srtm_rio_claro_log10, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
par(old_par)
Figura 15.31: Rasters de soma e log10 do mapa de elevação para Rio Claro/SP.
Além das operações aritméticas, a álgebra de rasters também permite operações lógicas, como criar um raster binário (composto por 1 quando a operação lógica é verdadeira, e 0 quanto é falsa). Em nosso caso, buscamos todos os pixels acima de 600 metros para o raster de elevação de Rio Claro/SP (Figura 15.32).
## Acima de 600
geo_raster_srtm_rio_claro_acima_600 <- geo_raster_srtm_rio_claro > 600
## Plot
plot(geo_raster_srtm_rio_claro_acima_600, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.32: Operação local lógica mostrando todos os pixels acima de 600 metros de elevação para Rio Claro/SP.
Além dos operadores aritméticos, também podemos usar as funções raster::calc() (uma camada) e raster::overlay() (duas ou mais camadas) para realizar operações em todas as células. Elas funcionam com a criação de uma função específica através da função function() (Capítulo 4), para que esta seja aplicada em todas as células do raster. Essas funções são muito eficientes, portanto, são preferíveis para grandes conjuntos de dados raster. Exemplificaremos essa operação calculando o produto de todos os pixels por eles mesmos do raster de elevação de Rio Claro/SP (Figura 15.33).
## Produto dos pixel - calc
geo_raster_srtm_rio_claro_prod <- raster::calc(x = geo_raster_srtm_rio_claro, fun = function(x){x * x})
geo_raster_srtm_rio_claro_prod
#> class : RasterLayer
#> dimensions : 370, 364, 134680 (nrow, ncol, ncell)
#> resolution : 0.0008333333, 0.0008333333 (x, y)
#> extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : layer
#> values : 241081, 970225 (min, max)
## Plot
plot(geo_raster_srtm_rio_claro_prod, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.33: Operação local de multiplicação de todos os pixels por eles mesmos do raster de elevação para Rio Claro/SP.
A predição de objetos raster (utilizando a função raster::predict()) é outra aplicação extremamente útil em operações locais, nós a veremos mais à frente neste capítulo. Essa função possui basicamente dois argumentos: object que é os rasters preditores e model com o modelo ajustado para o qual os valores serão preditos com base nos valores dos rasters. A partir da relação entre variáveis respostas (e.g, pontos no espaço, como ocorrência ou riqueza de espécies), e variáveis preditoras (rasters contínuos de elevação, pH, precipitação, temperatura, cobertura da terra ou classe de solo), criamos modelos usando funções como lm(), glm(), gam() ou uma técnica de aprendizado de máquina, e fazemos predições espaciais aplicando os coeficientes estimados aos valores dos raster preditores (consulte od Capítulos 7 e 8).
Por fim, a reclassificação de rasters é outra operação muito comum quando trabalhamos com rasters. Nela é realizada a classificação de intervalos de valores numéricos em grupos, e.g. agrupar um modelo digital de elevação em classes de valores. A função que faz essa operação é a raster::reclassify(). Ela possui dois argumentos: x que é o raster a ser reclassificado, e o segundo rcl para o qual devemos construir uma matriz de reclassificação, onde a primeira coluna é a extremidade inferior, a segunda coluna é a extremidade superior, e a terceira coluna representa o novo valor para os intervalos das colunas um e dois. Vamos reclassificar o raster de elevação de Rio Claro/SP para os intervalos 400–600, 600–800 e 800–1000 que são reclassificados para os valores 1, 2 e 3, respectivamente (Figura 15.34).
## Matriz de reclassificação
rcl <- matrix(c(400,600,1,
600,800,2,
800,1000,3),
ncol = 3, byrow = TRUE)
## Reclassifição
geo_raster_srtm_rio_claro_rcl <- raster::reclassify(x = geo_raster_srtm_rio_claro, rcl = rcl)
## Plot
plot(geo_raster_srtm_rio_claro_rcl, col = viridis::viridis(3))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.34: Operação local de reclassificação para três classes de elevação para Rio Claro/SP.
Operações focais
As operações focais levam em consideração uma célula central e seus vizinhos. A vizinhança (também chamada de janela móvel - moving window) tipicamente é composta de células de 3 por 3 (célula central e seus oito vizinhos), mas pode assumir outra forma. A operação focal aplica uma função de agregação a todas as células dentro da vizinhança especificada, e usa a saída correspondente como o novo valor para a célula central, e segue para a próxima célula central e seus vizinhos. Essa operação é realizada através da função raster::focal(). O parâmetro x especifica o raster de entrada, o parâmetro w define a janela móvel por uma matriz cujos valores correspondem a pesos, e por fim, o parâmetro fun especifica a função que desejamos aplicar às células, como min(), max(), sum(), mean(), sd() ou var(). Existem diversas aplicações dessa operação para dados raster, como no processamento de imagens de satélite (ver mais em Wegmann et al. (2016)). Outra utilidade é para o cálculo de características topográficas, como declividade, aspecto e direções de fluxo. Para calcular essas métricas específicas, podemos utilizar a função raster::terrain().
Para nosso exemplo, vamos realizar o cálculo do desvio padrão da elevação e a métrica de aspecto (orientação da vertente) para o raster de elevação em Rio Claro/SP (Figura 15.35).
## Janela móvel - moving window
geo_raster_srtm_rio_claro_focal_sd <- raster::focal(
x = geo_raster_srtm_rio_claro,
w = matrix(data = 1, nrow = 3, ncol = 3),
fun = sd)
## Declividade
geo_raster_srtm_rio_claro_asp <- raster::terrain(x = geo_raster_srtm_rio_claro, opt = "aspect")
## Plot
old_par <- par(mfrow = c(1, 2))
plot(geo_raster_srtm_rio_claro_focal_sd, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
plot(geo_raster_srtm_rio_claro_asp, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
par(old_par)
Figura 15.35: Cálculo do desvio padrão da elevação para uma janela de 3x3 e do aspecto para Rio Claro/SP.
Operações zonais
As operações zonais aplicam uma função de agregação para várias células de um raster. Geralmente usa-se um segundo raster categórico para definir as zonas, de modo que as células raster que definem a zona não precisam ser vizinhas, como na operação focal. O resultado de uma operação zonal é uma tabela de resumo agrupada por zona, explicando porque essa operação também é conhecida como estatística zonal. Isso é um contraste com as operações focais que retornam um objeto raster.
A operação zonal é realizada através da função raster::zonal(), que recebe de entrada no x o raster contínuo, em z o raster categórico, e em fun a função que resumirá as células. Em nosso exemplo, vamos calcular diversas medidas resumo da elevação com a função summary() para cada classe de elevação que criamos anteriormente.
## Estatística zonal
geo_raster_srtm_rio_claro_zonal <- data.frame(raster::zonal(geo_raster_srtm_rio_claro, geo_raster_srtm_rio_claro_rcl, fun = "summary"))
colnames(geo_raster_srtm_rio_claro_zonal) <- c("zona", "min", "1qt", "mediana", "media", "3qt", "max")
geo_raster_srtm_rio_claro_zonal
#> zona min 1qt mediana media 3qt max
#> 1 1 491 552 574.0 567.5995 589 600
#> 2 2 601 620 640.0 650.6829 670 800
#> 3 3 801 817 832.5 834.2732 846 985Operações globais
As operações globais usam todo o conjunto de dados raster representando uma única zona. As operações globais mais comuns são estatísticas descritivas para todos os pixels do raster, utilizando a função raster::cellStats() ou raster::freq(), já vistas. Além das estatísticas descritivas, podemos gerar rasters de distância, que calcula a distância de cada célula a uma ou um grupo células-alvo específica, utilizando a função raster::distance().
Em nosso exemplo, vamos selecionar os pixels abaixo de 500 m do raster de elevação e calcular a Distância Euclidiana (Figura 15.36).
## Distância euclidiana
geo_raster_srtm_rio_claro_abaixo_500 <- raster::calc(
x = geo_raster_srtm_rio_claro,
fun = function(x) ifelse(x < 500, 1, NA))
geo_raster_srtm_rio_claro_global_dist <- raster::distance(geo_raster_srtm_rio_claro_abaixo_500)
## Plot
plot(geo_raster_srtm_rio_claro_global_dist, col = viridis::viridis(10))
plot(geo_raster_srtm_rio_claro_abaixo_500, add = TRUE, col = "white", legend = FALSE)
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.36: Raster de distância Euclidiana dos pixels abaixo de 500 m de elevação para Rio Claro/SP.
15.9.3 Operações geométricas
As operações geométricas realizam modificações em objetos geoespaciais baseado na geometria do vetor ou do raster e na interação e conversão entre vetor-raster. As operações geométricas vetoriais podem ser unárias (funcionam em uma única geometria) ou binárias (modificam uma geometria com base na forma de outra geometria). Ainda podemos fazer transformações para alterar os tipos vetores, que refletirá se as feições são únicas ou múltiplas, inclusive na tabela de atributos. As operações geométricas em rasters envolvem mudar a posição, tamanho e número dos pixels subjacentes e atribuir-lhes novos valores. Por fim, podemos ainda fazer operações de interações e conversões entre raster-vetor para ajustar rasters a vetores, assim como converter um objeto espacial vetorial para raster e vice-versa.
Vetor
Como dissemos, as operações geométricas em vetores criam ou alteraram a geometria de objetos da classe sf, podendo fazer alterações em única geometria (unárias): i) simplificação, ii) centroides, iii) pontos aleatórios, iv) buffers, v) polígono convexo, vi) polígonos de Voronoi, vii) quadrículas e hexágonos; ou modificar uma geometria com base na forma de outra geometria (binárias), viii) união e ix) recortes; ou ainda fazer transformações de tipo de geometrias.
Para exemplificar as operações geométricas com vetores, vamos utilizar os dados do limite, nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Simplificação
A simplificação possui o intuito de generalizar linhas ou polígonos, diminuindo assim suas complexidades em relação ao número de vértices. É utilizada para representação em mapas menores ou mapas interativos ou ainda quando um objeto vetorial é muito grande. A função utilizada é a sf::st_simplify(), que usa o argumento x para uma geometria de entrada e dTolerance para controlar o nível de generalização nas unidades do mapa. Em nosso exemplo, simplificaremos a hidrografia de Rio Claro/SP (Figura 15.37).
## Simplificação
geo_vetor_hidrografia_simplificado <- sf::st_simplify(x = geo_vetor_hidrografia, dTolerance = 1000)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_hidrografia$geometry, col = "steelblue", lwd = 2, add = TRUE)
plot(geo_vetor_hidrografia_simplificado$geometry, col = adjustcolor("black", .7), add = TRUE)
Figura 15.37: Simplificação da hidrografia para Rio Claro/SP. As linhas em azul é a hidrografia original e as linhas em preto mostram a simplificação.
Centroides
A operação de centroides identifica o centro de objetos geoespaciais, geralmente o centro de massa das feições. É utilizado para gerar um ponto simples para representações complexas ou para estimar a distância entre polígonos utilizando esse centroide. Podemos calculá-los com a função sf::st_centroids() ou com a função sf::st_point_on_surface() para garantir que esses centroides caiam dentro das geometrias. Aqui calcularemos o centroide do município de Rio Claro/SP (Figura 15.38).
## Centroides
geo_vetor_rio_claro_sirgas2000_utm23s_cent <- sf::st_centroid(geo_vetor_rio_claro_sirgas2000_utm23s)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_cent$geom, cex = 3, pch = 20, add = TRUE)
Figura 15.38: Centroide do limite do município de Rio Claro/SP.
Pontos aleatórios
Por vezes precisamos criar algum padrão aleatório dentro de um contexto espacial. Isso pode ser realizado de diversas formas. Uma delas é a criação de pontos aleatórios dentro de um polígono. Podemos realizar essa operação com a função sf::st_sample(). Para essa função, dois argumentos são utilizados: x uma geometria de entrada e o size indicando o número de pontos à seres criados. Outro argumento bastante interessante é o type, indicando o tipo de amostragem espacial (aleatório, regular ou triangular). Para nosso exemplo, vamos fixar a amostragem utilizando a função set.seed() e sortear 30 pontos para o limite do município de Rio Claro/SP (Figura 15.39). Para mais detalhes da função set.seed(), consultar o Capítulo 4.
## Fixar amostragem
set.seed(42)
## Pontos aleatórios
geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios <- sf::st_sample(geo_vetor_rio_claro_sirgas2000_utm23s, size = 30)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, pch = 20, add = TRUE)
Figura 15.39: Sorteio de 30 pontos aleatório para Rio Claro/SP.
Buffer
Buffers são polígonos que representam a área dentro de uma determinada distância de um elemento geométrico, independentemente de ser um ponto, linha ou polígono. O buffer é comumente utilizado para análise de dados geoespaciais, geralmente sendo entendio como uma unidade amostral, delimitando uma porção no entorno de algum elemento ou evento, como as condições climáticas ou da estrutura da paisagem para uma amostragem, ou as características de cobertura da terra ao longo de um corpo d’água, e.g., Áreas de Preservação Permanente (APPs).
A função utilizada para criar buffers é a sf::st_buffer(), que requer pelo menos dois argumentos: x uma geometria de entrada e o dist uma distância para o buffer, fornecido nas unidades do CRS da geometria de entrada. Em nosso exemplo, vamos criar buffers circulares de 1000 metros para os 30 pontos aleatórios criados anteriormente para o município de Rio Claro/SP (Figura 15.40).
## Buffer
geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer <- sf::st_buffer(
x = geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, dist = 1000)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer, col = NA, lwd = 2, border = "red", add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Figura 15.40: Buffers de 1000 metros para os 30 pontos aleatórios no município de Rio Claro/SP.
Podemos ainda criar buffers quadrados acrescentando o argumento uendCapStyle = "SQUARE" (Figura 15.41).
## Buffer
geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_quad <- sf::st_buffer(
x = geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, dist = 1000, endCapStyle = "SQUARE")
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_quad, col = NA, lwd = 2, border = "red", add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Figura 15.41: Buffers quadrados de 1000 metros para os 30 pontos aleatórios no município de Rio Claro/SP.
Polígono convexo
Uma análise bastante comum, principalmente realizada pela IUCN, é a criação de polígonos convexos, para definir a extensão de ocorrência de uma espécie (Extent of occurrence - EOO). Nesse sentido, essa operação liga os pontos externos de um conjunto de pontos e cria um polígono a partir deles. Podemos criar esse polígono com a função sf::st_convex_hull(). Um único passo que precisamos adiantar é utilizar a função sf::st_union() para unir todos os pontos e criar um objeto sf MULTIPOINT, já iremos explicar com mais detalhes adiante. Vamos utilizar os pontos aleatórios que criamos anteriormente para criar o polígono convexo (Figura 15.42).
## Polígono convexo
geo_vetor_rio_claro_sirgas2000_utm23s_convexo <- geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios %>%
sf::st_union() %>%
sf::st_convex_hull()
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_convexo, col = NA, lwd = 2, border = "red", add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Figura 15.42: Polígono convexo para os 30 pontos criados aleatoriamente para Rio Claro/SP.
Polígonos de Voronoi
Uma outra forma de criar polígonos para resumir dados espaciais é através dos Polígonos de Voronoi ou Diagrama de Voronoi. Nele, polígonos irregulares são criados a partir da proximidade de pontos, de modo a estimar uma área de abrangência no entorno dos mesmos (Okabe 2000). Esses polígonos podem ser criados com a função sf::st_voronoi(), mas precisamos novamente utilizar a função sf::st_union() para unir todos os pontos e criar um objeto sf do tipo MULTIPOINT. Vamos utilizar os pontos aleatórios que criamos anteriormente para criar o polígono de Voronoi (Figura 15.43).
## Polígonos de Voronoi
geo_vetor_rio_claro_sirgas2000_utm23s_voronoi <- geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios %>%
sf::st_union() %>%
sf::st_voronoi()
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_voronoi, col = NA, lwd = 2, border = "red", add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Figura 15.43: Polígonos de Voronoi para os 30 pontos criados aleatoriamente para Rio Claro/SP.
Quadrículas e hexágonos
Muitas vezes precisamos criar unidades espaciais idênticas e igualmente espaçadas para resumir informações dispersas por toda a nossa área de estudo. Uma prática muito comum é a criação de um gride de pontos ou quadrículas em toda a área de estudo, e depois utilizar essas geometrias para associar ou resumir informações espacializadas, como a IUCN utiliza para a análise de área de ocupação (Are of occupancy - AOO). Além das quadrículas, uma outra geometria que se tornou bastante comum para as finalidades descritas, é a criação de hexágonos, que além de serem mais esteticamente atraentes, possuem uma explicação matemática de sua melhor funcionalidade para análises espaciais em Ecologia (Birch, Oom, and Beecham 2007).
A função utilizada para criar esses grides é a sf::st_make_grid(), que requer pelo menos dois argumentos: x uma geometria de entrada e o cellsize indicando o tamanho do gride a ser criado, fornecido nas unidades do CRS da geometria de entrada. Há diversos outros argumentos, mas os mais importantes são o square que definirá se o gride será de quadriculas ou de hexágonos, e o what que definirá se geraremos polígonos, cantos ou centroides.
Em nosso exemplo, vamos criar quadrículas e hexágonos de 2000 metros (i.e. 4000000 metros quadrados) para o município de Rio Claro/SP (Figura 15.44 e Figura 15.45). Podemos ainda utilizar as funções de filtros espaciais (Tabela 15.12) para definir como selecionaremos esses elementos para a área de estudo. Aqui utilizamos a função sf::st_intersects().
## Quadrículas
geo_vetor_rio_claro_sirgas2000_utm23s_grid <- sf::st_make_grid(
x = geo_vetor_rio_claro_sirgas2000_utm23s, cellsize = 2000, what = "polygons") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = geo_vetor_rio_claro_sirgas2000_utm23s, sparse = FALSE))
## Centroides das quadrículas
geo_vetor_rio_claro_sirgas2000_utm23s_grid_cent <- geo_vetor_rio_claro_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, what = "centers") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = sf::st_union(geo_vetor_rio_claro_sirgas2000_utm23s_grid), sparse = FALSE))
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_grid, col = NA, border = "red", lwd = 2, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_grid_cent, pch = 20, add = TRUE)
Figura 15.44: Quadrículas de 2000 metros de arestas e centroides para Rio Claro/SP.
## Hexágonos
geo_vetor_rio_claro_sirgas2000_utm23s_hex <- geo_vetor_rio_claro_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = .,
y = geo_vetor_rio_claro_sirgas2000_utm23s,
sparse = FALSE))
## Centroides de hexágonos
geo_vetor_rio_claro_sirgas2000_utm23s_hex_cent <- geo_vetor_rio_claro_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, square = FALSE, what = "centers") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = .,
y = sf::st_union(geo_vetor_rio_claro_sirgas2000_utm23s_hex),
sparse = FALSE))
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_hex, col = NA, border = "red", lwd = 2, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_hex_cent, pch = 20, add = TRUE)
Figura 15.45: Hexágonos equivalentes a quadrículas de 2000 metros de arestas e centroides para Rio Claro/SP.
União (“dissolver”)
Como vimos, na agregação por atributos podemos dissolver as geometrias de polígonos do mesmo grupo pelos valores da tabela de atributos, onde, naquele exemplo, contabilizamos quantas nascentes haviam por polígono de cobertura da terra para o município de Rio Claro/SP (Figura 15.29).
Nesta seção, vamos utilizar a função sf::st_union() para unir diversas feições em uma só, dissolvendo os limites entre elas. Vamos utilizar de exemplo os buffers que criamos a partir dos 30 pontos aleatórios (Figura 15.46).
## União
geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao <- sf::st_union(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
Figura 15.46: União - dissolução - dos buffers criados a partir dos 30 pontos aleatórios para Rio Claro/SP.
Recorte e apagar (“clip” e “erase”)
O recorte realiza um subconjunto espacial envolvendo dois objetos geoespaciais. O recorte é aplicado somente a linhas e polígonos, ou seja, usaremos linhas e polígonos para recortar linhas ou polígonos. Esse recorte pode ser realizado de três formas: i) intersecção (subconjunto das geometrias sobrepostas entre os dois objetos), ii) diferença (subconjunto das geometrias do primeiro objeto sem sobreposição com o segundo objeto), e iii) diferença simétrica (apenas as geometrias não sobrepostas entre os dois objetos). Respectivamente para cada uma dessas operações temos funções específicas: sf::st_intersection(), sf::st_difference() e sf::st_sym_difference().
Para nosso exemplo, faremos o recorte da hidrografia em relação aos buffers criados e unidos para os 30 pontos aleatórios em Rio Claro/SP. Primeiramente, faremos o recorte para dentro dos buffers com a função sf::st_intersection() (Figura 15.47).
## Recorte - intersecção
geo_vetor_hidrografia_interseccao <- sf::st_intersection(x = geo_vetor_hidrografia, y = geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
plot(geo_vetor_hidrografia_interseccao$geometry, col = "blue", add = TRUE)
Figura 15.47: Recorte da hidrografia para fora dos buffers dos 30 pontos aleatórios para Rio Claro/SP.
Para nosso segundo exemplo, realizamos o recorte da hidrografia em relação aos buffers, mas agora para fora dos buffers utilizando a função sf::st_difference() (Figura 15.48), que seria semelhando a operação de apagar (“erase”).
## Recorte - diferença
geo_vetor_hidrografia_diferenca <- sf::st_difference(x = geo_vetor_hidrografia, y = geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao)
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
plot(geo_vetor_hidrografia_diferenca$geometry, col = "blue", add = TRUE)
Figura 15.48: Recorte da hidrografia para fora dos buffers dos 30 aleatórios para Rio Claro/SP.
Transformações de tipo
Esse tópico possui muitas funcionalidades, que são exploradas no tópico “5.2.7 Type transformations” de Lovelace et al. (2019). Aqui, nosso interesse principal é em relação à transformação dos tipos de objetos geoespaciais da classe sf: MULTIPOINT, MULTILINESTRING e MULTIPOLYGON, para POINT, LINESTRING e POLYGON. Muitas vezes as feições de nossos objetos, i.e., as linhas da tabela de atributos estão agrupadas em apenas uma linha da tabela. Quando o objeto espacial está nesse formato, geralmente em alguma classe dessas (MULTIPOINT, MULTILINESTRING e MULTIPOLYGON), não temos como realizar operações espaciais ou geométricas para cada feição, e precisamos separá-las em linhas diferentes para que operações como o cálculo de comprimento ou área seja possível para cada feição.
Dessa forma, podemos utilizar a função sf::st_cast() para fazer essas transformações e atribuir cada feição a uma linha da tabela de atributos. Como exemplo, vamos separar os fragmentos de floresta e calcular a área para cada feição em hectares (Figura 15.49).
## Transformação de tipo
geo_vetor_cobertura_floresta_polygon <- geo_vetor_cobertura_floresta %>%
sf::st_cast("POLYGON") %>%
dplyr::mutate(area_ha = sf::st_area(.)/1e4 %>% round(2))
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(geo_vetor_cobertura_floresta_polygon["area_ha"], col = viridis::viridis(100), add = TRUE)
Figura 15.49: Transformação do vetor de florestas em POLYGON e cálculo da área para cada feição para Rio Claro/SP.
Raster
As operações geométricas em rasters envolvem mudar a posição, tamanho e número dos pixels e atribuir novos valores, geralmente aumentando ou diminuindo o tamanho desses pixels. Essas operações permitem alinhar rasters de diversas fontes, fazendo com que compartilhem uma correspondência entre seus pixels, permitindo que eles sejam processados todos juntos, ou simplesmente permite a realização de análises que demorariam muito, caso os rasters possuam um tamanho de pixel muito pequeno.
📝 Importante
Essas operações funcionam para as três classes dos objetos raster: RasterLayer, RasterBrick e RasterStack.
Para exemplificar as operações geométricas com rasters, vamos utilizar os dados de elevação para o município de Rio Claro/SP e bioclimáticos para o mundo.
Agregação
Na agregação de rasters, aumentamos o tamanho dos pixels (diminuindo a resolução), agregando os valores dos pixels em um pixel maior. Podemos realizar essa operação com a função raster::aggregate(), que possui três argumentos: x corresponde ao objeto raster de entrada, fact é o fator de agregação e corresponde ao número que definirá o novo tamanho do pixel (e.g., se um raster tem resolução de 90 m, um fator de agregação de 10 fará com o novo raster tenha a resolução de 900 m), e fun é a função utilizada para realizar a agregação dos pixels (Figura 15.50).
Em nosso exemplo, vamos aumentar o tamanho dos pixels para 900 metros do raster de elevação para Rio Claro/SP.
## Agregação - aumentar o tamanho do pixel
geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media <- raster::aggregate(x = geo_raster_srtm_rio_claro_sirgas2000_utm23s, fact = 10, fun = "mean")
geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media
#> class : RasterLayer
#> dimensions : 40, 37, 1480 (nrow, ncol, ncell)
#> resolution : 900, 900 (x, y)
#> extent : 214554.4, 247854.4, 7502625, 7538625 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
#> source : memory
#> names : srtm_27_17
#> values : 506.0024, 922.8709 (min, max)
## Plot
plot(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media, col = viridis::viridis(10))
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.50: Agregação (aumento do pixel para 900 metros) utilizando a média para o raster de elevação para Rio Claro/SP.
Desagregação
De modo contrário, na desagregação de rasters, diminuímos o tamanho dos pixels (aumentando a resolução), preenchendo com novos valores. Podemos realizar essa operação com a função raster::desaggregate(), que assim como a função anterior, possui três argumentos: x corresponde ao objeto raster de entrada, fact é o fator de desagregação e corresponde ao número que definirá o novo tamanho do pixel (e.g., se um raster tem resolução de 90 m, um fator de desagregação de 2 fará com que o novo raster tenha a resolução de 45 m), e method é a função utilizada para realizar a desagregação dos pixels (Figura 15.51).
Nesse exemplo, vamos diminuir o tamanho dos pixels para 45 metros do raster de elevação para Rio Claro/SP.
## Desagregação - diminuir o tamanho do pixel
geo_raster_srtm_rio_claro_desg_bil <- raster::disaggregate(x = geo_raster_srtm_rio_claro_sirgas2000_utm23s, fact = 2, method = "bilinear")
geo_raster_srtm_rio_claro_desg_bil
#> class : RasterLayer
#> dimensions : 792, 728, 576576 (nrow, ncol, ncell)
#> resolution : 45, 45 (x, y)
#> extent : 214554.4, 247314.4, 7502985, 7538625 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
#> source : memory
#> names : srtm_27_17
#> values : 493.7046, 986.5187 (min, max)
## Plot
plot(geo_raster_srtm_rio_claro_desg_bil, col = viridis::viridis(10))
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.51: Desagregação (diminuição do pixel para 45 metros) utilizando o método bilinear para o raster de elevação para Rio Claro/SP.
Alinhamento de rasters
Muitas vezes queremos ir além de ajustar o tamanho do pixel, ajustando também a extensão, número e origem dos pixels para várias camadas rasters, principalmente se precisamos criar objetos das classes RasterBrick ou RasterStack. Dessa forma, podemos utilizar a função raster::compareRaster() para comparar os rasters em relação a extensão, número de linhas e colunas, projeção, resolução e origem (ou um subconjunto dessas comparações).
Podemos utilizar a função raster::resample() para fazer esse alinhamento, ou ainda a função gdalUtils::align_rasters(). Para a primeira função, os argumentos são x para o raster de entrada, y para o raster de alinhamento e method para o método utilizado no alinhamento. Para nosso exemplo, vamos ajustar uma camada bioclimática (BIO01) à camada de elevação para Rio Claro/SP (Figura 15.52).
## Reamostragem
geo_raster_bioclim_rc <- raster::resample(x = geo_raster_bioclim$bio01,
y = geo_raster_srtm_rio_claro,
method = "bilinear")
geo_raster_bioclim_rc
#> class : RasterLayer
#> dimensions : 370, 364, 134680 (nrow, ncol, ncell)
#> resolution : 0.0008333333, 0.0008333333 (x, y)
#> extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : bio01
#> values : 19.8383, 20.60492 (min, max)
## Plot
plot(geo_raster_bioclim_rc$bio01, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.52: Reamostragem (alinhamento dos rasters) utilizando o método bilinear para alinhar o raster bioclimático (BIO01) ao de elevação para Rio Claro/SP.
Interações raster-vetor
Podemos fazer operações da interação entre objetos vetoriais e raster, como ajustes da extensão e limite do raster para vetores (corte e máscara), extração dos valores dos pixels para vetores (pontos, linhas e polígonos), e estatísticas zonais dos valores dos pixels dos raster para um vetor (linhas e polígonos).
Cortes e máscaras
Muitas vezes precisamos ajustar o tamanho de um objeto raster a uma área menor de interesse, geralmente definido por um objeto vetorial. Para realizar essa operação, dispomos de duas funções: raster::crop() e raster::mask().
📝 Importante
É fundamental que ambos os objetos raster a ser reduzido e o vetor como molde precisam estar no mesmo CRS.
A função raster::crop() ajusta o raster à extensão do vetor. Como exemplo, vamos retomar o raster de elevação original baixado e importado anteriormente (Figura 15.14). Primeiramente, vamos usar a função raster::crop() para ajustar esse raster à extensão do limite do município de Rio Claro/SP (Figura 15.53).
## Crop - adjuste da extensão
geo_raster_srtm_rio_claro_crop <- raster::crop(geo_raster_srtm, geo_vetor_rio_claro)
## Plot
plot(geo_raster_srtm_rio_claro_crop, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.53: Ajuste da extensão do raster de elevação para a extensão de Rio Claro/SP.
Para ajustar o raster ao limite do município de Rio Claro/SP, vamos usar a função raster::mask(). É importante notar que essa função preenche com NAs os pixels que estão fora do limite do polígono e não ajusta a extensão (Figura 15.53).
## Mask - adjuste ao limite
geo_raster_srtm_rio_claro_mask <- raster::mask(geo_raster_srtm, geo_vetor_rio_claro)
## Plot
plot(geo_raster_srtm_rio_claro_mask, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.54: Ajuste do raster de elevação para o limite de Rio Claro/SP.
Para ajustar o raster à extensão e ao limite do município de Rio Claro/SP, precisamos utilizar conjuntamente as funções raster::crop() e raster::mask() (Figura 15.55).
## Crop e mask - ajuste da extensão e do limite
geo_raster_srtm_rio_claro_crop_mask <- geo_raster_srtm %>%
raster::crop(geo_vetor_rio_claro) %>%
raster::mask(geo_vetor_rio_claro)
## Plot
plot(geo_raster_srtm_rio_claro_crop_mask, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.55: Ajuste da extensão e do limite do raster de elevação para Rio Claro/SP.
A função raster::mask() possui ainda um argumento chamado inverse, que cria uma máscara inversa ao limite, preenchendo com NA o pixels internos ao limite do polígono, como podemos ver para o raster de elevação e o limite de Rio Claro/SP (Figura 15.56).
## Crop e mask inversa - ajuste da extensão e do limite inverso
geo_raster_srtm_rio_claro_crop_mask_inv <- geo_raster_srtm %>%
raster::crop(geo_vetor_rio_claro) %>%
raster::mask(geo_vetor_rio_claro, inverse = TRUE)
## Plot
plot(geo_raster_srtm_rio_claro_crop_mask_inv, col = viridis::viridis(10))
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.56: Ajuste da extensão e do limite externo do raster de elevação para Rio Claro/SP.
Extração
A interação entre raster-vetor de extração é o processo que identifica e retorna valores associados de pixels de um raster com base em um objeto vetorial. É uma operação extremamente comum em análises espaciais, principalmente para associar valores de raster ambientais (contínuos ou categóricos) a pontos de ocorrência ou amostragem. Os valores retornados dependerão do tipo vetor (pontos, linhas ou polígonos) e de argumentos da função raster::extract() que alteram o funcionamento da extração.
Em nosso exemplo, vamos extrair os valores do raster de elevação para as nascentes do município de Rio Claro/SP (Figura 15.57).
## Extração
geo_vetor_nascentes_ele <- geo_vetor_nascentes %>%
dplyr::mutate(elev = raster::extract(x = geo_raster_srtm_rio_claro_sirgas2000_utm23s, y = .))
## Plot
plot(geo_vetor_nascentes_ele["elev"],
pch = 20, main = NA, axes = TRUE, graticule = TRUE)
Figura 15.57: Extração dos valores de elevação para as nascentes de Rio Claro/SP.
Além da extração dos valores totais, podemos resumir os valores dos pixels com a mesma operação de extração, utilizando ainda a função raster::extract(), mas utilizando uma outra função para resumir os valores dos pixels para um polígono, operação também denominada de estatística zonal (agora para vetores). Já vimos que ela pode ser realizada entre rasters na seção de operações zonais, mas aqui a realizaremos para rasters e vetores.
Para o exemplo, vamos calcular a elevação média dos valores para os hexágonos que criamos para o limite de Rio Claro/SP( Figura 15.58).
## Extração - estatística por zonas
geo_vetor_rio_claro_sirgas2000_utm23s_hex_alt <- geo_vetor_rio_claro_sirgas2000_utm23s_hex %>%
dplyr::mutate(elev_mean = raster::extract(x = geo_raster_srtm_rio_claro_sirgas2000_utm23s,
y = geo_vetor_rio_claro_sirgas2000_utm23s_hex,
fun = mean,
na.rm = TRUE))
## Plot
plot(geo_vetor_rio_claro_sirgas2000_utm23s_hex_alt["elev_mean"],
pch = 20, main = NA, axes = TRUE, graticule = TRUE)
Figura 15.58: Extração dos valores de elevação e resumo pela média para os hexágonos de Rio Claro/SP.
Conversões raster-vetor
Por fim, podemos ainda fazer operações de conversão entre objetos vetoriais para raster e vice-versa. Nessas operações, podemos resumir ou transformar objetos vetoriais (pontos, linhas ou polígonos) para rasters, escolhendo um raster previamente existente, processo denominado rasterização. Também podemos realizar o processo inverso, i.e., transformar o raster em um vetor, podendo esse vetor ser pontos, linhas ou polígonos, operação chamada de vetorização.
Rasterização
A conversão de vetor para raster pode ser realizada de pontos, linhas ou polígonos para rasters. Nesse processo, podemos utilizar uma função para resumir os dados pontuais para os pixels do raster que criaremos. Para essa operação, podemos utilizar a função raster::rasterize(), com o argumento x sendo o vetor de entrada, y o raster base, field a coluna ou campo da tabela de atributos do objeto vetorial para os quais os valores serão utilizados e fun a função utilizada para agregação dos dados.
Aqui, vamos contabilizar a quantidade de nascentes por pixel, utilizando como base o raster para o qual mudamos a resolução para 900 metros (Figura 15.59).
## Rasterizar pontos
geo_vetor_nascentes_rasterizacao <- raster::rasterize(x = geo_vetor_nascentes, y = geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media, field = 1, fun = "count")
## Plot
plot(geo_vetor_nascentes_rasterizacao, col = viridis::viridis(10))
plot(geo_vetor_nascentes$geometry, pch = 20, cex = .5, col = adjustcolor("gray", .5), add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.59: Rasterização das nascentes, com a operação de contabilização de pontos para Rio Claro/SP.
Além de pontos, podemos também rasterizar linhas. Aqui vamos contabilizar as linhas da hidrografia simplificada para Rio Claro/SP (Figura 15.60).
## Rasterizar linhas
geo_vetor_hidrografia_rasterizacao <- raster::rasterize(
x = geo_vetor_hidrografia_simplificado,
y = geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media,
field = 1, fun = "count")
## Plot
plot(geo_vetor_hidrografia_rasterizacao, col = viridis::viridis(10))
plot(geo_vetor_hidrografia_simplificado$geom, col = "gray", add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.60: Rasterização da hidrografia, com a operação de contabilização de linhas para Rio Claro/SP.
Podemos ainda rasterizar polígonos, de modo que cada pixel do raster a ser criado receberá o valor da tabela de atributos, ou uma análise pelo vizinho mais próximo no caso de um campo categórico, como a cobertura da terra, que também vai depender da resolução do raster base e do tamanho da feição do polígono. Para nosso exemplo, antes de criar o raster vamos transforma a coluna de classe de cobertura da terra em factor (Figura 15.61). Entretanto, essa operação de rasterização tende a demorar muito no caso de polígonos muito detalhados ou um raster com pixels muito pequenos, sendo que dois pacotes aceleram esse processamento (fasterize e gdalUtils), com suas respectivas funções: fasterize::fasterize() e gdalUtils::gdal_rasterize().
## Rasterizar polígonos
geo_vetor_cobertura_rasterizacao <- geo_vetor_cobertura %>%
dplyr::mutate(classe = as.factor(CLASSE_USO)) %>%
raster::rasterize(x = ., y = geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media, field = "classe")
## Plot
plot(geo_vetor_cobertura_rasterizacao, col = viridis::viridis(10))
plot(geo_vetor_cobertura$geom, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.61: Rasterização da cobertura da terra para Rio Claro/SP.
Vetorização
A operação inversa à rasterização é a vetorização, na qual convertemos um raster em um vetor, sendo que esse vetor receberá os valores dos pixels. O vetor em questão pode ser pontos (geralmente um gride de pontos), linhas (geralmente isolinhas ou linhas de contorno), ou polígonos (podendo esses polígonos ser ou não dissolvidos pelos valores dos pixels). Existem funções específicas para cada uma dessas conversões, sendo elas: raster::rasterToPoints(), raster::rasterToContour() e raster::rasterToPolygons(), respectivamente. Para a última função, ainda dispomos de uma alternativa mais veloz spex::polygonize().
Em nosso exemplo, vamos vetorizar o raster de elevação para Rio Claro/SP, criando um gride de pontos, sendo os pontos os centroides de cada pixels (Figura 15.62).
## Vetorização de pontos
geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media_pontos <- raster::rasterToPoints(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media, spatial = TRUE) %>%
sf::st_as_sf()
## Plot
plot(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media,
col = viridis::viridis(10, alpha = .8))
plot(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media_pontos,
pch = 20, cex = .7, main = FALSE, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom, col = NA,
border = "red", lwd = 2, add = TRUE)
Figura 15.62: Vetorização do raster de elevação criando um gride de pontos para Rio Claro/SP.
Neste outro exemplo, vamos vetorizar o raster de elevação para Rio Claro/SP novamente, mas agora criando isolinhas (Figura 15.63).
## Vetorização de linhas
geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media_linhas <- raster::rasterToContour(
x = geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media) %>%
sf::st_as_sf()
## Plot
plot(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media,
col = viridis::viridis(10, alpha = .8))
contour(geo_raster_srtm_rio_claro_sirgas2000_utm23s_agre_media,
labcex = 1, main = FALSE, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom,
col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.63: Vetorização do raster de elevação criando isolinhas para Rio Claro/SP.
Por fim, vamos vetorizar o raster de cobertura da terra criado anteriormente para Rio Claro/SP, criando polígonos não dissolvendo e dissolvidos (Figura 15.64).
## Vetorização de polígonos
geo_vetor_cobertura_rasterizacao_poligonos <- raster::rasterToPolygons(geo_vetor_cobertura_rasterizacao) %>%
sf::st_as_sf()
## Vetorização de polígonos dissolvendo
geo_vetor_cobertura_rasterizacao_poligonos_dissolvidos <- raster::rasterToPolygons(geo_vetor_cobertura_rasterizacao, dissolve = TRUE) %>%
sf::st_as_sf()
## Plot
old_par <- par(mfrow = c(1, 2))
plot(geo_vetor_cobertura_rasterizacao, col = viridis::viridis(10))
plot(geo_vetor_cobertura_rasterizacao_poligonos$geometry,
col = NA, border = "gray", lwd = 1, main = FALSE, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom,
col = NA, border = "red", lwd = 2, add = TRUE)
plot(geo_vetor_cobertura_rasterizacao, col = viridis::viridis(10))
plot(geo_vetor_cobertura_rasterizacao_poligonos_dissolvidos$geometry,
col = NA, border = "gray", lwd = 1, main = FALSE, add = TRUE)
plot(geo_vetor_rio_claro_sirgas2000_utm23s$geom,
col = NA, border = "red", lwd = 2, add = TRUE)
par(old_par)
Figura 15.64: Vetorização do raster de cobertura da terra para Rio Claro/SP, não dissolvendo e dissolvendos os polígonos gerados.
15.10 Visualização de dados geoespaciais
Um dos pontos finais de toda a análise envolvendo a manipulação de dados geoespaciais será a apresentação de um mapa com as informações de interesse geoespacializadas. Mas antes, é necessário ter conhecimento de alguns dos elementos principais para a composição de um mapa relativamente bem informativo. Além disso, o R nos permite criar tipos diferentes de mapas: estáticos, animados e interativos. Os mais comuns são os estáticos, mas podemos por vezes melhorar a apresentação dos dados geoespaciais criando mapas animados e/ou interativos, com o auxílio de páginas web. Por fim, veremos as melhores formas de exportar mapas para diferentes formatos.
15.10.1 Principais elementos de um mapa
Um mapa pode ser composto de vários elementos, tendo estes o intuito de auxiliar a visualização e entendimento de seu conteúdo. Apesar disso, nem todos os elementos necessitam estar presentes em todos os layouts de mapas, sendo que os mesmos devem atendem à necessidade das representações, podendo ser muitas vezes omitidos ou outros podem ser acrescentados.
Os principais elementos de um mapa geralmente são compostos por:
- Mapa principal (ocupando quase toda a área da figura)
- Mapa secundário (geralmente muito menor que o mapa principal e com o intuito de mostrar a localização do mapa principal num contexto mais amplo, como país ou continente)
- Título (para resumir o intuito do mapa)
- Legenda (apresentando as informações detalhadas das classes ou escala de valores, geralmente identificando as cores e/ou texturas),
- Barra de escala (representando a relação entre as unidades do mapa e o mundo real)
- Indicador de orientação (Norte) (indicando o norte geográfico, podendo ser representado por uma flecha, bússola ou compasso)
- Gride de coordenadas (coordenadas presentes nas laterais)
- Descrição do CRS (indicando qual o Sistema de Referência de Coordenadas)
- Informações de origem (informações sobre a fonte dos dados representados no mapa)
- Outros elementos auxiliares (como elementos textuais e figuras extras)
Podemos visualizar todos esses elementos resumidos na Figura 15.65.
Figura 15.65: Principais elementos de um mapa.
15.10.2 Principais pacotes para a composição de mapas
Há uma grande quantidade de pacotes para a composição de mapas no R. Aqui listamos os principais (Tabela 15.13).
| Pacote | Descrição |
|---|---|
ggplot2
|
Cria visualizações de dados elegantes usando a gramática de gráficos |
ggspatial
|
Estrutura de dados espaciais para ggplot2 |
ggmap
|
Visualização espacial com ggplot2 |
tmap
|
Mapas temáticos |
leaflet
|
Cria mapas da web interativos com a biblioteca JavaScript ‘Leaflet’ |
plotly
|
Cria gráficos interativos da Web por meio de ‘plotly.js’ |
cartography
|
Cartografia temática |
googleway
|
Cartografia temática |
mapview
|
Acessa APIs do Google Maps para recuperar dados e mapas de plotagem |
rasterVis
|
Visualização interativa de dados espaciais em R |
cartogram
|
Métodos de visualização para dados raster |
mapsf
|
Crie cartogramas com R |
geogrid
|
Transforme polígonos geoespaciais em grades regulares ou hexagonais |
geofacet
|
‘ggplot2’ Utilitários de facetação para dados geoespaciais |
globe
|
Mapas de visualização 2D e 3D da Terra, incluindo as linhas de costa |
linemap
|
Mapas de linhas |
15.10.3 Mapas estáticos
Mapas estáticos são mapas simples e fixos para visualização de dados, sendo o tipo mais comum de saída visual. No início da composição de mapas no R, esse era o único tipo de mapa que a linguagem permitia produzir, principalmente utilizando o pacote sp (E. J. Pebesma and Bivand 2005). No entanto, com o advento de ferramentas de visualização dinâmicas no R, como componentes HTML, os mapas puderam ser compostos de forma dinâmica (animados e interativos).
Neste tópico abordaremos funções simples para composição de mapas estáticos, como o plot(), além de pacotes para composição de mapas mais elaborados, como os pacotes ggplot2 (Wickham et al. 2020) e tmap (Tennekes 2021).
Função plot()
A função genérica plot() é a maneira mais rápida de compor mapas estáticos utilizando objetos geoespaciais vetoriais e raster, funcionando para ambos os pacotes que apresentamos anteriormente (sf e raster). Apesar da simplicidade, essa função geralmente tende a criar mapas com relativa velocidade, nos auxiliando principalmente em fases iniciais de desenvolvimento de um projeto. Essa função oferece dezenas de argumentos em R Base, permitindo alguns ajustes limitados, com resultados bastante interessantes.
Como dito anteriormente, a função plot() vai funcionar diferentemente dependendo da classe do objeto geoespacial. Para objetos geoespaciais sf, a função vai plotar um mapa para cada coluna da tabela de atributos. Vamos usar de exemplo nosso mapa de biomas mostrado com os principais elementos de um mapa, podendo inclusive selecionar apenas a coluna de características geoespaciais (geom).
Primeiramente, vamos fazer o download dos dados de limites de biomas, retirando os sistemas costeiros, usando o pacote geobr (Pereira and Goncalves 2021).
## Download de polígonos dos geo_vetor_biomas Brasileiros
geo_vetor_biomas <- geobr::read_biomes(showProgress = FALSE) %>%
dplyr::filter(name_biome != "Sistema Costeiro") %>%
dplyr::rename(nome_bioma = name_biome,
codigo_bioma = code_biome,
ano = year)Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_biomasAgora, quando utilizamos a função plot() para um objeto da classe sf, temos os três mapas, cada um indicando uma coluna da tabela de atritos (Figura 15.66).
## Plot
plot(geo_vetor_biomas)
Figura 15.66: Mapa feito com a função plot() de um objeto sf para os Biomas do Brasil.
Selecionando as colunas desse objeto, podemos escolher a informação que queremos plotar, por exemplo, apenas a geometria geom. Além disso, podemos acrescentar os argumentos col para colorir e main para o título, além dos argumentos axes e graticule para adicionar as coordenadas e quadrículas, respectivamente. A legenda pode ser adicionada com a função legend() (Figura 15.67).
## Plot
plot(geo_vetor_biomas$geom,
col = c("darkgreen", "orange", "orange4", "forestgreen",
"yellow", "yellow3"),
main = "Biomas do Brasil", axes = TRUE, graticule = TRUE)
legend(x = -75, y = -20, pch = 15, cex = .7, pt.cex = 2.5,
legend = geo_vetor_biomas$nome_bioma,
col = c("darkgreen", "orange", "orange4", "forestgreen",
"yellow", "yellow3"))
Figura 15.67: Mapa feito com a função plot() de um objeto vetor.
Para a classe dos objetos geoespaciais raster, a função plot() vai plotar um mapa para o tipo RasterLayer e quantos mapas houverem no objeto e couberem no espaço de plot para RasterBrick e RasterStack. Além disso, para essas classes do objeto raster, essa função provê também uma legenda e uma escala de cores automática (terrain). Vamos fazer o mapa da camada raster de elevação para o limite do município de Rio Claro/SP (Figura 15.68).
## Plot
plot(geo_raster_srtm_rio_claro)
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Figura 15.68: Mapa feito com a função plot() de um objeto raster com uma camada.
Agora vamos plotar objetos da classe RasterStack, alterando a cor para viridis, usando a função viridis::viridis() do pacote homônimo. Vamos fazer o mapa de duas camadas raster bioclimáticas para o mundo (Figura 15.69).
Figura 15.69: Mapa feito com a função plot() de um objeto raster com várias camadas.
Para exportar esses mapas podemos utilizar as funções png() ou pdf(), indicando os argumentos para ter as configurações que desejamos, e finalizando com a função dev.off(). Vamos exportar, a título de exemplo, a última figura (mais detalhes no Capítulo 6).
## Criar diretório
dir.create(here::here("dados", "mapas"))
## Exportar mapa
png(filename = here::here("dados", "mapas", "elev_rc.png"),
width = 20, height = 20, units = "cm", res = 300)
plot(geo_raster_srtm_rio_claro)
plot(geo_vetor_rio_claro$geom, col = NA, border = "red", lwd = 2, add = TRUE)
dev.off()Pacotes ggplot2 e ggspatial
Como discutimos no Capítulo 6 sobre gráficos, o pacote ggplot2 utiliza a gramática de gráficos para composição de figuras no R (Wilkinson and Wills 2005; Wickham 2016). Para cada classe de objeto geográfico há funções específicas para os dados: para objetos sf geom_sf() e para objetos raster geom_raster() e geom_tile().
Além do pacote ggplot2, podemos utilizar o pacote ggspatial para acrescentar elementos geoespaciais como a barra de escala e o indicador de orientação (Norte), através das funções annotation_scale() e annotation_north_arrow(), respectivamente, além de outras funções específicas que não abordaremos nesta seção.
A estrutura de composição das funções do pacote ggplot2 vai funcionar parecido com a estruturação de gráficos já vista no Capítulo 6, de modo que a cada função iremos utilizando o sinal de + para acrescentar outra camada. Indicaremos os dados com a função ggplot() e a coluna da tabela de atributos que queremos representar com a função aes(). Em seguida, utilizamos a função geom_sf() para indicar que trata-se de um objeto sf.
Além dessas funções, podemos ainda fazer alterações nos mapas através das funções: scale_*() que vai alterar as características indicadas em aes(), coord_*() que vai alterar construção do mapa em relação às coordenadas, facet_*() que altera a disposição de vários mapas, e theme_*() e theme() que alterarão características relacionadas ao tema, como cores de fundo, fontes e legenda. Podemos ainda utilizar as funções annotate() para adicionar textos e labs() para alterar o texto do título, legenda e eixos.
Vamos demonstrar esse funcionamento para compor o mapa de biomas, apresentado no início desta seção (Figura 15.70).
## Plot
mapa_vetor_biomas_ggplot2 <- ggplot(data = geo_vetor_biomas) +
aes(fill = nome_bioma) +
geom_sf(color = "black") +
scale_fill_manual(values = c("darkgreen", "orange", "orange4",
"forestgreen", "yellow", "yellow3")) +
annotation_scale(location = "br") +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(0, "cm"), pad_y = unit(.5, "cm"),
style = north_arrow_fancy_orienteering) +
annotate(geom = "text", label = "CRS: SIRGAS2000/Geo", x = -38, y = -31, size = 2.5) +
annotate(geom = "text", label = "Fonte: IBGE (2019)", x = -39, y = -32.5, size = 2.5) +
labs(title = "Biomas do Brasil", fill = "Legenda", x = "Longitude", y = "Latitude") +
theme_bw() +
theme(title = element_text(size = 15, face = "bold"),
legend.title = element_text(size = 10, face = "bold"),
legend.position = c(.15, .25),
legend.background = element_rect(colour = "black"),
axis.title = element_text(size = 10, face = "plain"))
mapa_vetor_biomas_ggplot2
Figura 15.70: Mapa de Biomas do Brasil com o pacote ggplot2.
Para objetos raster, o uso do pacote ggplot2 para compor mapas requer um passo preliminar. Primeiramente, vamos criar um data frame com os dados do raster, com as linhas sendo os pixels e as colunas sendo as coordenadas centrais da longitude e latitude, além dos valores de cada camada em cada coluna. Esse passo pode ser realizado com a função raster::rasterToPoints() ou raster::as.data.frame().
Uma vez que temos esses dados organizados, podemos utilizar as funções ggplot() para indicar o data frame, e as colunas com a função aes(). Em seguida, utilizamos a função geom_raster() para indicar que trata-se de um objeto raster. Além dessas funções, podemos ainda utilizar as demais funções para alterar as características do mapa, como comentamos acima. Entretanto, devemos nos atentar para a função coord_*() e escolher aquela que vai fazer a construção do mapa em relação à coordenadas e resolução das células.
Como exemplo, vamos compor o mapa de elevação para Rio Claro/SP, adicionando também o limite do município (Figura 15.71).
## Dados
geo_raster_srtm_rio_claro_dados <- raster::rasterToPoints(geo_raster_srtm_rio_claro) %>%
tibble::as_tibble()
head(geo_raster_srtm_rio_claro_dados)
#> # A tibble: 6 × 3
#> x y elevacao
#> <dbl> <dbl> <dbl>
#> 1 -47.8 -22.2 859
#> 2 -47.8 -22.2 856
#> 3 -47.8 -22.2 856
#> 4 -47.8 -22.2 856
#> 5 -47.8 -22.2 853
#> 6 -47.8 -22.2 852
## Mapa
mapa_srtm_rio_claro_ggplot2 <- ggplot() +
geom_raster(data = geo_raster_srtm_rio_claro_dados, aes(x = x, y = y, fill = elevacao)) +
geom_sf(data = geo_vetor_rio_claro, color = "red", fill = NA, size = 1.3) +
scale_fill_viridis_c() +
coord_sf() +
annotation_scale(location = "br",
pad_x = unit(.5, "cm"), pad_y = unit(.7, "cm"),) +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(.4, "cm"), pad_y = unit(1.3, "cm"),
style = north_arrow_fancy_orienteering) +
annotate(geom = "text", label = "CRS: WGS84/Geo",
x = -47.51, y = -22.53, size = 3) +
labs(title = "Elevação de Rio Claro/SP", fill = "Elevação (m)",
x = "Longitude", y = "Latitude") +
theme_bw() +
theme(title = element_text(size = 15, face = "bold"),
legend.title = element_text(size = 10, face = "bold"),
legend.position = c(.2, .25),
legend.background = element_rect(colour = "black"),
axis.title = element_text(size = 10, face = "plain"),
axis.text.y = element_text(angle = 90, hjust = .4))
mapa_srtm_rio_claro_ggplot2
Figura 15.71: Mapa raster com o pacote ggplot2.
Para exportar mapas criados com o pacote ggplot2, podemos utilizar a função ggplot2::ggsave(), indicando os argumentos para ter as configurações que desejamos. Vamos exportar, a título de exemplo, a última figura.
## Exportar mapa ggplot2
ggsave(filename = here::here("dados", "mapas", "srtm_rio_claro_ggplot2.png"),
plot = mapa_srtm_rio_claro_ggplot2, width = 20, height = 20, units = "cm", dpi = 300)Pacote tmap
O pacote tmap é um pacote direcionado à criação de mapas temáticos, com uma sintaxe concisa que permite a criação de mapas com o mínimo de códigos, mas muito similar à sintaxe do pacote ggplot2 (Tennekes 2018). Ele também pode gerar mapas estáticos ou interativos usando o mesmo código, apenas mudando a forma de visualização com a função tmap_mode(), com o argumento mode igual a “plot” para estático e “view” para interativo. Por fim, o pacote tmap aceita diversas classes espaciais, incluindo objetos raster, de forma bastante simples. Mais sobre o pacote pode ser lido aqui.
Todas as funções do pacote tmap iniciam-se com tm_*, facilitando seu uso. A cada função iremos utilizar o sinal de + para acrescentar outra camada, da mesma forma que o pacote ggplot2. A principal função, em que todos os objetos geoespaciais são dados de entrada, é tm_shape(). A partir dela, podemos seguir com funções específicas para visualização de objetos sf, como tm_polygons(), tm_borders(), tm_fill(), tm_lines(), tm_dots() ou tm_bubbles(), ou com funções para objetos raster como tm_raster(). Ainda há funções como tm_text() para representação de textos das colunas da tabela de atributos, e tm_scale_bar(), tm_compass() e tm_graticules(), para adicionar barra de escala, indicador de orientação (Norte) e gride de coordenadas, respectivamente. Por fim, a função tm_credits() adiciona um texto descritivo e a função tm_layout() faz diversas mudanças nos detalhes e apresentação do mapa.
Uma funcionalidade muito interessante do pacote tmap é o uso da função tmaptools::palette_explorer() para escolher as paletas de cores disponíveis. Essa função requer que os pacotes shiny e shinyjs estejam instalados, e quando executada, retorna uma aba onde é possível editar e escolher algumas paletas de cores nativas do tmap.
Diversos parâmetros podem ser acrescentados às funções de composição do tmap, mas não as detalharemos aqui, pois todas são descritas nos vignettes do pacote: tmap: get started! e tmap: version changes.
Vamos seguir com a composição do mapa de biomas para o Brasil apresentado no início dessa seção (Figura 15.72).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
mapa_vetor_biomas_tmap <- tm_shape(geo_vetor_biomas, bbox = c(-74, -35, -27, 10)) +
tm_polygons(col = "nome_bioma",
pal = c("darkgreen", "orange", "orange4", "forestgreen", "yellow", "yellow3"),
border.col = "black",
title = "Legenda") +
tm_compass() +
tm_scale_bar(text.size = .6) +
tm_graticules(lines = FALSE) +
tm_credits("CRS: SIRGAS2000/Geo", position = c(.63, .13)) +
tm_credits("Fonte: IBGE (2019)", position = c(.63, .09)) +
tm_layout(title = "Biomas do Brasil",
title.position = c(.25, .95),
title.size = 1.8,
title.fontface = "bold",
legend.frame = TRUE,
legend.position = c("left", "bottom"),
legend.title.fontface = "bold")
mapa_vetor_biomas_tmap
Figura 15.72: Mapa de Biomas do Brasil com o pacote tmap.
Além disso, o pacote tmap nos permite adicionar de forma simples um mapa secundário, provendo uma localização regional de interesse (Figura 15.73).
## Dados
geo_vetor_am_sul <- rnaturalearth::ne_countries(continent = "South America")
geo_vetor_brasil <- rnaturalearth::ne_countries(country = "Brazil")
geo_vetor_biomas <- geobr::read_biomes(showProgress = FALSE) %>%
dplyr::filter(name_biome != "Sistema Costeiro")Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_am_sul
ecodados::geo_vetor_brasil
ecodados::geo_vetor_biomas
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa secundário
mapa_am_sul <- tm_shape(geo_vetor_am_sul) +
tm_polygons() +
tm_shape(geo_vetor_brasil) +
tm_polygons(col = "gray50")
## Juntando os mapas
mapa_vetor_biomas_tmap
print(mapa_am_sul, vp = grid::viewport(.815, .875, wi = .2, he = .2))
Figura 15.73: Mapa vetorial primário e secundário com o pacote tmap.
Como exemplo de mapa raster, vamos compor novamente o mapa de elevação para Rio Claro/SP, adicionando também o limite do município (Figura 15.74).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
mapa_srtm_rio_claro_tmap <- tm_shape(geo_raster_srtm_rio_claro) +
tm_raster(pal = "viridis", title = "Elevação (m)") +
tm_shape(geo_vetor_rio_claro) +
tm_borders(col = "red", lwd = 2) +
tm_compass(position = c(.9, .08)) +
tm_scale_bar(text.size = .6, position = c(.67, 0)) +
tm_graticules(lines = FALSE) +
tm_credits("CRS: WGS84/Geo", position = c(.67, .06)) +
tm_layout(title = "Elevação Rio Claro/SP",
title.size = 1,
title.fontface = "bold",
legend.title.size = .7,
legend.text.size = .6,
legend.frame = TRUE,
legend.position = c(.01, .01),
legend.title.fontface = "bold")
mapa_srtm_rio_claro_tmap
Figura 15.74: Mapa raster de elevação com o pacote tmap.
Para exportar mapas criados com o pacote tmap podemos utilizar a função tmap::tmap_save(), indicando os argumentos para ter as configurações que desejamos. Vamos exportar, a título de exemplo, a última figura.
15.10.4 Mapas animados
Podemos montar mapas facetados para mostrar como padrões espaciais variam ao longo do tempo, como por exemplo, os limites do Brasil ao longo dos anos (Figura 15.75). Entretanto, essa abordagem possui algumas desvantagens, de modo que as facetas podem ficar muito pequenas quando há muitas delas.
## Dados
geo_vetor_brasil_anos <- NULL
for(i in c(1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970,
1980, 1991, 2001, 2010, 2019)){
geo_vetor_brasil_anos <- geobr::read_state(code_state = "all",
year = i, showProgress = FALSE) %>%
sf::st_geometry() %>%
sf::st_as_sf() %>%
dplyr::mutate(year = i) %>%
dplyr::bind_rows(geo_vetor_brasil_anos, .)
}Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_brasil_anos
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa facetado
mapa_brasil_tmap <- tm_shape(geo_vetor_brasil_anos) +
tm_polygons() +
tm_facets(by = "year", nrow = 4)
mapa_brasil_tmap
Figura 15.75: Mapa vetor facetado dos estados brasileiros ao longo do tempo com o pacote tmap.
Uma solução é a composição de mapas animados. Apesar de dependerem da publicação digital, os mapas animados podem aprimorar relatórios físicos à medida que o vínculo a uma página da web contendo a versão animada torna-se simples. Existem várias maneiras de gerar animações em R, por exemplo, com o pacote gganimate. Entretanto, aqui veremos a criação de mapas animados com tmap.
Podemos criar mapas animados alterando dois argumentos da função tm_facets():
- trocando o
by = yearporalong = year - indicando o
free.coords = FALSE
Por fim, podemos exportar o mapa animado no formato de .gif utilizando a função tmap::tmap_animation(), indicando a taxa de atualização com o argumento delay (Figura 15.76). Para visualização do .gif animado, acessar o livro on-line. Alguns pacotes extras são requeridos dependendo do sistema operacional utilizado.
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa animado
mapa_brasil_tmap_ani <- tm_shape(geo_vetor_brasil_anos) +
tm_polygons() +
tm_facets(along = "year", free.coords = FALSE)
## Exportar mapa tmap animado
tmap::tmap_animation(tm = mapa_brasil_tmap_ani,
filename = here::here("dados", "mapas",
"srtm_rio_claro_tmap_ani.gif"),
delay = 30)
Figura 15.76: Mapa vetorial animado mostrando os estados brasileiros ao longo do tempo com o pacote tmap.
15.10.5 Mapas interativos
Mapas interativos podem assumir muitas formas, sendo que a mais comum e útil é a capacidade de deslocar e ampliar qualquer parte de um conjunto de dados geoespaciais sobreposto em um “mapa da web”. Diversos pacotes nos permitem criar esse tipo de mapa, sendo os mais comuns o tmap, mapview e leaflet. Para visualização dos mapas interativos, acessar o livro on-line.
📝 Importante
Destacamos que esses mapas irão ser compostos numa janela especial de “Viewer”.
pacote tmap
Um recurso exclusivo do tmap é sua capacidade de criar mapas estáticos e interativos usando o mesmo código. Os mapas podem ser visualizados interativamente em qualquer ponto mudando para o modo de visualização, usando a função tmap::tmap_mode(mode = "view") (Figura 15.77).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "view")
## Atribuir novo mapa interativo
mapa_srtm_rio_claro_tmap_int <- mapa_srtm_rio_claro_tmap
Figura 15.77: Mapa vetorial interativo com o pacote tmap.
Para exportar mapas interativos criados com o pacote tmap, podemos utilizar novamente a função tmap::tmap_save(), indicando a extensão como .html.
## Exportar mapa tmap interativo
tmap::tmap_save(tm = mapa_srtm_rio_claro_tmap_int,
filename = here::here("dados", "mapas", "srtm_rio_claro_tmap_int.html"))Pacote mapview
O pacote mapview cria rapidamente mapas interativos simples com a função mapvew::mapview() (Figura 15.78). Entretanto, outras características podem ser mudadas para criar mapas mais elaborados, como pode ser visto através do site do pacote.
## Mapa
mapa_srtm_rio_claro_mapview_int <- mapview::mapview(
geo_raster_srtm_rio_claro, col.regions = viridis::viridis(100))
Figura 15.78: Mapa vetorial interativo com o pacote mapview.
Para exportar mapas interativos criados com o pacote mapview, podemos utilizar a função mapivew::mapshot(), indicando a extensão como .html.
## Exportar mapa mapview interativo
mapview::mapshot(x = mapa_srtm_rio_claro_mapview_int,
url = here::here("dados", "mapas", "srtm_rio_claro_mapview_int.html"))Pacote leaflet
O leaflet é um dos pacotes de mapeamento interativo mais utilizados e completos em R. Esse pacote fornece uma interface utilizando a biblioteca JavaScript e muitos argumentos podem ser compreendidos lendo a documentação da biblioteca original.
Mapas interativos usando esse pacote são criados utilizando a função leaflet::leaflet(). O resultado dessa função é um objeto da classe leaflet, que pode ser alterado por outras funções deste mesmo pacote, permitindo que várias camadas e configurações de controle sejam adicionadas interativamente (Figura 15.79).
Mais sobre o pacote leaflet pode ser consultado em seu site e CheatSheet.
## Paleta de cores
pal <- colorNumeric(viridis::viridis(10), raster::values(geo_raster_srtm_rio_claro))
## Mapa
mapa_srtm_rio_claro_leaflet_int <- leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(geo_raster_srtm_rio_claro, colors = pal, opacity = .8) %>%
addLegend(pal = pal, values = raster::values(geo_raster_srtm_rio_claro),
title = "Elevação (m)") %>%
addPolygons(data = geo_vetor_rio_claro, col = "red", fill = NA)
knitr::include_url("img/cap15_fig78.png")
Figura 15.79: Mapa vetorial interativo com o pacote leaflet.
Para exportar mapas interativos criados com o pacote leaflet, podemos utilizar novamente a função mapivew::mapshot(), indicando a extensão como .html.
15.11 Exemplos de aplicações de análises geoespaciais para dados ecológicos
Agora que vimos os principais conceitos e aplicações do manejo e visualização de dados geoespaciais, podemos avançar para realizar quatro exemplos de aplicações para dados ecológicos. Para isso, usaremos novamente os dados de comunidades de anfíbios da Mata Atlântica (Vancine et al. 2018). Primeiramente, veremos como resumir informações de biodiversidade (número de ocorrências e riqueza) para hexágonos. Num segundo momento, veremos como associar dados ambientais a coordenadas de espécies ou comunidades. Depois, como resumir dados de rasters para buffers. Por fim, realizaremos predições espaciais contínuas de adequabilidade de habitat para uma espécie e do número de espécies.
15.11.1 Resumir informações de biodiversidade para unidades espaciais
Resumir informações para unidades espaciais é um passo muito frequente em análises de Macroecologia, Biogeografia ou Ecologia da Paisagem. Nesta seção, contabilizaremos o número de ocorrências e riqueza de anfíbios para hexágonos na Mata Atlântica.
Primeiramente, vamos importar e preparar os dados de biodiversidade que usaremos nesses exemplos. Vamos começar importando os locais de amostragens de anfíbios na Mata Atlântica e selecionando apenas as colunas de interesse.
## Importar locais
geo_anfibios_locais <- readr::read_csv(
here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_sites.csv"),
col_types = cols()) %>%
dplyr::select(id, longitude, latitude, species_number)Agora vamos importar as espécies das comunidades, selecionando apenas as espécies com nomes válidos e transformando a coluna de indivíduos para 1, para compor posteriormente uma matriz de comunidade de espécies.
## Importar espécies
geo_anfibios_especies <- readr::read_csv(
here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_species.csv"),
col_types = cols()) %>%
tidyr::drop_na(valid_name) %>%
dplyr::select(id, valid_name, individuals) %>%
dplyr::distinct(id, valid_name, .keep_all = TRUE) %>%
dplyr::mutate(individuals = tidyr::replace_na(individuals, 1),
individuals = ifelse(individuals > 0, 1, 1))Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_anfibios_locais
ecodados::geo_anfibios_especiesPodemos agora juntar a tabela de locais, que possui as coordenadas à tabela de espécies. Em seguida convertemos essa única tabela na classe vetor sf.
## Junção das coordenadas e conversão para classe sf
geo_anfibios_especies_locais_vetor <- geo_anfibios_especies %>%
dplyr::left_join(geo_anfibios_locais, by = "id") %>%
dplyr::relocate(longitude, latitude, .after = 1) %>%
dplyr::mutate(lon = longitude, lat = latitude) %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)Agora vamos baixar o limite do Bioma da Mata Atlântica para o Brasil, converter o CRS para WGS84/Geo e ajustar sua extensão para remover as ilhas no Oceano Atlântico.
## Download do Bioma da Mata Atlântica
geo_vetor_mata_atlantica <- geobr::read_biomes(year = 2019, showProgress = FALSE) %>%
dplyr::filter(name_biome == "Mata Atlântica") %>%
sf::st_transform(crs = 4326) %>%
sf::st_crop(xmin = -55, ymin = -30, xmax = -34, ymax = -5)Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_vetor_mata_atlanticaPodemos verificar se as coordenadas e o limite do bioma estão todos corretos compondo um mapa preliminar, usando o pacote tmap (Figura 15.80).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_vetor_mata_atlantica,
bbox = geo_anfibios_especies_locais_vetor) +
tm_polygons() +
tm_shape(geo_anfibios_especies_locais_vetor) +
tm_dots(size = .1, col = "forestgreen")
Figura 15.80: Mapa dos locais do Atlantic Amphibians e do limite da Mata Atlântica.
Como o limite utilizado para reunir informações das comunidades de anfíbios foi o mais abrangente possível (Muylaert et al. 2018; Vancine et al. 2018), selecionaremos apenas os locais que caem dentro do limite da Mata Atlântica que estamos utilizando aqui.
## Selecionar os locais dentro do limite
geo_anfibios_especies_locais_vetor_mata_atlantica <- geo_anfibios_especies_locais_vetor[geo_vetor_mata_atlantica, ]Podemos refazer o mapa mostrando as coordenadas retiradas em vermelho e as que ficaram em verde (Figura 15.81).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_vetor_mata_atlantica,
bbox = geo_anfibios_especies_locais_vetor) +
tm_polygons() +
tm_shape(geo_anfibios_especies_locais_vetor) +
tm_bubbles(size = .1, col = "red") +
tm_shape(geo_anfibios_especies_locais_vetor_mata_atlantica) +
tm_bubbles(size = .1, col = "forestgreen")
Figura 15.81: Mapa dos locais do Atlantic Amphibians que caem dentro do limite da Mata Atlântica.
O próximo passo é criar um gride de hexágonos para o Bioma da Mata Atlântica. Usaremos a função sf::st_make_grid() que pode criar quadrículas ou hexágonos. Esses hexágonos terão a área equivalente à quadrículas de 1º de tamanho (aproximadamente 10000 km²). Usaremos a função sf::st_area() para calcular as áreas dos hexágonos e a função tibble::rowid_to_column() para criar uma identificação para cada feição.
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Criar hexágonos de ~110 km
geo_vetor_mata_atlantica_hex <- sf::st_make_grid(
x = geo_vetor_mata_atlantica, cellsize = 1, square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::mutate(areakm2 = sf::st_area(.)/1e6) %>%
tibble::rowid_to_column("id_hex")
## Selecionar os hexágonos dentro do limite da Mata Atlântica
geo_vetor_mata_atlantica_hex <- geo_vetor_mata_atlantica_hex[geo_vetor_mata_atlantica, ]Podemos conferir os hexágonos criados fazendo um mapa preliminar (Figura 15.82).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_vetor_mata_atlantica,
bbox = geo_vetor_mata_atlantica_hex) +
tm_polygons() +
tm_shape(geo_vetor_mata_atlantica_hex) +
tm_borders()
Figura 15.82: Mapa dos hexágonos para o limite da Mata Atlântica.
Podemos ser mais restritos e selecionar apenas os hexágonos dentro do limite do Bioma da Mata Atlântica utilizando o operador st_within() (Figura 15.83).
## Selecionar os hexágonos totalmente dentro do limite da Mata Atlântica
geo_vetor_mata_atlantica_hex_total <- geo_vetor_mata_atlantica_hex[geo_vetor_mata_atlantica, , op = st_within]
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_vetor_mata_atlantica,
bbox = geo_vetor_mata_atlantica_hex) +
tm_polygons() +
tm_shape(geo_vetor_mata_atlantica_hex_total) +
tm_borders()
Figura 15.83: Mapa dos hexágonos totalmente dentro do limite da Mata Atlântica.
Podemos agora associar as espécies aos hexágonos fazendo um “join” espacial, utilizando a função sf::st_join().
## Junção espacial dos locais com os hexágonos
geo_vetor_mata_atlantica_hex_especies <- sf::st_join(
x = geo_vetor_mata_atlantica_hex,
y = geo_anfibios_especies_locais_vetor_mata_atlantica,
left = TRUE)Por fim, podemos agregar os dados para ter o número de ocorrências e de espécies por hexágono.
## Agregar dados de ocorrências e número de espécies por hexágono
geo_vetor_mata_atlantica_hex_especies_oco_riq <- geo_vetor_mata_atlantica_hex_especies %>%
dplyr::group_by(id_hex) %>%
dplyr::summarise(ocorrencias = length(valid_name[!is.na(valid_name)]),
riqueza = n_distinct(valid_name, na.rm = TRUE))Finalmente podemos compor os mapas finais, mostrando os hexágonos com cores e valores do número de ocorrências e de espécies (Figura 15.84).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa de ocorrências
mapa_oco <- tm_shape(geo_vetor_mata_atlantica_hex_especies_oco_riq) +
tm_polygons(title = "Ocorrência de anfíbios", col = "ocorrencias",
pal = "viridis", style = "pretty") +
tm_text("ocorrencias", size = .4) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
## Mapa de riqueza
mapa_riq <- tm_shape(geo_vetor_mata_atlantica_hex_especies_oco_riq) +
tm_polygons(title = "Riqueza de anfíbios", col = "riqueza",
pal = "viridis", style = "pretty") +
tm_text("riqueza", size = .4) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
## União dos mapas
tmap_arrange(mapa_oco, mapa_riq)
Figura 15.84: Mapa com o número de ocorrências e riqueza de anfíbios para hexágonos no limite da Mata Atlântica.
15.11.2 Extrair dados raster para pontos
Atribuir informações ambientais a ocorrências é um passo fundamental para diversas análises. Nesta seção, atribuiremos os valores das variáveis bioclimáticas aos locais de amostragem de anfíbios na Mata Atlântica.
Já realizamos o download das variáveis bioclimáticas na seção de raster. Vamos importar novamente esses dados, primeiramente listando as camadas e depois importando com a função raster:stack().
## Listar arquivos
arquivos_raster <- dir(path = here::here("dados", "raster"), pattern = "wc") %>%
grep(".tif", ., value = TRUE)
## Importar rasters
geo_raster_bioclim <- raster::stack(here::here("dados", "raster", arquivos_raster))
## Renomear rasters
names(geo_raster_bioclim) <- c("bio01", paste0("bio", 10:19), paste0("bio0", 2:9))Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_raster_bioclimDa seção anterior, já temos o objeto com a tabela de coordenadas dos locais de amostragem das comunidades de anfíbios. Vamos agora criar um objeto vetorial das coordenadas e em seguida selecionar os locais dentro do limite do bioma da Mata Atlântica.
## Importar locais e converter em sf
geo_anfibios_locais_vetor <- geo_anfibios_locais %>%
dplyr::mutate(lon = longitude, lat = latitude) %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)
geo_anfibios_locais_vetor
#> Simple feature collection with 1163 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -56.74194 ymin: -33.51083 xmax: -34.79667 ymax: -3.51525
#> Geodetic CRS: WGS 84
#> # A tibble: 1,163 × 5
#> id longitude latitude species_number geometry
#> * <chr> <dbl> <dbl> <dbl> <POINT [°]>
#> 1 amp1001 -43.4 -8.68 19 (-43.42194 -8.68)
#> 2 amp1002 -38.9 -3.55 16 (-38.85783 -3.545527)
#> 3 amp1003 -38.9 -3.57 14 (-38.88869 -3.574194)
#> 4 amp1004 -38.9 -3.52 13 (-38.9188 -3.51525)
#> 5 amp1005 -38.9 -4.28 30 (-38.91083 -4.280556)
#> 6 amp1006 -36.4 -9.23 42 (-36.42806 -9.229167)
#> 7 amp1007 -40.9 -3.85 23 (-40.89444 -3.846111)
#> 8 amp1008 -40.9 -3.83 19 (-40.91944 -3.825833)
#> 9 amp1009 -40.9 -3.84 13 (-40.91028 -3.8375)
#> 10 amp1010 -35.2 -6.14 1 (-35.22944 -6.136944)
#> # … with 1,153 more rowsUsaremos agora a função raster::extract() para extrair e associar os valores das variáveis bioclimáticas para os locais de amostragem.
## Extrair valores das variáveis para os locais
geo_anfibios_locais_vetor_bioclim <- geo_anfibios_locais_vetor %>%
dplyr::mutate(raster::extract(geo_raster_bioclim, ., df = TRUE)) %>%
dplyr::select(-ID) %>%
dplyr::relocate(bio02:bio09, .after = bio01)Podemos ver esses dados na Tabela 15.14.
| id | longitude | latitude | species_number | bio01 | bio02 | bio03 | bio04 | bio05 |
|---|---|---|---|---|---|---|---|---|
| amp1001 | -43.42194 | -8.680000 | 19 | 24.88622 | 12.842188 | 76.28720 | 84.54608 | 33.17875 |
| amp1002 | -38.85783 | -3.545527 | 16 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1003 | -38.88869 | -3.574194 | 14 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1004 | -38.91880 | -3.515250 | 13 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1005 | -38.91083 | -4.280556 | 30 | 22.52411 | 8.434667 | 73.76507 | 64.62517 | 28.08925 |
| amp1006 | -36.42806 | -9.229167 | 42 | 21.61951 | 8.110271 | 67.74083 | 147.67953 | 27.98150 |
Podemos ainda fazer alguns mapas para espacializar essas variáveis (Figura 15.85).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
geo_anfibios_locais_vetor_bioclim %>%
dplyr::select(bio01:bio06) %>%
tidyr::gather(var, val, -geometry) %>%
tm_shape() +
tm_bubbles(size = .1, col = "val", pal = "viridis") +
tm_facets("var", free.scales = TRUE) +
tm_layout(legend.outside = FALSE)
Figura 15.85: Mapa com os valores das variáveis bioclimáticas (BIO01:BIO06) para os locais amostrados de comunidades de anfíbios para o limite da Mata Atlântica.
15.11.3 Resumir dados de rasters para buffers
Muitas análises requerem que façamos um resumo da composição da paisagem para buffers, sendo o buffer uma unidade de análise espacial no entorno de um ponto de amostragem.
Aqui, usaremos os dados do GlobCover v.2.3 de 2009 (Arino et al. 2012) como raster de cobertura da terra. O arquivo é grande (~400 Mb) e pode demorar muito, dependendo da velocidade da sua internet.
## Aumentar o tempo de download
options(timeout = 1e5)
## Download dos dados do GlobCover
download.file(url = "http://due.esrin.esa.int/files/Globcover2009_V2.3_Global_.zip",
destfile = here::here("dados", "raster", "Globcover2009_V2.3_Global.zip"),
mode = "wb", extra = "c")
## Unzip
unzip(zipfile = here::here("dados", "raster", "Globcover2009_V2.3_Global.zip"),
exdir = here::here("dados", "raster"))Depois de fazer o download, vamos importar e ajustar esse raster para o limite da Mata Atlântica (Figura 15.86).
## Importar raster do GlobCover
geo_raster_globcover <- raster::raster(
here::here("dados", "raster", "GLOBCOVER_L4_200901_200912_V2.3.tif"))
## Ajustar para o limite do bioma da Mata Atlântica
geo_raster_globcover_mata_atlantica <- geo_raster_globcover %>%
raster::crop(geo_vetor_mata_atlantica) %>%
raster::mask(geo_vetor_mata_atlantica)
geo_raster_globcover_mata_atlantica
#> class : RasterLayer
#> dimensions : 8940, 7274, 65029560 (nrow, ncol, ncell)
#> resolution : 0.002777778, 0.002777778 (x, y)
#> extent : -54.99861, -34.79306, -29.98194, -5.148611 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : GLOBCOVER_L4_200901_200912_V2.3
#> values : 14, 220 (min, max)Caso o download não funcione ou haja problemas com a importação, disponibilizamos os dados também no pacote ecodados.
## Importar os dados pelo pacote ecodados
ecodados::geo_raster_globcover_mata_atlantica
Figura 15.86: Camada raster do GlobCover 2.3 para o Bioma da Mata Atlântica.
Vamos agora transformar a tabela de locais em vetor, selecionar aleatoriamente 50 amostragens das comunidades de anfíbios e criar buffers de ~10 km.
## Fixar a amostragem
set.seed(42)
## Pontos
geo_anfibios_especies_locais_vetor_mata_atlantica <- geo_anfibios_locais %>%
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
dplyr::filter(sf::st_intersects(x = ., y = geo_vetor_mata_atlantica, sparse = FALSE)) %>%
dplyr::slice_sample(n = 50)
## Buffers de ~10 km
geo_anfibios_especies_locais_vetor_mata_atlantica_buffer10km <-
sf::st_buffer(geo_anfibios_especies_locais_vetor_mata_atlantica,
dist = 10000)Podemos conferir no mapa da Figura 15.87, atente para o aumento artificial do tamanho dos buffers para permitir a visualização.
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_vetor_mata_atlantica) +
tm_polygons() +
tm_shape(geo_anfibios_especies_locais_vetor_mata_atlantica_buffer10km) +
tm_bubbles(size = .3, border.col = "red", alpha = 0) +
tm_shape(geo_anfibios_especies_locais_vetor_mata_atlantica) +
tm_dots(size = .01, col = "forestgreen")
Figura 15.87: Distribuição de 50 localidades aleatórios e buffers de ~1 km (fora de escala).
Agora podemos utilizar a função raster::extract() para fazer a contabilização, já em porcentagem, de pixels de cada classe para cada buffer.
## Estatística zonal
geo_anfibios_locais_vetor_ma_buffer10km_ext <- raster::extract(
x = geo_raster_globcover_mata_atlantica,
y = geo_anfibios_especies_locais_vetor_mata_atlantica_buffer10km,
df = TRUE, na.rm = TRUE) %>%
dplyr::rename(id = ID, classe = 2) %>%
tidyr::drop_na(classe) %>%
dplyr::mutate(classe = paste0("classe_", classe)) %>%
dplyr::group_by(id, classe) %>%
dplyr::summarise(n = n()) %>%
tidyr::pivot_wider(id_cols = id, names_from = classe, values_from = n) %>%
dplyr::mutate(across(everything(), ~replace_na(.x, 0))) %>%
janitor::adorn_totals("col") %>%
janitor::adorn_percentages("row") %>%
janitor::adorn_pct_formatting(rounding = "half up", digits = 1)
head(geo_anfibios_locais_vetor_ma_buffer10km_ext)
#> id classe_100 classe_110 classe_120 classe_130 classe_14 classe_150 classe_170 classe_180 classe_20 classe_210 classe_30 classe_40 classe_50
#> 1 0.1% 0.2% 0.7% 12.0% 0.5% 0.1% 8.7% 0.6% 0.6% 24.8% 5.0% 43.7% 3.1%
#> 2 0.0% 1.7% 0.3% 8.4% 26.6% 0.5% 0.0% 0.0% 14.7% 1.5% 19.2% 17.9% 8.9%
#> 3 0.0% 0.0% 0.0% 4.9% 20.2% 0.0% 0.0% 0.0% 17.6% 0.2% 28.6% 27.2% 0.6%
#> 4 0.0% 2.3% 0.0% 26.0% 3.6% 0.0% 0.0% 0.0% 23.0% 0.2% 41.2% 3.6% 0.0%
#> 5 0.0% 0.0% 0.0% 13.5% 0.3% 0.0% 0.0% 0.1% 14.2% 0.0% 5.2% 53.6% 12.9%
#> 6 0.0% 0.4% 0.3% 7.2% 7.0% 0.0% 0.0% 0.0% 21.2% 0.0% 21.1% 35.8% 1.4%
#> classe_190 classe_60 classe_200 classe_220 classe_140 Total
#> 0.0% 0.0% 0.0% 0.0% 0.0% 100.0%
#> 0.4% 0.0% 0.0% 0.0% 0.0% 100.0%
#> 0.0% 0.6% 0.0% 0.0% 0.0% 100.0%
#> 0.0% 0.0% 0.0% 0.0% 0.0% 100.0%
#> 0.0% 0.1% 0.0% 0.0% 0.0% 100.0%
#> 5.2% 0.2% 0.1% 0.0% 0.0% 100.0%Agora podemos combinar esses dados aos dados dos buffers.
## Combinação
geo_anfibios_locais_vetor_ma_buffer10km_ext_bind <- dplyr::bind_cols(
x = geo_anfibios_especies_locais_vetor_mata_atlantica_buffer10km,
y = geo_anfibios_locais_vetor_ma_buffer10km_ext[, -1]) %>%
sf::st_drop_geometry()
geo_anfibios_locais_vetor_ma_buffer10km_ext_bind
#> # A tibble: 50 × 21
#> id species_number classe_100 classe_110 classe_120 classe_130 classe_14 classe_150 classe_170 classe_180 classe_20 classe_210 classe_30
#> * <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 amp1716 10 0.1% 0.2% 0.7% 12.0% 0.5% 0.1% 8.7% 0.6% 0.6% 24.8% 5.0%
#> 2 amp1351 16 0.0% 1.7% 0.3% 8.4% 26.6% 0.5% 0.0% 0.0% 14.7% 1.5% 19.2%
#> 3 amp1168 19 0.0% 0.0% 0.0% 4.9% 20.2% 0.0% 0.0% 0.0% 17.6% 0.2% 28.6%
#> 4 amp1085 20 0.0% 2.3% 0.0% 26.0% 3.6% 0.0% 0.0% 0.0% 23.0% 0.2% 41.2%
#> 5 amp1258 3 0.0% 0.0% 0.0% 13.5% 0.3% 0.0% 0.0% 0.1% 14.2% 0.0% 5.2%
#> 6 amp1160 9 0.0% 0.4% 0.3% 7.2% 7.0% 0.0% 0.0% 0.0% 21.2% 0.0% 21.1%
#> 7 amp1843 14 0.0% 0.0% 0.1% 3.9% 4.3% 0.0% 0.0% 0.0% 22.6% 0.6% 6.2%
#> 8 amp1060 8 0.0% 0.8% 0.0% 8.4% 11.5% 0.0% 0.0% 0.0% 53.3% 0.5% 22.7%
#> 9 amp1141 5 0.0% 0.2% 0.1% 2.6% 24.5% 0.0% 0.0% 0.1% 42.9% 0.0% 21.9%
#> 10 amp1333 7 0.0% 1.2% 0.2% 9.9% 8.5% 0.0% 0.0% 0.2% 21.9% 1.6% 21.0%
#> # … with 40 more rows, and 8 more variables: classe_40 <chr>, classe_50 <chr>, classe_190 <chr>, classe_60 <chr>, classe_200 <chr>,
#> # classe_220 <chr>, classe_140 <chr>, Total <chr>15.11.4 Predições espaciais de objetos raster
O pacote raster além de permitir realizar a manipulação e visualização de dados raster no R, também permite a extrapolação do ajuste de análises, como LMs, GLMs, GAMs dentre outras (Capítulos 7 e 8). Aqui, faremos uma pequena demonstração utilizando a função raster::predict(), predizendo o resultado de dois ajustes de GLMs para a presença/ausência de uma espécie de anuro e a extrapolação do número de espécies de anfíbios para o Bioma da Mata Atlântica.
Para ajustar um GLM para dados de presença/ausência, podemos usar o conjunto de dados já criado anteriormente, com as espécies e as coordenadas, e fazer uma junção com a última tabela que criamos com os dados bioclimáticos.
## Junção dos dados ambientais aos dados de espécies
geo_anfibios_locais_especies_vetor_bioclim <- geo_anfibios_especies %>%
dplyr::left_join(., sf::st_drop_geometry(geo_anfibios_locais_vetor_bioclim), by = "id")Agora, vamos selecionar ocorrências da espécie Haddadus binotatus, atribuindo 1 quando ela ocorre e 0 quando ela não ocorre. Essa espécie é relativamente comum na serrapilheira de fragmentos florestais da Mata Atlântica, e recebe esse nome em homenagem a um grande pesquisador de anfíbios da Mata Atlântica, o Prof. Célio Fernando Baptista Haddad.
## Seleção da espécie Haddadus binotatus
geo_anfibios_locais_especies_vetor_bioclim_hb <- geo_anfibios_locais_especies_vetor_bioclim %>%
dplyr::mutate(pa = ifelse(valid_name == "Haddadus binotatus", 1, 0), .after = individuals) %>%
dplyr::distinct(id, .keep_all = TRUE)Vamos utilizar apenas as variáveis não correlacionadas para o índice de correlação de Pearson para r < 0,7 (Capítulo 7).
## Correlação entre as variáveis
corr <- geo_anfibios_locais_especies_vetor_bioclim_hb %>%
dplyr::select(bio01:bio19) %>%
cor() %>%
caret::findCorrelation(.7, names = TRUE)
## Seleção das variáveis não correlacionadas
geo_anfibios_locais_especies_vetor_bioclim_hb_cor <- geo_anfibios_locais_especies_vetor_bioclim_hb %>%
dplyr::select(pa, bio01:bio19) %>%
dplyr::select(-c(corr))Agora sim, podemos ajustar um modelo simples da presença e ausência dessa espécie, utilizando as variáveis não correlacionadas, através de um GLM para a família binomial. Nosso intuito não é analisar se o modelo atende à todos os pressupostos, e sim exemplificar a predição espacial, para esses detalhes, consulte o Capítulo 8.
## Ajustar um modelo GLM binomial
modelo_pa <- glm(formula = pa ~ .,
data = geo_anfibios_locais_especies_vetor_bioclim_hb_cor,
family = binomial("logit"))Antes de fazermos a predição da distribuição potencial da espécie é fundamental que o objeto raster esteja ajustado para o limite da Mata Atlântica. Para isso vamos utilizar as funções raster::crop() e raster::mask() para fazer esse ajuste (Figura 15.88).
## Ajuste da extensão e limite
geo_raster_bioclim_mata_atlantica <- geo_raster_bioclim %>%
raster::crop(geo_vetor_mata_atlantica) %>%
raster::mask(geo_vetor_mata_atlantica)
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Mapa
tm_shape(geo_raster_bioclim_mata_atlantica[[c(1, 4)]]) +
tm_raster(pal = "viridis", title = c("bio01", "bio12")) +
tm_facets(free.scales.raster = TRUE)
Figura 15.88: Mapa de dois rasters (BIO01 e BIO12) ajustados ao limite da Mata Atlântica.
Agora podemos fazer a predição desse modelo para todo o Bioma da Mata Atlântica. Essa função vai utilizar os coeficientes do modelo ajustado para gerar um raster de predição para todos os pixels da Mata Atlântica. Vamos usar o argumento type = "response" para que os valores da predição sejam ajustados a 0 a 1.
📝 Importante
Para que a predição funcione, os nomes das camadas raster geo_raster_bioclim_mata_atlantica devem possuir o mesmo nome das colunas das variáveis preditoras ajustadas no modelo modelo_pa.
## Predições
modelo_pa_pred <- raster::predict(
object = geo_raster_bioclim_mata_atlantica,
model = modelo_pa,
type = "response")
modelo_pa_pred
#> class : RasterLayer
#> dimensions : 149, 121, 18029 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -55, -34.83333, -30, -5.166667 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : layer
#> values : 1.255833e-12, 0.1983049 (min, max)Por fim, no último passo podemos tornar esse modelo binário, ou seja, apenas com valores 0 ou 1. Para isso vamos adotar arbitrariamente o valor de 0,01 como ponto de corte. A partir desse valor consideraremos os pixels acima como 1 e abaixo como 0.
## Seleção dos pixels de presença/ausência potencial
modelo_pa_pred_corte <- modelo_pa_pred >= .01Por fim, vamos produzir dois mapas mostrando os valores das predições e o mapa binário da espécie (Figura 15.89).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Combinar os dois raster
modelo_pa_pred <- raster::stack(modelo_pa_pred, modelo_pa_pred_corte)
names(modelo_pa_pred) <- c("Contínuo", "Binário")
## Mapas de predição contínua e binária
tm_shape(modelo_pa_pred) +
tm_raster(pal = "viridis", title = c("Predição", "Predição"),
style = "fisher") +
tm_facets(free.scales.raster = TRUE)
Figura 15.89: Predição contínua e binária do modelo ajustado para a presença/ausência da espécie Haddadus binotatus na Mata Atlântica.
📝 Importante
Essa análise foi realizada com o intuito de exemplificar o funcionamento da função raster::predict(), para mais detalhes, consultar livros específicos da área de Modelagem de Distribuição de Espécies ou Modelagem de Nicho Ecológico Fletcher and Fortin (2018).
Em nossa segunda análise, vamos predizer os dados de riqueza de anfíbios (i.e. número de espécies) para todo o bioma da Mata Atlântica. Para isso, temos de retirar novamente as variáveis correlacionadas.
## Correlação
corr <- geo_anfibios_locais_vetor_bioclim %>%
sf::st_drop_geometry() %>%
dplyr::select(bio01:bio19) %>%
cor() %>%
caret::findCorrelation(.7, names = TRUE)
## Seleção das variáveis não correlacionadas
geo_anfibios_locais_bioclim_cor <- geo_anfibios_locais_vetor_bioclim %>%
sf::st_drop_geometry() %>%
dplyr::select(species_number, bio01:bio19) %>%
dplyr::select(-c(corr))Agora sim, podemos criar os GLMs com famílias de distribuição apropriadas para dados de contagem como Poisson e Binomial Negativa. Novamente, nosso intuito não é analisar se o modelo atende a todos os pressupostos, e sim exemplificar a predição espacial; para esses detalhes, consulte o Capítulo 8.
## Modelo Poisson
modelo_riq_pois <- glm(
formula = species_number ~ .,
data = geo_anfibios_locais_bioclim_cor,
family = poisson)
## Modelo Binomial Negativo
modelo_riq_nb <- MASS::glm.nb(
formula = species_number ~ .,
data = geo_anfibios_locais_bioclim_cor)Com os modelos ajustados, podemos fazer as predições utilizando os objetos raster com as variáveis ambientais.
## Predição do modelo Poisson
modelo_riq_pois_pred <- raster::predict(
object = geo_raster_bioclim_mata_atlantica,
model = modelo_riq_pois,
type = "response")
modelo_riq_pois_pred
#> class : RasterLayer
#> dimensions : 149, 121, 18029 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -55, -34.83333, -30, -5.166667 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : layer
#> values : 8.249131, 24.76082 (min, max)
## Predição do modelo Binomial Negativo
modelo_riq_nb_pred <- raster::predict(
object = geo_raster_bioclim_mata_atlantica,
model = modelo_riq_nb,
type = "response")
modelo_riq_nb_pred
#> class : RasterLayer
#> dimensions : 149, 121, 18029 (nrow, ncol, ncell)
#> resolution : 0.1666667, 0.1666667 (x, y)
#> extent : -55, -34.83333, -30, -5.166667 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : layer
#> values : 9.262095, 24.45297 (min, max)Por fim, podemos compor os dois mapas de predições (Figura 15.90).
## Mudar o modo de exibição do tmap
tmap::tmap_mode(mode = "plot")
## Combinar os dois raster
modelo_riq_pred <- raster::stack(modelo_riq_pois_pred, modelo_riq_nb_pred)
names(modelo_riq_pred) <- c("Poisson", "Binomial Negativa")
## Mapas da predição Poisson e Binomial Negativa
tm_shape(modelo_riq_pred) +
tm_raster(pal = "viridis", title = c("Número de espécies", "Número de espécies"),
style = "fisher") +
tm_facets(free.scales.raster = TRUE)
Figura 15.90: Mapa da predição de riqueza utilizando o modelo Poisson e Binomial Negativa para a Mata Atlântica.
📝 Importante
Essa análise foi realizada com o intuito de exemplificar o funcionamento da função raster::predict(), para mais detalhes, consultar livros específicos da área de Ecologia Espacial (Fletcher and Fortin 2018).
15.12 Para se aprofundar
15.12.1 Livros
Listamos aqui as principais referências sobre manipulação, visualização de dados geoespaciais e análises geoespaciais no R. Recomendamos aos (às) interessados(as) os livros: i) Lovelace, Nowosad & Muenchow (2019) Geocomputation with R, ii) Mas et al. (2019) Análise espacial com R, iii) Olaya (2020) Sistemas de Información Geográfica, iv) Moraga (2019) Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny, v) Brunsdon & Comber (2015) An Introduction to Spatial Analysis and Mapping in R, vi) Wegmann, Leutner & Dech (2016) Remote Sensing and GIS for Ecologists: Using Open Source Software, vii) Wegmann, Schwalb-Willmann & Dech (2020) An Introduction to Spatial Data Analysis Remote Sensing and GIS with Open Source Software, viii) Fletcher & Fortin (2018) Spatial ecology and conservation modeling: Applications with R, ix) Lapaine, Miljenko, Usery, E. Lynn (2017) Choosing a Map Projection.
15.13 Exercícios
15.1
Importe o limite dos estados brasileiros no formato sf com o nome br. Para isso, use a função ne_states do pacote rnaturalearth. Crie um mapa simples cinza utilizando a função plot(), selecionando a coluna geometry com o operador $ e com os argumentos axes e graticule verdadeiros.
15.2
Dados vetoriais podem ser criados com diversos erros de topologia, e.g., sobreposição de linhas ou polígonos ou buracos. Algumas funções exigem que os objetos vetoriais aos quais são atribuídos esses dados não possuam esses erros para que o algoritmo funcione. Para verificar se há erros, podemos usar a função st_is_valid() do pacote sf. Há diversas forma de correções desses erros, mas vamos usar uma correção simples do R, com a função st_make_valid(). Vamos fazer essa correção para o br importado anteriormente e atribuindo ao objeto br_valid. Podemos conferir para saber se há erros e fazer um plot.
15.3
Crie um objeto RasterLayer vazio chamado ra com resolução de 5º (~600 km). Atribua um sistema de referência de coordenadas com o código 4326. Atribua valores aleatórios de uma distribuição normal e plote o mesmo.
15.4
Reprojete o limite dos estados brasileiros do exercício anterior para o CRS SIRGAS 2000/Brazil Polyconic, utilizando o código EPSG:5880 e chamando de br_poly. Faça um mapa simples como no exercício 1. Atente para as curvaturas das linhas.
15.5
Utilizando a função st_centroid do pacote sf, crie um vetor chamado br_valid_cen que armazenará o centroide de cada estado brasileiro do objeto br_valid do exercício 2 e plot o resultado.
15.6
Ajuste o limite e máscara do objeto raster criado no exercício 3 para o limite do Brasil, atribuindo ao objeto ra_br. Depois reprojete esse raster para a mesma projeção utilizada no exercício 4 com o nome ra_br_poly e plote o mapa resultante.
15.7
Extraia os valores de cada pixel do raster criado no exercício 6 para os centroides dos estados do Brasil criado no exercício 5, atribuindo à coluna val do objeto espacial chamado br_valid_poly_cent_ra.
15.8
Crie um mapa final usando os resultados dos exercícios 4, 5 e 6. Utilize o pacote tmap e inclua todos os principais elementos de um mapa.