Type:	Package
Title:	Utilities to Support Lidar Applications at the Landscape, Forest, and Tree Scale
Version:	1.0.2
Date:	2026-01-19
Description:	Implements algorithms for terrestrial, mobile, and airborne lidar processing, tree detection, segmentation, and attribute estimation (Donager et al., 2021) <doi:10.3390/rs13122297>, and a hierarchical patch delineation algorithm 'PatchMorph' (Girvetz & Greco, 2007) <doi:10.1007/s10980-007-9104-8>. Tree detection uses rasterized point cloud metrics (relative neighborhood density and verticality) combined with RANSAC cylinder fitting to locate tree boles and estimate diameter at breast height. Tree segmentation applies graph-theory approaches inspired by Tao et al. (2015) <doi:10.1016/j.isprsjprs.2015.08.007> with cylinder fitting methods from de Conto et al. (2017) <doi:10.1016/j.compag.2017.07.019>. PatchMorph delineates habitat patches across spatial scales using organism-specific thresholds. Built on 'lidR' (Roussel et al., 2020) <doi:10.1016/j.rse.2020.112061>.
URL:	https://github.com/bi0m3trics/spanner
License:	GPL-3
Encoding:	UTF-8
RoxygenNote:	7.3.3
RdMacros:	mathjaxr
LinkingTo:	lidR, RcppArmadillo, Rcpp (≥ 1.0.13), RcppEigen, BH,
Imports:	Rcpp (≥ 1.0.13), conicfit, FNN, RANN, cppRouting, sf, terra, sfheaders, Rfast, geometry, dplyr, mathjaxr, data.table
Depends:	R (≥ 4.0.0), lidR (≥ 4.2.0),
Suggests:	testthat (≥ 3.0.0), magick, rgl, rstac
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2026-01-29 02:44:22 UTC; ajsm
Author:	Andrew Sanchez Meador [aut, cre, ctb], Jonathon Donager [aut, ctb], Blackburn Ryan [aut, ctb], Cannon Jeffery [ctb], Tiago de Conto [ctb, cph] (Author of included TreeLS code), Keith O'Hara [ctb, cph] (Author of included OptimLib code)
Maintainer:	Andrew Sanchez Meador <Andrew.SanchezMeador@nau.edu>
Repository:	CRAN
Date/Publication:	2026-02-03 10:30:02 UTC

Colorize a LAS object using multiple methods

Description

Colors a LAS object using one of three methods: attribute-based coloring, raster-based RGB coloring, or PCV (Portion de Ciel Visible) ambient occlusion.

Usage

colorize_las(
  las,
  method = "attr",
  attribute_name = NULL,
  palette = c("black", "white"),
  raster_path = NULL,
  radius = 1,
  num_directions = 60,
  kernel_size = 5,
  pixel_size = 0.1,
  num_samples = 16,
  ncpu = 4
)

Arguments

Details

The function supports four coloring methods:

attr

Attribute-based coloring: normalizes attribute values and maps them to colors using the palette.

rgb

Raster-based coloring: extracts RGB values from georeferenced raster(s) that align with the point cloud. Requires matching CRS between LAS and raster. Can use a single 3-band RGB raster or three separate rasters.

pcv

PCV (Portion de Ciel Visible): computes 3D ambient occlusion by calculating sky visibility for each point. Based on the algorithm from Duguet & Girardeau-Montaut (2004). More accurate but slower than SSAO.

ssao

SSAO (Screen Space Ambient Occlusion): fast ambient occlusion using 2D depth map techniques. Projects points to a depth buffer and calculates occlusion based on depth differences. Much faster than PCV.

Value

A LAS object with updated R, G, and B fields based on the selected method.

Examples

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

# Attribute-based coloring
las_colored <- colorize_las(las, method="attr", attribute_name="Z",
                            palette=c("blue", "green", "yellow", "red"))

# Raster-based coloring with RGB file
rgb_file <- system.file("extdata", "UAS_Clip_RGB.tif", package="spanner")
las_colored <- colorize_las(las, method="rgb", raster_path=rgb_file)

