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_atlantica

15.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].

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: epsg ou proj4string indicando o CRS (Coordinate Reference System)
  • tibble: tabela de atributos, com destaque para a coluna geom ou geometry que 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 rows

Podemos 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.

## Plot dos polígonos do mundo
plot(world[1], col = viridis::viridis(100), main = "Mapa do mundo")
Mapa vetorial do mundo.

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].

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].

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  105

Vamos 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).

## Plot raster layers
plot(raster_layer, col = viridis::viridis(n = 100))
Mapa simples de um `RasterLayer`.

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.290035

Ao utilizarmos a função plot() do pacote raster, podemos visualizar os raster contidos no objeto RasterBrick (Figura 15.6).

## Plot raster brick
plot(raster_brick, col = viridis::viridis(n = 25), main = "")
Mapas simples de um raster `RasterBrick`.

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.290035

Da mesma forma, ao utilizar a função plot() do pacote raster, podemos visualizar os raster contidos no objeto RasterStack (Figura 15.7).

## Plot raster stack
plot(raster_stack, col = viridis::viridis(n = 25), main = "")
Mapas simples de um raster `RasterStack`.

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)
Tabela 15.1: Listagem de Sistemas de Referências de Coordenadas disponíveis no R, com informações dos códigos epsg e proj4string
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.

Tabela 15.2: Principais bases de dados vetoriais para o Brasil e o Mundo.
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.

Tabela 15.3: Principais bases de dados raster para o Brasil e o Mundo.
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.

Tabela 15.4: Principais pacotes no R para download de dados vetoriais e raster.
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.

Tabela 15.5: Principais formatos de arquivos geoespaciais. Adaptado de: Lovelace et al. (2019).
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())
Tabela 15.6: Alguns formatos vetoriais importados e exportados pelo pacote sf.
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_cobertura

Agora 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)
Mapa de nascentes de Rio Claro/SP.

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)
Mapa da hidrografia de Rio Claro/SP.

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"))
Mapa de cobertura da terra de Rio Claro/SP.

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()
Tabela 15.7: Alguns dados disponíveis no pacote 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)
Limite do município de Rio Claro/SP.

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)
Limite do Brasil.

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_locais

Por 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)
Coordenadas das comunidades do Atlantic Amphibians.

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().

## Polígonos países sf
co110_sf <- sf::st_as_sf(co110_sp)
class(co110_sf)
#> [1] "sf"         "data.frame"

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_srtm

Agora 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)
Camada raster do DEM em relação ao limite do município de Rio Claro/SP.

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_bioclim

Para 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, ...
## Plot
plot(geo_raster_bioclim[[c(1, 4)]], col = viridis::viridis(10))
Camadas rasters do WorldClim (BIO01 e BIO12) para o mundo.

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.

Tabela 15.8: Principais formatos de arquivos raster exportados do R.
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.

Tabela 15.9: Tipos de dados suportados pelo pacote raster.
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     Sudeste

15.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 383

15.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] TRUE

As 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.

Tabela 15.10: Principais CRSs utilizados.
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)
Limites do município de Rio Claro/SP com CRS SIRGAS2000/geográfico e com CRS SIRGAS2000/UTM23S.

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)
Limite dos países do mundo com CRS geográfico e datum WGS84.

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)
Limite dos países do mundo com CRS Projeção de Mollweide.

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)
Limite dos países do mundo com CRS Projeção Azimutal de Lambert centrado no Brasil.

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)
Ajuste da extensão do raster de elevação para o município de Rio Claro/SP.

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)
Reprojeção do raster de elevação para SIRGAS2000/UTM23S especificado por um objeto e informando o tamanho da célula.

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)
Reprojeção do raster de média de temperatura anual (BIO01) para Projeção de Mollweide informando o tamanho da célula.

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.

Tabela 15.11: Principais funções para realizar operações de atributos e suas descrições.
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)
Filtro da classe floresta para o mapeamento de cobertura da terra para o município de Rio Claro/SP.

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 silvicultura

Agregaçã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)
Agregação e contagem das nascentes para o município de Rio Claro/SP.

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              NA

