sylva.canopy¶

sylva.canopy ¶

Canopy structure: voxels, plant area density, gap fraction, LAI.

VoxelGrid `dataclass` ¶

Point counts on a regular 3D grid; build with voxelize.

Attributes:

Name	Type	Description
`counts`	`ndarray`	`(nz, ny, nx)` integer counts; `counts[k, j, i]` indexes z, y, x.
`origin`	`ndarray`	Minimum corner of the grid, `(x, y, z)`.
`voxel_size`	`float`	Voxel edge length (m).

Notes

Occupancy is not density: an empty voxel may simply not have been seen. For plant area density use the ray-traced sylva.voxels.ray_voxelize.

occupied `property` ¶

occupied: ndarray

Boolean (nz, ny, nx), True where a voxel holds at least one point.

shape `property` ¶

shape: tuple[int, int, int]

Grid size as (nx, ny, nz) (the reverse of counts.shape).

z_levels ¶

z_levels() -> np.ndarray

Bottom z of each voxel layer.

Returns:

Type	Description
`ndarray`	Length `nz`, bottom first (m).

occupied_centers ¶

occupied_centers() -> np.ndarray

Centres of the occupied voxels, e.g. for plotting.

Returns:

Type	Description
`ndarray`	`(n, 3)` x, y, z of each voxel with at least one point.

vertical_profile ¶

vertical_profile() -> np.ndarray

Fraction of voxels occupied in each layer.

Returns:

Type	Description
`ndarray`	Length `nz`, bottom first, values 0 to 1.

DensityGrid `dataclass` ¶

Ray-traced voxel statistics (raycloudtools / Lowe et al. 2021 style).

Arrays are [k, j, i] (z, y, x). density is 2 (n-1)/n · hits / path_length per voxel (spherical leaf angles), NaN where fewer than min_hits hits. Build with density_grid.

For most work sylva.voxels.ray_voxelize is the fuller model (beam width, several estimators, occlusion); this grid is the simple, fast one.

Attributes:

Name	Type	Description
`n_rays`	`ndarray`	Pulses entering each voxel.
`n_hits`	`ndarray`	Echoes inside each voxel.
`path_length`	`ndarray`	Summed beam path length through each voxel (m).
`density`	`ndarray`	Plant area density (m² m⁻³).
`profile`	`ndarray`	Mean density per layer, bottom first.
`origin`	`ndarray`	Minimum corner `(x, y, z)`.
`voxel_size`	`float`	Voxel edge (m).

pai `property` ¶

pai: float

Plant area index: the layer means summed over height (m² m⁻²).

Includes ground layers unless mask_ground was applied.

z_levels ¶

z_levels() -> np.ndarray

Bottom z of each voxel layer.

Returns:

Type	Description
`ndarray`	Length `nz`, bottom first (m).

centers ¶

centers() -> tuple[np.ndarray, np.ndarray, np.ndarray]

Voxel-centre coordinates.

Returns:

Type	Description
`X, Y, Z : numpy.ndarray`	Each shaped like `density`, `(nz, ny, nx)`.

height_above ¶

height_above(dtm: Raster) -> np.ndarray

Height of each voxel centre above the terrain.

Parameters:

Name	Type	Description	Default
`dtm`	`Raster`	Terrain in the grid's frame.	required

Returns:

Type	Description
`ndarray`	Shaped like `density` (m).

mask_ground ¶

mask_ground(dtm: Raster, margin: float | None = None) -> DensityGrid

Remove ground voxels so terrain returns do not count as plant material.

Parameters:

Name	Type	Description	Default
`dtm`	`Raster`	Terrain in the grid's frame.	required
`margin`	`float \| None`	Voxels whose centre is less than this height (m) above the terrain are set to NaN; one voxel if None.	`None`

Returns:

Type	Description
`DensityGrid`	A copy with `density` masked and `profile` recomputed.

profile_above_ground ¶

profile_above_ground(dtm: Raster, bin_size: float | None = None, max_height: float | None = None, margin: float | None = None, pooled: bool = True) -> tuple[np.ndarray, np.ndarray]

Plant area density by height above the terrain.

