9. Ray-traced voxels¶

sylva.voxels traces every pulse through a voxel grid and estimates the attenuation coefficient λ of each voxel from how far pulses got, then plant area density as λ / G. It follows Vincent et al. 2017 and the rayvoxel tool; see the guide for the estimators.

In [1]:

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import sylva
from sylva import synthetic

DATA = Path("data")          # the Litchfield tile, cut by make_litch_subset.py
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

from sylva import Shots, voxels

from pathlib import Path import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap import sylva from sylva import synthetic DATA = Path("data") # the Litchfield tile, cut by make_litch_subset.py 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 from sylva import Shots, voxels

Tracing a real tile¶

The Litchfield pulses, with the echoes' classification so the terrain can be left out of the plant area. occlusion=True also records what the pulses never reached.

In [2]:

shots = Shots.load(DATA / "litch_tile_shots.parquet")
grid = voxels.ray_voxelize(
    shots, voxel_size=0.5, bounds=((0, 0, 0), (20, 20, 18)),
    ground_class=2, leaf_classes=[4], attenuation=["fpl", "ppl"], occlusion=True,
)
print(grid)
print("raw fields:", ", ".join(grid.fields[:10]), "...")

shots = Shots.load(DATA / "litch_tile_shots.parquet") grid = voxels.ray_voxelize( shots, voxel_size=0.5, bounds=((0, 0, 0), (20, 20, 18)), ground_class=2, leaf_classes=[4], attenuation=["fpl", "ppl"], occlusion=True, ) print(grid) print("raw fields:", ", ".join(grid.fields[:10]), "...")

RayVoxelGrid(40x40x36 @ 0.5 m)
raw fields: num_hits, num_beams, num_beams_occluded, num_unbound_rays, num_miss_rays, num_hit_leaf, num_hit_wood, num_hit_plant, num_beams_weighted, path_length ...

state is 0 unobserved, 1 occluded (only reached behind a last echo), 2 empty, 3 filled. Arrays are (nz, ny, nx). Here 97 % of the tile was observed and 3 % is occluded, which is what many scan positions buy you.

In [3]:

state = grid.state
names = ["unobserved", "occluded", "empty", "filled"]
print({n: f"{np.mean(state == i):.1%}" for i, n in enumerate(names)})
print({k: round(v, 3) for k, v in grid.occlusion_profile()["total"].items()})
from matplotlib.colors import LogNorm

j = 15                                        # the row of voxels at y = 7.5 to 8 m
ext = (grid.origin[0], grid.origin[0] + grid.voxel_size * state.shape[2],
       grid.origin[2], grid.origin[2] + grid.voxel_size * state.shape[0])
state_colours = ListedColormap(["white", "0.45", "#CFE3F0", LEAF])      # unobserved, occluded, empty, filled
fig, ax = plt.subplots(1, 3, figsize=(10, 3.3), sharey=True)
im = ax[0].imshow(state[:, j, :], origin="lower", extent=ext, cmap=state_colours, vmin=-0.5, vmax=3.5)
cb = fig.colorbar(im, ax=ax[0], ticks=range(4), shrink=0.9)
cb.ax.set_yticklabels(names)
im = ax[1].imshow(grid.num_beams[:, j, :], origin="lower", extent=ext, norm=LogNorm(1, 1e4))
fig.colorbar(im, ax=ax[1], shrink=0.9, label="pulses through the voxel")
im = ax[2].imshow(grid.num_hits[:, j, :], origin="lower", extent=ext, vmax=50)
fig.colorbar(im, ax=ax[2], shrink=0.9, label="echoes in the voxel", extend="max")
for a, title in zip(ax, ("Voxel state", "Pulses", "Echoes")):
    a.set(title=title, xlabel="x (m)")
ax[0].set_ylabel("z (m)")
fig.suptitle("A vertical slice of 0.5 m voxels at y = 7.5 to 8 m");

