10. A RIEGL project, end to end¶

This notebook runs the whole Sylva chain on one real plot, straight from the scanner's files: a RiSCAN PRO project of the TERN Litchfield Savanna SuperSite core hectare, scanned in July 2021 with a RIEGL VZ-2000i from 64 positions (36 on a 20 m grid inside the plot, 28 on a ring around it).

Step	Sylva	Output
1. Project	read_riscan_project	scan positions, SOPs, scan pattern
2. Read	io.read_rxp_shots, Shots.fill_missing, filters.voxel_downsample	2 cm plot cloud, pulse file with misses
3. Registration	coreg.prepare_scan, coreg.coregister_prepared	the ring registered into the inner grid, checked on the ground
4. Ground	ground.classify_ground_csf, make_dtm, normalize_height, make_chm	DTM, CHM, heights
5. Trees	trees.detect_stems, merge_branches, segment_trees, prune_trees, crown_metrics_all	tree table, segmented cloud
6. Scan quality	quality.stem_noise	range noise, per-scan horizontal registration
7. Wood and leaves	qsm.wood_points, build_qsm, QSM.metrics, leaves.*	cylinder models, leaf meshes
8. Canopy	canopy.GapProfile	PAI, PAVD profile, clumping
9. Voxels	voxels.ray_voxelize, occlusion_profile, tree_sampling	plant area density, what was seen

It needs RiVLib (see Point clouds and files) and about 40 GB of memory at its peak. The first run reads 35 GB from the project, which takes about 25 minutes; the thinned read is cached in OUT, so later runs start at step 3 within a minute. The data are TERN's.

Step 3 finds that the ring positions are not registered with the inner grid (metres off in height) and registers them into it with sylva.coreg, so everything after it uses all 64 scans. Checking registration before anything else is the lesson of this plot.

In [1]:

import os
from pathlib import Path
import time

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import sylva
import json
import pickle

from sylva import canopy, coreg, filters, ground, io, leaves, qsm, quality, trees, voxels

# The TERN Litchfield core plot, in the folder SYLVA_TERN_DIR names.
PROJECT = Path(os.environ.get("SYLVA_TERN_DIR", "TERN_TLS_RAW")) / "LSS_2021_07_core.RiSCAN"
OUT = Path(os.environ.get("SYLVA_RUNS_DIR", Path.home() / "sylva_runs")) / PROJECT.stem   # cached read, products
OUT.mkdir(parents=True, exist_ok=True)

PLOT = (0.0, -100.0, 100.0, 0.0)    # xmin, ymin, xmax, ymax of the core hectare (project frame, m)
BUFFER = 10.0                       # ground is classified this far beyond the plot: a cloth needs both sides
STRIDE = 8                          # keep every 8th pulse when reading
RAY_EVERY = 4                       # ... and every 4th of those for the gap profile and voxels
VOXEL = 0.02                        # point spacing of the plot cloud (m)

plt.style.use("sylva.mplstyle")     # shared figure style; C0 to C7 are the Okabe-Ito colours
GROUND, WOOD, LEAF, GRASS, CONTEXT = "#997A5C", "#4D2B12", "#009E73", "#E69F00", "0.8"   # the same in every notebook
DIVERGING = "RdBu_r"                # for signed differences, centred on zero
LABELS = ListedColormap(["#332288", "#88CCEE", "#44AA99", "#117733", "#999933", "#DDCC77",
                         "#CC6677", "#882255", "#AA4499"])   # tree ids and clusters: Paul Tol's muted set
print("sylva", sylva.__version__, "| RiVLib:", io.find_rivlib().name)

import os from pathlib import Path import time import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap import sylva import json import pickle from sylva import canopy, coreg, filters, ground, io, leaves, qsm, quality, trees, voxels # The TERN Litchfield core plot, in the folder SYLVA_TERN_DIR names. PROJECT = Path(os.environ.get("SYLVA_TERN_DIR", "TERN_TLS_RAW")) / "LSS_2021_07_core.RiSCAN" OUT = Path(os.environ.get("SYLVA_RUNS_DIR", Path.home() / "sylva_runs")) / PROJECT.stem # cached read, products OUT.mkdir(parents=True, exist_ok=True) PLOT = (0.0, -100.0, 100.0, 0.0) # xmin, ymin, xmax, ymax of the core hectare (project frame, m) BUFFER = 10.0 # ground is classified this far beyond the plot: a cloth needs both sides STRIDE = 8 # keep every 8th pulse when reading RAY_EVERY = 4 # ... and every 4th of those for the gap profile and voxels VOXEL = 0.02 # point spacing of the plot cloud (m) plt.style.use("sylva.mplstyle") # shared figure style; C0 to C7 are the Okabe-Ito colours GROUND, WOOD, LEAF, GRASS, CONTEXT = "#997A5C", "#4D2B12", "#009E73", "#E69F00", "0.8" # the same in every notebook DIVERGING = "RdBu_r" # for signed differences, centred on zero LABELS = ListedColormap(["#332288", "#88CCEE", "#44AA99", "#117733", "#999933", "#DDCC77", "#CC6677", "#882255", "#AA4499"]) # tree ids and clusters: Paul Tol's muted set print("sylva", sylva.__version__, "| RiVLib:", io.find_rivlib().name)

sylva 0.2.1 | RiVLib: libscanifc-mt.so

1. The project¶

read_riscan_project parses project.rsp: every scan position with its .rxp, its SOP (scanner to project transform) and the angular scan pattern. Nothing is read from the scans yet.

In [2]:

project = sylva.read_riscan_project(PROJECT)
positions = project.with_scans()
origins = np.array([p.sop[:3, 3] for p in positions])
tilt = np.array([np.degrees(np.arccos(p.sop[2, 2])) for p in positions])
pat = positions[0].pattern
print(f"{project.name}: {len(project)} positions, {len(positions)} with a scan and a SOP; {positions[0].instrument}")
print(f"zenith {pat['theta_start']:.0f}-{pat['theta_start'] + pat['theta_delta'] * pat['theta_count']:.0f} deg "
      f"in {pat['theta_count']} lines, {pat['phi_count']} azimuth steps: "
      f"{pat['theta_count'] * pat['phi_count'] / 1e6:.0f} M pulses per scan; tilt up to {tilt.max():.1f} deg")

in_plot = ((origins[:, 0] >= PLOT[0] - 1) & (origins[:, 0] <= PLOT[2] + 1)
           & (origins[:, 1] >= PLOT[1] - 1) & (origins[:, 1] <= PLOT[3] + 1))
fig, ax = plt.subplots(figsize=(6.5, 5))
ax.add_patch(plt.Rectangle(PLOT[:2], PLOT[2] - PLOT[0], PLOT[3] - PLOT[1], fill=False, lw=1.2))
ax.scatter(*origins[in_plot, :2].T, c="C0", s=20, label=f"inside the plot ({in_plot.sum()})")
ax.scatter(*origins[~in_plot, :2].T, c="C1", s=20, label=f"ring ({(~in_plot).sum()})")
ax.set(aspect="equal", xlabel="x, project (m)", ylabel="y, project (m)", title="Scan positions and the core hectare")
ax.legend(loc="upper left", bbox_to_anchor=(1.0, 1.0));

project = sylva.read_riscan_project(PROJECT) positions = project.with_scans() origins = np.array([p.sop[:3, 3] for p in positions]) tilt = np.array([np.degrees(np.arccos(p.sop[2, 2])) for p in positions]) pat = positions[0].pattern print(f"{project.name}: {len(project)} positions, {len(positions)} with a scan and a SOP; {positions[0].instrument}") print(f"zenith {pat['theta_start']:.0f}-{pat['theta_start'] + pat['theta_delta'] * pat['theta_count']:.0f} deg " f"in {pat['theta_count']} lines, {pat['phi_count']} azimuth steps: " f"{pat['theta_count'] * pat['phi_count'] / 1e6:.0f} M pulses per scan; tilt up to {tilt.max():.1f} deg") in_plot = ((origins[:, 0] >= PLOT[0] - 1) & (origins[:, 0] <= PLOT[2] + 1) & (origins[:, 1] >= PLOT[1] - 1) & (origins[:, 1] <= PLOT[3] + 1)) fig, ax = plt.subplots(figsize=(6.5, 5)) ax.add_patch(plt.Rectangle(PLOT[:2], PLOT[2] - PLOT[0], PLOT[3] - PLOT[1], fill=False, lw=1.2)) ax.scatter(*origins[in_plot, :2].T, c="C0", s=20, label=f"inside the plot ({in_plot.sum()})") ax.scatter(*origins[~in_plot, :2].T, c="C1", s=20, label=f"ring ({(~in_plot).sum()})") ax.set(aspect="equal", xlabel="x, project (m)", ylabel="y, project (m)", title="Scan positions and the core hectare") ax.legend(loc="upper left", bbox_to_anchor=(1.0, 1.0));

