Skip to content
FastFuels Docs
Searching in-progress documentation for v2 (Beta). Switch to the stable version for production results.

From lidar to a tree inventory

You are viewing in-progress documentation for v2 (Beta). Switch to the stable version for the current production release.

Somebody flew an aeroplane over the Blackfoot River valley in Montana in 2021, pointed a laser at the ground, and published the returns. In this tutorial we will turn those returns into a tree inventory: one row per tree, each with a position, a height, a diameter, a crown ratio, and a species.

We will work on about half a square kilometre of river valley, building it in six calls and then inspecting what we made:

StepWhat we do
1Create a domain — the patch of ground everything else hangs off
2Run a coverage check — is there lidar here, and whose?
3Fetch a point cloud — 10.5 million returns, clipped to the domain
4Build a canopy height model — the returns rasterized to a 1 m height surface
5Detect trees — a detected inventory, one tree per treetop
6Fill in the rest — a fully attributed inventory with diameter, crown ratio, and species
7Inspect the result, and correct or calibrate it

Every step after the first is asynchronous: one POST starts it, and you GET the resource until its status reaches completed. So the rhythm of the whole tutorial is create → poll → use the id in the next call. In the Python SDK the poll is resource.wait(), and the Python tabs form one continuous script — you import and authenticate once in Step 1, then each step keeps the object it just made and passes it to the next. The whole run takes about three minutes.

Prerequisites

Section titled “Prerequisites”

  1. An API key. Create one in the FastFuels web app under your account settings. Set it once here and it propagates to every code block on the page: my-api-key. The SDK reads it via ff.set_api_key(...) or the FASTFUELS_API_KEY environment variable.

  2. curl, or Python with the FastFuels SDK — pip install fastfuels-sdk (v0.21.0+ for v2).

That’s all — we create everything else as we go.

Step 1 — Create a domain

Section titled “Step 1 — Create a domain”

A domain is the georeferenced extent every other resource belongs to. We POST a GeoJSON FeatureCollection to /domains.

We give the coordinates in UTM zone 12N (EPSG:32612) rather than longitude and latitude, because we already know exactly which patch of valley we want and the lidar we are about to fetch is published for it.

POST /domains
curl -X 'POST' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key' \
  -H 'Content-Type: application/json' \
  -d '{
  "type": "FeatureCollection",
  "name": "Blackfoot River",
  "description": "About half a square kilometre of the Blackfoot River valley, Montana.",
  "crs": { "type": "name", "properties": { "name": "EPSG:32612" } },
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [[
          [294095.0, 5198982.0],
          [294785.0, 5198982.0],
          [294785.0, 5199750.0],
          [294095.0, 5199750.0],
          [294095.0, 5198982.0]
        ]]
      }
    }
  ]
}'

This one returns immediately — there is no job to wait for. Record the id, since every call below hangs off it: your-domain-id.

Notice the bbox in the response: 294095, 5198982 to 294785, 5199750. That is 690 m east–west by 768 m north–south, and those numbers reappear in every resource we build from here.

Domain
The domain outline over satellite imagery — 0.53 km² of the Blackfoot River valley, about 30 km east of Missoula, Montana.

Step 2 — Ask what lidar is there

Section titled “Step 2 — Ask what lidar is there”

Before fetching anything, we ask what the USGS 3D Elevation Program actually holds over this domain. The answer comes back immediately — there is no job to wait for — and asking first is a habit worth keeping: it tells us in advance whether the fetch will succeed, and what it will contain.

GET pointclouds/3dep/coverage
curl -X 'GET' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/pointclouds/3dep/coverage' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key'

One acquisition, MT_Statewide_P3_4_B21, covers the whole domain on its own — contribution_fraction is 1.0. That is the simplest case, and it is why this valley makes a good first lesson.

Note the acquisition’s name. We are going to pin it in the next call.

Step 3 — Fetch the point cloud

Section titled “Step 3 — Fetch the point cloud”

Now we pull the returns. Passing datasets pins the fetch to the acquisition we just found, so this tutorial produces the same cloud today and next year even if USGS publishes new lidar over the valley.

POST pointclouds/3dep
curl -X 'POST' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/pointclouds/3dep' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key' \
  -H 'Content-Type: application/json' \
  -d '{
  "name": "Blackfoot ALS (pinned)",
  "description": "Pinned to a single 3DEP acquisition so the fetch is reproducible.",
  "tags": ["blackfoot", "3dep"],
  "datasets": ["MT_Statewide_P3_4_B21"]
}'

