sylva.als: canopy structure

Re-exported by sylva.als; see the guide.

sylva.als.canopy

Ray-based canopy structure from airborne and UAV lidar.

The functions here are re-exported by sylva.als, and take either one ~sylva.PointCloud or a whole ~sylva.als.Catalog, which they process chunk by chunk with the engine of sylva.als:

1. Trajectory. read_trajectory reads an SBET or a text table of time, position and attitude; Trajectory.positions interpolates the sensor position at any GPS time. Without one, estimate_trajectory recovers an approximate trajectory from the pulses with several returns.
2. Pulses. pulses groups returns into pulses by GPS time and gives ~sylva.Shots whose origin is the sensor, so that sylva.voxels.ray_voxelize traces airborne pulses as it traces terrestrial ones.
3. Canopy. gap_profile gives gap-fraction profiles, plant area density and gridded plant area index by the method of MacArthur & Horn (1969), corrected for each return's beam angle; ray_voxelize traces the pulses of a whole catalogue through one voxel lattice.

Examples:

>>> from sylva import als
>>> traj = als.read_trajectory("flight.out", crs="EPSG:32755")
>>> prof = als.gap_profile("tiles/", traj, resolution=20.0)
>>> height, pad = prof.profile()
>>> vox = als.ray_voxelize("tiles/", traj, voxel_size=2.0)
>>> vox.pai(dtm=dtm).to_geotiff("pai.tif")

ALSProfile dataclass

Gap-fraction counts of an area, per grid cell and height layer.

From gap_profile. Layer 0 holds the returns below min_height (the pulses that got through the canopy), layers 1 to nz the bins of bin_size above it, and the last layer the returns above the top bin. weight sums each return's share of its pulse, weight_k that share times its extinction per unit vertical plant area, G(θ) / cos θ at its beam zenith θ. Because these are counts, cells and chunks add up exactly; the densities are derived from them. It is plain data and pickles.

Parameters:

Name	Type	Description	Default
xmin	float	South-west corner of the grid.	required
ymin	float	South-west corner of the grid.	required
resolution	float	Cell size (m).	required
min_height	float	Bottom of the first layer and layer thickness (m).	required
bin_size	float	Bottom of the first layer and layer thickness (m).	required
weight	ndarray	(nz + 2, ny, nx), row 0 at the south.	required
weight_k	ndarray	(nz + 2, ny, nx), row 0 at the south.	required
crs	str | None	CRS of the grid, if known.	None
n_skipped	int	Returns left out (no trajectory position, or beyond max_zenith).	0

shape property

shape: tuple[int, int, int]

(nz, ny, nx): layers and cells.

heights property

heights: ndarray

Bottom of each layer (m above ground).

pad

pad() -> np.ndarray

Plant area density of every layer and cell (m²/m³).

Layer i lets through the weight below it out of the weight at or below it, T; its density is -ln T / (k̄ dz), with k̄ the mean extinction of the returns at or below it (MacArthur & Horn 1969, with the beam angle of each return). A layer nothing passed gets T = 0.5 / max(E, 1) (half a pulse of the E that entered), a lower bound (not NaN).

Returns:

Type	Description
ndarray	(nz, ny, nx); NaN where no pulse reached a layer.

pai

pai() -> Raster

Plant area index of each cell above min_height.

-ln P / k̄, with P the share of the weight below min_height (the gap probability) and k̄ the mean extinction of all returns; saturated cells (no return below) as in pad. With one beam angle this is the sum of the layer densities times bin_size.

Returns:

Type	Description
Raster	NaN for cells without returns.

cover

cover() -> Raster

Canopy cover of each cell: 1 - P at min_height.

Returns:

Type	Description
Raster

profile

profile(mask=None) -> tuple[np.ndarray, np.ndarray]

Plant area density profile of an area, from its pooled counts.

The counts of the cells are added before the densities are derived, so this is the profile of the area seen as one column: an effective profile, lower than the mean of the cells' profiles where the canopy is clumped between cells.

