sylva.voxels¶

sylva.voxels ¶

Ray-traced voxel grids with Beer-Lambert statistics.

A port of rayvoxel (J. Rivory, unpublished), built on raycloudtools. Every pulse is traced through the grid, and each voxel accumulates beam counts, potential and free path lengths, beam sections and mean beam angles. From these the attenuation coefficient is estimated by free path length (fpl, bias corrected after Pimont et al. 2018), exact potential path length (ppl), transmittance or Bailey & Mahaffee's (2017) eq. 10 (bailey), and divided by a projection function G to give plant, leaf and wood area density. G comes from an analytic leaf angle distribution or, with inclination=True, from per-tree inclination angle distributions estimated on the echoes (Vicari et al. 2019).

sylva.canopy.density_grid is the lightweight alternative when only hits and path lengths are needed.

RayVoxelGrid ¶

Voxel statistics from ray_voxelize.

The grid lives in the Rust core; arrays are copied out on first use and cached, shaped (nz, ny, nx) like sylva.canopy.DensityGrid. Raw accumulators (fields) and derived quantities (metrics, plus attenuation_<m>, pad_<m>, lad_<m>, wad_<m> for any attenuation method m) are available by name or as attributes::

grid["num_hits"]; grid.path_length; grid.pad_fpl; grid.transmittance

Notes

Main raw fields (sums over the pulses crossing each voxel):

Field	Meaning
`num_beams`	pulses entering the voxel
`num_hits`	echoes in the voxel (`num_hit_leaf`, `_wood`, `_plant` by class)
`num_hits_weighted`	echoes weighted by their share of the pulse (`weighting`)
`num_beams_occluded`	pulses reaching it only after their last echo (`occlusion=True`)
`path_length`	potential path length: full chords of the entering pulses (m)
`free_path_length`	path actually travelled inside the voxel (m)
`bs_entering`, `bs_intercepted`	beam section entering and stopped (m², needs `beam`)
`wood_volume`	QSM wood (m³) after `add_wood_volume`

Derived metrics:

Metric	Meaning
`state`	0 unobserved, 1 occluded, 2 empty, 3 filled (`STATES`)
`attenuation_<m>`	attenuation λ (m⁻¹) by method `m` (see `attenuation`)
`pad_<m>`, `lad_<m>`, `wad_<m>`	plant, leaf, wood area density λ / G (m² m⁻³)
`pad_g0_5`	plant area density with G = 0.5
`pad_g_corrected`	PAD with G at the voxel's mean beam zenith, first attenuation method
`surface_area`	`pad_g0_5` × voxel volume (m²)
`transmittance`	beam-section transmittance (0-1)
`mean_zenith_angle`, `mean_azimuth_angle`	mean beam direction (degrees)
`azimuth_concentration`	0 (all azimuths) to 1 (one azimuth)
`mean_laser_dist`	mean distance from the scanner (m)
`sd_path_length`	spread of the chord lengths (m)
`distance_from_ground`	voxel centre above the DTM (m, needs `dtm`)
`exploration_rate`	share of sub-voxels crossed (`subvoxel_split`)
`g_plant`, `g_leaf`, `g_wood`	G used per voxel
`wood_volume_density`	QSM wood volume per voxel volume (m³ m⁻³)

origin `property` ¶

origin: ndarray

Minimum corner of the grid, (x, y, z).

voxel_size `property` ¶

voxel_size: float

Voxel edge (m).

shape `property` ¶

shape: tuple[int, int, int]

Grid size as (nx, ny, nz); arrays are (nz, ny, nx).

fields `property` ¶

fields: list[str]

Names of the raw accumulators this grid holds (depends on the options).

metrics `property` ¶

metrics: list[str]

Names of the derived metrics; see the class notes.

observed `property` ¶

observed: ndarray

Boolean (nz, ny, nx): voxels crossed by at least one pulse before its last echo.

tree_iad `property` ¶

tree_iad: dict[int, dict]

Per-tree inclination angle distributions: normalised liad / wiad / piad histograms over bin_centres (rad), the angle-integrated g_leaf / g_wood / g_plant and the closest de Wit type of each. Empty unless inclination=True.

Returns:

Type	Description
`dict`	`{tree_id: {...}}`.

occlusion_profile ¶

occlusion_profile(min_height: float = 0.0, max_height: float | None = None) -> dict

What the scan saw of the canopy space, layer by layer.