Record the id — your-point-cloud-id — and poll until status is completed. This one takes about twenty seconds.

GET the point cloud
curl -X 'GET' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/pointclouds/your-point-cloud-id' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key'

10,512,255 returns at 19.8 per square metre. Two things in that response are worth stopping on.

First, the coverage check in Step 2 estimated 7,847,291 points and we received 10.5 million — about a third more. The estimate comes from the acquisition’s published average density, and real flight lines overlap, so the estimate is a floor rather than a forecast. See estimate versus delivered for the detail.

Second, point_classes lists 1, 2, 7, 9, 18, 20. Class 2 is ground and class 1 is unclassified — mostly vegetation here. Those two are what the next step needs: without ground returns, there is no surface to measure height above.

The returns themselves, coloured by height above ground:

An oblique three-dimensional view of a 120 metre block of the lidar point cloud, with an inset map showing where that block sits inside the domain. Brown points form a continuous ground surface; individual conifer crowns rise out of it as distinct spikes, shading from pale green at their base to dark green at their tips.

A 120 m block of the domain, holding 358,942 returns — the inset shows where the block sits. Brown points are at ground level, green points are canopy. Thinned for display. The full cloud covers the whole domain and holds 10,512,255 returns.

Step 4 — Build the canopy height model

Section titled “Step 4 — Build the canopy height model”

A canopy height model is a raster where each cell holds the height of the tallest return above the ground beneath it. This is the step that turns a cloud of points into a surface we can look for trees in.

POST grids/canopy/point_cloud
curl -X 'POST' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/grids/canopy/point_cloud' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key' \
  -H 'Content-Type: application/json' \
  -d '{
  "source_point_cloud_id": "your-point-cloud-id",
  "name": "Canopy height from 3DEP lidar",
  "description": "CHM rasterized from the pinned Blackfoot point cloud.",
  "tags": ["blackfoot", "chm"]
}'

Record the grid id — your-chm-grid-id — and poll it to completed, about thirty seconds.

GET the grid
curl -X 'GET' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/grids/your-chm-grid-id' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key'

We asked for no resolution, so the grid defaulted to 1 m: a shape of [768, 690], one cell per metre of the domain.

Stop and read source.ground

Section titled “Stop and read source.ground”

This is the most important field on the page, and it is the reason to build a CHM from a cloud rather than upload one you rasterized elsewhere:

"ground": {
  "max_ground_distance_m": 35.2,
  "ground_source": "classification",
  "ground_coverage": 0.9053
}

A canopy height is a subtraction: the top of the canopy minus the ground under it. If the ground is wrong, every height above it is wrong by the same amount, and nothing downstream can tell. So the grid reports how the ground was established:

  • ground_source: "classification" — the ground came from returns the vendor labelled as ground (class 2), which is the trustworthy case. Had the cloud carried no usable ground returns, this would read derived and the ground would be an estimate.
  • ground_coverage: 0.9053 — just over 90% of cells had ground returns near enough to use directly. The rest were interpolated.

Look also at the chm band summary: count is 493,763 and nodata_count is 36,157 — about 7% of the domain has no height at all. That is the river. Water absorbs the pulse and returns nothing, so those cells are genuinely empty rather than zero. The tallest cell is 37.68 m.

For what each of these fields means in full, and what to do when ground_coverage is low, see Read source.ground.

The canopy height model of the Blackfoot River domain, rendered green on cream. Individual tree crowns appear as distinct round dots, densest in a band along both banks. A pale grey ribbon winds across the lower half where the river returns no data.

The finished CHM — 10.5 million returns rasterized to 1 m. Every dot is a crown. The grey ribbon is the Blackfoot River: water absorbs the pulse, so those cells are the 7% nodata counted above.

Step 5 — Detect the trees

Section titled “Step 5 — Detect the trees”

Now we look for treetops. A treetop is a local maximum in the height surface: a cell taller than everything around it. The lmf algorithm slides a window over the CHM and places one tree at each peak it finds.

Two parameters shape the result. footprint_size: 5 is the window, in pixels — on our 1 m grid, a 5 m square. min_height: 2 ignores anything under 2 m, so we detect trees rather than shrubs.

