sylva.waveform¶

sylva.waveform ¶

Full-waveform lidar: reading, Gaussian decomposition and conversion to pulses.

A full-waveform scanner records the received power of each pulse as a run of samples instead of (or besides) a few discrete returns. This module reads such records from LAS 1.3 / 1.4 files with waveform data packets (point formats 4, 5, 9 and 10, packets in the file or in a .wdp file next to it) and from PulseWaves (.pls with its .wvs), decomposes them into Gaussian echoes (Hofton et al. 2000; Wagner et al. 2006), and turns the echoes into sylva.shots.Shots for the ray-traced voxel and canopy tools. RIEGL SDF files need RIEGL's proprietary library and are not read; export them to LAS with waveform data packets first.

Waveform `dataclass` ¶

One digitised waveform: a view of a row of Waveforms.

Attributes:

Name	Type	Description
`samples`	`ndarray`	Received power (digitiser units, or physical units after a gain, offset or lookup table from the file).
`interval`	`float`	Sampling interval (ns).
`anchor, direction`		A point on the beam and the unit beam direction (away from the scanner).
`offset`	`float`	Time from the anchor to the first sample (ns).
`metres_per_ns`	`float`	Range per nanosecond of round-trip time (m/ns).
`origin`	`ndarray`	Scanner position, NaN when the file does not give it.
`gps_time, pulse`		Time of the pulse and its identifier.

times ¶

times() -> np.ndarray

Time of each sample after the first (ns).

Returns:

Type	Description
`ndarray`	`arange(n) * interval`.

positions ¶

positions() -> np.ndarray

Position of each sample.

Returns:

Type	Description
`ndarray`	`(n, 3)`: `anchor + direction * metres_per_ns * (offset + t)`.

Waveforms `dataclass` ¶

A set of waveforms, their geometry and their samples.

The samples of waveform i are samples[sample_start[i] : sample_start[i] + sample_count[i]]; sample k of it lies at anchor[i] + direction[i] * metres_per_ns[i] * (offset[i] + k * interval[i]). Rows that belong to one pulse (segments of a PulseWaves sampling, or its outgoing and returning waveforms) are consecutive and share pulse.

Build one with read, chunks or sylva.synthetic.waveforms.

Attributes:

Name	Type	Description
`pulse`	`ndarray`	Pulse identifier: the index of the pulse's first point record (LAS) or of the pulse (PulseWaves).
`gps_time`	`ndarray`	Time of the pulse (s, as stored in the file).
`origin`	`ndarray`	`(n, 3)` scanner positions; NaN where unknown (LAS files never give them).
`anchor, direction`		`(n, 3)` point on the beam and unit direction away from the scanner.
`offset`	`ndarray`	Time from the anchor to the first sample (ns).
`interval`	`ndarray`	Sampling interval (ns).
`metres_per_ns`	`ndarray`	Range per nanosecond of round-trip time (m/ns), about `C_HALF`.
`sample_start, sample_count`		CSR offsets into `samples`.
`samples`	`ndarray`	All samples, float32.
`attrs`	`dict[str, ndarray]`	Per-waveform attributes: from LAS `intensity`, `return_number`, `number_of_returns`, `classification`, `point_source_id`, `n_records` (point records sharing the packet) and `descriptor`; from PulseWaves `kind` (1 outgoing, 2 returning), `channel`, `segment`, `sampling`, `intensity`, `classification` and `descriptor`.

n_waveforms `property` ¶

n_waveforms: int

Number of waveforms (rows).

n_samples `property` ¶

n_samples: int

Number of samples over all waveforms.

sample_positions ¶

sample_positions() -> np.ndarray

Position of every sample, in the order of samples.

Returns:

Type	Description
`ndarray`	`(n_samples, 3)`.

sample_times ¶

sample_times() -> np.ndarray

Time of every sample after its waveform's first sample (ns).

Returns:

Type	Description
`ndarray`	Length `n_samples`.

subset ¶

subset(mask) -> Waveforms

Select waveforms.

Parameters:

Name	Type	Description	Default
`mask`		Boolean array of length `n_waveforms`, or row indices.	required

Returns:

Type	Description
`Waveforms`	The selected rows with their samples, repacked.

Raises:

Type	Description
`ValueError`	If a boolean mask has the wrong length.

concatenate `classmethod` ¶

concatenate(parts: list[Waveforms]) -> Waveforms

Stack waveform sets, e.g. the chunks of chunks.

Parameters:

Name	Type	Description	Default
`parts`	`list[Waveforms]`	Waveform sets.	required

Returns:

Type	Description
`Waveforms`	All rows in order; only attributes present in every part are kept.

Raises:

Type	Description
`ValueError`	If `parts` is empty.

write_las ¶

write_las(path: str | Path, scale: float = 0.001, bits: int = 16, external: bool = False, crs_wkt: str | None = None) -> None

Write as LAS 1.4 with waveform data packets; see write_las.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	As for `write_las`.	required
`scale`	`str \| Path`	As for `write_las`.	required
`bits`	`str \| Path`	As for `write_las`.	required
`external`	`str \| Path`	As for `write_las`.	required
`crs_wkt`	`str \| Path`	As for `write_las`.	required

write_pulsewaves ¶

write_pulsewaves(path: str | Path, scale: float = 0.0001) -> None

Write as PulseWaves; see write_pulsewaves.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	As for `write_pulsewaves`.	required
`scale`	`str \| Path`	As for `write_pulsewaves`.	required

decompose ¶

decompose(**kwargs) -> Echoes

Gaussian decomposition; see decompose.

Parameters:

Name	Type	Description	Default
`**kwargs`		Parameters of `decompose`.	`{}`

Returns:

Type	Description
`Echoes`

Echoes `dataclass` ¶

Echoes found by decompose (or the truth of sylva.synthetic.waveforms), one row per echo, grouped by waveform in ascending order.

Attributes:

Name	Type	Description
`waveform`	`ndarray`	Row of the waveform each echo belongs to.
`time`	`ndarray`	Time of the echo's peak after the waveform's first sample (ns).
`amplitude`	`ndarray`	Peak of the fitted Gaussian above the background (sample units).
`width`	`ndarray`	Standard deviation of the fitted Gaussian (ns). The full width at half maximum is `2.3548 * width`.
`xyz`	`ndarray`	`(n, 3)` echo positions.
`range`	`ndarray`	Range from the scanner along the beam (m); NaN where the origin is unknown.
`stats`	`dict[str, ndarray]`	Per-waveform results: `background`, `noise` (standard deviation), `rmse` of the fit, `n_echoes` and `iterations`.

energy `property` ¶

energy: ndarray

Area under each echo, amplitude * width * sqrt(2 pi) (sample units times ns).

cross_section ¶

cross_section(calibration: float = 1.0, range=None) -> np.ndarray

Backscatter cross-section of each echo; see backscatter_cross_section.

Parameters:

Name	Type	Description	Default
`calibration`	`float`	Calibration constant `C_cal` (see `calibration_constant`); 1 gives relative values.	`1.0`
`range`		Ranges (m) to use instead of `range`, e.g. from a trajectory when the file has no scanner positions.	`None`

Returns:

Type	Description
`ndarray`	`C_cal R^4 amplitude width`; NaN where the range is unknown.

to_pointcloud ¶

to_pointcloud() -> PointCloud

The echoes as points.

Returns:

Type	Description
`PointCloud`	One point per echo with `amplitude`, `width`, `time`, `range` and `waveform` attributes.

to_shots ¶

to_shots(waveforms: Waveforms, origin=None) -> Shots

Pulses with these echoes; see to_shots.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms the echoes were found in.	required
`origin`		Scanner positions per waveform, as for `to_shots`.	`None`

Returns:

Type	Description
`Shots`

info ¶

info(path: str | Path) -> dict

Describe a waveform file without reading its waveforms.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	LAS/LAZ with waveform data packets, or a PulseWaves `.pls`.	required

Returns:

Type	Description
`dict`	`format` (`"las"` or `"pulsewaves"`), `version`, `bounds`, `compressed`, the wave packet or pulse `descriptors`; for LAS also `point_format`, `n_records` and `packets` (`"internal"`, `"external"` with the `wdp` path, or `"missing"`); for PulseWaves `n_pulses`, `system`, `software`, `waves` (the `.wvs` path) and `lookup_tables`.