LSS_2021_07_core: 66 positions, 64 with a scan and a SOP; VZ-2000i
zenith 30-130 deg in 3311 lines, 12001 azimuth steps: 40 M pulses per scan; tilt up to 3.8 deg

2. One pass over the scans¶

Reading a scan takes about 20 s per 0.6 GB scan, whatever the thinning, so every scan is read once, as pulses, and everything is derived from that read:

• shot_stride=STRIDE keeps every 8th pulse with all its echoes.
• Points. The pulses' echoes, moved into the project frame with the SOP, cropped to the plot plus BUFFER and thinned to 2 cm.
• Pulses for the canopy. Every RAY_EVERY-th pulse. RiVLib only streams pulses that returned something, so the pulses that went to the sky are added back with fill_missing. Without them the free space above the canopy looks unsampled and plant area is overestimated, more and more with height. The number fired per zenith line is counted on the downward lines (100-125 deg), where every pulse hits the ground. This has to happen in the scanner frame, before the SOP.

Only the attributes used later are kept, which keeps the pass under about 25 GB.

In [3]:

cloud_path, rays_path, counts_path = OUT / "cloud_2cm.laz", OUT / "rays.parquet", OUT / "ray_counts.npy"
features_path = OUT / "coreg_features.pkl"
lo = np.array([PLOT[0] - BUFFER, PLOT[1] - BUFFER, -np.inf])
hi = np.array([PLOT[2] + BUFFER, PLOT[3] + BUFFER, np.inf])


def fired_per_line(s, pattern):
    """Median pulses per zenith line on the downward lines, where every pulse returns."""
    theta = pattern["theta_start"] + pattern["theta_delta"] * np.arange(pattern["theta_count"])   # line zeniths
    edges = np.r_[theta - pattern["theta_delta"] / 2, theta[-1] + pattern["theta_delta"] / 2]
    observed, _ = np.histogram(s.zenith_azimuth()[0], bins=edges)
    return int(round(np.median(observed[(theta >= 100) & (theta <= 125)])))


if not all(p.exists() for p in (cloud_path, rays_path, counts_path, features_path)):
    t0 = time.time()
    parts, rays, features = [], [], []
    for i, pos in enumerate(positions):
        s = io.read_rxp_shots(pos.rxp, shot_stride=STRIDE)            # scanner frame
        r = s.subset(np.arange(s.n_shots) % RAY_EVERY == 0)
        r = r.fill_missing(pos.pattern, pulses_per_line=fired_per_line(r, pos.pattern))
        r = r.transform(pos.sop)
        r.echo_attrs = {}
        rays.append(r)
        pts = s.transform(pos.sop).to_pointcloud()
        # Registration features from everything within 40 m of the scanner.
        near = pts[np.linalg.norm(pts.xyz - pos.sop[:3, 3], axis=1) < 40.0]
        features.append(coreg.prepare_scan(sylva.PointCloud(near.xyz), name=pos.name, origin=pos.sop[:3, 3]))
        pts = pts[np.all((pts.xyz >= lo) & (pts.xyz <= hi), axis=1)]
        pts = sylva.PointCloud(pts.xyz, {"reflectance": pts.attrs["reflectance"]})
        pts = filters.voxel_downsample(pts, VOXEL)
        parts.append(pts.with_attrs(scan_id=np.full(len(pts), i, np.int32)))
        if i % 16 == 0:
            print(f"{i:2d} {pos.name}: {len(pts):,} points, {r.n_shots:,} pulses ({time.time() - t0:.0f} s)")
    cloud = filters.voxel_downsample(sylva.PointCloud.concatenate(parts), VOXEL)
    del parts
    sylva.write(cloud, cloud_path)
    np.save(counts_path, np.array([r.n_shots for r in rays]))
    pickle.dump(features, open(features_path, "wb"))
    shots = sylva.Shots.concatenate(rays)
    del rays
    shots.save(rays_path)
    print(f"read {len(positions)} scans in {time.time() - t0:.0f} s")
else:
    cloud = sylva.read(cloud_path)
    shots = sylva.Shots.load(rays_path)
    features = pickle.load(open(features_path, "rb"))
ray_counts = np.load(counts_path)
print(f"plot cloud: {len(cloud):,} points at {VOXEL * 100:.0f} cm; pulses: {shots.n_shots:,} "
      f"({(shots.echo_count == 0).mean():.0%} misses), {shots.n_echoes:,} echoes")

cloud_path, rays_path, counts_path = OUT / "cloud_2cm.laz", OUT / "rays.parquet", OUT / "ray_counts.npy" features_path = OUT / "coreg_features.pkl" lo = np.array([PLOT[0] - BUFFER, PLOT[1] - BUFFER, -np.inf]) hi = np.array([PLOT[2] + BUFFER, PLOT[3] + BUFFER, np.inf]) def fired_per_line(s, pattern): """Median pulses per zenith line on the downward lines, where every pulse returns.""" theta = pattern["theta_start"] + pattern["theta_delta"] * np.arange(pattern["theta_count"]) # line zeniths edges = np.r_[theta - pattern["theta_delta"] / 2, theta[-1] + pattern["theta_delta"] / 2] observed, _ = np.histogram(s.zenith_azimuth()[0], bins=edges) return int(round(np.median(observed[(theta >= 100) & (theta <= 125)]))) if not all(p.exists() for p in (cloud_path, rays_path, counts_path, features_path)): t0 = time.time() parts, rays, features = [], [], [] for i, pos in enumerate(positions): s = io.read_rxp_shots(pos.rxp, shot_stride=STRIDE) # scanner frame r = s.subset(np.arange(s.n_shots) % RAY_EVERY == 0) r = r.fill_missing(pos.pattern, pulses_per_line=fired_per_line(r, pos.pattern)) r = r.transform(pos.sop) r.echo_attrs = {} rays.append(r) pts = s.transform(pos.sop).to_pointcloud() # Registration features from everything within 40 m of the scanner. near = pts[np.linalg.norm(pts.xyz - pos.sop[:3, 3], axis=1) < 40.0] features.append(coreg.prepare_scan(sylva.PointCloud(near.xyz), name=pos.name, origin=pos.sop[:3, 3])) pts = pts[np.all((pts.xyz >= lo) & (pts.xyz <= hi), axis=1)] pts = sylva.PointCloud(pts.xyz, {"reflectance": pts.attrs["reflectance"]}) pts = filters.voxel_downsample(pts, VOXEL) parts.append(pts.with_attrs(scan_id=np.full(len(pts), i, np.int32))) if i % 16 == 0: print(f"{i:2d} {pos.name}: {len(pts):,} points, {r.n_shots:,} pulses ({time.time() - t0:.0f} s)") cloud = filters.voxel_downsample(sylva.PointCloud.concatenate(parts), VOXEL) del parts sylva.write(cloud, cloud_path) np.save(counts_path, np.array([r.n_shots for r in rays])) pickle.dump(features, open(features_path, "wb")) shots = sylva.Shots.concatenate(rays) del rays shots.save(rays_path) print(f"read {len(positions)} scans in {time.time() - t0:.0f} s") else: cloud = sylva.read(cloud_path) shots = sylva.Shots.load(rays_path) features = pickle.load(open(features_path, "rb")) ray_counts = np.load(counts_path) print(f"plot cloud: {len(cloud):,} points at {VOXEL * 100:.0f} cm; pulses: {shots.n_shots:,} " f"({(shots.echo_count == 0).mean():.0%} misses), {shots.n_echoes:,} echoes")

plot cloud: 112,165,249 points at 2 cm; pulses: 80,575,741 (25% misses), 77,848,847 echoes

3. Registering the outer ring¶

