Speed-up the computations on a LAScatalog

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 What takes time when processing a LAScatalog is not necessarily the computation itself but the time required to read the files. In fact the read time (i.e the time needed to load the data in R) might be much longer than the actual computation time. This vignette explains why and how to speed-up the computation by a factor of 2 to 8 by carefully preparing the catalog. Generic considerations on LAScatalog processing When processing a LAScatalog the area covered is divided into chunks that are then processed sequentially. In the following examples we present the case where chunks are equal to tiles, i.e. when processing each file sequentially. This is the common way to process data but not the only one. In any case, the explanation remains valid even when chunks are not equal to tiles. So each file is processed sequentially. For a given processed file, the content is read and loaded into R. In the following figure chunk number 42 is currently read and processed (in blue): bbox = as.numeric(ctg@data[42,1:4]) plot(ctg) graphics::rect(bbox[1], bbox[3], bbox[2], bbox[4], border = "black", col = "blue") But the current processed file is not the only one that needs to be read. To properly process the catalog without edge artifacts we need to load an extra buffer around the processed file (in red). bbox = as.numeric(ctg@data[42,1:4]) plot(ctg) graphics::rect(bbox[1]+200, bbox[3]+200, bbox[2]-200, bbox[4]-200, border = "black", col = "red") graphics::rect(bbox[2], bbox[4], bbox[1], bbox[3], border = "black", col = "blue") To load a buffer the processing engine must not only read the processed file but also the 8 neighbouring tiles (in red) to selectively read a small buffer around the processed file. Thus, for each processed file it is not one file that is read but nine. This is one of the reasons why the read time is far from negligible compared to the actual computation time. neighbourg = ctg@data[c(19,20,23, 41, 43, 44, 45, 47),1:4] plot(ctg) for (i in 1:nrow(neighbourg)) { bbox = as.numeric(neighbourg[i,]) graphics::rect(bbox[2], bbox[4], bbox[1], bbox[3], border = "black", col = "red") } To process a LAScatalog faster we need to read the files faster. Read a las file vs read a laz file A laz file is a strongly compressed las file. Reading a laz file is thus slower than reading a las file because it must be un-compressed on the fly at reading time. The following graph shows a benchmark of read time for a single file. system.time(readLAS(file.las)) system.time(readLAS(file.laz)) LASfile <- system.file("extdata", "Megaplot.laz", package = "lidR") las = readLAS(LASfile) f = tempfile(fileext = ".las") writeLAS(las, f) t1 = system.time(for (i in 1:10) readLAS(LASfile)) t2 = system.time(for (i in 1:10) readLAS(f)) t1 = t1[3] t2 = t2[3] t = c(t1,t2) format = c("laz", "las") X = data.frame(t, format) op <- graphics::par(mar = c(4,4,1,1) + 0.1) barplot(X$t/min(X$t), names.arg = X$format, col = "darkred", xlab = "File format", ylab = "Relative read time", asp = 1) graphics::par(op)

So let's assume that the total computation time is 1 unit of time divided into 0.25 units of actual processing time and 0.75 units of read time (which is a fairly reasonable ratio). We can divide the read time by 3, and thus have 0.25 units of read time and 0.25 units of computing time, which gives a computation time of 0.5 instead of 1. We can therefore speed-up the computation time by a factor of 2 by using the las format instead of laz. Obviously the gain is less significant for more computationally demanding processes.

So for faster computation users can opt for las files instead of laz files. Obviously, there are good reasons to use laz files instead of las files. The strong compression brought by the laz format has a lot of advantages for storing data. It is up to the user to choose a format by considering the trade-offs between space and computation time. This section explains how it works only to help users make a decision that best suits their needs.

Indexation of the points with lax files

Another way to speed-up the total computation time is to avoid reading all 8 neighbouring tiles to load a buffer. Instead, we can read only parts of the neighbouring tiles. The gain comes from the fact that we read only a small portion of the neighbouring files to extract the buffer, skipping most of the file contents. Indeed, the buffer usually corresponds to only a very small percentage of the actual contents of a file (equivalent to a few square meters).

This is made possible by indexing the las or laz files with lax files. A lax file is a tiny file associated with a las or laz file that spatially indexes the points to make faster spatial queries. This file type was created by Martin Isenburg in LAStools. For a better understanding of how it works one can refer to a talk given by Martin Isenburg about lasindex.

By adding lax files along with your las/laz files the buffer can be added around the processing file by only partially reading the 8 neigbouring files (in red)