Para 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)
# Plot
plot(geo_raster_bioclim_bio01, col = viridis::viridis(10))
Camada BIO01 selecionada pelas operações de subconjunto.

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.3774489

Ou 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  7892

15.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).

Tabela 15.12: Principais pacotes para composição de mapas no R.
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)
Nascentes dentro de florestas no município de Rio Claro/SP.

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)
Nascentes fora de florestas no município de Rio Claro/SP.

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"))
Junção espacial da cobertura da terra para as nascentes no município de Rio Claro/SP.

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)
Agregação espacial contabilizando o número de nascentes para cada classe de cobertura da terra no município de Rio Claro/SP.

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)
Distância espacial das nascentes até o fragmento de floresta mais próxima no município de Rio Claro/SP.

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)
Rasters de soma e log10 do mapa de elevação para Rio Claro/SP.

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)
Operação local lógica mostrando todos os pixels acima de 600 metros de elevação para Rio Claro/SP.

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)
Operação local de multiplicação de todos os pixels por eles mesmos do raster de elevação para Rio Claro/SP.

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)
Operação local de reclassificação para três classes de elevação para Rio Claro/SP.

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)
Cálculo do desvio padrão da elevação para uma janela de 3x3 e do aspecto para Rio Claro/SP.

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 985

Operaçõ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)
Raster de distância Euclidiana dos pixels abaixo de 500 m de elevação para Rio Claro/SP.

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)
Simplificação da hidrografia para Rio Claro/SP. As linhas em azul é a hidrografia original e as linhas em preto mostram a simplificação.

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)
Centroide do limite do município de Rio Claro/SP.

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)
Sorteio de 30 pontos aleatório para Rio Claro/SP.

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)
Buffers de 1000 metros para os 30 pontos aleatórios no município de Rio Claro/SP.

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)
Buffers quadrados de 1000 metros para os 30 pontos aleatórios no município de Rio Claro/SP.

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)
Polígono convexo para os 30 pontos criados aleatoriamente para Rio Claro/SP.

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)
Polígonos de Voronoi para os 30 pontos criados aleatoriamente para Rio Claro/SP.

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)
Quadrículas de 2000 metros de arestas e centroides para Rio Claro/SP.

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)
Hexágonos equivalentes a quadrículas de 2000 metros de arestas e centroides para Rio Claro/SP.

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)
União - dissolução - dos buffers criados a partir dos 30 pontos aleatórios para Rio Claro/SP.

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)
Recorte da hidrografia para fora dos buffers dos 30 pontos aleatórios para Rio Claro/SP.

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)
Recorte da hidrografia para fora dos buffers dos 30 aleatórios para Rio Claro/SP.

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)
Transformação do vetor de florestas em `POLYGON` e cálculo da área para cada feição para Rio Claro/SP.

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)
Agregação (aumento do pixel para 900 metros) utilizando a média para o raster de elevação para Rio Claro/SP.

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)
Desagregação (diminuição do pixel para 45 metros) utilizando o método bilinear para o raster de elevação para Rio Claro/SP.

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)
Reamostragem (alinhamento dos rasters) utilizando o método bilinear para alinhar o raster bioclimático (BIO01) ao de elevação para Rio Claro/SP.

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)
Ajuste da extensão do raster de elevação para a extensão de Rio Claro/SP.

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)
Ajuste do raster de elevação para o limite de Rio Claro/SP.

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)
Ajuste da extensão e do limite do raster de elevação para Rio Claro/SP.

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)
Ajuste da extensão e do limite externo do raster de elevação para Rio Claro/SP.

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)
Extração dos valores de elevação para as nascentes de Rio Claro/SP.

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)
Extração dos valores de elevação e resumo pela média para os hexágonos de Rio Claro/SP.

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)
Rasterização das nascentes, com a operação de contabilização de pontos para Rio Claro/SP.

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)
Rasterização da hidrografia, com a operação de contabilização de linhas para Rio Claro/SP.

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)
Rasterização da cobertura da terra para Rio Claro/SP.

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)
Vetorização do raster de elevação criando um gride de pontos para Rio Claro/SP.

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)
Vetorização do raster de elevação criando isolinhas para Rio Claro/SP.

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)
Vetorização do raster de cobertura da terra para Rio Claro/SP, não dissolvendo e dissolvendos os polígonos gerados.

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:

  1. Mapa principal (ocupando quase toda a área da figura)
  2. 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)
  3. Título (para resumir o intuito do mapa)
  4. Legenda (apresentando as informações detalhadas das classes ou escala de valores, geralmente identificando as cores e/ou texturas),
  5. Barra de escala (representando a relação entre as unidades do mapa e o mundo real)
  6. Indicador de orientação (Norte) (indicando o norte geográfico, podendo ser representado por uma flecha, bússola ou compasso)
  7. Gride de coordenadas (coordenadas presentes nas laterais)
  8. Descrição do CRS (indicando qual o Sistema de Referência de Coordenadas)
  9. Informações de origem (informações sobre a fonte dos dados representados no mapa)
  10. Outros elementos auxiliares (como elementos textuais e figuras extras)