state = grid.state names = ["unobserved", "occluded", "empty", "filled"] print({n: f"{np.mean(state == i):.1%}" for i, n in enumerate(names)}) print({k: round(v, 3) for k, v in grid.occlusion_profile()["total"].items()}) from matplotlib.colors import LogNorm j = 15 # the row of voxels at y = 7.5 to 8 m ext = (grid.origin[0], grid.origin[0] + grid.voxel_size * state.shape[2], grid.origin[2], grid.origin[2] + grid.voxel_size * state.shape[0]) state_colours = ListedColormap(["white", "0.45", "#CFE3F0", LEAF]) # unobserved, occluded, empty, filled fig, ax = plt.subplots(1, 3, figsize=(10, 3.3), sharey=True) im = ax[0].imshow(state[:, j, :], origin="lower", extent=ext, cmap=state_colours, vmin=-0.5, vmax=3.5) cb = fig.colorbar(im, ax=ax[0], ticks=range(4), shrink=0.9) cb.ax.set_yticklabels(names) im = ax[1].imshow(grid.num_beams[:, j, :], origin="lower", extent=ext, norm=LogNorm(1, 1e4)) fig.colorbar(im, ax=ax[1], shrink=0.9, label="pulses through the voxel") im = ax[2].imshow(grid.num_hits[:, j, :], origin="lower", extent=ext, vmax=50) fig.colorbar(im, ax=ax[2], shrink=0.9, label="echoes in the voxel", extend="max") for a, title in zip(ax, ("Voxel state", "Pulses", "Echoes")): a.set(title=title, xlabel="x (m)") ax[0].set_ylabel("z (m)") fig.suptitle("A vertical slice of 0.5 m voxels at y = 7.5 to 8 m");

{'unobserved': '0.0%', 'occluded': '2.9%', 'empty': '80.1%', 'filled': '17.0%'}
{'observed': 0.971, 'occluded': 0.029, 'unobserved': 0.0, 'top': 17.75}

Attenuation, and where the tile can support an estimate¶

FPL and PPL agree where voxels are well sampled and diverge where they are not, which is the useful diagnostic. Before reading any profile, look at how many beams reached each layer: the tile holds one real scan position (its corner) plus rays from the rest of the plot clipped at the boundary, so sampling falls away with height. Layers whose voxels see only a few tens of beams produce large, meaningless λ -- a voxel with 20 beams, 2 echoes and short path lengths reads as dense -- and profile averages over whichever voxels pass min_beams, so a stricter threshold can raise the number rather than settle it.

In [4]:

beams = np.median(grid.num_beams, axis=(1, 2))
z = grid.z_levels() + 0.25
usable = z[beams >= 200].max()
fig, ax = plt.subplots(1, 3, figsize=(10, 3.6))
well = (grid.num_beams >= 200) & (grid.num_hits > 0)
ax[0].loglog(grid.attenuation_fpl[well], grid.attenuation_ppl[well], ".", ms=1.5, alpha=0.3, c="C0")
ax[0].plot([1e-3, 40], [1e-3, 40], "k", lw=0.8, label="1:1")
ax[0].set(xlabel="λ, FPL (m⁻¹)", ylabel="λ, PPL (m⁻¹)", title="The two estimators, per voxel\n(voxels with ≥ 200 pulses)")
ax[0].legend(loc="upper left")
for a in ax[1:]:
    a.axhspan(usable, z.max() + 0.25, color=CONTEXT, alpha=0.5, lw=0)
ax[1].plot(beams, z, "C0")
ax[1].axvline(200, color="k", ls="--", lw=0.8)
ax[1].text(230, 0.5, "200", fontsize=8.5)
ax[1].set(xscale="log", xlabel="median pulses per voxel", ylabel="z (m)", title="Sampling by layer\n(grey: median below 200)")
for mb, style in ((20, "-"), (200, "--")):
    ax[2].plot(grid.profile("pad_ppl", min_beams=mb), z, style, label=f"min_beams = {mb}")
ax[2].set(xlabel="PAD (m² m⁻³)", ylabel="z (m)", title="The PPL profile, two thresholds"); ax[2].legend(loc="upper right")
print(f"median beams per voxel stays above 200 up to {usable:.1f} m, and falls below 20 above "
      f"{z[beams >= 20].max():.1f} m")

