# Vetor contendo nomes dos pacotes; tidyverse inclui ggplot2
pac_nomes <- c("sf", "stars", "tidyverse", "ggspatial")
# Verifica se cada pacote está instalado e, se não, instala
for (item in pac_nomes) { if (!item %in% installed.packages()[,"Package"]) {install.packages(item)}
# Carrega na memória RAM cada pacote (cada item em pac_nomes)
library(item, character.only = TRUE)}Mapas e rotinas geo-espacials no R
1. Instalar R, seguido de RStudio
Baixar R-4.3.1 for Windows (79 megabytes, 64 bit) – Instalar primeiro
Baixar e instalar RStudio (ou atualize sua instalação)
2. Consulte livro online
3. Baixar pasta do projeto
Baixar a pasta ‘geospat’, neste link de GDrive.
Baixar a pasta inteira. Na página de destino do link clicar no diretório ‘geospat’ com o botão direito do mouse e selecione ‘download’. Deve cair na sua pasta de Downloads. Retire do zip e coloque a pasta ‘geospat’ em qualquer lugar no seu micro.
3. Ativar o projeto e abrir o script QMD
No RStudio: File - Open Project - navegue até o local de sua pasta ‘geospat’ e clicar duas vezes no arquivo ‘geospat.Rproj’. File - Open File - clicar duas vezes no arquivo ‘mapas_e_rotinas_geoespaciais_no_R.qmd’.
Agora seus caminhos dentro do script qmd serão relativos a esta pasta. A estrutura de pastas e subpastas acima da pasta de projeto não deve ser incluída nos caminhos.
Para rodar uma ou mais linhas no script, selecione-as e tecle Ctrl-Enter
4. Algumas rotinas básicas geo-espaciais
11 video-aulas do grupo datacarpentry.org; dados do projeto NEON.
Acessar aqui, mas siga as rotinas modificadas neste arquivo QMD
AULA 1. LER DADOS RASTER; PRIMEIRO MAPA
1.1. Instale e carregue pacotes
Rodar este primeiro bloco de código usando ‘Run current chunk’
1.2. Importar raster geotiff para o formato stars
O valor de cada célula do raster geotiff representa a elevação em metros do terreno, na Floresta de Harvard.
mdt_harv <- read_stars("data/HARV/HARV_dtmCrop.tif")1.2.1. Solicita descrição do objeto stars raster
# o que signfica "2 dimensions and 1 attribute"?
# o valor de cada pixel de um Modelo Digital do Terreno representa o que?
# Quantas colunas e quantas linhas tem no raster?
# Qual é a projeção ("refsys")?
mdt_harvstars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
Min. 1st Qu. Median Mean 3rd Qu. Max.
HARV_dtmCrop.tif 335.37 343.25 349.39 352.947 360.54 389.82
dimension(s):
from to offset delta refsys point x/y
x 1 1697 731453 1 WGS 84 / UTM zone 18N FALSE [x]
y 1 1367 4713838 -1 WGS 84 / UTM zone 18N FALSE [y]1.2.2. Solicite mais informação sobre o sistema de coordenadas de referência (projeção) do raster
# Na última linha do informe, qual é código EPSG da projeção?
st_crs(mdt_harv)Coordinate Reference System:
User input: WGS 84 / UTM zone 18N
wkt:
PROJCRS["WGS 84 / UTM zone 18N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
BBOX[0,-78,84,-72]],
ID["EPSG",32618]]1.3 Criar primeiro mapa com apenas duas linhas de código
ggplot() +
geom_stars(data = mdt_harv) 
- Renderize o QMD para ver o primeiro mapa na aba Viewer;
- geom_stars desenha o mapa com um mínimo de informação fornecida pelo usuário;
- Foi usada uma rampa de cores padrão (tons de azul).
1.3.1. Mude a rampa de cores para viridis_c, que detecta nuances nas elevações;
- Renderize de novo após rodar o próximo bloco.
ggplot() +
geom_stars(data = mdt_harv) +
scale_fill_viridis_c()
1.3.2. Ggplot2 usa ‘+’ ao acresentar informação incremental.
- Os incrementos ou são camadas novas, ou modificam camadas existentes.
- Muitas vezes a camada nova inicia com geom_. Exemplos comuns são geom_point(), geom_smooth(), geom_line(). Geom_sf() acrescenta camadas vetoriais.
- Scale_fill() é uma camada de coloração que preenche os pixels;
- Com a rampa viridis_c, a percepção de mudança gradual da cor acompanha a mudança linear dos valores de elevação das células.
1.4. Acrescente uma barra de escala e seta de norte
- São mais duas camadas. Cada uma segue um símbolo ‘+’
- bl e tr signficam bottom-left, e top-right,respectivamente.
- pad_x e pad_y empurram a escala e a seta para dentro da figura
ggplot() +
geom_stars(data = mdt_harv) +
scale_fill_viridis_c() +
ggspatial::annotation_scale(location = "bl",
pad_x=unit(1, "cm"), pad_y = unit(1, "cm")) +
ggspatial::annotation_north_arrow(location = "tr",
pad_x=unit(1, "cm"), pad_y = unit(1, "cm"))
1.5. Acrescente título, rótulos dos eixos, trocar rótulo da legenda, remover fundo cinza, encolher seta de norte
Uma consulta ao GPT4 e ao livro “ggplot2: Elegant Graphics for Data Analysis” indicou que basta alterar duas camadas e acrescentar outras quatro. O bloco de código fica assim:
ggplot() +
geom_stars(data = mdt_harv) +
scale_fill_viridis_c(
guide = guide_colourbar(title = "Elevação (m)")) +
theme_bw() +
labs (x ="UTM_18N Easting",
y ="UTM_18N Northing",
title = "Modelo Digital do Terreno, Harvard Forest") +
ggspatial::annotation_scale(location = "bl",
pad_x=unit(1, "cm"),
pad_y = unit(1, "cm")) +
ggspatial::annotation_north_arrow(location = "tr",
height = unit(1.0, "cm"), # encolhe a seta
width = unit(1.0, "cm"), # encolhe a seta
pad_x=unit(1, "cm"),
pad_y = unit(1, "cm"))
1.6. Salvar o mapa final como arquivo de imagem em disco
Para isso, atribuimos o mapa ao objeto ‘primeiro_mapa’ e rodamos ggsave.
primeiro_mapa <- ggplot() +
geom_stars(data = mdt_harv) +
scale_fill_viridis_c(guide = guide_colourbar(title = "Elevação (m)")) +
theme_bw() +
labs (x ="UTM_18N Easting",
y ="UTM_18N Northing",
title = "Modelo Digital do Terreno, Harvard Forest") +
ggspatial::annotation_scale(location = "bl",
pad_x=unit(1, "cm"),
pad_y = unit(1, "cm")) +
ggspatial::annotation_north_arrow(location = "tr",
height = unit(1.0, "cm"), # encolhe a seta
width = unit(1.0, "cm"), # encolhe a seta
pad_x=unit(1, "cm"),
pad_y = unit(1, "cm"))
ggsave(filename = "mapa_MDT_harvard_forest.png",
plot = primeiro_mapa,
width = 7,
height = 5,
dpi = 300,
bg = "white")Localize o arquivo png na pasta do projeto. Clicar duas vezes para visualizar. Caso necessário, altere ‘width’ e ‘height’ no código de ggsave() para obter escalas iguais nos sentidos X e Y, sem distorção.
FIM.
AULA 2. ÁLGEBRA COM RASTERS: IMAGEM-DIFERENÇA
2.1. Rodar todos os blocos de código acima deste ponto e renderizar
2.2. Importar dois modelos digitais de elevações
- modelo digital do terreno e
- modelo digital da superfície superior do dossel
mdt_harv <- read_stars("data/HARV/HARV_dtmCrop.tif")
mds_harv <- read_stars("data/HARV/HARV_dsmCrop.tif")2.3. Antes de realizar a subtração entre imagens:
2.3.1. Verifique se os dois rasters casem em termos de Sistema de Coordenadas de Referência (projeção), tamanho, resolução espacial e retângulo envolvente.
Selecione cada linha decomando dentro do bloco de código abaixo e rodar linha por linha usando Ctrl-Enter. O output de interesse vai aparecer no Console do R e no documento renderizado
# numero de colunas e linhas
dim(mds_harv) x y
1697 1367 dim(mdt_harv) x y
1697 1367 # retângulo envolvente, em unidades projetadas
st_bbox(mds_harv) xmin ymin xmax ymax
731453 4712471 733150 4713838 st_bbox(mdt_harv) xmin ymin xmax ymax
731453 4712471 733150 4713838 # Extrair EPSG de mds_harv
epsg_mds_harv <- st_crs(mds_harv)$epsg
print(paste("EPSG for mds_harv:", epsg_mds_harv))[1] "EPSG for mds_harv: 32618"# Extrair EPSG de mdt_harv
epsg_mdt_harv <- st_crs(mdt_harv)$epsg
print(paste("EPSG for mdt_harv:", epsg_mdt_harv))[1] "EPSG for mdt_harv: 32618"2.3.2 As duas imagens casam perfeitamente.
2.3.3. Qual é a resolução espacial?
Veja como calcular a resolução em metros, na direção horizontal (cell_size_x), com este código fornecido sugerido pelo GPT4
#Computar o retangulo envolvente
bbox <- st_bbox(mds_harv)
# Extrair xmax e xmin em metros (coordenadas UTM projetadas) como valores numéricos, removendo seus nomes
xmax_value <- as.numeric(bbox["xmax"])
xmin_value <- as.numeric(bbox["xmin"])
# Obter a dimensão X da imagem, em metros (coordenadas UTM)
x_dist_utm <- xmax_value - xmin_value
# Obter as dimensões em pixels (colunas e linhas)
dims <- dim(mds_harv)
# Extrair número de colunas
num_cols <- dims["x"]
# Dividir para obter resolução em metros
cell_size_x <- (x_dist_utm / num_cols)
message <- paste("Resolução em metros no sentido x é ", unname(cell_size_x))
print(message)[1] "Resolução em metros no sentido x é 1"2.4. Criar a imagem diferença ‘mav_harv’
A diferença entre a imagem contendo as elevações do topo do dossel e a imagem contendo as elevações do terreno representa a altura da vegetação em relação ao terreno.
# Modelo de Altura da Vegetação
mav_harv <- (mds_harv - mdt_harv)
ggplot()+
geom_stars(data = mav_harv)
2.5. Exercício prático, Mapa de Altura da Vegetação
(1). Usando um código que inicia com ggplot( ), visualize a imagem do Modelo de Altura da Vegetação. Renderize-o para que apareca na aba Viewer como a parte final da Aula 2.
(2). Siga o exemplo da Aula 1, inclua os seguintes elementos:
- Escala de cores viridis_c
- Tema que remove o fundo cinza
- Título do mapa
- Eixos rotulados
- Título “Altura da Vegetação (m)” para a legenda de cores
- Barra de escala em metros
- Seta para norte
(3). Cria outro bloco com o mesmo conteúdo, mas atribuindo o resultado a um objeto, como feito no final da Aula 1.
(4) Use ggsave() para criar uma cópia do mapa em disco, como arquivo de imagem ‘altura_floresta_harvard.png’