This book is in Open Review. I want your feedback to make the book better for you and other readers. To add your annotation, select some text and then click the on the pop-up menu. To see the annotations of others, click the in the upper right hand corner of the page

Cours 7 Work with LiDAR data

LiDAR : “light detection and ranging” ou “laser imaging detection and ranging”
Lumière d’un laser en général : visible, infrarouge ou ultraviolet.

Télémétrie : détermination de la distance d’un objet

Principe d’écholocalisation : La distance est donnée par la mesure du délai entre l’émission d’une impulsion et la détection d’une impulsion réfléchie, connaissant la vitesse de la lumière.

Mesure de la matière : isoler l’effet des différentes interactions entre la lumière et la matière le long du faisceau laser.

On appelle « équation lidar » le bilan de liaison du lidar entre son émission et sa réception, c’est-à-dire l’énergie lumineuse E (en J) de l’impulsion rétrodiffusée par une cible (supposée lambertienne, i.e. qui diffuse la lumière uniformément) située à une distance z et captée par le lidar.

Emet de plusieurs dizaines à plusieurs centaines de milliers d’impulsions chaque seconde.

Fichiers : - trajectoire : décrit les positions dans le temps - nuage de points (.las ou .laz en compréssé)

7.1 Tools

logiciels : cloudcompare, Fugroviewer

7.1.1 LAStools

Using LAStools, help here. In R :

# Define the path for lastools executables
LAStoolsDir =  "D:/Program/LAStools/bin/"

# Define the R function that calls lastools executables
LAStool <- function(tool, inputFile, ...){
  cmd = paste(paste(LAStoolsDir, tool, sep=''), '-i', inputFile , ...)
  cat(cmd)
  system(cmd)
  return(cmd)
}

# Define the directory for the ALS project
cloudPath <- "//amap-data.cirad.fr/work/users/VincyaneBadouard/Lidar/HovermapUAV2023/P16_C19C20/" 

# Define las/laz files to be processed
inFiles = paste(cloudPath, 'Output_laz1_4.laz', sep ='')
outFile = "FullDensity.laz"

# Define output directory (and creates it ... if doesn't exist)
outDir = paste(cloudPath, 'FullDens_split', sep ='')
dir.create(outDir) # , showWarnings = F (pas le droit de le faire dans le serveur safe)

# Call the lastools function to split in several files by gpstime
cores <- detectCores()

LASrun <- LAStool('lassplit64', inFiles,
                  '-cores', cores,
                  '-by_gps_time_interval 60',
                  '-odir', outDir,
                  '-odix _split',
                  '-o', outFile
)

# Merge files once corrected
outDir2 <- paste(cloudPath, 'FullDens_cor', sep='') # new folder: files rebinded and corrected
dir.create(outDir2, showWarnings = F)
inFiles <- paste(outDir, '/*.laz', sep='') # files corrected

LASrun <- LAStool('lasmerge64', inFiles,
                  '-odir', outDir2,
                  '-odix _cor',
                  '-o', outFile
)
# laz to dem
las2dem -i input/Paracou_284500_581500.laz -o output/dem.asc

7.2 Load laz

Le format binaire ASPRS LASer (.las) est un format ouvert permettant le stockage de données à 3 dimensions sous forme binaire. En plus des coordonnées x,y,z de chaque point, ce format permet de stocker des informations propres aux données LiDAR : - temps GPS - intensité - angle de scan - numéro du retour - classification - données utilisateurs - nombre de retours par impulsion - intensité du rouge - intensité du vert - intensité du bleu - indicateur de bord de ligne de vol - ID du point source - couleur RGB https://r-lidar.github.io/lidRbook/io.html

las1 <- lidR::readLAS(files = "D:/Mes Donnees/PhD/SIG_data/Lidar/P16_2022_ptf6.laz")
# ou
las2 <- rlas::read.las(files = "D:/Mes Donnees/PhD/SIG_data/Lidar/P16_2022_ptf6.laz")

7.2.1 Load several laz

# readLAScatalog() creates a las catalog from a folder to search in
ctg <- readLAScatalog(folder = "D:/Mes Donnees/PhD/Lidar", filter = "keep") # las catalog. filter only 1st returns
# readLAS(filter = "-help")
# "-first_only"
# , select = "xyz")  # load XYZ only
# filter = "-keep_first -drop_z_below 5 -drop_z_above 50")