# PCV ambient occlusion (slow, high quality)
las_colored <- colorize_las(las, method="pcv", radius=1.0,
                            num_directions=30, palette=c("black", "white"))

# SSAO ambient occlusion (faster alternative to PCV)
las_colored <- colorize_las(las, method="ssao", pixel_size=0.1,
                            kernel_size=5, num_samples=16, palette=c("black", "white"), ncpu=8)

Compute PCV (Portion de Ciel Visible) for point cloud

Description

Calculates ambient occlusion using sky visibility algorithm. Based on Duguet & Girardeau-Montaut (2004).

Usage

compute_pcv(las, radius = 1, num_directions = 60, ncpu = 4)

Arguments

Details

The PCV algorithm computes the visible portion of the sky from each point by:

• Casting rays in multiple directions around each point

• Computing the maximum elevation angle (horizon) in each direction

• Calculating the average sky visibility across all directions

This function uses an optimized C++ implementation with OpenMP parallelization for improved performance on large point clouds.

Value

Numeric vector of PCV values (sky visibility) for each point

Compute SSAO (Screen Space Ambient Occlusion) for point cloud

Description

Fast ambient occlusion using 2D depth buffer technique. Much faster than PCV as it works on projected depth maps.

Usage

compute_ssao(
  las,
  kernel_size = 5,
  pixel_size = 0.1,
  num_samples = 16,
  ncpu = 4
)

Arguments

Details

The SSAO algorithm computes ambient occlusion by:

• Projecting the point cloud to a 2D depth map (grid)

• For each point, sampling the depth buffer around it

• Calculating occlusion based on depth differences

• Applying distance and angle-based falloff

This is significantly faster than full 3D PCV because it only requires 2D image processing operations rather than 3D neighbor searches and ray tracing.

Value

Numeric vector of SSAO values (ambient occlusion) for each point

Create animated GIF of rotating 3D point cloud

Description

Generates a 360-degree rotating animation of a LAS point cloud using rgl and saves it as an animated GIF.

Usage

create_rotation_gif(
  las,
  output_path = "pointcloud_rotation.gif",
  duration = 12,
  rpm = 5,
  background = "white",
  axis = "z",
  show_axis = TRUE,
  show_legend = TRUE,
  screen_size = c(800, 600),
  overwrite = FALSE
)

Arguments

Details

This function creates a smooth 360-degree rotation animation by:

• Plotting the point cloud using lidR's plot function with RGB colors

• Using rgl's movie3d and spin3d to create smooth rotation

• Saving the result as an animated GIF

The rotation speed is controlled by the rpm (rotations per minute) parameter. The total duration determines how long the animation will be.

Requires the rgl package and lidR for plotting.

Value

Character string of the output file path (invisible)

Examples

# Load example LAS file
LASfile <- system.file("extdata", "ALS_Clip.laz", package="spanner")
las <- readLAS(LASfile)

# Create basic rotation GIF with attribute coloring
las_colored <- colorize_las(las, method="attr", attribute_name="Z")
create_rotation_gif(las_colored, output_path=tempfile(fileext = ".gif"))

# High quality with specific settings
create_rotation_gif(las_colored,
                    output_path=tempfile(fileext = ".gif"),
                    duration=15,
                    rpm=10,
                    background="black",
                    show_axis=FALSE,
                    show_legend=FALSE)

# Rotate around X axis for side-to-side view
create_rotation_gif(las_colored, output_path=tempfile(fileext = ".gif"), axis="x")

Point cloud cylinder fitting as per de Conto et al. 2017 as implemented here: https://github.com/tiagodc/TreeLS

Description

Fits a cylinder on a set of 3D points.

Usage

cylinderFit(
  las,
  method = "ransac",
  n = 5,
  inliers = 0.9,
  conf = 0.95,
  max_angle = 30,
  n_best = 20
)

Arguments

Value

vector of parameters

Examples

# Define the cylinder attributes
npts = 500
cyl_length = 0.5
radius = 0.2718

