library(lidR)
library(ggplot2)
library(dplyr)Análise dos dados - SABRINA
1 Pacotes carregados
2 Leitura dos dados
Primeira leitura dos dados completos.
dados <- readLAS('E:/Geoprocessamento/LIDAR/sabrina/lidars/terra_las/cloud281cf13d0ecc133c.las')2.1 Análise preliminar dos dados.
dadosclass : LAS (v1.2 format 3)
memory : 1.9 Gb
extent : 527354.8, 527778, 9695128, 9695529 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 22S
area : 126816 m²
points : 34.14 million points
density : 269.21 points/m²
density : 239.42 pulses/m²Assim podemos observar que esse documento lido é da classe .LAS, tem 1,9 Gb, a varredura foi feita de uma área de 126.816 m², gerando uma nuvem com 34,14 milhões de pontos, cuja densidade ficou em torno de 269,21 pontos/m² com uma densidade de varedura laser emitindo 239,42 pulsos/m².
Com a leitura do cabeçalho dos dados, podemos ver mais informações relevantes.
dados@headerFile signature: LASF
File source ID: 0
Global encoding:
- GPS Time Type: Standard GPS 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: REpJLUwxLUFDQy4=
Generating software: DJI TERRA 3.11.13.0
File creation d/y: 109/2023
header size: 227
Offset to point data: 1089
Num. var. length record: 3
Point data format: 3
Point data record length: 34
Num. of point records: 34140305
Num. of points by return: 30361913 3623422 154970 0 0
Scale factor X Y Z: 0.0001 0.0001 0.0001
Offset X Y Z: 527566.4 9695328 74.48956
min X Y Z: 527354.8 9695128 7.76886
max X Y Z: 527778 9695529 141.2103
Variable Length Records (VLR):
Variable Length Record 1 of 3
Description: GeoTIFF GeoKeyDirectoryTag
Tags:
Key 1024 value 1
Key 1025 value 1
Key 1026 value 0
Key 2049 value 22
Key 2054 value 9102
Key 3072 value 32722
Key 3076 value 9001
Variable Length Record 2 of 3
Description: GeoTIFF GeoAsciiParamsTag
data: WGS 84 / UTM zone 22S|WGS 84|
Variable Length Record 3 of 3
Description: OGR variant of OpenGIS WKT SRS
WKT OGC COORDINATE SYSTEM: [...] (truncated)
Extended Variable Length Records (EVLR): voidFazendo a leitura da carga útil em sí, temos a noção do formato da tabela de atributos dos dados.
dados@data X Y Z gpstime Intensity ReturnNumber
1: 527759.2 9695238 51.56856 366049825 17 1
2: 527757.7 9695238 52.45726 366049825 17 1
3: 527757.2 9695238 52.33056 366049825 22 1
4: 527758.2 9695238 50.07036 366049825 35 1
5: 527756.6 9695238 51.51956 366049825 22 1
---
34140301: 527459.9 9695410 49.74826 366050036 0 3
34140302: 527446.4 9695410 51.81346 366050036 0 3
34140303: 527446.8 9695410 51.50966 366050036 0 3
34140304: 527428.3 9695411 44.19926 366050036 0 3
34140305: 527428.9 9695411 55.33526 366050036 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
---
34140301: 3 0 0 0
34140302: 3 0 0 0
34140303: 3 0 0 0
34140304: 3 0 0 0
34140305: 3 0 0 0
Synthetic_flag Keypoint_flag Withheld_flag ScanAngleRank UserData
1: FALSE FALSE FALSE 35 0
2: FALSE FALSE FALSE 34 0
3: FALSE FALSE FALSE 34 0
4: FALSE FALSE FALSE 34 0
5: FALSE FALSE FALSE 34 0
---
34140301: FALSE FALSE FALSE 3 0
34140302: FALSE FALSE FALSE -4 0
34140303: FALSE FALSE FALSE -4 0
34140304: FALSE FALSE FALSE -13 0
34140305: FALSE FALSE FALSE -14 0
PointSourceID R G B
1: 0 19200 28672 14080
2: 0 25856 36608 17664
3: 0 21248 26880 17408
4: 0 18432 24320 15104
5: 0 14848 19200 11008
---
34140301: 0 34048 37120 29952
34140302: 0 2304 4096 2304
34140303: 0 10496 9728 7168
34140304: 0 18176 23296 11008
34140305: 0 31744 31488 21760Podemos observar que a tabela de atributos pussui 34.140.305 linhas e 19 colunas, onde cada linha se refere a um ponto da nuvem e cada coluna especifica um tipo de informação para o ponto.
plot(dados, color = 'Z', bg = 'white', axis = T, legend = T)
Podemos observar que alguns pontos estão muito alto, o que pode ser um erro ou um objeto estranho como um pássaro por exemplo.
3 Validação dos dados
Essa etapa é importate para garantir que seus dados estejam completos e válidos de acordo com as especificações.
las_check(dados)
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...[0;37m skipped[0m
- Checking attribute population...
[0;32m 🛈 'PointSourceID' attribute is not populated[0m
[0;32m 🛈 'ScanDirectionFlag' attribute is not populated[0m
[0;32m 🛈 'EdgeOfFlightline' attribute is not populated[0m
- Checking gpstime incoherances[0;32m ✓[0m
- Checking flag attributes...[0;32m ✓[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;31m no[0m
- Checking normalization...[0;31m no[0m
- Checking negative outliers...[0;32m ✓[0m
- Checking flightline classification...[0;31m no[0m
Checking compression
- Checking attribute compression...[0;31m no[0m4 Filtragem da base e seleção das colunas de dados
Com esse código eu realizei a leitura dos dados e LAS e armazenei em uma variável chamada dados2, selecionando com o parâmetro select apenas as colunas das coordenadas dos pontos (X,Y,Z), a coluna de intensidade, número retornado e número de retorno. E filtrei as linhas com o parâmetro filter = ‘-drop_z_above 100’ para filtrar apenas os dados que tem coordenadas Z abaixo de 100, dessa forma eliminando os dados outliers.
dados2 <- readLAS('E:/Geoprocessamento/LIDAR/sabrina/lidars/terra_las/cloud281cf13d0ecc133c.las',
select = 'xyzirn',
filter = '-drop_z_above 100')Ao plotar esses dados, vemos que os eixos estão mais certos com os dados da realidade e que o processo é executado mais rapidamente, pois a base está mais leve depois da seleção dos dados.
4.1 Avaliando os dados filtrados
dados2class : LAS (v1.2 format 3)
memory : 1.1 Gb
extent : 527354.8, 527778, 9695128, 9695529 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 22S
area : 126816 m²
points : 34.14 million points
density : 269.21 points/m²
density : 239.42 pulses/m²Pode fazer a mesma interpretação feita durante o capítulo 1.
dados2@data X Y Z Intensity ReturnNumber NumberOfReturns
1: 527759.2 9695238 51.56856 17 1 1
2: 527757.7 9695238 52.45726 17 1 1
3: 527757.2 9695238 52.33056 22 1 1
4: 527758.2 9695238 50.07036 35 1 1
5: 527756.6 9695238 51.51956 22 1 1
---
34140286: 527459.9 9695410 49.74826 0 3 3
34140287: 527446.4 9695410 51.81346 0 3 3
34140288: 527446.8 9695410 51.50966 0 3 3
34140289: 527428.3 9695411 44.19926 0 3 3
34140290: 527428.9 9695411 55.33526 0 3 3Para avaliar a distribuição das intensidades dos retornos dos pontos, pode aplicar esse código.
ggplot(dados2@data)+
geom_histogram(aes(Intensity), fill = 'tomato', color = 'white', binwidth = 10)+
theme_bw()+
scale_x_continuous(name = 'Intensity', n.breaks = 20)+
labs(title = 'Distribuição da intensidade dos pontos',y = 'Contagem')
Para avaliar a distribuição da altura dos pontos, pode aplicar esse código.
ggplot(dados2@data)+
geom_histogram(aes(Z), fill = 'tomato', color = 'white', binwidth = 5)+
theme_bw()+
scale_x_continuous(name = 'Altura', n.breaks = 20)+
labs(title = 'Distribuição da altura dos pontos',y = 'Contagem')
Pode avaliar, também, como estão se comportando a coluna dos números de retorno.
table(dados2@data$NumberOfReturns)
1 2 3
26708355 6963470 468465 Vemos que 30361898 pontos foram de um retorno e assim por diante.
Agora pode plotar o dados2.
plot(dados2, color = 'Z', bg = 'white', axis = T, legend = T)