lazer.structure¶
LiDAR data is typically delivered in point format, where a series of xyz locations describe all of the discrete light interactions a laser beam makes between the aircraft and the surface.
By aggregating these point returns over a large horizontal area, say 10x10 meters, and binning all of the returns into discrete vertical classes (0-1m, 1-2m, 2-3m, ...), we can construct a continuous profile of vegetation structure. The structure module provides tools for working with these profile signals.
The profile we use to calculate these metrics is derived from the slicer dataset, which we create using the laser slicer command line tool. Each band of the slicer rasters represents the percent of canopy material present in each vertical bin (b1 maps to 0-1m, b2 maps to 1-2m, b3 maps to 2-3m, ...).
The routines below compute a series of forest vertical structure metrics from slicer data.
Functions for computing tree crowns and canopy structural metrics from lidar vertical profiles.
compute_base_height(layers)
¶
Returns the lowest layer location.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
layers |
ndarray |
array of layer position indices from |
required |
Returns:
| Type | Description |
|---|---|
int |
Base height. |
Source code in lazer/structure.py
def compute_base_height(layers: np.ndarray) -> int:
"""Returns the lowest layer location.
Args:
layers: array of layer position indices from `compute_layers()`.
Returns:
Base height.
"""
return np.min(layers)
compute_flame_gap(understory, overstory)
¶
Computes the distance between overstory and understory base heights.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
understory |
int |
understory base height from |
required |
overstory |
int |
overstory base height from |
required |
Returns:
| Type | Description |
|---|---|
int |
Flame gap. |
Source code in lazer/structure.py
def compute_flame_gap(understory: int, overstory: int) -> int:
"""Computes the distance between overstory and understory base heights.
Args:
understory: understory base height from `compute_understory_base()`
overstory: overstory base height from `compute_overstory_base()`.
Returns:
Flame gap.
"""
return np.max((overstory - understory, 0))
compute_layer_count(layers)
¶
Returns the number of distinct vertical layers.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
layers |
ndarray |
array of layer position indices from |
required |
Returns:
| Type | Description |
|---|---|
int |
Layer count. |
Source code in lazer/structure.py
def compute_layer_count(layers: np.ndarray) -> int:
"""Returns the number of distinct vertical layers.
Args:
layers: array of layer position indices from `compute_layers()`.
Returns:
Layer count.
"""
return len(layers)
compute_layers(profile, threshold=3)
¶
Identifies the locations of each unique vertical layer in a canopy profile.
Analyzes the vertical profile like a spectrum, looking for changes in signs between positive and negative changes in the shape of the profile.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
profile |
ndarray |
a 1-d array with the vertical profile information derived from slicer data. |
required |
threshold |
int |
the layer-to-layer change in density required to count as a new layer. |
3 |
Returns:
| Type | Description |
|---|---|
ndarray |
the indices of profile locations where the start of a layer was identified. |
Source code in lazer/structure.py
def compute_layers(
profile: np.ndarray, threshold: int = defaults.vertical_layer_threshold
) -> np.ndarray:
"""Identifies the locations of each unique vertical layer in a canopy profile.
Analyzes the vertical profile like a spectrum, looking for changes in signs between
positive and negative changes in the shape of the profile.
Args:
profile: a 1-d array with the vertical profile information derived from slicer data.
threshold: the layer-to-layer change in density required to count as a new layer.
Returns:
the indices of profile locations where the start of a layer was identified.
"""
# calculate the change in sign as the difference between each vertical stratum
n_spec = len(profile)
strata = np.arange(n_spec)
offset = np.zeros(n_spec)
offset[0:-1] = profile[1:]
vchange = offset - profile
# calculate zero and positive change between strata to evaluate multi-layer structure
pchange = vchange >= 0
zchange = vchange == 0
# evaluate cumulative increases in canopy volume to handle continuous growth
for i in strata[1:]:
if pchange[i]:
vchange[i] += np.max((0, vchange[i - 1]))
vchange[i - 1] = 0
# set changes below a threshold to 0 to reduce noise
vchange[np.abs(vchange) < threshold] = 0
# if the layer is calculated above a 0 change layer, move the index down one
zoffset = (vchange > 0) * zchange
while zoffset.any():
zs = (np.where(zoffset))[0]
for z in zs:
vchange[z - 1] = vchange[z]
vchange[z] = 0
zchange[z] = False
zoffset = (vchange > 0) * zchange
# and return a layer as only areas where there is positive signal change
layers = np.where(vchange > 0)[0] + 1
return layers
compute_overstory_base(layers, over_under=4)
¶
Returns the lowest layer location in the canopy overstory.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
layers |
ndarray |
array of layer position indices from |
required |
over_under |
int |
height that distinguishes between under and overstory. |
4 |
Returns:
| Type | Description |
|---|---|
int |
Overstory base height. |
Source code in lazer/structure.py
def compute_overstory_base(
layers: np.ndarray, over_under: int = defaults.overstory_threshold
) -> int:
"""Returns the lowest layer location in the canopy overstory.
Args:
layers: array of layer position indices from `compute_layers()`.
over_under: height that distinguishes between under and overstory.
Returns:
Overstory base height.
"""
over = layers[layers > over_under]
ilayers = np.insert(over, 0, 0)
return np.max(ilayers)
compute_understory_base(layers, over_under=4)
¶
Returns the lowest layer location in the canopy understory.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
layers |
ndarray |
array of layer position indices from |
required |
over_under |
int |
height that distinguishes between under and overstory. |
4 |
Returns:
| Type | Description |
|---|---|
int |
Understory base height. |
Source code in lazer/structure.py
def compute_understory_base(
layers: np.ndarray, over_under: int = defaults.overstory_threshold
) -> int:
"""Returns the lowest layer location in the canopy understory.
Args:
layers: array of layer position indices from `compute_layers()`.
over_under: height that distinguishes between under and overstory.
Returns:
Understory base height.
"""
ilayers = np.insert(layers, 0, 0)
return np.max(ilayers[ilayers <= over_under])
crown_seg_dalponte_coomes(height, tree_top_x, tree_top_y, min_inclusion_height=2, height_decrease_rate=0.5, height_grouping_rate=0.5, max_diameter=8)
¶
Delineates crown vectors using the Dalponte/Coomes approach.
besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12575
Method described in the Supporting Information Appendix S2.
This routine calls a numba function to run the operation, which
doesn't appear to allow specifying default keyword arguments.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
height |
ndarray |
2d array with tree height |
required |
tree_top_x |
ndarray |
x index locations for tree tops |
required |
tree_top_y |
ndarray |
y index locations for tree tops |
required |
min_inclusion_height |
float |
minimum height threshold for inclusion in crown |
2 |
height_decrease_rate |
float |
relative pixel-to-pixel decrease allowed for inclusion in a neighboring crown |
0.5 |
height_grouping_rate |
float |
relative pixel-to-pixel increase allowed for inclusion in a neighboring crown |
0.5 |
max_diameter |
float |
maximum diameter allowed for single crown growth |
8 |
Returns:
| Type | Description |
|---|---|
ndarray |
2-d array with each pixel labeled with a unique crown ID. |
Source code in lazer/structure.py
def crown_seg_dalponte_coomes(
height: np.ndarray,
tree_top_x: np.ndarray,
tree_top_y: np.ndarray,
min_inclusion_height: float = defaults.min_istree_height,
height_decrease_rate: float = defaults.crown_height_decrease_rate,
height_grouping_rate: float = defaults.crown_height_grouping_rate,
max_diameter: float = defaults.max_crown_diameter,
) -> np.ndarray:
"""Delineates crown vectors using the Dalponte/Coomes approach.
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12575
Method described in the Supporting Information Appendix S2.
This routine calls a `numba` function to run the operation, which
doesn't appear to allow specifying default keyword arguments.
Args:
height: 2d array with tree height
tree_top_x: x index locations for tree tops
tree_top_y: y index locations for tree tops
min_inclusion_height: minimum height threshold for inclusion in crown
height_decrease_rate: relative pixel-to-pixel decrease allowed for
inclusion in a neighboring crown
height_grouping_rate: relative pixel-to-pixel increase allowed for
inclusion in a neighboring crown
max_diameter: maximum diameter allowed for single crown growth
Returns:
2-d array with each pixel labeled with a unique crown ID.
"""
return _crown_seg_dalponte_coomes(
height,
tree_top_x,
tree_top_y,
min_inclusion_height,
height_decrease_rate,
height_grouping_rate,
max_diameter,
)
crowns_to_vector(crowns, height, vector, transform=None, crs=None, epsg=None, wkt=None, driver='ESRI Shapefile', nodata=0)
¶
Converts labeled tree crown arrays to a polygon vector file.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
crowns |
ndarray |
2-d array with unique crown ID labels for each pixel value. typically generated by one of the crown_seg_() algorithms. |
required |
height |
ndarray |
2-d array of tree height values used to create the crown data. |
required |
vector |
ndarray |
path to output vector file. |
required |
transform |
Affine |
the rasterio affine transform for applying georeferencing. you typically get this via src.transform when opening a file with rasterio. the format is [xres, xoffset, xmin, yoffset, yres, ymax]. uses pixel coordinates if this isn't passed. |
None |
crs |
CRS |
the coordinate reference system for these data. A crs object is created if you instead pass epsg or wkt arguments. |
None |
epsg |
int |
the raster EPSG code. Use instead of passing crs or wkt. |
None |
wkt |
str |
the raster wkt string. Use instead of passing crs or epsg. |
None |
driver |
str |
the output vector format driver. available choices can be read with lazer.geo.get_fiona_drivers(). |
'ESRI Shapefile' |
nodata |
float |
the background value to ignore in labeling unique regions. |
0 |
Returns:
| Type | Description |
|---|---|
None |
None. Point-format data are written to the vector file. |
Source code in lazer/structure.py
def crowns_to_vector(
crowns: np.ndarray,
height: np.ndarray,
vector: np.ndarray,
transform: rio.transform.Affine = None,
crs: pyproj.CRS = None,
epsg: int = None,
wkt: str = None,
driver: str = defaults.vector_format,
nodata: float = 0,
) -> None:
"""Converts labeled tree crown arrays to a polygon vector file.
Args:
crowns: 2-d array with unique crown ID labels for each pixel value.
typically generated by one of the crown_seg_() algorithms.
height: 2-d array of tree height values used to create the crown data.
vector: path to output vector file.
transform: the rasterio affine transform for applying georeferencing.
you typically get this via src.transform when opening a file with rasterio.
the format is [xres, xoffset, xmin, yoffset, yres, ymax].
uses pixel coordinates if this isn't passed.
crs: the coordinate reference system for these data.
A crs object is created if you instead pass epsg or wkt arguments.
epsg: the raster EPSG code. Use instead of passing crs or wkt.
wkt: the raster wkt string. Use instead of passing crs or epsg.
driver: the output vector format driver.
available choices can be read with lazer.geo.get_fiona_drivers().
nodata: the background value to ignore in labeling unique regions.
Returns:
None. Point-format data are written to the vector file.
"""
if not crs:
if epsg:
crs = pyproj.CRS.from_epsg(epsg)
elif wkt:
crs = pyproj.CRS.from_wkt(wkt)
# get the geojson coordinates for each labeled region
crown_shapes = raster_to_shapes(
crowns, mask=(crowns != nodata), transform=transform
)
# output vector format definition
crown_schema = {
"geometry": "Polygon",
"properties": OrderedDict(
[
("ht_mean", "float"),
("ht_max", "float"),
("area", "float"),
("circ", "float"),
]
),
}
# write to the output file feature-by-feature
with fio.open(vector, "w", driver=driver, crs=crs, schema=crown_schema) as dst:
features = []
for polygon, value in crown_shapes:
crown_height = height[crowns == value]
height_mean = float(crown_height.mean())
height_max = float(crown_height.max())
crown_polygon = Polygon(polygon["coordinates"][0])
area = crown_polygon.area
circ = crown_polygon.length
feature = {
"geometry": polygon,
"properties": {
"ht_mean": height_mean,
"ht_max": height_max,
"area": area,
"circ": circ,
},
}
features.append(feature)
dst.writerecords(features)
get_circular_kernel(radius=3)
¶
Returns a circular filtering kernel.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
radius |
int |
kernel radius, which defines the output filter size. |
3 |
Returns:
| Type | Description |
|---|---|
ndarray |
a 2d ndarray of shape (radius2, radius2). |
Source code in lazer/structure.py
def get_circular_kernel(radius: int = defaults.crown_search_radius) -> np.ndarray:
"""Returns a circular filtering kernel.
Args:
radius: kernel radius, which defines the output filter size.
Returns:
a 2d ndarray of shape (radius*2, radius*2).
"""
y, x = np.ogrid[-radius : radius + 1, -radius : radius + 1]
kernel = x ** 2 + y ** 2 <= radius ** 2
return kernel
identify_tree_tops(height, crown_radius=3, min_height=5, radius_units='pixels', resolution=1.0)
¶
Identifies tree top locations from a 2d array of tree height values.
Uses an image-based maximum filter, meaning the seach radius for what constitutes a tree must resolve to integer units. This is of particular importance if you use radius_units='meters'.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
height |
ndarray |
2d array of shape (ysize, xsize) with height (in meters) for each element |
required |
crown_radius |
int |
the maximum radius of an individual tree. determines how many unique tree tops will be identified. |
3 |
min_height |
float |
the minimum height for a tree top to be labeled as a tree. |
5 |
radius_units |
str |
whether to compute radius by 'pixels' or by 'meters'. must also set the resolution keyword if radius_units='meters'. |
'pixels' |
resolution |
float |
the grid size for the height array. |
1.0 |
Returns:
| Type | Description |
|---|---|
(x, y, tree_height) |
1-d numpy arrays with the tree top locations in array coordinates and the height of the trees at each location. |
Source code in lazer/structure.py
def identify_tree_tops(
height: np.ndarray,
crown_radius: int = defaults.crown_search_radius,
min_height: float = defaults.min_treetop_height,
radius_units: str = "pixels",
resolution: float = defaults.resolution,
) -> _Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Identifies tree top locations from a 2d array of tree height values.
Uses an image-based maximum filter, meaning the seach radius for what constitutes
a tree must resolve to integer units. This is of particular importance if you use
radius_units='meters'.
Args:
height: 2d array of shape (ysize, xsize) with height (in meters) for each element
crown_radius: the maximum radius of an individual tree.
determines how many unique tree tops will be identified.
min_height: the minimum height for a tree top to be labeled as a tree.
radius_units: whether to compute radius by 'pixels' or by 'meters'.
must also set the resolution keyword if radius_units='meters'.
resolution: the grid size for the height array.
Returns:
(x, y, tree_height): 1-d numpy arrays with the tree top locations in array
coordinates and the height of the trees at each location.
"""
if radius_units == "meters":
crown_radius = round(crown_radius / resolution)
max_kernel = get_circular_kernel(crown_radius)
max_filter = ndimage.filters.maximum_filter(height, footprint=max_kernel)
max_height = height == max_filter
max_height[height < min_height] = 0
tree_tops, tree_counts = ndimage.label(max_height)
if tree_counts == 0:
return None, None, None
elif tree_counts == 1:
tree_ids = 1
else:
tree_ids = range(1, tree_counts + 1)
yx = (
np.array(
ndimage.center_of_mass(height, tree_tops, tree_ids),
dtype=np.float32,
)
+ 0.5
)
x, y = yx[:, 1], yx[:, 0]
tree_height = height[y.astype(int), x.astype(int)]
return x, y, tree_height
tree_tops_to_vector(tree_tops, vector, bounds, resolution, crs=None, epsg=None, wkt=None, driver='ESRI Shapefile')
¶
Converts tree top locations from array indices to a point-format vector file.
Many of the parameters required are used to translate from pixel coordinates to geographic coordinates.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
tree_tops |
Tuple[numpy.ndarray, numpy.ndarray, <built-in function array>] |
(x, y, tree_height) tuple of locations/heights from |
required |
vector |
str |
path to output vector file. |
required |
bounds |
Tuple[float, float, float, float] |
tuple of [xmin, ymin, xmax, ymax] coordinates for the input raster. |
required |
resolution |
Tuple[float, float] |
tuple of [xres, yres] for the input spatial resolution. |
required |
crs |
CRS |
the coordinate reference system for these data. A crs object is created if you instead pass epsg or wkt arguments. |
None |
epsg |
int |
the raster EPSG code. Use instead of passing crs or wkt. |
None |
wkt |
str |
the raster wkt string. Use instead of passing crs or epsg. |
None |
driver |
str |
the output vector format driver. available choices can be read with lazer.geo.get_fiona_drivers(). |
'ESRI Shapefile' |
Returns:
| Type | Description |
|---|---|
None |
None. Point-format data are written to the vector file. |
Source code in lazer/structure.py
def tree_tops_to_vector(
tree_tops: _Tuple[np.ndarray, np.ndarray, np.array],
vector: str,
bounds: _Tuple[float, float, float, float],
resolution: _Tuple[float, float],
crs: pyproj.CRS = None,
epsg: int = None,
wkt: str = None,
driver: str = defaults.vector_format,
) -> None:
"""Converts tree top locations from array indices to a point-format vector file.
Many of the parameters required are used to translate from pixel coordinates
to geographic coordinates.
Args:
tree_tops: (x, y, tree_height) tuple of locations/heights from `identify_tree_tops()`.
vector: path to output vector file.
bounds: tuple of [xmin, ymin, xmax, ymax] coordinates for the input raster.
resolution: tuple of [xres, yres] for the input spatial resolution.
crs: the coordinate reference system for these data.
A crs object is created if you instead pass epsg or wkt arguments.
epsg: the raster EPSG code. Use instead of passing crs or wkt.
wkt: the raster wkt string. Use instead of passing crs or epsg.
driver: the output vector format driver.
available choices can be read with lazer.geo.get_fiona_drivers().
Returns:
None. Point-format data are written to the vector file.
"""
if not crs:
if epsg:
crs = pyproj.CRS.from_epsg(epsg)
elif wkt:
crs = pyproj.CRS.from_wkt(wkt)
else:
raise KeyError("One of crs, epsg, or wkt must be set")
xp, yp, tree_heights = tree_tops
xmin, ymin, xmax, ymax = bounds
xres, yres = resolution
xs = (xp * xres) + xmin
ys = (yp * -yres) + ymax
# output vector format definition
tree_schema = {
"geometry": "Point",
"properties": OrderedDict([("height", "float")]),
}
with fio.open(vector, "w", driver=driver, crs=crs, schema=tree_schema) as dst:
features = []
for x, y, height in zip(xs, ys, tree_heights):
feature = {
"geometry": {"type": "Point", "coordinates": (x, y)},
"properties": {"height": float(height)},
}
features.append(feature)
dst.writerecords(features)