# Generate the X,Y,Z values
Z=runif(n = npts, min = 0, max = cyl_length)
angs = runif(npts, 0, 2*pi)
X = sin(angs)*radius
Y = cos(angs)*radius

# Creation of a LAS object out of external data
cloud <- LAS(data.frame(X,Y,Z))

# Fit a cylinder and retrun the information
cyl_par = spanner::cylinderFit(cloud, method = 'ransac', n=5, inliers=.9,
                               conf=.95, max_angle=30, n_best=20)

Download NAIP Imagery for LiDAR Extent

Description

Downloads NAIP (National Agriculture Imagery Program) imagery from Microsoft Planetary Computer STAC API for the extent of a LAS/LAZ point cloud.

Usage

download_naip_for_las(
  las,
  output_path = NULL,
  year_range = c("2018-01-01", "2023-12-31"),
  buffer = 0,
  overwrite = FALSE
)

Arguments

Details

This function queries the Microsoft Planetary Computer STAC API to find and download NAIP imagery that overlaps with the extent of the input LAS file. The imagery is automatically cropped to match the LiDAR extent and saved as a GeoTIFF.

NAIP imagery is typically 4-band (RGB + NIR) with 0.6m or 1m resolution, collected annually or biannually across the continental United States.

Requires the rstac package for STAC API access.

Value

Character string of the output file path, or NULL if download failed.

Examples

# Load example LAS file
LASfile <- system.file("extdata", "ALS_Clip.laz", package="spanner")
las <- readLAS(LASfile)

# Download NAIP for a LAS file
naip_path <- download_naip_for_las(las, output_path = tempfile(fileext = ".tif"))

# Download with buffer and specific year range
naip_path2 <- download_naip_for_las(las, buffer = 10,
                                    output_path = tempfile(fileext = ".tif"),
                                    year_range = c("2020-01-01", "2023-12-31"))

# Then use with colorize_las
las_colored <- colorize_las(las, method = "rgb", raster_path = naip_path)

Calculates eigen decomposition metrics for fixed neighborhood point cloud data

Description

This function calculates twelve (plus the first and second PCA) for several point geometry-related metrics (listed below) in parallel using C++ for a user-specified radius.

Usage

eigen_metrics(las = las, radius = 0.1, ncpu = 8)

Arguments

Value

A labeled data.table of point metrics for each point in the LAS object

List of available point metrics

• eLargest: first eigenvalue, \(\lambda_{1}\)

• eMedium: second eigenvalue, \(\lambda_{2}\)

• eSmallest: third eigenvalue, \(\lambda_{3}\)

• eSum: sum of eigenvalues, \(\sum_{i=1}^{n=3} \lambda_{i}\)

• Curvature: surface variation, \(\lambda_{3} / \sum_{i=1}^{n=3} \lambda_{i}\)

• Omnivariance: high values correspond to spherical features and low values to planes or linear features, \((\lambda_{1} * \lambda_{2} * \lambda_{3})^{1/3}\)

• Anisotropy: relationships between the directions of the point distribution, \((\lambda_{1} - \lambda_{3}) / \lambda_{1}\)

• Eigentropy: entropy in the eigenvalues, \(- \sum_{i=1}^{n=3} \lambda_{i} * ln(\lambda_{i})\)

• Linearity: linear saliency, \((\lambda_{1} - \lambda_{2}) / \lambda_{1}\)

• Verticality: vertical saliency, \(1-abs(\langle (0,0,1),e_3\rangle)\)

• Planarity: planar saliency, \((\lambda_{2} - \lambda_{3}) / \lambda_{1}\)

• Sphericity: spherical saliency, \(\lambda_{3} / \lambda_{1}\)

• Nx,Ny,Nz: 3 components of the normal vector (smallest eigenvector)

• SurfaceVariation: surface variation (change of curvature), same as Curvature

• ChangeCurvature: alternative name for surface variation

• SurfaceDensity: 2D point density using circle area, \(k / (\pi R^{2})\)

• VolumeDensity: 3D point density using sphere volume, \(k / (\frac{4}{3}\pi R^{3})\)