Bins voxel centres by their height above dtm (default bin = voxel size), ignoring voxels whose centre is below margin (default one voxel, which drops the ground-containing layer). With pooled=True each bin is treated as one big voxel — 2 Σhits / Σpath over all its rays — which weights voxels by how well they were sampled and is robust to the sparsely-sampled, occluded voxels of a single scan. pooled=False averages the per-voxel density instead.

Parameters:

Name	Type	Description	Default
`dtm`	`Raster`	Terrain in the grid's frame.	required
`bin_size`	`float \| None`	Height bin (m); the voxel size if None.	`None`
`max_height`	`float \| None`	Top of the profile; the highest sampled voxel if None.	`None`
`margin`	`float \| None`	Minimum voxel-centre height (m); one voxel if None.	`None`
`pooled`	`bool`	Pool hits and path lengths per bin (recommended).	`True`

Returns:

Type	Description
`bin_bottoms, pad : numpy.ndarray`	PAD in m² m⁻³; `np.nansum(pad) * bin_size` is PAI.

GapProfile `dataclass` ¶

Gap probability by zenith ring, azimuth sector and height above ground, pooled over scans (Jupp et al. 2009), and what follows from it.

Build it with empty and add_scan (one call per scan position), then read report. Returns are weighted 1 / n per echo of an n-echo pulse (equal weighting, Armston et al. 2013); shots are pulses fired.

Validated against pylidar on TERN plots (within 4 % where unsaturated) and against hemispherical photographs; see Benchmarks > Canopy gap profiles and the Pulse data guide for the full workflow.

Examples:

>>> prof = GapProfile.empty()
>>> for pos in project.with_scans():
...     s = io.read_rxp_shots(pos.rxp, shot_stride=4)
...     fired = fired_pulses_per_ring(s, pos.pattern, prof.zenith_edges, shot_stride=4)
...     s = s.transform(pos.sop)
...     xyz = s.echo_xyz()
...     a, b, c = fit_ground_plane(xyz, centre=pos.sop[:2, 3], radius=25)
...     prof.add_scan(s, xyz[:, 2] - (a * xyz[:, 0] + b * xyz[:, 1] + c), fired_per_ring=fired)
>>> prof.report()["pai_hinge"]

heights `property` ¶

heights: ndarray

Top of each height bin (m): the profile at z counts returns below it.

zenith `property` ¶

zenith: ndarray

Ring centres (deg).

empty `classmethod` ¶

empty(zenith_edges=np.arange(5.0, 75.0, 5.0), n_azimuth: int = 36, height_bin: float = 0.5, max_height: float = 80.0) -> 'GapProfile'

Start an empty profile.

Parameters:

Name	Type	Description	Default
`zenith_edges`		Ring edges (degrees from up). The default, 5-75 degrees in 5 degree rings, avoids the near-horizontal rings that see mostly stems.	`arange(5.0, 75.0, 5.0)`
`n_azimuth`	`int`	Azimuth sectors per ring, used for the clumping index.	`36`
`height_bin`	`float`	Height resolution (m).	`0.5`
`max_height`	`float`	Top of the profile (m); returns above it go in the top bin.	`80.0`

Returns:

Type	Description
`GapProfile`

add_scan ¶

add_scan(shots: Shots, echo_heights, fired_per_ring=None, min_height: float = -np.inf) -> None

Add one scan position's pulses; echo_heights above ground, one per echo. Without fired_per_ring the shots must include the pulses that returned nothing (e.g. Shots.fill_missing or a ray cloud).

Returns below min_height are not counted. By default all count, and those below zero go in the lowest bin: a pulse fired upwards cannot hit the ground, so a negative height means the ground model is off there, not that the return is not vegetation.

Parameters:

Name	Type	Description	Default
`shots`	`Shots`	One scan position's pulses. Only their directions and echo counts are used, so the scanner or project frame both work.	required
`echo_heights`		Height above ground of each echo (length `shots.n_echoes`), e.g. from a DTM or `fit_ground_plane`.	required
`fired_per_ring`		Pulses fired into each ring (`fired_pulses_per_ring`) for RiVLib streams, which lack the misses.	`None`
`min_height`	`float`	Echoes below this height (m) are dropped entirely.	`-inf`

Notes

The profile is updated in place; add every scan of the plot before reading it.

pgap ¶

pgap() -> np.ndarray

Gap probability by ring and height.

Returns:

Type	Description
`ndarray`	`(rings, heights)`: 1 minus the returns between `min_height` and the top of each height bin, over the pulses fired, clipped to 0-1. NaN for rings without pulses.

pai_profile ¶

pai_profile(method: str = 'hinge') -> np.ndarray

Cumulative plant area index below each height (effective, not corrected for clumping).

hinge: -1.1 ln P(57.5 deg) from the ring holding 57.5 deg. linear: Jupp et al. (2009): -ln P(theta) = PAI_h + (2 / pi) tan(theta) PAI_v (horizontal and vertical foliage) fitted over the rings, PAI = PAI_h + PAI_v; the mean leaf angle is atan(PAI_v / PAI_h). weighted: 2 cos(theta) (-ln P) averaged over rings with weights sin(theta) (Miller 1967 over the rings measured, spherical leaves).

Parameters:

Name	Type	Description	Default
`method`	`('hinge', 'linear', 'weighted')`	Estimator. Hinge is the most robust; linear also gives leaf angle.	`"hinge"`

Returns:

Type	Description
`ndarray`	PAI below the top of each height bin (m² m⁻²), non-decreasing.

Raises:

Type	Description
`ValueError`	For an unknown method, or `"hinge"` with no pulses in the ring holding 57.5 degrees.

pavd_profile ¶

pavd_profile(method: str = 'hinge') -> np.ndarray

Plant area volume density: the height derivative of pai_profile.

Parameters:

Name	Type	Description	Default
`method`	`('hinge', 'linear', 'weighted')`	Estimator, as for `pai_profile`.	`"hinge"`

Returns:

Type	Description
`ndarray`	PAVD per height bin (m² m⁻³).

clumping ¶

clumping(zenith: float = 57.5) -> float

Lang & Xiang (1986) clumping index at the ring holding zenith: ln(mean P) / mean(ln P) over every (scan, azimuth sector) segment. 1 is random foliage; below 1 clumped.

Parameters:

Name	Type	Description	Default
`zenith`	`float`	Zenith angle (degrees) whose ring is used.	`57.5`

Returns:

Type	Description
`float`	Clumping index; NaN without scans or without any gap. Divide effective PAI by it for true PAI.

report ¶

report(top_fraction: float = 0.99, saturation_gap: float = 0.005) -> dict

Plot summary: effective PAI (hinge, linear, weighted), mean leaf angle from the linear fit, clumping index and clumping-corrected hinge PAI, canopy height (where the hinge PAI profile reaches top_fraction of its total), canopy closure at 57.5 deg, cover at the steepest ring measured, and the profiles.

saturated is set when less than saturation_gap of the pulses at 57.5 deg got through the canopy (hinge PAI above about 5.8): the PAI is then bounded by the pulse count, not measured. Dense rainforest seen from the ground reaches it.

Parameters:

Name	Type	Description	Default
`top_fraction`	`float`	Fraction of total PAI that defines canopy height.	`0.99`
`saturation_gap`	`float`	Gap fraction at 57.5 degrees below which `saturated` is set.	`0.005`

Returns:

Type	Description
`dict`	Scalars `saturated`, `gap_57`, `pai_hinge`, `pai_linear`, `pai_weighted`, `mla_linear` (degrees), `clumping`, `pai_hinge_corrected`, `canopy_height` (m), `closure_57`, `cover`, `cover_zenith`, `n_scans`, `pulses`; and arrays `height`, `pai_hinge_profile`, `pavd_hinge`, `pai_linear_profile`, `pavd_linear`. Check `saturated` before reporting PAI.

voxelize ¶

voxelize(cloud: PointCloud, voxel_size: float, origin=None, shape=None) -> VoxelGrid

Count points per voxel.

Parameters:

Name	Type	Description	Default
`cloud`	`PointCloud`	Input points.	required
`voxel_size`	`float`	Voxel edge length (m).	required
`origin`		Minimum corner `(x, y, z)`; the cloud's minimum if None. Fix it (with `shape`) to compare grids across dates or plots.	`None`
`shape`		`(nx, ny, nz)`; enough to cover the cloud if None. Points outside the grid are ignored.	`None`

Returns:

Type	Description
`VoxelGrid`

vertical_profile ¶

