,
J. Antonio Guzmán Q. [cph] (Author of original version of voxelize()
function)| Title: | Fast 3D Lacunarity for Voxel Data |
|---|---|
| Description: | Calculates 3D lacunarity from voxel data. It is designed for use with point clouds generated from Light Detection And Ranging (LiDAR) scans in order to measure the spatial heterogeneity of 3-dimensional structures such as forest stands. It provides fast 'C++' functions to efficiently bin point cloud data into voxels and calculate lacunarity using different variants of the gliding-box algorithm originated by Allain & Cloitre (1991) <doi:10.1103/PhysRevA.44.3552>. |
| Authors: | Elliott Smeds [aut, cre, cph] ,
J. Antonio Guzmán Q. [cph] (Author of original version of voxelize()
function) |
| Maintainer: | Elliott Smeds <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.1.9000 |
| Built: | 2025-03-02 05:06:58 UTC |
| Source: | https://github.com/elliottsmeds/lacunr |
bounding_box() takes a table of voxel data and converts it into a
3-dimensional array, with the original voxels arranged in their correct
spatial positions inside of a 3-D "box" of empty voxels. This array can be
input directly into lacunarity() to generate lacunarity curves.
bounding_box(x, threshold = 0, edge_length = NULL)bounding_box(x, threshold = 0, edge_length = NULL)
x |
A ' |
threshold |
An integer specifying the minimum number of points to use
when determining if a voxel is occupied. The default is |
edge_length |
a numeric vector of length 3, specifying the X, Y, and Z
dimensions of each voxel. This argument should only be necessary when
supplying voxel data generated by a function other than |
bounding_box() relies on the spatial coordinates of the input data
to determine the dimensions of the resulting array. Noisy point cloud data
will often produce "outlier" voxels surrounding points that are far removed
from the bulk of the point cloud. These can drastically alter the output,
creating an array in which the occupied voxels are surrounded by a large
region of empty space. It is highly recommended that users supply a cutoff
value for threshold to ideally remove these outliers. More elaborate
filtering tools for trimming the point data before voxelization are
available in the lidR package.
A 3-dimensional integer array containing values of 0
or 1, representing the occupancy of a given voxel. Occupied voxels (1)
are arranged according to their relative positions in 3-dimensional space,
and fully encapsulated within a rectangular volume of unoccupied voxels
(0). The XYZ positions of the voxels are retained in the array
dimnames.
# Basic usage --------------------------------------------------------------- # simulate a diagonal line of points with XYZ coordinates pc <- data.frame(X = as.numeric(0:24), Y = as.numeric(0:24), Z = as.numeric(0:24)) # convert point data to cubic voxels of length 5 vox <- voxelize(pc, edge_length = c(5,5,5)) # convert to voxel array box <- bounding_box(vox) # Using lidR::voxel_metrics ------------------------------------------------- if (require("lidR")){ # reformat point data into rudimentary LAS object las <- suppressMessages(lidR::LAS(pc)) # convert to voxels of length 5 vox <- lidR::voxel_metrics(las, ~list(N = length(Z)), res = 5) # convert to voxel array box <- bounding_box(vox) }# Basic usage --------------------------------------------------------------- # simulate a diagonal line of points with XYZ coordinates pc <- data.frame(X = as.numeric(0:24), Y = as.numeric(0:24), Z = as.numeric(0:24)) # convert point data to cubic voxels of length 5 vox <- voxelize(pc, edge_length = c(5,5,5)) # convert to voxel array box <- bounding_box(vox) # Using lidR::voxel_metrics ------------------------------------------------- if (require("lidR")){ # reformat point data into rudimentary LAS object las <- suppressMessages(lidR::LAS(pc)) # convert to voxels of length 5 vox <- lidR::voxel_metrics(las, ~list(N = length(Z)), res = 5) # convert to voxel array box <- bounding_box(vox) }
This dataset contains point cloud data from two terrestrial LiDAR scans of a Northern California oak forest shortly before and after the Glass Fire, which burned some 27000 hectares of land in Sonoma and Napa counties between September 27 and October 20, 2020. The scans encompass an identical 24m by 24m rectangular plot at a study site within the Saddle Mountain Open Space Preserve in Sonoma County.
glassfireglassfire
A data table with 1,000,000 rows and 4 columns: X, Y, Z, and
Year
X,Y,Z
The XYZ spatial positions of each point, in meters. X and Y denote the East-West and North-South horizontal positions, respectively, while Z denotes the vertical position
YearThe year each LiDAR scan was taken, either 2020, immediately before the Glass Fire, or 2021, a few months after
The original terrain topography has been removed using digital elevation model (DEM)-based height normalization, and ground points removed by clipping all points below 0.25m. The raw point cloud data were normalized via voxelization at a resolution of 0.125m, and the results further down-sampled to make the dataset more compact. The X, Y, and Z coordinates were generated from the original Easting, Northing, and elevation by subtracting their minimum values.
glassfire is technically encoded as a lasmetrics3d object from the lidR
package. This class inherits from data.table(), but has the added benefit
that it can be rendered as a 3D rgl plot using lidR::plot.lasmetrics3d().
Plot lacunarity curve(s)
lac_plot(..., log = TRUE, group_names = NULL) lacnorm_plot(..., log = TRUE, group_names = NULL) hr_plot(..., group_names = NULL)lac_plot(..., log = TRUE, group_names = NULL) lacnorm_plot(..., log = TRUE, group_names = NULL) hr_plot(..., group_names = NULL)
... |
One or more |
log |
A Boolean. |
group_names |
A character |
A ggplot object displaying the lacunarity or H(r) curve(s). If
multiple curves are supplied, their ordering in the plot legend will
reflect the order they were listed in the function call.
# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # calculate lacunarity at all box sizes lac_curve <- lacunarity(a, box_sizes = "all") # plot raw lacunarity curve lac_plot(lac_curve)# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # calculate lacunarity at all box sizes lac_curve <- lacunarity(a, box_sizes = "all") # plot raw lacunarity curve lac_plot(lac_curve)
Generates lacunarity curves for a specified set of box sizes, using one
of two versions of the gliding-box algorithm
lacunarity(x, box_sizes = "twos", periodic = FALSE)lacunarity(x, box_sizes = "twos", periodic = FALSE)
x |
A 3-dimensional |
box_sizes |
Which box sizes to use for calculating lacunarity:
|
periodic |
A Boolean. Determines which boundary algorithm to use, the classic fixed boundary by Allain and Cloitre (default) or the periodic boundary algorithm introduced by Feagin et al. 2007. The latter is slightly slower but is more robust to edge effects. |
The raw values depend on the proportion of occupied voxels
within the data space. As a result, it is difficult to compare two spatial
patterns with different occupancy proportions because the curves will begin
at different y-intercepts. This is rectified by normalizing the curve,
typically by log-transforming it and dividing by the lacunarity value at
the smallest box size (i.e. ). lacunarity() outputs
both normalized and non-normalized curves for convenience.
The function also computes , a transformed lacunarity curve introduced
by Feagin 2003. rescales normalized in terms of the Hurst
exponent, where values greater than 0.5 indicate heterogeneity and values
less than 0.5 indicate homogeneity. Where describes a pattern's
deviation from translational invariance, describes its deviation from
standard Brownian motion.
A data.frame containing box sizes and their
corresponding , normalized , and values. Lacunarity is
always computed for box size 1, even if the user supplies a custom
box_sizes vector that omits it, as this value is required to calculate
normalized lacunarity.
Allain, C., & Cloitre, M. (1991). Characterizing the lacunarity of random and deterministic fractal sets. Physical Review A, 44(6), 3552–3558. doi:10.1103/PhysRevA.44.3552.
Feagin, R. A. (2003). Relationship of second-order lacunarity, Hurst exponent, Brownian motion, and pattern organization. Physica A: Statistical Mechanics and its Applications, 328(3-4), 315-321. doi:10.1016/S0378-4371(03)00524-7.
Feagin, R. A., Wu, X. B., & Feagin, T. (2007). Edge effects in lacunarity analysis. Ecological Modelling, 201(3–4), 262–268. doi:10.1016/j.ecolmodel.2006.09.019.
# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # calculate lacunarity with default options lacunarity(a) # supply custom vector of box sizes lacunarity(a, box_sizes = c(1,3,5)) # calculate lacunarity at all box sizes using the periodic boundary algorithm lacunarity(a, box_sizes = "all", periodic = TRUE)# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # calculate lacunarity with default options lacunarity(a) # supply custom vector of box sizes lacunarity(a, box_sizes = c(1,3,5)) # calculate lacunarity at all box sizes using the periodic boundary algorithm lacunarity(a, box_sizes = "all", periodic = TRUE)
pad_array() adds additional rows, columns, or slices to a 3-dimensional
array, increasing the array's dimensions by the desired amount and filling
the new space with a uniform value. It is intended for adding empty or
occupied space to the edges of a 3D spatial map.
pad_array(a, x = 0, y = 0, z = 0, fill = 0)pad_array(a, x = 0, y = 0, z = 0, fill = 0)
a |
A 3-dimensional |
x |
A positive or negative integer, denoting the number of rows to add to the array. The sign dictates which side of the array to pad. Default is zero. |
y |
A positive or negative integer, denoting the number of columns to add to the array. The sign dictates which side of the array to pad. Default is zero. |
z |
A positive or negative integer, denoting the number of "slices" to add to the array. The sign dictates which side of the array to pad. Default is zero. |
fill |
The desired value to fill the array padding with. Default is zero. |
pad_array() uses the signs of x, y, and z to determine where
to add padding. Negative values are prepended before the array's lower
indices (what one might call the "top", "front", or "beginning" of the
array), while positive values are appended after the upper indices (the
"bottom", "back", or "end" of the array).
A 3-dimensional array with the desired padding added.
The padded portions are labelled using their dimnames. If
no padding has been specified, the function returns the original array.
# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # add two rows of zeroes to top of array pad <- pad_array(a, x = -2) # add one row of zeroes to bottom of array, and two columns to beginning pad <- pad_array(a, x = 1, y = -2)# generate array a <- array(data = rep(c(1,0), 125), dim = c(5,5,5)) # add two rows of zeroes to top of array pad <- pad_array(a, x = -2) # add one row of zeroes to bottom of array, and two columns to beginning pad <- pad_array(a, x = 1, y = -2)
Bins point cloud data into 3D pixels, otherwise known as 'voxels'.
voxelize(x, edge_length, threads = 1L)voxelize(x, edge_length, threads = 1L)
x |
A |
edge_length |
A numeric |
threads |
The number of threads to use for computing the voxel data. Default is 1. |
A data object of class 'lac_voxels', which inherits from
data.table. The output contains 4 columns: X, Y, Z, and
N. The first three columns encode the spatial coordinates of each voxel
while the fourth denotes the total number of points they contain.
# simulate a diagonal line of points with XYZ coordinates pc <- data.frame(X = 0:99, Y = 0:99, Z = 0:99) # convert point data to cubic voxels of length 5 voxelize(pc, edge_length = c(5,5,5))# simulate a diagonal line of points with XYZ coordinates pc <- data.frame(X = 0:99, Y = 0:99, Z = 0:99) # convert point data to cubic voxels of length 5 voxelize(pc, edge_length = c(5,5,5))