Parameters:

Name	Type	Description	Default
mask		(ny, nx) booleans selecting the cells; all if None.	None

Returns:

Name	Type	Description
height	ndarray	Bottom of each layer (m).
pad	ndarray	Plant area density (m²/m³); NaN above the highest return.

pgap

pgap(mask=None) -> tuple[np.ndarray, np.ndarray]

Gap probability profile of an area: the share of the pooled weight below each layer boundary.

Parameters:

Name	Type	Description	Default
mask		(ny, nx) booleans selecting the cells; all if None.	None

Returns:

Name	Type	Description
height	ndarray	min_height and the top of every layer (m).
pgap	ndarray	Gap probability at each height (not corrected for beam angle).

pooled_pai

pooled_pai(mask=None) -> float

Plant area index of an area from its pooled counts (see profile).

Parameters:

Name	Type	Description	Default
mask		(ny, nx) booleans selecting the cells; all if None.	None

Returns:

Type	Description
float

metrics

metrics(strata: float = 5.0) -> dict[str, Raster]

Summary metrics of each cell's profile, as rasters, for use beside the height metrics of sylva.als.grid_metrics in a model of height, cover or biomass.

• pulses: the cell's weight (its number of pulses with the "equal" weighting); a check on how well the cell is sampled.
• pai and cover: as pai and cover.
• fhd: foliage height diversity, the Shannon index -sum(p ln p) of the shares p of the plant area in each layer (MacArthur & MacArthur 1961; GEDI L2B fhd_normal). 0 where there is no plant area.
• pad_max and height_pad_max: the density of the densest layer and the height of its middle.
• height_pad_mean and height_pad_sd: the mean and standard deviation of height weighted by plant area: where the canopy is and how deep it is.
• For each stratum [a, b) of strata m from the ground (the strata of GEDI L2B with strata=5): pavd_<a>_<b>, its mean plant area density; pai_above_<h>, the plant area index above h, the bottom of the stratum's lowest layer; cover_above_<h>, the canopy cover at h.

A layer belongs to the stratum its middle falls in, and strata below min_height are left out (so with min_height=1 the lowest stratum is pavd_0_5 over 1 to 5 m). Layers no pulse reached are left out of every sum. pai corrects each return for its own beam angle; the layer sums behind pai_above_<h> correct each layer by the mean angle of the returns that reached it, so the two can differ slightly under a wide swath.

Parameters:

Name	Type	Description	Default
strata	float	Thickness of the strata (m); at least bin_size.	5.0

Returns:

Type	Description
dict[str, Raster]	One raster per metric, in a fixed order. Cells without returns are NaN, but for pulses (0).

area_metrics

area_metrics(mask=None, strata: float = 5.0) -> dict[str, float]

The metrics of metrics for an area, from its pooled counts (the area seen as one column, as profile).

Parameters:

Name	Type	Description	Default
mask		(ny, nx) booleans selecting the cells; all if None.	None
strata	float	Thickness of the strata (m).	5.0

Returns:

Type	Description
dict[str, float]

plot_cells

plot_cells(plots, radius=None) -> list[np.ndarray]

The cells of each plot: those whose centre lies inside it.

Parameters:

Name	Type	Description	Default
plots		As for sylva.als.plot_metrics: (N, 2) centres with a radius, or polygons in any form sylva.geo.masks accepts.	required
radius		As for sylva.als.plot_metrics: (N, 2) centres with a radius, or polygons in any form sylva.geo.masks accepts.	required

Returns:

Type	Description
list of numpy.ndarray	Row-major cell indices (row * nx + column) of each plot.

plot_metrics

plot_metrics(plots, radius=None, ids=None, strata: float = 5.0)

The metrics of metrics for field plots, each from the pooled counts of the cells whose centre lies inside it, as a table to join with sylva.als.plot_metrics and the plot data.

