print(las)
#> class : LAS (v1.2 format 1)
#> memory : 4.4 Mb
#> extent : 684766.4, 684993.3, 5017773, 5018007 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83 / UTM zone 17N
#> area : 51572 m²
#> points : 81.6 thousand points
#> density : 1.58 points/m²
#> density : 1.08 pulses/m²2 Reading, Plotting, Querying & Validating
2.1 Reading LiDAR data using readLAS
Discrete return ALS sensors record various types of data. Primarily, they capture positional data in three dimensions (X, Y, Z), followed by additional information like the intensity for each point, the position of each point in the return sequence, and the beam incidence angle of each point. Reading, writing, and efficient storage of these ALS data are critical steps prior to any subsequent analysis.
ALS data are most commonly distributed in LAS format, which is specifically designed to store ALS data in a standardized way. These data are officially documented and maintained by the ASPRS. However, LAS files require a large amount of memory because they are not compressed. The LAZ format has become the standard compression scheme because it is free and open-source.
The widespread use, standardization, and open-source nature of the LAS and LAZ formats promoted the development of the lidR package. This package is designed to process LAS and LAZ files both as input and output, leveraging the LASlib and LASzip C++ libraries via the rlas package.
The function readLAS() reads a LAS or LAZ file and returns an object of class LAS. The LAS formal class is documented in detail in a dedicated vignette. To briefly summarize, a LAS file consists of two parts:
- The header, which stores summary information about its content, including the bounding box of the file, coordinate reference system, and point format.
- The payload, i.e., the point cloud itself.
The function readLAS() reads and creates an object that contains both the header and the payload.
las <- readLAS("files.las")When printed it displays a summary of its content.
For a more in-depth print out of the data use the function summary() instead of print().
Parameter select
A LAS file stores the X Y Z coordinates of each point as well as many other data such as intensity, incidence angle, and return sequence position. These data are called attributes. In practice, many attributes are not actually useful but are loaded by default. This can consume a lot of processing memory because R does not allow for choosing data storage modes (see this vignette for more details).
To save memory, readLAS() can take an optional parameter select, which enables the user to selectively load the attributes of interest. For example, one can choose to load only the X Y Z attributes.
las <- readLAS("file.las", select = "xyz") # load XYZ only
las <- readLAS("file.las", select = "xyzi") # load XYZ and intensity onlyExamples of other attribute abbreviations are: t - gpstime, a - scan angle, n - number of returns, r - return number, c - classification, s - synthetic flag, k - keypoint flag, w - withheld flag, o - overlap flag (format 6+), u - user data, p - point source ID, e - edge of flight line flag, d - direction of scan flag
Parameter filter
While select enables the user to choose “columns” (or attributes) while reading files, filter allows selection of “rows” (or points) during the reading process. Removing superfluous data at read time saves memory and increases computation speed. For example, it’s common practice in forestry to process only the first returns.
las <- readLAS("file.las", filter = "-keep_first") # Read only first returnsIt is important to understand that the filter option in readLAS() keeps or discards points at read time, i.e., while reading at the C++ level, without involving any R code. For example, the R function filter_poi() may return the same output as the filter option in readLAS():
las1 <- readLAS("file.las", filter = "-keep_first")
las2 <- readLAS("file.las")
las2 <- filter_poi(las2, ReturnNumber == 1L)In the example above, we are (1) reading only the first returns or (2) reading all the points and then filtering the first returns in R. Both outputs are strictly identical, but the first method is faster and more memory-efficient because it doesn’t load the entire file into R and avoids using extra processing memory. It should always be preferred when possible. Multiple filter commands can be used simultaneously to, for example, read only the first returns between 5 and 50 meters.
las <- readLAS("file.las", filter = "-keep_first -drop_z_below 5 -drop_z_above 50")The full list of available commands can be obtained by using readLAS(filter = "-help"). Users of LAStools may recognize these commands, as both LAStools and lidR use the same libraries (LASlib and LASzip) to read and write LAS and LAZ files.
2.2 Validating LiDAR Data
An important first step in ALS data processing is ensuring that your data is complete and valid according to the ASPRS LAS specifications. Users commonly report bugs arising from invalid data. This is why we introduced the las_check() function to perform a thorough inspection of LAS objects. This function checks whether a LAS object meets the ASPRS LAS specifications and whether it is valid for processing, providing warnings if it does not.
A common issue is that a LAS file contains duplicate points. This can lead to problems such as trees being detected twice, invalid metrics, or errors in DTM generation. We may also encounter invalid return numbers, incoherent return numbers and number of returns attributes, and invalid coordinate reference systems, among other issues. Always make sure to run the las_check() function before delving deeply into your data.
las_check(las)
#> Checking the data
#> - Checking coordinates... ✓
#> - Checking coordinates type... ✓
#> - Checking coordinates range... ✓
#> - Checking coordinates quantization... ✓
#> - Checking attributes type... ✓
#> - Checking ReturnNumber validity...
#> ⚠ Invalid data: 1 points with a return number equal to 0 found.
#> [...]A check is performed at read time regardless, but the read time check is not as thorough as las_check() for computation time reasons. For example duplicated points are not checked at read time.
las <- readLAS("data/chap1/corrupted.laz")
#> Warning: Invalid data: 174638 points with a 'return number' greater than the
#> 'number of returns'.2.3 Plotting
The lidR package takes advantage of the rgl package to provide a versatile and interactive 3D viewer with points colored by Z coordinates on a black background as default.
Basic 3D rendering
The very basic way to render a point cloud is the function plot().
plot(las)
Users can change the attributes used for coloring by providing the name of the attribute used to colorize the points. The background color of the viewer can also be changed by assigning a color using the bg argument. Axes can also be added and point sizes can be changed.
# Plot las object by scan angle,
# make the background white,
# display XYZ axis and scale colors
plot(las, color = "ScanAngleRank", bg = "white", axis = TRUE, legend = TRUE)
Note that if your file contains RGB data the string "RGB" is supported:
plot(las, color = "RGB")The argument breaks enables to defined more adequate breaks in the color palette for example when intensity contains large outliers. Otherwise the palette range would be too large and most of the values would be considered as “very low”, so everything would appear in the same color.
plot(las, color = "Intensity", breaks = "quantile", bg = "white")
Overlays
The package also provides some easy to use functions for common overlay. For example add_dtm3d() to add a digital terrain model (section Chapter 4)) and add_treetops3d() to visualize the output of an individual tree detection (section Section 7.1))
x <- plot(las, bg = "white", size = 3)
add_dtm3d(x, dtm)
x <- plot(las, bg = "white", size = 3)
add_treetops3d(x, ttops)
It is also possible to combine two point clouds with different color palettes. In the following example we are using a previously classified point cloud. We first separate the vegetation and non vegetation points using filter_poi() and then plot both on top of each other with different colour schemes using add options in plot()
nonveg <- filter_poi(las, Classification != LASHIGHVEGETATION)
veg <- filter_poi(las, Classification == LASHIGHVEGETATION)
x <- plot(nonveg, color = "Classification", bg = "white", size = 3)
plot(veg, add = x)
Advanced 3D Rendering
Since lidR is based on rgl, it is easy to add objects to the main rendering using rgl functions such as rgl::point3d(), rgl::text(), rgl::surface3d(), and so on to produce publication-ready renderings. However, lidR introduces an additional challenge: it does not display the points with their actual coordinates. Instead, the points are shifted to be rendered close to (0, 0) due to accuracy issues, as rgl uses float (32-bit decimal numbers) rather than double (64-bit decimal numbers). When plot() is used, it invisibly returns the shift values, which can later be used to realign other objects.
offsets <- plot(las)
print(offsets)
#> [1] 391867.8 3901019.3The coordinates of the objects must be corrected to align with the point cloud. In the following we will add lines to render the trunks. We read a file, we locate the trees (see Section 7.1)), we extract the coordinates and sizes of the trees and plot lines with rgl::segment3d().
Show the code
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzc")
# get the location of the trees
ttops <- locate_trees(las, lmf(ws = 5))
# plot the point cloud
offsets <- plot(las, bg = "white", size = 3)
add_treetops3d(offsets, ttops)
# extract the coordinates of the trees and
# apply the shift to display the lines
# in the rendering coordinate system
x <- sf::st_coordinates(ttops)[,1] - offsets[1]
y <- sf::st_coordinates(ttops)[,2] - offsets[2]
z <- ttops$Z
# Build a GL_LINES matrix for fast rendering
x <- rep(x, each = 2)
y <- rep(y, each = 2)
tmp <- numeric(2*length(z))
tmp[2*1:length(z)] <- z
z <- tmp
M <- cbind(x,y,z)
# Display lines
rgl::segments3d(M, col = "black", lwd = 2)
../../../../../../tmp/RtmpNMWMOm/file1737024096ea6.png
Voxel rendering
It is possible to render voxels. This is useful to render the output of the function voxelise_points() or voxel_metrics() for examples.
vox <- voxelize_points(las, 6)
plot(vox, voxel = TRUE, bg = "white")
Cross Sections 2D Rendering
To better visualize the vertical structure of a point cloud, investigate classification results, or compare the results of different interpolation routines, a cross section can be plotted. To do this, we first need to decide where the cross section will be located (i.e., define the beginning and end) and specify its width. The point cloud can then be clipped, and the X and Z coordinates used to create the plot.
For example, to create a 200 m long cross section, we might define the beginning and end, and then use the clip_transect() function to subset the point cloud.
p1 <- c(273457, 5274357)
p2 <- c(273542, 5274542)
las_tr <- clip_transect(las, p1, p2, width = 5, xz = TRUE)Rendering can be achieved with base plot or ggplot2. Notice the use of payload() to extract the data.frame from the LAS object.
library(ggplot2)
ggplot(payload(las_tr), aes(X,Z, color = Z)) +
geom_point(size = 0.5) +
coord_equal() +
theme_minimal() +
scale_color_gradientn(colours = height.colors(50))