Package 'lidR'

Description

lidR provides a set of tools to manipulate airborne LiDAR data in forestry contexts. The package works with .las or .laz files. The toolbox includes algorithms for DSM, CHM, DTM, ABA, normalisation, tree detection, tree segmentation, tree delineation, colourization, validation and other tools, as well as a processing engine to process broad LiDAR coverage split into many files.

Details

To learn more about lidR, start with the vignettes: browseVignettes(package = "lidR"). Users can also find unofficial supplementary documentation in the lidR book. To ask "how to" questions please ask on gis.stackexchange.com with the tag lidr.

Package options

lidR.progress

Several functions have a progress bar for long operations (but not all). Should lengthy operations show a progress bar? Default: TRUE

lidR.progress.delay

The progress bar appears only for long operations. After how many seconds of computation does the progress bar appear? Default: 2

lidR.raster.default

The functions that return a raster are raster agnostic meaning that they work either with rasters from packages 'raster', 'stars' or 'terra'. By default they return rasters from 'stars'. Can be one of "raster", "stars" or "terra". Default: "terra"

lidR.check.nested.parallelism

The catalog processing engine (catalog_apply) checks the parallel strategy chosen by the user and verify if C++ parallelization with OpenMP should be disabled to avoid nested parallel loops. Default: TRUE. If FALSE the catalog processing engine will not check for nested parallelism and will respect the settings of set_lidr_threads.

Author(s)

Maintainer: Jean-Romain Roussel [email protected] [copyright holder]

Authors:

• David Auty (Reviews the documentation) [contributor]

Other contributors:

• Florian De Boissieu (Fixed bugs and improved catalog features) [contributor]

• Andrew Sánchez Meador (Implemented wing2015() for segment_snags()) [contributor]

• Bourdon Jean-François (Contributed to Roussel2020() for track_sensor()) [contributor]

• Gatziolis Demetrios (Implemented Gatziolis2019() for track_sensor()) [contributor]

• Leon Steinmeier (Contributed to parallelization management) [contributor]

• Stanislaw Adaszewski (Author of the C++ concaveman code) [copyright holder]

• Benoît St-Onge (Author of the 'chm_prep' function) [copyright holder]

See Also

Useful links:

• https://github.com/r-lidar/lidR

• Report bugs at https://github.com/r-lidar/lidR/issues

Description

A LAS object represents a las file in R. According to the LAS specifications a las file contains a core of defined attributes, such as XYZ coordinates, intensity, return number, and so on, for each point. It is possible to add supplementary attributes.

Usage

add_attribute(las, x, name)

add_lasattribute(las, x, name, desc)

add_lasattribute_manual(
  las,
  x,
  name,
  desc,
  type,
  offset = NULL,
  scale = NULL,
  NA_value = NULL
)

add_lasrgb(las, R, G, B)

add_lasnir(las, NIR)

remove_lasattribute(las, name)

Arguments

Details

Users cannot assign names that are the same as the names of the core attributes. These functions are dedicated to adding data that are not part of the LAS specification. For example, add_lasattribute(las, x, "R") will fail because R is a name reserved for the red channel of a .las file that contains RGB attributes. Use add_lasrgb instead.

add_attribute

Simply adds a new column in the data but does not update the header. Thus the LAS object is not strictly valid. These data will be temporarily usable at the R level but will not be written in a las file with writeLAS.

add_lasattribute

Does the same as add_attribute but automatically updates the header of the LAS object. Thus, the LAS object is valid and the new data is considered as "extra bytes". This new data will be written in a las file with writeLAS.

add_lasattribute_manual

Allows the user to manually write all the extra bytes metadata. This function is reserved for experienced users with a good knowledge of the LAS specifications. The function does not perform tests to check the validity of the information. When using add_lasattribute and add_lasattribute_manual, x can only be of type numeric, (integer or double). It cannot be of type character or logical as these are not supported by the LAS specifications. The types that are supported in lidR are types 0 to 10 (Table 24 on page 25 of the specification). Types greater than 10 are not supported.

add_lasrgb

Adds 3 columns named RGB and updates the point format of the LAS object for a format that supports RGB attributes. If the RGB values are ranging from 0 to 255 they are automatically scaled on 16 bits.

Value

An object of class LAS

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las <- readLAS(LASfile, select = "xyz")

print(las)
print(header(las))

x <- 1:30

las <- add_attribute(las, x, "mydata")
print(las)        # The las object has a new attribute called "mydata"
print(header(las)) # But the header has not been updated. This new data will not be written

las <- add_lasattribute(las, x, "mydata2", "A new data")
print(las)        # The las object has a new attribute called "mydata2"
print(header(las)) # The header has been updated. This new data will be written

# Optionally if the data is already in the LAS object you can update the header skipping the
# parameter x
las <- add_attribute(las, x, "newattr")
las <- add_lasattribute(las, name = "newattr", desc = "Amplitude")
print(header(las))

# Remove an extra bytes attribute
las <- remove_lasattribute(las, "mydata2")
print(las)
print(header(las))

las <- remove_lasattribute(las, "mydata")
print(las)
print(header(las))

Description

template_metrics() computes a series of user-defined descriptive statistics for a LiDAR dataset within each element of a template. Depending on the template it can be for each pixel of a raster (area-based approach), or each polygon, or each segmented tree, or on the whole point cloud. Other functions are convenient and simplified wrappers around template_metrics() and are expected to be the actual functions used. See Details and Examples.

Usage

cloud_metrics(las, func, ...)

crown_metrics(
  las,
  func,
  geom = "point",
  concaveman = c(3, 0),
  attribute = "treeID",
  ...
)

hexagon_metrics(las, func, area = 400, ...)

pixel_metrics(las, func, res = 20, start = c(0, 0), ...)

plot_metrics(las, func, geometry, ..., radius)

polygon_metrics(las, func, geometry, ...)

template_metrics(las, func, template, filter = NULL, by_echo = "all", ...)

voxel_metrics(las, func, res = 1, ..., all_voxels = FALSE)

Arguments

Details

pixel_metrics

Area-based approach. Computes metrics in a square tessellation. The output is a raster.

hexagon_metrics

Computes metrics in an hexagon tessellation. The output is a sf/sfc_POLYGON

plot_metrics

Computes metrics for each plot of a ground inventory by 1. clipping the plot inventories with clip_roi, 2. computing the user's metrics for each plot with cloud_metrics, and 3. combining spatial data and metrics into one data.frame ready for statistical modelling with cbind. The output is of the class of the input.

cloud_metrics

Computes a series of user-defined descriptive statistics for an entire point cloud. The output is a list

crown_metrics

Once the trees are segmented, i.e. attributes exist in the point cloud that reference each tree, computes a set of user-defined descriptive statistics for each individual tree. The output can be spatial points or spatial polygons (sf/sfc_POINT or sf/sfc_POLYGON)

voxel_metrics

Is a 3D version of pixel_metrics. It creates a 3D matrix of voxels with a given resolution. It creates a voxel from the cloud of points if there is at least one point. The output is a data.frame

point_metrics

Is a bit more complex and is documented in point_metrics

Value

Depends on the function, the template and the number of metrics. Can be a RasterLayer, a RasterBrick, a stars, a SpatRaster a sf/sfc, a list, a SpatialPolygonDataFrame, or a data.table. Functions are supposed to return an object that is best suited for storing the level of regularization needed.

Parameter func

The function to be applied to each cell is a classical function (see examples) that returns a labelled list of metrics. For example, the following function f is correctly formed.

f = function(x) {list(mean = mean(x), max = max(x))}

And could be applied either on the Z coordinates or on the intensities. These two statements are valid:

pixel_metrics(las, f(Z), res = 20)
voxel_metrics(las, f(Intensity), res = 2)

The following existing functions allow the user to compute some predefined metrics: stdmetrics entropy, VCI, LAD. But usually users must write their own functions to create metrics. template_metrics will dispatch the point cloud in the user's function.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile, filter = "-keep_random_fraction 0.5")
col <- sf::sf.colors(15)
fun1 <- ~list(maxz = max(Z))
fun2 <- ~list(q85 = quantile(Z, probs = 0.85))

set_lidr_threads(1) ; data.table::setDTthreads(1) # for cran only

# ================
# CLOUD METRICS
# ================

cloud_metrics(las, .stdmetrics_z)

# ================
# PIXEL METRICS
# ================

m <- pixel_metrics(las, fun1, 20)
#plot(m, col = col)

# ================
# PLOT METRICS
# ================

shpfile <- system.file("extdata", "efi_plot.shp", package="lidR")
inventory <- sf::st_read(shpfile, quiet = TRUE)
inventory # contains an ID and a Value Of Interest (VOI) per plot

m <- plot_metrics(las, fun2, inventory, radius = 11.28)
#plot(header(las))
#plot(m["q85"], pch = 19, cex = 3, add = TRUE)


# Works with polygons as well
inventory <- sf::st_buffer(inventory, 11.28)
#plot(header(las))
#plot(sf::st_geometry(inventory), add = TRUE)
m <- plot_metrics(las, .stdmetrics_z, inventory)
#plot(m["zq85"], pch = 19, cex = 3, add = TRUE)


# ================
# VOXEL METRICS
# ================

m <- voxel_metrics(las, length(Z), 8)
m <- voxel_metrics(las, mean(Intensity), 8)
#plot(m, color = "V1", colorPalette = heat.colors(50), trim = 60)
#plot(m, color = "V1", colorPalette = heat.colors(50), trim = 60, voxel = TRUE)

# ================
# CROWN METRICS
# ================

# Already tree-segmented point cloud
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
trees <- readLAS(LASfile, filter = "-drop_z_below 0")

metrics <- crown_metrics(trees, .stdtreemetrics)
#plot(metrics["Z"], pch = 19)

metrics <- crown_metrics(trees, .stdtreemetrics, geom = "convex")
#plot(metrics["Z"])

metrics <- crown_metrics(trees, .stdtreemetrics, geom = "bbox")
#plot(metrics["Z"])


metrics <- crown_metrics(trees, .stdtreemetrics, geom = "concave")
#plot(metrics["Z"])

# ================
# ARGUMENT FILTER
# ================

# Compute using only some points: basic
first = filter_poi(las, ReturnNumber == 1)
metrics = pixel_metrics(first, mean(Z), 20)

# Compute using only some points: optimized
# faster and uses less memory. No intermediate object
metrics = pixel_metrics(las, mean(Z), 20, filter = ~ReturnNumber == 1)

# Compute using only some points: best
# ~50% faster and uses ~10x less memory
las = readLAS(LASfile, filter = "-keep_first")
metrics = pixel_metrics(las, mean(Z), 20)

# ================
# ARGUMENT BY_ECHO
# ================

func = ~list(avgI = mean(Intensity))
echo = c("all", "first","multiple")

# func defines one metric but 3 are computed respectively for: (1) all echo types,
# (2) for first returns only and (3) for multiple returns only
metrics <- pixel_metrics(las, func, 20, by_echo = echo)
#plot(metrics, col = heat.colors(25))

cloud_metrics(las, func, by_echo = echo)

## Not run: 
# ================
# TEMPLATE METRICS
# ================

# a raster as template
template <- raster::raster(extent(las), nrow = 15, ncol = 15)
raster::crs(template) <- crs(las)
m <- template_metrics(las, fun1, template)
#plot(m, col = col)

# a sfc_POLYGON as template
sfc <- sf::st_as_sfc(st_bbox(las))
template <- sf::st_make_grid(sfc, cellsize = 20, square = FALSE)
m <- template_metrics(las, fun1, template)
#plot(m)

# a bbox as template
template <- st_bbox(las) + c(50,30,-50,-70)
plot(sf::st_as_sfc(st_bbox(las)), col = "gray")
plot(sf::st_as_sfc(template), col = "darkgreen", add = TRUE)
m <- template_metrics(las, fun2, template)
print(m)

# ================
# CUSTOM METRICS
# ================

# Define a function that computes custom metrics
# in an R&D perspective.
myMetrics = function(z, i) {
  metrics = list(
     zwimean = sum(z*i)/sum(i), # Mean elevation weighted by intensities
     zimean  = mean(z*i),       # Mean products of z by intensity
     zsqmean = sqrt(mean(z^2))) # Quadratic mean

   return(metrics)
}

# example with a stars template
template <- stars::st_as_stars(st_bbox(las), dx = 10, dy = 10)
m <- template_metrics(las, myMetrics(Z, Intensity), template)
#plot(m, col = col)

## End(Not run)

Description

Functions to construct, coerce and check for both kinds of R lists.

Usage

## S3 method for class 'LASheader'
as.list(x, ...)