• MomentOrder1: 1st order moment from CloudCompare, projection onto 2nd eigenvector, \(m_{1}^{2} / m_{2}\)

• NormalChangeRate: normal change rate, same as Curvature, \(\lambda_{3} / \sum_{i=1}^{n=3} \lambda_{i}\)

• Roughness: distance from query point to fitted plane, \(|\vec{d} \cdot \vec{n}|\)

• MeanCurvature: mean curvature from quadric surface fitting, \(H = \frac{(1+f_{y}^{2})f_{xx} - 2f_{x}f_{y}f_{xy} + (1+f_{x}^{2})f_{yy}}{2(1+f_{x}^{2}+f_{y}^{2})^{3/2}}\)

• GaussianCurvature: Gaussian curvature from quadric surface fitting, \(K = \frac{f_{xx}f_{yy} - f_{xy}^{2}}{(1+f_{x}^{2}+f_{y}^{2})^{2}}\)

• PCA1: eigenvector projection variance normalized by eigensum, \(\sigma_{PC1}^{2} / \sum_{i=1}^{n=3} \lambda_{i}\)

• PCA2: eigenvector projection variance normalized by eigensum, \(\sigma_{PC2}^{2} / \sum_{i=1}^{n=3} \lambda_{i}\)

• NumNeighbors: number of points in the spherical neighborhood, \(k\)

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
eigen = eigen_metrics(las, radius=2, ncpu=4)

Obtain tree information by rasterizing point cloud values of relative neighborhood density and verticality within a slice of a normalized point cloud

Description

get_raster_eigen_treelocs returns a data.frame containing TreeID, X, Y, Z, Radius and Error in the same units as the .las

Usage

get_raster_eigen_treelocs(
  las = las,
  res = 0.05,
  pt_spacing = 0.0254,
  dens_threshold = 0.2,
  neigh_sizes = c(0.333, 0.166, 0.5),
  eigen_threshold = 0.6666,
  grid_slice_min = 0.6666,
  grid_slice_max = 2,
  minimum_polygon_area = 0.025,
  cylinder_fit_type = "ransac",
  max_dia = 0.5,
  SDvert = 0.25,
  n_best = 25,
  n_pts = 20,
  inliers = 0.9,
  conf = 0.99,
  max_angle = 20
)

Arguments

Details

For terrestrial and mobile lidar datasets, tree locations and estimates of DBH are provided by rasterizing individual point cloud values of relative neighborhood density (at 0.3 and 1 m radius) and verticality within a slice of the normalized point cloud around breast height (1.34 m). The algorithim then uses defined threshold values to classify the resulting rasters and create unique polygons from the resulting classified raster. These point-density and verticality polygons were selected by their intersection with one another, resulting in a final set of polygons which were used to clip out regions of the point cloud that were most likely to represent tree boles. A RANSAC cylinder fitting algorithm was then used to estimate the fit of a cylinder to individual bole points. Cylinder centers and radius were used as inputs to an individual tree segmentation

Value

sf A sf object containing the following tree seed information: TreeID, Radius, and Error in the same units as the .las, as well as the point geometry

Examples

# Set the number of threads to use in lidR
set_lidr_threads(8)

LASfile = system.file("extdata", "TLS_Clip.laz", package="spanner")
las = readTLSLAS(LASfile, select = "xyzcr", "-filter_with_voxel 0.01")
# Don't forget to make sure the las object has a projection
sf::st_crs(las) <- 26912

# Pre-process the example lidar dataset by classifying the ground  and noise points
# using lidR::csf(), normalizing it, and removing outlier points
# using lidR::ivf()
# las = classify_ground(las, csf(sloop_smooth = FALSE,
#                                 class_threshold = 0.5,
#                                cloth_resolution = 0.5, rigidness = 1L,
#                                 iterations = 500L, time_step = 0.65))
# las = normalize_height(las, tin())
# las = classify_noise(las, ivf(0.25, 3))
# las = filter_poi(las, Classification != LASNOISE)