vertical_profile(cloud: PointCloud, bin_size: float = 0.5, height_attr: str = 'height', max_height: float | None = None) -> tuple[np.ndarray, np.ndarray]

Point counts by height above ground.

Parameters:

Name	Type	Description	Default
`cloud`	`PointCloud`	A height-normalised cloud.	required
`bin_size`	`float`	Bin height (m).	`0.5`
`height_attr`	`str`	Attribute holding heights; z is used if absent.	`'height'`
`max_height`	`float \| None`	Top of the last bin; the highest point if None.	`None`

Returns:

Type	Description
`bin_bottoms, counts : numpy.ndarray`	Point counts depend on scanner distance and occlusion, so this is a description of the data, not of the canopy.

pad_profile_voxel ¶

pad_profile_voxel(cloud: PointCloud, voxel_size: float = 0.5, height_attr: str = 'height', max_height: float | None = None, clumping: float = 1.0) -> tuple[np.ndarray, np.ndarray]

Plant area density profile by the vertical contact-frequency method (simplified Hosoi & Omasa 2006, vertical beams, G = 0.5): PAD_i = -ln(1 - N_i) / (G dz), where N_i is the fraction of voxel columns not yet intercepted above layer i that are occupied in it.

Quick and needs only points, but it ignores beam direction and occlusion. For a path-length estimate use sylva.voxels.ray_voxelize or GapProfile with pulse data.

Parameters:

Name	Type	Description	Default
`cloud`	`PointCloud`	A height-normalised cloud.	required
`voxel_size`	`float`	Voxel edge and layer thickness (m); results depend strongly on it.	`0.5`
`height_attr`	`str`	Attribute holding heights.	`'height'`
`max_height`	`float \| None`	Top of the profile; the highest point if None.	`None`
`clumping`	`float`	Factor the PAD is multiplied by (1 = none).	`1.0`

Returns:

Type	Description
`layer_bottoms, pad : numpy.ndarray`	Heights (m) and plant area density (m² m⁻³).

gap_fraction_zenith ¶

gap_fraction_zenith(shots: Shots, echo_heights: ndarray, min_height: float = 0.0, zenith_edges: ndarray | None = None) -> tuple[np.ndarray, np.ndarray]

Directional gap fraction by zenith ring from pulses that include misses.

A pulse is a gap when none of its echoes is above min_height. Needs every fired pulse, so use shots from a ray cloud or after sylva.Shots.fill_missing; for raw RiVLib streams use gap_fraction_pattern or GapProfile.

Parameters:

Name	Type	Description	Default
`shots`	`Shots`	Pulses of one scan position.	required
`echo_heights`	`ndarray`	Height above ground of each echo (length `shots.n_echoes`).	required
`min_height`	`float`	Echoes at or below this height (m) do not block the pulse.	`0.0`
`zenith_edges`	`ndarray \| None`	Ring edges (degrees from up); 0-90 in 5 degree rings if None.	`None`

Returns:

Type	Description
`centres_deg, gap : numpy.ndarray`	Ring centres and gap fraction; NaN for rings without pulses.

gap_fraction_pattern ¶

gap_fraction_pattern(shots: Shots, echo_heights: ndarray, pattern: dict, min_height: float = 0.0, zenith_edges: ndarray | None = None, pulses_per_line: int | None = None) -> tuple[np.ndarray, np.ndarray]

Gap fraction by zenith ring for a RIEGL scan whose no-return pulses are absent from the stream (RiVLib). The number of pulses fired per ring comes from the angular scan pattern (see Shots.expected_per_zenith); a fired pulse with no echo above min_height is a gap. Pass the shots in the scanner frame, or pulses_per_line from the scanner-frame shots.

Parameters:

Name	Type	Description	Default
`shots`	`Shots`	Pulses of one scan position.	required
`echo_heights`	`ndarray`	Height above ground of each echo.	required
`pattern`	`dict`	Scan pattern from `sylva.riscan.ScanPosition.pattern`.	required
`min_height`	`float`	Echoes at or below this height (m) do not block the pulse.	`0.0`
`zenith_edges`	`ndarray \| None`	Ring edges (degrees); 0-90 in 5 degree rings if None.	`None`
`pulses_per_line`	`int \| None`	Pulses fired per zenith line; estimated from `shots` if None.	`None`

