sylva.change

sylva.change

Change between two epochs of one plot.

The principle throughout is to separate real change from noise and from not having seen something. Every change carries its uncertainty or level of detection, and whatever the data cannot support is labelled (below detection, unmeasured, ambiguous) rather than reported as change.

Trees and plots::

from sylva import change

al = change.align_epochs(cloud_2019, cloud_2024)          # on stems and ground
m = change.match_trees(stems_2019, stems_2024, transform=al)
inc = change.tree_increments(m, cloud_2019, cloud_2024, labels_2019, labels_2024,
                             noise_a=0.003, noise_b=0.004)
s = change.plot_summary(inc, area=np.pi * 20**2, years=5)
change.provenance(settings_2019, settings_2024)           # warns if they differ

Points, surfaces and voxels: distances (C2C and M3C2), dod and occupancy. QSMs: compare_qsms and compare_plot_qsms.

sylva.synthetic.forest_epochs makes two scanned epochs of a plot with known changes, against which all of this is validated.

Airborne surveys (catalogues of tiles, sylva.change.als): align_surveys, chm_change and dtm_change (surface_change), harmonise, canopy_gaps and gap_change, tree_change, metric_change and pai_change.

sylva.change.epochs

Two epochs in one frame, and a record of how each was processed.

EpochAlignment dataclass

Result of align_epochs: the later epoch in the frame of the earlier one, and how well that is known.

Attributes:

Name	Type	Description
transform	ndarray	(4, 4) transform taking the new epoch onto the reference.
registration_sigma	float	One-sigma uncertainty (m) of where a point of the new epoch lands, as the RMS over the reference stems (or evaluate) of the 3-D displacement uncertainty that the covariance implies; NaN if a direction is not constrained (stems alone leave z free).
sigma_horizontal, sigma_vertical		Its horizontal and vertical parts. The vertical part is what moves stem slices along the stem when increments are measured.
sigma_xyz	ndarray	One-sigma uncertainty of x, y, z at centre.
covariance	ndarray	(6, 6) covariance of a small rotation (rad, about x, y, z through centre) and a translation (m).
centre	ndarray	Point the parameters are expressed about.
stem_pairs	ndarray	(n, 2) indices of the matched stems in the two stem maps (reference.stem_map, new.stem_map).
stem_residuals	ndarray	(n, 2) horizontal residual (m) of each matched stem after alignment.
ground_residuals	ndarray	Height (m) of each terrain sample of the new epoch above the reference terrain after alignment.
stem_sigma, ground_sigma		Robust spread (m) of one stem residual component and of one terrain residual: the precision of a single stable feature.
coarse	PairResult | None	The coregistration pair that gave the starting transform (None when initial was given).
reference, new		The prepared epochs (sylva.coreg.ScanFeatures).
stable	str	The features used.
settings	dict	The settings of the call, for provenance.

n_stems property

n_stems: int

Stems matched between the epochs.

apply

apply(cloud: PointCloud) -> PointCloud

The new epoch's cloud moved into the reference frame.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	A cloud in the new epoch's frame.	required

Returns:

Type	Description
PointCloud

report

report() -> str

Plain-text summary for a log or a QC record.

ProvenanceWarning

Bases: UserWarning

Two epochs were processed with different software or settings.

Provenance dataclass

Result of provenance.

Attributes:

Name	Type	Description
records	list[dict]	One per epoch: sylva_version, settings and whatever else was recorded (created, info).
differences	list[tuple]	(setting, value in the first epoch, value in the other) for every setting that differs, with dotted keys (settings.detect_stems.min_slices) and values as text (None where an epoch lacks the setting). With more than two epochs, each is compared with the first and the key is prefixed epoch <k>:.

consistent property

consistent: bool

Were all epochs processed with the same version and settings?

to_json

to_json() -> str

The records as JSON text, to store beside each epoch's products.

align_epochs

align_epochs(ref, new, stable: str = 'stems+ground', config: CoregConfig | None = None, initial: ndarray | None = None, stem_tolerance: float = 0.3, dbh_tolerance: float = 0.35, ground_spacing: float = 1.0, ground_block: float = 5.0, stem_floor: float = 0.002, ground_floor: float = 0.005, huber: float = 2.0, iterations: int = 30, evaluate=None) -> EpochAlignment

Align a later epoch onto an earlier one on features that do not change.

Between surveys crowns grow, lose limbs and move in the wind, so a registration over all points is pulled by the change it should reveal. Here a coarse transform comes from sylva.coreg.register_pair (a global stem-map match refined by ICP), and is then refined on the stable features only: the horizontal position of every stem matched in both epochs (a stem thickens, its axis stays put) and the height of the new terrain above the reference terrain. Robust Gauss-Newton with Huber (1964) weights, each feature type weighted by its own spread, gives the transform and the covariance of its six parameters, from which registration_sigma follows.

Parameters:

Name	Type	Description	Default
ref		The two epochs: ~sylva.PointCloud, (n, 3) arrays, file paths, or scans already prepared with sylva.coreg.prepare_scan. For a multi-scan project, pass its registered, merged cloud (sylva.coreg.merge_clouds).	required
new		The two epochs: ~sylva.PointCloud, (n, 3) arrays, file paths, or scans already prepared with sylva.coreg.prepare_scan. For a multi-scan project, pass its registered, merged cloud (sylva.coreg.merge_clouds).	required
stable	str	"stems+ground" (default), "stems" (x, y and rotation about z only; z and the tilts stay at the coarse transform's) or "ground" (z and the tilts only).	'stems+ground'
config	CoregConfig | None	Settings of the coarse coregistration and of the stem maps and terrain models.	None
initial	ndarray | None	(4, 4) starting transform (new onto reference); skips the coarse coregistration.	None
stem_tolerance	float	Farthest (m) a stem may lie from its match under the current transform to count as stable.	0.3
dbh_tolerance	float	Largest relative DBH difference of a stable stem pair.	0.35
ground_spacing	float	Spacing (m) of the terrain samples.	1.0
ground_block	float	Terrain errors are taken as correlated within blocks of this side (m): all terrain samples together weigh as much as one independent sample per block, so the smooth, spatially coherent errors of a terrain model do not pass for precision.	5.0
stem_floor	float	Smallest error (m) assumed for one stem position (per axis) and one terrain sample. No bark surface or terrain model is known better, whatever the residuals of an unusually clean plot suggest; on synthetic plots without these floors the alignment error exceeded its sigma.	0.002
ground_floor	float	Smallest error (m) assumed for one stem position (per axis) and one terrain sample. No bark surface or terrain model is known better, whatever the residuals of an unusually clean plot suggest; on synthetic plots without these floors the alignment error exceeded its sigma.	0.002
huber	float	Residuals beyond this many robust standard deviations are down-weighted.	2.0
iterations	int	Most Gauss-Newton iterations.	30
evaluate		(n, 3) points (reference frame) over which registration_sigma is summarised; the reference stems if None.	None

Returns:

Type	Description
EpochAlignment

Raises:

Type	Description
ValueError	If stable is unknown, or the features do not constrain the alignment (too few stems matched and no shared terrain).

Notes

registration_sigma is the uncertainty of the transform, not the residual of one feature: a stem-map residual of 5 mm over 15 stems constrains the translation to about 5 / sqrt(15) mm. Stem errors are taken as independent, terrain errors as shared within ground_block. A bias common to all stems or to the whole terrain of one epoch (a different terrain filter, say) is not visible in the residuals and not included.

Examples:

>>> ep = synthetic.forest_epochs(seed=1)
>>> al = change.align_epochs(*ep.clouds)
>>> new_in_ref = al.apply(ep.clouds[1])
>>> al.registration_sigma

provenance

provenance(*epochs, warn: bool = True, rtol: float = 1e-09, ignore=('created', 'info')) -> Provenance

Record how each epoch was processed, and compare the epochs.

A difference between two epochs is change only if both were measured the same way: another stem-detection setting, voxel size or Sylva version can move a DBH or a profile by more than a year's growth. Record each epoch when it is processed (provenance(settings).records[0], stored with its products), then compare the records of the epochs to be compared.

Parameters:

Name	Type	Description	Default
*epochs		Per epoch, the settings it was processed with (a dict, possibly nested, or a dataclass such as sylva.coreg.CoregConfig; the settings of this module's results can be included), or a record made earlier by this function.	()
warn	bool	Emit a ProvenanceWarning when the epochs differ.	True
rtol	float	Numbers within this relative tolerance count as equal.	1e-09
ignore		Top-level record keys left out of the comparison.	('created', 'info')

Returns:

Type	Description
Provenance

Examples:

>>> rec_a = change.provenance({"detect_stems": {"min_slices": 3}}).records[0]
>>> rec_b = change.provenance({"detect_stems": {"min_slices": 4}}).records[0]
>>> change.provenance(rec_a, rec_b).differences   # warns
[('settings.detect_stems.min_slices', '3', '4')]

sylva.change.trees

Trees matched between epochs, and their increments.

TreeMatch dataclass

Result of match_trees.

Attributes:

Name	Type	Description
trees_a, trees_b		The trees of each epoch, as given (each in its own frame).
pairs	ndarray	(n, 2) indices into trees_a and trees_b of every survivor, sorted by the first.
distance	ndarray	Horizontal distance (m) of each pair in the reference frame.
cost	ndarray	Assignment cost of each pair.
status_a	ndarray	Per tree of the first epoch: survivor, dead, merged (found in the second epoch only as part of a neighbour's stem) or split (found as several stems, none of them matched to it).
status_b	ndarray	Per tree of the second epoch: survivor, recruit, split (part of a stem the first epoch saw as one with a neighbour) or merged (one stem standing for several of the first epoch, none of them matched to it).
related_a, related_b		The partner (for a survivor) or the stem merged into / split from, as an index into the other epoch's trees; -1 for deaths and recruits.
transform	ndarray	The transform applied to the second epoch's positions (identity if none was given).
alignment	EpochAlignment | None	The EpochAlignment given as transform, if any.
settings	dict	The settings of the call, for provenance.

survivors property

survivors: list[tuple]

(tree in epoch a, tree in epoch b) of every survivor.

deaths property

deaths: list

Trees of the first epoch not found in the second.

recruits property

recruits: list

Trees of the second epoch not found in the first.

merged property

merged: list[tuple]

(tree in epoch a, stem in epoch b it went into) of the trees of the first epoch that merged into a neighbour's stem.

split property

split: list[tuple]

(stem in epoch a it came from, tree in epoch b) of the trees of the second epoch that split off a stem of the first.

ambiguous_pairs

ambiguous_pairs() -> np.ndarray

Survivor pairs involved in a merge or a split, whose increments compare a stem with more (or less) than itself.

Returns:

Type	Description
ndarray	One bool per row of pairs.

counts

counts() -> dict

Numbers of survivors, deaths and recruits, and of the trees of either epoch in a merge or a split.

TreeIncrements dataclass

Result of tree_increments: one row per survivor.

Columns (columns[name], or increments[name]), lengths in m:

tree_id_a, tree_id_b, x, y The tree in each epoch and its position in the first. dbh_a, dbh_a_se, dbh_b, dbh_b_se DBH of each epoch at 1.3 m from a weighted line through the stem slices, and its standard error. d_dbh, d_dbh_se, d_dbh_mdi, n_slices DBH increment from the paired slices, its standard error, the minimum detectable increment at the chosen confidence, and the number of slices measured in both epochs. dbh_change growth (at least the MDI), below_detection, decrease (a significant decrease) or unmeasured. height_a, height_b, d_height, d_height_se, d_height_mdi, height_change Top height of each epoch, and the same for its increment. crown_area_a, crown_area_b, d_crown_area, crown_volume_a, crown_volume_b, d_crown_volume Crown projection area (m²) and stacked-hull volume (m³); descriptive, without a detection level: occlusion changes them as much as growth. implausible A significant DBH decrease, or an increase above max_dbh_increment. ambiguous The pair is part of a merge or split (TreeMatch.ambiguous_pairs). flags The row's warnings as text: unmeasured, implausible, ambiguous, height_decrease.

Attributes:

Name	Type	Description
columns	dict	The table, a dict of equal-length arrays.
match	TreeMatch	The TreeMatch the rows come from.
measures	list	Per epoch, the measurement of every tree (not only survivors): a dict of arrays with dbh, dbh_se, taper, height, height_se, crown_area, crown_volume, n_points and the (n, k) slice diameter and diameter_se.
slice_heights	ndarray	Heights (m) of the stem slices.
settings	dict	The settings of the call, for provenance.

as_dict

as_dict() -> dict

The table as a dict of arrays (a copy).

to_pandas

to_pandas()

The table as a pandas.DataFrame (needs pandas).

match_trees

match_trees(trees_a, trees_b, max_distance: float = 1.0, transform=None, dbh_tolerance: float = 0.35, max_shrink: float = 0.15, dbh_weight: float = 1.0, height_weight: float = 0.25, merge_factor: float = 0.8) -> TreeMatch

Match the trees of two epochs: survivors, deaths and recruits.

The second epoch's stem positions are moved by transform into the first epoch's frame, and the pairs are chosen by an optimal assignment (Kuhn 1955; Munkres 1957) that first matches as many trees as it can and then minimises the total cost (d / max_distance)² + dbh_weight * ddbh² + height_weight * dh², with ddbh and dh the relative differences ((b - a) / max(a, b)). A pair is allowed only within max_distance, and only if the DBH grew by at most dbh_tolerance or shrank by at most max_shrink, so a felled tree and a young one grown beside it are not taken for one tree. Unlike nearest-neighbour matching, a neighbour's closer stem cannot take a tree's partner when the pair would leave the neighbour unmatched. The assignment is solved in each group of trees that could be matched to one another, so large plots cost little.

Trees left over are checked for a merge before being called dead: when a stem of the second epoch is wide enough (merge_factor times the quadratic sum of the DBHs) to stand for its partner and an unmatched neighbour of the first epoch, the neighbour is merged into it, as happens when the detector sees two touching stems as one. Splits are the mirror case.

Parameters:

Name	Type	Description	Default
trees_a		Trees of the two epochs (sylva.trees.Tree, e.g. from sylva.trees.detect_stems and sylva.trees.tree_heights), or tables (dicts of arrays with x, y and optionally tree_id, dbh, height), each in its own frame.	required
trees_b		Trees of the two epochs (sylva.trees.Tree, e.g. from sylva.trees.detect_stems and sylva.trees.tree_heights), or tables (dicts of arrays with x, y and optionally tree_id, dbh, height), each in its own frame.	required
max_distance	float	Farthest apart (m) the two stems of one tree can be.	1.0
transform		EpochAlignment or (4, 4) matrix taking the second epoch onto the first; the epochs are taken to share a frame if None.	None
dbh_tolerance	float	Largest relative DBH increase of one tree; NaN DBHs are not tested.	0.35
max_shrink	float	Largest relative DBH decrease of one tree. Stems do not shrink, so this only absorbs the DBH error of the two detections.	0.15
dbh_weight	float	Weights of the relative DBH and height differences in the cost.	1.0
height_weight	float	Weights of the relative DBH and height differences in the cost.	1.0
merge_factor	float	Test for merges and splits, as above.	0.8

Returns:

Type	Description
TreeMatch

Raises:

Type	Description
ValueError	If a position is not finite or a setting is out of range.

Examples:

>>> m = change.match_trees(stems_a, stems_b, transform=alignment)
>>> m.counts()
{'survivors': 14, 'deaths': 2, 'recruits': 2, 'merged': 0, 'split': 0}

tree_increments

tree_increments(match: TreeMatch, cloud_a: PointCloud, cloud_b: PointCloud, labels_a, labels_b, noise_a=None, noise_b=None, registration_sigma=None, height_attr: str = 'height', slice_heights=None, slice_thickness: float = 0.1, search_radius: float = 0.75, min_slices: int = 3, confidence: float = 0.95, max_dbh_increment: float | None = None, top_points: int = 5, height_error: float = 0.02, top_radius: float = 0.3, top_gap: float = 0.5, slice_correlation: float = 0.3, dbh_accuracy: float = 0.002) -> TreeIncrements

DBH, height and crown increments of the survivors, each with the smallest increment the data can detect.

Each epoch's stems are cut into thin slices (1 to 3 m every 0.25 m by default) and a circle fitted to each: RANSAC, then least squares on the inliers, with one inlier distance for both epochs (three times the larger range noise, at least 1 cm). Slices far off the taper line of the others (a branch junction) are dropped. The DBH increment is the weighted mean of the paired differences at the same heights, which cancels the stem's own irregularity (a flattened or fluted stem has the same shape in both epochs). Its standard error combines

• each slice's circle-fit precision, 2 s / sqrt(n) / sqrt(c) with s the fit residual (at least the epoch's range noise), n the inlier points and c the share of the circumference seen;
• the scatter of the differences between slices, where it exceeds that (Birge ratio);
• a share slice_correlation of the slice errors that does not average out, because the same scanners see every slice of a stem from the same directions;
• the vertical registration uncertainty times the stem taper.

The minimum detectable increment (MDI) is that error times the normal quantile of confidence: an increment smaller than it is reported as below_detection, not as growth.

Height is the highest point of the tree above its epoch's own terrain, continued upwards through points the segmentation left unassigned within top_radius of the stem (segmentation often loses a thin leader); vertical misregistration cancels because each epoch is normalised by its own terrain. Its uncertainty combines how far the next highest points (top_points) lie below the top with an assumed error of height_error times the height for a top that was missed: a thin leader that no pulse hit, or a top hidden behind the crown, leaves no trace in the epoch's own data. Crown area and volume are reported as measured.

Parameters:

Name	Type	Description	Default
match	TreeMatch	From match_trees; its trees' x, y place the slices.	required
cloud_a	PointCloud	The segmented clouds of the two epochs, each in its own frame, with height above ground in height_attr.	required
cloud_b	PointCloud	The segmented clouds of the two epochs, each in its own frame, with height above ground in height_attr.	required
labels_a		Tree label of each point (sylva.trees.segment_trees), equal to the trees' tree_id; -1 for none.	required
labels_b		Tree label of each point (sylva.trees.segment_trees), equal to the trees' tree_id; -1 for none.	required
noise_a		Range noise of each epoch (m), or a sylva.quality.StemNoise result whose sigma_local is used; a floor on the circle-fit residual. None uses the fit residual alone.	None
noise_b		Range noise of each epoch (m), or a sylva.quality.StemNoise result whose sigma_local is used; a floor on the circle-fit residual. None uses the fit residual alone.	None
registration_sigma		Vertical one-sigma alignment uncertainty (m), or an EpochAlignment (its sigma_vertical); by default the alignment match was made with, else 0.	None
height_attr	str	Attribute holding height above ground.	'height'
slice_heights		Heights (m) of the stem slices; 1 to 3 m every 0.25 m if None.	None
slice_thickness	float	Slice thickness (m).	0.1
search_radius	float	Only points this close (m) to the stem centre enter a slice.	0.75
min_slices	int	Fewest slices measured in both epochs for an increment.	3
confidence	float	Two-sided confidence level of the MDI.	0.95
max_dbh_increment	float | None	Increments above this (m) are flagged implausible; None for no limit.	None
top_points	int	Highest points used for the height uncertainty.	5
height_error	float	One-sigma error of a top height as a fraction of the height. The default of 2 % covers the missed tops of the synthetic validation (multi-scan, 0.25 degree); calibrate it against repeated or destructive measurements for other scanners and stands.	0.02
top_radius	float	Unassigned points (label -1) within top_radius (m) of the stem extend the tree upwards while each lies within top_gap (m) of the top reached so far; 0 disables.	0.3
top_gap	float	Unassigned points (label -1) within top_radius (m) of the stem extend the tree upwards while each lies within top_gap (m) of the top reached so far; 0 disables.	0.3
slice_correlation	float	Correlation of the slice errors of one stem in one epoch, 0 to 1.	0.3
dbh_accuracy	float	Absolute one-sigma accuracy (m) of a single-epoch DBH beyond its precision (bark roughness and edge points widen every circle of a stem alike). It cancels in increments and is added only to dbh_a_se and dbh_b_se, which plot_summary uses for dead trees and recruits.	0.002

Returns:

Type	Description
TreeIncrements

Raises:

Type	Description
ValueError	If the labels do not fit the clouds or a setting is out of range.

Notes

On twenty sylva.synthetic.forest_epochs plots (range noise 3 and 5 mm, five scans per epoch) the MDI of the DBH increment was 1.4 mm (median) and every one of the 280 true increments lay within the MDI of its measurement; with two echoes per pulse, whose displaced edge points widen the circles unevenly, 93 % did. See Change detection in the guide.

Examples:

>>> inc = change.tree_increments(m, cloud_a, cloud_b, labels_a, labels_b,
...                              noise_a=0.003, noise_b=0.005)
>>> inc["dbh_change"]

sylva.change.summary

Plot-level growth, mortality and recruitment between two epochs.

Estimate

Bases: NamedTuple

One plot-level quantity with its interval.

PlotSummary dataclass

Result of plot_summary.

Quantities (summary[name] gives an Estimate):

n_survivors, n_deaths, n_recruits, n_ambiguous Tree counts; ambiguous trees (merges and splits) are left out of everything else. mortality_rate, recruitment_rate Annual rates, 1 - (S / N0) ** (1 / years) and 1 - (S / N1) ** (1 / years) (Sheil et al. 1995), with S survivors and N0, N1 the trees of each epoch. dbh_increment Mean DBH increment of the survivors (m/yr). basal_area_growth, basal_area_mortality, basal_area_recruitment, basal_area_net m²/ha/yr: the survivors' increase, the basal area of the dead trees (first epoch), of the recruits (second epoch), and growth plus recruitment minus mortality. volume_growth, volume_mortality, volume_recruitment, volume_net The same for stem volume (form_factor x basal area x height), m³/ha/yr. biomass_* Volume times wood_density (Mg/ha/yr), when a density is given.

Attributes:

Name	Type	Description
quantities	dict	{name: Estimate} in the order above.
settings	dict	The settings of the call, for provenance.

as_dict

as_dict() -> dict

{name: (estimate, low, high, unit)} as plain tuples.

report

report() -> str

Plain-text table of the quantities and their intervals.

plot_summary

plot_summary(increments: TreeIncrements, area: float, years: float = 1.0, n_draws: int = 2000, confidence: float = 0.95, form_factor: float = 0.5, wood_density: float | None = None, seed: int = 0) -> PlotSummary

Growth, mortality and recruitment of a plot, with intervals.

Every quantity is computed from the measured trees, and again for n_draws Monte Carlo draws in which each tree's DBH, height and increments are perturbed by their standard errors (the registration uncertainty is part of the increments' errors). Survivors whose increment could not be measured take the mean measured increment in the estimate and a randomly drawn measured survivor's increment in each draw. The interval is the central confidence range of the draws.

Parameters:

Name	Type	Description	Default
increments	TreeIncrements	From tree_increments; its match gives the deaths and recruits, its measurements their DBH and height. A tree without a stem profile falls back to the DBH and height of its tree object, with no measurement error.	required
area	float	Plot area (m²). Restrict the trees to the plot beforehand so that trees and area cover the same ground.	required
years	float	Time between the epochs.	1.0
n_draws	int	Monte Carlo draws.	2000
confidence	float	Central probability of the intervals.	0.95
form_factor	float	Stem volume over basal area times height.	0.5
wood_density	float | None	Oven-dry wood density (t/m³) for biomass; biomass is left out if None.	None
seed	int	Random seed of the draws.	0

Returns:

Type	Description
PlotSummary

Raises:

Type	Description
ValueError	If a setting is out of range or a tree has no DBH at all.

Notes

The intervals cover measurement noise and registration only. The counts of deaths and recruits are taken as exact (a wrong match is not modelled), and the sampling error of a plot as a sample of a stand is not included. On repeated sylva.synthetic.forest_epochs plots the 95 % intervals of the basal-area terms cover the true values in about 95 % of the runs; see Change detection in the guide.

Examples:

>>> s = change.plot_summary(inc, area=30 * 30, years=5)
>>> s["basal_area_net"]

sylva.change.points

Surface change between two epochs: point distances and rasters of difference.

• distances compares two point clouds, either by the distance to the nearest point (cloud-to-cloud) or by M3C2 (Lague et al. 2013), which also gives each distance its 95 % level of detection.
• dod subtracts two rasters (a DTM or CHM of difference) and marks the cells whose change exceeds the level of detection.

Both epochs must be in one frame (align_epochs in this package registers them); the registration error left after alignment enters the level of detection through registration_sigma.

PointDistances dataclass

Distances between two epochs at a set of core points.

Attributes:

Name	Type	Description
method	{'c2c', 'm3c2'}	How the distances were measured.
core_points	ndarray	(N, 3) positions the distances refer to.
distance	ndarray	(N,) distance (m). For "c2c" the unsigned distance to the nearest point of the reference epoch; for "m3c2" the signed distance from the reference to the compared surface along normal (positive where the surface moved the way the normal points). NaN where it could not be measured.
lod	ndarray | None	(N,) 95 % level of detection (m), M3C2 only; NaN where either cylinder held fewer than min_points points.
significant	ndarray | None	(N,) bool, abs(distance) > lod; M3C2 only.
normal	ndarray | None	(N, 3) unit normals used by M3C2.
n_a, n_b		Points of each epoch in the projection cylinder (M3C2).
spread_a, spread_b		Standard deviation (m, n - 1 denominator) of each epoch's positions along the normal inside the cylinder: the local roughness plus the range noise (M3C2).

to_cloud

to_cloud() -> PointCloud

The core points with the results as attributes.

Returns:

Type	Description
PointCloud	distance and, for M3C2, lod, significant, nx, ny, nz (the normal), n_a, n_b, spread_a and spread_b, ready for sylva.write.

DoD dataclass

A raster of difference with its level of detection.

Attributes:

Name	Type	Description
difference	Raster	raster_b - raster_a on the overlap of the two rasters (m); NaN where either is NaN.
lod	Raster	Level of detection per cell (m); NaN where it is unknown.
significant	ndarray	(rows, cols) bool, abs(difference) > lod.
volume_gained, volume_lost		Volume (m³) raised and lowered over the significant cells (both positive).
net_volume	float	volume_gained - volume_lost (m³).
area_changed	float	Area of the significant cells (m²).
area_compared	float	Area of the cells with a finite difference and level of detection (m²).

thresholded

thresholded() -> Raster

The difference where it is significant, 0 elsewhere and NaN where it could not be assessed.

Returns:

Type	Description
Raster

distances

distances(a, b, method: str = 'm3c2', core_points=None, *, normal_scale: float | None = None, projection_scale: float | None = None, max_depth: float | None = None, registration_sigma: float = 0.0, normals=None, orientation='up', min_points: int = 4, max_distance: float | None = None) -> PointDistances

Distances between two epochs of a point cloud.

"c2c" gives, for each point of b, the distance to the nearest point of a: quick, unsigned, and biased upwards by point spacing and noise, so it suits a first look rather than a test of change.

"m3c2" is the Multiscale Model to Model Cloud Comparison of Lague, Brodu and Leroux (2013). At each core point a plane is fitted to the points of a within normal_scale / 2; both clouds are then projected onto its normal inside a cylinder of diameter projection_scale reaching max_depth either way, and the distance is the difference of the two mean positions. Its 95 % level of detection is (their eq. 1)::

LoD95 = 1.96 * (sqrt(s_a**2 / n_a + s_b**2 / n_b) + reg)

where s is the spread of each cloud's positions along the normal, n their count and reg the registration error registration_sigma. A distance is significant when it exceeds this level. Every core point is computed on its own in parallel (k-d trees on both clouds), so tens of millions of points are practical and the result does not depend on the thread count; use a spatially thinned copy of a as core points to save time.

Parameters:

Name	Type	Description	Default
a		Reference (earlier) epoch, a ~sylva.PointCloud or an (N, 3) array. Points with a non-finite coordinate are ignored.	required
b		Compared (later) epoch, in the same frame.	required
method	('m3c2', 'c2c')	Comparison.	"m3c2"
core_points		Where distances are reported (cloud or (M, 3) array). By default the points of a for M3C2 and of b for C2C; C2C accepts no other core points.	None
normal_scale	float | None	M3C2: diameter D (m) of the neighbourhood the normal is fitted to. Lague et al. suggest about 20 to 25 times the roughness of the surface; for stems a scale below the stem diameter keeps the normal radial.	None
projection_scale	float | None	M3C2: diameter d (m) of the projection cylinder; large enough to hold at least about 30 points of each cloud.	None
max_depth	float | None	M3C2: how far along the normal, either way, points are searched (m); at least the largest change expected. Defaults to normal_scale.	None
registration_sigma	float	M3C2: registration error (m, one standard deviation), for example EpochAlignment.registration_sigma.	0.0
normals		M3C2: (M, 3) normals at the core points to use instead of fitting them (they are normalised, not re-oriented).	None
orientation		M3C2: sign of the fitted normals. "up" (+z; ground and other near-horizontal surfaces), a direction (dx, dy, dz) given as ("direction", (dx, dy, dz)), or a location such as the scanner position given as ("towards", (x, y, z)), which makes stem normals point outwards to it so that growth is positive.	'up'
min_points	int	M3C2: fewest points of each epoch in the cylinder for a level of detection (Lague et al. regard the estimate as reliable from about 4); below it lod is NaN and the point is never significant, though a distance is still given.	4
max_distance	float | None	C2C: distances above this are NaN (m).	None

Returns:

Type	Description
PointDistances

Raises:

Type	Description
ValueError	For an unknown method, missing or non-positive M3C2 scales, a negative registration_sigma, an invalid orientation, normals not matching the core points, or core points given to C2C.

Examples:

>>> d = sylva.change.distances(epoch1, epoch2, "m3c2", core_points=core,
...                            normal_scale=0.5, projection_scale=0.3,
...                            max_depth=1.0, registration_sigma=0.005)
>>> d.distance[d.significant]

dod

dod(raster_a: Raster, raster_b: Raster, min_detectable=None, sigma_a=None, sigma_b=None) -> DoD

Raster of difference between two epochs (a DTM or CHM of difference).

The difference is raster_b - raster_a. A cell's change is significant when it exceeds the level of detection, which is either given (min_detectable) or propagated from the standard deviations of the two surfaces as::

LoD95 = 1.96 * sqrt(sigma_a**2 + sigma_b**2)

(Brasington et al. 2003; Wheaton et al. 2010), which assumes the errors of the two surfaces are independent and normal. Include the vertical registration error in the sigmas, for example sigma_a = np.hypot(dtm_sigma, registration_sigma).

Parameters:

Name	Type	Description	Default
raster_a	Raster	Earlier and later surface, with the same resolution on the same lattice (corners a whole number of cells apart); the result covers their overlap.	required
raster_b	Raster	Earlier and later surface, with the same resolution on the same lattice (corners a whole number of cells apart); the result covers their overlap.	required
min_detectable		Level of detection (m): a number, or a ~sylva.Raster on the same lattice.	None
sigma_a		Standard deviation of each surface (m): numbers or rasters on the same lattice. Give one and the other is taken as equal to it.	None
sigma_b		Standard deviation of each surface (m): numbers or rasters on the same lattice. Give one and the other is taken as equal to it.	None

Returns:

Type	Description
DoD

Raises:

Type	Description
ValueError	If the rasters differ in resolution or lattice or do not overlap, if neither or both of min_detectable and the sigmas are given, or if a threshold is negative.

Examples:

>>> change = sylva.change.dod(chm_2020, chm_2024, sigma_a=0.15)
>>> change.net_volume, change.area_changed

sylva.change.voxels

Occupancy change between two ray-traced voxel grids.

A voxel that holds no echoes in the later epoch has only lost its contents if pulses went through it; if none did, the later scan simply did not see it. occupancy uses the ray-traced grids of sylva.voxels.ray_voxelize to tell the two apart.

Occupancy dataclass

Voxel occupancy change between two epochs, from occupancy.

Arrays are (nz, ny, nx) like those of ~sylva.voxels.RayVoxelGrid.

Attributes:

Name	Type	Description
classes	ndarray	uint8 class code per voxel (OCCUPANCY_CLASSES).
pad_a, pad_b		Plant area density of each epoch (m² m⁻³), the pad metric of the two grids.
pad_change	ndarray	pad_b - pad_a where both epochs sent at least min_pulses pulses into the voxel, NaN elsewhere.
layers	dict[str, ndarray]	Per voxel layer, bottom first: z (layer centre, grid frame), the voxel count of every class (unobserved, stable_empty, stable_occupied, gained, lost), n_compared (voxels in the PAD comparison), and the mean pad_a, pad_b and their difference pad_change over those voxels (NaN where there are none). Multiply pad_change by the voxel size and sum for the change in plant area index over the commonly observed space.
origin	ndarray	Minimum corner of the grids.
voxel_size	float	Voxel edge (m).
min_pulses, min_hits, alpha		Criteria used.

mask

mask(name: str) -> np.ndarray

Boolean (nz, ny, nx) mask of one class.

Parameters:

Name	Type	Description	Default
name	str	A key of OCCUPANCY_CLASSES, e.g. "lost".	required

Returns:

Type	Description
ndarray

counts

counts() -> dict[str, int]

Number of voxels in each class.

Returns:

Type	Description
dict	{class name: count}.

volume

volume(name: str) -> float

Volume (m³) of the voxels in one class.

Parameters:

Name	Type	Description	Default
name	str	A key of OCCUPANCY_CLASSES.	required

Returns:

Type	Description
float

centers

centers(name: str) -> np.ndarray

Centres of the voxels in one class, e.g. to plot what was lost.

Parameters:

Name	Type	Description	Default
name	str	A key of OCCUPANCY_CLASSES.	required

Returns:

Type	Description
ndarray	(K, 3) coordinates.

occupancy

occupancy(grid_a: RayVoxelGrid, grid_b: RayVoxelGrid, pad: str = 'pad_fpl', min_pulses: int = 10, min_hits: int = 1, alpha: float = 0.05) -> Occupancy

Voxel occupancy change between two epochs of a plot.

In each epoch a voxel is occupied when it holds at least min_hits echoes, empty when it holds none and at least min_pulses pulses entered it before their last echo, and not observed otherwise: no pulse reached it, only pulses already stopped did (occluded), or too few passed to call it empty. A voxel occupied in one epoch and empty in the other is only called lost (or gained) if the empty epoch sent enough pulses to have found its contents: if the occupied epoch saw a share p of the pulses entering the voxel stopped there, n pulses all pass through contents as dense with probability (1 - p) ** n, which must not exceed alpha. Sparse contents grazed by a few pulses are thus not reported as change. The classes are then

Class	Epoch a	Epoch b
stable_empty	empty	empty
stable_occupied	occupied	occupied
gained	empty	occupied
lost	occupied	empty
unobserved	not observed in one or both epochs, or too few pulses to tell

so foliage the later scan could not see is unobserved, never lost. Plant area density change is reported per voxel and per layer over the voxels that both epochs sampled with at least min_pulses pulses, so that the layer means compare like with like.

Parameters:

Name	Type	Description	Default
grid_a	RayVoxelGrid	Earlier and later grids from sylva.voxels.ray_voxelize, of pulses registered to one frame and traced with the same voxel_size and bounds. occlusion=True is not needed for the classes (occluded and unreached voxels are both unobserved) but lets ~sylva.voxels.RayVoxelGrid.occlusion_profile say which.	required
grid_b	RayVoxelGrid	Earlier and later grids from sylva.voxels.ray_voxelize, of pulses registered to one frame and traced with the same voxel_size and bounds. occlusion=True is not needed for the classes (occluded and unreached voxels are both unobserved) but lets ~sylva.voxels.RayVoxelGrid.occlusion_profile say which.	required
pad	str	Plant area density metric compared, e.g. "pad_fpl" or "pad_ppl".	'pad_fpl'
min_pulses	int	Fewest pulses entering a voxel for it to count as empty, and for its density to be compared.	10
min_hits	int	Fewest echoes for a voxel to count as occupied; 2 or more discounts isolated noise echoes.	1
alpha	float	Largest probability that the empty epoch missed contents that were still there, for a voxel to be called lost or gained (1 disables the test).	0.05

Returns:

Type	Description
Occupancy

Raises:

Type	Description
ValueError	If the grids differ in origin, voxel size or shape, if pad is not a metric of the grids, if min_pulses or min_hits is below 1, or if alpha is not in (0, 1].

Examples:

>>> ga = sylva.voxels.ray_voxelize(shots_2020, 0.25, bounds)
>>> gb = sylva.voxels.ray_voxelize(shots_2024, 0.25, bounds)
>>> occ = sylva.change.occupancy(ga, gb, min_pulses=10)
>>> occ.volume("lost"), occ.layers["pad_change"]

QSM change

sylva.change.qsm

Change between the quantitative structure models of two epochs.

compare_qsms compares two cylinder models of one tree, the epochs aligned beforehand: the stem radius change by height (the taper increment), branch matching (matched branches with their growth, lost and new branches), volume change by branch order, and height, DBH and crown change.

Every quantity is labelled trusted only where both models fitted that part to points. A cylinder counts as measured when its n_points is positive, as in sylva.qsm.QSM.metrics; the rest of a model comes from the taper and pipe-model priors, and a change there is reported, but as untrusted. With a ray-traced grid of the later epoch (sylva.voxels.ray_voxelize), a branch missing from the later model is lost only where the later scans saw its space empty; where they did not see it, it is unobserved.

compare_plot_qsms does the same for every tree of a plot, given a tree match (survivor pairs, deaths and recruits), and totals the changes.

QSMChange dataclass

Change between two models of one tree; from compare_qsms.

Model a is the earlier epoch and b the later; every change is b - a. Lengths m, volumes m³, angles degrees.

Attributes:

Name	Type	Description
volume_a, volume_b		Whole-tree wood volume of each model.
trusted_change, untrusted_change		volume_b - volume_a split in two: the part both models measured (stem bins and branches fitted to points in both, lost branches confirmed empty, new branches measured) and the rest. They sum to change.
trusted_sigma	float	Standard uncertainty of trusted_change, from the scatter of the measured radii (stem bins and branches) and radius_sigma.
measured_volume_a, measured_volume_b		Share of each model's volume in measured cylinders.
orders	dict[str, ndarray]	Volume change by branch order (0 the stem): arrays order, volume_a, volume_b, trusted_change and untrusted_change. The stem's trusted part is the sum over the trusted taper bins; a branch counts under its order in the earlier model (the later one for new branches).
taper	dict[str, ndarray]	The stem profile, one entry per height bin: z0 and z1 (bin above the base of model a), radius_a and radius_b (a line through the measured radii, at the bin centre), increment (radius_b - radius_a), sigma (its standard uncertainty: the residual scatter about both lines and radius_sigma, not divided by the number of cylinders, since neighbouring fits are smoothed together), n_fits_a / n_fits_b (measured cylinders), measured_a / measured_b (measured share of the stem length), volume_a / volume_b, fitted (both models fitted the bin to points: min_fits and min_measured) and trusted (fitted, and within clip of its uncertainty from the mean increment; a fitted bin far out of line is almost always a radius the QSM regularised, such as a stem tip tapered to its apex).
taper_increment, taper_sigma, n_taper_bins		The mean radius increment of the stem over the trusted bins, weighted by 1 / sigma**2, its standard uncertainty (the weighted-mean error, scaled up by the Birge ratio when the bins scatter more than their uncertainties allow) and the number of bins used; NaN without trusted bins.
matched	dict[str, ndarray]	One entry per matched branch pair (the stem pair included): id_a, id_b, order_a, order_b, base_distance, angle, parent_consistent, length_a, length_b, volume_a, volume_b, mean_radius_a, mean_radius_b, tip_shift (distance between the tips), measured_a, measured_b, volume_sigma and trusted.
lost, new		Branches of a without a match in b, and of b without one in a: id, order, parent, base_x/y/z, length, volume, measured, status, observed_share and filled_share (NaN without a grid), volume_sigma and trusted. status is "lost" (or "new"); "unobserved" when the other epoch's grid did not observe the branch's space; "present" when that space holds returns (the branch is there but the other model lacks it).
height_a, height_b, height_trusted		Tree heights, and whether a measured cylinder reaches within top_band of the top in both models.
dbh_a, dbh_b, dbh_trusted		Stem diameter at 1.3 m, and whether its bin is trusted.
crown_area_a, crown_area_b, crown_volume_a, crown_volume_b, crown_trusted		Projected area and stacked-hull volume of the crown outlined by the branches (sylva.qsm.QSM.metrics); trusted when both models measured at least min_measured of their branch length and no unmatched branch is unobserved.
base_z	float	Height of the base both stem profiles are measured from.

change property

change: float

Total volume change, volume_b - volume_a (m³).

branches

branches(status: str) -> dict[str, np.ndarray]

Unmatched branches of one status.

Parameters:

Name	Type	Description	Default
status	str	"lost", "unobserved", "present" (branches of model a) or "new".	required

Returns:

Type	Description
dict	The columns of lost (or new) for those branches.

Raises:

Type	Description
ValueError	For an unknown status.

summary

summary() -> dict

Headline numbers with units in the keys.

Returns:

Type	Description
dict	volume_a_m3, volume_b_m3, change_m3, trusted_change_m3, untrusted_change_m3, trusted_sigma_m3, taper_increment_m, taper_sigma_m, height_change_m, dbh_change_m, n_matched, and the unmatched branches by status (n_lost, n_unobserved, n_present, n_new).

PlotQSMChange dataclass

QSM change over a plot; from compare_plot_qsms.

Attributes:

Name	Type	Description
rows	list[dict]	One dict per tree, in the order survivors, deaths, recruits (each as given): fate ("survivor", "death", "recruit"), tree_id_a, tree_id_b (None where the tree has none), volume_a_m3, volume_b_m3, change_m3, trusted_change_m3, untrusted_change_m3, trusted_sigma_m3, stem_change_m3, branch_change_m3, taper_increment_m, taper_sigma_m, height_a_m, height_b_m, n_matched, n_lost, n_unobserved, n_present, n_new and note (why a row has no comparison). A dead tree counts its whole volume as lost, a recruit its whole volume as new; the measured part is the trusted one.
changes	dict[tuple[int, int], QSMChange]	{(tree_id_a, tree_id_b): QSMChange} for the compared survivors.
totals	dict	n_survivors, n_deaths, n_recruits, n_unmodelled (rows without a model, left out of the sums), growth_m3 (survivors' change), growth_trusted_m3, growth_sigma_m3, mortality_m3 and recruitment_m3 (positive volumes), net_m3 (growth - mortality + recruitment) and net_trusted_m3.

table

table() -> list[dict]

The rows, ready for pandas.DataFrame or a CSV.

Returns:

Type	Description
list of dict	rows, copied.

to_csv

to_csv(path: str | Path) -> None

Write table as a CSV (empty cells for NaN and missing ids).

Parameters:

Name	Type	Description	Default
path	str | Path	Output file.	required

compare_qsms

compare_qsms(qsm_a: QSM, qsm_b: QSM, grid_a=None, grid_b=None, *, height_step: float = 1.0, max_base_distance: float = 0.5, max_angle: float = 35.0, parent_penalty: float = 1.0, direction_reach: float = 1.0, min_measured: float = 0.5, min_fits: int = 3, radius_sigma: float = 0.001, clip: float = 3.0, min_observed: float = 0.5, max_filled: float = 0.5, top_band: float = 1.0, crown_branch_length: float = 1.0, crown_slice: float = 0.5) -> QSMChange

Compare two models of one tree from two epochs.

Both models must be in one frame (align the epochs first). The stem is compared by height bins above the base of qsm_a: in each bin a line through each model's measured stem radii gives the radius at the bin centre, and their difference is the radius increment. Bins both models fitted and in line with the others (within clip uncertainties of the mean) are trusted, and their weighted mean is the taper increment. Branches (one branch_id chain each, as sylva.qsm.QSM.branches) are matched one order at a time, parents first: a pair is admissible when the bases are within max_base_distance and the directions (over direction_reach past the first cylinder) within max_angle, costs distance / max_base_distance + angle / max_angle, plus parent_penalty when the parents are not matched to each other, and the cheapest pairs are taken first. Unmatched branches of qsm_a are lost and of qsm_b new; with a grid of the other epoch, one whose space the grid did not observe is "unobserved" instead.

Parameters:

Name	Type	Description	Default
qsm_a	QSM	The earlier and the later model (sylva.qsm.QSM, or its (n, 12) cylinder array).	required
qsm_b	QSM	The earlier and the later model (sylva.qsm.QSM, or its (n, 12) cylinder array).	required
grid_a		Optional sylva.voxels.RayVoxelGrid of each epoch, from sylva.voxels.ray_voxelize over that epoch's pulses (misses included; occlusion=True separates occluded space). grid_b checks lost branches and grid_a new ones: each is sampled every half voxel along the branch's cylinder axes, leaving out the part inside its parent.	None
grid_b		Optional sylva.voxels.RayVoxelGrid of each epoch, from sylva.voxels.ray_voxelize over that epoch's pulses (misses included; occlusion=True separates occluded space). grid_b checks lost branches and grid_a new ones: each is sampled every half voxel along the branch's cylinder axes, leaving out the part inside its parent.	None
height_step	float	Height of the stem bins (m).	1.0
max_base_distance	float	Branch matching, as above (m, degrees, cost, m).	0.5
max_angle	float	Branch matching, as above (m, degrees, cost, m).	0.5
parent_penalty	float	Branch matching, as above (m, degrees, cost, m).	0.5
direction_reach	float	Branch matching, as above (m, degrees, cost, m).	0.5
min_measured	float	Share of a stem bin's length, or of a branch's length, that must be measured (fitted to points) in a model for that part to be trusted.	0.5
min_fits	int	Measured stem cylinders a bin needs in each model to be trusted.	3
radius_sigma	float	Radius uncertainty of one fit (m), added in quadrature to the scatter-based uncertainties. The scatter about a line understates it, since the QSM smooths radii along the stem; 1 mm suits circle fits on TLS stems.	0.001
clip	float	Bins further than this many of their own uncertainties from the mean increment are not trusted (see QSMChange.taper).	3.0
min_observed	float	Share of a branch's samples the grid must have observed (crossed by a pulse or holding a return); below it the branch is unobserved.	0.5
max_filled	float	Share of the observed samples holding returns above which the branch is "present" rather than lost or new.	0.5
top_band	float	A model's height is trusted when a measured cylinder reaches within this distance of its top (m).	1.0
crown_branch_length	float	Crown base and slice height, as in sylva.qsm.QSM.metrics.	1.0
crown_slice	float	Crown base and slice height, as in sylva.qsm.QSM.metrics.	1.0

Returns:

Type	Description
QSMChange

Raises:

Type	Description
ValueError	For a model that is not a QSM or (n, 12) array, non-finite cylinders, a grid that is not a ray-traced grid, or a setting out of range.

Notes

On a synthetic tree scanned from three positions in each epoch (the change detection guide), a stem thickened by 10 mm was recovered as a taper increment of 10.1 to 10.3 mm with a stated uncertainty of 0.6 mm over five noise draws; a cut limb was lost, a limb extended by 1 m was matched 0.95 to 1.00 m longer, and a limb hidden by foliage in the later epoch was unobserved in every draw.

A cylinder is measured when its n_points is positive; the QSM sets it to 0 where it replaced the cylinder's own circle by its allometric prior, so a prior radius is never read as trusted change.

Examples:

>>> grid_b = voxels.ray_voxelize(shots_b, 0.1, occlusion=True)
>>> c = change.compare_qsms(model_2020, model_2025, grid_b=grid_b)
>>> c.taper_increment, c.taper_sigma
>>> c.branches("lost")["volume"].sum()

compare_plot_qsms

compare_plot_qsms(qsms_a, qsms_b, pairs: Iterable, deaths: Iterable = (), recruits: Iterable = (), grid_a=None, grid_b=None, **settings) -> PlotQSMChange

Compare the QSMs of every tree of a plot between two epochs.

Survivors are compared with compare_qsms; a dead tree counts its whole earlier volume as mortality and a recruit its whole later volume as recruitment. Trees run in parallel; the result does not depend on the number of threads.

Parameters:

Name	Type	Description	Default
qsms_a		Models of each epoch: sylva.qsm.PlotQSMs or {tree_id: QSM}. A tree without a model gives a row with a note and no change.	required
qsms_b		Models of each epoch: sylva.qsm.PlotQSMs or {tree_id: QSM}. A tree without a model gives a row with a note and no change.	required
pairs	Iterable	Survivors as (tree_id_a, tree_id_b) pairs, e.g. from a tree match between the epochs.	required
deaths	Iterable	Tree ids of epoch a that died.	()
recruits	Iterable	Tree ids of epoch b that are new.	()
grid_a		Optional ray-traced grids of the whole plot in each epoch, as for compare_qsms.	None
grid_b		Optional ray-traced grids of the whole plot in each epoch, as for compare_qsms.	None
**settings		Keyword settings of compare_qsms.	{}

Returns:

Type	Description
PlotQSMChange

Raises:

Type	Description
ValueError	For a tree id used twice in one epoch, a malformed pair, an unknown setting, or anything compare_qsms rejects.

Examples:

>>> plot = change.compare_plot_qsms(qsms_2020, qsms_2025, pairs=[(1, 4), (2, 5)],
...                                 deaths=[3], recruits=[9], grid_b=grid_2025)
>>> plot.totals["net_trusted_m3"]; plot.to_csv("qsm_change.csv")