# Plot the non-ground points, colored by height
# plot(filter_poi(las, Classification != 2), color = "Z")

# find tree locations and attribute data
myTreeLocs = get_raster_eigen_treelocs(las = las, res = 0.025, pt_spacing = 0.0254,
                                       dens_threshold = 0.25,
                                       neigh_sizes = c(0.25, 0.15, 0.66),
                                       eigen_threshold = 0.75,
                                       grid_slice_min = 1,
                                       grid_slice_max = 2,
                                       minimum_polygon_area = 0.005,
                                       cylinder_fit_type = "ransac",
                                       max_dia = 1,
                                       SDvert = 0.33,
                                       n_pts = 20,
                                       n_best = 25,
                                       inliers = 0.9,
                                       conf = 0.99,
                                       max_angle = 20)

# Plot results if trees were found
if (!is.null(myTreeLocs) && nrow(myTreeLocs) > 0) {
  plot(lidR::rasterize_canopy(las, res = 0.2, p2r()))
  symbols(sf::st_coordinates(myTreeLocs)[,1], sf::st_coordinates(myTreeLocs)[,2],
          circles = myTreeLocs$Radius^2*3.14, inches = FALSE, add = TRUE, bg = 'black')
} else {
  message("No tree locations were found. Try adjusting the parameters.")
}

Convert LAS object to XYZ matrix

Description

Extracts the X, Y, and Z coordinates from a LAS object and returns them as a matrix.

Usage

las2xyz(las)

Arguments

Value

A numeric matrix with three columns (X, Y, Z) containing the point coordinates

Examples

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
xyz_matrix <- las2xyz(las)
head(xyz_matrix)

Merge RGB colors from two colorized LAS objects

Description

Blends the RGB values from two LAS objects to create a new composite coloring. Useful for combining different coloring methods (e.g., ambient occlusion with raster RGB).

Usage

merge_las_colors(las1, las2, alpha = 0.5, method = "alpha")

Arguments

Details

Blending methods:

alpha

Simple linear interpolation: (1-alpha)las1 + alphalas2

multiply

Multiplicative blend (darkens): (las1 * las2) / 255

screen

Screen blend (lightens): 255 - ((255-las1) * (255-las2)) / 255

overlay

Overlay blend: combines multiply and screen based on base color

Common use cases:

• Combine ambient occlusion (PCV/SSAO) with aerial RGB for realistic shading

• Blend attribute coloring with terrain colors

• Overlay multiple visualization layers

Value

A LAS object (copy of las1) with merged R, G, and B fields

Examples

# Load example LAS file
LASfile <- system.file("extdata", "ALS_Clip.laz", package="spanner")
las <- readLAS(LASfile)

# Combine SSAO ambient occlusion with aerial RGB
las_ao <- colorize_las(las, method="ssao", palette=c("black", "white"))
rgb_file <- system.file("extdata", "UAS_Clip_RGB.tif", package="spanner")
las_rgb <- colorize_las(las, method="rgb", raster_path=rgb_file)
las_merged <- merge_las_colors(las_ao, las_rgb, alpha=0.3, method="multiply")

# Blend attribute coloring with RGB at 50/50
las_height <- colorize_las(las, method="attr", attribute_name="Z",
                           palette=c("blue", "red"))
las_merged <- merge_las_colors(las_height, las_rgb, alpha=0.5)

Plot a raster by its name

Description

plot_raster_by_name plots a raster from a list of rasters based on the provided raster name.

Usage

plot_raster_by_name(rasters, raster_name)

Arguments

Value

NULL This function does not return a value. It plots the raster if found.

Examples

# Define input parameters
las <- lidR::readLAS(system.file("extdata", "MixedConifer.laz", package="lidR"))
input_raster <- lidR::rasterize_canopy(las, res = 1, lidR::pitfree(c(0,2,5,10,15), c(0, 2)))
suitList <- c(0, 2, 32)
gapList <- seq(1, 8, by = 1)
spurList <- seq(1, 8, by = 1)