Returns:

Type	Description
`centres_deg, gap : numpy.ndarray`	Rings outside the scanned zenith range are NaN.

lai_from_gap_fraction ¶

lai_from_gap_fraction(zenith_deg, gap_fraction, method: str = 'hinge') -> float

Effective PAI from directional gap fraction.

"hinge": -ln P(57.5°) cos(57.5°) / 0.5 (Wilson 1963). "miller": 2 ∫ -ln P(θ) cos θ sin θ dθ (Miller 1967).

NaN rings are skipped; a gap fraction of 0 (saturated ring, common in dense forest) is floored at 1e-5, so the result is large but finite.

Parameters:

Name	Type	Description	Default
`zenith_deg`		Ring centres (degrees).	required
`gap_fraction`		Gap fraction per ring.	required
`method`	`('hinge', 'miller')`	`"hinge"` uses the ring nearest 57.5 degrees, where the result is almost independent of leaf angle; `"miller"` needs rings covering 0-90 degrees.	`"hinge"`

Returns:

Type	Description
`float`	Effective PAI (m² m⁻²), not corrected for clumping.

Raises:

Type	Description
`ValueError`	For an unknown method or no valid rings.

canopy_cover ¶

canopy_cover(chm_data: ndarray, threshold: float = 2.0) -> float

Canopy cover from a canopy height model.

Parameters:

Name	Type	Description	Default
`chm_data`	`ndarray`	CHM values, e.g. `make_chm(...).data`; NaN cells are ignored.	required
`threshold`	`float`	Height (m) at or above which a cell counts as canopy.	`2.0`

Returns:

Type	Description
`float`	Fraction of valid cells that are canopy; NaN if none are valid.

density_grid ¶

density_grid(shots: Shots, voxel_size: float, origin=None, shape=None, min_hits: int = 2) -> DensityGrid

Ray-trace pulses through a voxel grid and estimate plant area density.

Pulses without echoes run to the grid edge.

Parameters:

Name	Type	Description	Default
`shots`	`Shots`	Pulses from one or more scan positions in a common frame.	required
`voxel_size`	`float`	Voxel edge (m).	required
`origin`		Grid minimum corner and `(nx, ny, nz)`; both None covers the echoes.	`None`
`shape`		Grid minimum corner and `(nx, ny, nz)`; both None covers the echoes.	`None`
`min_hits`	`int`	Voxels with fewer echoes get NaN density.	`2`

Returns:

Type	Description
`DensityGrid`

fit_ground_plane ¶

fit_ground_plane(points, cell: float = 1.0, centre=None, radius: float | None = None, iterations: int = 20) -> np.ndarray

Ground plane z = a x + b y + c through the lowest point of every cell (m) grid cell, optionally within radius of centre (xy), fitted with Huber-weighted least squares so tree bases and pits pull little (after Calders et al. 2014).

Used for the single-scan gap profiles, where a plane around the scanner is enough and a DTM may not exist.

Parameters:

Name	Type	Description	Default
`points`		A `~sylva.PointCloud` or `(N, 3)` array; for a scan, downward echoes work best.	required
`cell`	`float`	Grid cell size (m) for the lowest-point selection.	`1.0`
`centre`		Only use points within `radius` (m) of `centre` (xy).	`None`
`radius`		Only use points within `radius` (m) of `centre` (xy).	`None`
`iterations`	`int`	Reweighting iterations.	`20`

Returns:

Type	Description
`ndarray`	`[a, b, c]`; ground height at (x, y) is `a * x + b * y + c`.

Raises:

Type	Description
`ValueError`	With fewer than 3 points.

fired_pulses_per_ring ¶

fired_pulses_per_ring(shots_scanner: Shots, pattern: dict, zenith_edges, shot_stride: int = 1, ground_zenith=(100.0, 125.0)) -> np.ndarray

Pulses a RIEGL scan fired into each zenith ring, for streams without the pulses that returned nothing (RiVLib). shots_scanner in the scanner frame, read with shot_stride.

Every azimuth step fires one pulse per zenith line, so a ring gets (lines in the ring) x (pulses per line). Pulses per line are counted on the downward lines in ground_zenith, where every pulse hits the ground and so every fired pulse is in the stream: the median there. The nominal phi_count of the pattern is not used, because the scanner fires about 1 % more pulses than nominal and a percentile over all lines overshoots, since the mirror's zenith angles do not sit exactly on the nominal lines. In dense canopy, where nearly every pulse returns, a few per cent too many fired pulses read as gaps and cap the PAI.