neighbourg = ctg@data[c(19,20,23, 41, 43, 44, 45, 47),1:4] J = list(c(1200,1000,0,0), c(0,1000,0,0), c(0,1000,-1200,0), c(1200,0,0,0), c(1200,0,0,-1000), c(0,0,0,-1000), c(0,0,-1200,0), c(0,0,-1200,-1000)) plot(ctg) for (i in 1:nrow(neighbourg)) { j = J[[i]] bbox = as.numeric(neighbourg[i,]) graphics::rect(bbox[2] + j[1], bbox[4] + j[2], bbox[1] + j[3], bbox[3] + j[4], border = "black", col = "red") }

The best way to create a lax file is to use laxindex from LAStools. It is a free and open-source part of LAStools. If you cannot or do not want to use LAStools the lidR package has an (undocumented) internal function that creates lax files using the rlas package:

lidR:::catalog_laxindex(ctg)

The gain is really significant and allows an additional 2- to 3-fold saving in terms of read time, which significantly speeds up the computation time. Changing from laz to las format has a cost because it implies storing more data. However, using lax files provides a significant gain for free, so there is no incentive not to create lax files.

Changing the hardware

We demonstrated the importance of decreasing the time taken to read files to improve the overall computation time. The faster you read the files the faster you perform the computation because the read time is non-negligible. Opting for a fast SSD disk instead of a slow HDD disk may significantly speed-up the computation time independently of the power of your processor. Hardware matters!

Use more cores

This section comes at the end for a reason. lidR allows you to gain speed using multicore processing. However, users will not necessarily have a significant gain with multicore processing. Multicore is not a magic trick. Processing a LAScatalog implies reading a file on disk, a task that is not really parallelizable. As a consequence, multicore is not necessarily faster than single core. The multicore engine implemented in lidR reads several files at a time to process them in parallel. It uses more memory (four cores means four chunks read at one time) and the overhead of reading four files at a time may be more penalizing than the gain.

In conclusion, multicore processing may or may not speed-up the computation time and it has a significant memory cost.

Benchmarks

The following are benchmarks for some functions

Simple canopy height model

A simple point-to-raster canopy height model using 25 files of 150 x 150 m with 30 pts/m² on a laptop with an SSD and an intel core i7 processor.

chm <- grid_canopy(ctg, 1, p2r())
Cores = c(1,1,1,1,2,2) Format = c("laz", "laz + lax", "las", "las + lax", "laz + lax", "las + lax") Runtime = c("40 sec", "20 sec", "10 sec", "7 sec", "20 sec", "5 sec") Timing = data.frame(Format, Cores, Runtime) knitr::kable(Timing)

Here we found an almost 8-fold increase in speed simply by changing the file types. With 2 cores we gained an additional 30%, reaching an actual 8-fold speed-up compared to the basic use of laz files. One can try laz files only + 32 cores on a supercomputer with a lot of RAM, but the gain will probably not reach an 8-fold speed-up!

Area-based Approach

Computation of a single metric on 360 files of 1 x 1 km with 3 pts/m² (~300 km² and 900 millions points) on a laptop with an SSD and an intel core i7 processor.

hmean <- grid_metric(ctg, mean(Z), 20)
Cores = c(1,1,1,4) Format = c("laz", "las", "las + lax", "las + lax") Runtime = c("45 min", "15 min", "8 min", "7 min") Timing = data.frame(Format, Cores, Runtime) knitr::kable(Timing)

Here we found an almost 6-fold speed-up by changing only the file types. With 4 cores we gained an additional 20%, reaching almost a 7-fold speed-up compared to the basic use of laz files. Using 2 cores instead of 4 is likely to be faster. Again, one can try laz files only + 32 cores on a supercomputer with a lot of RAM, but the gain will probably not reach an 8-fold speed-up!

Clip ground inventories

Extraction of 140 plots of 12 m radius from 1 x 1 km files with 3 pts/m² on a laptop with an SSD and an intel core i7 processor.

# Streaming way opt_output_files(ctg) <- "/path/to/output/file_{ID}" new_ctg <- lasclip(ctg, shapefile, radius = 12) # Load in memory opt_output_files(ctg) <- "" new_ctg <- lasclip(ctg, shapefile, radius = 12)
Cores = c(1,4,1,4, 1,4,1,4) Method = c("Streaming", "Streaming", "Streaming", "Streaming", "Memory", "Memory", "Memory", "Memory") Format = c("las", "las", "las + lax", "las + lax", "las", "las", "las + lax", "las + lax") Runtime = c("45 sec", "35 sec", "4 sec", "6 sec", "45 sec", "35 sec", "4 sec", "17 sec") Timing = data.frame(Format, Method, Cores, Runtime) knitr::kable(Timing)

Here using 4 cores instead of 1 considerably increases the computation time. Indeed the runtime is very small and the overhead of the parallelization does not compensate for the potential gain.