The scans are only as good as their SOPs. The 36 inner-grid scans were registered together in RiSCAN; the 28 on the ring around the plot carry little more than the scanner's own GNSS fix (the project gives them 0.8 m horizontal and 1.3 m vertical accuracy). sylva.coreg registers them into the inner grid from the trees themselves:

• Stems. Each scan's stems are matched against the others' with no initial guess. Stem pairs at equal separation are candidate matches, which gives yaw and translation.
• ICP. Point-to-plane ICP on locally planar points (stems, ground, logs) refines each pair. A pair is accepted only if it also fits above the ground.
• Pose graph. The accepted pairs are solved together, with the inner scans held fixed.
• Prior. The SOPs serve as a prior on where each scanner stood. Their headings are not trusted, and two turn out to be wrong by 123 and 155 degrees.

The check is the ground: build a DTM from the inner scans, take each scan's lowest point per 0.5 m cell within 25 m of it, and read the median height above that DTM. A scan registered with the others sits a few centimetres above it (grass and litter).

In [4]:

t0 = time.time()
fixed = {k: np.eye(4) for k in np.flatnonzero(in_plot)}
reg = coreg.coregister_prepared(features, coreg.CoregConfig(max_pair_distance=40.0, workers=8, verbose=False),
                                approximate_positions=origins, fixed=fixed, priors=[np.eye(4)] * len(positions))
corr = reg.poses                                  # correction applied on top of each SOP
registered = np.asarray(reg.registered)
print(f"{registered.sum()} of {len(positions)} scans registered in {time.time() - t0:.0f} s; "
      f"{sum(p.success for p in reg.pairs)} accepted pairs")
ring = pd.DataFrame({
    "scan": [positions[k].name for k in np.flatnonzero(~in_plot)],
    "shift_m": [np.linalg.norm(features[k].location(corr[k]) - origins[k]) for k in np.flatnonzero(~in_plot)],
    "dz_m": [(features[k].location(corr[k]) - origins[k])[2] for k in np.flatnonzero(~in_plot)],
    "rotation_deg": [np.degrees(np.linalg.norm(coreg.se3_log(corr[k])[:3])) for k in np.flatnonzero(~in_plot)],
}).round(2)
json.dump({p.name: (T @ p.sop).tolist() for p, T in zip(positions, corr)}, open(OUT / "sop_corrected.json", "w"))
ring.sort_values("shift_m", ascending=False).head(10)

t0 = time.time() fixed = {k: np.eye(4) for k in np.flatnonzero(in_plot)} reg = coreg.coregister_prepared(features, coreg.CoregConfig(max_pair_distance=40.0, workers=8, verbose=False), approximate_positions=origins, fixed=fixed, priors=[np.eye(4)] * len(positions)) corr = reg.poses # correction applied on top of each SOP registered = np.asarray(reg.registered) print(f"{registered.sum()} of {len(positions)} scans registered in {time.time() - t0:.0f} s; " f"{sum(p.success for p in reg.pairs)} accepted pairs") ring = pd.DataFrame({ "scan": [positions[k].name for k in np.flatnonzero(~in_plot)], "shift_m": [np.linalg.norm(features[k].location(corr[k]) - origins[k]) for k in np.flatnonzero(~in_plot)], "dz_m": [(features[k].location(corr[k]) - origins[k])[2] for k in np.flatnonzero(~in_plot)], "rotation_deg": [np.degrees(np.linalg.norm(coreg.se3_log(corr[k])[:3])) for k in np.flatnonzero(~in_plot)], }).round(2) json.dump({p.name: (T @ p.sop).tolist() for p, T in zip(positions, corr)}, open(OUT / "sop_corrected.json", "w")) ring.sort_values("shift_m", ascending=False).head(10)

64 of 64 scans registered in 860 s; 116 accepted pairs

Out[4]:

	scan	shift_m	dz_m	rotation_deg
26	ScanPos065	3.24	1.64	1.79
18	ScanPos057	3.21	2.80	3.12
3	ScanPos040	3.04	2.95	123.44
15	ScanPos054	2.82	-2.78	0.74
27	ScanPos066	2.58	1.03	1.31
24	ScanPos063	2.56	2.10	1.23
0	ScanPos037	2.15	0.95	0.49
25	ScanPos064	2.12	1.32	2.10
16	ScanPos055	2.00	1.62	0.64
17	ScanPos056	1.70	1.24	1.38

In [5]:

c5 = filters.voxel_downsample(cloud, 0.05)
sid = c5.attrs["scan_id"]
g = ground.classify_ground_csf(c5[in_plot[sid]], cloth_resolution=0.5, rigidness=2)
ref = ground.make_dtm(g, resolution=0.5, bounds=(lo[0], lo[1], hi[0], hi[1]))
del g


def ground_offset(xyz, origin):
    """Median height above the inner-grid DTM of a scan's lowest point per 0.5 m cell within 25 m."""
    r = np.hypot(xyz[:, 0] - origin[0], xyz[:, 1] - origin[1])
    xyz = xyz[(r > 2) & (r < 25)]
    if len(xyz) < 1000:
        return np.nan
    key = np.floor(xyz[:, :2] / 0.5).astype(np.int64)
    k = key[:, 0] * 1_000_000 + key[:, 1]
    order = np.lexsort((xyz[:, 2], k))
    low = xyz[order][np.r_[True, k[order][1:] != k[order][:-1]]]
    return float(np.median(low[:, 2] - ref.sample(low[:, 0], low[:, 1])))


before = np.array([ground_offset(c5.xyz[sid == k], origins[k]) for k in range(len(positions))])
after = np.array([ground_offset(coreg.transform_points(corr[k], c5.xyz[sid == k]), features[k].location(corr[k]))
                  for k in range(len(positions))])
del c5
print(f"ring scans: {np.nanmin(before[~in_plot]):+.2f} to {np.nanmax(before[~in_plot]):+.2f} m before, "
      f"{np.nanmin(after[~in_plot]):+.2f} to {np.nanmax(after[~in_plot]):+.2f} m after; "
      f"inner scans {np.nanmin(after[in_plot]):+.2f} to {np.nanmax(after[in_plot]):+.2f} m")
fig, axes = plt.subplots(1, 2, figsize=(10, 4.6), sharey=True)
for ax, v, title in [(axes[0], before, "As delivered"), (axes[1], after, "After sylva.coreg")]:
    sc_ = ax.scatter(*origins[:, :2].T, c=v, cmap=DIVERGING, vmin=-3, vmax=3, s=40, edgecolor="k", lw=0.3)
    ax.add_patch(plt.Rectangle(PLOT[:2], 100, 100, fill=False, lw=1))
    ax.set(aspect="equal", title=title, xlabel="x, project (m)")
axes[0].set_ylabel("y, project (m)")
fig.suptitle("Height of each scan's ground above the inner-grid DTM (white: in agreement)")
fig.colorbar(sc_, ax=axes, label="median ground height above the DTM (m)", shrink=0.9, extend="both");

c5 = filters.voxel_downsample(cloud, 0.05) sid = c5.attrs["scan_id"] g = ground.classify_ground_csf(c5[in_plot[sid]], cloth_resolution=0.5, rigidness=2) ref = ground.make_dtm(g, resolution=0.5, bounds=(lo[0], lo[1], hi[0], hi[1])) del g def ground_offset(xyz, origin): """Median height above the inner-grid DTM of a scan's lowest point per 0.5 m cell within 25 m.""" r = np.hypot(xyz[:, 0] - origin[0], xyz[:, 1] - origin[1]) xyz = xyz[(r > 2) & (r < 25)] if len(xyz) < 1000: return np.nan key = np.floor(xyz[:, :2] / 0.5).astype(np.int64) k = key[:, 0] * 1_000_000 + key[:, 1] order = np.lexsort((xyz[:, 2], k)) low = xyz[order][np.r_[True, k[order][1:] != k[order][:-1]]] return float(np.median(low[:, 2] - ref.sample(low[:, 0], low[:, 1]))) before = np.array([ground_offset(c5.xyz[sid == k], origins[k]) for k in range(len(positions))]) after = np.array([ground_offset(coreg.transform_points(corr[k], c5.xyz[sid == k]), features[k].location(corr[k])) for k in range(len(positions))]) del c5 print(f"ring scans: {np.nanmin(before[~in_plot]):+.2f} to {np.nanmax(before[~in_plot]):+.2f} m before, " f"{np.nanmin(after[~in_plot]):+.2f} to {np.nanmax(after[~in_plot]):+.2f} m after; " f"inner scans {np.nanmin(after[in_plot]):+.2f} to {np.nanmax(after[in_plot]):+.2f} m") fig, axes = plt.subplots(1, 2, figsize=(10, 4.6), sharey=True) for ax, v, title in [(axes[0], before, "As delivered"), (axes[1], after, "After sylva.coreg")]: sc_ = ax.scatter(*origins[:, :2].T, c=v, cmap=DIVERGING, vmin=-3, vmax=3, s=40, edgecolor="k", lw=0.3) ax.add_patch(plt.Rectangle(PLOT[:2], 100, 100, fill=False, lw=1)) ax.set(aspect="equal", title=title, xlabel="x, project (m)") axes[0].set_ylabel("y, project (m)") fig.suptitle("Height of each scan's ground above the inner-grid DTM (white: in agreement)") fig.colorbar(sc_, ax=axes, label="median ground height above the DTM (m)", shrink=0.9, extend="both");