POST inventories/tree/chm
curl -X 'POST' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/inventories/tree/chm' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key' \
  -H 'Content-Type: application/json' \
  -d '{
  "name": "Detected treetops",
  "source_chm_grid_id": "your-chm-grid-id",
  "algorithm": {
    "name": "lmf",
    "min_height": 2,
    "footprint_size": 5
  }
}'

Record this inventory’s id — inventory-to-complete — and poll it to completed, under a minute.

GET the inventory
curl -X 'GET' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/inventories/inventory-to-complete' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key'

7,317 trees. We have gone from a cloud of points to a list of individual stems, and we can already see the shape of the stand: heights run from the 2 m floor we set to 37.68 m, averaging 12.4 m.

Notice that the tallest tree is exactly as tall as the tallest CHM cell — 37.684879 m in both. It has to be: a detected treetop is a cell of the height surface, so detection can never invent a height the CHM did not already hold.

Now notice what is missing. The columns are x, y, and height — and forestry_metrics is null. A height surface seen from above cannot tell us how thick a trunk is, how far down the crown reaches, or what species we are looking at, so the API does not guess. Without a diameter there is no basal area, and without basal area there are no stand metrics.

That is the honest limit of the measurement, and Step 6 is how we get past it.

A 200 metre square of the canopy height model at high magnification, with a small white circle drawn on the peak of each detected crown. The circles sit one per crown across dense forest, thinning out toward the bare river edge along the top.

Detected treetops over a 200 m window of dense forest on the south bank — 693 of the 7,317 trees. Each circle is one row of the inventory, placed at a local maximum of the height surface.

Step 6 — Fill in what the CHM could not see

Section titled “Step 6 — Fill in what the CHM could not see”

GDAM — a Generalized Dendro Allometric Model — takes each tree’s position and height and imputes the attributes lidar could not observe. We point it at the inventory we just made.

POST inventories/tree/allometry/gdam
curl -X 'POST' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/inventories/tree/allometry/gdam' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key' \
  -H 'Content-Type: application/json' \
  -d '{
  "name": "Fully attributed tree inventory",
  "source_tree_inventory_id": "inventory-to-complete"
}'

This creates a new inventory rather than modifying the detected one — record its id as your-inventory-id and poll it to completed, about a minute.

GET the inventory
curl -X 'GET' \
  'https://api-v2-prod-nyvjyh5ywa-uw.a.run.app/domains/your-domain-id/inventories/your-inventory-id' \
  -H 'accept: application/json' \
  -H 'api-key: my-api-key'

Same 7,317 trees, in the same places, at the same heights — with three new columns: dbh (diameter at breast height, in cm), crown_ratio (the fraction of the tree’s height that is crown), and fia_species_code.

And forestry_metrics is no longer null. It now carries the numbers a forester would ask for — note that these follow FIA convention rather than the metric columns above: basal area in ft²/acre, tree density in trees per acre, and quadratic mean diameter in inches. Our stand comes out at roughly 56 trees per acre, a basal area of 20.5 ft²/acre, and a quadratic mean diameter of 8.2 in, with Douglas-fir holding about 66% of the basal area and true fir/hemlock most of the rest.

Nothing about the stand changed between Step 5 and Step 6. What changed is that the diameters now exist, and every metric that depends on them followed.

Three histograms side by side across all 7,317 trees: canopy height in metres with a long tail out to 37 m, diameter at breast height in centimetres with a hard floor near 11 cm, and crown ratio spanning only 0.42 to 0.56.

The three columns across all 7,317 trees. Height is measured, and its long tail is the real spread of the stand. Diameter and crown ratio are imputed, and it shows: diameter has a floor near 11 cm, and crown ratio spans barely a tenth. A model asked for a value it cannot observe returns a typical one.

Step 7 — Inspect the trees

Section titled “Step 7 — Inspect the trees”

Every column we just added came out of a model, and a model returns a number whether or not that number is right. Before building anything on this inventory, look at it.

The map below draws all 7,317 trees. Click any one to see its stored values, and zoom in where the canopy is dense:

Domain Trees
Every tree in the inventory. Click a marker for its species, height, diameter, and crown ratio.

Fetch the same values for the whole inventory with the data endpoints, as Parquet, CSV, or GeoJSON.