Arguments

Description

A set of global variables corresponding to the point classification defined by the ASPRS for the LAS format. Instead of remembering the classification table of the specification it is possible to use one of these global variables.

Usage

LASNONCLASSIFIED

LASUNCLASSIFIED

LASGROUND

LASLOWVEGETATION

LASMEDIUMVEGETATION

LASHIGHVEGETATION

LASBUILDING

LASLOWPOINT

LASKEYPOINT

LASWATER

LASRAIL

LASROADSURFACE

LASWIREGUARD

LASWIRECONDUCTOR

LASTRANSMISSIONTOWER

LASBRIGDE

LASNOISE

Format

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

An object of class integer of length 1.

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las = readLAS(LASfile)
las2 = filter_poi(las, Classification %in% c(LASGROUND, LASWATER))

print(LASGROUND)

Description

This function gives users access to the LAScatalog processing engine. It allows the application of a user-defined routine over a collection of LAS/LAZ files. The LAScatalog processing engine is explained in the LAScatalog class

catalog_apply() is the core of the lidR package. It drives every single function that can process a LAScatalog. It is flexible and powerful but also complex. catalog_map() is a simplified version of catalog_apply() that suits for 90% of use cases.

catalog_sapply() is a previous attempt to provide simplified version of catalog_apply(). Use catalog_map() instead.

Usage

catalog_apply(ctg, FUN, ..., .options = NULL)

catalog_sapply(ctg, FUN, ..., .options = NULL)

catalog_map(ctg, FUN, ..., .options = NULL)

Arguments

Edge artifacts

It is important to take precautions to avoid 'edge artifacts' when processing wall-to-wall tiles. If the points from neighbouring tiles are not included during certain processes, this could create 'edge artifacts' at the tile boundaries. For example, empty or incomplete pixels in a rasterization process, or dummy elevations in a ground interpolation. The LAScatalog processing engine provides internal tools to load buffered data 'on-the-fly'. catalog_map() takes care of removing automatically the results computed in the buffered area to avoid unexpected output with duplicated entries or conflict between values computed twice. It does that in predefined way that may not suit all cases. catalog_apply() does not remove the buffer and leave users free to handle this in a custom way. This is why catalog_apply() is more complex but gives more freedom to build new applications.

Buffered data

The LAS objects loaded in theses functions have a special attribute called 'buffer' that indicates, for each point, if it comes from a buffered area or not. Points from non-buffered areas have a 'buffer' value of 0, while points from buffered areas have a 'buffer' value of 1, 2, 3 or 4, where 1 is the bottom buffer and 2, 3 and 4 are the left, top and right buffers, respectively. This allows for filtering of buffer points if required.

Function template

The parameter FUN of catalog_apply expects a function with a first argument that will be supplied automatically by the LAScatalog processing engine. This first argument is a LAScluster. A LAScluster is an internal undocumented class but the user needs to know only three things about this class:

• It represents a chunk of the file collection

• The function readLAS can be used with a LAScluster

• The function extent or st_bbox can be used with a LAScluster and they return the bounding box of the chunk without the buffer. It must be used to clip the output and remove the buffered region (see examples).

A user-defined function must be templated like this:

myfun <- function(chunk, ...) {
   # Load the chunk + buffer
   las <- readLAS(chunk)
   if (is.empty(las)) return(NULL)

   # do something
   output <- do_something(las, ...)

   # remove the buffer of the output
   bbox <- bbox(chunk)
   output <- remove_buffer(output, bbox)
   return(output)
}

The line if (is.empty(las)) return(NULL) is important because some clusters (chunks) may contain 0 points (we can't know this before reading the file). In this case an empty point cloud with 0 points is returned by readLAS() and this may fail in subsequent code. Thus, exiting early from the user-defined function by returning NULL indicates to the internal engine that the chunk was empty.

catalog_map is much simpler (but less versatile). It automatically takes care of reading the chunk and checks if the point cloud is empty. It also automatically crop the buffer. The way it crops the buffer suits for most cases but for some special cases it may be advised to handle this in a more specific way i.e. using catalog_apply(). For catalog_map() the first argument is a LAS and the template is:

myfun <- function(las, ...) {
   # do something
   output <- do_something(las, ...)
   return(output)
}

.options

Users may have noticed that some lidR functions throw an error when the processing options are inappropriate. For example, some functions need a buffer and thus buffer = 0 is forbidden. Users can add the same constraints to protect against inappropriate options. The .options argument is a list that allows users to tune the behaviour of the processing engine.

• drop_null = FALSE: not intended to be used by regular users. The engine does not remove NULL outputs

• need_buffer = TRUE: the function complains if the buffer is 0.

• need_output_file = TRUE the function complains if no output file template is provided.

• raster_alignment = ... the function checks the alignment of the chunks. This option is important if the output is a raster. See below for more details.

• automerge = TRUE by default the engine returns a list with one item per chunk. If automerge = TRUE, it tries to merge the outputs into a single object: a ⁠Raster*⁠, a ⁠Spatial*⁠, a ⁠LAS*⁠, an sf, a stars similarly to other functions of the package. This is a fail-safe option so in the worst case, if the merging fails, the list is returned. This is activated by catalog_map making it simpler.

• autoread = TRUE. Introduced in v3.0.0 this option enables to get rid of the first steps of the function i.e readLAS() and ⁠if (is.empty())⁠. In this case the function must take two objects as input, first a LAS object and second a bbox from sf. This is activated by catalog_map making it simpler.

• autocrop = TRUE. Introduced in v4.0.0 this option enables to get rid of the buffer crop steps of the function i.e something <- remove_buffer(something, bbox). In this case the function must take one LAS object as input. This is activated by catalog_map making it simpler.

When the function FUN returns a raster it is important to ensure that the chunks are aligned with the raster to avoid edge artifacts. Indeed, if the edge of a chunk does not correspond to the edge of the pixels, the output will not be strictly continuous and will have edge artifacts (that might not be visible). Users can check this with the options raster_alignment, that can take the resolution of the raster as input, as well as the starting point if needed. The following are accepted:

# check if chunks are aligned with a raster of resolution 20
raster_alignment = 20
raster_alignment = list(res = 20)

# check if chunks are aligned with a raster of resolution 20
# that starts at (0,10)
raster_alignment = list(res = 20, start = c(0,10))

Examples

# More examples might be avaible in the official lidR vignettes or
# on the github book <https://jean-romain.github.io/lidRbook/>

## =========================================================================
## Example 1: detect all the tree tops over an entire catalog
## (this is basically a reproduction of the existing function 'locate_trees')
## =========================================================================

# 1. Build the user-defined function that analyzes each chunk of the catalog.
# The function's first argument is a LAScluster object. The other arguments can be freely
# chosen by the users.
my_tree_detection_method <- function(chunk, ws)
{
  # The chunk argument is a LAScluster object. The users do not need to know how it works.
  # readLAS will load the region of interest (chunk) with a buffer around it, taking advantage of
  # point cloud indexation if possible. The filter and select options are propagated automatically
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)

  # Find the tree tops using a user-developed method
  # (here simply a LMF for the example).
  ttops <- locate_trees(las, lmf(ws))

  # ttops is an sf object that contains the tree tops in our region of interest
  # plus the trees tops in the buffered area. We need to remove the buffer otherwise we will get
  # some trees more than once.
  bbox  <- st_bbox(chunk)
  ttops <- sf::st_crop(ttops, bbox)
  return(ttops)
}

# 2. Build a collection of file
# (here, a single file LAScatalog for the purposes of this simple example).
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
ctg <- readLAScatalog(LASfile)
plot(ctg)

# 3. Set some processing options.
# For this small single file example, the chunk size is 100 m + 10 m of buffer
opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg)   <- 100            # Small because this is a dummy example.
opt_chunk_alignment(ctg) <- c(-50, -35) # Align such as it creates 2 chunks only.
opt_select(ctg)       <- "xyz"          # Read only the coordinates.
opt_filter(ctg)       <- "-keep_first"  # Read only first returns.

# 4. Apply a user-defined function to take advantage of the internal engine
opt    <- list(need_buffer = TRUE,   # catalog_apply will throw an error if buffer = 0
               automerge   = TRUE)   # catalog_apply will merge the outputs into a single object
output <- catalog_apply(ctg, my_tree_detection_method, ws = 5, .options = opt)

plot(output)


## =========================================================================
## Example 1: simplified. There is nothing that requires special data
## manipulation in the previous example. Everything can be handled automatically
##=========================================================================

# 1. Build the user-defined function that analyzes a point cloud.
my_tree_detection_method <- function(las, ws)
{
  # Find the tree tops using a user-developed method
  # (here simply a LMF for the example).
  ttops <- locate_trees(las, lmf(ws))
  return(ttops)
}

# 2. Build a project
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
ctg <- readLAScatalog(LASfile)
plot(ctg)

# 3. Set some processing options.
# For this dummy example, the chunk size is 100 m and the buffer is 10 m
opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg)   <- 100            # small because this is a dummy example.
opt_chunk_alignment(ctg) <- c(-50, -35) # Align such as it creates 2 chunks only.
opt_select(ctg)       <- "xyz"          # Read only the coordinates.
opt_filter(ctg)       <- "-keep_first"  # Read only first returns.

# 4. Apply a user-defined function to take advantage of the internal engine
opt    <- list(need_buffer = TRUE)   # catalog_apply will throw an error if buffer = 0
output <- catalog_map(ctg, my_tree_detection_method, ws = 5, .options = opt)


## Not run: 
## ===================================================
## Example 2: compute a rumple index on surface points
## ===================================================

rumple_index_surface = function(las, res)
{
  las    <- filter_surfacepoints(las, 1)
  rumple <- pixel_metrics(las, ~rumple_index(X,Y,Z), res)
  return(rumple)
}

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg <- readLAScatalog(LASfile)

opt_chunk_buffer(ctg) <- 1
opt_chunk_size(ctg)   <- 140     # small because this is a dummy example.
opt_select(ctg)       <- "xyz"   # read only the coordinates.

opt    <- list(raster_alignment = 20) # catalog_apply will adjust the chunks if required
output <- catalog_map(ctg, rumple_index_surface, res = 20, .options = opt)

plot(output, col = height.colors(25))

## End(Not run)

Description

Computes the polygon that encloses the points. It reads all the files one by one and computes a concave hull using the st_concave_hull function. When all the hulls are computed it updates the LAScatalog to set the true polygons instead of the bounding boxes.

Usage

catalog_boundaries(ctg, ...)

Arguments

Value

A LAScatalog with true boundaries

Non-supported LAScatalog options

The options 'select', 'output files', 'chunk size', 'chunk buffer', 'chunk alignment' are not supported and not respected in 'catalog_boundaries*' because the function must always process by file, without buffer and knows which attributes to load.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg <- readLAScatalog(LASfile, filter = "-drop_z_below 0.5")
ctg2 <- catalog_boundaries(ctg, concavity = 1, length_threshold = 15)
plot(ctg)
plot(ctg2, add = TRUE)

Description

Subset a LAScatalog interactively using the mouse. Subset a LAScatalog with a spatial object to keep only the tiles of interest. It can be used to select tiles of interest that encompass spatial objects.

Usage

catalog_intersect(
  ctg,
  y,
  ...,
  subset = c("subset", "flag_unprocessed", "flag_processed")
)

catalog_select(
  ctg,
  mapview = TRUE,
  subset = c("subset", "flag_unprocessed", "flag_processed")
)

Arguments

Value

A LAScatalog object

Examples

## Not run: 
ctg = readLAScatalog("<Path to a folder containing a set of .las files>")
new_ctg = catalog_select(ctg)

## End(Not run)

Description

Splits or merges files to reshape the original files collection (.las or .laz) into smaller or larger files. It also enables the addition or removal of a buffer around the tiles. Internally, the function reads and writes the chunks defined by the internal processing options of a LAScatalog. Thus, the function is flexible and enables the user to retile the dataset, retile while adding or removing a buffer (negative buffers are allowed), or optionally to compress the data by retiling without changing the pattern but by changing the format (las/laz). This function is does not load the point cloud into R memory It streams from input file(s) to output file(s) and can be applied to large point-cloud with low memory computer.

Note that this function is not actually very useful because lidR manages everything (clipping, processing, buffering, ...) internally using the proper options. Thus, retiling may be useful for working in other software, for example, but not in lidR

Usage

catalog_retile(ctg)

Arguments

Value

A new LAScatalog object

Non-supported LAScatalog options

The option select is not supported and not respected because it always preserves the file format and all the attributes. select = "*" is imposed internally.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)
plot(ctg)