ring scans: -2.72 to +2.86 m before, -0.08 to +0.24 m after; inner scans +0.01 to +0.09 m

The corrections are applied to the cached points and pulses, so nothing has to be read again. With the ring left as delivered, the DTM stepped down by metres in wedges behind the misregistered positions, and ghost stems appeared at the plot edge. From here on all 64 scans are used.

In [6]:

use = registered
for k in range(len(positions)):
    if use[k] and not np.allclose(corr[k], np.eye(4)):
        m = cloud.attrs["scan_id"] == k
        cloud.xyz[m] = coreg.transform_points(corr[k], cloud.xyz[m])
cloud = cloud[use[cloud.attrs["scan_id"]]]
start = np.r_[0, np.cumsum(ray_counts)]
origin, direction = shots.origin.copy(), shots.direction.copy()
keep = np.zeros(shots.n_shots, bool)
for k in np.flatnonzero(use):
    s_, e_ = start[k], start[k + 1]
    origin[s_:e_] = coreg.transform_points(corr[k], origin[s_:e_])
    direction[s_:e_] = direction[s_:e_] @ corr[k][:3, :3].T
    keep[s_:e_] = True
shots = sylva.Shots(origin, direction, shots.echo_start, shots.echo_count, shots.echo_range).subset(keep)
used_counts = ray_counts[use]
inner_used = in_plot[use]                       # the gap profile uses the scans inside the plot
origins = np.array([features[k].location(corr[k]) for k in range(len(positions))])
rays_path = OUT / "rays_registered.parquet"
shots.save(rays_path)
print(f"{use.sum()} scans: {len(cloud):,} points, {shots.n_shots:,} pulses")

use = registered for k in range(len(positions)): if use[k] and not np.allclose(corr[k], np.eye(4)): m = cloud.attrs["scan_id"] == k cloud.xyz[m] = coreg.transform_points(corr[k], cloud.xyz[m]) cloud = cloud[use[cloud.attrs["scan_id"]]] start = np.r_[0, np.cumsum(ray_counts)] origin, direction = shots.origin.copy(), shots.direction.copy() keep = np.zeros(shots.n_shots, bool) for k in np.flatnonzero(use): s_, e_ = start[k], start[k + 1] origin[s_:e_] = coreg.transform_points(corr[k], origin[s_:e_]) direction[s_:e_] = direction[s_:e_] @ corr[k][:3, :3].T keep[s_:e_] = True shots = sylva.Shots(origin, direction, shots.echo_start, shots.echo_count, shots.echo_range).subset(keep) used_counts = ray_counts[use] inner_used = in_plot[use] # the gap profile uses the scans inside the plot origins = np.array([features[k].location(corr[k]) for k in range(len(positions))]) rays_path = OUT / "rays_registered.parquet" shots.save(rays_path) print(f"{use.sum()} scans: {len(cloud):,} points, {shots.n_shots:,} pulses")

64 scans: 112,165,249 points, 80,575,741 pulses

4. Ground and heights¶

The cloth simulation filter runs on a 5 cm copy of the plot and its buffer, and the DTM covers the same area, so the plot edge is not a cloth edge. Heights are then added to the full 2 cm cloud and the plot is cut out.

A CHM takes the highest point in each cell, and 64 scans collect a few returns far above the canopy (birds, insects, mixed-pixel ghosts), up to 90 m here. Isolated points are removed before the CHM is built.

The CHM counts a 0.5 m cell as canopy if any point lies above the threshold, so its "cover" is crown extent, gaps within crowns included. In this open eucalypt forest that is more than half the plot, while only about a third of the pulses going up at 30 deg are stopped (step 8). The two answer different questions.

In [7]:

t0 = time.time()
g = ground.classify_ground_csf(filters.voxel_downsample(cloud, 0.05), cloth_resolution=0.5, rigidness=2)
dtm = ground.make_dtm(g, resolution=0.5, bounds=(lo[0], lo[1], hi[0], hi[1]))
del g
cloud = ground.normalize_height(cloud, dtm)
inside = ((cloud.x >= PLOT[0]) & (cloud.x <= PLOT[2]) & (cloud.y >= PLOT[1]) & (cloud.y <= PLOT[3]))
plot = cloud[inside]
del cloud
# A CHM takes the highest point per cell, so one stray return (a bird, a ghost) sets a cell:
# drop isolated points first.
tops = filters.voxel_downsample(plot, 0.05)
tops = filters.radius_outlier_removal(tops, radius=0.25, min_neighbors=5)
chm = ground.make_chm(tops, resolution=0.5, bounds=PLOT)
print(f"highest point {plot.attrs['height'].max():.0f} m, after removing isolated points {tops.attrs['height'].max():.0f} m")
del tops
dtm.to_ascii_grid(OUT / "dtm.asc")
chm.to_ascii_grid(OUT / "chm.asc")
print(f"{len(plot):,} points in the plot; terrain {np.nanmin(dtm.data):.1f} to {np.nanmax(dtm.data):.1f} m; "
      f"CHM cells with anything above 5 m: {canopy.canopy_cover(chm.data, 5.0):.0%}; {time.time() - t0:.0f} s")

X, Y = dtm.cell_centers()
dtm_plot = np.where((X >= PLOT[0]) & (X <= PLOT[2]) & (Y >= PLOT[1]) & (Y <= PLOT[3]), dtm.data, np.nan)
fig, axes = plt.subplots(1, 2, figsize=(10, 4.4), sharey=True)
for ax, data, r, title, label in [(axes[0], dtm_plot, dtm, "Terrain (DTM)", "elevation (m)"),
                                  (axes[1], chm.data, chm, "Canopy height (CHM)", "height above ground (m)")]:
    lo_, hi_ = np.nanpercentile(data, [1, 99.5])
    im = ax.imshow(data, origin="lower", extent=(r.xmin, r.xmax, r.ymin, r.ymax), vmin=lo_, vmax=hi_)
    ax.set(title=title, aspect="equal", xlim=PLOT[::2], ylim=PLOT[1::2], xlabel="x, project (m)")
    fig.colorbar(im, ax=ax, shrink=0.9, label=label)
axes[0].set_ylabel("y, project (m)");

t0 = time.time() g = ground.classify_ground_csf(filters.voxel_downsample(cloud, 0.05), cloth_resolution=0.5, rigidness=2) dtm = ground.make_dtm(g, resolution=0.5, bounds=(lo[0], lo[1], hi[0], hi[1])) del g cloud = ground.normalize_height(cloud, dtm) inside = ((cloud.x >= PLOT[0]) & (cloud.x <= PLOT[2]) & (cloud.y >= PLOT[1]) & (cloud.y <= PLOT[3])) plot = cloud[inside] del cloud # A CHM takes the highest point per cell, so one stray return (a bird, a ghost) sets a cell: # drop isolated points first. tops = filters.voxel_downsample(plot, 0.05) tops = filters.radius_outlier_removal(tops, radius=0.25, min_neighbors=5) chm = ground.make_chm(tops, resolution=0.5, bounds=PLOT) print(f"highest point {plot.attrs['height'].max():.0f} m, after removing isolated points {tops.attrs['height'].max():.0f} m") del tops dtm.to_ascii_grid(OUT / "dtm.asc") chm.to_ascii_grid(OUT / "chm.asc") print(f"{len(plot):,} points in the plot; terrain {np.nanmin(dtm.data):.1f} to {np.nanmax(dtm.data):.1f} m; " f"CHM cells with anything above 5 m: {canopy.canopy_cover(chm.data, 5.0):.0%}; {time.time() - t0:.0f} s") X, Y = dtm.cell_centers() dtm_plot = np.where((X >= PLOT[0]) & (X <= PLOT[2]) & (Y >= PLOT[1]) & (Y <= PLOT[3]), dtm.data, np.nan) fig, axes = plt.subplots(1, 2, figsize=(10, 4.4), sharey=True) for ax, data, r, title, label in [(axes[0], dtm_plot, dtm, "Terrain (DTM)", "elevation (m)"), (axes[1], chm.data, chm, "Canopy height (CHM)", "height above ground (m)")]: lo_, hi_ = np.nanpercentile(data, [1, 99.5]) im = ax.imshow(data, origin="lower", extent=(r.xmin, r.xmax, r.ymin, r.ymax), vmin=lo_, vmax=hi_) ax.set(title=title, aspect="equal", xlim=PLOT[::2], ylim=PLOT[1::2], xlabel="x, project (m)") fig.colorbar(im, ax=ax, shrink=0.9, label=label) axes[0].set_ylabel("y, project (m)");