The plot is represented by whole cells, so use a resolution well below the plot size (a 2 m grid for a 0.1 ha plot). Because the counts add up exactly, pooling small cells gives the profile of the plot itself. area gives the area of the cells used.

Parameters:

Name	Type	Description	Default
plots		As for sylva.als.plot_metrics.	required
radius		As for sylva.als.plot_metrics.	required
ids		Optional identifiers, one per plot.	None
strata	float	Thickness of the strata (m).	5.0

Returns:

Type	Description
PlotMetrics	Columns plot, id (with ids), area (m²), then the metrics. A plot covering no cell centre has area 0 and pulses 0, NaN elsewhere.

ALSVoxels dataclass

Ray-traced voxel fields of an airborne survey, from ray_voxelize.

The lattice is anchored at multiples of voxel_size; every array is (nz, ny, nx) with row 0 at the south and layer 0 at the bottom. Columns no chunk traced are NaN. It is plain data and pickles.

Parameters:

Name	Type	Description	Default
origin	ndarray	(x, y, z) of the lattice's lower corner.	required
voxel_size	float	Voxel edge (m).	required
fields	dict[str, ndarray]	The fields and metrics asked for, by name (see sylva.voxels.RayVoxelGrid); num_beams is always there.	required
reach	float	The largest horizontal distance (m) a pulse covered between the top of the grid and its last return: the buffer that brings every pulse crossing a chunk into it.	0.0
n_pulses	int	Pulses traced (a pulse within reach of several chunks counts in each).	0
crs	str | None	CRS, if known.	None

shape property

shape: tuple[int, int, int]

(nz, ny, nx).

z_levels property

z_levels: ndarray

Bottom z of each layer.

profile

profile(name: str | None = None, dtm: Raster | None = None, bin_size: float | None = None, min_beams: float = 1.0, mask=None) -> tuple[np.ndarray, np.ndarray]

Mean of a field per height above ground.

Parameters:

Name	Type	Description	Default
name	str | None	Field; the first pad_ field if None.	None
dtm	Raster | None	Terrain: heights are those of voxel centres above it. Without one, above the grid floor.	None
bin_size	float | None	Height bin (m); the voxel size if None.	None
min_beams	float	Voxels crossed by fewer pulses are left out.	1.0
mask		(ny, nx) booleans selecting the columns; all if None.	None

Returns:

Name	Type	Description
height	ndarray	Bottom of each bin, from 0.
mean	ndarray	NaN for bins without voxels.

pai

pai(name: str | None = None, dtm: Raster | None = None, min_height: float = 0.0, min_beams: float = 1.0) -> Raster

Plant area index of each column: the sum of a density field times the voxel size.

Parameters:

Name	Type	Description	Default
name	str | None	Density field; the first pad_ field if None.	None
dtm	Raster | None	Terrain, for min_height.	None
min_height	float	Voxels whose centre is lower above the DTM are left out (the ground itself, understorey).	0.0
min_beams	float	Voxels crossed by fewer pulses count as 0.	1.0

Returns:

Type	Description
Raster	(ny, nx) on the voxel columns; NaN where no voxel was seen.

gap_profile

gap_profile(source, trajectory=None, resolution: float = 10.0, bin_size: float = 1.0, min_height: float = 1.0, max_height: float | None = None, top_quantile: float = 0.99999, drop_noise: bool = True, weighting: str = 'equal', angles: str = 'auto', lad: str = 'spherical', lad_params: Sequence[float] = (), g: float | None = None, max_zenith: float = 90.0, anchor: str = 'ground', dtm='auto', dtm_resolution: float = 1.0, bounds=None, max_gap: float | None = None, time_offset: float = 0.0, chunk_size: float | None = None, buffer: float = 20.0, workers: int | None = None) -> ALSProfile

Gap-fraction profiles, plant area density and plant area index.