# Create a new set of 200 x 200 m.las files with first returns only
opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 200
opt_filter(ctg) <- "-keep_first"
opt_chunk_alignment(ctg) <- c(275, 90)
opt_output_files(ctg) <- paste0(tempdir(), "/retile_{XLEFT}_{YBOTTOM}")

# preview the chunk pattern
plot(ctg, chunk = TRUE)

newctg = catalog_retile(ctg)

plot(newctg)

# Create a new set of 200 x 200 m.las files
# but extended with a 50 m buffer in the folder

opt_chunk_buffer(ctg) <- 25
opt_chunk_size(ctg) <- 200
opt_filter(ctg) <- ""
opt_chunk_alignment(ctg) <- c(275, 90)
opt_output_files(ctg) <- paste0(tempdir(), "/{XLEFT}_{YBOTTOM}_buffered")
newctg = catalog_retile(ctg)

plot(newctg)


## Not run: 
# Create a new set of compressed .laz file equivalent to the original, keeping only
# first returns above 2 m

opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 0
opt_laz_compression(ctg) <- TRUE
opt_filter(ctg) <- "-keep_first -drop_z_below 2"
opt_output_files(ctg) <- paste0(tempdir(), "/{ORIGINALFILENAME}_first_2m")
newctg = catalog_retile(ctg)

## End(Not run)

Description

Classify points that meet some criterion and/or that belong in a region of interest. The functions updates the attribute Classification of the LAS object according to las specifications

Usage

classify_ground(las, algorithm, last_returns = TRUE)

classify_noise(las, algorithm)

classify_poi(
  las,
  class,
  poi = NULL,
  roi = NULL,
  inverse_roi = FALSE,
  by_reference = FALSE
)

Arguments

Details

classify_noise

Classify points as 'noise' (outliers) with several possible algorithms. lidR has: sor, ivf. The points classified as 'noise' are assigned a value of 18.

classify_ground

Classify points as 'ground' with several possible algorithms. lidR has pmf, csf and mcc. The points classified as 'ground' are assigned a value of 2

classify_poi

Classify points that meet some logical criterion and/or that belong in a region of interest with class of choice.

Non-supported LAScatalog options

The option select is not supported and not respected because it always preserves the file format and all the attributes. select = "*" is imposed internally.

Examples

# ===============
# Classify ground
# ===============

if (require(RCSF, quietly = TRUE))
{
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn", filter = "-inside 273450 5274350 273550 5274450")

# (Parameters chosen mainly for speed)
mycsf <- csf(TRUE, 1, 1, time_step = 1)
las <- classify_ground(las, mycsf)
#plot(las, color = "Classification")
}

# ===============
# Classify noise
# ===============

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, filter = "-inside 273450 5274350 273550 5274450")

# Add 20 artificial outliers
set.seed(314)
id = round(runif(20, 0, npoints(las)))
set.seed(42)
err = runif(20, -50, 50)
las$Z[id] = las$Z[id] + err

# Using IVF
las <- classify_noise(las, ivf(5,2))
#plot(las, color = "Classification")

# Remove outliers using filter_poi()
las_denoise <- filter_poi(las, Classification != LASNOISE)

# ===============
# Classify POI
# ===============

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
shp <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")

las  <- readLAS(LASfile, filter = "-keep_random_fraction 0.1")
lake <- sf::st_read(shp, quiet = TRUE)

# Classifies the points that are NOT in the lake and that are NOT ground points as class 5
poi <- ~Classification != LASGROUND
las <- classify_poi(las, LASHIGHVEGETATION, poi = poi, roi = lake, inverse = TRUE)

# Classifies the points that are in the lake as class 9
las <- classify_poi(las, LASWATER, roi = lake, inverse = FALSE)

#plot(las, color = "Classification")

Description

Clip points within a given region of interest (ROI) from a point cloud (LAS object) or a collection of files (LAScatalog object).

Usage

clip_roi(las, geometry, ...)

clip_rectangle(las, xleft, ybottom, xright, ytop, ...)

clip_polygon(las, xpoly, ypoly, ...)

clip_circle(las, xcenter, ycenter, radius, ...)

clip_transect(las, p1, p2, width, xz = FALSE, ...)

Arguments

Value

If the input is a LAS object: an object of class LAS, or a list of LAS objects if the query implies several regions of interest.

If the input is a LAScatalog object: an object of class LAS, or a list of LAS objects if the query implies several regions of interest, or a LAScatalog if the queries are immediately written into files without loading anything in R.

Non-supported LAScatalog options

The option ⁠chunk size⁠, buffer, ⁠chunk alignment⁠ and select are not supported by ⁠clip_*⁠ because they are meaningless in this context.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")

# Load the file and clip the region of interest
las = readLAS(LASfile, select = "xyz", filter = "-keep_first")
subset1 = clip_rectangle(las, 684850, 5017850, 684900, 5017900)

# Do not load the file(s), extract only the region of interest
# from a bigger dataset
ctg = readLAScatalog(LASfile, progress = FALSE, filter = "-keep_first")
subset2 = clip_rectangle(ctg, 684850, 5017850, 684900, 5017900)

# Extract all the polygons from a shapefile
f <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")
lakes <- sf::st_read(f, quiet = TRUE)
subset3 <- clip_roi(las, lakes)

# Extract the polygons for a catalog, write them in files named
# after the lake names, do not load anything in R
opt_output_files(ctg) <- paste0(tempfile(), "_{LAKENAME_1}")
new_ctg = clip_roi(ctg, lakes)
plot(new_ctg)

# Extract a transect
p1 <- c(684800, y = 5017800)
p2 <- c(684900, y = 5017900)
tr1 <- clip_transect(las, p1, p2, width = 4)

## Not run: 
plot(subset1)
plot(subset2)
plot(subset3)

plot(tr1, axis = TRUE, clear_artifacts = FALSE)

## End(Not run)

Description

Reduce the number of points using several possible algorithms.

Usage

decimate_points(las, algorithm)

Arguments

Value

If the input is a LAS object, returns a LAS object. If the input is a LAScatalog, returns a LAScatalog.

Non-supported LAScatalog options

The option 'select' is not supported and not respected because it always preserves the file format and all the attributes. 'select = "*"' is imposed internally.
The options 'chunk buffer' is not supported and not respected because it is not needed.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile, select = "xyz")

# Select points randomly to reach an overall density of 1
thinned1 <- decimate_points(las, random(1))
#plot(rasterize_density(las))
#plot(rasterize_density(thinned1))

# Select points randomly to reach an homogeneous density of 1
thinned2 <- decimate_points(las, homogenize(1,5))
#plot(rasterize_density(thinned2))

# Select the highest point within each pixel of an overlayed grid
thinned3 = decimate_points(las, highest(5))
#plot(thinned3)

Description

These functions are provided for compatibility with older versions of lidR but are deprecated. They will progressively print a message, throw a warning and eventually be removed. The links below point to the documentation of the new names. In version 4 they now throw an error. In version 4.1 they ill be removed definitively.

lasadd lascheck lasclip lasdetectshape lasfilter lasfiltersurfacepoints lasflightline lasground lasmergespatial lasnormalize laspulse lasrangecorrection lasflightline lasreoffset lasrescale lasscanlines lassmooth lassnags lastrees lasvoxelize sensor_tracking tree_detection tree_hull

Usage

lascheck(las)

lasclip(las, geometry, ...)

lasclipRectangle(las, xleft, ybottom, xright, ytop, ...)

lasclipPolygon(las, xpoly, ypoly, ...)

lasclipCircle(las, xcenter, ycenter, radius, ...)

lasdetectshape(las, algorithm, attribute = "Shape", filter = NULL)

lasfilter(las, ...)

lasfilterfirst(las)

lasfilterfirstlast(las)

lasfilterfirstofmany(las)

lasfilterground(las)

lasfilterlast(las)

lasfilternth(las, n)

lasfiltersingle(las)

lasfilterdecimate(las, algorithm)

lasfilterduplicates(las)

lasfiltersurfacepoints(las, res)

lasground(las, algorithm, last_returns = TRUE)

laspulse(las)

lasflightline(las, dt = 30)

lasscanline(las)

lasmergespatial(las, source, attribute = NULL)

lasnormalize(
  las,
  algorithm,
  na.rm = FALSE,
  use_class = c(2L, 9L),
  ...,
  add_lasattribute = FALSE
)

lasunnormalize(las)

lasrangecorrection(
  las,
  sensor,
  Rs,
  f = 2.3,
  gpstime = "gpstime",
  elevation = "Z"
)

lasrescale(las, xscale, yscale, zscale)

lasreoffset(las, xoffset, yoffset, zoffset)

lassmooth(
  las,
  size,
  method = c("average", "gaussian"),
  shape = c("circle", "square"),
  sigma = size/6
)

lasunsmooth(las)

lassnags(las, algorithm, attribute = "snagCls")

lastrees(las, algorithm, attribute = "treeID", uniqueness = "incremental")

lasadddata(las, x, name)

lasaddextrabytes(las, x, name, desc)

lasaddextrabytes_manual(
  las,
  x,
  name,
  desc,
  type,
  offset = NULL,
  scale = NULL,
  NA_value = NULL
)

lasremoveextrabytes(las, name)

lasvoxelize(las, res)

sensor_tracking(
  las,
  interval = 0.5,
  pmin = 50,
  extra_check = TRUE,
  thin_pulse_with_time = 0.001
)

tree_detection(las, algorithm)

tree_hulls(
  las,
  type = c("convex", "concave", "bbox"),
  concavity = 3,
  length_threshold = 0,
  func = NULL,
  attribute = "treeID"
)

hexbin_metrics(...)

filter_surfacepoints(las, res)

## S3 method for class 'LAS'
filter_surfacepoints(las, res)

## S3 method for class 'LAScatalog'
filter_surfacepoints(las, res)

Arguments

Description

This function is made to be used in rasterize_canopy. It implements the pit-free algorithm developed by Khosravipour et al. (2014), which is based on the computation of a set of classical triangulations at different heights (see references). The subcircle tweak replaces each point with 8 points around the original one. This allows for virtual 'emulation' of the fact that a lidar point is not a point as such, but more realistically a disc. This tweak densifies the point cloud and the resulting canopy model is smoother and contains fewer 'pits' and empty pixels.

Usage

pitfree(
  thresholds = c(0, 2, 5, 10, 15),
  max_edge = c(0, 1),
  subcircle = 0,
  highest = TRUE
)

Arguments

References

Khosravipour, A., Skidmore, A. K., Isenburg, M., Wang, T., & Hussin, Y. A. (2014). Generating pit-free canopy height models from airborne lidar. Photogrammetric Engineering & Remote Sensing, 80(9), 863-872.

See Also

Other digital surface model algorithms: dsm_point2raster, dsm_tin

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
poi = "-drop_z_below 0 -inside 481280 3812940 481330 3812990"
las <- readLAS(LASfile, filter = poi)
col <- height.colors(50)

# Basic triangulation and rasterization of first returns
chm <- rasterize_canopy(las, res = 0.5, dsmtin())
plot(chm, col = col)

# Khosravipour et al. pitfree algorithm
chm <- rasterize_canopy(las, res = 0.5, pitfree(c(0,2,5,10,15), c(0, 1.5)))
plot(chm, col = col)

## Not run: 
# Potentially complex concave subset of point cloud
x = c(481340, 481340, 481280, 481300, 481280, 481340)
y = c(3812940, 3813000, 3813000, 3812960, 3812940, 3812940)
las2 = clip_polygon(las,x,y)
plot(las2)

# Because the TIN interpolation is done within the convex hull of the point cloud
# dummy pixels are interpolated that are correct according to the interpolation
# method used, but meaningless in our CHM
chm <- rasterize_canopy(las2, res = 0.5, pitfree(max_edge = c(0, 1.5)))
plot(chm, col = col)

chm = rasterize_canopy(las2, res = 0.5, pitfree(max_edge = c(3, 1.5)))
plot(chm, col = col)

## End(Not run)

Description

This function is made to be used in rasterize_canopy. It implements an algorithm for digital surface model computation based on a points-to-raster method: for each pixel of the output raster the function attributes the height of the highest point found. The subcircle tweak replaces each point with 8 points around the original one. This allows for virtual 'emulation' of the fact that a lidar point is not a point as such, but more realistically a disc. This tweak densifies the point cloud and the resulting canopy model is smoother and contains fewer 'pits' and empty pixels.

Usage

p2r(subcircle = 0, na.fill = NULL)

Arguments

See Also

Other digital surface model algorithms: dsm_pitfree, dsm_tin

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
col <- height.colors(50)