beams = np.median(grid.num_beams, axis=(1, 2)) z = grid.z_levels() + 0.25 usable = z[beams >= 200].max() fig, ax = plt.subplots(1, 3, figsize=(10, 3.6)) well = (grid.num_beams >= 200) & (grid.num_hits > 0) ax[0].loglog(grid.attenuation_fpl[well], grid.attenuation_ppl[well], ".", ms=1.5, alpha=0.3, c="C0") ax[0].plot([1e-3, 40], [1e-3, 40], "k", lw=0.8, label="1:1") ax[0].set(xlabel="λ, FPL (m⁻¹)", ylabel="λ, PPL (m⁻¹)", title="The two estimators, per voxel\n(voxels with ≥ 200 pulses)") ax[0].legend(loc="upper left") for a in ax[1:]: a.axhspan(usable, z.max() + 0.25, color=CONTEXT, alpha=0.5, lw=0) ax[1].plot(beams, z, "C0") ax[1].axvline(200, color="k", ls="--", lw=0.8) ax[1].text(230, 0.5, "200", fontsize=8.5) ax[1].set(xscale="log", xlabel="median pulses per voxel", ylabel="z (m)", title="Sampling by layer\n(grey: median below 200)") for mb, style in ((20, "-"), (200, "--")): ax[2].plot(grid.profile("pad_ppl", min_beams=mb), z, style, label=f"min_beams = {mb}") ax[2].set(xlabel="PAD (m² m⁻³)", ylabel="z (m)", title="The PPL profile, two thresholds"); ax[2].legend(loc="upper right") print(f"median beams per voxel stays above 200 up to {usable:.1f} m, and falls below 20 above " f"{z[beams >= 20].max():.1f} m")

median beams per voxel stays above 200 up to 11.2 m, and falls below 20 above 16.2 m

Read the profile up to that height and no further. Summed over the well-sampled layers it comes to about 0.7 -- the grass layer and the lower canopy, and in the same range as the plot's hinge PAI of 0.87. Summed over every layer it comes to 6.3, which is the sparsely sampled top of the tile inventing plant area. A 20 x 20 m cut-out cannot give a plot's plant area index either way: its rays were clipped at the boundary and the upper canopy is barely sampled. The canopy benchmark has the plot values from whole scans, where Sylva's profiles match pylidar-tls-canopy to within 4 %.

The slice below shows the same voxel by voxel: most of the dense-looking voxels in the upper canopy are ones that fewer than 200 pulses reached.

In [5]:

ok = beams >= 200
print("PAI over the well-sampled layers only:",
      round(float(np.nansum(grid.profile("pad_ppl", min_beams=200)[ok]) * grid.voxel_size), 2),
      " over every layer:",
      round(float(np.nansum(grid.profile("pad_ppl", min_beams=200)) * grid.voxel_size), 2))
pad_slice = grid.pad_ppl[:, j, :]
fig, ax = plt.subplots(1, 2, figsize=(10, 3.8), sharey=True)
ax[0].imshow(pad_slice, origin="lower", extent=ext, vmin=0, vmax=3)
ax[1].imshow(np.where(grid.num_beams[:, j, :] >= 200, np.nan, 1.0), origin="lower", extent=ext,
             cmap=ListedColormap([CONTEXT]))
im = ax[1].imshow(np.where(grid.num_beams[:, j, :] >= 200, pad_slice, np.nan), origin="lower", extent=ext, vmin=0, vmax=3)
fig.colorbar(im, ax=ax, shrink=0.9, label="PAD, PPL (m² m⁻³)", extend="max")
ax[0].set(title="Every voxel", xlabel="x (m)", ylabel="z (m)")
ax[1].set(title="Voxels with ≥ 200 pulses (grey: fewer)", xlabel="x (m)")
fig.suptitle("Plant area density in the slice at y = 7.5 to 8 m");