Parameters:

Name	Type	Description	Default
`shots_scanner`	`Shots`	Pulses of one scan in the scanner frame.	required
`pattern`	`dict`	Scan pattern from `sylva.riscan.ScanPosition.pattern`.	required
`zenith_edges`		Ring edges (degrees).	required
`shot_stride`	`int`	The `shot_stride` used when reading (for the fallback estimate).	`1`
`ground_zenith`		Zenith range (degrees) of the lines used to count pulses per line. If fewer than 10 lines fall in it, the estimate from `sylva.Shots.pulses_per_line` is used instead.	`(100.0, 125.0)`

Returns:

Type	Description
`ndarray`	Pulses fired per ring (float, length `len(zenith_edges) - 1`), for `GapProfile.add_scan`.

fired_pulses_from_points ¶

fired_pulses_from_points(shots_scanner: Shots, zenith_edges, ground_zenith=(100.0, 125.0), limit_quantile: float = 1e-05, field_of_view: float | None = 100.0) -> np.ndarray

Pulses a scan fired into each zenith ring, from its returns alone.

For a scan known only by its points (a LAS/LAZ export and its SOP), with no .rxp or scan pattern. A RIEGL scanner fires at a constant angular step down each mirror sweep and sweeps every azimuth the same number of times, so pulses fired per degree of zenith are the same at every zenith. On the downward lines in ground_zenith nearly every pulse hits the ground (0.3-0.8 % did not on a VZ-2000i scan in savanna), so the returns there, per degree, give that rate; a ring gets the rate times the part of its width inside the scan's zenith limits. The lower limit is read from the returns, since every downward pulse hits the ground; the upper one is field_of_view degrees above it, because in open vegetation the most upward pulses return nothing.

Parameters:

Name	Type	Description	Default
`shots_scanner`	`Shots`	Pulses of one scan in the scanner frame, e.g. `Shots.from_pointcloud(cloud.transform(inv(sop)))`. Echoes of a pulse must share `gps_time`, or the pulse counts more than once.	required
`zenith_edges`		Ring edges (degrees).	required
`ground_zenith`		Zenith range (degrees) where every pulse is taken to return. It must lie inside the scan's field of view and look at open ground.	`(100.0, 125.0)`
`limit_quantile`	`float`	The lower zenith limit is this upper quantile of the returns' zenith.	`1e-05`
`field_of_view`	`float \| None`	The scanner's vertical field of view (degrees): 100 for RIEGL's VZ scanners. None reads the upper limit from the returns too, which is right only where the canopy returns the most upward pulses.	`100.0`

Returns:

Type	Description
`ndarray`	Pulses fired per ring (float, length `len(zenith_edges) - 1`), for `GapProfile.add_scan`.

Raises:

Type	Description
`ValueError`	If `ground_zenith` lies outside the scan's zenith limits.

sylva.canopy¶

sylva.canopy ¶

VoxelGrid dataclass ¶

occupied property ¶

shape property ¶

z_levels ¶

occupied_centers ¶

vertical_profile ¶

DensityGrid dataclass ¶

pai property ¶

z_levels ¶

centers ¶

height_above ¶

mask_ground ¶

profile_above_ground ¶

GapProfile dataclass ¶

heights property ¶

zenith property ¶

empty classmethod ¶

add_scan ¶

pgap ¶

pai_profile ¶

pavd_profile ¶

clumping ¶

report ¶

voxelize ¶

vertical_profile ¶

pad_profile_voxel ¶

gap_fraction_zenith ¶

gap_fraction_pattern ¶

lai_from_gap_fraction ¶

canopy_cover ¶

density_grid ¶

fit_ground_plane ¶

fired_pulses_per_ring ¶

fired_pulses_from_points ¶

VoxelGrid `dataclass` ¶

occupied `property` ¶

shape `property` ¶

DensityGrid `dataclass` ¶

pai `property` ¶

GapProfile `dataclass` ¶

heights `property` ¶

zenith `property` ¶

empty `classmethod` ¶