The canopy space is every voxel from min_height above the ground (distance_from_ground, needs a DTM; else above the grid floor) up to max_height -- by default the highest layer holding a filled voxel. Returns per layer the height of its centre, the voxel count and the shares observed (a pulse went through or ended in it), occluded (only pulses already stopped reached it; needs occlusion=True) and unobserved, the mean pulses entering a voxel (mean_beams), and plot totals under "total".

Parameters:

Name	Type	Description	Default
`min_height`	`float`	Bottom of the canopy space (m above ground).	`0.0`
`max_height`	`float \| None`	Top of the canopy space; the highest filled voxel if None.	`None`

Returns:

Type	Description
`dict`	Arrays `height`, `n_voxels`, `observed`, `occluded`, `unobserved`, `mean_beams` (one value per layer) and `total`: `{"observed", "occluded", "unobserved", "top"}`.

observed_map ¶

observed_map(min_height: float = 0.0, max_height: float | None = None) -> np.ndarray

Map of how much of each column's canopy space was observed.

Parameters:

Name	Type	Description	Default
`min_height`	`float`	Canopy space, as for `occlusion_profile`.	`0.0`
`max_height`	`float`	Canopy space, as for `occlusion_profile`.	`0.0`

Returns:

Type	Description
`ndarray`	`(ny, nx)` share observed; NaN for columns with no canopy space. Useful to find the parts of a plot to rescan.

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 `(nz, ny, nx)`, aligned with the voxel arrays.

profile ¶

profile(name: str = 'pad_fpl', min_beams: int = 1) -> np.ndarray

Mean of a metric per vertical layer over voxels crossed by at least min_beams pulses (unobserved voxels do not count as empty). Multiply by the voxel size and sum for a plant area index.

Parameters:

Name	Type	Description	Default
`name`	`str`	Field or metric, e.g. `"pad_fpl"` or `"pad_ppl"`.	`'pad_fpl'`
`min_beams`	`int`	Minimum pulses for a voxel to count. Raise it (5-10) to keep poorly sampled voxels out of the mean.	`1`

Returns:

Type	Description
`ndarray`	One value per layer, bottom first; NaN where no voxel qualifies. Layers are in grid z, not height above ground.

add_wood_volume ¶

add_wood_volume(qsms: QSM | Iterable[QSM]) -> None

Add QSM wood volume to the grid.

Rasterises QSM cylinders into the wood_volume (m³) field and the wood_volume_density (m³ m⁻³) metric. Repeated calls add up.

Parameters:

Name	Type	Description	Default
`qsms`	`QSM \| Iterable[QSM]`	One QSM or several, in the grid's frame.	required

to_dict ¶

to_dict(names: Iterable[str] | None = None) -> dict[str, np.ndarray]

Copy arrays out by name.

Parameters:

Name	Type	Description	Default
`names`	`Iterable[str] \| None`	Fields and metrics to include; every raw field if None.	`None`

Returns:

Type	Description
`dict`	`{name: (nz, ny, nx) array}`, e.g. for `np.savez`.

write ¶

write(path: str | Path, format: str | None = None, include_unobserved: bool = False, filled_only: bool = False) -> int

Write a .vox voxel-space file, or a space-delimited .txt table with voxel centres (format "vox" / "text", by default from the extension). Observed and occluded voxels are written unless filled_only or include_unobserved.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	Output file; overwritten.	required
`format`	`('vox', 'text')`	From the extension if None (`.vox` is a voxel-space file, anything else text).	`"vox"`
`include_unobserved`	`bool`	Also write voxels no pulse reached.	`False`
`filled_only`	`bool`	Only write voxels holding echoes.	`False`

Returns:

Type	Description
`int`	Number of voxels written.

write_iad_csv ¶

write_iad_csv(path: str | Path) -> None

Write per-tree inclination distributions as CSV.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	Output file: one row per tree with de Wit types and inclination histograms. Needs `inclination=True`.	required

laser_spec ¶

laser_spec(name: str) -> tuple[float, float]

Beam geometry of a known scanner.

Parameters:

Name	Type	Description	Default
`name`	`str`	Scanner name, e.g. `"VZ-400"`, `"LMS-Q780"`, `"FARO-FOCUS-X330"`.	required

Returns:

Type	Description
`diameter, divergence : float`	Beam diameter at exit (m) and full divergence (rad).

Raises:

Type	Description
`ValueError`	For an unknown scanner; pass `beam=(diameter, divergence)` to `ray_voxelize` instead.

leaf_projection ¶