highest point 93 m, after removing isolated points 23 m
75,337,114 points in the plot; terrain -2.8 to 11.2 m; CHM cells with anything above 5 m: 58%; 93 s

5. Trees¶

Two settings differ from the defaults, both for a plot this size:

• min_arc_deg=130. A stem circle must cover 130 degrees of arc unbroken. With a scan every 20 m, every real stem is seen from most sides, while grass tussocks, termite mounds and shrubs give short arcs. Without it, detection returned over 5 000 candidates on this plot, many of them 1.4 m wide at the diameter limit.
• voxel_size=0.1 for the merge and segmentation graphs. On 70-76 M points the default 5 cm graph took 17 minutes and 38 GB. At 10 cm the whole step takes about 4 minutes and 15 GB.

Some stems come out wider than 40 cm yet lower than 8 m: solid, round and short. On this plot these are most likely termite mounds, or stumps. The table flags them as wide_and_short instead of dropping them; check them in the segmented cloud before using the table.

In [8]:

t0 = time.time()
cands = trees.detect_stems(plot, min_arc_deg=130)
stems, _ = trees.merge_branches(plot, cands, voxel_size=0.1)
labels = trees.segment_trees(plot, stems, voxel_size=0.1)
trees.tree_heights(plot, labels, stems, percentile=99)
stems, labels = trees.prune_trees(stems, labels, min_height=2.0)
crowns = trees.crown_metrics_all(plot, labels)
table = pd.DataFrame([{**t.as_dict(), **crowns.get(t.tree_id, {})} for t in stems])
table["wide_and_short"] = (table.dbh > 0.4) & (table.height < 8)
table.to_csv(OUT / "trees.csv", index=False)
print(f"{len(cands)} candidates -> {len(stems)} trees in {time.time() - t0:.0f} s; "
      f"{table.wide_and_short.sum()} are wider than 40 cm but lower than 8 m")
table[["dbh", "height", "crown_area", "crown_base_height", "quality"]].describe().round(2)

t0 = time.time() cands = trees.detect_stems(plot, min_arc_deg=130) stems, _ = trees.merge_branches(plot, cands, voxel_size=0.1) labels = trees.segment_trees(plot, stems, voxel_size=0.1) trees.tree_heights(plot, labels, stems, percentile=99) stems, labels = trees.prune_trees(stems, labels, min_height=2.0) crowns = trees.crown_metrics_all(plot, labels) table = pd.DataFrame([{**t.as_dict(), **crowns.get(t.tree_id, {})} for t in stems]) table["wide_and_short"] = (table.dbh > 0.4) & (table.height < 8) table.to_csv(OUT / "trees.csv", index=False) print(f"{len(cands)} candidates -> {len(stems)} trees in {time.time() - t0:.0f} s; " f"{table.wide_and_short.sum()} are wider than 40 cm but lower than 8 m") table[["dbh", "height", "crown_area", "crown_base_height", "quality"]].describe().round(2)

782 candidates -> 638 trees in 204 s; 10 are wider than 40 cm but lower than 8 m

Out[8]:

	dbh	height	crown_area	crown_base_height	quality
count	638.00	638.00	638.00	638.00	638.00
mean	0.12	9.24	10.33	4.32	0.47
std	0.10	3.70	16.03	2.79	0.11
min	0.03	2.13	0.03	0.50	0.13
25%	0.07	6.81	2.84	2.00	0.46
50%	0.09	8.94	4.85	4.00	0.52
75%	0.12	11.13	9.04	6.50	0.54
max	1.06	22.21	118.27	13.00	0.60

In [9]:

fig, axes = plt.subplots(1, 3, figsize=(10, 3.6), gridspec_kw={"width_ratios": [1.25, 1, 1]})
order = np.argsort(table.height.to_numpy())
t_ = table.iloc[order]
sc_ = axes[0].scatter(t_.x, t_.y, s=t_.dbh * 150, c=t_.height, alpha=0.8)
fig.colorbar(sc_, ax=axes[0], shrink=0.9, label="tree height (m)")
axes[0].set(aspect="equal", title="Stems (size by DBH)", xlim=PLOT[::2], ylim=PLOT[1::2],
            xlabel="x, project (m)", ylabel="y, project (m)")
axes[1].hist(table.dbh * 100, bins=np.arange(0, 80, 2.5), color="C0")
axes[1].set(xlabel="DBH (cm)", ylabel="trees", title=f"{len(table)} trees, {(table.dbh >= 0.1).sum()} of 10 cm or more")
w = table.wide_and_short
axes[2].scatter(table.dbh[~w] * 100, table.height[~w], s=6, alpha=0.6, c="C0", label="trees")
axes[2].scatter(table.dbh[w] * 100, table.height[w], s=12, c="C3", label="wide and short")
axes[2].set(xlabel="DBH (cm)", ylabel="height (m)", title="Height against DBH")
axes[2].legend(loc="center right")

sub = np.random.default_rng(0).choice(len(plot), 2_000_000, replace=False)
lab = labels[sub]
shuffle = np.random.default_rng(3).permutation(lab.max() + 2)       # neighbouring trees get unlike colours
cols = np.where(lab[:, None] >= 0, LABELS(shuffle[lab] % LABELS.N)[:, :3], 0.85)
fig, ax = plt.subplots(figsize=(6.5, 6.5))
order = np.argsort(plot.attrs["height"][sub])
ax.scatter(plot.x[sub][order], plot.y[sub][order], c=cols[order], s=0.05)
ax.set(aspect="equal", title="Trees from above, one colour each (grey: no tree)",
       xlabel="x, project (m)", ylabel="y, project (m)");

fig, axes = plt.subplots(1, 3, figsize=(10, 3.6), gridspec_kw={"width_ratios": [1.25, 1, 1]}) order = np.argsort(table.height.to_numpy()) t_ = table.iloc[order] sc_ = axes[0].scatter(t_.x, t_.y, s=t_.dbh * 150, c=t_.height, alpha=0.8) fig.colorbar(sc_, ax=axes[0], shrink=0.9, label="tree height (m)") axes[0].set(aspect="equal", title="Stems (size by DBH)", xlim=PLOT[::2], ylim=PLOT[1::2], xlabel="x, project (m)", ylabel="y, project (m)") axes[1].hist(table.dbh * 100, bins=np.arange(0, 80, 2.5), color="C0") axes[1].set(xlabel="DBH (cm)", ylabel="trees", title=f"{len(table)} trees, {(table.dbh >= 0.1).sum()} of 10 cm or more") w = table.wide_and_short axes[2].scatter(table.dbh[~w] * 100, table.height[~w], s=6, alpha=0.6, c="C0", label="trees") axes[2].scatter(table.dbh[w] * 100, table.height[w], s=12, c="C3", label="wide and short") axes[2].set(xlabel="DBH (cm)", ylabel="height (m)", title="Height against DBH") axes[2].legend(loc="center right") sub = np.random.default_rng(0).choice(len(plot), 2_000_000, replace=False) lab = labels[sub] shuffle = np.random.default_rng(3).permutation(lab.max() + 2) # neighbouring trees get unlike colours cols = np.where(lab[:, None] >= 0, LABELS(shuffle[lab] % LABELS.N)[:, :3], 0.85) fig, ax = plt.subplots(figsize=(6.5, 6.5)) order = np.argsort(plot.attrs["height"][sub]) ax.scatter(plot.x[sub][order], plot.y[sub][order], c=cols[order], s=0.05) ax.set(aspect="equal", title="Trees from above, one colour each (grey: no tree)", xlabel="x, project (m)", ylabel="y, project (m)");