Raises:

Type	Description
`OSError`	If the file is missing, is not a waveform file, or is a RIEGL SDF file (which only RIEGL's proprietary library reads).

read ¶

read(path: str | Path, start: int = 0, count: int | None = None, kind: str = 'returning', dedupe: bool = True, lookup: bool = True) -> Waveforms

Read waveforms from a LAS/LAZ or PulseWaves file.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	LAS 1.3 / 1.4 (or LAZ) with point format 4, 5, 9 or 10, with its waveform packets inside or in a `.wdp` file of the same name; or a PulseWaves `.pls` with its `.wvs`.	required
`start`	`int`	Point records (LAS) or pulses (PulseWaves) to read, from `start`; all by default. Use `chunks` to stream a large file.	`0`
`count`	`int`	Point records (LAS) or pulses (PulseWaves) to read, from `start`; all by default. Use `chunks` to stream a large file.	`0`
`kind`	`str`	PulseWaves only: `"returning"` (default), `"outgoing"` or `"all"` samplings.	`'returning'`
`dedupe`	`bool`	LAS only: consecutive point records that point to the same packet (the returns of one pulse) give one waveform; False gives one per record.	`True`
`lookup`	`bool`	PulseWaves only: turn sample values into physical units with the file's lookup tables.	`True`

Returns:

Type	Description
`Waveforms`

Raises:

Type	Description
`OSError`	If the file is missing or unreadable, its packets cannot be found, or it uses an unsupported variant (compressed packets, PulseWaves compression).
`ValueError`	For a bad `kind`, `start` or `count`.

chunks ¶

chunks(path: str | Path, size: int = 100000, kind: str = 'returning', dedupe: bool = True, lookup: bool = True) -> Iterator[Waveforms]

Stream a waveform file in chunks.

Each chunk holds about size point records (LAS) or pulses (PulseWaves) and the samples they refer to, so memory stays bounded whatever the file size. A pulse belongs to the chunk it starts in: in a LAS file the records that continue a chunk's last packet are read with it, so no waveform is split or repeated.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	As for `read`.	required
`size`	`int`	Records or pulses per chunk.	`100000`
`kind`	`str`	As for `read`.	`'returning'`
`dedupe`	`str`	As for `read`.	`'returning'`
`lookup`	`str`	As for `read`.	`'returning'`

Yields:

Type	Description
`Waveforms`

Raises:

Type	Description
`ValueError`	If `size` is not positive.

write_las ¶

write_las(waveforms: Waveforms, path: str | Path, scale: float = 0.001, bits: int = 16, external: bool = False, crs_wkt: str | None = None) -> None

Write waveforms as LAS 1.4 point format 9 with waveform data packets.

Each waveform becomes one point at its anchor, with the LAS waveform fields that put its samples where they were: the packet offset and size, the return point location -offset (ps) and the vector -direction * metres_per_ns (m/ps). gps_time and the attributes intensity, return_number, number_of_returns, classification and point_source_id are written where present. One wave packet descriptor is written per distinct (number of samples, interval) pair, at most 255.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms. The origin and `pulse` are not stored (LAS has no field for them).	required
`path`	`str \| Path`	Output `.las` file.	required
`scale`	`float`	Coordinate resolution (m).	`0.001`
`bits`	`int`	8 or 16 bits per sample. Integer samples that fit are stored as they are; others with a digitiser gain and offset spanning their range (the rounding error is at most half the gain).	`16`
`external`	`bool`	Put the packets in a `.wdp` file next to `path` instead of an extended VLR.	`False`
`crs_wkt`	`str \| None`	Coordinate reference system as OGC WKT, stored in a VLR.	`None`

Raises:

Type	Description
`ValueError`	For bad arguments, more than 255 descriptors, or an interval that is not a whole number of picoseconds.

write_pulsewaves ¶

write_pulsewaves(waveforms: Waveforms, path: str | Path, scale: float = 0.0001) -> None

Write waveforms as PulseWaves 0.3 (path and a .wvs next to it).

Consecutive rows of one pulse on the same beam become one pulse whose outgoing (kind == 1) and returning rows are the segments of an outgoing and a returning sampling. The scanner position becomes the anchor where it is known (the file then says the optical centre and the anchor coincide). Samples are stored as 16-bit integers, segment starts to a thousandth of a sample.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms; samples must round to integers from 0 to 65535.	required
`path`	`str \| Path`	Output `.pls` file.	required
`scale`	`float`	Resolution of the anchor and target coordinates (m).	`0.0001`

Raises:

Type	Description
`ValueError`	If a sample does not fit 16 bits or a segment has more than 65535 samples.

estimate_noise ¶

estimate_noise(waveforms: Waveforms) -> tuple[np.ndarray, np.ndarray]

Background level and noise of each waveform.

First robustly: the median of the samples and 1.4826 times their median absolute deviation (the standard deviation, for Gaussian noise), recomputed without the samples more than three standard deviations above the median until these settle, so that echoes do not inflate them. Then the mean and standard deviation of the samples within three standard deviations of that background (divided by 0.9866 for the clipped tails), iterated; this keeps the noise right for digitised samples, whose median absolute deviation moves in whole counts.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms.	required

Returns:

Type	Description
`background, noise : numpy.ndarray`	Per waveform, in sample units.

smooth ¶

smooth(waveforms: Waveforms, sigma: float) -> Waveforms

Smooth the samples with a Gaussian kernel.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms.	required
`sigma`	`float`	Standard deviation of the kernel (ns), truncated at four; near the ends of a waveform the kernel is renormalised over the samples it covers.	required

Returns:

Type	Description
`Waveforms`	A copy with smoothed samples.

Raises:

Type	Description
`ValueError`	If `sigma` is negative.

decompose ¶

decompose(waveforms: Waveforms, smooth: float = 1.0, threshold: float = 4.0, min_amplitude: float = 0.0, peaks: str = 'inflection', max_echoes: int = 10, min_width: float = 0.3, max_width: float = 20.0, noise: float | None = None, background: float | None = None, max_iter: int = 100, refine: bool = True) -> Echoes

Gaussian decomposition of waveforms into echoes (Wagner et al. 2006).

Each waveform is modelled as a constant background plus one Gaussian per echo. The background and noise come from estimate_noise (or the values given). The waveform is smoothed with a Gaussian of smooth ns, and candidate echoes are taken from the smoothed signal where it exceeds the threshold: with peaks="derivative" at the zero crossings of its first derivative (the local maxima, as Wagner et al. 2006), with peaks="inflection" at the local minima of its second derivative, the width from the inflection points around them (Hofton et al. 2000), which also separates echoes that overlap too much to make a maximum of their own. The amplitude, position and width of every candidate are then fitted together to the raw samples by Levenberg-Marquardt least squares. Fitted echoes weaker than the threshold are dropped, coincident ones merged, and with refine an echo is added where the residual still shows one, until the fit settles. Waveforms are processed in parallel; the result does not depend on the number of threads.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms (outgoing ones too, if they are in the set).	required
`smooth`	`float`	Standard deviation of the smoothing kernel (ns) used to find the candidates; 0 for none. The fit always uses the raw samples.	`1.0`
`threshold`	`float`	Detection threshold in noise standard deviations.	`4.0`
`min_amplitude`	`float`	Smallest echo amplitude (sample units), whatever the noise.	`0.0`
`peaks`	`str`	`"inflection"` or `"derivative"`.	`'inflection'`
`max_echoes`	`int`	Most echoes per waveform.	`10`
`min_width`	`float`	Bounds on each echo's standard deviation (ns).	`0.3`
`max_width`	`float`	Bounds on each echo's standard deviation (ns).	`0.3`
`noise`	`float \| None`	Values to use instead of the estimates (sample units).	`None`
`background`	`float \| None`	Values to use instead of the estimates (sample units).	`None`
`max_iter`	`int`	Levenberg-Marquardt iterations per fit.	`100`
`refine`	`bool`	Look for echoes in the residual of the fit.	`True`

Returns:

Type	Description
`Echoes`	With per-waveform `stats`.

Raises:

Type	Description
`ValueError`	For bad parameters.

to_shots ¶

to_shots(waveforms: Waveforms, echoes: Echoes, origin=None) -> Shots

Pulses from decomposed waveforms, for ray tracing and canopy metrics.

One shot per pulse (the consecutive rows sharing pulse; outgoing rows are ignored), with the echoes of all its returning rows sorted by range. A pulse in which no echo was found stays as a shot without an echo, which ray-traced products count as free space.

Parameters:

Name	Type	Description	Default
`waveforms`	`Waveforms`	The waveforms.	required
`echoes`	`Echoes`	Echoes found in them by `decompose`.	required
`origin`		`(n_waveforms, 3)` scanner positions (e.g. from a trajectory), used where finite. Otherwise the waveform's own origin is used, and where that is unknown too the shot starts at the pulse's first sample: its free space is then only what the record covers.	`None`

Returns:

Type	Description
`Shots`	Echo attributes `amplitude`, `width` (ns), `time` (ns after the first sample of its waveform) and `waveform` (row).

Raises:

Type	Description
`ValueError`	If `origin` does not have one row per waveform, or the echoes refer to rows that are not there.

backscatter_cross_section ¶

backscatter_cross_section(range, amplitude, width, calibration: float = 1.0) -> np.ndarray

Backscatter cross-section of echoes (Wagner et al. 2006).

By the radar equation for a Gaussian system pulse, the cross-section of an echo is sigma = C_cal R^4 P s: the fourth power of the range times the echo's amplitude P and width s, scaled by a calibration constant that folds in the transmitted power, the receiver and the atmosphere.

Parameters:

Name	Type	Description	Default
`range`		Ranges (m).	required
`amplitude`		Echo amplitudes.	required
`width`		Echo widths (ns, standard deviation).	required
`calibration`	`float`	`C_cal`; 1 gives relative values, comparable within a survey.	`1.0`

Returns:

Type	Description
`ndarray`	Cross-sections (m² once calibrated).

Raises:

Type	Description
`ValueError`	If the arrays differ in length.

calibration_constant ¶

calibration_constant(range, amplitude, width, reflectance, beam_divergence: float, incidence=0.0) -> float

Calibration constant C_cal from echoes of reference targets.

For an extended Lambertian target of reflectance rho hit at incidence angle alpha by a beam of divergence beta (full angle, rad), sigma = pi rho R^2 beta^2 cos(alpha) (Wagner 2010), so each echo gives C_cal = pi rho beta^2 cos(alpha) / (R^2 P s); the median over the echoes is returned.

Parameters:

Name	Type	Description	Default
`range`		Range (m), amplitude and width (ns) of the reference echoes.	required
`amplitude`		Range (m), amplitude and width (ns) of the reference echoes.	required
`width`		Range (m), amplitude and width (ns) of the reference echoes.	required
`reflectance`		Reflectance of each target (0 to 1), or one value for all.	required
`beam_divergence`	`float`	Full beam divergence (rad).	required
`incidence`		Incidence angle (rad), per echo or one value.	`0.0`

Returns:

Type	Description
`float`

Raises:

Type	Description
`ValueError`	If no echo gives a finite, positive value.

sylva.waveform¶

sylva.waveform ¶

Waveform dataclass ¶

times ¶

positions ¶

Waveforms dataclass ¶

n_waveforms property ¶

n_samples property ¶

sample_positions ¶

sample_times ¶

subset ¶

concatenate classmethod ¶

write_las ¶

write_pulsewaves ¶

decompose ¶

Echoes dataclass ¶

energy property ¶

cross_section ¶

to_pointcloud ¶

to_shots ¶

info ¶

read ¶

chunks ¶

write_las ¶

write_pulsewaves ¶

estimate_noise ¶

smooth ¶

decompose ¶

to_shots ¶

backscatter_cross_section ¶

calibration_constant ¶

Waveform `dataclass` ¶

Waveforms `dataclass` ¶

n_waveforms `property` ¶

n_samples `property` ¶

concatenate `classmethod` ¶

Echoes `dataclass` ¶

energy `property` ¶