Calculating Forest Metrics¶
This notebook will explain how to calculate and visualize forest metrics. Code snippets for how to get started using PyForestScan can be found in the notebook getting-started-importing-preprocessing-dtm-chm, and in the documentation.
The example dataset is a one-square-kilometer tile derived from a 2019 aerial LiDAR survey of the Big Island of Hawaii. The data is captured over a dry forest and has a nominal pulse spacing of 0.35 meters. The data has been preprocessed to classify ground and vegetation points.
First we will load the pyforestscan functions that will be used.
from pyforestscan.handlers import read_lidar, create_geotiff
from pyforestscan.visualize import plot_pad, plot_metric
from pyforestscan.filters import filter_hag
from pyforestscan.calculate import generate_dtm, assign_voxels, calculate_pad, calculate_pai, calculate_fhd, calculate_chm
Import Data¶
We will begin by importing and preprocessing the data. It is important for a HeightAboveGround dimension to be defined, which can be accomplished by setting hag=True when importing data.
file_path = "../example_data/20191210_5QKB020880.laz"
arrays = read_lidar(file_path, "EPSG:32605", hag=True)
Next, we will use filter_hag() to remove any points below ground.
arrays = filter_hag(arrays)
points = arrays[0]
Create Voxels¶
We need to create voxels, and assign points to them. Voxel resolution must be given as a tuple with the format (x_res, y_res, z_res) where x_res, y_res, and z_res are the resolutions along the x, y, and z axes. Points within voxel bounds are assigned to them.
voxel_resolution = (5, 5, 1)
voxels, extent = assign_voxels(points, voxel_resolution)
Calculate Forest Metrics¶
Canopy Height Model¶
We will calculate and plot a canopy height model (CHM), showing the highest point above ground in each voxel column. calculate_chm will return the canopy height model along with its extent. Forest metrics are calculated in the same units as the data. In this case, since the data is in EPSG:32605, output will be in meters.
chm, extent = calculate_chm(points, voxel_resolution)
plot_metric('Canopy Height Model', chm, extent, metric_name='Height (m)', cmap='viridis', fig_size=None)
Plant Area Density¶
Next we will calculate and plot plant area density (PAD), which shows density of plant matter within each voxel.
pad = calculate_pad(voxels, voxel_resolution[-1])
plot_pad(pad, 5, axis='y', cmap='viridis')
Plant Area Index¶
You can calculate and plot plant area index (PAI), which quantifies the total plant surface area of above-ground plant matter (leaves, branches, and stems) per unit of ground area, which corresponds to the vertical summation of PAD.
pai = calculate_pai(pad)
plot_metric('Plant Area Index', pai, extent, metric_name='PAI', cmap='viridis', fig_size=None)
Foliage Height Diversity¶
We can calculate foliage height diversity (FHD). FHD is a measure of how plant material is distributed across the vertical layers of a canopy, as calculated using Shannon entropy.
fhd = calculate_fhd(voxels)
plot_metric('Foliage Height Diversity', fhd, extent, metric_name='FHD', cmap='viridis', fig_size=None)