# Process the rasters
processed_rasters <- process_rasters_patchmorph(input_raster, suitList, gapList, spurList)

# Plot a raster by its name
plot_raster_by_name(processed_rasters, "suit_2_gap_2_spur_6")

Process rasters based on suitability, gap, and spur parameters

Description

process_rasters_patchmorph processes an input raster by reclassifying it based on suitability levels, and then applying gap and spur distance transformations to generate a list of processed rasters.

Usage

process_rasters_patchmorph(input_raster, suitList, gapList, spurList)

Arguments

Value

list A list of processed rasters with names indicating the suitability, gap, and spur parameters used.

Examples

# Define input parameters
las <- lidR::readLAS(system.file("extdata", "MixedConifer.laz", package="lidR"))
input_raster <- lidR::rasterize_canopy(las, res = 1, lidR::pitfree(c(0,2,5,10,15), c(0, 2)))
suitList <- c(0, 2, 32)
gapList <- seq(1, 8, by = 1)
spurList <- seq(1, 8, by = 1)

# Process the rasters
processed_rasters <- process_rasters_patchmorph(input_raster, suitList, gapList, spurList)

# Plot the first processed raster
plot(processed_rasters[[1]])

Obtain tree information by processing point cloud data

Description

process_tree_data processes the output of get_raster_eigen_treelocs and segment_graph to add information about the height, crown area, and diameter for each unique TreeID. It also has an optional parameter to return an sf object representing the convex hulls for each tree.

Usage

process_tree_data(treeData, segmentedLAS, return_sf = FALSE)

Arguments

Details

For terrestrial and mobile lidar datasets, tree locations and estimates of DBH are provided by rasterizing individual point cloud values of relative neighborhood density (at 0.3 and 1 m radius) and verticality within a slice of the normalized point cloud around breast height (1.34 m). The algorithm then uses defined threshold values to classify the resulting rasters and create unique polygons from the resulting classified raster. These point-density and verticality polygons were selected by their intersection with one another, resulting in a final set of polygons which were used to clip out regions of the point cloud that were most likely to represent tree boles. A RANSAC cylinder fitting algorithm was then used to estimate the fit of a cylinder to individual bole points. Cylinder centers and radius were used as inputs to an individual tree segmentation.

Value

sf object An updated sf object with the original columns plus:

height

numeric: Height of the highest point for each TreeID.

crown_area

numeric: Area of the convex hull for each TreeID.

crown_base_height

numeric: Height to the base of the live crown for each TreeID.

crown_volume

numeric: Volume of the convex hull for the crown of each TreeID.

diameter

numeric: Diameter of the tree, calculated as twice the Radius.

If return_sf is TRUE, returns an sf object where the geometry is the convex hulls for each tree. If return_sf is FALSE, returns an sf object with point geometries using treeData.

Examples

# Set the number of threads to use in lidR
set_lidr_threads(8)

LASfile = system.file("extdata", "TLS_Clip.laz", package="spanner")
las = readTLSLAS(LASfile, select = "xyzcr", "-filter_with_voxel 0.01")
# Don't forget to make sure the las object has a projection
sf::st_crs(las) <- 26912

# Pre-process the example lidar dataset by classifying the ground  and noise points
# using lidR::csf(), normalizing it, and removing outlier points
# using lidR::ivf()
# las = classify_ground(las, csf(sloop_smooth = FALSE,
#                                 class_threshold = 0.5,
#                                cloth_resolution = 0.5, rigidness = 1L,
#                                 iterations = 500L, time_step = 0.65))
# las = normalize_height(las, tin())
# las = classify_noise(las, ivf(0.25, 3))
# las = filter_poi(las, Classification != LASNOISE)

# Plot the non-ground points, colored by height
# plot(filter_poi(las, Classification != 2), color = "Z")