6. Scan quality¶

Stems between 1 and 3 m are fitted from all scans together, and each point's radial residual is read per scan: range noise with the stem's shape removed, and each scan's horizontal offset from where the others put the stem. The best-supported stems are enough. summary() leaves scans that were measured in no stem slice out of the registration figures and reports how many were used (n_scans_registered).

In [10]:

t0 = time.time()
good = [t for t in stems if t.dbh >= 0.1 and t.quality >= np.quantile([s.quality for s in stems], 0.5)]
q = quality.stem_noise(plot, scan_id="scan_id", stems=good)
summary = q.summary()
print(f"{time.time() - t0:.0f} s;", {k: (round(v * 1000, 1) if isinstance(v, float) and k != "tail_fraction" else v)
                                   for k, v in summary.items()}, "(sigmas and offsets in mm)")
sc = pd.DataFrame(q.scans)
sc = sc[sc.n_slices >= 1]
fig, ax = plt.subplots(figsize=(5.5, 5.5))
o = origins[sc.scan]
arrows = ax.quiver(o[:, 0], o[:, 1], sc.tx * 1000, sc.ty * 1000, angles="xy", scale_units="xy", scale=0.5, width=0.004)
ax.quiverkey(arrows, 0.8, 1.02, 5, "5 mm", labelpos="E")
ax.add_patch(plt.Rectangle(PLOT[:2], 100, 100, fill=False, lw=1))
ax.set(aspect="equal", title="Horizontal offset of each scan, from the stems\n(1 mm drawn as 2 m)",
       xlabel="x, project (m)", ylabel="y, project (m)");

t0 = time.time() good = [t for t in stems if t.dbh >= 0.1 and t.quality >= np.quantile([s.quality for s in stems], 0.5)] q = quality.stem_noise(plot, scan_id="scan_id", stems=good) summary = q.summary() print(f"{time.time() - t0:.0f} s;", {k: (round(v * 1000, 1) if isinstance(v, float) and k != "tail_fraction" else v) for k, v in summary.items()}, "(sigmas and offsets in mm)") sc = pd.DataFrame(q.scans) sc = sc[sc.n_slices >= 1] fig, ax = plt.subplots(figsize=(5.5, 5.5)) o = origins[sc.scan] arrows = ax.quiver(o[:, 0], o[:, 1], sc.tx * 1000, sc.ty * 1000, angles="xy", scale_units="xy", scale=0.5, width=0.004) ax.quiverkey(arrows, 0.8, 1.02, 5, "5 mm", labelpos="E") ax.add_patch(plt.Rectangle(PLOT[:2], 100, 100, fill=False, lw=1)) ax.set(aspect="equal", title="Horizontal offset of each scan, from the stems\n(1 mm drawn as 2 m)", xlabel="x, project (m)", ylabel="y, project (m)");

5 s; {'n_stems': 139, 'n_slices': 820, 'n_scans': 63, 'sigma_total': 7.7, 'sigma_corrected': 6.9, 'sigma_within': 6.1, 'sigma_local': 4.9, 'tail_fraction': 0.02919919850756581, 'n_scans_registered': 40, 'registration_rms': 5.4, 'registration_max': 16.6, 'worst_scan': 35} (sigmas and offsets in mm)

7. Wood models and leaves¶

For the largest well-supported trees: the wood filter, a cylinder model and its architecture. measured_volume_fraction is the share of the volume fitted to points rather than filled in by the taper and pipe-model priors, so read the volumes with it.

In [11]:

t0 = time.time()
for f in list(OUT.glob("qsm_[0-9]*.*")) + list(OUT.glob("tree_*.obj")):
    f.unlink()                                          # models of an earlier run
q_med = np.median([t.quality for t in stems])
big = sorted([t for t in stems if t.quality >= q_med and t.height >= 10], key=lambda t: -t.dbh)[:4]  # largest real trees
models, rows = {}, []
for t in big:
    tree = plot[labels == t.tree_id]
    model = qsm.build_qsm(qsm.wood_points(tree), base_xy=(t.x, t.y))
    model.to_csv(OUT / f"qsm_{t.tree_id:04d}.csv")
    model.to_ply(OUT / f"qsm_{t.tree_id:04d}.ply")
    m = model.metrics()
    models[t.tree_id] = (tree, model)
    rows.append({"tree_id": t.tree_id, "dbh_stem": t.dbh, "dbh_qsm": m["dbh"], "height": m["height"],
                 "volume_m3": m["total_volume"], "measured_volume": m["measured_volume_fraction"],
                 "max_order": m["max_order"], "insertion_deg": m["median_insertion_angle"],
                 "crown_area": m["crown"]["projected_area"]})
print(f"{len(big)} QSMs in {time.time() - t0:.0f} s")
pd.DataFrame(rows).round(3)

t0 = time.time() for f in list(OUT.glob("qsm_[0-9]*.*")) + list(OUT.glob("tree_*.obj")): f.unlink() # models of an earlier run q_med = np.median([t.quality for t in stems]) big = sorted([t for t in stems if t.quality >= q_med and t.height >= 10], key=lambda t: -t.dbh)[:4] # largest real trees models, rows = {}, [] for t in big: tree = plot[labels == t.tree_id] model = qsm.build_qsm(qsm.wood_points(tree), base_xy=(t.x, t.y)) model.to_csv(OUT / f"qsm_{t.tree_id:04d}.csv") model.to_ply(OUT / f"qsm_{t.tree_id:04d}.ply") m = model.metrics() models[t.tree_id] = (tree, model) rows.append({"tree_id": t.tree_id, "dbh_stem": t.dbh, "dbh_qsm": m["dbh"], "height": m["height"], "volume_m3": m["total_volume"], "measured_volume": m["measured_volume_fraction"], "max_order": m["max_order"], "insertion_deg": m["median_insertion_angle"], "crown_area": m["crown"]["projected_area"]}) print(f"{len(big)} QSMs in {time.time() - t0:.0f} s") pd.DataFrame(rows).round(3)

4 QSMs in 10 s

Out[11]:

	tree_id	dbh_stem	dbh_qsm	height	volume_m3	measured_volume	max_order	insertion_deg	crown_area
0	9	0.480	0.482	21.466	3.113	0.684	7	45.284	65.487
1	13	0.404	0.396	20.688	2.079	0.781	6	53.313	91.278
2	14	0.385	0.384	22.469	1.989	0.761	7	52.269	91.279
3	18	0.351	0.330	19.323	1.976	0.630	8	55.520	111.510

In [12]:

from matplotlib.collections import PolyCollection


def side_view(model, axis=0):
    # Each cylinder as its true-width outline seen from the side (axis 0: x-z, 1: y-z).
    a, b = model.start[:, [axis, 2]], model.end[:, [axis, 2]]
    d = b - a
    n = np.c_[-d[:, 1], d[:, 0]] / np.maximum(np.hypot(*d.T), 1e-9)[:, None] * model.column("radius")[:, None]
    return PolyCollection(np.stack([a + n, b + n, b - n, a - n], axis=1), facecolor=WOOD, edgecolor=WOOD, lw=0.3)


tid = big[0].tree_id
tree, model = models[tid]
s = tree[np.random.default_rng(0).choice(len(tree), min(len(tree), 80_000), replace=False)]
fig, axes = plt.subplots(1, 3, figsize=(10, 5.5), sharey=True)
axes[0].scatter(s.x, s.z, s=0.1, c="0.5")
axes[0].set(title=f"Tree {tid}: points, x-z", xlabel="x, project (m)", ylabel="z (m)")
for ax, k, lab in [(axes[1], 0, "x"), (axes[2], 1, "y")]:
    ax.scatter(s.xyz[:, k], s.z, s=0.1, c=CONTEXT)
    ax.add_collection(side_view(model, k))
    ax.set(title=f"QSM, {lab}-z ({len(model)} cylinders)", xlabel=f"{lab}, project (m)")
