# Instalar o pacote
#install.packages('lidR')
# Carregando o pacote
library(lidR)lidR e suas funcionalidades
1 INTRODUÇÃO
lidRé um pacote R para manipulação e visualização de dados de varredura a laser aerotransportada (ALS) com ênfase em aplicações florestais.
O pacote é totalmente de código aberto e está integrado ao ecossistema R geoespacial (ou seja, raster/terra/stars e sp/sf).
A principal funcionalidade lidR inclui funções para:
Ler, gravar e renderizar exibição de nuvem de pontos personalizada (
.lase.laz);Processar nuvens de pontos, incluindo classificação de pontos, modelos digitais de terreno(MDT), normalização e modelos digitais de superfície(MDS);
Executa a segmentação de árvore individual;
Calcula métricas padrão em diferentes níveis de regularização;
Gerenciar o processamento de conjuntos de arquivos de nuvem de pontos - referidos como
LAScatalog;Facilitar fluxos de processamento definidos pelo usuário para pesquisa e desenvolvimento.
2 LENDO DADOS LIDAR
Os sensores ALS de retorno discreto registram vários dados. Em primeiro lugar, dados posicionais em três dimensões (X,Y,Z), seguidos de informações adicionais como a intensidade de cada ponto, a posição de cada ponto na sequência de retorno ou o ângulo de incidência do feixe de cada ponto.
Dados de varedura a laser aerotransportados são armazenados no formato .las , mas por serem muito pesados utlizam-se o formato de compactação .laz.
A função readLAS() executa a leitura tanto de arquivo .las como .laz e retorna um objeto da classe .las dentro do R.
Um arquivo LAS é composto de duas partes:
O cabeçalho que armazena informações resumidas sobre seu conteúdo, o sistema de referência de coordenadas e o formato do ponto.
A carga útil – ou seja, a própria nuvem de pontos.
Para mais informações sobre o formato
.lasacesse: https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAS-class.html
Executando a leitura com a função readLAS()
las <- readLAS('data/MixedEucaNat.laz')
# Inspecionar informações do cabeçalho
las@headerFile signature: LASF
File source ID: 0
Global encoding:
- GPS Time Type: GPS Week Time
- Synthetic Return Numbers: no
- Well Know Text: CRS is GeoTIFF
- Aggregate Model: false
Project ID - GUID: 00000000-0000-0000-0000-000000000000
Version: 1.2
System identifier:
Generating software: rlas R package
File creation d/y: 0/2013
header size: 227
Offset to point data: 297
Num. var. length record: 1
Point data format: 0
Point data record length: 20
Num. of point records: 553544
Num. of points by return: 403345 126738 21779 1637 45
Scale factor X Y Z: 0.01 0.01 0.01
Offset X Y Z: 2e+05 7300000 0
min X Y Z: 203830 7358900 731.83
max X Y Z: 203980 7359050 806.34
Variable Length Records (VLR):
Variable Length Record 1 of 1
Description: by LAStools of rapidlasso GmbH
Tags:
Key 3072 value 31983
Extended Variable Length Records (EVLR): void# Inspecionar atributos da nuvem de pontos
las@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 203851.6 7359049 766.84 285 1 1
2: 203922.2 7359048 768.11 343 1 1
3: 203942.9 7359045 770.18 104 2 2
4: 203830.0 7359045 765.83 284 1 1
5: 203841.2 7359047 766.29 290 1 1
---
553540: 203902.5 7359050 766.85 259 2 2
553541: 203907.1 7359050 766.81 206 1 1
553542: 203956.0 7359050 771.03 309 1 1
553543: 203962.5 7359050 771.28 100 2 2
553544: 203972.6 7359050 771.64 46 2 2
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 2 FALSE
2: 0 0 2 FALSE
3: 0 0 2 FALSE
4: 0 0 2 FALSE
5: 0 0 2 FALSE
---
553540: 0 0 2 FALSE
553541: 0 0 2 FALSE
553542: 0 0 2 FALSE
553543: 0 0 2 FALSE
553544: 0 0 2 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID
1: FALSE TRUE -21 0 14
2: FALSE TRUE -21 0 14
3: FALSE TRUE -21 0 14
4: FALSE TRUE -21 0 14
5: FALSE TRUE -21 0 14
---
553540: FALSE TRUE -3 0 15
553541: FALSE TRUE -3 0 15
553542: FALSE TRUE -21 0 14
553543: FALSE TRUE -21 0 14
553544: FALSE TRUE -1 0 15# Verificar o tamanho do arquivo dos dados LiDAR carregados
format(object.size(las), "Mb")[1] "38 Mb"2.1 Parâmetro de seleção
Um arquivo LAS armazena as coordenadas X Y Z de cada ponto, bem como muitos outros dados, como intensidade, ângulo de incidência e posição da sequência de retorno.
Na prática, muitos atributos não são realmente úteis, mas são carregados de qualquer maneira por padrão, isso pode ocupar muita memória de processamento.
Para economizar memória, readLAS()pode-se levar um parâmetro opcional select que permite ao usuário carregar seletivamente os atributos de interesse.
las <- readLAS("data/MixedEucaNat.laz", select = "xyz") # lê somente XYZ
las <- readLAS("data/MixedEucaNat.laz", select = "xyzinr") # lê somente XYZ, intensity, número retornado e número de retorno
lasclass : LAS (v1.2 format 0)
memory : 19 Mb
extent : 203830, 203980, 7358900, 7359050 (xmin, xmax, ymin, ymax)
coord. ref. : SIRGAS 2000 / UTM zone 23S
area : 22500 m²
points : 553.5 thousand points
density : 24.6 points/m²
density : 17.93 pulses/m²format(object.size(las), 'Mb')[1] "19 Mb"las@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 203851.6 7359049 766.84 285 1 1
2: 203922.2 7359048 768.11 343 1 1
3: 203942.9 7359045 770.18 104 2 2
4: 203830.0 7359045 765.83 284 1 1
5: 203841.2 7359047 766.29 290 1 1
---
553540: 203902.5 7359050 766.85 259 2 2
553541: 203907.1 7359050 766.81 206 1 1
553542: 203956.0 7359050 771.03 309 1 1
553543: 203962.5 7359050 771.28 100 2 2
553544: 203972.6 7359050 771.64 46 2 2Exemplos de outras abreviações de atributos são: t- gpstime, a- ângulo de varredura, n- número retornado, r- número de retorno, c- classificação, s- sinalizador sintético, k- sinalizador de ponto-chave, w- sinalizador retido, o- sinalizador de sobreposição (formato 6+), u- dados do usuário , p- ID da fonte do ponto, e- borda da bandeira da linha de voo, d- direção da bandeira de varredura.
2.2 Filtragem dos dados
filter permite a seleção de “linhas” (ou pontos) durante a leitura. A remoção de dados supérfluos no momento da leitura economiza memória e aumenta a velocidade de computação. Por exemplo, é prática comum na silvicultura processar usando os primeiros retornos.
las <- readLAS("data/MixedEucaNat.laz", filter = "-keep_first") # Leitura somente dos primeiros retornos
las@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 203851.6 7359049 766.84 285 1 1
2: 203922.2 7359048 768.11 343 1 1
3: 203830.0 7359045 765.83 284 1 1
4: 203841.2 7359047 766.29 290 1 1
5: 203851.5 7359049 766.87 242 1 1
---
403341: 203854.9 7359050 767.48 276 1 1
403342: 203868.2 7359050 768.39 340 1 1
403343: 203881.0 7359049 768.35 401 1 1
403344: 203907.1 7359050 766.81 206 1 1
403345: 203956.0 7359050 771.03 309 1 1
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 2 FALSE
2: 0 0 2 FALSE
3: 0 0 2 FALSE
4: 0 0 2 FALSE
5: 0 0 2 FALSE
---
403341: 0 0 2 FALSE
403342: 0 0 2 FALSE
403343: 0 0 1 FALSE
403344: 0 0 2 FALSE
403345: 0 0 2 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID
1: FALSE TRUE -21 0 14
2: FALSE TRUE -21 0 14
3: FALSE TRUE -21 0 14
4: FALSE TRUE -21 0 14
5: FALSE TRUE -21 0 14
---
403341: FALSE TRUE -21 0 14
403342: FALSE TRUE -3 0 15
403343: FALSE TRUE -3 0 15
403344: FALSE TRUE -3 0 15
403345: FALSE TRUE -21 0 14Você pode usar tanto o parâmetro filter usado acima, ou pode usar as técnicas de filtragem comuns no R, como usando a função filter_poi(), a seguir:
las2 <- readLAS("data/MixedEucaNat.laz")
las2 <- filter_poi(las2, ReturnNumber == 1L)
las2@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 203851.6 7359049 766.84 285 1 1
2: 203922.2 7359048 768.11 343 1 1
3: 203830.0 7359045 765.83 284 1 1
4: 203841.2 7359047 766.29 290 1 1
5: 203851.5 7359049 766.87 242 1 1
---
403341: 203854.9 7359050 767.48 276 1 1
403342: 203868.2 7359050 768.39 340 1 1
403343: 203881.0 7359049 768.35 401 1 1
403344: 203907.1 7359050 766.81 206 1 1
403345: 203956.0 7359050 771.03 309 1 1
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 2 FALSE
2: 0 0 2 FALSE
3: 0 0 2 FALSE
4: 0 0 2 FALSE
5: 0 0 2 FALSE
---
403341: 0 0 2 FALSE
403342: 0 0 2 FALSE
403343: 0 0 1 FALSE
403344: 0 0 2 FALSE
403345: 0 0 2 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID
1: FALSE TRUE -21 0 14
2: FALSE TRUE -21 0 14
3: FALSE TRUE -21 0 14
4: FALSE TRUE -21 0 14
5: FALSE TRUE -21 0 14
---
403341: FALSE TRUE -21 0 14
403342: FALSE TRUE -3 0 15
403343: FALSE TRUE -3 0 15
403344: FALSE TRUE -3 0 15
403345: FALSE TRUE -21 0 14Nos exemplos acima (1) lendo apenas os primeiros retornos ou (2) lendo todos os pontos e depois filtrando os primeiros retornos em R. Ambas as saídas são estritamente idênticas, mas a primeira é mais rápida e mais eficiente em termos de memória porque não carrega o arquivo inteiro em R e não usa memória de processamento extra. Vários comandos de filtro podem ser usados ao mesmo tempo para, por exemplo, primeiro retorno somente leitura entre 5 e 50 m.
las <- readLAS("data/MixedEucaNat_normalized.laz", filter = "-keep_first -drop_z_below 5 -drop_z_above 50")
las@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 203830.7 7358900 6.63 284 1 1
2: 203830.0 7358900 6.78 317 1 1
3: 203830.4 7358901 6.27 256 1 1
4: 203830.4 7358900 6.82 239 1 1
5: 203830.7 7358900 5.66 191 1 2
---
238244: 203979.2 7359049 9.27 66 1 2
238245: 203978.6 7359050 16.74 56 1 2
238246: 203978.7 7359050 5.41 143 1 2
238247: 203979.8 7359049 5.82 105 1 2
238248: 203979.9 7359050 18.23 90 1 3
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 1 FALSE
2: 0 0 1 FALSE
3: 0 0 1 FALSE
4: 0 0 1 FALSE
5: 0 0 1 FALSE
---
238244: 0 0 1 FALSE
238245: 0 0 1 FALSE
238246: 0 0 1 FALSE
238247: 0 0 1 FALSE
238248: 0 0 1 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID
1: FALSE TRUE 6 0 16
2: FALSE TRUE 6 0 16
3: FALSE TRUE 6 0 16
4: FALSE TRUE 6 0 16
5: FALSE TRUE 6 0 16
---
238244: FALSE TRUE -1 0 15
238245: FALSE TRUE -1 0 15
238246: FALSE TRUE -21 0 14
238247: FALSE TRUE -1 0 15
238248: FALSE TRUE -1 0 15A lista completa de comandos disponíveis é fornecida por readLAS(filter = "-help")
readLAS(filter = "-help")3 VALIDAÇÃO DOS DADOS
Um primeiro passo importante no processamento de dados do ALS é garantir que seus dados estejam completos e válidos de acordo com as especificações.
Um exemplo simples que acontece com bastante frequência é que um LAS arquivo contendo pontos duplicados. Isso pode levar a problemas como árvores sendo detectadas duas vezes, métricas inválidas ou erros na geração do MDT. Também podemos encontrar números de retorno inválidos, números de retorno incoerentes, sistema de referência de coordenadas inválido, etc. Sempre certifique-se de executar a função las_check() antes de se aprofundar em seus dados.
las <- readLAS('data/MixedEucaNat.laz')
las_check(las)
Checking the data
- Checking coordinates...[0;32m ✓[0m
- Checking coordinates type...[0;32m ✓[0m
- Checking coordinates range...[0;32m ✓[0m
- Checking coordinates quantization...[0;32m ✓[0m
- Checking attributes type...[0;32m ✓[0m
- Checking ReturnNumber validity...[0;32m ✓[0m
- Checking NumberOfReturns validity...[0;32m ✓[0m
- Checking ReturnNumber vs. NumberOfReturns...[0;32m ✓[0m
- Checking RGB validity...[0;32m ✓[0m
- Checking absence of NAs...[0;32m ✓[0m
- Checking duplicated points...[0;32m ✓[0m
- Checking degenerated ground points...
[1;33m ⚠ There were 37 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates[0m
- Checking attribute population...
[0;32m 🛈 'ScanDirectionFlag' attribute is not populated[0m
[0;32m 🛈 'EdgeOfFlightline' attribute is not populated[0m
- Checking gpstime incoherances[0;37m skipped[0m
- Checking flag attributes...
[0;32m 🛈 127974 points flagged 'withheld'[0m
- Checking user data attribute...[0;32m ✓[0m
Checking the header
- Checking header completeness...[0;32m ✓[0m
- Checking scale factor validity...[0;32m ✓[0m
- Checking point data format ID validity...[0;32m ✓[0m
- Checking extra bytes attributes validity...[0;32m ✓[0m
- Checking the bounding box validity...[0;32m ✓[0m
- Checking coordinate reference system...[0;32m ✓[0m
Checking header vs data adequacy
- Checking attributes vs. point format...[0;32m ✓[0m
- Checking header bbox vs. actual content...[0;32m ✓[0m
- Checking header number of points vs. actual content...[0;32m ✓[0m
- Checking header return number vs. actual content...[0;32m ✓[0m
Checking coordinate reference system...
- Checking if the CRS was understood by R...[0;32m ✓[0m
Checking preprocessing already done
- Checking ground classification...[0;32m yes[0m
- Checking normalization...[0;31m no[0m
- Checking negative outliers...[0;32m ✓[0m
- Checking flightline classification...[0;32m yes[0m
Checking compression
- Checking attribute compression...
- ScanDirectionFlag is compressed
- EdgeOfFlightline is compressed
- Synthetic_flag is compressed
- Keypoint_flag is compressed
- UserData is compressed4 PLOTAGEM DOS DADOS
O pacote lidR aproveita o pacote rgl para fornecer um visualizador 3D versátil e interativa com pontos coloridos por coordenadas Z em um fundo preto como padrão.
4.1 Renderização 3D básica
A maneira básica de renderizar uma nuvem de pontos é a função plot().
plot(las)4.2 Renderização estética
Os usuários podem alterar os atributos usados para colorir fornecendo o nome do atributo usado para colorir os pontos. A cor de fundo do visualizador também pode ser alterada atribuindo uma cor usando o bgargumento. Os eixos também podem ser adicionados e os tamanhos dos pontos podem ser alterados e a escala da legenda também.
plot(las, color = "Z", bg = "white", axis = TRUE, legend = TRUE)
plot(las, color = "Z", breaks = "quantile", bg = "white", axis = T, legend = T, size = 2)Observe que se o seu arquivo contiver dados RGB, a string "RGB"será suportada:
lasRGB <- readLAS('data/exemple_rgb.las')
lasRGB@data X Y Z gpstime Intensity ReturnNumber NumberOfReturns
1: 350825.8 5742066 1.890 0 0 0 0
2: 350825.8 5742066 1.956 0 0 0 0
3: 350825.8 5742066 1.921 0 0 0 0
4: 350825.8 5742066 1.949 0 0 0 0
5: 350825.8 5742066 1.989 0 0 0 0
---
545359: 350902.3 5742102 0.386 0 0 0 0
545360: 350902.3 5742102 0.382 0 0 0 0
545361: 350902.3 5742102 0.192 0 0 0 0
545362: 350902.3 5742102 0.201 0 0 0 0
545363: 350902.4 5742102 0.461 0 0 0 0
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 1 FALSE
2: 0 0 1 FALSE
3: 0 0 1 FALSE
4: 0 0 1 FALSE
5: 0 0 1 FALSE
---
545359: 0 0 1 FALSE
545360: 0 0 1 FALSE
545361: 0 0 1 FALSE
545362: 0 0 1 FALSE
545363: 0 0 1 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID R
1: FALSE FALSE 0 0 0 52224
2: FALSE FALSE 0 0 0 34304
3: FALSE FALSE 0 0 0 41216
4: FALSE FALSE 0 0 0 38656
5: FALSE FALSE 0 0 0 33536
---
545359: FALSE FALSE 0 0 0 18688
545360: FALSE FALSE 0 0 0 23296
545361: FALSE FALSE 0 0 0 18176
545362: FALSE FALSE 0 0 0 19200
545363: FALSE FALSE 0 0 0 17408
G B
1: 50688 46336
2: 33536 22528
3: 38912 33536
4: 38144 27648
5: 33536 21504
---
545359: 19200 23552
545360: 23296 26368
545361: 19200 23040
545362: 20224 24576
545363: 18944 23296Realizando a plotagem do arquivo que contém dados RGB.
#plot(lasRGB, color = "RGB", bg = "white", axis = TRUE, legend = TRUE)4.3 Renderização de voxel
vox <- voxelize_points(las,res = 6)
plot(vox, voxel = TRUE, bg = "white", axis = TRUE, legend = TRUE)4.4 Renderização 2D de secções transversais
Para melhor visualizar a estrutura vertical de uma nuvem de pontos, investigar os resultados da classificação ou comparar os resultados de diferentes rotinas de interpolação, é possível traçar uma seção transversal. Para realizar esse procedimento, primeiro precisamos determinar onde a seção transversal será localizada, ou seja, definir o ponto de início e o ponto de fim, além de especificar a sua largura. Em seguida, a nuvem de pontos pode ser cortada por meio das função clip_transect, e as coordenadas X e Z podem ser usadas para criar o gráfico.
las <- readLAS('data/MixedEucaNat.laz')
lasclass : LAS (v1.2 format 0)
memory : 27.5 Mb
extent : 203830, 203980, 7358900, 7359050 (xmin, xmax, ymin, ymax)
coord. ref. : SIRGAS 2000 / UTM zone 23S
area : 22500 m²
points : 553.5 thousand points
density : 24.6 points/m²
density : 17.93 pulses/m²p1 <- c(203830, 7358900)
p2 <- c(203980, 7359050)
las_tr <- clip_transect(las, p1, p2, width = 4, xz = TRUE)
las_tr@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 208.55 -1.11 771.82 337 2 2
2: 203.00 -0.40 771.96 36 4 4
3: 208.62 0.90 772.05 446 1 1
4: 206.50 -1.18 771.91 37 2 2
5: 206.26 0.08 771.78 229 2 2
---
21522: 152.16 -0.96 771.30 145 3 3
21523: 163.22 -0.22 772.56 59 2 2
21524: 165.74 0.22 772.58 469 1 1
21525: 171.58 -0.67 772.86 485 1 1
21526: 173.16 1.70 772.72 204 1 1
ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag
1: 0 0 2 FALSE
2: 0 0 2 FALSE
3: 0 0 2 FALSE
4: 0 0 2 FALSE
5: 0 0 2 FALSE
---
21522: 0 0 2 FALSE
21523: 0 0 2 FALSE
21524: 0 0 2 FALSE
21525: 0 0 2 FALSE
21526: 0 0 2 FALSE
Keypoint_flag Withheld_flag ScanAngleRank UserData PointSourceID
1: FALSE TRUE -21 0 14
2: FALSE TRUE -21 0 14
3: FALSE TRUE -21 0 14
4: FALSE TRUE -21 0 14
5: FALSE TRUE -21 0 14
---
21522: FALSE FALSE 19 0 16
21523: FALSE FALSE 1 0 15
21524: FALSE FALSE 1 0 15
21525: FALSE FALSE 1 0 15
21526: FALSE FALSE 0 0 15A renderização pode ser realizada com auxilio do pacote ggplot2.
library('ggplot2')
ggplot(las_tr@data, aes(X,Z, color = Z)) +
geom_point(size = 0.5) +
coord_equal() +
theme_minimal() +
scale_color_gradientn(colours = height.colors(50))
5 CLASSIFICAÇÃO DO SOLO
A classificação dos pontos terrestres é uma etapa importante no processamento de dados de nuvem de pontos. A distinção entre pontos terrestres e não terrestres permite a criação de um modelo contínuo de elevação do terreno. Muitos algoritmos foram relatados na literatura e o lidR atualmente fornece três deles, utilizáveis com a função classify_ground():
Filtro Morfológico Progressivo (PMF)
Função de Simulação de Tecido (CSF)
Classificação de Curvatura Multiescala (MCC)
5.1 Filtro Morfológico Progressivo
A implementação do algoritmo PMF lidRé baseada no método descrito em (Keqi Zhang et al. 2003), com algumas modificações técnicas.
O método original é baseado em raster, enquanto lidR realiza operações morfológicas baseadas em pontos, pois lidR é um software orientado a nuvem de pontos.
A etapa principal do método está resumida na figura abaixo:
A função pmf() requer a definição dos seguintes parâmetros de entrada:
ws(tamanho da janela ou sequência de tamanhos de janela)th(tamanho do limite ou sequência de alturas de limite).
Usuários mais experientes podem experimentar esses parâmetros para obter a melhor precisão de classificação; no entanto,lidR contém util_makeZhangParam() uma função que inclui os valores de parâmetro padrão descritos em (Keqi Zhang et al. 2003), no trabalho a seguir: https://ieeexplore.ieee.org/document/1202973
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn")
las <- classify_ground(las, algorithm = pmf(ws = 5, th = 3))
las@data X Y Z ReturnNumber NumberOfReturns Classification
1: 273357.1 5274360 806.5340 1 1 2
2: 273357.2 5274359 806.5635 2 2 2
3: 273357.2 5274358 806.0248 2 2 2
4: 273357.2 5274510 809.6303 2 2 2
5: 273357.2 5274509 809.3880 1 1 2
---
73399: 273642.8 5274579 806.2622 1 2 1
73400: 273642.8 5274578 806.9623 1 1 1
73401: 273642.8 5274577 807.2855 1 1 1
73402: 273642.8 5274577 807.1540 1 3 1
73403: 273642.8 5274576 806.0433 1 3 1Veja o resultado:
plot(las, color = "Classification", size = 3, bg = "white") Para ilustrar melhor os resultados da classificação podemos gerar e traçar uma secção transversal da nuvem de pontos.
p1 <- c(273420, 5274455)
p2 <- c(273570, 5274460)
las_tr <- clip_transect(las, p1, p2, width = 4, xz = TRUE)
las_tr@data X Y Z ReturnNumber NumberOfReturns Classification
1: 0.11950 3.58400 810.3363 1 1 2
2: 0.09175 2.01100 810.1248 1 1 2
3: 0.73725 2.69800 809.9230 1 1 2
4: 1.74300 3.33025 810.5065 1 1 2
5: 1.68300 0.19925 810.1843 1 1 2
---
734: 149.02375 0.11825 805.6260 1 1 2
735: 149.76875 2.99500 809.8782 1 1 1
736: 149.82350 2.23975 807.5643 1 2 1
737: 149.75125 0.66700 809.1952 1 2 1
738: 149.79275 0.70100 807.7798 2 2 2Construindo o gráfico.
library(ggplot2)
las_tr@data$Classification <- as.factor(las_tr@data$Classification)
ggplot(las_tr@data, aes(X,Z, color = Classification)) +
geom_point(size = 1) +
coord_equal() +
theme_minimal()
Como podemos ver, alguns pontos acima do terreno foram classificados como 2, o que de acordo com as especificações ASPRS indica solo.
Isso indica claramente que os parâmetros ws e th devem ser ajustados.
5.2 Funções de simulação de tecido
A filtragem de simulação de tecido (CSF) usa o algoritmo Zhang et al 2016 desenvolvido no estudo https://www.mdpi.com/2072-4292/8/6/501/htm e consiste em simular um pedaço de tecido enrolado sobre uma nuvem de pontos invertida. Neste método a nuvem de pontos é virada de cabeça para baixo e então um pano é jogado sobre a superfície invertida.
Os pontos de solo são determinados analisando as interações entre os nós do tecido e a superfície invertida. A simulação do tecido em si é baseada em uma grade que consiste em partículas com massa e interconexões que juntas determinam a posição tridimensional e a forma do tecido.
As csf()funções utilizam os valores padrão propostos por Zhang et al 2016 e podem ser utilizadas sem fornecer quaisquer argumentos.
library(RCSF)
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn")
mycsf <- csf(sloop_smooth = TRUE, class_threshold = 1, cloth_resolution = 1, time_step = 1)
las <- classify_ground(las, mycsf)
las <- classify_ground(las, algorithm = csf())
las@data X Y Z ReturnNumber NumberOfReturns Classification
1: 273357.1 5274360 806.5340 1 1 2
2: 273357.2 5274359 806.5635 2 2 2
3: 273357.2 5274358 806.0248 2 2 2
4: 273357.2 5274510 809.6303 2 2 2
5: 273357.2 5274509 809.3880 1 1 2
---
73399: 273642.8 5274579 806.2622 1 2 1
73400: 273642.8 5274578 806.9623 1 1 1
73401: 273642.8 5274577 807.2855 1 1 1
73402: 273642.8 5274577 807.1540 1 3 1
73403: 273642.8 5274576 806.0433 1 3 1Veja o resultado.
library(ggplot2)
las@data$Classification <- as.factor(las@data$Classification)
ggplot(las@data, aes(X,Z, color = Classification)) +
geom_point(size = 1) +
coord_equal() +
theme_minimal()
Podemos filtrar esses dados para exibir em 3D apenas o resultado do solo.
las_ground <- filter_ground(las)
las_ground@data X Y Z ReturnNumber NumberOfReturns Classification
1: 273357.1 5274360 806.5340 1 1 2
2: 273357.2 5274359 806.5635 2 2 2
3: 273357.2 5274358 806.0248 2 2 2
4: 273357.2 5274510 809.6303 2 2 2
5: 273357.2 5274509 809.3880 1 1 2
---
19364: 273642.7 5274625 790.5415 1 1 2
19365: 273642.7 5274623 790.7990 2 2 2
19366: 273642.8 5274614 791.9695 1 1 2
19367: 273642.8 5274609 793.1893 4 4 2
19368: 273642.8 5274600 795.5983 1 1 2plot(las_ground, bg = 'white')5.3 Classificação de Curvatura Multiescala (MCC)
A Classificação de Curvatura Multiescala (MCC) usa o algoritmo Evans e Hudak 2016 (Evans and Hudak 2007), no artigo entitulado Um algoritmo de curvatura multiescala para classificação LiDAR de retorno discreto em ambientes florestais; https://ieeexplore.ieee.org/document/4137852.
library(RMCC)
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn")
las <- classify_ground(las, mcc(1.5,0.3))Veja o resultado.
library(ggplot2)
las@data$Classification <- as.factor(las@data$Classification)
ggplot(las@data, aes(X,Z, color = Classification)) +
geom_point(size = 1) +
coord_equal() +
theme_minimal()
Resultado em 3D da filtragem do solo.
las_ground <- filter_ground(las)
plot(las_ground, bg = 'white')5.4 O problema da borda
Não importa qual algoritmo é usado em lidR ou outro software, a classificação do terreno sempre será mais fraca nas bordas das nuvens de pontos, pois o algoritmo deve analisar a vizinhança local.
6 MODELO DIGITAL DO TERRENO
A geração de um modelo digital de terreno (MDT) é geralmente a segunda etapa do processamento que segue a classificação dos pontos do terreno.
A construção de um DTM começa com pontos terrestres conhecidos ou amostrados e utiliza várias técnicas de interpolação espacial para inferir pontos terrestres em locais não amostrados.
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzc")
plot(las, size = 3, bg = "white")6.1 Rede irregular triangular
Este método é baseado em rede irregular triangular (TIN) de dados de pontos terrestres para derivar uma função bivariada para cada triângulo, que é então usada para estimar os valores em locais não amostrados (entre pontos terrestres conhecidos).
Para gerar um modelo DTM com o algoritmo TIN usamos a função rasterize_terrain() com o algorithm = tin().
dtm_tin <- rasterize_terrain(las, res = 1, algorithm = tin())
plot_dtm3d(dtm_tin, bg = "white")
#raster::writeRaster(dtm_tin, filename = 'E:/Geoprocessamento/LIDAR/Apresentação/data/dtm_tin.tif')6.2 Ponderação de distância invertida (IDW)
Baseia-se na suposição de que o valor em um ponto não amostrado pode ser aproximado como uma média ponderada de valores em pontos dentro de uma certa distância de corte d, ou a partir de um determinado número k de vizinhos mais próximos.
Para gerar um modelo DTM com o algoritmo IDW usamos rasterize_terrain() com o algorithm = knnidw()
dtm_idw <- rasterize_terrain(las, algorithm = knnidw(k = 10L, p = 2))
plot_dtm3d(dtm_idw, bg = "white") 6.3 Krigagem
A krigagem é a abordagem mais avançada e utiliza métodos avançados de interpolação geoestatística que levam em consideração as relações entre os retornos e suas respectivas distâncias entre si.
lidRusa o pacote gstatpara realizar a krigagem.
Para gerar um modelo DTM com o algoritmo de krigagem usamos rasterize_terrain()where algorithm = kriging().
dtm_kriging <- rasterize_terrain(las, algorithm = kriging(k = 40))
plot_dtm3d(dtm_kriging, bg = "white") 6.4 Executando o MDT e exibindo a plantação.
las <- readLAS(LASfile, select = "xyzc")
las@data X Y Z Classification
1: 273357.1 5274360 806.5340 1
2: 273357.2 5274359 806.5635 1
3: 273357.2 5274358 806.0248 2
4: 273357.2 5274510 809.6303 1
5: 273357.2 5274509 809.3880 2
---
73399: 273642.8 5274579 806.2622 1
73400: 273642.8 5274578 806.9623 1
73401: 273642.8 5274577 807.2855 1
73402: 273642.8 5274577 807.1540 1
73403: 273642.8 5274576 806.0433 1table(las@data$Classification)
1 2 9
61347 8159 3897 dtm_euc <- rasterize_terrain(las, res = 1, algorithm = tin())
las2 <- readLAS(LASfile, select = "xyzc", filter = "-keep_class 1")
table(las2@data$Classification)
1
61347 x <- plot(las2, bg = 'white', axis = T, legend = T)
add_dtm3d(x, dtm_euc)7 NORMALIZAÇÃO DA ALTURA
A normalização da nuvem de pontos remove a influência do terreno nas medições acima do solo. Isso torna possível a comparação das alturas da vegetação acima do solo e simplifica as análises nas áreas de aquisição.
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile)
plot(las, size = 3, bg = "white")Para ter uma ideia melhor da aparência do terreno, vamos remover todos os pontos que não são do solo
gnd <- filter_ground(las)
plot(gnd, size = 3, bg = "white", color = "Classification")7.1 Normalização do MDT
Para normalizar pontos usando um MDT, primeiro precisamos criar-lo. Para isso usamos a rasterize_terrain(). Para este exemplo vou utilizar uma resolução de grade de 1 m e utilizar o knnidw() algoritmo com parâmetros padrão.
dtm <- rasterize_terrain(las, 1, knnidw())
plot(dtm, col = gray(1:50/50))
Agora que temos a nossa superfície e estamos satisfeitos com ela, podemos utilizá-la para normalizar a nossa nuvem de pontos através da subtração.
nlas <- las - dtm
plot(nlas, size = 4, bg = "white")Avaliando como ficou a altura dos pontos após a normalização.
library(ggplot2)
dados_solo <- filter_ground(nlas)
dados_solo <- dados_solo@data
ggplot(data = dados_solo)+
geom_histogram(aes(Z), bins = 30, fill = 'tomato', color = 'white')+
labs(title = 'Distribuição das alturas dos pontos do solo', x = 'Elevação', y = 'Frequência')+
theme_bw()+
scale_x_continuous(n.breaks = 20)
Avaliando com um Boxplot.
ggplot(dados_solo)+
geom_boxplot(aes(x = Z))+
theme_bw()
7.2 Normalização de novem de pontos
Para calcular a normalização contínua, podemos alimentar normalize_height()com um algoritmo de interpolação espacial em vez de um raster.
nlas <- normalize_height(las, knnidw())Avaliando como ficou a altura dos pontos após a normalização.
library(ggplot2)
dados_solo <- filter_ground(nlas)
dados_solo <- dados_solo@data
table(dados_solo$Z)
0
8159 ggplot(dados_solo)+
geom_boxplot(aes(x = Z))+
theme_bw()
8 MODELO DIGITAL DE SUPERFÍCIE E MODELO DE ALTURA DO DOSSEL
Modelos Digitais de Superfície (DSM) e Modelos de Altura de Canopy (CHM) são camadas raster que representam - mais ou menos - a maior elevação dos retornos de ALS.
Se a nuvem de pontos for normalizada a superficie derivada representa a altura da copa e é chamada de CHM;
Quando a nuvem de pontos original é usada, a camada derivada representa a elevação do topo da copa acima do nível do mar e é chamada de DSM;
Nesta seção usaremos o MixedConifer.laz conjunto de dados normalizados, que é incluído internamente lidR para criar exemplos reproduzíveis.
LASfile <- system.file("extdata", "MixedConifer.laz", package ="lidR")
las <- readLAS(LASfile)
plot(las, size = 3, bg = "white")8.1 Ponto para raster
Consiste em estabelecer uma grade com uma resolução definida pelo usuário e atribuir a elevação do ponto mais alto a cada pixel.
chm <- rasterize_canopy(las, res = 1, algorithm = p2r())
col <- height.colors(25)
plot(chm, col = col)
Para eliminar pontos vazios pode-se utilizar métodos de interpolação como tin(), knnidw() e kriging().
chm <- rasterize_canopy(las, res = 0.5, p2r(0.2, na.fill = tin()))
plot(chm, col = col)
9 DETECÇÃO E SEGMENTAÇÃO DE ÁRVORES INDIVIDUAIS
A detecção individual de árvores (ITD) é o processo de localização espacial de árvores e extração de informações de altura.
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzr", filter = "-drop_z_below 0")
chm <- rasterize_canopy(las, 0.5, pitfree(subcircle = 0.2))
plot(las, bg = "white", size = 4)As copas das árvores podem ser detectadas aplicando um Filtro Máximo Local (LMF) no conjunto de dados carregado.
ttops <- locate_trees(las, lmf(ws = 5))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
Os resultados da detecção de árvores também podem ser visualizados em 3D!
x <- plot(las, bg = "white", size = 4)
add_treetops3d(x, ttops)Extraindo a localização das árvores.
LASfile <- system.file('extdata', 'MixedConifer.laz', package = 'lidR')
las <- readLAS(LASfile, select = "xyzc")
# Extrair a localização das árvores
ttops <- locate_trees(las, lmf(ws = 5))
# plotagem
offsets <- plot(las, bg = "white", size = 3)
add_treetops3d(offsets, ttops)
# Extração das coordenadas do topo das árvores
x <- sf::st_coordinates(ttops)[,1] - offsets[1]
y <- sf::st_coordinates(ttops)[,2] - offsets[2]
z <- ttops$Z
localizacao <- data.frame(x,y,z)References
Evans, Jeffrey S., and Andrew T. Hudak. 2007. “A Multiscale Curvature Algorithm for Classifying Discrete Return LiDAR in Forested Environments.” IEEE Transactions on Geoscience and Remote Sensing 45 (4): 1029–38. https://doi.org/10.1109/tgrs.2006.890412.
Keqi Zhang, Shu-Ching Chen, D. Whitman, Mei-Ling Shyu, Jianhua Yan, and Chengcui Zhang. 2003. “A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LIDAR Data.” IEEE Transactions on Geoscience and Remote Sensing 41 (4): 872–82. https://doi.org/10.1109/tgrs.2003.810682.