# To apply a function to a LAScatalog
catalog_apply(ctg,fct)

7.2.2 Load/Create a subset of the laz to developpe the code with smaller data

range(las$gpstime) #  1697733332 1697734174
readLAS(filter = "-help")
las <- readLAS("Y:/users/VincyaneBadouard/Lidar/HovermapUAV2023/Translation/LAZ/UAV_C19C20_translated.las",
               filter = "-keep_gps_time 1697733332 1697733334")

sampletowork <- decimate_points(las, random(1)) # reduce point density (13 MB)
plot(sampletowork)

7.3 Write laz

writeLAS(las, "Z:/users/VincyaneBadouard/ALS2022/P16_2022_v2_4ha.laz")

7.4 Check if the las is complete

lidR::las_check(las1)

7.5 See their summaries

class(las1)
# [1] "LAS"
# attr(,"package")
# [1] "lidR"
class(las2) # "data.table" "data.frame"

print(las1) 
print(las2)
# area         : 294704 m²
# points       : 116.93 million points
# density      : 396.77 points/m²
# density      : 180.46 pulses/m²

7.6 Plot

lidR::plot(las1)
# persp(x = las2$X, y = las2$Y, z = las2$Z)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)

lidR::plot(las)

Ploter une raster layer : d’abord en faire un data.frame

TWI_df <- as.data.frame(TWI, xy = TRUE) %>%
  na.omit() %>%
  dplyr::rename('TWI' = 'P13_2022_TWI')

ggplot() + 
  geom_raster(data = TWI_df, aes(x = x, y = y, fill = TWI)) +
  scale_fill_viridis_c() +
  ggtitle("Paracou P13 - TWI") + 
  coord_sf()

7.7 Trajectory

Permet :
- de calculer les hauteurs de vol - d’estimer les incertitudes de positionnement des points (en fonction de la hauteur) - de normaliser les intensités des échos (si cela n’a pas été fait par le fournisseur)

Contient :
- The X axis is tangent to the geoid and is pointing to the north, the Z axis is pointing to the center of the earth, and the Y axis is tangent to the geoid and is pointing to the East. - temps GPS (seconde) - Latitude et longitude (radian) - Altitude (mètre) - Vitesse dans l’axe X (est-ouest), Y (nord-sud), Z (verticale) (m/s) - Quaternions q0, q1, q2, q3 - “r”, “g” ,“b” : red, green, blue - “nx” , “ny” , “nz” : coordonnées du normal vector - “roll” (roulis) : roll angle φ which is defined around the xb axis (radian) - “pitch” (tangage) : the pitch angle θ which is defined about the intermediate axis y1 (radian) - “yaw” (lacet) : the yaw angle ψ which is defined around the vertical axis zi (radian) - Angle du système d’axes par rapport au Nord (radian) - Accélération sur les axes X, Y, Z (m/s2) - Vitesse angulaire autour des axes X, Y, Z (rad/s)

traj <- fread("Z:/users/VincyaneBadouard/ZebHorizon/LiDAR_Manu/2022-09-07_15-46-34.gs-traj") # or .txt
traj 

plot3d(traj[,2:4], aspect=F) # plot the trajectory (x,y,z)

7.8 Clip the point cloud

  • selon une géométrie : clip_roi(las = ctg, geometry = ROI)
  • selon des coordonées : clip_rectangle()/clip_polygon()/clip_circle()
  • sur un LAScatalog (ça ne coupera pas dans une dalle mais ça réduit le nombre de dalle à la région d’intéret) : catalog_intersect(ctg, ROI)

7.9 Classification

7 = bruit 2 = sol 9 = eau vegetation = 0,3,4,5

table(ST@data$Classification) # 7 = bruit, 2 = sol, veg = 3,4,5
ST@data <- ST@data[Classification != 7,] # enlever le bruit

Classifier le sol

# Modèle numérique de terrain
mnt <- rast("Z:/users/VincyaneBadouard/ALS2022/mnt_roi36ha_1m.asc")
plot(mnt)

# Normaliser la hauteur du sol
ST_norm <- lidR::normalize_height(ST, mnt) # applati le relief
plot(ST_norm)# sol plat