ok = beams >= 200 print("PAI over the well-sampled layers only:", round(float(np.nansum(grid.profile("pad_ppl", min_beams=200)[ok]) * grid.voxel_size), 2), " over every layer:", round(float(np.nansum(grid.profile("pad_ppl", min_beams=200)) * grid.voxel_size), 2)) pad_slice = grid.pad_ppl[:, j, :] fig, ax = plt.subplots(1, 2, figsize=(10, 3.8), sharey=True) ax[0].imshow(pad_slice, origin="lower", extent=ext, vmin=0, vmax=3) ax[1].imshow(np.where(grid.num_beams[:, j, :] >= 200, np.nan, 1.0), origin="lower", extent=ext, cmap=ListedColormap([CONTEXT])) im = ax[1].imshow(np.where(grid.num_beams[:, j, :] >= 200, pad_slice, np.nan), origin="lower", extent=ext, vmin=0, vmax=3) fig.colorbar(im, ax=ax, shrink=0.9, label="PAD, PPL (m² m⁻³)", extend="max") ax[0].set(title="Every voxel", xlabel="x (m)", ylabel="z (m)") ax[1].set(title="Voxels with ≥ 200 pulses (grey: fewer)", xlabel="x (m)") fig.suptitle("Plant area density in the slice at y = 7.5 to 8 m");

PAI over the well-sampled layers only: 0.74  over every layer: 6.31

Checking the estimators against a known scene¶

To see whether the numbers are right, the scene has to be known. A synthetic plot (synthetic.plot, 20 m square: 20 broadleaf and eucalypt trees with shrubs, grass and dead wood on rough ground) is scanned from five positions with the model of a RIEGL VZ-400, its beam footprint, 5 mm of range noise and mixed pixels. The echoes keep their true class, so the leaf area density can be estimated from the leaf echoes alone and compared with the leaf area the generator placed.

In [6]:

p = synthetic.plot(size=20, density=500, archetypes={"broadleaf": 1, "eucalypt": 1}, seed=1)
positions = [(4, 4), (16, 4), (4, 16), (16, 16), (10, 10)]
sim = Shots.concatenate([
    synthetic.scan(p.points, origin=(x, y, p.ground_height(x, y) + 1.5), resolution_deg=0.1,
                   scanner="vz400", range_noise=0.005, seed=k)
    for k, (x, y) in enumerate(positions)])
box = ((0, 0, float(p.points.z.min()) - 0.5), (20, 20, float(p.points.z.min()) + 29.5))
sim_grid = voxels.ray_voxelize(sim, 0.5, box, ground_class=2, leaf_classes=[4], wood_classes=[5],
                               attenuation=["fpl", "ppl"])
pai = float(np.nansum(sim_grid.profile("pad_ppl", min_beams=20)) * 0.5)
lai = float(np.nansum(sim_grid.profile("lad_ppl", min_beams=20)) * 0.5)
print(sim)
print(f"PAI {pai:.2f}  LAI {lai:.2f}  true LAI of the trees {p.trees['leaf_area'].sum() / 20**2:.2f}")

p = synthetic.plot(size=20, density=500, archetypes={"broadleaf": 1, "eucalypt": 1}, seed=1) positions = [(4, 4), (16, 4), (4, 16), (16, 16), (10, 10)] sim = Shots.concatenate([ synthetic.scan(p.points, origin=(x, y, p.ground_height(x, y) + 1.5), resolution_deg=0.1, scanner="vz400", range_noise=0.005, seed=k) for k, (x, y) in enumerate(positions)]) box = ((0, 0, float(p.points.z.min()) - 0.5), (20, 20, float(p.points.z.min()) + 29.5)) sim_grid = voxels.ray_voxelize(sim, 0.5, box, ground_class=2, leaf_classes=[4], wood_classes=[5], attenuation=["fpl", "ppl"]) pai = float(np.nansum(sim_grid.profile("pad_ppl", min_beams=20)) * 0.5) lai = float(np.nansum(sim_grid.profile("lad_ppl", min_beams=20)) * 0.5) print(sim) print(f"PAI {pai:.2f} LAI {lai:.2f} true LAI of the trees {p.trees['leaf_area'].sum() / 20**2:.2f}")

Shots(n_shots=18,000,000, n_echoes=9,919,519)
PAI 4.00  LAI 2.46  true LAI of the trees 2.29

The leaf area index comes out within about a tenth of the truth: with five positions and the misses traced, the voxels see most of the crowns. The PAI is higher, as it should be: it also counts the wood, the shrubs and the grass, which the truth here leaves out.

Leaf angles¶

With inclination=True, normals of the echoes give an inclination angle distribution per tree, and G is integrated over it and over the tree's beam zeniths instead of assuming a spherical distribution. The synthetic leaves have known orientations: the broadleaf trees' are planophile (mostly horizontal blades, a mean inclination of 27 degrees), the eucalypts' erectophile (hanging leaves, 63 degrees).

