Canopy Height Estimation¶
Overview¶
This notebook demonstrates how to estimate canopy height from RGB aerial/satellite imagery using Meta's HighResCanopyHeight model integrated into geoai.
The model uses a DINOv2 backbone with a DPT (Dense Prediction Transformer) decoder to predict per-pixel canopy height in meters from standard RGB imagery. It was trained on NAIP aerial imagery with Aerial LiDAR ground truth.
Reference: Tolan et al., "Very high resolution canopy height maps from RGB imagery using self-supervised vision transformer and convolutional decoder trained on Aerial Lidar," Remote Sensing of Environment, 2024. https://doi.org/10.1016/j.rse.2023.113888
Install packages¶
Uncomment the following line to install the required packages.
# %pip install geoai-py
Import libraries¶
import os
import numpy as np
import matplotlib.pyplot as plt
import geoai
from geoai.canopy import CanopyHeightEstimation, list_canopy_models
List available models¶
Several model variants are available with different trade-offs between accuracy and computational requirements.
models = list_canopy_models()
for name, desc in models.items():
print(f"{name:30s} -> {desc}")
Download sample data¶
We'll use a NEON aerial image of a forested area (WLOU site) with a corresponding LiDAR-derived Canopy Height Model (CHM) for comparison.
image_url = "https://huggingface.co/datasets/giswqs/geospatial/resolve/main/2017_WLOU_1_NEON_D13_WLOU_DP3_420000_4415000_RGB.tif"
chm_url = "https://huggingface.co/datasets/giswqs/geospatial/resolve/main/2017_WLOU_1_NEON_D13_WLOU_DP3_420000_4415000_CHM.tif"
image_path = geoai.download_file(image_url)
chm_path = geoai.download_file(chm_url)
Visualize input image and ground-truth CHM¶
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
import rasterio
# Input image
with rasterio.open(image_path) as src:
img = src.read()[:3]
img_display = np.moveaxis(img, 0, 2)
axes[0].imshow(img_display)
axes[0].set_title("Input Aerial Image")
axes[0].axis("off")
# Ground-truth CHM
with rasterio.open(chm_path) as src:
chm_gt = src.read(1)
chm_gt[chm_gt < 0] = 0
vmax = np.percentile(chm_gt[chm_gt > 0], 98) if np.any(chm_gt > 0) else 1
im = axes[1].imshow(chm_gt, cmap="viridis", vmin=0, vmax=vmax)
axes[1].set_title("Ground-Truth CHM (LiDAR)")
axes[1].axis("off")
plt.colorbar(im, ax=axes[1], fraction=0.046, pad=0.04, label="Height (m)")
plt.tight_layout()
plt.show()
Initialize the model¶
Create a CanopyHeightEstimation instance. The model checkpoint will be downloaded automatically on first use (~749 MB for the default compressed model).
For aerial imagery (e.g., NEON, NAIP), use compressed_SSLhuge_aerial which was fine-tuned on aerial data. For satellite imagery (e.g., Maxar), use SSLhuge_satellite (requires GPU).
estimator = CanopyHeightEstimation(model_name="compressed_SSLhuge_aerial")
Run canopy height prediction¶
The predict() method processes the input image in 256×256 tiles and outputs per-pixel canopy height in meters. Adjacent tiles overlap by 128 pixels (default) and are blended using raised-cosine weights for seamless output.
output_path = "canopy_height_output.tif"
height_map = estimator.predict(
image_path,
output_path=output_path,
batch_size=4,
)
Examine results¶
print(f"Height map shape: {height_map.shape}")
print(f"Height range: {height_map.min():.2f} - {height_map.max():.2f} meters")
print(f"Mean height: {height_map.mean():.2f} meters")
print(f"Non-zero pixels: {(height_map > 0.1).sum()} / {height_map.size}")
Visualize results¶
Show the input image, predicted canopy height, and ground-truth CHM side by side.
fig, axes = plt.subplots(1, 3, figsize=(20, 6))
# Input image
with rasterio.open(image_path) as src:
img = src.read()[:3]
img_display = np.moveaxis(img, 0, 2)
axes[0].imshow(img_display)
axes[0].set_title("Input Aerial Image")
axes[0].axis("off")
# Predicted canopy height
vmax = np.percentile(height_map[height_map > 0], 98) if np.any(height_map > 0) else 1
im1 = axes[1].imshow(height_map, cmap="viridis", vmin=0, vmax=vmax)
axes[1].set_title("Predicted Canopy Height")
axes[1].axis("off")
plt.colorbar(im1, ax=axes[1], fraction=0.046, pad=0.04, label="Height (m)")
# Ground-truth CHM
with rasterio.open(chm_path) as src:
chm_gt = src.read(1)
chm_gt[chm_gt < 0] = 0
im2 = axes[2].imshow(chm_gt, cmap="viridis", vmin=0, vmax=vmax)
axes[2].set_title("Ground-Truth CHM (LiDAR)")
axes[2].axis("off")
plt.colorbar(im2, ax=axes[2], fraction=0.046, pad=0.04, label="Height (m)")
plt.tight_layout()
plt.show()
Visualize the saved GeoTIFF¶
The output height map is saved as a GeoTIFF with the same georeferencing as the input.
fig = estimator.visualize(output_path, cmap="viridis")
Using the convenience function¶
For quick one-off predictions, use the canopy_height_estimation() function directly.
from geoai.canopy import canopy_height_estimation
height_map2 = canopy_height_estimation(
image_path,
output_path="canopy_height_output2.tif",
model_name="compressed_SSLhuge_aerial",
)
Clean up¶
# import glob
# for f in ["canopy_height_output.tif", "canopy_height_output2.tif"]:
# if os.path.exists(f):
# os.remove(f)
# for f in glob.glob("*WLOU*"):
# os.remove(f)