Very high resolution canopy height maps from RGB imagery using self-supervised vision transformer and convolutional decoder trained on Aerial Lidar

In plain words

Satellite photos show where trees are, but not how tall they are. Airborne lidar, a laser scanner flown over an area, measures height well, but it is scarce and does not cover the whole world. The GEDI laser measures height across much of the globe, but each reading covers a patch about 25 m wide, too coarse to pick out single trees.

This work combines the two. A vision model is first trained without labels on a very large set of satellite image tiles, so it learns general features of vegetation. A second, smaller network is then trained on those features to predict height, with aerial lidar as the answer key. A third step uses GEDI to adjust the output, to correct for the fact that the lidar covers only a few places. The result is a height map at the resolution of the satellite photos, fine enough in some places to show individual trees.

TL;DR

The paper produces canopy height maps for California and São Paulo state from RGB Maxar satellite imagery. A self-supervised DINOv2 ViT-H is pretrained on 18 million unlabeled tiles, then frozen while a DPT decoder learns heights from NEON aerial lidar. A GEDI-trained network rescales the output per pixel. The authors report an average MAE of 2.8 m over three lidar test sets tab. 2. Only the inference half is in the public code; training and the GEDI correction are described in the paper alone. Several numbers in the paper disagree with each other (see Pitfalls), so check against the table cited.

Problem

Canopy height and crown area at high resolution matter for carbon accounting and for monitoring forest restoration and agroforestry §1. Existing global canopy maps are around 10 m, built from Sentinel-2 and GEDI §1. Airborne lidar is the accurate, high-resolution reference, but it is scarce, and models trained on it may not transfer to places without lidar §1. GEDI measures canopy height with about 25 m beams and cannot resolve individual trees §3. The method assumes very high resolution optical imagery (under 1 m ground sample distance) and lidar-derived canopy height models (CHMs) for training sites; the training CHMs come from NEON sites in the USA, about 5,800 of them §2.4. The authors state that ALS coverage is limited to US sites, so they expect the ALS model to work best in ecosystems similar to the training set §3.

Method

The method has three parts fig. 1 fig. 2 fig. 3:

1. Self-supervised encoder §3.1. DINOv2 is run on 18 million randomly sampled 256×256 Maxar tiles, with no labels §2.3.2. The encoder is a huge ViT with 16×16 patches §3.1. Two versions of the network (teacher and student) receive different crops, some masked, and are trained to agree. The SSL procedure itself follows DINOv2, with only the normalisation constants changed app. D.
2. ALS decoder §3.2. The encoder weights are frozen (the teacher's weights). A DPT decoder is trained from scratch on 256×256 random crops of paired RGB and 1 m ALS CHM tiles §2.4 §3.2. The decoder outputs a histogram over 256 height bins per pixel; the prediction is the expected height over the bins §3.2. The loss is Sigloss, borrowed from monocular depth estimation §3.2 eq. 1.
3. GEDI model §3.3. A separate CNN (five convolution layers and five fully connected layers) is trained on 128×128 Maxar patches to predict the GEDI RH95 value, with an L1 loss. Its inputs also include latitude, longitude, off-nadir angle, sun zenith angle and terrain slope app. B.3. Terrain slope is used in training and set to zero at inference app. B.4.
4. Rescaling §3.4. A dense factor γ(x) is computed from the smoothed GEDI prediction and the smoothed 95th percentile of the ALS map over the same region eq. 2 eq. 3. The ALS map is multiplied by γ, and γ is clipped to [0.5, 2] §3.4. The ratio form is chosen to down-weight noisy estimates near zero height §3.4.

Key equations, numbered in the order they appear in the source eq. 1:

\mathcal{L} = \alpha \sqrt{\frac{1}{T}\sum_i \delta_i^2 - \frac{\lambda}{T^2}\Big(\sum_i \delta_i\Big)^2}, \quad \delta_i = \log \hat{c}_i - \log c_i eq. 1

C'(x) = \gamma(x) \odot C(x) eq. 2

\gamma(x) = \frac{1 + s_{\sigma}(G(x))}{1 + s_{\sigma}(Q(x)_{95})} eq. 3

Here T is the number of pixels with valid ground truth, \lambda = 0.85 and \alpha = 10 §3.2. G(x) is the GEDI model output, Q(x)_{95} is the upsampled 95th percentile of the ALS CHM, and s_\sigma is a Gaussian smoothing with σ = 20 pixels §3.4.

The classification head and loss, in pseudocode. Names and the bin range come from the paper and code; see Paper vs code for where they differ.

# Expected height from 256 bin logits  [§3.2] [code: models/dpt_head.py#L521-L531]
# logits: (B, 256, H, W) from the DPT decoder
# bin_values: 256 height values (paper: [0, B]; code: linspace 0.001 to 10)
def expected_height(logits, bin_values):
    # 1. Turn logits into weights over bins; the code applies ReLU, adds
    #    0.1 and divides by the sum  [code: models/dpt_head.py#L527-L530]
    weights = normalise(relu(logits) + 0.1, axis=bins)  # (B, 256, H, W)
    # 2. Height is the weighted sum of the bin values  [§3.2]
    height = sum(weights * bin_values, axis=bins)       # (B, H, W)
    return height

# Sigloss over valid pixels  [§3.2] [eq. 1]
# lam = 0.85 and alpha = 10 are fixed by the paper
def sigloss(pred, truth, lam=0.85, alpha=10.0):
    # 1. Keep only pixels with a valid ground-truth height  [§3.2]
    pred = pred[valid]                                  # (T,)
    truth = truth[valid]                                # (T,)
    # 2. Per-pixel difference of log heights  [§3.2] [eq. 1]
    delta = log(pred) - log(truth)                      # (T,)
    # 3. Mean square minus lam times square of mean  [eq. 1]
    #    (lam / T^2) * (sum delta)^2 equals lam * mean(delta)^2
    spread = mean(delta ** 2) - lam * mean(delta) ** 2
    # 4. Scale and take the square root  [eq. 1]
    return alpha * sqrt(spread)

# GEDI-based rescaling of the ALS map  [§3.4] [eq. 2] [eq. 3]
# als_chm: (H, W) heights in m; gedi_pred and als_p95: lower resolution
def rescale(als_chm, gedi_pred, als_p95):
    # 1. Smooth both lower-resolution maps, sigma = 20 px  [§3.4]
    smooth_gedi = gaussian_blur(gedi_pred, sigma=20)
    smooth_p95 = gaussian_blur(als_p95, sigma=20)
    # 2. Ratio form, to down-weight noise near zero  [§3.4] [eq. 3]
    gamma = (1 + smooth_gedi) / (1 + smooth_p95)
    # 3. Clip the factor to [0.5, 2]  [§3.4]
    gamma = clip(gamma, 0.5, 2.0)
    # 4. Upsample to ALS resolution and multiply  [§3.4] [eq. 2]
    return upsample(gamma, als_chm.shape) * als_chm

Hyperparameters

Results

Author-reported. Setting: ALS-trained models, Maxar satellite input at 0.59 m GSD; set-aside test sets NEON test, CA-Brande and São Paulo §2.4. Headline comparisons from tab. 2:

Metric definitions app. C:

R^2_{\mathrm{block}} = 1 - \frac{\sum_i (c^{(b)}_i - \hat{c}^{(b)}_i)^2}{\sum_i (c^{(b)}_i - \bar{c}^{(b)})^2} eq. 4

\mathrm{ME} = \frac{1}{|\mathcal{D}|} \sum_i (\hat{c}_i - c_i) eq. 5

The Edge Error (EE) is defined in alg. 1 app. C.3: Sobel edge maps of the predicted and truth CHMs are compared, and the summed absolute difference is divided by the summed edge magnitude of both maps. Lower is better.

On the 20,000-sample GEDI test, SSL+GEDI has the highest R² in 13 of 15 subregions and the GEDI rescaling raises R² in 42 of 45 model–subregion pairs §4.1.3 tab. 3.

Cost of the SSL model and the aerial result are described in app. D and §4.5.

Paper vs code

The public repository is facebookresearch/HighResCanopyHeight (Apache-2.0), at the commit in the frontmatter. The paper does not link it. Its README links the paper and lists the same authors, and the owner is Meta's facebookresearch organisation README.md.

• Scope. The release contains inference only. It has no Sigloss, no SSL or decoder training loop and no GEDI correction. The README says the GEDI correction is not included README.md. All training details therefore come from the paper alone.
• Classification bins. The paper takes a scalar product with a "linear vector ranging [0, B]" with B = 256 §3.2. The code builds 256 bins with linspace(0.001, 10, 256) dpt_head.py [code: models/dpt_head.py#L525]. Before normalising it applies ReLU and adds 0.1 [code: models/dpt_head.py#L527-L529]. The released wrapper multiplies the head output by 10 inference.py. Output scale in the code therefore differs from the paper's description.
• Input normalisation. The paper gives no constants app. D. The inference script uses a red-channel mean of 0.420 inference.py. The src/transforms.py file uses 0.430 transforms.py. Which value the released checkpoints were trained with is not stated.
• Evaluation clipping and ReLU. Negative truth heights are set to zero [code: inference.py#L116], and predictions are passed through ReLU [code: inference.py#L205]. The paper does not describe either step.
• Block R². The paper uses 50×50 blocks §4.1 app. C.1. The code crops a 3-pixel border and then averages with a 50-pixel pool on the 253 remaining pixels, giving a 5×5 grid inference.py. The paper does not mention the border crop.
• Aerial normalisation. The paper describes matching the 5th and 95th percentiles to the paired satellite image app. G. The code's default (normtype=2) uses a learned quantile predictor instead [code: inference.py#L176-L177] inference.py. The paired-satellite variant is normtype=1 inference.py. That branch refers to s, which is not defined in the code shown inference.py; it appears to fail as written.
• Released checkpoint vs paper. The README says the INT8 compressed huge model (compressed_SSLhuge.pth) is "used in the evaluations of the paper" README.md. The loader applies dynamic quantisation [code: inference.py#L64-L70]. The paper does not mention quantisation.
• Aerial numbers. The README table for aerial input gives MAE 3.15, R² 0.51 and bias −1.6 for the huge satellite model README.md. Table 6 gives MAE 3.0, R² 0.55 and ME +1.7 for the satellite-trained SSL model [tab. 6]. The bias sign and size differ. The cause is not established.
• Feature grid. Paper: 8×8×1280 per layer for the huge model app. D. With a 256-pixel input and 16-pixel patches the code produces a 16×16 grid [code: models/backbone.py#L773] backbone.py. The decoder configuration matches the paper otherwise: depth 32, 1280 channels, 20 heads, layers 9/16/22/29 and post-process channels 160/320/640/1280 [code: inference.py#L31-L43].

Pitfalls

• Numbering. Equation numbers above are inferred from source order in the LaTeX; the compiled PDF was not checked. Section, table and appendix labels come directly from the source. Table numbers are in source order.
• Abstract ME. The abstract reports a Mean Error of 0.6 m. Table 2 has no ME column. The ME values in tab. 1 and tab. 7 do not give that figure for the SSL+GEDI model. Its source is not identified.
• Upper vs lower accuracy. Section 4.2.1 states user's accuracy 0.88 and producer's accuracy 0.82 for the global set. tab. 5 lists the SSL+GEDI row as 0.82 / 0.88 under the U/P column, which is the reverse. The table caption defines U/P as user's/producer's.
• Thresholds. The set-aside segmentation uses a 5 m canopy threshold §4.2. The human-labelled set uses a 1 m threshold §4.2.1. Results are not comparable across the two.
• Sample counts. §2.3.3 says 9,000 thumbnails were annotated; §4.2.1 uses 8,903. The difference is not explained.
• Learning-rate schedule. §3.2 says one-cycle with warmup for the encoder training. app. D gives a cosine schedule with 12k warmup for the decoder. Optimiser, weight decay and whether the decoder schedule is the same are not stated.
• Log of zero. Sigloss takes logs of heights. The paper says T counts pixels with "valid ground truth" but does not say how zero-height truth or zero predictions are handled. Clipping to a small positive value is not described.
• Square root. The term under the root is non-negative for λ ≤ 1, but it is zero when all δ are equal. The gradient of √x is unbounded at zero. No stabilising constant is stated.
• Bins and scale. See Paper vs code. The bin range, the +0.1 offset before normalising and the ×10 output factor are not in the paper. An implementation that follows only the paper will not match the released model's output scale.
• Rescaling detail. The paper does not say whether the upsampling happens before or after the Gaussian smoothing, or whether the clip applies to γ or to the output.
• GEDI model input size. The ALS model takes 256×256 tiles, and each tile is split into four 128×128 crops for the GEDI model §3.4. The paper does not spell out the mapping of the four crop values back to the tile.
• Terrain slope. Slope is an input during training and zero at inference app. B.4. Without the GEDI code there is no way to check this.
• Saturation. Predictions saturate above about 30 m for GEDI RH95 §5. Field validation covers trees under 25 m only §5.
• Sampling bias. The lidar training sites are in the USA §3. The São Paulo field check uses Espírito Santo because São Paulo NFI data were not available §4.1.4.
• Model cost. The huge SSL model takes about 16,896 GPU-hours and an estimated 1.8 t CO₂ app. D. The code runs compressed checkpoints on CPU inference.py, and the README says they were not tested on GPU README.md.
• Data access. Training data are not in the repository. The NEON list is a CSV URL app. A. The map outputs are on AWS Open Data and a Google Earth Engine viewer §4. This card does not redistribute the Maxar, NEON or GEDI data.
• Text vs table (§4.5). The text compares the satellite-trained SSL model's NEON MAE with 2.7. Table 1 gives 2.6 for that model without GEDI, and 2.7 is the SSL+GEDI value in Table 2. The text does not say which row it means.

Lineage

• Builds on: self-supervised DINOv2 (Oquab et al., 2023) for the encoder §3.1; DPT decoder (Ranftl et al., 2021) §3.2; Sigloss from Eigen et al. (2014) §3.2; GEDI simulator (Hancock et al., 2019) §3.4.
• Compared against: ResUNet-18 and a SWAG ViT-L encoder §3.5; Lang et al. (2022) and Potapov et al. (2021) canopy maps §1 §4; Wagner et al. (2023), a submeter California map trained on NAIP imagery §1 §4.5.
• Superseded by: not stated in this card's sources.