# LAScatalog formal class

knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(lidR)  data = data.table::data.table( Max.X = c(885228.88, 886993.96, 885260.93, 887025.96, 885292.94, 887056.88, 892199.94, 893265.54, 892229.99, 893295.15, 888759.96, 890524.95, 892259.98, 894025.98, 892289.96, 894055.93, 888790.91, 890554.98, 888820.95, 890585.99, 892319.96, 894084.97, 892349.89, 894114.29, 895250.23, 895094.78, 895044.96, 895053.55, 885323.96, 887087.95, 885355.95, 887119.96, 883657.85, 885387.95, 887150.97, 885419.98, 887182.95, 883688.44, 885442.91, 887193.9, 888851.96, 890615.97, 888882.94, 890646.97, 892379.94, 894127.84, 892409.97, 892676.58, 888913.92, 890676.93, 888944.86, 890707.98, 892439.95, 894124.59, 892469.86, 894232.94, 894786.68, 888958.83, 890713.51, 892476.43, 894239.97, 894786.07), Min.X = c(885022.37, 885204.73, 885027.52, 885229.03, 885040.86, 885261.03, 891503.09, 892198.69, 891501.42, 892200.07, 886970.07, 888735.55, 891499.96, 892230.05, 890521.99, 892260.01, 886994.05, 888760.09, 887026.07, 888791.01, 890525.05, 892290.04, 890555.01, 892320.12, 894002.98, 894026.02, 894056.02, 894085.03, 885051.45, 885293.03, 885063.29, 885324.03, 883166.09, 885072.16, 885356.09, 883642.36, 885388.15, 883180.23, 883658.11, 885420.02, 887057.07, 888821.02, 887088.11, 888852.03, 890586.03, 892350.02, 890616.07, 892380.01, 887120.07, 888883.03, 887151.11, 888914.02, 890647.06, 892410.06, 890677.07, 892440.07, 894209.19, 887183.07, 888945.12, 890708.03, 892470.16, 894233.07), Max.Y = c(630219.48, 630214.96, 631609.95, 631604.97, 633001.65, 632995.99, 625898.35, 625882.94, 627289.82, 627273.89, 630174.88, 630134.94, 628681.66, 628664.99, 630094.95, 630057.95, 631564.98, 631524.94, 632955.82, 632915.99, 631486.9, 631447.96, 632876.93, 632838.96, 628627.89, 630019.93, 631410.97, 631740.88, 634393.05, 634386.96, 635786.24, 635779.96, 638613.36, 637176.84, 637169.92, 638601.99, 638560.96, 639938.36, 639926.95, 639558.31, 634346.93, 634307.92, 635739.92, 635699.92, 634268.97, 634229.95, 635659.89, 635622.88, 637129.84, 637089.93, 638520.94, 638481.91, 637051.99, 637012.92, 638442.98, 638403.94, 638366.87, 639177.04, 639133.74, 638702.56, 638702.56, 638702.56), Min.Y = c(629157.18, 629099.31, 630215.04, 630175.05, 631605.02, 631565.05, 625816.52, 625793.6, 625883.01, 625860.81, 629036.82, 629017.72, 627274.01, 627251.36, 628665.04, 628628.01, 630135.08, 630095.02, 631525.01, 631487.19, 630058.02, 630020.05, 631448.08, 631411.03, 627506.32, 628612.41, 629999.84, 631390.38, 632996.06, 632956.04, 634387.01, 634347.01, 637939.24, 635780.07, 635740.05, 637170.11, 637130.14, 638602.13, 638561.04, 638521.07, 632916.05, 632877.04, 634308.06, 634269.04, 632839.06, 632801.04, 634230.04, 634223.9, 635700.07, 635660.11, 637090.03, 637052.15, 635623.06, 635619.13, 637013.1, 636979.71, 637259.33, 638482.01, 638443.02, 638404.08, 638367.11, 638355.37), Min.Z = c(325.12, 251.48, 244.68, 286.7, 338.86, 320.68, 118.08, 60.69, -5.01, -7.58, 225.29, 252, 87.3, 41.7, 115.01, 28.77, 205.11, 200.85, 200.54, 169.5, 90.64, 19.72, 126.04, 28.6, 41.98, 43.15, 7.74, 6.7, 199.26, 190.02, 284.92, 216.16, 218.14, 318.93, 220.21, 218.04, 137.31, 218.13, 217.34, 190.42, 207.62, 113.84, 118.18, 141.52, 92, 52.5, 91.2, 77.92, 156.57, 89.53, 108.83, 93.98, 23.45, 1.64, 33.22, 3.29, 0.61, 108.03, 208.38, 121.18, 58.83, 0.95), Max.Z = c(418.46, 990.54, 409.06, 1021.87, 996.42, 1005.02, 173.77, 393.97, 836.52, 820.98, 414.2, 936.47, 792.95, 822.51, 777.31, 837.87, 419.15, 741.84, 907.2, 872.27, 898.53, 822.53, 846.77, 740.65, 826.61, 890.21, 828.86, 680.32, 390.31, 997.2, 965.55, 969.24, 249.34, 849.5, 950.2, 848.64, 904.1, 880, 827.92, 888.34, 462.88, 906.61, 440.83, 887.34, 860.37, 747.1, 808.75, 194.76, 734.21, 838.34, 834.76, 758.91, 771.76, 670.1, 810.94, 761.53, 109.26, 303.94, 349.94, 799.8, 737.01, 593.91), filename = paste0("path/to/las/files/file", 1:62, ".las") ) pgeom <- lapply(1:nrow(data), function(i) { mtx <- matrix(c(data$Min.X[i], data$Max.X[i], data$Min.Y[i], data$Max.Y[i])[c(1, 1, 2, 2, 1, 3, 4, 4, 3, 3)], ncol = 2) sp::Polygons(list(sp::Polygon(mtx)),as.character(i)) }) Sr = sp::SpatialPolygons(pgeom, proj4string = sp::CRS("+init=epsg:3005")) ctg <- new("LAScatalog") ctg@bbox <- Sr@bbox ctg@proj4string <- Sr@proj4string ctg@plotOrder <- Sr@plotOrder ctg@data <- data ctg@polygons <- Sr@polygons A LAScatalog class is a representation in R of a file or a set of las files not loaded in memory. Indeed, a regular computer cannot load the entire point cloud in R if it covers a wide area. For very high density datasets it can even fail loading a single file (see also the "LAS formal class"" vignette). In lidR, we use a LAScatalog to process datasets that cannot fit in memory. ## Build a LAScatalog object reading a folder of las files ctg <- catalog("path/to/las/files/") # or ctg <- catalog("path/to/las/files/big_file.las") print(ctg) ## Basic structure of a LAScatalog object A LAScatalog inherits a SpatialPolygonsDataFrame. Thus, it has the structure of a SpatialPolygonsDataFrame plus some extra slots that store information relative to how the LAScatalog will be processed. The slot data of a LAScatalog object contains a data.frame with the most important information read from the header of .las or .laz files. Reading only the header of the file provides an overview of the content of the files very quickly without actually loading the point cloud. The columns of the table are named after the LAS specification version 1.4 The other slots are well documented in the documentation ?LAScatalog-class so they are not described in this vignette. ## Allowed and non-allowed manipulation of a LAScatalog object A LAScatalog is formally a SpatialPolygonsDataFrame although its purpose is not to manipulate spatial data in R. The purpose of a LAScatalog is to represent a set of existing las/laz files. Thus a LAScatalog cannot be modified because it must be related to the actual content of the files. The following throws an error: ctg$Min.Z <- 0

Obviously it is always possible to modify an R object by bypassing such simple restrictions. In this case the user will break something internally and a correct output is not guaranteed.

## Validation of LAScatalog object

Users commonly report bugs arising from the fact that the point cloud is invalid. This is why we introduced the function lascheck to perform an inspection of the LAScatalog objects. This function checks if a LAScatalog object is consistent (files are all of the same type for example). For example, it may happen that a catalog mixes files of type 1 with files of type 3 or files with different scale factors.

lascheck(ctg)

The function lascheck when applied to a catalog does not perform a deep inspection of the point cloud unlike when applied to a LAS object. Indeed the point cloud is not actually read.

## Display a LAScatalog object

lidR provides a simple plot function to plot a LAScatalog object:

plot(ctg)

The option mapview = TRUE displays the catalog on an interactive map with pan and zoom and allows the addition of a satellite map in the background. It uses the package mapview internally. It is often useful to check if the CRS of the file are properly registered. The epsg code of the las file appears to be often incorrect, according to our own experience.

plot(ctg, mapview = TRUE, map.type = "Esri.WorldImagery")

.

Being a SpatialPolygonsDataFrame, the function spplot from sp can be used on a LAScatalog to check the content of the files. In the following we can immediately see that the catalog is not normalized and is likely to contain outliers:

spplot(ctg, "Min.Z")

## Apply lidR functions on a LAScatalog

Most of lidR functions are compatible with a LAScatalog and work almost like with a single point cloud loaded in memory. In the following example we use the function grid_metrics to compute the mean elevation of the points. The output is a continuous wall-to-wall RasterLayer. It works exactly as if the input was a LAS object.

hmean <- grid_metrics(ctg, mean(Z), 20)

However, processing a LAScatalog usually requires some tuning of the processing options to get better control of the computation. Indeed, if the catalog is huge the output is likely to be huge as well, and maybe the output cannot fit in the R memory. For example, lasnormalize throws an error if used 'as is' without tuning the processing options. Using lasnormalize like in the following example the expected output would be a huge point cloud loaded in memory. The lidR package forbids such a call:

output <- lasnormalize(ctg, tin())

Instead, one can use the processing option opt_output_files. Processing options drive how the big files are split in small chunks and how the outputs are either returned into R or written on disk into files.

opt_output_file(ctg) <- "folder/where/to/store/outputs/{ORIGINALFILENAME}_normalized" output <- lasnormalize(ctg, tin())

Here the output is not a point cloud but a LAScatalog pointing to the newly created files. The user can check how the catalog will be processed by calling summary

summary(ctg)

Also the plot function can displays the chunks pattern i.e. how the dataset is split into small chunks that will be sequentially processed

opt_chunk_size(ctg) <- 0 plot(ctg, chunk = TRUE) opt_chunk_size(ctg) <- 900 plot(ctg, chunk = TRUE)

## Some practical examples

#### Example 1 - Raster

Load a catalog. Process each file sequentially. Returns a RasterLayer into R.

ctg <- catalog("path/to/las/files/") hmean <- grid_metrics(ctg, mean(Z), 20)

#### Example 2 - Raster

Load a catalog. Process each file sequentially. For each file write a raster on disk named after the name of the processed files. Returns a lightweight virtual raster mosaic.

ctg <- catalog("path/to/las/files/") opt_output_files(ctg) <- "folder/where/to/store/outputs/dtm_{ORIGINALFILENAME}" dtm <- grid_terrain(ctg, tin())

#### Example 3 - Raster

Load a single big file too big to be actually loaded in memory. Process small chunks of 100 x 100 meters at a time. Returns a RasterLayer into R.

ctg <- catalog("path/to/las/files/bigfile.las") opt_chunk_size(ctg) <- 100 chm <- grid_canopy(ctg, p2r())

#### Example 4 - Tree detection

Load a catalog. Process small chunks of 200 x 200 meters at a time. Each chunk is loaded with an extra 20 m buffer. Use 3 cores. Returns a Shapefile into R.

ctg <- catalog("path/to/las/files/") opt_chunk_size(ctg) <- 200 opt_chunk_buffer(ctg) <- 20 opt_cores(ctg) <- 3 ttops <- tree_detection(ctg, lmf(5))

#### Example 5 - Decimate

This is forbidden. The output would be too big.

ctg <- catalog("path/to/las/files/") decimated <- lasfilterdecimate(ctg, homogenize(4))

#### Example 6 - Decimate

Load a catalog. Process small chunks of 500 x 500 meter sequentially. For each chunk write a laz file on disk named after the coordinates of the chunk. Returns a lightweight LAScatalog. Note that the original catalog has been retiled.

ctg <- catalog("path/to/las/files/") opt_chunk_size(ctg) <- 500 opt_output_files(ctg) <- "folder/where/to/store/outputs/project_{XLEFT}_{YBOTTOM}_decimated" opt_laz_compression(ctg) <- TRUE decimated <- lasfilterdecimate(ctg, homogenize(4))

#### Example 7 - Clip

Load a catalog. Load a shapefile of plot centers. Extract the plots. Returns a list of extracted point clouds in R.

ctg <- catalog("path/to/las/files/") shp <- shapefile("plot_center.shp") plots <- lasclip(ctg, shp, radius = 11.2)

#### Example 8 - Clip

Load a catalog. Load a shapefile of plot centers. Extract the plots and immediately write them into a file named after the coordinates of the plot and an attributes of the shapefile (here PLOTID if such an attribute exists in the shapefile). Returns a lightweight LAScatalog.

ctg <- catalog("path/to/las/files/") shp <- shapefile("plot_center.shp") opt_output_files(ctg) <- "folder/where/to/store/outputs/plot_{XCENTER}_{YCENTER}_{PLOTID}" plots <- lasclip(ctg, plot_centers, radius = 11.2)