# Points-to-raster algorithm with a resolution of 1 meter
chm <- rasterize_canopy(las, res = 1, p2r())
plot(chm, col = col)

# Points-to-raster algorithm with a resolution of 0.5 meters replacing each
# point by a 20 cm radius circle of 8 points
chm <- rasterize_canopy(las, res = 0.5, p2r(0.2))
plot(chm, col = col)

## Not run: 
chm <- rasterize_canopy(las, res = 0.5, p2r(0.2, na.fill = tin()))
plot(chm, col = col)

## End(Not run)

Description

This function is made to be used in rasterize_canopy. It implements an algorithm for digital surface model computation using a Delaunay triangulation of first returns with a linear interpolation within each triangle.

Usage

dsmtin(max_edge = 0, highest = TRUE)

Arguments

See Also

Other digital surface model algorithms: dsm_pitfree, dsm_point2raster

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
col <- height.colors(50)

# Basic triangulation and rasterization of first returns
chm <- rasterize_canopy(las, res = 1, dsmtin())
plot(chm, col = col)

## Not run: 
# Potentially complex concave subset of point cloud
x = c(481340, 481340, 481280, 481300, 481280, 481340)
y = c(3812940, 3813000, 3813000, 3812960, 3812940, 3812940)
las2 = clip_polygon(las,x,y)
plot(las2)

# Because the TIN interpolation is done within the convex hull of the point cloud
# dummy pixels are interpolated that are correct according to the interpolation method
# used, but meaningless in our CHM
chm <- rasterize_canopy(las2, res = 0.5, dsmtin())
plot(chm, col = col)

# Use 'max_edge' to trim dummy triangles
chm = rasterize_canopy(las2, res = 0.5, dsmtin(max_edge = 3))
plot(chm, col = col)

## End(Not run)

Description

This function is made to be used in rasterize_terrain or normalize_height. It implements an algorithm for spatial interpolation. Interpolation is done using a k-nearest neighbour (KNN) approach with an inverse-distance weighting (IDW).

Usage

knnidw(k = 10, p = 2, rmax = 50)

Arguments

See Also

Other dtm algorithms: dtm_kriging, dtm_tin

Examples

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)

#plot(las)

dtm = rasterize_terrain(las, algorithm = knnidw(k = 6L, p = 2))

#plot(dtm)
#plot_dtm3d(dtm)

Description

This function is made to be used in rasterize_terrain or normalize_height. It implements an algorithm for spatial interpolation. Spatial interpolation is based on universal kriging using the krige() function from gstat. This method combines the KNN approach with the kriging approach. For each point of interest it kriges the terrain using the k-nearest neighbour ground points. This method is more difficult to manipulate but it is also the most advanced method for interpolating spatial data.

Usage

kriging(model = gstat::vgm(0.59, "Sph", 874), k = 10L)

Arguments

See Also

Other dtm algorithms: dtm_idw, dtm_tin

Examples

## Not run: 
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)

plot(las)

dtm = rasterize_terrain(las, algorithm = kriging())

plot(dtm)
plot_dtm3d(dtm)

## End(Not run)

Description

This function is made to be used in rasterize_terrain or normalize_height. It implements an algorithm for spatial interpolation. Spatial interpolation is based on a Delaunay triangulation, which performs a linear interpolation within each triangle. There are usually a few points outside the convex hull, determined by the ground points at the very edge of the dataset, that cannot be interpolated with a triangulation. Extrapolation can be performed with another algorithm.

Usage

tin(..., extrapolate = knnidw(3, 1, 50))

Arguments

See Also

Other dtm algorithms: dtm_idw, dtm_kriging

Examples

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile, filter = "-inside 273450 5274350 273550 5274450")

#plot(las)

dtm = rasterize_terrain(las, algorithm = tin())

#plot(dtm)
#plot_dtm3d(dtm)

Description

Functions for the LAScatalog processing engine not meant to be called directly by users. They are exported for debugging and to simplify export of internal functions when processing in parallel

Usage

engine_apply(
  .CHUNKS,
  .FUN,
  .PROCESSOPT,
  .OUTPUTOPT,
  .GLOBALS = NULL,
  .AUTOREAD = FALSE,
  .AUTOCROP = FALSE,
  ...
)

engine_chunks(ctg, realignment = FALSE, plot = opt_progress(ctg))

engine_crop(x, bbox)

engine_merge(ctg, any_list, ...)

engine_write(x, path, drivers)

Arguments

See Also

Other LAScatalog processing engine: engine_options

Other LAScatalog processing engine: engine_options

Description

The names of the options and their roles are documented in LAScatalog. The options are used by all the functions that support a LAScatalog as input. Most options are easy to understand and to link to the documentation of LAScatalog but some need more details. See section 'Details'.

Usage

opt_chunk_buffer(ctg)

opt_chunk_buffer(ctg) <- value

opt_chunk_size(ctg)

opt_chunk_size(ctg) <- value

opt_chunk_alignment(ctg)

opt_chunk_alignment(ctg) <- value

opt_restart(ctg) <- value

opt_progress(ctg)

opt_progress(ctg) <- value

opt_stop_early(ctg)

opt_stop_early(ctg) <- value

opt_wall_to_wall(ctg)

opt_wall_to_wall(ctg) <- value

opt_independent_files(ctg)

opt_independent_files(ctg) <- value

opt_output_files(ctg)

opt_output_files(ctg) <- value

opt_laz_compression(ctg)

opt_laz_compression(ctg) <- value

opt_merge(ctg)

opt_merge(ctg) <- value

opt_select(ctg)

opt_select(ctg) <- value

opt_filter(ctg)

opt_filter(ctg) <- value

Arguments

Details

• opt_restart() automatically sets the chunk option named "drop" in such a way that the engine will restart at a given chunk and skip all previous chunks. Useful for restarting after a crash.

• opt_independent_file() automatically sets the chunk size, chunk buffer and wall-to-wall options to process files that are not spatially related to each other, such as plot inventories.

• opt_laz_compression() automatically modifies the drivers to write LAZ files instead of LAS files.

See Also

Other LAScatalog processing engine: engine

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)

plot(ctg, chunk_pattern = TRUE)

opt_chunk_size(ctg) <- 150
plot(ctg, chunk_pattern = TRUE)

opt_chunk_buffer(ctg) <- 10
plot(ctg, chunk_pattern = TRUE)

opt_chunk_alignment(ctg) <- c(270,250)
plot(ctg, chunk_pattern = TRUE)

summary(ctg)

opt_output_files(ctg) <- "/path/to/folder/templated_filename_{XBOTTOM}_{ID}"
summary(ctg)

Description

Operators acting on LAS* objects. However, some have modified behaviors to prevent some irrelevant modifications. Indeed, a LAS* object cannot contain anything, as the content is restricted by the LAS specifications. If a user attempts to use one of these functions inappropriately an informative error will be thrown.

Usage

## S4 method for signature 'LAS'
x$name

## S4 method for signature 'LAS,ANY,missing'
x[[i, j, ...]]

## S4 replacement method for signature 'LAS'
x$name <- value

## S4 replacement method for signature 'LAS,ANY,missing'
x[[i, j]] <- value

## S4 method for signature 'LAS,numeric,ANY'
x[i]

## S4 method for signature 'LAS,logical,ANY'
x[i]

## S4 method for signature 'LAS,sf,ANY'
x[i]

## S4 method for signature 'LAS,sfc,ANY'
x[i]

## S4 method for signature 'LAScatalog'
x$name

## S4 method for signature 'LAScatalog,ANY,missing'
x[[i, j, ...]]

## S4 method for signature 'LAScatalog,ANY,ANY'
x[i, j, ..., drop = FALSE]

## S4 method for signature 'LAScatalog,logical,ANY'
x[i]

## S4 method for signature 'LAScatalog,sf,ANY'
x[i]

## S4 method for signature 'LAScatalog,sfc,ANY'
x[i]

## S4 replacement method for signature 'LAScatalog,ANY,ANY'
x[[i, j]] <- value

## S4 replacement method for signature 'LAScatalog'
x$name <- value

## S4 method for signature 'LASheader'
x$name

## S4 replacement method for signature 'LASheader'
x$name <- value

## S4 method for signature 'LASheader,ANY,missing'
x[[i, j, ...]]

## S4 replacement method for signature 'LASheader,character,missing'
x[[i]] <- value

Arguments

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las = readLAS(LASfile)

las$Intensity
las[["Z"]]
las[["Number of points by return"]]

## Not run: 
las$Z = 2L
las[["Z"]] = 1:10
las$NewCol = 0
las[["NewCol"]] = 0

## End(Not run)

Description

Filter points of interest (POI) from a LAS object where conditions are true.

Usage

filter_poi(las, ...)

filter_first(las)

filter_firstlast(las)

filter_firstofmany(las)

filter_ground(las)

filter_last(las)

filter_nth(las, n)

filter_single(las)

filter_duplicates(las)

## S3 method for class 'LAS'
filter_duplicates(las)

## S3 method for class 'LAScatalog'
filter_duplicates(las)

Arguments

Details

• filter_poi Select points of interest based on custom logical criteria.

• filter_first Select only the first returns.

• filter_firstlast Select only the first and last returns.

• filter_ground Select only the returns classified as ground according to LAS specification.

• filter_last Select only the last returns i.e. the last returns and the single returns.

• filter_nth Select the returns from their position in the return sequence.

• filter_firstofmany Select only the first returns from pulses which returned multiple points.

• filter_single Select only the returns that return only one point.

• filter_duplicates Removes the duplicated points (duplicated by XYZ)

Value

An object of class LAS

Non-supported LAScatalog options

The option select is not supported and not respected because it always preserves the file format and all the attributes. select = "*" is imposed internally.

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
lidar = readLAS(LASfile)

# Select the first returns classified as ground
firstground = filter_poi(lidar, Classification == 2L & ReturnNumber == 1L)

# Multiple arguments are equivalent to &
firstground = filter_poi(lidar, Classification == 2L, ReturnNumber == 1L)

# Multiple criteria
first_or_ground = filter_poi(lidar, Classification == 2L | ReturnNumber == 1L)

Description

This function is made to be used in classify_ground. It implements an algorithm for segmentation of ground points base on a Cloth Simulation Filter. This method is a strict implementation of the CSF algorithm made by Zhang et al. (2016) (see references) that relies on the authors' original source code written and exposed to R via the the RCSF package.

Usage

csf(
  sloop_smooth = FALSE,
  class_threshold = 0.5,
  cloth_resolution = 0.5,
  rigidness = 1L,
  iterations = 500L,
  time_step = 0.65
)

Arguments

References

W. Zhang, J. Qi*, P. Wan, H. Wang, D. Xie, X. Wang, and G. Yan, “An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation,” Remote Sens., vol. 8, no. 6, p. 501, 2016. (http://www.mdpi.com/2072-4292/8/6/501/htm)

See Also

Other ground segmentation algorithms: gnd_mcc, gnd_pmf

Examples

if (require(RCSF, quietly = TRUE))
{
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn")

mycsf <- csf(TRUE, 1, 1, time_step = 1)
las <- classify_ground(las, mycsf)
#plot(las, color = "Classification")
}

Description

This function is made to be used in classify_ground. It implements an algorithm for segmentation of ground points base on a Multiscale Curvature Classification. This method is a strict implementation of the MCC algorithm made by Evans and Hudak. (2007) (see references) that relies on the authors' original source code written and exposed to R via the the RMCC package.

Usage

mcc(s = 1.5, t = 0.3)

Arguments

Details

There are two parameters that the user must define, the scale (s) parameter and the curvature threshold (t). The optimal scale parameter is a function of 1) the scale of the objects (e.g., trees) on the ground, and 2) the sampling interval (post spacing) of the lidar data. Current lidar sensors are capable of collecting high density data (e.g., 8 pulses/m2) that translate to a spatial sampling frequency (post spacing) of 0.35 m/pulse (1/sqrt(8 pulses/m2) = 0.35 m/pulse), which is small relative to the size of mature trees and will oversample larger trees (i.e., nominally multiple returns/tree). Sparser lidar data (e.g., 0.25 pulses/m2) translate to a post spacing of 2.0 m/pulse (1/sqrt(0.25 pulses/m2) = 2.0 m/pulse), which would undersample trees and fail to sample some smaller trees (i.e., nominally <1 return/tree).