What to look for

Section titled “What to look for”
  • Detections on things that are not trees. Detection runs on a surface model, which records what is physically there rather than what is vegetation. This domain has a clear example — see below.
  • Counts that do not match the stand. Our 7,317 detections over 0.53 km² are 56 trees per acre. If your knowledge of the site says otherwise, footprint_size and min_height from Step 5 are the first things to revisit; Tune CHM detection walks through exporting the CHM and the treetops together and reading the overlay.
  • Attribute ranges. Across the inventory the imputed columns span 11.0–53.0 cm dbh and 0.42–0.56 crown_ratio. A range far narrower than a real stand’s is expected of an imputed value, but check the centre of it against local data before trusting stand totals.
  • Species composition. See Species below.

The power line

Section titled “The power line”

A transmission line crosses the north-east of this domain, and detection has placed a chain of trees along it:

A 315 metre square of the canopy height model with detected treetops drawn as circles. Two thin parallel lines run diagonally across it, and the detections lying on them are circled in red, forming a continuous chain, while the grey circles on the surrounding round crowns are real trees.

The conductors appear in the CHM as two thin parallel lines, and 191 of the 7,317 detections — 2.6% of the inventory — sit on them.

A conductor holds its elevation while the ground falls away beneath it, so it stands well above the ground surface and reads as canopy. In this cloud the vendor left the wires in ASPRS class 1, so they survive into the CHM.

Today this is yours to correct. The detection step has no notion of infrastructure, and nothing in the pipeline will flag it for you — which is the reason to inspect. Remove them with a spatial modification rule: draw a polygon over the corridor and use it as the condition, exactly as Remove trees with road and water features does with a feature. Both the inline-geometry and feature-reference forms are described in About modifications.

Automatic detection of power lines and other structures is something we are looking at for a future release. If it matters for your work, tell us at support.fastfuels@silvxlabs.com — it helps us prioritise.

Species

Section titled “Species”

GDAM assigns species from position and height, and this is where its limits are easiest to see. The inventory comes back:

SpeciesTrees
Douglas-fir4,658
western hemlock2,626
red alder19
Pacific silver fir10
mountain hemlock4

Western hemlock is unlikely on this site, and ponderosa pine — the species most likely to be here — does not appear at all. GDAM is an imputation and it has inaccuracies; treat its species assignment as a starting point to check against what you know of the site, not as an observation.

fia_species_code is a modifiable attribute, so a reassignment is a modification rule like any other. Conditions on it accept eq and ne only, since it is categorical. To reassign every western hemlock (FIA code 263) to ponderosa pine (122):

{
  "modifications": [
    {
      "conditions": [
        { "attribute": "fia_species_code", "operator": "eq", "value": 263 }
      ],
      "actions": [
        { "attribute": "fia_species_code", "modifier": "replace", "value": 122 }
      ]
    }
  ]
}

POST that to .../inventories/{inventory_id}/modifications and poll the inventory back to completed. See Modify an inventory for the full request and the polling loop.

Correcting and calibrating

Section titled “Correcting and calibrating”

Nothing here requires rebuilding from the point cloud. Rules are appended to the inventory and applied to its stored trees under the same id, so the order you apply them in is the order they take effect.

To do thisUseWhere it is documented
Change the tree countRe-run detection with different parametersTune CHM detection
Drop trees by size or attributeAn attribute or expression condition with a remove actionRemove small trees, Remove by expression
Drop trees in a placeA spatial condition — inline polygon or a featureRemove trees under roads or water, Remove trees with features
Reassign speciesfia_species_code with a replace actionSpecies above
Calibrate a numeric columnmultiply, divide, add, subtract, or replaceScale an attribute
Keep the original to compareDuplicate before modifyingBranch a scenario

For the concepts behind all of it — conditions, actions, buffer_m, and how create-time rules differ from in-place ones — read About modifications.

What we built, and what to question

Section titled “What we built, and what to question”

We started from published laser returns and ended with 7,317 trees we can voxelize, treat, or export. It is worth being clear about which parts we measured and which parts we modelled:

ColumnWhere it came from
x, yMeasured — the position of a peak in the height surface
heightMeasured — a return, above a ground surface built from other returns
dbh, crown_ratio, fia_species_codeImputed by GDAM from position and height