The returns of each grid cell are binned by height above ground. A pulse passes a height if it has weight below it, so the gap probability at height z is the share of the cell's weighted returns below z (Armston et al. 2013), and the transmittance of a layer the weight below it over the weight at or below it. Inverting Beer-Lambert layer by layer gives the plant area density profile of MacArthur & Horn (1969), as computed by Bouvier et al. (2015); that is weighting="all", angles="none", g=0.5 here.

Scan angle. A beam at zenith θ crosses dz / cos θ of a layer of thickness dz and meets foliage in proportion to G(θ), the projection of the leaf angle distribution, so each return carries the extinction k = G(θ) / cos θ per unit of vertical plant area, and a layer's density divides -ln T by the mean k of the returns that reached it. Without this, returns from the edge of a 30° swath read 15 % denser than at nadir. θ comes from the trajectory (the angle of the line from each return to the sensor), else from the LAS scan_angle (which includes roll but not pitch).

Which cell. An oblique pulse stopped in the canopy leaves its return to one side of where it would have reached the ground. Counted where they are, the returns of a cell come from pulses picked partly by their outcome, which biases the gap probability wherever the pulse density varies across the swath (by 5 % on the synthetic layer of the guide). With a trajectory each return is therefore counted, by default, in the cell where its beam meets the ground (anchor="ground"): the pulses of a cell are then chosen by their geometry alone. Without a trajectory the returns stay where they are.

Every return is counted once, so a catalogue gives the same counts, cell for cell, whatever the chunks, as long as the buffer is wider than the drift of a beam between its returns and the ground (canopy height times the tangent of the scan angle); the buffer also serves the "auto" DTM.

Parameters:

Name	Type	Description	Default
source		A ~sylva.PointCloud, or a catalogue (anything sylva.als.catalog accepts).	required
trajectory		Trajectory (or mapping, or file) for the beam angles.	None
resolution	float	Cell size (m). Cells should hold a few hundred pulses; 10-30 m is usual.	10.0
bin_size	float	Layer thickness (m).	1.0
min_height	float	Bottom of the lowest layer (m): returns below it, ground included, are pulses that went through.	1.0
max_height	float | None	Top of the highest layer; returns above it are intercepted above every layer. By default the layers stop where top_quantile of the weight above min_height lies below.	None
top_quantile	float	With max_height None, the top layers that together hold no more than 1 - top_quantile of the weight above min_height are dropped (their returns count as above every layer), so that a few stray returns far above the canopy (birds, haze, a mast) do not stretch the profile with empty layers. The layers kept are as they would be without it, since a layer's transmittance depends only on the returns at and below it. 1 reaches the highest return.	0.99999
drop_noise	bool	Leave out returns classified as noise (7 or 18).	True
weighting	('equal', 'first', 'all')	Share of a pulse each return stands for: 1 / number_of_returns, first returns only (the original MacArthur-Horn), or one per return (Bouvier et al. 2015).	"equal"
angles	('auto', 'trajectory', 'scan_angle', 'none')	Source of the beam zenith; "auto" uses the trajectory if given, else scan_angle. "none" treats every beam as vertical.	"auto"
lad	str	Leaf angle distribution for G(θ) (see sylva.voxels.leaf_projection).	'spherical'
lad_params	str	Leaf angle distribution for G(θ) (see sylva.voxels.leaf_projection).	'spherical'
g	float | None	A constant G instead of lad (0.5 in Bouvier et al. 2015).	None
max_zenith	float	Returns whose beam is further than this from nadir (degrees) are left out.	90.0
anchor	('ground', 'return')	Count each return in the cell where its beam meets the ground (needs the trajectory; see above), or where the return is.	"ground"
dtm	('auto', None or Raster)	Heights: "auto" makes a DTM from the ground points (classification 2; per chunk with its buffer for a catalogue), a ~sylva.Raster is subtracted, None uses z as it is (normalised tiles).	'auto'
dtm_resolution	float	Cell size (m) of the "auto" DTM.	1.0
bounds		(xmin, ymin, xmax, ymax) of the grid for a point cloud; its extent by default. A catalogue uses the catalogue grid.	None
max_gap	float | None	As for pulses.	None
time_offset	float | None	As for pulses.	None
chunk_size	float | None	As for sylva.als.apply.	None
buffer	float | None	As for sylva.als.apply.	None
workers	float | None	As for sylva.als.apply.	None