# Find individual tree locations and attribute data
# find tree locations and attribute data
myTreeLocs = get_raster_eigen_treelocs(las = las, res = 0.025, pt_spacing = 0.0254,
                                       dens_threshold = 0.25,
                                       neigh_sizes = c(0.25, 0.15, 0.66),
                                       eigen_threshold = 0.75,
                                       grid_slice_min = 1,
                                       grid_slice_max = 2,
                                       minimum_polygon_area = 0.005,
                                       cylinder_fit_type = "ransac",
                                       max_dia = 1,
                                       SDvert = 0.33,
                                       n_pts = 20,
                                       n_best = 25,
                                       inliers = 0.9,
                                       conf = 0.99,
                                       max_angle = 20)

# Plot results if trees were found
if (!is.null(myTreeLocs) && nrow(myTreeLocs) > 0) {
  plot(lidR::rasterize_canopy(las, res = 0.2, p2r()))
  symbols(sf::st_coordinates(myTreeLocs)[,1], sf::st_coordinates(myTreeLocs)[,2],
          circles = myTreeLocs$Radius^2*3.14, inches = FALSE, add = TRUE, bg = 'black')
} else {
  message("No tree locations were found. Try adjusting the parameters.")
}

# Segment the point cloud
# For areas with interlocking crowns and trees of different sizes,
# enable metabolic scaling to prevent height overestimation
myTreeGraph = segment_graph(las = las, tree.locations = myTreeLocs, k = 50,
                             distance.threshold = 0.5,
                             use.metabolic.scale = FALSE,
                             ptcloud_slice_min = 1,
                             ptcloud_slice_max = 2,
                             subsample.graph = 0.1,
                             return.dense = TRUE)

# Plot it in 3D colored by treeID
plot(myTreeGraph, color = "treeID", pal=spanner_pal())

# Process the data
processed_data <- process_tree_data(myTreeLocs, myTreeGraph, return_sf = TRUE)

# Print the processed data
print(processed_data$data)
# Print the sf object if return_sf is TRUE
if (!is.null(processed_data$sf)) {
  print(processed_data$sf)
}

Segment a terrestrial point cloud using graph theory.

Description

segment_graph returns a .las object with a new column "treeID".

Usage

segment_graph(
  las,
  tree.locations,
  k = 50,
  distance.threshold = 0.33,
  use.metabolic.scale = FALSE,
  ptcloud_slice_min = 0.5,
  ptcloud_slice_max = 2,
  metabolic.scale.function = NULL,
  subsample.graph = 0.1,
  return.dense = FALSE
)

Arguments

Details

Preforms Individual tree segmentation following ecological principles for “growing” trees based on these input locations in a graph-theory approach inspired by work of Tao and others (2015). Point coordinates are linked together based on proximity and turned into a connected graph object, using the estimated tree bole locations as origin points, connecting individual points back to those tree bole origins based on shortest paths within the graph network, and finally assigning those points a unique tree identification based on the bole coordinate for which they are connected. Input point cloud is subsampled to a lower resolution before processing to increase processing efficiency. However, graph objects can still get large quite rapidly. Take this into consideration when choosing the extent of the input las object.

Value

a sparse/dense normalized .las with the column treeID added.

References

Tao, S., Wu, F., Guo, Q., Wang, Y., Li, W., Xue, B., ... & Fang, J. (2015). Segmenting tree crowns from terrestrial and mobile LiDAR data by exploring ecological theories. ISPRS Journal of Photogrammetry and Remote Sensing, 110, 66-76.

Examples

# Set the number of threads to use in lidR
set_lidr_threads(8)

LASfile = system.file("extdata", "TLS_Clip.laz", package="spanner")
las = readTLSLAS(LASfile, select = "xyzcr", "-filter_with_voxel 0.01")
# Don't forget to make sure the las object has a projection
sf::st_crs(las) <- 26912

# Pre-process the example lidar dataset by classifying the ground  and noise points
# using lidR::csf(), normalizing it, and removing outlier points
# using lidR::ivf()
# las = classify_ground(las, csf(sloop_smooth = FALSE,
#                                 class_threshold = 0.5,
#                                cloth_resolution = 0.5, rigidness = 1L,
#                                 iterations = 500L, time_step = 0.65))
# las = normalize_height(las, tin())
# las = classify_noise(las, ivf(0.25, 3))
# las = filter_poi(las, Classification != LASNOISE)