Podemos visualizar todos esses elementos resumidos na Figura 15.65.

Principais elementos de um mapa.

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).

Tabela 15.13: Principais pacotes para composição de mapas no R.
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_biomas

Agora, 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)
Mapa feito com a função `plot()` de um objeto `sf` para os Biomas do Brasil.

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"))
Mapa feito com a função `plot()` de um objeto vetor.

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)
Mapa feito com a função `plot()` de um objeto raster com uma camada.

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).

## Plot
plot(geo_raster_bioclim[[1:2]], col = viridis::viridis(10))
Mapa feito com a função `plot()` de um objeto raster com várias camadas.

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
Mapa de Biomas do Brasil com o pacote `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
Mapa raster com o pacote `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.

📝 Importante
Atentar para a instalação extra em Sistemas Operacionais GNU/Linux e MacOS.

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
Mapa de Biomas do Brasil com o pacote `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))
Mapa vetorial primário e secundário com o pacote `tmap`.

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
Mapa raster de elevação com o pacote `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.

## Exportar mapa tmap
tmap::tmap_save(tm = mapa_srtm_rio_claro_tmap, 
                filename = here::here("dados", "mapas", "srtm_rio_claro_tmap.png"),
                width = 20, height = 20, units = "cm", dpi = 300)

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
Mapa vetor facetado dos estados brasileiros ao longo do tempo com o pacote `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 = year por along = 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)
Mapa vetorial animado mostrando os estados brasileiros ao longo do tempo com o pacote `tmap`.

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.

## Exportar mapa leaflet interativo
mapview::mapshot(x = mapa_srtm_rio_claro_leaflet_int, 
                 url = here::here("dados", "mapas", "srtm_rio_claro_leaflet_int.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_especies

Podemos 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_atlantica

Podemos 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")
Mapa dos locais do Atlantic Amphibians e do limite da Mata Atlântica.

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")
Mapa dos locais do Atlantic Amphibians que caem dentro do limite da Mata Atlântica.

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()
Mapa dos hexágonos para o limite da Mata Atlântica.

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()
Mapa dos hexágonos totalmente dentro do limite da Mata Atlântica.

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)
Mapa com o número de ocorrências e riqueza de anfíbios para hexágonos no limite da Mata Atlântica.

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_bioclim

Da 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 rows

Usaremos 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.

Tabela 15.14: Dados extraídos e atribuídos aos locais de amostragens de comunidades de anfíbios na Mata Atlântica
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)
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.

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
## Plot
plot(geo_raster_globcover_mata_atlantica, col = viridis::viridis(n = 200))
Camada raster do GlobCover 2.3 para o Bioma da Mata Atlântica.

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")
Distribuição de 50 localidades aleatórios e buffers de ~1 km (fora de escala).

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)
Mapa de dois rasters (BIO01 e BIO12) ajustados ao limite da Mata Atlântica.

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 >= .01

Por 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)
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.

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)
Mapa da predição de riqueza utilizando o modelo Poisson e Binomial Negativa para a Mata Atlântica.

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.

Soluções dos exercícios.