Returns:

Type	Description
ALSProfile

Raises:

Type	Description
ValueError	For bad settings, a weighting whose attributes are missing, or angles that cannot be had.

ray_voxelize

ray_voxelize(source, trajectory, voxel_size: float = 1.0, fields: Sequence[str] | None = None, *, dtm: Raster | None = None, ground_class: int | None = 2, ground_distance: float = 0.2, class_attr: str = 'classification', leaf_classes: Sequence[int] = (), wood_classes: Sequence[int] = (), tree_attr: str = '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] = (), average_leaf_area: float = 0.005, occlusion: bool = False, unbounded_range: float = np.inf, bounds=None, z_range=None, max_gap: float | None = None, time_offset: float = 0.0, fill_missing: bool = False, max_fill: int = 8, drop_incomplete: bool = False, chunk_size: float | None = None, buffer: float = 20.0, workers: int | None = None) -> ALSVoxels

Trace the pulses of a survey through one voxel lattice.

The returns are grouped into pulses from the sensor as pulses does, and traced by the ray tracer of sylva.voxels.ray_voxelize along their actual beam directions, so every voxel's path lengths, hits and mean beam angle, and so its plant area density, account for the scan angle.

On a catalogue, each chunk reconstructs the pulses of its core and buffer points and traces them through the columns whose centres lie in its core; the columns are put together into one lattice anchored at the catalogue's minimum corner snapped down to a multiple of voxel_size. A pulse that crosses a chunk's core from outside it is traced there as long as one of its returns lies in the buffer. That needs a buffer at least the height of the grid above the pulse's last return times the tangent of its zenith (for a 30 m canopy and a 30° swath, 17 m); ALSVoxels.reach reports the largest such distance among the pulses traced, and a warning is given when it exceeds the buffer. With a sufficient buffer the result does not depend on the chunks, up to floating-point rounding.

Use z as measured, not heights above ground: normalising moves each return by a different amount and bends the rays.

Parameters:

Name	Type	Description	Default
source		A ~sylva.PointCloud, or a catalogue.	required
trajectory		Trajectory (or mapping, or file).	required
voxel_size	float	Voxel edge (m). Memory is about 0.4 kB per voxel of a chunk while tracing, and 8 bytes per voxel and field of the result.	1.0
fields	Sequence[str] | None	Fields and metrics to keep (see sylva.voxels.RayVoxelGrid); by default pad_<attenuation>, transmittance, num_hits_weighted and free_path_length. num_beams is always added.	None
dtm	Raster | None	Echo labels, as for sylva.voxels.ray_voxelize; by default returns of class 2 are ground (traced up to, never hits).	None
ground_class	Raster | None	Echo labels, as for sylva.voxels.ray_voxelize; by default returns of class 2 are ground (traced up to, never hits).	None
ground_distance	Raster | None	Echo labels, as for sylva.voxels.ray_voxelize; by default returns of class 2 are ground (traced up to, never hits).	None
class_attr	Raster | None	Echo labels, as for sylva.voxels.ray_voxelize; by default returns of class 2 are ground (traced up to, never hits).	None
leaf_classes	Sequence[int]	As for sylva.voxels.ray_voxelize.	()
wood_classes	Sequence[int]	As for sylva.voxels.ray_voxelize.	()
tree_attr	Sequence[int]	As for sylva.voxels.ray_voxelize.	()
intensity_attr	Sequence[int]	As for sylva.voxels.ray_voxelize.	()
weighting	str	As for sylva.voxels.ray_voxelize.	'equal'
attenuation	str	As for sylva.voxels.ray_voxelize.	'equal'
laser	str	As for sylva.voxels.ray_voxelize.	'equal'
beam	str	As for sylva.voxels.ray_voxelize.	'equal'
lad	str	As for sylva.voxels.ray_voxelize.	'equal'
lad_params	str	As for sylva.voxels.ray_voxelize.	'equal'
average_leaf_area	float	As for sylva.voxels.ray_voxelize.	0.005
occlusion	float	As for sylva.voxels.ray_voxelize.	0.005
unbounded_range	float	As for sylva.voxels.ray_voxelize.	0.005
bounds		(xmin, ymin, zmin, xmax, ymax, zmax) for a point cloud; its extent by default.	None
z_range		(zmin, zmax) of a catalogue's lattice; the header z range by default.	None
max_gap	float | None	As for pulses.	None
time_offset	float | None	As for pulses.	None
fill_missing	float | None	As for pulses.	None
max_fill	float | None	As for pulses.	None
drop_incomplete	float | None	As for pulses.	None
chunk_size	float | None	As for sylva.als.apply. Chunks are sized to memory with their voxels counted.	None
buffer	float | None	As for sylva.als.apply. Chunks are sized to memory with their voxels counted.	None
workers	float | None	As for sylva.als.apply. Chunks are sized to memory with their voxels counted.	None

Returns:

Type	Description
ALSVoxels

Raises:

Type	Description
ValueError	For bad settings, a lattice too large for the memory budget, or returns outside the trajectory.

Warns:

Type	Description
UserWarning	When a catalogue's buffer is narrower than the pulses' reach.

sylva.als.trajectory

The flight trajectory of an airborne survey, and the pulses recovered from it.

Trajectory dataclass

The path of the scanner: positions, and optionally attitude, in time.

Build one with read_trajectory or estimate_trajectory, or from arrays. Samples are sorted by time and a sample repeating an earlier time is dropped. It is plain data and pickles.

Parameters:

Name	Type	Description	Default
time	ndarray	GPS time of each sample (s), on the clock of the returns' gps_time (see time_offset of pulses otherwise).	required
x	ndarray	Sensor position (m), in the frame of the points.	required
y	ndarray	Sensor position (m), in the frame of the points.	required
z	ndarray	Sensor position (m), in the frame of the points.	required
roll	ndarray | None	Attitude (degrees; heading clockwise from north), if known.	None
pitch	ndarray | None	Attitude (degrees; heading clockwise from north), if known.	None
heading	ndarray | None	Attitude (degrees; heading clockwise from north), if known.	None
rms	ndarray | None	For an estimated trajectory: the root mean square distance (m) of the pulse lines from each sample, the pulses used, and the flight line.	None
n_pulses	ndarray | None	For an estimated trajectory: the root mean square distance (m) of the pulse lines from each sample, the pulses used, and the flight line.	None
line	ndarray | None	For an estimated trajectory: the root mean square distance (m) of the pulse lines from each sample, the pulses used, and the flight line.	None

Raises:

Type	Description
ValueError	Fewer than two samples at different times, arrays of different lengths, or non-finite values.

xyz property

xyz: ndarray

(N, 3) positions.

has_attitude property

has_attitude: bool

Whether roll, pitch and heading are known.

from_dict classmethod

from_dict(d: Mapping) -> Trajectory

From a mapping of arrays, such as sylva.synthetic.ALSFlight.trajectory.

Parameters:

Name	Type	Description	Default
d	Mapping	time, x, y, z and optionally roll, pitch, heading, rms, n_pulses, line.	required

Returns:

Type	Description
Trajectory

positions

positions(gps_time, max_gap: float | None = None, time_offset: float = 0.0) -> np.ndarray

Sensor position at each time, linearly interpolated.

Parameters:

Name	Type	Description	Default
gps_time		Times, e.g. cloud.attrs["gps_time"].	required
max_gap	float | None	Longest gap between samples (s) to interpolate across; ten median sample intervals if None.	None
time_offset	float	Seconds added to gps_time first.	0.0

Returns:

Type	Description
ndarray	(N, 3); NaN outside the trajectory or within a longer gap.

attitude

attitude(gps_time, max_gap: float | None = None, time_offset: float = 0.0) -> np.ndarray

Roll, pitch and heading at each time, linearly interpolated.

Parameters:

Name	Type	Description	Default
gps_time		Times.	required
max_gap	float | None	As for positions.	None
time_offset	float | None	As for positions.	None

Returns:

Type	Description
ndarray	(N, 3) degrees; the heading is interpolated along the shorter arc and wrapped to [0, 360).

Raises:

Type	Description
ValueError	If the trajectory has no attitude.

week_seconds

week_seconds(gps_time)

GPS seconds of the week of adjusted standard GPS times.

LAS files flagged for it store adjusted standard GPS time (GPS time minus 10⁹ s); SBET trajectories count seconds from the start of the GPS week. Within one week the two differ by a constant, which week_seconds(t) - t gives, for the time_offset arguments.

Parameters:

Name	Type	Description	Default
gps_time		Adjusted standard GPS times (scalar or array).	required

Returns:

Type	Description
float or ndarray	(t + 1e9) mod 604800.

read_trajectory

read_trajectory(path: str | Path, format: str | None = None, crs: str | None = None, columns: Sequence[str] | None = None, z_offset: float = 0.0) -> Trajectory

Read a trajectory file: an SBET or a delimited text table.

SBET (format="sbet", or a .out, .sbet, .sbt or .bin file): the binary "smoothed best estimate of trajectory" of Applanix POSPac and most GNSS/INS software, records of 17 little-endian doubles (time in GPS seconds of the week, latitude, longitude and platform heading in radians, ellipsoidal height, ...). Latitude and longitude are projected to crs, which is required; the heading is the platform heading minus the wander angle.

Text (any other file): columns separated by commas, semicolons or white space, with a header line naming them or with columns given. Names are matched case-insensitively: time (or gps_time, t, timestamp), x (easting), y (northing), z (height, altitude, elevation), roll, pitch and heading (yaw); other columns are ignored. With crs, x and y are longitude and latitude in degrees and are projected.

Parameters:

Name	Type	Description	Default
path	str | Path	The file.	required
format	(None, 'sbet', 'text')	None to decide from the suffix.	None
crs	str | None	The CRS of the points (e.g. "EPSG:32755"), for geographic trajectories.	None
columns	Sequence[str] | None	Names of the text columns, in order, instead of a header line.	None
z_offset	float	Added to every height. SBET heights are ellipsoidal: for points in orthometric heights, give minus the geoid separation. A wrong height shows as a large line_offset_median in the report of pulses.	0.0

Returns:

Type	Description
Trajectory

Raises:

Type	Description
OSError	If the file cannot be read or is not a whole number of SBET records.
ValueError	Without crs for an SBET, for a projection that is not in metres, or for a text table without time, x, y and z columns.

estimate_trajectory

estimate_trajectory(source, interval: float = 0.5, min_pulses: int = 30, min_separation: float = 2.0, max_pulses: int = 4000, extend: float = 2.0, workers: int | None = None) -> Trajectory

An approximate trajectory from the returns themselves.

Every pulse with two or more returns defines a line through the scanner: from its last return through its first. Within a short time window the platform flies a nearly straight line at nearly constant speed, so its position and velocity there follow from a linear least squares fit of a moving point a + b (t - t_c) to those lines. This is the idea of Roussel et al. (2020), after Gatziolis & McGaughey (2019), which fits a fixed point per window; the motion within the window is modelled here, so longer windows can be used without the position lagging.

Windows of interval seconds are anchored at multiples of interval in GPS time, per flight line (point_source_id). Each window with min_pulses usable pulses (first and last returns at least min_separation apart, the max_pulses widest used) gives a sample at the mean time of its pulses; the first and last windows of a line also give samples towards the line's first and last returns, extrapolated with the fitted velocity by at most extend seconds (over open ground a line has no multiple returns to fit). Lines more than three median distances (and 5 cm) from the first fit are left out of a second one, and a window whose sensor comes out below its returns is rejected.

The result is approximate: it depends on the multiple returns (few over open ground or sparse canopy), carries no attitude, and cannot follow the platform's wobble within a window. Check rms and use a measured trajectory where one exists.

Parameters:

Name	Type	Description	Default
source		A ~sylva.PointCloud with gps_time (and return_number, point_source_id if present), or a catalogue, whose tiles are read one at a time (each tile's own points only).	required
interval	float	Window length (s).	0.5
min_pulses	int	Fewest usable pulses in a window.	30
min_separation	float	Least first-to-last distance (m) of a usable pulse.	2.0
max_pulses	int	Most pulses fitted per window.	4000
extend	float	Longest extrapolation (s) at each end of a line.	2.0
workers	int | None	Tiles read at once for a catalogue.	None

Returns:

Type	Description
Trajectory	With rms (m), n_pulses and line per sample.

Raises:

Type	Description
ValueError	For bad settings, points without gps_time, or too few multiple returns for any window.

pulses

pulses(cloud: PointCloud, trajectory, max_gap: float | None = None, time_offset: float = 0.0, fill_missing: bool = False, max_fill: int = 8, drop_incomplete: bool = False, report: bool = False)

Group discrete returns into pulses fired from the sensor.

Returns sharing a gps_time (and point_source_id, when present) are one pulse; a group that repeats a return_number (two channels firing at the same instant) is split where it repeats. Each pulse starts at the sensor position interpolated from trajectory at its time, and points at its farthest return; each return's range is its distance from the sensor. The result traces with sylva.voxels.ray_voxelize like a terrestrial scan.

What discrete returns cannot give back:

• Returns missing from a pulse (a return below the detection threshold, or clipped away at a tile or chunk edge) are counted from number_of_returns; the pulse keeps the returns it has (drop_incomplete leaves it out). If its last return is missing, the tracer ends the pulse early.
• Pulses with no return (into water, onto a dark roof, lost to absorption) leave no record. With fill_missing, a hole of at most max_fill pulses in a line's regular firing is filled with pulses without returns, their directions interpolated between the pulses on either side, provided the mirror swept steadily across the hole. Larger holes, and pulses at the edge of the data, are not recovered.
• Returns closer than the range resolution were merged by the receiver, and the energy each return carried is not recorded, so every return stands for an equal share of its pulse.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	Returns with gps_time (return_number, number_of_returns and point_source_id are used when present). All attributes become echo attributes.	required
trajectory		A Trajectory, a mapping of arrays (such as sylva.synthetic.ALSFlight.trajectory) or a trajectory file.	required
max_gap	float | None	Longest gap between trajectory samples (s) to interpolate across; ten median sample intervals if None. Returns in longer gaps, or outside the trajectory, are left out.	None
time_offset	float	Seconds added to gps_time to reach the trajectory's clock (see week_seconds).	0.0
fill_missing	bool	Infer pulses without a return in holes of at most max_fill pulses.	False
max_fill	bool	Infer pulses without a return in holes of at most max_fill pulses.	False
drop_incomplete	bool	Leave out pulses with fewer returns than number_of_returns.	False
report	bool	Also return what was found.	False

Returns:

Type

Description

Shots or (Shots, dict)

Pulses ordered by flight line and time. The report has n_returns, n_pulses, n_incomplete, n_missing_returns, n_dropped, n_split, n_unpositioned (returns left out), n_filled, pulse_interval (s) and line_offset_median / line_offset_p95: the distance (m) from the sensor to the line through each pulse's first and last returns. Centimetres to decimetres mean the trajectory fits the points; metres mean a time or height offset.

Raises:

Type	Description
ValueError	Without gps_time, or when no return falls within the trajectory (the message gives the offset if the clocks look like adjusted standard time against seconds of the week).