In [7]:

import pandas as pd

inc = voxels.ray_voxelize(sim, 0.5, box, ground_class=2, leaf_classes=[4], wood_classes=[5], inclination=True)
leaf_tree = np.asarray(p.leaves["tree"])
leaf_inc = np.degrees(np.arccos(np.abs(np.asarray(p.leaves["normal"])[:, 2])))
leaf_w = np.asarray(p.leaves["area"])
rows = []
for tid, t in inc.tree_iad.items():
    if tid == 0:
        continue
    k = leaf_tree == tid
    centres = np.degrees(t["bin_centres"])
    rows.append({"tree": tid, "archetype": p.trees["archetype"][tid - 1],
                 "mean, true (deg)": np.average(leaf_inc[k], weights=leaf_w[k]),
                 "mean, estimated (deg)": np.sum(centres * t["liad"]) / np.sum(t["liad"]),
                 "estimated type": t["liad_de_wit"], "G_leaf": t["g_leaf"]})
rows = pd.DataFrame(rows)
summary = rows.groupby("archetype").agg(trees=("tree", "size"), true_mean=("mean, true (deg)", "mean"),
                                        estimated_mean=("mean, estimated (deg)", "mean"),
                                        estimated_range=("mean, estimated (deg)", lambda v: f"{v.min():.0f} to {v.max():.0f}"),
                                        G_leaf=("G_leaf", "mean"))
summary.round(2)

import pandas as pd inc = voxels.ray_voxelize(sim, 0.5, box, ground_class=2, leaf_classes=[4], wood_classes=[5], inclination=True) leaf_tree = np.asarray(p.leaves["tree"]) leaf_inc = np.degrees(np.arccos(np.abs(np.asarray(p.leaves["normal"])[:, 2]))) leaf_w = np.asarray(p.leaves["area"]) rows = [] for tid, t in inc.tree_iad.items(): if tid == 0: continue k = leaf_tree == tid centres = np.degrees(t["bin_centres"]) rows.append({"tree": tid, "archetype": p.trees["archetype"][tid - 1], "mean, true (deg)": np.average(leaf_inc[k], weights=leaf_w[k]), "mean, estimated (deg)": np.sum(centres * t["liad"]) / np.sum(t["liad"]), "estimated type": t["liad_de_wit"], "G_leaf": t["g_leaf"]}) rows = pd.DataFrame(rows) summary = rows.groupby("archetype").agg(trees=("tree", "size"), true_mean=("mean, true (deg)", "mean"), estimated_mean=("mean, estimated (deg)", "mean"), estimated_range=("mean, estimated (deg)", lambda v: f"{v.min():.0f} to {v.max():.0f}"), G_leaf=("G_leaf", "mean")) summary.round(2)

Out[7]:

	trees	true_mean	estimated_mean	estimated_range	G_leaf
archetype
broadleaf	12	26.79	36.55	34 to 40	0.52
eucalypt	8	63.17	56.28	52 to 63	0.50

In [8]:

fig, ax = plt.subplots(1, 2, figsize=(10, 3.6), sharey=True)
for a, kind in zip(ax, ("broadleaf", "eucalypt")):
    tid = int(rows.tree[rows.archetype == kind].iloc[0])
    t = inc.tree_iad[tid]
    edges = np.linspace(0, 90, len(t["bin_centres"]) + 1)
    k = leaf_tree == tid
    true_hist = np.histogram(leaf_inc[k], edges, weights=leaf_w[k])[0]
    a.stairs(true_hist / true_hist.sum(), edges, color="k", lw=1.5, label="true leaves")
    a.stairs(t["liad"] / np.sum(t["liad"]), edges, color=LEAF, lw=2, label="estimated, leaf echoes")
    a.stairs(t["wiad"] / np.sum(t["wiad"]), edges, color=WOOD, lw=1.5, ls="--", label="estimated, wood echoes")
    a.set(xlabel="inclination of the surface normal from vertical (deg)", xlim=(0, 90), title=f"Tree {tid}, {kind}")