leaf_projection(theta, lad: str = 'spherical', lad_params: Sequence[float] = ()) -> np.ndarray

Leaf projection function G(θ) of an analytic leaf angle distribution.

G is the mean projected leaf area per unit leaf area in the beam direction; area density is attenuation over G.

Parameters:

Name	Type	Description	Default
`theta`		Beam zenith angle(s) in radians.	required
`lad`	`str`	`spherical` (G = 0.5 everywhere), `uniform`, `planophile`, `erectophile`, `plagiophile`, `extremophile` (de Wit 1965), `ellipsoidal` (Campbell 1990; `lad_params=[chi]`) or `twoParamBeta` (Goel & Strebel 1984; `[mu, nu]`).	`'spherical'`
`lad_params`	`Sequence[float]`	Parameters of the ellipsoidal or beta distribution. The beta distribution follows Goel and Strebel: on `t = 2θ/π` its density is proportional to `(1 - t)^(mu - 1) t^(nu - 1)`, so their fits apply as published: planophile `[2.770, 1.172]`, erectophile `[1.172, 2.770]`, plagiophile `[3.326, 3.326]`, extremophile `[0.433, 0.433]`, uniform `[1, 1]` and spherical `[1.101, 1.930]`. A fit from `sylva.leaves.leaf_angle_distribution` is given as `sylva.leaves.LeafAngleDistribution.goel_strebel`.	`()`

Returns:

Type	Description
`ndarray`	G for each angle.

Raises:

Type	Description
`ValueError`	For an unknown distribution or wrong parameters.

ray_voxelize ¶

ray_voxelize(shots: Shots | str | Path, voxel_size: float = 0.1, bounds=None, *, dtm: Raster | None = None, ground: ndarray | None = None, ground_class: int | None = None, ground_distance: float = 0.2, foliage: ndarray | None = None, leaf_classes: Sequence[int] = (), wood_classes: Sequence[int] = (), class_attr: str = 'classification', tree_attr: str | None = 'tree_id', intensity_attr: str = 'intensity', weighting: str = 'equal', attenuation: str | Sequence[str] = 'fpl', laser: str | None = None, beam: tuple[float, float] | None = None, lad: str = 'spherical', lad_params: Sequence[float] = (), inclination: bool = False, n_iad_bins: int = 18, knn_normal: int = 10, triangle_lmax: float = 0.05, occlusion: bool = False, flat_top: bool = False, neighbour_prior_min_rays: int = 0, subvoxel_split: int = 0, subvoxel_min_beams: int = 10, average_leaf_area: float = 0.005, unbounded_range: float = np.inf, block_size: int | Sequence[int] | None = None, out: str | Path | None = None, max_memory: float | None = None, workers: int | None = None)

Trace shots through a voxel grid (see the module docstring).

shots may be the path of a shots file (sylva.Shots.save). It is then streamed a few row groups at a time, so memory stays at the size of the grid plus a few million pulses; the grid defaults to the echo bounds recorded in the file, and echo labels must come from attributes (ground_class / dtm, leaf_classes / wood_classes, tree_attr) rather than from ground / foliage arrays.

Parameters:

Name	Type	Description	Default
`shots`	`Shots \| str \| Path`	Pulses in one frame (every scan of a plot together), or the path of a shots file. Must include the misses: from `sylva.Shots.fill_missing`, a ray cloud or a shots file made from either.	required
`voxel_size`	`float`	Voxel edge (m). 0.1-0.5 m is usual; memory is about 0.25 kB per voxel while tracing, so a 50 x 50 x 40 m plot at 0.1 m needs ~25 GB unless it is traced in blocks (`block_size`, `out`).	`0.1`
`bounds`		`(min_xyz, max_xyz)`; by default the extent of the echoes. The max corner is snapped up to a whole number of voxels.	`None`
`dtm`	`Raster \| None`	Ground echoes are traced up to but never count as hits. Give them as a boolean `ground` mask per echo, as a `ground_class` code of the `class_attr` echo attribute, or as echoes no more than `ground_distance` above the `dtm` (or below it). The DTM also stops occlusion rays below ground and gives `distance_from_ground`.	`None`
`ground`	`Raster \| None`	Ground echoes are traced up to but never count as hits. Give them as a boolean `ground` mask per echo, as a `ground_class` code of the `class_attr` echo attribute, or as echoes no more than `ground_distance` above the `dtm` (or below it). The DTM also stops occlusion rays below ground and gives `distance_from_ground`.	`None`
`ground_class`	`Raster \| None`	Ground echoes are traced up to but never count as hits. Give them as a boolean `ground` mask per echo, as a `ground_class` code of the `class_attr` echo attribute, or as echoes no more than `ground_distance` above the `dtm` (or below it). The DTM also stops occlusion rays below ground and gives `distance_from_ground`.	`None`
`ground_distance`	`Raster \| None`	Ground echoes are traced up to but never count as hits. Give them as a boolean `ground` mask per echo, as a `ground_class` code of the `class_attr` echo attribute, or as echoes no more than `ground_distance` above the `dtm` (or below it). The DTM also stops occlusion rays below ground and gives `distance_from_ground`.	`None`
`foliage`	`ndarray \| None`	Per-echo codes (`EXCLUDED`, `PLANT`, `LEAF`, `WOOD`), or the `class_attr` codes that mean leaf and wood (other codes below 3 are excluded, the rest are plant, as in rayvoxel). With neither, every non-ground echo is plant.	`None`
`leaf_classes`	`ndarray \| None`	Per-echo codes (`EXCLUDED`, `PLANT`, `LEAF`, `WOOD`), or the `class_attr` codes that mean leaf and wood (other codes below 3 are excluded, the rest are plant, as in rayvoxel). With neither, every non-ground echo is plant.	`None`
`wood_classes`	`ndarray \| None`	Per-echo codes (`EXCLUDED`, `PLANT`, `LEAF`, `WOOD`), or the `class_attr` codes that mean leaf and wood (other codes below 3 are excluded, the rest are plant, as in rayvoxel). With neither, every non-ground echo is plant.	`None`
`class_attr`	`ndarray \| None`	Per-echo codes (`EXCLUDED`, `PLANT`, `LEAF`, `WOOD`), or the `class_attr` codes that mean leaf and wood (other codes below 3 are excluded, the rest are plant, as in rayvoxel). With neither, every non-ground echo is plant.	`None`
`tree_attr`	`str \| None`	Echo attribute grouping echoes into trees for the inclination distributions; all echoes are pooled as tree 0 when it is missing.	`'tree_id'`
`weighting`	`str`	Share of a pulse carried by each echo: `equal` (`1 / n`), `full` (last), `first`, `relative` or `strongest` (by `intensity_attr`).	`'equal'`
`attenuation`	`str \| Sequence[str]`	One or more of `fpl`, `ppl`, `transmittance`, `bailey`; the first is used for `pad_g_corrected`. `ppl` enables the exact solve; `bailey` needs leaf and wood echoes and turns on `inclination`.	`'fpl'`
`laser`	`str \| None`	A scanner name for `laser_spec`, or `(exit diameter [m], divergence [rad])`. Enables the beam-section metrics (`bs_*`, `transmittance`) and weights FPL / PPL by beam section.	`None`
`beam`	`str \| None`	A scanner name for `laser_spec`, or `(exit diameter [m], divergence [rad])`. Enables the beam-section metrics (`bs_*`, `transmittance`) and weights FPL / PPL by beam section.	`None`
`lad`	`str`	Analytic leaf angle distribution, see `leaf_projection`.	`'spherical'`
`lad_params`	`str`	Analytic leaf angle distribution, see `leaf_projection`.	`'spherical'`
`inclination`	`bool`	Estimate inclination angle distributions from echo normals (PCA over `knn_normal` neighbours); `triangle_lmax` caps Bailey facet edges.	`False`
`n_iad_bins`	`bool`	Estimate inclination angle distributions from echo normals (PCA over `knn_normal` neighbours); `triangle_lmax` caps Bailey facet edges.	`False`
`knn_normal`	`bool`	Estimate inclination angle distributions from echo normals (PCA over `knn_normal` neighbours); `triangle_lmax` caps Bailey facet edges.	`False`
`triangle_lmax`	`bool`	Estimate inclination angle distributions from echo normals (PCA over `knn_normal` neighbours); `triangle_lmax` caps Bailey facet edges.	`False`
`occlusion`	`bool`	Also trace beyond each pulse's last echo (`num_beams_occluded`).	`False`
`flat_top`	`bool`	Start the path in each column's top voxel at the highest echo.	`False`
`neighbour_prior_min_rays`	`int`	Top up voxels crossed by fewer weighted beams from their 26 neighbours.	`0`
`subvoxel_split`	`int`	`N` (2-4) for an `N³` sub-voxel grid giving `exploration_rate`.	`0`
`subvoxel_min_beams`	`int`	`N` (2-4) for an `N³` sub-voxel grid giving `exploration_rate`.	`0`
`average_leaf_area`	`float`	Mean leaf area (m²) of the effective free path correction for leaves of finite size (Pimont et al. 2018, 2019); 0 disables.	`0.005`
`unbounded_range`	`float`	How far pulses without an echo are traced (default: to the grid edge).	`inf`
`block_size`	`int \| Sequence[int] \| None`	Trace the grid block by block, a few blocks at a time: voxels per block edge, or `(bx, by, bz)`. Every voxel receives exactly the sums of a whole-grid trace (see the guide, "Block tracing"); memory while tracing is `max_memory` rather than 0.25 kB per voxel of the grid. 64 when only `out` is given.	`None`
`out`	`str \| Path \| None`	Directory to write the blocks to as they finish; the grid is then returned as a `~sylva.voxels.blocks.BlockedVoxelGrid` on disk. Without it the blocks are assembled into one grid in memory (0.13 kB per voxel). Neighbour priors, `inclination` and `bailey` need the assembled grid.	`None`
`max_memory`	`float \| None`	GB the tracing accumulators of one pass may take (default: half of `sylva.util.limits.budget`). Fewer GB means more passes over the pulses.	`None`
`workers`	`int \| None`	Threads for a blocked trace (default: all). Blocked results do not depend on it.	`None`