Therefore, a bit of trial and error is warranted to determine the best scale and curvature parameters to use. Select a las tile containing a good variety of canopy and terrain conditions, as it's likely the parameters that work best there will be applicable to the rest of your project area tiles. As a starting point: if the scale (post spacing) of the lidar survey is 1.5 m, then try 1.5. Try varying it up or down by 0.5 m increments to see if it produces a more desirable digital terrain model (DTM) interpolated from the classified ground returns in the output file. Use units that match the units of the lidar data. For example, if your lidar data were delivered in units of feet with a post spacing of 3 ft, set the scale parameter to 3, then try varying it up or down by 1 ft increments and compare the resulting interpolated DTMs. If the trees are large, then it's likely that a scale parameter of 1 m (3 ft) will produce a cleaner DTM than a scale parameter of 0.3 m (1 ft), even if the pulse density is 0.3 m (1 ft). As for the curvature threshold, a good starting value to try might be 0.3 (if data are in meters; 1 if data are in feet), and then try varying this up or down by 0.1 m increments (if data are in meters; 0.3 if data are in feet).

References

Evans, Jeffrey S.; Hudak, Andrew T. 2007. A multiscale curvature algorithm for classifying discrete return LiDAR in forested environments. IEEE Transactions on Geoscience and Remote Sensing. 45(4): 1029-1038.

See Also

Other ground segmentation algorithms: gnd_csf, gnd_pmf

Examples

if (require(RMCC, quietly = TRUE))
{
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn", filter = "-inside 273450 5274350 273550 5274450")

las <- classify_ground(las, mcc(1.5,0.3))
#plot(las, color = "Classification")
}

Description

This function is made to be used in classify_ground. It implements an algorithm for segmentation of ground points based on a progressive morphological filter. This method is an implementation of the Zhang et al. (2003) algorithm (see reference). Note that this is not a strict implementation of Zhang et al. This algorithm works at the point cloud level without any rasterization process. The morphological operator is applied on the point cloud, not on a raster. Also, Zhang et al. proposed some formulas (eq. 4, 5 and 7) to compute the sequence of windows sizes and thresholds. Here, these parameters are free and specified by the user. The function util_makeZhangParam enables computation of the parameters according to the original paper.

Usage

pmf(ws, th)

util_makeZhangParam(
  b = 2,
  dh0 = 0.5,
  dhmax = 3,
  s = 1,
  max_ws = 20,
  exp = FALSE
)

Arguments

Details

The progressive morphological filter allows for any sequence of parameters. The 'util_makeZhangParam' function enables computation of the sequences using equations (4), (5) and (7) from Zhang et al. (see reference and details).

In the original paper the windows size sequence is given by eq. 4 or 5:

$w_k = 2kb + 1$

or

$w_k = 2b^k + 1$

In the original paper the threshold sequence is given by eq. 7:

$th_k = s*(w_k - w_{k-1})*c + th_0$

Because the function classify_ground applies the morphological operation at the point cloud level the parameter $c$ is set to 1 and cannot be modified.

References

Zhang, K., Chen, S. C., Whitman, D., Shyu, M. L., Yan, J., & Zhang, C. (2003). A progressive morphological filter for removing nonground measurements from airborne LIDAR data. IEEE Transactions on Geoscience and Remote Sensing, 41(4 PART I), 872–882. http:#doi.org/10.1109/TGRS.2003.810682.

See Also

Other ground segmentation algorithms: gnd_csf, gnd_mcc

Examples

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzrn", filter = "-inside 273450 5274350 273550 5274450")

ws <- seq(3,12, 3)
th <- seq(0.1, 1.5, length.out = length(ws))

las <- classify_ground(las, pmf(ws, th))
#plot(las, color = "Classification")

Description

Full waveform can be difficult to manipulate and visualize in R. This function converts a LAS object with full waveform data into a regular point cloud. Each waveform record becomes a point with XYZ coordinates and an amplitude (units: volts) and an ID that records each original pulse. Notice that this has the effect of drastically inflating the size of the object in memory, which is likely already very large

Usage

interpret_waveform(las)

Arguments

Value

An object of class LAS 1.2 format 0 with one point per records

Full waveform

With most recent versions of the rlas package, full waveform (FWF) can be read and lidR provides some compatible functions. However, the support of FWF is still a work-in-progress in the rlas package. How it is read, interpreted and represented in R may change. Consequently, tools provided by lidR may also change until the support of FWF becomes mature and stable in rlas. See also rlas::read.las.

Remember that FWF represents an insanely huge amount of data. It terms of memory it is like having between 10 to 100 times more points. Consequently, loading FWF data in R should be restricted to relatively small point clouds.

Examples

## Not run: 
LASfile <- system.file("extdata", "fwf.laz", package="rlas")
fwf <- readLAS(LASfile)
las <- interpret_waveform(fwf)
x <- plot(fwf, size = 3, pal = "red")
plot(las, color = "Amplitude", bg = "white", add = x, size = 2)

## End(Not run)

Description

is.empty tests if a LAS object is a point cloud with 0 points.
is.overlapping tests if a LAScatalog has overlapping tiles.
is.indexed tests if the points of a LAScatalog are indexed with .lax files.
is.algorithm tests if an object is an algorithm of the lidR package.
is.parallelised tests if an algorithm of the lidR package is natively parallelised with OpenMP. Returns TRUE if the algorithm is at least partially parallelised i.e. if some portion of the code is computed in parallel.

Usage

## S4 method for signature 'LAS'
is.empty(x)

is.overlapping(catalog)

is.indexed(catalog)

is.algorithm(x)

is.parallelised(algorithm)

Arguments

Value

TRUE or FALSE

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las = readLAS(LASfile)
is.empty(las)

las = new("LAS")
is.empty(las)

f <- lmf(2)
is.parallelised(f)

g <- pitfree()
is.parallelised(g)

ctg <- readLAScatalog(LASfile)
is.indexed(ctg)

Description

This function is made to be used in locate_trees. It implements an algorithm for tree detection based on a local maximum filter. The windows size can be fixed or variable and its shape can be square or circular. The internal algorithm works either with a raster or a point cloud. It is deeply inspired by Popescu & Wynne (2004) (see references).

Usage

lmf(ws, hmin = 2, shape = c("circular", "square"), ws_args = "Z")

Arguments

References

Popescu, Sorin & Wynne, Randolph. (2004). Seeing the Trees in the Forest: Using Lidar and Multispectral Data Fusion with Local Filtering and Variable Window Size for Estimating Tree Height. Photogrammetric Engineering and Remote Sensing. 70. 589-604. 10.14358/PERS.70.5.589.

See Also

Other individual tree detection algorithms: itd_manual

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzi", filter = "-inside 481250 3812980 481300 3813050")

# =================
# point-cloud-based
# =================

# 5x5 m fixed window size
ttops <- locate_trees(las, lmf(5))

#plot(las) |> add_treetops3d(ttops)

# variable windows size
f <- function(x) { x * 0.07 + 3}
ttops <- locate_trees(las, lmf(f))

#plot(las) |> add_treetops3d(ttops)

# Very custom variable windows size
f <- function(x, y, z) { x * 0.07 + y * 0.01 + z}
ws_args <- list(x = "Z", y = "Intensity", z = 3)
ttops <- locate_trees(las, lmf(f, ws_args = ws_args))

# ============
# raster-based
# ============

chm <- rasterize_canopy(las, res = 1, p2r(0.15), pkg = "terra")
ttops <- locate_trees(chm, lmf(5))

plot(chm, col = height.colors(30))
plot(sf::st_geometry(ttops), add = TRUE, col = "black", cex = 0.5, pch = 3)

# variable window size
f <- function(x) { x * 0.07 + 3 }
ttops <- locate_trees(chm, lmf(f))

plot(chm, col = height.colors(30))
plot(sf::st_geometry(ttops), add = TRUE, col = "black", cex = 0.5, pch = 3)

Description

This function is made to be used in locate_trees. It implements an algorithm for manual tree detection. Users can pinpoint the tree top positions manually and interactively using the mouse. This is only suitable for small-sized plots. First the point cloud is displayed, then the user is invited to select a rectangular region of interest in the scene using the mouse button. Within the selected region the highest point will be flagged as 'tree top' in the scene. Once all the trees are labelled the user can exit the tool by selecting an empty region. Points can also be unflagged. The goal of this tool is mainly for minor correction of automatically-detected tree outputs.
This algorithm does not preserve tree IDs from detected and renumber all trees. It also looses all attributes

Usage

manual(detected = NULL, radius = 0.5, color = "red", button = "middle", ...)

Arguments

See Also

Other individual tree detection algorithms: itd_lmf

Examples

## Not run: 
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las = readLAS(LASfile)

# Full manual tree detection
ttops = locate_trees(las, manual())

# Automatic detection with manual correction
ttops = locate_trees(las, lmf(5))
ttops = locate_trees(las, manual(ttops))

## End(Not run)

Description

This function is made to be used in segment_trees. It implements an algorithm for tree segmentation based on Dalponte and Coomes (2016) algorithm (see reference). This is a seeds + growing region algorithm. This algorithm exists in the package itcSegment. This version has been written from the paper in C++. Consequently it is hundreds to millions times faster than the original version. Note that this algorithm strictly performs a segmentation, while the original method as implemented in itcSegment and described in the manuscript also performs pre- and post-processing tasks. Here these tasks are expected to be done by the user in separate functions.

Usage

dalponte2016(
  chm,
  treetops,
  th_tree = 2,
  th_seed = 0.45,
  th_cr = 0.55,
  max_cr = 10,
  ID = "treeID"
)

Arguments

Details

Because this algorithm works on a CHM only there is no actual need for a point cloud. Sometimes the user does not even have the point cloud that generated the CHM. lidR is a point cloud-oriented library, which is why this algorithm must be used in segment_trees to merge the result with the point cloud. However the user can use this as a stand-alone function like this:

 chm <- raster("chm.tif")
 ttops <- locate_trees(chm, lmf(3))
 crowns <- dalponte2016(chm, ttops)()

References

Dalponte, M. and Coomes, D. A. (2016), Tree-centric mapping of forest carbon density from airborne laser scanning and hyperspectral data. Methods Ecol Evol, 7: 1236–1245. doi:10.1111/2041-210X.12575.

See Also

Other individual tree segmentation algorithms: its_li2012, its_silva2016, its_watershed

Other raster based tree segmentation algorithms: its_silva2016, its_watershed

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
poi <- "-drop_z_below 0 -inside 481280 3812940 481320 3812980"
las <- readLAS(LASfile, select = "xyz", filter = poi)
col <- pastel.colors(200)

chm <- rasterize_canopy(las, 0.5, p2r(0.3))
ker <- matrix(1,3,3)
chm <- terra::focal(chm, w = ker, fun = mean, na.rm = TRUE)

ttops <- locate_trees(chm, lmf(4, 2))
las   <- segment_trees(las, dalponte2016(chm, ttops))
#plot(las, color = "treeID", colorPalette = col)

Description

This functions is made to be used in segment_trees. It implements an algorithm for tree segmentation based on Li et al. (2012) (see reference). This method is a growing region method working at the point cloud level. It is an implementation by lidR authors, from the original paper, as close as possible from the original description. However we added a parameter hmin to prevent over-segmentation for objects that are too low. This algorithm is known to be slow because it has an algorithmic complexity worst that O(n^2).

Usage

li2012(dt1 = 1.5, dt2 = 2, R = 2, Zu = 15, hmin = 2, speed_up = 10)

Arguments

References

Li, W., Guo, Q., Jakubowski, M. K., & Kelly, M. (2012). A new method for segmenting individual trees from the lidar point cloud. Photogrammetric Engineering & Remote Sensing, 78(1), 75-84.

See Also

Other individual tree segmentation algorithms: its_dalponte2016, its_silva2016, its_watershed

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
poi <- "-drop_z_below 0 -inside 481280 3812940 481320 3812980"
las <- readLAS(LASfile, select = "xyz", filter = poi)
col <- pastel.colors(200)

las <- segment_trees(las, li2012(dt1 = 1.4))
#plot(las, color = "treeID", colorPalette = col)

Description

This functions is made to be used in segment_trees. It implements an algorithm for tree segmentation based on Silva et al. (2016) (see reference). This is a simple method based on seed + voronoi tesselation (equivalent to nearest neighbour). This algorithm is implemented in the package rLiDAR. This version is not the version from rLiDAR. It is code written from the original article by the lidR authors and is considerably (between 250 and 1000 times) faster.

Usage

silva2016(chm, treetops, max_cr_factor = 0.6, exclusion = 0.3, ID = "treeID")

Arguments

Details

Because this algorithm works on a CHM only there is no actual need for a point cloud. Sometimes the user does not even have the point cloud that generated the CHM. lidR is a point cloud-oriented library, which is why this algorithm must be used in segment_trees to merge the result into the point cloud. However, the user can use this as a stand-alone function like this:

 chm <- raster("chm.tif")
 ttops <- locate_trees(chm, lmf(3))
 crowns <- silva2016(chm, ttops)()