# Classifier les points
ST_norm@data[(Z>(-0.2) & Z< 0.5), Classification := 2] 

table(ST_norm@data$Classification) # 7 = bruit, 2 = sol, vegetation = 3,4,5
# plot(ST_norm, color = "Classification") # trop long

# Dénormaliser la hauteur du sol
ST <- lidR::unnormalize_height(ST_norm)

7.10 Ajouter des points à un las, combiner plusieurs las

# transform in a LAS object
las2 <- LAS(vector_coord[,.(XB,YB,ZB)], las1@header)

# This triggers some warning because we input some incorrectly quantized coordinates. Check the output of las_check()
las_check(las2)
# We can fix that
las2 <- las_quantize(las2, TRUE)
las2 <- las_update(las2)

las_check(las2)

# To combine the las1 and las2 they must have the same columns.
# We need to make them manually.
# Here if you load only XYZ in las the job is easier.

las2@data$gpstime = 0
las2@data$Intensity = 0L
las2@data$ReturnNumber = 1L
las2@data$NumberOfReturns = 1L
las2@data$ScanDirectionFlag = 0L
las2@data$EdgeOfFlightline = 0L
las2@data$Classification = 0L
las2@data$Synthetic_flag = FALSE
las2@data$Keypoint_flag = FALSE
las2@data$Withheld_flag = FALSE
las2@data$ScanAngleRank = 0L
las2@data$UserData = 0L
las2@data$PointSourceID = 0L

# Bind
las3 <- lidR::rbind(las1, las2)

7.11 Caluler des métriques

7.11.1 Generate a CHM (Canopy Height Model) and DSM (Digital Surface Model)

In the case of a normalized point cloud, the derived surface represents the canopy height (for vegetated areas) and is referred to as CHM. When the original (non-normalized) point cloud with absolute elevations is used, the derived layer represents the elevation of the top of the canopy above sea level, and is referred to as DSM

CHM <- rasterize_canopy(normalize_height(ST, knnidw()), res = 1, algorithm = p2r()) 
# p2r() : for each pixel of the output raster the function attributes the height of the highest point found
# tin() : triangulation algorithm
plot(CHM, col = height.colors(25))

DSM <- rasterize_canopy(ST, res = 0.5, algorithm = dsmtin())
# dsmtin() : algorithm for digital surface model computation using a Delaunay triangulation of first returns with a linear interpolation within each triangle.

ggplot() + 
  geom_spatraster(data = CHM, aes(fill = Z)) +
  scale_fill_gradientn(name = "Canopy height (m)",
                       colors = height.colors(25), na.value="white") + 
  geom_sf(data = sf::st_cast(ROI, "LINESTRING")) +
  ggtitle("Paracou P16 4ha - CHM ALS 2023 - 1m res") + 
  theme_classic() +
  coord_sf()

7.11.2 Generate a DTM (MNT)

https://r-lidar.github.io/lidRbook/dtm.html l’élévation est dans la variable Z (dtm$Z)

# Triangulation (généralement utilisée)
dtm_tin <- rasterize_terrain(ST, res = 1, algorithm = tin())
plot_dtm3d(dtm_tin, bg = "black") 

# Invert distance weighting
dtm_idw <- rasterize_terrain(ST, algorithm = knnidw(k = 10L, p = 2))
plot_dtm3d(dtm_idw, bg = "black") 

# Kriging
dtm_kriging <- rasterize_terrain(ST, algorithm = kriging(k = 40))
plot_dtm3d(dtm_kriging, bg = "black") 

mnt <- dtm_tin

ggplot() + 
  geom_spatraster(data = mnt, aes(fill = Z)) +
  scale_fill_gradientn(name = "Elevation (m)",
                       colors = terrain.colors(20), na.value="white") +
  theme_classic()

7.11.3 La réflectance (albedo)

Calculer la réflectance apparente en ratio (albedo)

Réflectance (ou Albedo) = énergie lumineuse réfléchie/énergie lumineuse incidente

VP$initial_intensity <- VP$Intensity 

# Intensité = réflectance apparente en ratio (albedo)
VP$Intensity = 10^(VP$Reflectance/10) # Réflectance initialment en decibel (et selon une référence connue)