Returns:

Type	Description
`RayVoxelGrid or BlockedVoxelGrid`	A `~sylva.voxels.blocks.BlockedVoxelGrid` when `out` is given. After a blocked trace, `block_stats` holds the blocks, passes and the most voxels held at once (`peak_voxels`).

Raises:

Type	Description
`ValueError`	For inconsistent options (`laser` and `beam`; arrays with a shots file; arrays of the wrong length; unknown methods; a block that does not fit in `max_memory`).
`KeyError`	If a named echo attribute is missing.

tree_sampling ¶

tree_sampling(grid: RayVoxelGrid, cloud, labels, min_beams: float = 10.0, above: float = 2.0) -> dict[str, np.ndarray]

How well each tree was seen, from a ray-traced grid built with occlusion=True (so occluded voxels are told from unreached ones).

A tree's crown envelope is, voxel layer by voxel layer, the convex hull of its points (labels per point, negative ignored). That envelope is observed almost by construction -- the parts of a crown nobody saw left no points -- so the informative columns are the pulses that reached it and whether its top is real:

above_observed_fraction share of the voxels up to above metres over the tree's highest point, within the footprint of its top metre, that were observed (empty, or filled by a neighbour). Low means the top may be hidden and the tree taller than its points. median_beams, p10_beams, beams_by_quarter pulses entering an envelope voxel (unseen voxels count 0): overall, the 10th percentile, and the median per quarter of the envelope's height, bottom first. well_sampled_fraction observed envelope voxels with at least min_beams pulses.

Also tree_id, n_voxels and volume of the envelope and its observed_fraction, occluded_fraction, unobserved_fraction.

Parameters:

Name	Type	Description	Default
`grid`	`RayVoxelGrid`	Grid from `ray_voxelize` with `occlusion=True`, covering the trees; a `~sylva.voxels.blocks.BlockedVoxelGrid` is read a slab of blocks at a time (it needs 5 bytes per voxel).	required
`cloud`		A `~sylva.PointCloud` or `(N, 3)` array in the grid's frame.	required
`labels`		Tree id per point (e.g. from `sylva.trees.segment_trees`).	required
`min_beams`	`float`	Pulses for a voxel to count as well sampled.	`10.0`
`above`	`float`	Height (m) above each tree's top that is checked.	`2.0`

Returns:

Type	Description
`dict`	The columns above as arrays, one row per tree; ready for `pandas.DataFrame`.

sylva.voxels¶

sylva.voxels ¶

RayVoxelGrid ¶

origin property ¶

voxel_size property ¶

shape property ¶

fields property ¶

metrics property ¶

observed property ¶

tree_iad property ¶

occlusion_profile ¶

observed_map ¶

z_levels ¶

centers ¶

profile ¶

add_wood_volume ¶

to_dict ¶

write ¶

write_iad_csv ¶

laser_spec ¶

leaf_projection ¶

ray_voxelize ¶

tree_sampling ¶

origin `property` ¶

voxel_size `property` ¶

shape `property` ¶

fields `property` ¶

metrics `property` ¶

observed `property` ¶

tree_iad `property` ¶