References

Silva, C. A., Hudak, A. T., Vierling, L. A., Loudermilk, E. L., O’Brien, J. J., Hiers, J. K., Khosravipour, A. (2016). Imputation of Individual Longleaf Pine (Pinus palustris Mill.) Tree Attributes from Field and LiDAR Data. Canadian Journal of Remote Sensing, 42(5), 554–573. https://doi.org/10.1080/07038992.2016.1196582.

See Also

Other individual tree segmentation algorithms: its_dalponte2016, its_li2012, its_watershed

Other raster based tree segmentation algorithms: its_dalponte2016, its_watershed

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
poi <- "-drop_z_below 0 -inside 481280 3812940 481320 3812980"
las <- readLAS(LASfile, select = "xyz", filter = poi)
col <- pastel.colors(200)

chm <- rasterize_canopy(las, res = 0.5, p2r(0.3))
ker <- matrix(1,3,3)
chm <- terra::focal(chm, w = ker, fun = mean, na.rm = TRUE)

ttops <- locate_trees(chm, lmf(4, 2))
las   <- segment_trees(las, silva2016(chm, ttops))
#plot(las, color = "treeID", colorPalette = col)

Description

This function is made to be used in segment_trees. It implements an algorithm for tree segmentation based on a watershed. It is based on the bioconductor package EBIimage. You need to install this package to run this method (see its github page). Internally, the function EBImage::watershed is called.

Usage

watershed(chm, th_tree = 2, tol = 1, ext = 1)

Arguments

Details

Because this algorithm works on a CHM only there is no actual need for a point cloud. Sometimes the user does not even have the point cloud that generated the CHM. lidR is a point cloud-oriented library, which is why this algorithm must be used in segment_trees to merge the result into the point cloud. However, the user can use this as a stand-alone function like this:

 chm <- raster("chm.tif")
 crowns <- watershed(chm)()

See Also

Other individual tree segmentation algorithms: its_dalponte2016, its_li2012, its_silva2016

Other raster based tree segmentation algorithms: its_dalponte2016, its_silva2016

Examples

## Not run: 
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
poi <- "-drop_z_below 0 -inside 481280 3812940 481320 3812980"
las <- readLAS(LASfile, select = "xyz", filter = poi)
col <- pastel.colors(250)

# Using raster because focal does not exist in stars
chm <- rasterize_canopy(las, res = 0.5, p2r(0.3), pkg = "raster")
ker <- matrix(1,3,3)
chm <- raster::focal(chm, w = ker, fun = mean, na.rm = TRUE)
las <- segment_trees(las, watershed(chm))

plot(las, color = "treeID", colorPalette = col)

## End(Not run)

Description

Performs a deep inspection of a LAS or LAScatalog object and prints a report.

For a LAS object it checks:

• if the point cloud is valid according to las specification

• if the header is valid according to las specification

• if the point cloud is in accordance with the header

• if the point cloud has duplicated points and degenerated ground points

• if gpstime and pulses are consistent

• it the coordinate reference sytem is correctly recorded

• if some pre-processing, such as normalization or ground filtering, is already done.

• and much more

For a LAScatalog object it checks:

• if the headers are consistent across files

• if the files are overlapping

• if some pre-processing, such as normalization, is already done.

For the pre-processing tests the function only makes an estimation and may not be correct.

Usage

las_check(las, print = TRUE, ...)

Arguments

Value

A list with three elements named message, warnings and errors. This list is returned invisibly if print = TRUE. If deep = TRUE a nested list is returned with one element per file.

See Also

Other las utilities: las_utilities

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)
las_check(las)

Description

Package rlas 1.6.0 supports compact representation of non populated attributes. For example UserData is usually populated with zeros (not populated). Yet it takes 32 bits per point to store each 0. With rlas 1.6.0 it can now use 644 bits no matter the number of points loaded if it is not populated or populated with a unique value.

Usage

las_is_compressed(las)

las_size(las)

Arguments

Details

las_is_compressed test each attributes and returns a named vector with TRUE if the attribute is compressed FALSE otherwise.

las_size returns the true size of a LAS object by considering the compression. object.size from base R does not account for ALTREP and consequently cannot measure properly the size of a LAS object

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
las_is_compressed(las)

format(object.size(las), units = "MB")
format(las_size(las), units = "MB")

Description

Tools to manipulate LAS objects maintaining compliance with ASPRS specification

Usage

las_rescale(las, xscale, yscale, zscale)

las_reoffset(las, xoffset, yoffset, zoffset)

las_quantize(las, by_reference = TRUE)

las_update(las)

quantize(x, scale, offset, by_reference = TRUE, ...)

is.quantized(x, scale, offset, ...)

count_not_quantized(x, scale, offset)

storable_coordinate_range(scale, offset)

header(las)

payload(las)

phb(las)

vlr(las)

evlr(las)

Arguments

Details

In the specification of the LAS format the coordinates are expected to be given with a certain precision e.g. 0.01 for a millimeter precision (or millifeet), meaning that a file records e.g. 123.46 not 123.45678. Also, coordinates are stored as integers. This is made possible with a scale and offset factor. For example, 123.46 with an offset of 100 and a scale factor of 0.01 is actually stored as (123.46 - 100)/0.01 = 2346. Storing 123.45678 with a scale factor of 0.01 and an offset of 100 is invalid because it does not convert to an integer: (123.45678-100)/0.01 = 2345.678. Having an invalid LAS object may be critical in some lidR applications. When writing into a LAS file, users will loose the extra precision without warning and some algorithms in lidR use the integer conversion to make integer-based computation and thus speed-up some algorithms and use less memory. Creation of an invalid LAS object may cause problems and incorrect outputs.

See Also

Other las utilities: las_check()

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las = readLAS(LASfile)

# Manual modification of the coordinates (e.g. rotation, re-alignment, ...)
las@data$X <- las@data$X + 2/3
las@data$Y <- las@data$Y - 5/3

# The point cloud is no longer valid
las_check(las)

# It is important to fix that
las_quantize(las)

# Now the file is almost valid
las_check(las)

# Update the object to set up-to-date header data
las <- las_update(las)
las_check(las)

# In practice the above code is not useful for regular users because the operators
# $<- already perform such operations on-the-fly. Thus the following
# syntax must be preferred and returns valid objects. Previous tools
# were only intended to be used in very specific cases.
las$X <- las$X + 2/3
las$Y <- las$Y - 5/3

# Rescale and reoffset recompute the coordinates with
# new scales and offsets according to LAS specification
las <- las_rescale(las, xscale = 0.01, yscale = 0.01)
las <- las_reoffset(las, xoffset = 300000, yoffset = 5248000)

Description

Class LAS is the representation of a las/laz file according to the LAS file format specifications.

Usage

LAS(data, header = list(), crs = sf::NA_crs_, check = TRUE, index = NULL, ...)

Arguments

Details

A LAS object contains a data.table with the data read from a las/laz file and a LASheader (see the ASPRS documentation for the LAS file format for more information). Because las files are standardized the table of attributes read from the las/laz file is also standardized. Columns are named:

• X, Y, Z (numeric)

• gpstime (numeric)

• Intensity (integer)

• ReturnNumber, NumberOfReturns (integer)

• ScanDirectionFlag (integer)

• EdgeOfFlightline (integer)

• Classification (integer)

• Synthetic_flag,Keypoint_flag, Withheld_flag (logical)

• ScanAngleRank/ScanAngle (integer/numeric)

• UserData (integer)

• PointSourceID (integer)

• R,G,B, NIR (integer)

Value

An object of class LAS

Functions

• LAS(): creates objects of class LAS. The original data is updated by reference to quantize the coordinates according to the scale factor of the header if no header is provided. In this case the scale factor is set to 0.001

Slots

crs

Object of class crs from sf.

data

Object of class data.table. Point cloud data according to the LAS file format

header

Object of class LASheader. LAS file header according to the LAS file format

index

list. See spatial indexing.

See Also

readLAS

Examples

# Read a las/laz file
LASfile <- system.file("extdata", "example.laz", package="rlas")
las <- readLAS(LASfile)
las

# Creation of a LAS object out of external data
data <- data.frame(X = runif(100, 0, 100),
                   Y = runif(100, 0, 100),
                   Z = runif(100, 0, 20))

# 'data' has many decimal digits
data

# Create a default header and quantize *by reference*
# the coordinates to fit with offset and scale factors
cloud <- LAS(data)

# 'data' has been updated and coordinates were quantized
data
cloud

# Be careful when providing a header the function assumes that
# it corresponds to the data and won't quantize the coordinates
data <- data.frame(X = runif(100, 0, 100),
                   Y = runif(100, 0, 100),
                   Z = runif(100, 0, 20))
header <- header(las)

# This works but triggers warnings and creates an invalid LAS object
cloud <- LAS(data, header)

las_check(cloud)

Description

A LAScatalog object is a representation of a collection of las/laz files. A LAScatalog is a way to manage and batch process a lidar coverage. It allows the user to process a large area, or to selectively clip data from a large area without loading all the data into computer memory. A LAScatalog can be built with the function readLAScatalog.

Details

A LAScatalog contains an sf object to store the geometry and metadata. It is extended with slots that contain processing options. In lidR, each function that supports a LAScatalog as input will respect these processing options. Internally, processing a catalog is almost always the same and relies on a few steps:

1. Define chunks. A chunk is an arbitrarily-defined region of interest (ROI) of the collection. Altogether, the chunks are a wall-to-wall set of ROIs that encompass the whole dataset.

2. Loop over each chunk (in parallel or not).

3. For each chunk, load the points inside the ROI into R, run some R functions, return the expected output.

4. Merge the outputs of the different chunks once they are all processed to build a continuous (wall-to-wall) output.

So basically, a LAScatalog is an object that allows for batch processing but with the specificity that lidR does not loop through LAS or LAZ files, but loops seamlessly through chunks that do not necessarily match with the file pattern. This way lidR can sequentially process tiny ROIs even if each file may be individually too big to fit in memory. This is also why point cloud indexation with lax files may significantly speed-up the processing.

It is important to note that catalogs with files that overlap each other are not natively supported by lidR. When encountering such datasets the user should always filter any overlaps if possible. This is possible if the overlapping points are flagged, for example in the 'withheld' attribute. Otherwise lidR will not be able to process the dataset correctly.

Slots

data

sf. An sf data.frame with the bounding box of each file as well as all the information read from the header of each LAS/LAZ file.

processing_options

list. A list that contains some settings describing how the collection will be processed (see dedicated section).

chunk_options

list. A list that contains some settings describing how the collection will be sub-divided into chunks to be processed (see dedicated section).

output_options

list. A list that contains some settings describing how the collection will return the outputs (see dedicated section).

input_options

list. A list of parameters to pass to readLAS (see dedicated section).

index

list. See spatial indexing.

Processing options

The slot ⁠@processing_options⁠ contains a list of options that determine how chunks (the sub-areas that are sequentially processed) are processed.

• progress: boolean. Display a progress bar and a chart of progress. Default is TRUE. Progress estimation can be enhanced by installing the package progress. See opt_progress.

• stop_early: boolean. Stop the processing if an error occurs in a chunk. If FALSE the process can run until the end, removing chunks that failed. Default is TRUE and the user should have no reason to change this. See opt_stop_early.

• wall.to.wall logical. The catalog processing engine always guarantees to return a continuous output without edge effects, assuming that the catalog is a wall-to-wall catalog. To do so, some options are checked internally to guard against bad settings, such as buffer = 0 for an algorithm that requires a buffer. In rare cases it might be useful to disable these controls. If wall.to.wall = FALSE controls are disabled and wall-to-wall outputs cannot be guaranteed. See opt_wall_to_wall

Chunk options

The slot ⁠@chunk_options⁠ contains a list of options that determine how chunks (the sub-areas that are sequentially processed) are made.

• chunk_size: numeric. The size of the chunks that will be sequentially processed. A small size allows small amounts of data to be loaded at once, saving computer memory. With big chunks the computation is usually faster but uses much more memory. If chunk_size = 0 the chunk pattern is build using the file pattern. The chunks are expecting to be wall-to-wall coverage, which means that chunk_size = 0 has meaning only if the files are not overlapping. Default is 0 i.e. by default the processing engine respects the existing tiling pattern. See opt_chunk_size.

• buffer: numeric. Each chunk can be read with an extra buffer around it to ensure there are no edge effects between two independent chunks and that the output is continuous. This is mandatory for some algorithms. Default is 30. See opt_chunk_buffer.