for ax in axes:
    ax.set_aspect("equal")
    ax.autoscale_view()

from matplotlib.collections import PolyCollection def side_view(model, axis=0): # Each cylinder as its true-width outline seen from the side (axis 0: x-z, 1: y-z). a, b = model.start[:, [axis, 2]], model.end[:, [axis, 2]] d = b - a n = np.c_[-d[:, 1], d[:, 0]] / np.maximum(np.hypot(*d.T), 1e-9)[:, None] * model.column("radius")[:, None] return PolyCollection(np.stack([a + n, b + n, b - n, a - n], axis=1), facecolor=WOOD, edgecolor=WOOD, lw=0.3) tid = big[0].tree_id tree, model = models[tid] s = tree[np.random.default_rng(0).choice(len(tree), min(len(tree), 80_000), replace=False)] fig, axes = plt.subplots(1, 3, figsize=(10, 5.5), sharey=True) axes[0].scatter(s.x, s.z, s=0.1, c="0.5") axes[0].set(title=f"Tree {tid}: points, x-z", xlabel="x, project (m)", ylabel="z (m)") for ax, k, lab in [(axes[1], 0, "x"), (axes[2], 1, "y")]: ax.scatter(s.xyz[:, k], s.z, s=0.1, c=CONTEXT) ax.add_collection(side_view(model, k)) ax.set(title=f"QSM, {lab}-z ({len(model)} cylinders)", xlabel=f"{lab}, project (m)") for ax in axes: ax.set_aspect("equal") ax.autoscale_view()

In [13]:

wood = leaves.classify_leaf_wood(tree)
foliage = tree[~wood]
angles = leaves.leaf_angle_distribution(foliage)
area = leaves.leaf_area_density(foliage, voxel_size=0.25)
mesh = leaves.add_leaves(model, area, angles, leaf_points=foliage)
leaves.write_tree_obj(OUT / f"tree_{tid:04d}.obj", model, mesh)
print(f"wood {wood.mean():.0%} of points; mean leaf angle {angles.mean_deg:.0f} deg (nearest de Wit: {angles.de_wit}); "
      f"leaf area seen {area.total_area:.1f} m2 as {len(mesh)} leaves -> tree_{tid:04d}.obj")

wood = leaves.classify_leaf_wood(tree) foliage = tree[~wood] angles = leaves.leaf_angle_distribution(foliage) area = leaves.leaf_area_density(foliage, voxel_size=0.25) mesh = leaves.add_leaves(model, area, angles, leaf_points=foliage) leaves.write_tree_obj(OUT / f"tree_{tid:04d}.obj", model, mesh) print(f"wood {wood.mean():.0%} of points; mean leaf angle {angles.mean_deg:.0f} deg (nearest de Wit: {angles.de_wit}); " f"leaf area seen {area.total_area:.1f} m2 as {len(mesh)} leaves -> tree_{tid:04d}.obj")

wood 42% of points; mean leaf angle 57 deg (nearest de Wit: spherical); leaf area seen 242.4 m2 as 134773 leaves -> tree_0009.obj

8. Canopy gap profile¶

Gap probability by zenith ring and height, pooled over the 36 scans inside the plot (the ring looks at the plot from outside, across its edge) (Jupp et al. 2009). The VZ-2000i scans from 30 deg zenith down, so the rings run from 30 to 70 deg. The pulses already hold their misses, so no fired-pulse count is needed. Heights come from the DTM.

In [14]:

t0 = time.time()
prof = canopy.GapProfile.empty(zenith_edges=np.arange(30.0, 75.0, 5.0))
start = np.r_[0, np.cumsum(used_counts)]
xyz = shots.echo_xyz()
echo_h = xyz[:, 2] - dtm.sample(xyz[:, 0], xyz[:, 1])
del xyz
soe = shots.shot_of_echo()
for j in np.flatnonzero(inner_used):
    keep = np.zeros(shots.n_shots, bool)
    keep[start[j]:start[j + 1]] = True
    prof.add_scan(shots.subset(keep), echo_h[keep[soe]])
del soe, echo_h
rep = prof.report()
print(f"{rep['n_scans']} scans, {time.time() - t0:.0f} s")
print({k: round(v, 3) if isinstance(v, float) else v for k, v in rep.items() if np.isscalar(v)})

fig, axes = plt.subplots(1, 2, figsize=(8, 4.2), sharey=True)
axes[0].plot(rep["pai_hinge_profile"], rep["height"], label="hinge")
axes[0].plot(rep["pai_linear_profile"], rep["height"], "--", label="linear")
axes[0].set(xlabel="cumulative PAI below the height", ylabel="height above ground (m)",
            ylim=(0, rep["canopy_height"] + 3), title="Plant area index")
axes[0].legend(loc="lower right")
axes[1].plot(rep["pavd_hinge"], rep["height"])
axes[1].set(xlabel="PAVD (m² m⁻³)", title="Plant area volume density, hinge");

t0 = time.time() prof = canopy.GapProfile.empty(zenith_edges=np.arange(30.0, 75.0, 5.0)) start = np.r_[0, np.cumsum(used_counts)] xyz = shots.echo_xyz() echo_h = xyz[:, 2] - dtm.sample(xyz[:, 0], xyz[:, 1]) del xyz soe = shots.shot_of_echo() for j in np.flatnonzero(inner_used): keep = np.zeros(shots.n_shots, bool) keep[start[j]:start[j + 1]] = True prof.add_scan(shots.subset(keep), echo_h[keep[soe]]) del soe, echo_h rep = prof.report() print(f"{rep['n_scans']} scans, {time.time() - t0:.0f} s") print({k: round(v, 3) if isinstance(v, float) else v for k, v in rep.items() if np.isscalar(v)}) fig, axes = plt.subplots(1, 2, figsize=(8, 4.2), sharey=True) axes[0].plot(rep["pai_hinge_profile"], rep["height"], label="hinge") axes[0].plot(rep["pai_linear_profile"], rep["height"], "--", label="linear") axes[0].set(xlabel="cumulative PAI below the height", ylabel="height above ground (m)", ylim=(0, rep["canopy_height"] + 3), title="Plant area index") axes[0].legend(loc="lower right") axes[1].plot(rep["pavd_hinge"], rep["height"]) axes[1].set(xlabel="PAVD (m² m⁻³)", title="Plant area volume density, hinge");

36 scans, 95 s
{'saturated': False, 'gap_57': 0.455, 'pai_hinge': 0.866, 'pai_linear': 0.768, 'pai_weighted': 0.8, 'mla_linear': 70.819, 'clumping': 0.93, 'pai_hinge_corrected': 0.931, 'canopy_height': 20.0, 'closure_57': 0.545, 'cover': 0.354, 'cover_zenith': 32.5, 'n_scans': 36, 'pulses': 17710145.0}

9. Ray-traced voxels, and what was seen¶

Every pulse is traced through a 0.5 m grid over the plot. With occlusion=True each voxel is observed (a pulse went through or ended in it), occluded (only pulses already stopped reached it) or unreached.

• Beam geometry. beam= weights each pulse by its cross-section in the voxel, which grows with range. The values are approximate (7 mm exit, 0.27 mrad). Only relative sections enter a pooled layer estimate, so they matter little: without beam the PAI differs by about 3 %.
• Profile. Plant area density per height layer is pooled over the layer (intercepted section over effective free path), and only layers where the median voxel saw at least 50 pulses are counted. With every 32nd pulse, that is where the numbers can be trusted.

The voxel PAI comes out close to the gap profile's hinge PAI (0.88 against 0.87, or 0.93 with the clumping correction). They are independent estimates from the same pulses: one traces every pulse through the grid, the other counts gaps by zenith ring. Both depend on the misses being there. Traced without them, the grid gave about twice the plant area, rising with height.

With 64 positions and the misses traced, a pulse crosses every voxel of this open canopy, so nothing is occluded or unreached. In denser plots, and with fewer scans, occlusion_profile and observed_map show where to rescan. The right panel shows how many pulses reached each layer, which sets how far up the profile can be trusted.

In [15]:

t0 = time.time()
bounds = ((PLOT[0], PLOT[1], float(np.nanmin(dtm.data)) - 1.0), (PLOT[2], PLOT[3], float(np.nanmax(dtm.data)) + 30.0))
grid = voxels.ray_voxelize(rays_path, 0.5, bounds, dtm=dtm, occlusion=True, beam=(0.007, 0.00027))
occ = grid.occlusion_profile(min_height=0.5)
print(f"{grid} in {time.time() - t0:.0f} s; canopy space:",
      {k: round(v, 3) for k, v in occ["total"].items()})