Both halves are legitimate; they just carry different kinds of uncertainty, and three things are worth carrying forward:

  • A CHM sees the overstory only. Trees beneath the dominant canopy return no visible peak, so 7,317 is a count of what can be seen from above, not of what is standing (why).
  • A surface model sees whatever is physically there. The power line in Step 7 is 2.6% of this inventory, on a domain that is otherwise almost all vegetation. A site with buildings would carry more.
  • The ground was 90.5% classified, not 100%. The remaining cells rest on an interpolated surface, and their heights inherit that.

The whole pipeline in one script

Section titled “The whole pipeline in one script”

Every call above, in order, with the polling folded in:

build_inventory.py — 3DEP lidar to an attributed tree inventory
"""Blackfoot River: 3DEP lidar to a fully attributed tree inventory.


Runs the whole tutorial end to end with the FastFuels v2 SDK. Takes about
three minutes. `resource.wait()` replaces the manual poll loop — it blocks
until the resource reaches `completed`, or raises if it `failed`.
"""


import fastfuels_sdk.v2 as ff
from fastfuels_sdk.v2.client_library.models import StemIsolationLmf


ff.set_api_key("my-api-key")


# 1 — the domain (Blackfoot River valley, EPSG:32612)
domain = ff.Domain.from_geojson(
    {
        "type": "FeatureCollection",
        "crs": {"type": "name", "properties": {"name": "EPSG:32612"}},
        "features": [
            {
                "type": "Feature",
                "properties": {},
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [
                        [
                            [294095.0, 5198982.0],
                            [294785.0, 5198982.0],
                            [294785.0, 5199750.0],
                            [294095.0, 5199750.0],
                            [294095.0, 5198982.0],
                        ]
                    ],
                },
            }
        ],
    },
    name="Blackfoot River",
    description="About half a square kilometre of the Blackfoot River valley, Montana.",
)
print(f"domain {domain.id}")


# 2 — is there lidar here, and whose?
coverage = ff.point_clouds.check_3dep_coverage(domain).to_dict()
print(f"coverage {coverage['coverage_fraction']:.0%}, "
      f"~{coverage['estimated_point_count']:,} points estimated")


# 3 — fetch it, pinned to one acquisition so the run is reproducible
point_cloud = ff.point_clouds.create_point_cloud_from_3dep(
    domain,
    datasets=["MT_Statewide_P3_4_B21"],
    name="Blackfoot ALS (pinned)",
    description="Pinned to a single 3DEP acquisition so the fetch is reproducible.",
    tags=["blackfoot", "3dep"],
)
point_cloud.wait()
summary = point_cloud.to_dict()["summary"]
print(f"point cloud {point_cloud.id}: "
      f"{summary['point_count']:,} points, {summary['density']:.1f} pts/m2")


# 4 — rasterize the canopy surface
grid = ff.grids.create_canopy_height_grid_from_point_cloud(
    point_cloud,
    name="Canopy height from 3DEP lidar",
    description="CHM rasterized from the pinned Blackfoot point cloud.",
    tags=["blackfoot", "chm"],
)
grid.wait()
ground = grid.to_dict()["source"]["ground"]
print(f"chm {grid.id}: ground from {ground['ground_source']}, "
      f"{ground['ground_coverage']:.1%} covered")


# 5 — one tree per treetop
detected = ff.inventories.create_tree_inventory_from_chm_grid(
    domain,
    grid,
    algorithm=StemIsolationLmf(min_height=2, footprint_size=5),
    name="Detected treetops",
)
detected.wait()
print(f"detected {detected.id}: "
      f"{[c['key'] for c in detected.to_dict()['columns']]}")


# 6 — fill in what the CHM could not see
inventory = ff.inventories.create_tree_inventory_from_gdam(
    domain,
    detected,
    name="Fully attributed tree inventory",
)
inventory.wait()
doc = inventory.to_dict()
metrics = doc["forestry_metrics"]
print(f"inventory {inventory.id}: {[c['key'] for c in doc['columns']]}")
print(f"  {metrics['tree_count']:,} trees, "
      f"{metrics['tree_density']:.0f} per acre, "
      f"QMD {metrics['quadratic_mean_diameter']:.1f} in")
for group in metrics["dominant_species_groups"]:
    print(f"  {group['basal_area_share']:.1%} {group['name']}")

Where to go next

Section titled “Where to go next”