ax[0].set_ylabel("fraction of the surface")
ax[1].legend(loc="upper left");

fig, ax = plt.subplots(1, 2, figsize=(10, 3.6), sharey=True) for a, kind in zip(ax, ("broadleaf", "eucalypt")): tid = int(rows.tree[rows.archetype == kind].iloc[0]) t = inc.tree_iad[tid] edges = np.linspace(0, 90, len(t["bin_centres"]) + 1) k = leaf_tree == tid true_hist = np.histogram(leaf_inc[k], edges, weights=leaf_w[k])[0] a.stairs(true_hist / true_hist.sum(), edges, color="k", lw=1.5, label="true leaves") a.stairs(t["liad"] / np.sum(t["liad"]), edges, color=LEAF, lw=2, label="estimated, leaf echoes") a.stairs(t["wiad"] / np.sum(t["wiad"]), edges, color=WOOD, lw=1.5, ls="--", label="estimated, wood echoes") a.set(xlabel="inclination of the surface normal from vertical (deg)", xlim=(0, 90), title=f"Tree {tid}, {kind}") ax[0].set_ylabel("fraction of the surface") ax[1].legend(loc="upper left");

The estimates tell the two kinds of tree apart, but both are pulled towards the middle: the broadleaf leaves come out steeper than they are, the eucalypt leaves flatter. The normals are fitted to a few noisy echoes on each leaf, and their scatter flattens any distribution towards a uniform one. The wood is mostly stem and limbs, near-vertical surfaces whose normals lie close to 90 degrees from vertical, hence its peak there.

Wood volume, files and streaming¶

QSM cylinders can be rasterised into the same grid; write produces a .vox voxel-space file or a text table; and a shots file is voxelised without being loaded. The cylinders here are the largest tree's own, from the generator.

In [9]:

tid = int(p.trees["tree_id"][np.argmax(p.trees["dbh"])])
model = p.qsm(tid)
sim_grid.add_wood_volume(model)
reach = np.r_[model.start[:, :2], model.end[:, :2]]
print(f"wood volume in the grid {sim_grid.wood_volume.sum():.3f} m3 of the {model.total_volume:.3f} m3 in tree {tid}'s "
      f"cylinders; the rest lies beyond the grid, as the crown reaches from x = {reach[:, 0].min():.1f} to "
      f"{reach[:, 0].max():.1f} m and y = {reach[:, 1].min():.1f} to {reach[:, 1].max():.1f} m")

n = grid.write("tile.vox")
print(n, "voxels written to tile.vox")
streamed = voxels.ray_voxelize(DATA / "litch_tile_shots.parquet", 0.5, ((0, 0, 0), (20, 20, 18)),
                               ground_class=2, leaf_classes=[4])
print("streamed from the file:", streamed,
      "- same echo count:", int(streamed.num_hits.sum()) == int(grid.num_hits.sum()))

tid = int(p.trees["tree_id"][np.argmax(p.trees["dbh"])]) model = p.qsm(tid) sim_grid.add_wood_volume(model) reach = np.r_[model.start[:, :2], model.end[:, :2]] print(f"wood volume in the grid {sim_grid.wood_volume.sum():.3f} m3 of the {model.total_volume:.3f} m3 in tree {tid}'s " f"cylinders; the rest lies beyond the grid, as the crown reaches from x = {reach[:, 0].min():.1f} to " f"{reach[:, 0].max():.1f} m and y = {reach[:, 1].min():.1f} to {reach[:, 1].max():.1f} m") n = grid.write("tile.vox") print(n, "voxels written to tile.vox") streamed = voxels.ray_voxelize(DATA / "litch_tile_shots.parquet", 0.5, ((0, 0, 0), (20, 20, 18)), ground_class=2, leaf_classes=[4]) print("streamed from the file:", streamed, "- same echo count:", int(streamed.num_hits.sum()) == int(grid.num_hits.sum()))

wood volume in the grid 0.896 m3 of the 1.047 m3 in tree 1's cylinders; the rest lies beyond the grid, as the crown reaches from x = 6.2 to 13.7 m and y = 16.8 to 24.9 m
57600 voxels written to tile.vox

streamed from the file: RayVoxelGrid(40x40x36 @ 0.5 m) - same echo count: True