h, beams = grid.distance_from_ground, grid.num_beams
edges = np.arange(0.0, 30.5, 0.5)
k = np.digitize(h, edges) - 1
layer = []
for j in range(len(edges) - 1):
    m = k == j
    bi, be = grid.bs_intercepted[m].sum(), grid.bs_effective_free_path[m].sum()
    layer.append((edges[j], np.median(beams[m]) if m.any() else 0, 2 * bi / be if be > 0 else np.nan))
layer = pd.DataFrame(layer, columns=["height", "median_beams", "pad"])
ok = (layer.median_beams >= 50) & (layer.height >= 0.5)
top = layer.height[ok].max() + 0.5
print(f"voxel PAI from 0.5 to {top:.1f} m (sampled layers): {layer.pad[ok].sum() * 0.5:.2f}; "
      f"gap profile hinge PAI {rep['pai_hinge']:.2f} (clumping-corrected {rep['pai_hinge_corrected']:.2f})")

fig, axes = plt.subplots(1, 2, figsize=(8, 4.2), sharey=True)
for ax in axes:
    ax.axhspan(top, 30.5, color=CONTEXT, alpha=0.5, lw=0)          # layers with too few pulses
axes[0].plot(layer.pad.where(ok), layer.height + 0.25, "o-", ms=3, label="layers used")
axes[0].plot(layer.pad.where(~ok), layer.height + 0.25, "o", ms=3, color="0.5", label="not counted")
axes[0].set(xlabel="PAD (m² m⁻³)", ylabel="height above ground (m)", ylim=(0, 30.5), title="Voxel plant area density")
axes[0].legend(loc="upper right")
axes[1].semilogx(layer.median_beams.clip(lower=1), layer.height + 0.25)
axes[1].axvline(50, color="k", ls="--", lw=0.8)
axes[1].text(55, 1.0, "50 pulses", fontsize=8.5)
axes[1].set(xlabel="median pulses per voxel", xlim=(10, None), title="Sampling of each layer\n(grey: median below 50)");

t0 = time.time() bounds = ((PLOT[0], PLOT[1], float(np.nanmin(dtm.data)) - 1.0), (PLOT[2], PLOT[3], float(np.nanmax(dtm.data)) + 30.0)) grid = voxels.ray_voxelize(rays_path, 0.5, bounds, dtm=dtm, occlusion=True, beam=(0.007, 0.00027)) occ = grid.occlusion_profile(min_height=0.5) print(f"{grid} in {time.time() - t0:.0f} s; canopy space:", {k: round(v, 3) for k, v in occ["total"].items()}) h, beams = grid.distance_from_ground, grid.num_beams edges = np.arange(0.0, 30.5, 0.5) k = np.digitize(h, edges) - 1 layer = [] for j in range(len(edges) - 1): m = k == j bi, be = grid.bs_intercepted[m].sum(), grid.bs_effective_free_path[m].sum() layer.append((edges[j], np.median(beams[m]) if m.any() else 0, 2 * bi / be if be > 0 else np.nan)) layer = pd.DataFrame(layer, columns=["height", "median_beams", "pad"]) ok = (layer.median_beams >= 50) & (layer.height >= 0.5) top = layer.height[ok].max() + 0.5 print(f"voxel PAI from 0.5 to {top:.1f} m (sampled layers): {layer.pad[ok].sum() * 0.5:.2f}; " f"gap profile hinge PAI {rep['pai_hinge']:.2f} (clumping-corrected {rep['pai_hinge_corrected']:.2f})") fig, axes = plt.subplots(1, 2, figsize=(8, 4.2), sharey=True) for ax in axes: ax.axhspan(top, 30.5, color=CONTEXT, alpha=0.5, lw=0) # layers with too few pulses axes[0].plot(layer.pad.where(ok), layer.height + 0.25, "o-", ms=3, label="layers used") axes[0].plot(layer.pad.where(~ok), layer.height + 0.25, "o", ms=3, color="0.5", label="not counted") axes[0].set(xlabel="PAD (m² m⁻³)", ylabel="height above ground (m)", ylim=(0, 30.5), title="Voxel plant area density") axes[0].legend(loc="upper right") axes[1].semilogx(layer.median_beams.clip(lower=1), layer.height + 0.25) axes[1].axvline(50, color="k", ls="--", lw=0.8) axes[1].text(55, 1.0, "50 pulses", fontsize=8.5) axes[1].set(xlabel="median pulses per voxel", xlim=(10, None), title="Sampling of each layer\n(grey: median below 50)");

RayVoxelGrid(200x200x90 @ 0.5 m) in 76 s; canopy space: {'observed': 1.0, 'occluded': 0.0, 'unobserved': 0.0, 'top': 29.203}

voxel PAI from 0.5 to 30.0 m (sampled layers): 0.88; gap profile hinge PAI 0.87 (clumping-corrected 0.93)

Per tree, the voxels show whether the space just above the tree was seen. Where it was not, the tree may continue where no scan reached, and its height is a lower bound. Here every tree top was seen, so the heights in the table are not limited by occlusion. How many pulses reached each crown is the more telling figure: it sets how well crown shape and leaf area are known.

In [16]:

samp = voxels.tree_sampling(grid, plot, labels)
samp = pd.DataFrame({k: v for k, v in samp.items() if np.ndim(v) == 1})     # beams_by_quarter is (trees, 4)
table = table.merge(samp[["tree_id", "above_observed_fraction", "median_beams"]], on="tree_id", how="left")
table.to_csv(OUT / "trees.csv", index=False)
flag = table.above_observed_fraction < 0.5
print(f"{flag.sum()} of {len(table)} trees have less than half of the space above them observed")
fig, ax = plt.subplots()
ax.scatter(table.height, table.median_beams, s=8, alpha=0.6, c="C0")
ax.set(xlabel="tree height (m)", ylabel="median pulses per crown voxel", yscale="log",
       title="Sampling of each crown (every 32nd pulse traced)");

samp = voxels.tree_sampling(grid, plot, labels) samp = pd.DataFrame({k: v for k, v in samp.items() if np.ndim(v) == 1}) # beams_by_quarter is (trees, 4) table = table.merge(samp[["tree_id", "above_observed_fraction", "median_beams"]], on="tree_id", how="left") table.to_csv(OUT / "trees.csv", index=False) flag = table.above_observed_fraction < 0.5 print(f"{flag.sum()} of {len(table)} trees have less than half of the space above them observed") fig, ax = plt.subplots() ax.scatter(table.height, table.median_beams, s=8, alpha=0.6, c="C0") ax.set(xlabel="tree height (m)", ylabel="median pulses per crown voxel", yscale="log", title="Sampling of each crown (every 32nd pulse traced)");

0 of 638 trees have less than half of the space above them observed

Outputs¶

Everything lands in OUT:

In [17]:

sylva.write(plot.with_attrs(tree_id=labels.astype(np.int32)), OUT / "plot_segmented.laz")
for f in sorted(OUT.iterdir()):
    print(f"{f.name:28s} {f.stat().st_size / 1e6:9.1f} MB")

sylva.write(plot.with_attrs(tree_id=labels.astype(np.int32)), OUT / "plot_segmented.laz") for f in sorted(OUT.iterdir()): print(f"{f.name:28s} {f.stat().st_size / 1e6:9.1f} MB")

chm.asc                            0.3 MB
cloud_2cm.laz                    814.4 MB
coreg_features.pkl               548.3 MB
dtm.asc                            0.5 MB
plot_segmented.laz              1081.4 MB
qsm_0009.csv                       0.6 MB
qsm_0009.ply                       5.4 MB
qsm_0013.csv                       0.5 MB
qsm_0013.ply                       4.4 MB
qsm_0014.csv                       0.7 MB
qsm_0014.ply                       6.9 MB
qsm_0018.csv                       1.3 MB
qsm_0018.ply                      12.2 MB
ray_counts.npy                     0.0 MB
rays.parquet                     747.8 MB
rays_registered.parquet          747.8 MB
sop_corrected.json                 0.0 MB
tree_0009.obj                     42.7 MB
trees.csv                          0.2 MB