• alignment: numeric. A vector of size 2 (x and y coordinates, respectively) to align the chunk pattern. By default the alignment is made along (0,0), meaning that the edge of the first chunk will belong on x = 0 and y = 0 and all the the other chunks will be multiples of the chunk size. Not relevant if chunk_size = 0. See opt_chunk_alignment.

• drop: integers. A vector of integers that specify the IDs of the chunks that should not be created. This is designed to enable users to restart a computation that failed without reprocessing everything. See opt_restart<-. Technically, this option may be used for partial processing of a collection, but it generally should not be. Partial processing is already a feature of the engine. See this vignette

Output options

The slot ⁠@output_options⁠ contains a list of options that determine how chunks (the sub-areas that are sequentially processed) are written. By "written" we mean written to files or written in R memory.

• output_files: string. If output_files = "" outputs are returned in R. Otherwise, if output_files is a string the outputs will be written to files. This is useful if the output is too big to be returned in R. A path to a filename template without a file extension (the engine guesses it for you) is expected. When several files are going to be written a single string is provided with a template that is automatically filled. For example, the following file names are possible:

  "/home/user/als/normalized/file_{ID}_segmented"
  "C:/user/document/als/zone52_{XLEFT}_{YBOTTOM}_confidential"
  "C:/user/document/als/{ORIGINALFILNAME}_normalized"
  

  This option will generate as many filenames as needed with custom names for each file. The allowed templates are {XLEFT}, {XRIGHT}, {YBOTTOM}, {YTOP}, {ID}, {XCENTER}, {YCENTER}, {ORIGINALFILENAME}. See opt_output_files.

• drivers: list. This contains all the drivers required to seamlessly write ⁠Raster*⁠, SpatRaster, stars, ⁠Spatial*⁠, sf, and LAS objects. It is recommended that only advanced users change this option. A dedicated page describes the drivers in lidR-LAScatalog-drivers.

• merge: boolean. Multiple objects are merged into a single object at the end of the processing. See opt_merge.

Input options

The slot ⁠@input_options⁠ contains a list of options that are passed to the function readLAS. Indeed, the readLAS function is not called directly by the user but by the internal processing engine. Users can propagate these options through the LAScatalog settings.

• select: string. The option select. Usually this option is not respected because each function knows which data must be loaded or not. This is documented in each function. See opt_select.

• filter: string. The option filter. See opt_filter.

Examples

## Not run: 
# Build a catalog
ctg <- readLAScatalog("filder/to/las/files/", filter =  "-keep_first")


# Summary gives a summary of how the catalog will be processed
summary(ctg)

# We can seamlessly use lidR functions
hmean <- pixel_metrics(ctg, ~~mean(Z), 20)
ttops <- tree_detection(ctg, lmf(5))

# For low memory config it is probably advisable not to load entire files
# and process chunks instead
opt_chunk_size(ctg) <- 500

# Sometimes the output is likely to be very large
# e.g. large coverage and small resolution
dtm <- rasterize_terrain(ctg, 1, tin())

# In that case it is advisable to write the output(s) to files
opt_output_files(ctg) <- "path/to/folder/DTM_chunk_{XLEFT}_{YBOTTOM}"

# Raster will be written to disk. The list of written files is returned
# or, in this specific case, a virtual raster mosaic.
dtm <- rasterize_terrain(ctg, 1, tin())

# When chunks are files the original names of the las files can be preserved
opt_chunk_size(ctg) <- 0
opt_output_files(ctg) <- "path/to/folder/DTM_{ORIGINALFILENAME}"
dtm <- rasterize_terrain(ctg, 1, tin())

# For some functions, files MUST be written to disk. Indeed, it is certain that R cannot
# handle the entire output.
opt_chunk_size(ctg) <- 0
opt_output_files(ctg) <- "path/to/folder/{ORIGINALFILENAME}_norm"
opt_laz_compression(ctg) <- TRUE
new_ctg <- normalize_height(ctg, tin())

# The user has access to the catalog engine through the functions catalog_apply
# and catalog_map
output <- catalog_apply(ctg, FUN, ...)

## End(Not run)

Description

Creates a LASheader object either from a raw list containing all the elements named according to the rlas package or creates a header from a data.frame or data.table containing a point cloud. In the latter case it will generate a header according to the data using rlas::header_create(). It will guess the LAS file format, the point data format, and initialize the scale factors and offsets, but these may not suit a user's needs. Users are advised to manually modify the results to fit their specific needs.

Usage

LASheader(data = list())

Arguments

Value

An object of class LASheader

Examples

data = data.frame(X = c(339002.889, 339002.983, 339002.918),
                  Y = c(5248000.515, 5248000.478, 5248000.318),
                  Z = c(975.589, 974.778, 974.471),
                  gpstime = c(269347.28141, 269347.28142, 269347.28143),
                  Intensity = c(82L, 54L, 27L),
                  ReturnNumber = c(1L, 1L, 2L),
                  NumberOfReturns = c(1L, 1L, 2L),
                  ScanDirectionFlag = c(1L, 1L, 1L),
                  EdgeOfFlightline = c(1L, 0L, 0L),
                  Classification = c(1L, 1L, 1L),
                  ScanAngleRank = c(-21L, -21L, -21L),
                  UserData = c(32L, 32L, 32L),
                  PointSourceID = c(17L, 17L, 17L))

header = LASheader(data)
header

# Record an EPSG code
epsg(header) <- 32618
header

las <- LAS(data, header)
las

# The function inferred a LAS 1.2 format 1 which is correct
# Upgrade to LAS 1.4 for the example
header@VLR <- list() # Erase VLR previously written
header@PHB[["Global Encoding"]][["WKT"]] <- TRUE
header@PHB[["Version Minor"]] <- 4L
header@PHB[["Header Size"]] <- 375L
header@PHB[["Offset to point data"]] <- 375L
wkt(header) <- sf::st_crs("EPSG:32618")$wkt
header
las1.4 <- LAS(data, header)
las1.4

Description

An S4 class to represent the header of .las or .laz files according to the LAS file format specifications. A LASheader object contains a list in the slot ⁠@PHB⁠ with the data read from the Public Header Block, a list in the slot ⁠@VLR⁠ with the data read from the Variable Length Records and a list in the slot EVLR with the data read from the Extended Variable Lenght Records.

Slots

PHB

list. Represents the Public Header Block

VLR

list. Represents the Variable Length Records

EVLR

list. Represents the Extended Variable Length Records

Description

This document explains how objects are written on disk when processing a LAScatalog. As mentioned in LAScatalog-class, users can set a templated filename to store the outputs on disk instead of in R memory. By defaut LAS objects are stored in .las files with writeLAS, raster objects are stored in .tif files using native write function for raster, stars or terra, ⁠Spatial*⁠ objects are stored in .shp files with after coercion to sf with st_write, sf object are stored to .gpkg files with st_write.data.frame objects are stored in .csv files with fwrite, and other objects are not supported. However, users can modify all these default settings and even add new drivers. This manual page explain how. One may also refer to some unofficial documentation here or here.

Generic form of a driver

A driver is stored in the slot ⁠@output_options⁠ of a LAScatalog. It is a list that contains:

write

A function that receives an object and a path, and writes the object into a file using the path. The function can also have extra options.

extension

A string that gives the file extension.

object

A string that gives the name of the argument used to pass the object to write in the function used to write the object.

path

A string that gives the name of the argument used to pass the path of the file to write in the function used to write the object.

param

A labelled list of extra parameters for the function used to write the object

For example, the driver to write a ⁠Raster*⁠ is

list(
 write = raster::writeRaster,
 extension = ".tif",
 object = "x",
 path = "filename",
 param = list(format = "GTiff"))

And the driver to write a LAS is

list(
 write = lidR::writeLAS,
 extension = ".las",
 object = "las",
 path = "file",
 param = list())

Modify a driver (1/2)

Users can modify the drivers to write different file types than the default. For example, to write in shapefile instead of a GeoPackage, one must change the sf driver:

ctg@output_options$drivers$sf$extension <- ".shp"

To write a ⁠Raster*⁠ in .grd files instead of .tif files one must change the Raster driver:

ctg@output_options$drivers$Raster$extension <- ".grd"
ctg@output_options$drivers$Raster$param$format <- "raster"

To write in .laz files instead of .las files one must change the LAS driver:

ctg@output_options$drivers$LAS$extension <- ".laz"

Add a new driver

The drivers allow LAS, Spatial,⁠sf,⁠ Raster, stars, SpatRaster and data.frame objects to be written. When using the engine (catalog_apply) to build new tools, users may need to be able to write other objects such as a list. To do that users need to add a list element into ⁠@output_options⁠:

ctg@output_options$drivers$list = list(
 write = base::saveRDS,
 object = "object",
 path = "file",
 extension = ".rds",
 param = list(compress = TRUE))

The LAScatalog now has a new driver capable of writing a list.

Modify a driver (2/2)

It is also possible to completely overwrite an existing driver. By default sf objects are written into GeoPackage with st_write. st_write can also wite in GeoJSON and even in SQLlite database objects. But it cannot add data into an existing SQLlite database. Let's create our own driver for a sf. First we need a function able to write and append a sf into a SQLlite database from the object and the path.

dbWrite_sf = function(x, path, name)
{
 x <- sf::st_drop_geometry(x)
 con <- RSQLite::dbConnect(RSQLite::SQLite(), path)
 RSQLite::dbWriteTable(con, name, x, append = TRUE)
 RSQLite::dbDisconnect(con)
}

Then we create the driver. User-defined drivers supersede default drivers:

ctg@output_options$drivers$sf = list(
 write = dbWrite_sf,
 extension = ".sqlite",
 object = "x",
 path = "path",
 param = list(name = "layername"))

Then to be sure that we do not write several .sqlite files, we don't use templated filename.

opt_output_files(ctg) <- paste0(tempdir(), "/mysqlitefile")

And all the sf will be appended in a single database. To preserve the geometry one can

Description

This document explains how to process point clouds taking advantage of parallel processing in the lidR package. The lidR package has two levels of parallelism, which is why it is difficult to understand how it works. This page aims to provide users with a clear overview of how to take advantage of multicore processing even if they are not comfortable with the parallelism concept.

Algorithm-based parallelism

When processing a point cloud we are applying an algorithm on data. This algorithm may or may not be natively parallel. In lidR some algorithms are fully computed in parallel, but some are not because they are not parallelizable, while some are only partially parallelized. It means that some portions of the code are computed in parallel and some are not. When an algorithm is natively parallel in lidR it is always a C++ based parallelization with OpenMP. The advantage is that the computation is faster without any consequence for memory usage because the memory is shared between the processors. In short, algorithm-based parallelism provides a significant gain without any cost for your R session and your system (but obviously there is a greater workload for the processors). By default lidR uses half of your cores but you can control this with set_lidr_threads. For example, the lmf algorithm is natively parallel. The following code is computed in parallel:

las  <- readLAS("file.las")
tops <- locate_trees(las, lmf(2))

However, as stated above, not all algorithms are parallelized or even parallelizable. For example, li2012 is not parallelized. The following code is computed in serial:

las <- readLAS("file.las")
dtm <- segment_trees(las, li2012())

To know which algorithms are parallelized users can refer to the documentation or use the function is.parallelised.

is.parallelised(lmf(2))   #> TRUE
is.parallelised(li2012()) #> FALSE

Chunk-based parallelism

When processing a LAScatalog, the internal engine splits the dataset into chunks and each chunk is read and processed sequentially in a loop. This loop can be parallelized with the future package. By default the chunks are processed sequentially, but they can be processed in parallel by registering an evaluation strategy. For example, the following code is evaluated sequentially:

ctg <- readLAScatalog("folder/")
out <- pixel_metrics(ctg, ~mean(Z))

But this one is evaluated in parallel with two cores:

library(future)
plan(multisession, workers = 2L)
ctg <- readLAScatalog("folder/")
out <- pixel_metrics(ctg, ~mean(Z))

With chunk-based parallelism any algorithm can be parallelized by processing several subsets of a dataset. However, there is a strong cost associated with this type of parallelism. When processing several chunks at a time, the computer needs to load the corresponding point clouds. Assuming the user processes one square kilometer chunks in parallel with 4 cores, then 4 chunks are loaded in the computer memory. This may be too much and the speed-up is not guaranteed since there is some overhead involved in reading several files at a time. Once this point is understood, chunk-based parallelism is very powerful since all the algorithms can be parallelized whether or not they are natively parallel. It also allows to parallelize the computation on several machines on the network or to work on a HPC.

Nested parallelism - part 1

