rm(list = ls())
options(scipen = 9999)
mapview::mapviewOptions(basemaps = 'Esri.WorldImagery', platform = 'leaflet')PROCESSAMENTO E EXTRAÇÃO DE MÉTRICAS LiDAR
FAZENDA APARÁ
1 PARTE
1.1 PRÉ-CONFIGURAÇÕES
Inicialmente, precisamos limpar a memória do computador com a função rm() do R. Além disso, definiremos uma opção para evitar a exibição de valores altos em notação científica com options(scipen = 9999). Essas configurações são úteis, pois trabalharemos com dados com grandes quantidades de observações.
1.2 PACOTES
Nesta etapa, vamos carregar o pacote pacman(Rinker and Kurkiewicz 2018) para gerenciar os pacotes que usaremos durante o trabalho. A vantagem de usar pacman é que sua função p_load() carrega os pacotes especificados, se já estiverem instalados; caso contrário, ele os instala primeiro e depois os carrega. Se o pacote pacman não estiver instalado, você pode usar a função install.packages('pacman') para instalá-lo. A linha correspondente para instalação está comentada no bloco abaixo com um # no início.
#install.packages('pacman')
library(pacman)
pacman::p_load(lidR, tidyverse, tmap, sf, terra, mapview)1.3 LEITURA INICIAL DOS DADOS
Realizaremos a leitura inicial dos dados com a função readLAS() do pacote lidR. (Roussel et al. 2020a).
las <- readLAS('Fazenda Apara/lidars/terra_las/cloud31e1d8fb1c24cc10.las')
lasclass : LAS (v1.2 format 3)
memory : 139.7 Mb
extent : -393421.7, -393271.4, 9877640, 9877815 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 24S
area : 21836 m²
points : 2.44 million points
density : 111.84 points/m²
density : 107.52 pulses/m²O arquivo ocupa um tamanho de 139,7 MB na memória. O voo foi realizado em uma área de 21.836 m², gerando 2,44 milhões de pontos. Com uma densidade de 107,52 pulsos por m², foi possível detectar 111,84 pulsos/m².
Outra informação importante é o sistema de projeção de coordenadas, que está em WGS 84 / UTM zone 24S.
1.4 PRÉ-VISUALIZAÇÃO
A pré-visualização é importante para avaliar a presença de outliers nos dados, como pássaros ou objetos voadores não identificados.
Para isso, usaremos a função plot(), que já vem pré-instalada com o R. Acrescentaremos alguns parâmetros, como bg = 'white' para deixar o plano de fundo branco, axis = TRUE para a plotagem dos eixos cartesianos, legend = TRUE para gerar uma legenda, e color = "Z" para graduar os dados com base na altura dos pontos.
plot(las, bg = 'white', axis = T, legend = T, color = "Z")
Nesse caso, os dados não apresentam nenhum outlier aparente, portanto, não será necessário realizar nenhuma filtragem.
1.5 AVALIAÇÃO DOS ATRIBUTOS
Cada coluna da tabela representa uma característica do ponto. Observamos que, com este voo, foram detectados 2.442.038 pontos e, para cada ponto, há 9 atributos.
las@data X Y Z gpstime Intensity ReturnNumber
1: -393317.2 9877811 17.09978 396444786 33 1
2: -393317.8 9877811 16.68088 396444786 21 1
3: -393318.2 9877811 15.39658 396444786 31 1
4: -393319.0 9877811 15.19328 396444786 23 1
5: -393319.7 9877811 14.86018 396444786 17 1
---
2442034: -393326.4 9877650 -14.74972 396444801 0 3
2442035: -393344.1 9877649 -15.89532 396444801 0 3
2442036: -393345.8 9877649 -17.63802 396444801 0 2
2442037: -393289.3 9877648 -10.71892 396444802 0 3
2442038: -393289.3 9877648 -10.60792 396444802 0 3
NumberOfReturns ScanDirectionFlag EdgeOfFlightline Classification
1: 1 0 0 0
2: 1 0 0 0
3: 1 0 0 0
4: 1 0 0 0
5: 1 0 0 0
---
2442034: 3 0 0 0
2442035: 3 0 0 0
2442036: 2 0 0 0
2442037: 3 0 0 0
2442038: 3 0 0 0
Synthetic_flag Keypoint_flag Withheld_flag ScanAngleRank UserData
1: FALSE FALSE FALSE -20 0
2: FALSE FALSE FALSE -19 0
3: FALSE FALSE FALSE -19 0
4: FALSE FALSE FALSE -18 0
5: FALSE FALSE FALSE -18 0
---
2442034: FALSE FALSE FALSE -16 0
2442035: FALSE FALSE FALSE -16 0
2442036: FALSE FALSE FALSE -16 0
2442037: FALSE FALSE FALSE -16 0
2442038: FALSE FALSE FALSE -16 0
PointSourceID R G B
1: 0 25088 31488 20736
2: 0 34816 41472 29440
3: 0 21248 28672 16384
4: 0 39168 48384 29440
5: 0 15360 19968 9216
---
2442034: 0 33536 41216 25856
2442035: 0 43520 50688 24576
2442036: 0 33024 41728 22784
2442037: 0 34048 41984 26624
2442038: 0 43776 54272 389121.6 SELEÇÃO DOS ATRIBUTOS DOS DADOS
Grande parte dos atributos do voo não será utilizada; portanto, vamos selecionar apenas os atributos importantes para nossa análise, que são as coordenadas dos pontos (xyz), intensidade de retorno (i), número de retorno (r), número de retornos (n) e os parâmetros RGB (RGB). Para isso, utilizaremos o parâmetro select() da função readLAS().
las <- readLAS('Fazenda Apara/lidars/terra_las/cloud31e1d8fb1c24cc10.las', select = 'xyzirnRGB')
lasclass : LAS (v1.2 format 3)
memory : 111.8 Mb
extent : -393421.7, -393271.4, 9877640, 9877815 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 24S
area : 21836 m²
points : 2.44 million points
density : 111.84 points/m²
density : 107.52 pulses/m²plot(las, mapview = T)1.7 ALTERANDO O SISTEMA DE PROJEÇÃO
Por um motivo desconhecido, o sistema de coordenadas veio incorreto no momento da exportação dos dados em formato .las. Para corrigir esse erro, vamos reprojetar os dados para o sistema SIRGAS 2000 / UTM zona 23S, cujo código EPSG é 31983. Faremos isso com a função st_transform() do pacote sf (Pebesma and Bivand 2023). Em seguida, plotaremos novamente os dados para reavaliar.
las <- st_transform(las, crs = 31983)
plot(las, mapview = T)lasclass : LAS (v1.2 format 3)
memory : 111.8 Mb
extent : 276989.5, 277138.4, 9878765, 9878938 (xmin, xmax, ymin, ymax)
coord. ref. : SIRGAS 2000 / UTM zone 23S
area : 21440 m²
points : 2.44 million points
density : 113.9 points/m²
density : 109.51 pulses/m²Agora o dado está no sistema de projeção correto.
1.8 CHECAGEM DOS DADOS
Nesta etapa, realizaremos uma checagem para avaliar se os dados estão livres de imperfeições.
las_check(las)
1.9 CLASSIFICAÇÃO
Nesta parte, classificaremos os dados em duas classes (solo e vegetação) usando a função classify_ground() do pacote lidR com o algoritmo csf(), desenvolvido por (Roussel and Qi 2020).
Após a classificação, podemos filtrar apenas os dados de solo usando a função filter_ground(), também do pacote lidR. Para melhorar a classificação, excluiremos os dados com elevação maior que -17, pois foram visualmente identificados como não pertencentes ao solo.
las <- classify_ground(las = las, algorithm = csf())
solo <- filter_ground(las) %>% filter_poi(Z <= -17)podemos visualizar o resultado com a função plot().
plot(solo, bg = 'white', axis = T, legend = T)
Poucos pontos foram detectados como solo. Três hipóteses foram levantadas para explicar essa situação:
A alta densidade da vegetação dificultou a penetração dos pulsos do LiDAR.
A área estava muito alagada, o que impediu o retorno dos pulsos que atravessaram a vegetação.
Eu sou muito ruim de programação e não consegui adaptar os parâmetros para uma boa classificação. (acho que é isso)
1.10 GERAR MDT
Para gerar o modelo digital do terreno, foi desenvolvida uma aplicação do algoritmo Multilevel B-spline Approximation (MBA) (Finley, Banerjee, and Hjelle 2022). Esse algoritmo busca aproximar ao máximo o terreno gerado, mantendo a suavidade da elevação.
mba <- function(n = 1, m = 1, h = 8, extend = TRUE) {
f <- function(gnd, where) {
res <- MBA::mba.points(gnd@data, where, n, m , h, extend)
return(res$xyz.est[,3])
}
f <- plugin_dtm(f)
return(f)
}Para aplicar esse algoritmo, utilizamos a função rasterize_terrain() do pacote lidR na nossa variável solo.
mdt <- rasterize_terrain(las = solo,algorithm = mba())Podemos visualizar o resultado com a função Plot_dtm3d();
plot_dtm3d(mdt)
Se for necessário exportar esse modelo como um raster para utilizá-lo no QGIS, por exemplo, para criar um layout personalizado ou aplicá-lo em alguma modelagem, podemos usar a função writeLAS() do pacote terra (Hijmans 2023). Estou utilizando o parâmetro overwrite() para reescrever o arquivo caso já exista um com o mesmo nome no diretório.
terra::writeRaster(x = mdt, filename = 'produtosR/MDT.tif', overwrite = T)1.11 NORMALIZAR A NUVEM
A normalização é necessária antes da extração das métricas, pois remove a interferência do terreno sobre as características morfométricas da vegetação. Isso faz com que as métricas representem com mais precisão a realidade.
Para realizar essa etapa, usaremos a função normalize_height(), mais uma vez do pacote lidR, com o algoritmo tin() (Turner 2023).
las <- normalize_height(las = las, algorithm = tin())1.12 EXPORTAÇÃO DA NUVEM DE PONTOS PROCESSADA
Dependendo do tamanho dos dados, esses processos podem ser extremamente lentos. Por isso, é importante exportar a nova nuvem de pontos já classificada e normalizada para usar nas próximas etapas. Para isso, utilizaremos a função writeLAS() do pacote lidR.
writeLAS(las = las, file = 'produtosR/nuvem_normalizada.las')2 PARTE
Para iniciar a segunda parte, precisaremos esvaziar a memória do computador e pré-configurar novamente o ambiente R, como foi feito no início da parte 1. Também será necessário instalar os pacotes necessários.
rm(list = ls())
options(scipen = 9999)library(pacman)
pacman::p_load(lidR, tidyverse, tmap, sf, terra)Agora podemos dar início a nova fase.
2.1 LENDO A NUVEM NORMALIZADA
Vamos realizar a leitura da nuvem de pontos processada que exportamos ao final da parte 1. Em seguida, faremos uma plotagem para visualização prévia.
las <- readLAS('produtosR/nuvem_normalizada.las', select = 'xyzirncRGB')
plot(las, mapview = T)2.2 GERAR O CHM
Para gerar o modelo digital do dossel, utilizaremos a função grid_canopy() do pacote lidR, com os parâmetros adicionais pitfree() para remoção de ruídos e subcircle() para simular a pegada do pulso e proporcionar mais realismo ao processo físico, conforme descrito por (Roussel et al. 2020b).
chm <- grid_canopy(las = las, res = 0.5, pitfree(c(0,20,40,60,80), c(0,1), subcircle = 0.2))Visualização do chm em escala de cor viridis().
plot(chm, col = viridis::viridis(15),
main = 'Modelo digital de superficie',
xlab = 'Coordenadas X',
ylab = 'Coordenadas Y')
Caso seja necessário exportar esse raster para uso em um software externo, podemos utilizar:
terra::writeRaster(x = chm, filename = 'produtosR/chm.tif', overwrite = T)2.3 AVALIANDO AS EXTENÇÕES DO VOO COM A ÁREA PROJETADA PREVIAMENTE
Previamente, estabelecemos uma área de sobrevoo (area1.gpkg). Vamos ler essa área no R usando a função st_read() do pacote sf e salvaremos com o nome roi (região de interesse).
roi <- st_read('Fazenda Apara/areas de voo/area1.gpkg')Reading layer `area_florestal__area_florestalkml' from data source
`E:\Geoprocessamento\LIDAR\fazenda APARA\Fazenda Apara\areas de voo\area1.gpkg'
using driver `GPKG'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension: XYZ
Bounding box: xmin: 276929.9 ymin: 9878583 xmax: 277545.2 ymax: 9879121
z_range: zmin: 0 zmax: 0
Projected CRS: SIRGAS 2000 / UTM zone 23SAgora, vamos plotar a região de interesse (ROI) com a área e o modelo digital do dossel sobrepostos.
plot(roi$geom, axes = T)
plot(chm, col = height.colors(10),
main = 'Modelo digital de superficie',
xlab = 'Coordenadas X',
ylab = 'Coordenadas Y',
legend = F,
add = T)
Podemos observar que a área sobrevoada foi muito menor do que o planejado, ou o voo foi dividido em várias parcelas e nos foi disponibilizada apenas uma parte.
2.4 EXTRAINDO AS MÉTRICAS A NÍVEL DE PIXEL
O ideal seria que a área sobrevoada fosse maior do que a área de interesse, para que as bordas pudessem ser recortadas, uma vez que a falta de contexto espacial pode prejudicar a extração das métricas.
No entanto, vamos supor que isso não irá influenciar o resultado e proceder com a geração das métricas.
Para isso, utilizaremos a função pixel_metrics() do pacote lidR com o parâmetro .stdmetrics (que calcula automaticamente 56 métricas padrão) e uma resolução de 10m, apenas a título de exemplo.
metrics <- pixel_metrics(las, .stdmetrics, res = 10)E vamos ver os resultados:
Da métrica 1 à 16;
plot(metrics[[1:16]], col = height.colors(50))
Da métrica 17 à 33;
plot(metrics[[17:33]], col = height.colors(50))
Da métrica 34 à 50;
plot(metrics[[34:50]], col = height.colors(50))
Da métrica 51 à 56;
plot(metrics[[51:56]], col = height.colors(50))
References
Finley, Andrew, Sudipto Banerjee, and Øyvind Hjelle. 2022. “MBA: Multilevel b-Spline Approximation.” https://CRAN.R-project.org/package=MBA.
Hijmans, Robert J. 2023. “Terra: Spatial Data Analysis.” https://CRAN.R-project.org/package=terra.
Pebesma, Edzer, and Roger Bivand. 2023. “Spatial Data Science: With Applications in r.” https://doi.org/10.1201/9780429459016.
Rinker, Tyler W., and Dason Kurkiewicz. 2018. “Pacman: Package Management for r.” http://github.com/trinker/pacman.
Roussel, Jean-Romain, David Auty, Nicholas C. Coops, Piotr Tompalski, Tristan R. H. Goodbody, Andrew Sánchez Meador, Jean-François Bourdon, Florian de Boissieu, and Alexis Achim. 2020a. “lidR: An r Package for Analysis of Airborne Laser Scanning (ALS) Data” 251: 112061. https://doi.org/10.1016/j.rse.2020.112061.
———. 2020b. “lidR: An R Package for Analysis of Airborne Laser Scanning (ALS) Data.” Remote Sensing of Environment 251 (December): 112061. https://doi.org/10.1016/j.rse.2020.112061.
Roussel, Jean-Romain, and Jianbo Qi. 2020. “RCSF: Airborne LiDAR Filtering Method Based on Cloth Simulation.” https://CRAN.R-project.org/package=RCSF.
Turner, Rolf. 2023. “Deldir: Delaunay Triangulation and Dirichlet (Voronoi) Tessellation.” https://CRAN.R-project.org/package=deldir.