# Plot the non-ground points, colored by height
# plot(filter_poi(las, Classification != 2), color = "Z")

# Perform a deep inspection of the las object. If you see any
# red text, you may have issues!
las_check(las)

# Find individual tree locations and attribute data
# find tree locations and attribute data
myTreeLocs = get_raster_eigen_treelocs(las = las, res = 0.025, pt_spacing = 0.0254,
                                       dens_threshold = 0.25,
                                       neigh_sizes = c(0.25, 0.15, 0.66),
                                       eigen_threshold = 0.75,
                                       grid_slice_min = 1,
                                       grid_slice_max = 2,
                                       minimum_polygon_area = 0.005,
                                       cylinder_fit_type = "ransac",
                                       max_dia = 1,
                                       SDvert = 0.33,
                                       n_pts = 20,
                                       n_best = 25,
                                       inliers = 0.9,
                                       conf = 0.99,
                                       max_angle = 20)

# Plot results if trees were found
if (!is.null(myTreeLocs) && nrow(myTreeLocs) > 0) {
  plot(lidR::rasterize_canopy(las, res = 0.2, p2r()))
  symbols(sf::st_coordinates(myTreeLocs)[,1], sf::st_coordinates(myTreeLocs)[,2],
          circles = myTreeLocs$Radius^2*3.14, inches = FALSE, add = TRUE, bg = 'black')
} else {
  message("No tree locations were found. Try adjusting the parameters.")
}

# Segment the point cloud
# For areas with interlocking crowns and trees of different sizes,
# enable metabolic scaling to prevent height overestimation
myTreeGraph = segment_graph(las = las, tree.locations = myTreeLocs, k = 50,
                             distance.threshold = 0.5,
                             use.metabolic.scale = FALSE,
                             ptcloud_slice_min = 1,
                             ptcloud_slice_max = 2,
                             subsample.graph = 0.1,
                             return.dense = TRUE)

# Plot it in 3D colored by treeID
plot(myTreeGraph, color = "treeID", pal=spanner_pal())

# Optional: Use a custom metabolic scaling function
# myTreeGraph = segment_graph(las = las, tree.locations = myTreeLocs, k = 50,
#                              distance.threshold = 0.5,
#                              use.metabolic.scale = TRUE,
#                              metabolic.scale.function = '1/((2*x)^(1/8))',
#                              ptcloud_slice_min = 1,
#                              ptcloud_slice_max = 2,
#                              subsample.graph = 0.1,
#                              return.dense = TRUE)

Spanner color palette

Description

Returns a named vector of colors for use in spanner visualizations. The palette includes 10 distinct colors suitable for categorical data visualization.

Usage

spanner_pal()

Value

A named character vector of hex color codes

Examples

# Get the palette
colors <- spanner_pal()

# Use in a plot
barplot(1:10, col = spanner_pal(), names.arg = names(spanner_pal()), las = 2)

Sum rasters by suitability level

Description

sum_rasters_by_suitability sums rasters from a list based on their suitability levels.

Usage

sum_rasters_by_suitability(rasters, suitList)

Arguments

Value

list A list of summed rasters for each suitability level.

Examples

# Define input parameters
las <- lidR::readLAS(system.file("extdata", "MixedConifer.laz", package="lidR"))
input_raster <- lidR::rasterize_canopy(las, res = 1, lidR::pitfree(c(0,2,5,10,15), c(0, 2)))
suitList <- c(0, 2, 32)
gapList <- seq(1, 8, by = 1)
spurList <- seq(1, 8, by = 1)

# Process the rasters
processed_rasters <- process_rasters_patchmorph(input_raster, suitList, gapList, spurList)

# Sum rasters by suitability level
summed_rasters <- sum_rasters_by_suitability(processed_rasters, suitList)

# Call the plot_raster_by_name function to plot the raster named "suit_2_sum"
plot_raster_by_name(summed_rasters, "suit_2_sum")