Previous sections stated that some algorithms are natively parallel, such as lmf, and some are not, such as li2012. Anyway, users can split the dataset into chunks to process them simultaneously with the LAScatalog processing engine. Let's assume that the user's computer has four cores, what happens in this case:

library(future)
plan(multisession, workers = 4L)
set_lidr_threads(4L)
ctg <- readLAScatalog("folder/")
out <- locate_trees(ctg, lmf(2))

Here the catalog will be split into chunks that will be processed in parallel. And each computation itself implies a parallelized task. This is a nested parallelism task and it is dangerous! Hopefully the lidR package handles such cases and chooses by default to give precedence to chunk-based parallelism. In this case chunks will be processed in parallel and the points will be processed serially by disabling OpenMP.

Nested parallelism - part 2

We explained rules of precedence. But actually the user can tune the engine more accurately. Let's define the following function:

myfun = function(las, ws, ...)
{
  las  <- normalize_height(las, tin())
  tops <- locate_tree(las, lmf(ws))
  return(tops)
}

out <- catalog_map(ctg, myfun, ws = 5)

This function used two algorithms, one is partially parallelized (tin) and one is fully parallelized lmf. The user can manually use both OpenMP and future. By default the engine will give precedence to chunk-based parallelism because it works in all cases but the user can impose something else. In the following 2 workers are attributed to future and 2 workers are attributed to OpenMP.

plan(multisession, workers = 2L)
set_lidr_threads(2L)
catalog_map(ctg, myfun, ws = 5)

The rule is simple. If the number of workers needed is greater than the number of available workers then OpenMP is disabled. Let suppose we have a 4 cores:

# 2 chunks 2 threads: OK
plan(multisession, workers = 2L)
set_lidr_threads(2L)

# 4 chunks 1 threads: OK
plan(multisession, workers = 4L)
set_lidr_threads(1L)

# 1 chunks 4 threads: OK
plan(sequential)
set_lidr_threads(4L)

# 3 chunks 2 threads: NOT OK
# Needs 6 workers, OpenMP threads are set to 1 i.e. sequential processing
plan(multisession, workers = 3L)
set_lidr_threads(2L)

Complex computing architectures

For more complex processing architectures such as multiple computers controlled remotely or HPC a finer tuning might be necessary. Using

options(lidR.check.nested.parallelism = FALSE)

lidR will no longer check for nested parallelism and will never automatically disable OpenMP.

Description

This document explains how to process point-clouds taking advantage of different spatial indexes available in the lidR package. lidR can use several types of spatial indexes to apply algorithms (that need a spatial indexing) as fast as possible. The choice of the spatial index depends on the type of point-cloud that is processed and the algorithm that is performed. lidR can use a grid partition, a voxel partition, a quadtree or an octree. See details.

Usage

sensor(las, h = FALSE)

sensor(las) <- value

index(las, h = FALSE)

index(las) <- value

Arguments

Details

From lidR (>= 3.1.0), a LAS object records the sensor used to sample the point-cloud (ALS, TLS, UAV, DAP) as well as the spatial index that must be used for processing the point cloud. This can be set manually by the user but the simplest is to use one of the read*LAS() functions. By default a point-cloud is associated to a sensor and the best spatial index is chosen on-the-fly depending on the algorithm applied. It is possible to force the use of a specific spatial index.

Information relative to the spatial indexing is stored in slot ⁠@index⁠ that contains a list with two elements:

• sensor: an integer that records the sensor type

• index: an integer that records the spatial index to be used

By default the spatial index code is 0 ("automatic") meaning that each function is free to choose a different spatial index depending on the recorded sensor. If the code is not 0 then each function will be forced to used the spatial index that is imposed. This, obviously, applies only to functions that use spatial indexing.

LAScatalog objects also record such information that is automatically propagated to the LAS objects when processing.

Note: before version 3.1.0, point-clouds were all considered as ALS because lidR was originally designed for ALS. Consequently, for legacy and backwards-compatibility reasons, readLAS() and readALSLAS() are actually equivalent. readLAS() tags the point cloud with "unknown" sensor while readALSLAS() tags it with 'ALS'. Both behave the same and this is especially true compared with versions < 3.1. As a consequence, using readLAS() provides the same performance (no degradation) than in previous versions, while using one of the read*LAS() functions may improve the performance.

Examples

LASfile <- system.file("extdata", "example.laz", package="rlas")
las <- readLAS(LASfile)

# By default the sensor and spatial index codes are 0
sensor(las)
index(las)

# Codes are used internally and not intended to be known by users
# Use h option for human readable output
sensor(las, h = TRUE)
index(las, h = TRUE)

# Modification of the sensor enables users to select a better spatial index
# when processing the point-cloud.
sensor(las) <- "tls"
sensor(las, h = TRUE)
index(las, h = TRUE)

# Modification of the spatial index forces users to choose one of the available
# spatial indexes.
index(las) <- "quadtree"
sensor(las, h = TRUE)
index(las, h = TRUE)

# The simplest way to take advantage of appropriate spatial indexing is
# to use one of the read*LAS() functions.
las <- readTLSLAS(LASfile)
sensor(las, h = TRUE)
index(las, h = TRUE)

# But for some specific point-clouds / algorithms it might be advisable to force
# the use of a specific spatial index to perform the computation faster
index(las) <- "voxelpartition"
index(las, h = TRUE)

# With a LAScatalog, spatial indexing information is propagated to the
# different chunks
ctg = readTLSLAScatalog(LASfile)
index(ctg) <- "voxelpartition"
sensor(ctg, h = TRUE)
index(ctg, h = TRUE)

# ==================
# PERFORMANCE TESTS
# ==================

## Not run: 
# Performance tests on TLS
# ------------------------

# The package does not include TLS data
# so we can generate something that looks TLS-ish
# >>>>>>>>>>>
X <- runif(50, -25, 25)
Y <- runif(50, -25, 25)
X <- as.numeric(sapply(X, function(x) rnorm(2000, x, 2)))
Y <- as.numeric(sapply(Y, function(x) rnorm(2000, x, 2)))
Z <- abs(rnorm(length(Y), 10, 5))
veg <- data.frame(X,Y,Z)
X <- runif(5000, -30, 30)
Y <- runif(5000, -30, 30)
Z <- runif(5000, 0, 1)
ground <- data.frame(X,Y,Z)
X <- runif(30, -30, 30)
Y <- runif(30, -30, 30)
Z <- runif(30, 0, 30)
noise <- data.frame(X,Y,Z)
las <- LAS(rbind(ground, veg, noise))
# <<<<<<<<<<<<<

plot(las)

# If read with readALSLAS()
sensor(las) <- "als"
system.time(classify_noise(las, sor(20, 8)))
#> 1.5 sec

# If read with readTLSLAS()
sensor(las) <- "tls"
system.time(classify_noise(las, sor(20, 8)))
#> 0.6 sec

# Performance tests on ALS
# ------------------------

# The package does not include large ALS data
# so we can generate something that looks ALS-ish
# >>>>>>>>>>>
X <- runif(4e5, 0, 1000)
Y <- runif(4e5, 0, 1000)
Z <- 40*sin(0.01*X) + 50*cos(0.005*Y) + abs(rnorm(length(Y), 10, 5))
veg <- data.frame(X,Y,Z)
X <- runif(100, 0, 1000)
Y <- runif(100, 0, 1000)
Z <- 40*sin(0.01*X) + 50*cos(0.005*Y) + abs(rnorm(length(Y), 10, 5)) + runif(100, 30, 70)
noise <- data.frame(X,Y,Z)
las <- LAS(rbind(veg, noise))
# <<<<<<<<<<<<<

plot(las)

# If read with readALSLAS()
sensor(las) <- "als"
system.time(classify_noise(las, sor(15, 8)))
#> 3 sec

# If read with readTLSLAS()
sensor(las) <- "tls"
system.time(classify_noise(las, sor(15, 8)))
#> 4.3 sec

## End(Not run)

Description

Individual tree detection function that find the position of the trees using several possible algorithms.

Usage

locate_trees(las, algorithm, uniqueness = "incremental")

Arguments

Value

locate_trees returns an sf object with POINT Z geometries. The table of attributes contains a column treeID with an individual ID for each tree. The height of the trees (Z) are also repeated in the table of attribute to be analysed as an attribute and not as a coordinate.

Uniqueness

By default the tree IDs are numbered from 1 to n, n being the number of trees found. The problem with such incremental numbering is that, while it ensures a unique ID is assigned for each tree in a given point-cloud, it also guarantees duplication of tree IDs in different tiles or chunks when processing a LAScatalog. This is because each chunk/file is processed independently of the others and potentially in parallel on different computers. Thus, the index always restarts at 1 on each chunk/file. Worse, in a tree segmentation process, a tree that is located exactly between 2 chunks/files will have two different IDs for its two halves.

This is why we introduced some uniqueness strategies that are all imperfect and that should be seen as experimental. Please report any troubleshooting. Using a uniqueness-safe strategy ensures that trees from different files will not share the same IDs. It also ensures that two halves of a tree on the edge of a processing chunk will be assigned the same ID.

incremental

Number from 0 to n. This method does not ensure uniqueness of the IDs. This is the legacy method.

gpstime

This method uses the gpstime of the highest point of a tree (apex) to create a unique ID. This ID is not an integer but a 64-bit decimal number, which is suboptimal but at least it is expected to be unique if the gpstime attribute is consistent across files. If inconsistencies with gpstime are reported (for example gpstime records the week time and was reset to 0 in a coverage that takes more than a week to complete), there is a (low) probability of getting ID attribution errors.

bitmerge

This method uses the XY coordinates of the highest point (apex) of a tree to create a single 64-bit number with a bitwise operation. First, XY coordinates are converted to 32-bit integers using the scales and offsets of the point cloud. For example, if the apex is at (10.32, 25.64) with a scale factor of 0.01 and an offset of 0, the 32-bit integer coordinates are X = 1032 and Y = 2564. Their binary representations are, respectively, (here displayed as 16 bits) 0000010000001000 and 0000101000000100. X is shifted by 32 bits and becomes a 64-bit integer. Y is kept as-is and the binary representations are unionized into a 64-bit integer like (here displayed as 32 bit) 00000100000010000000101000000100 that is guaranteed to be unique. However R does not support 64-bit integers. The previous steps are done at C++ level and the 64-bit binary representation is reinterpreted into a 64-bit decimal number to be returned in R. The IDs thus generated are somewhat weird. For example, the tree ID 00000100000010000000101000000100 which is 67635716 if interpreted as an integer becomes 3.34164837074751323479078607289E-316 if interpreted as a decimal number. This is far from optimal but at least it is guaranteed to be unique if all files have the same offsets and scale factors.

All the proposed options are suboptimal because they either do not guarantee uniqueness in all cases (inconsistencies in the collection of files), or they imply that IDs are based on non-integers or meaningless numbers. But at least it works and deals with some of the limitations of R.

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyz", filter = "-inside 481250 3812980 481300 3813030")

ttops <- locate_trees(las, lmf(ws = 5))

#plot(las) |> add_treetops3d(ttops)

Description

Merge a point cloud with a source of spatial data. It adds an attribute along each point based on a value found in the spatial data. Sources of spatial data can be a ⁠SpatialPolygons*⁠, an sf/sfc, a ⁠Raster*⁠, a stars, or a SpatRaster.

• ⁠SpatialPolygons*⁠, sf and sfc: it checks if the points belongs within each polygon. If the parameter attribute is the name of an attribute in the table of attributes it assigns to the points the values of that attribute. Otherwise it classifies the points as boolean. TRUE if the points are in a polygon, FALSE otherwise.

• RasterLayer, single band stars or single layer SpatRaster: it attributes to each point the value found in each pixel of the raster.

• RasterStack, RasterBrick, multibands stars or multilayer SpatRaster must have 3 layers for RGB colors. It colorizes the point cloud with RGB values.

Usage

merge_spatial(las, source, attribute = NULL)

Arguments

Value

a LAS object

Examples

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
shp     <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")

las   <- readLAS(LASfile, filter = "-keep_random_fraction 0.1")
lakes <- sf::st_read(shp, quiet = TRUE)

# The attribute "inlake" does not exist in the shapefile.
# Points are classified as TRUE if in a polygon
las    <- merge_spatial(las, lakes, "inlakes")     # New attribute 'inlakes' is added.
names(las)

forest <- filter_poi(las, inlakes == FALSE)
#plot(forest)

# The attribute "LAKENAME_1" exists in the shapefile.
# Points are classified with the values of the polygons
las <- merge_spatial(las, lakes, "LAKENAME_1")     # New column 'LAKENAME_1' is added.
names(las)