Machine Learning for Precision Forestry: A Case Study in Satellite-Based Forest Monitoring

Forestry is one of those industries where data science can create enormous value, but the gap between a promising prototype and a production system is wider than most practitioners expect. This post walks through a real project — anonymized but technically accurate — where we built an ML-based forest monitoring system for a large plantation operator in Southeast Asia.

The goal was deceptively simple: detect individual tree crowns from satellite imagery, classify their health status, and track changes over time. The reality involved coordinate system headaches, class imbalance that would make a fraud detection engineer weep, and the constant tension between model accuracy and operational utility.

The Problem Context

Our client managed approximately 200,000 hectares of plantation forest across multiple concessions. Their existing monitoring process was manual: field teams would sample plots on a rotating basis, visiting each hectare roughly once every 3-5 years. This meant that disease outbreaks, illegal logging, and fire damage could go undetected for months.

The ask was straightforward: “Can you use satellite imagery to tell us which trees are healthy, which are stressed, and which are dead — across the entire estate — every two weeks?”

The answer was yes, but the path from question to production was anything but straight.

Satellite Data Pipeline

Data Sources

We used a combination of optical and radar satellite data:

• Sentinel-2 (optical): 10m resolution, 5-day revisit, free. Provided spectral information for vegetation health assessment.
• Sentinel-1 (SAR): 10m resolution, 6-day revisit, free. Radar data that penetrates clouds — critical in tropical environments where cloud cover exceeds 70% on most days.
• Planet Labs (optical): 3m resolution, daily revisit, commercial. Used for validation and high-resolution analysis in priority areas.

The cloud cover problem cannot be overstated. In tropical forestry, you might get 2-3 usable optical images per month. SAR data becomes essential for temporal continuity.

Preprocessing Pipeline

The preprocessing pipeline was the most engineering-intensive component:

import rioxarray
import xarray as xr
import numpy as np
from pathlib import Path

class Sentinel2Preprocessor:
    """Production Sentinel-2 preprocessing for forestry applications."""

    def __init__(self, output_resolution: float = 10.0):
        self.output_resolution = output_resolution
        self.bands = {
            "B02": "blue",
            "B03": "green",
            "B04": "red",
            "B05": "rededge1",
            "B06": "rededge2",
            "B07": "rededge3",
            "B08": "nir",
            "B8A": "nir_narrow",
            "B11": "swir1",
            "B12": "swir2",
        }

    def process_safe(self, safe_path: Path) -> xr.Dataset:
        """Process a Sentinel-2 SAFE directory to analysis-ready data."""
        # Read and stack bands
        band_datasets = {}
        for band_id, band_name in self.bands.items():
            band_path = self._find_band(safe_path, band_id)
            ds = rioxarray.open_rasterio(band_path)
            ds = ds.rio.reproject(ds.rio.crs, resolution=self.output_resolution)
            band_datasets[band_name] = ds

        stacked = xr.concat(list(band_datasets.values()), dim="band")
        stacked = stacked.assign_coords(band=list(band_datasets.keys()))

        # Apply SCL-based cloud mask
        scl_path = self._find_band(safe_path, "SCL")
        scl = rioxarray.open_rasterio(scl_path)
        cloud_mask = self._create_cloud_mask(scl)

        # Apply mask
        stacked = stacked.where(~cloud_mask)

        # Calculate vegetation indices
        indices = self._calculate_indices(stacked)

        return xr.merge([stacked, indices])

    def _create_cloud_mask(self, scl: xr.DataArray) -> xr.DataArray:
        """Create cloud mask from Scene Classification Layer."""
        # SCL classes: 3=cloud shadow, 8=cloud medium, 9=cloud high, 10=cirrus
        cloud_classes = [3, 8, 9, 10]
        mask = scl.isin(cloud_classes)
        return mask.squeeze()

    def _calculate_indices(self, ds: xr.Dataset) -> xr.Dataset:
        """Calculate forestry-relevant spectral indices."""
        nir = ds.sel(band="nir").astype(float)
        red = ds.sel(band="red").astype(float)
        rededge1 = ds.sel(band="rededge1").astype(float)
        swir1 = ds.sel(band="swir1").astype(float)

        indices = {}

        # NDVI — general vegetation health
        indices["ndvi"] = (nir - red) / (nir + red + 1e-10)

        # NDRE — red edge NDVI, more sensitive to chlorophyll content
        indices["ndre"] = (nir - rededge1) / (nir + rededge1 + 1e-10)

        # NDMI — moisture content, critical for drought stress detection
        indices["ndmi"] = (nir - swir1) / (nir + swir1 + 1e-10)

        # LAI proxy — simplified
        indices["lai_proxy"] = np.log((0.7 - indices["ndvi"]) / (0.1 - indices["ndvi"] + 1e-10)) / 2.5

        return indices

The Scene Classification Layer (SCL) is Sentinel-2’s built-in classification of each pixel. It’s imperfect — I’ve seen it misclassify haze as clear sky — but it’s the best automated cloud mask available without running your own cloud detection model.

Temporal Compositing

Because of cloud cover, we couldn’t rely on single images. Instead, we created temporal composites — combining the best pixels from a 14-day window:

def create_temporal_composite(
    images: list[xr.Dataset],
    composite_period_days: int = 14,
    method: str = "max_ndvi",
) -> xr.Dataset:
    """Create a cloud-free composite from multiple images."""
    # Stack all images along time dimension
    stacked = xr.concat(images, dim="time")

    if method == "max_ndvi":
        # For each pixel, select the observation with highest NDVI
        # This naturally selects cloud-free, sunlit observations
        ndvi = stacked["ndvi"]
        best_idx = ndvi.argmax(dim="time")

        # Index into original data using the best time index
        composite = stacked.isel(time=best_idx)

    elif method == "median":
        # Median is more robust to outliers
        composite = stacked.median(dim="time")

    elif method == "most_recent_cloud_free":
        # Use the most recent observation that passes cloud threshold
        cloud_threshold = 0.1  # Maximum acceptable cloud probability
        # ... implementation depends on cloud probability band

    return composite

The “max NDVI” compositing method works surprisingly well for forestry applications. Cloud pixels have low NDVI, so selecting the maximum NDVI observation naturally picks cloud-free pixels. It fails in areas of actual low vegetation (recent harvest, bare soil), but for monitoring existing forest, it’s robust.

Tree Crown Detection

The Instance Segmentation Approach

Individual tree crown detection is fundamentally an instance segmentation problem — each tree is a distinct object that needs to be localized and delineated. We approached this with a two-stage pipeline:

Stage 1: Crown detection — locate tree tops using a local maximum filter on the canopy height model (CHM).

Stage 2: Crown delineation — segment individual crowns using a watershed algorithm seeded from detected treetops.

import numpy as np
from scipy.ndimage import maximum_filter, label
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from skimage.morphology import disk

def detect_tree_crowns(
    chm: np.ndarray,
    min_height: float = 5.0,
    min_distance: int = 5,
    smoothing_sigma: float = 1.5,
) -> np.ndarray:
    """Detect and delineate individual tree crowns from CHM."""
    from scipy.ndimage import gaussian_filter

    # Smooth CHM to reduce noise
    chm_smooth = gaussian_filter(chm, sigma=smoothing_sigma)

    # Mask out non-tree areas
    tree_mask = chm_smooth > min_height

    # Find local maxima (treetops)
    coordinates = peak_local_max(
        chm_smooth,
        min_distance=min_distance,
        threshold_abs=min_height,
        labels=tree_mask.astype(int),
    )

    # Create marker image for watershed
    markers = np.zeros_like(chm, dtype=int)
    for i, (row, col) in enumerate(coordinates, 1):
        markers[row, col] = i

    # Watershed segmentation
    # Use inverted CHM so watershed fills from peaks
    crowns = watershed(-chm_smooth, markers, mask=tree_mask)

    return crowns

The min_distance parameter is the single most impactful hyperparameter. For plantation forests with regular spacing (e.g., 3m x 3m oil palm), you can set this to the known spacing. For natural forest, you need to estimate from the data. Too small and you get over-segmentation (one tree becomes many); too large and you get under-segmentation (many trees become one).

Deep Learning Approach

For natural forest where tree spacing is irregular and crowns overlap, classical methods struggle. We experimented with a Mask R-CNN approach:

import torch
import torchvision
from torchvision.models.detection import maskrcnn_resnet50_fpn

def build_tree_detector(num_classes: int = 2):
    """Build Mask R-CNN for tree crown detection."""
    model = maskrcnn_resnet50_fpn(pretrained=True)

    # Replace the classifier head
    in_features = model.roi_heads.box_predictor.cls_score.in_features
    model.roi_heads.box_predictor = torchvision.models.detection.faster_rcnn.FastRCNNPredictor(
        in_features, num_classes
    )

    # Replace the mask predictor
    in_features_mask = model.roi_heads.mask_predictor.conv5_mask.in_channels
    model.roi_heads.mask_predictor = torchvision.models.detection.mask_rcnn.MaskRCNNPredictor(
        in_features_mask, 256, num_classes
    )

    return model

Training data was the bottleneck. We had field inventory data for ~500 plots, each with GPS-tagged tree positions. Converting these to polygon masks for Mask R-CNN required careful alignment between field GPS coordinates (typically ±3m accuracy) and the satellite imagery.

The deep learning approach improved crown detection accuracy from ~65% (watershed) to ~82% (Mask R-CNN) in dense natural forest, but at significantly higher computational cost and complexity. For plantation forest with regular spacing, the watershed approach was sufficient and much simpler to maintain.

Tree Health Classification

Feature Engineering

Health classification relied on spectral signatures combined with temporal patterns:

def extract_tree_features(
    crown_mask: np.ndarray,
    composite: xr.Dataset,
    time_series: list[xr.Dataset],
) -> dict:
    """Extract features for a single tree crown."""
    # Spatial features from current composite
    crown_pixels = composite.where(crown_mask > 0)

    features = {
        # Spectral means
        "ndvi_mean": float(crown_pixels["ndvi"].mean()),
        "ndre_mean": float(crown_pixels["ndre"].mean()),
        "ndmi_mean": float(crown_pixels["ndmi"].mean()),

        # Spectral variability (within-crown heterogeneity)
        "ndvi_std": float(crown_pixels["ndvi"].std()),
        "red_std": float(crown_pixels.sel(band="red").std()),

        # Crown shape features
        "crown_area_pixels": int((crown_mask > 0).sum()),
        "crown_aspect_ratio": calculate_aspect_ratio(crown_mask),

        # Temporal features (from time series)
        "ndvi_trend": calculate_ndvi_trend(time_series, crown_mask),
        "ndvi_seasonal_amplitude": calculate_seasonal_amplitude(time_series, crown_mask),
        "days_since_ndvi_drop": detect_ndvi_drop(time_series, crown_mask),
    }

    return features

The temporal features were the most predictive. A single snapshot can’t distinguish between a deciduous tree in dry season and a dead tree. But the trend over 6 months tells you everything: a healthy tree shows seasonal variation, a stressed tree shows gradual decline, and a dead tree shows a sudden drop followed by consistently low values.

Classification Model

We used a gradient-boosted tree model (LightGBM) for health classification:

import lightgbm as lgb
from sklearn.model_selection import StratifiedKFold

def train_health_classifier(features_df, labels):
    """Train tree health classifier with proper validation."""
    # Classes: healthy, stressed, dead, gap (no tree)
    class_weights = {
        "healthy": 1.0,
        "stressed": 3.0,  # Upweight rare class
        "dead": 5.0,      # Most important to detect
        "gap": 1.0,
    }

    params = {
        "objective": "multiclass",
        "num_class": 4,
        "metric": "multi_logloss",
        "learning_rate": 0.05,
        "num_leaves": 31,
        "min_child_samples": 50,
        "class_weight": class_weights,
        "verbose": -1,
    }

    # Stratified k-fold with spatial blocking
    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    scores = []
    for train_idx, val_idx in cv.split(features_df, labels):
        train_data = lgb.Dataset(
            features_df.iloc[train_idx],
            label=labels.iloc[train_idx],
        )
        val_data = lgb.Dataset(
            features_df.iloc[val_idx],
            label=labels.iloc[val_idx],
        )

        model = lgb.train(
            params,
            train_data,
            num_boost_round=1000,
            valid_sets=[val_data],
            callbacks=[lgb.early_stopping(50)],
        )

        scores.append(evaluate_model(model, features_df.iloc[val_idx], labels.iloc[val_idx]))

    return model, scores

The Class Imbalance Problem

In a typical plantation, 85-90% of trees are healthy, 5-10% are stressed, and 1-3% are dead. This extreme imbalance means a naive classifier achieves 90% accuracy by just predicting “healthy” for everything.

We addressed this through multiple strategies:

1. Class weights in the loss function — as shown above, upweighting stressed and dead classes
2. Focal loss variant — for neural network approaches, focal loss down-weights easy examples
3. Oversampling with SMOTE — synthetic minority oversampling for training data
4. Threshold tuning — optimizing classification thresholds per class rather than using argmax
5. Evaluation metrics that matter — focusing on recall for stressed/dead classes, not overall accuracy

from sklearn.metrics import classification_report, confusion_matrix

def evaluate_for_business_impact(y_true, y_pred):
    """Evaluate with focus on operational utility."""
    report = classification_report(y_true, y_pred, output_dict=True)

    # Key metric: what percentage of stressed/dead trees do we catch?
    stressed_recall = report["stressed"]["recall"]
    dead_recall = report["dead"]["recall"]

    # Key metric: what percentage of alerts are false alarms?
    stressed_precision = report["stressed"]["precision"]
    dead_precision = report["dead"]["precision"]

    # Business-adjusted score
    # Missing a dead tree costs $X in delayed response
    # A false alarm costs $Y in unnecessary field visit
    dead_miss_cost = 500  # USD per missed dead tree
    false_alarm_cost = 50  # USD per unnecessary field visit

    n_dead = (y_true == "dead").sum()
    n_false_dead = ((y_pred == "dead") & (y_true != "dead")).sum()
    n_missed_dead = ((y_pred != "dead") & (y_true == "dead")).sum()

    cost = n_missed_dead * dead_miss_cost + n_false_dead * false_alarm_cost

    return {
        "stressed_recall": stressed_recall,
        "dead_recall": dead_recall,
        "total_cost_usd": cost,
        "cost_per_hectare": cost / n_hectares,
    }

Change Detection

The most operationally valuable output wasn’t the health classification itself — it was the change in health status over time. A tree that transitions from healthy to stressed between two observation periods is far more actionable than a tree that’s been stressed for a year.

def detect_health_changes(
    current_classification: np.ndarray,
    previous_classification: np.ndarray,
    min_change_confidence: float = 0.7,
) -> np.ndarray:
    """Detect significant health status changes between periods."""
    # Define transition significance
    transitions = {
        ("healthy", "stressed"): "early_warning",
        ("healthy", "dead"): "rapid_decline",
        ("stressed", "dead"): "progression",
        ("stressed", "healthy"): "recovery",
        ("dead", "healthy"): "regrowth",  # Possible misclassification or replanting
    }

    changes = np.full_like(current_classification, "no_change", dtype=object)

    for (from_class, to_class), change_type in transitions.items():
        mask = (previous_classification == from_class) & (current_classification == to_class)
        changes[mask] = change_type

    return changes

“Rapid decline” transitions (healthy → dead) were the highest-priority alerts. These typically indicate fire, illegal logging, or acute disease — events that require immediate response. We routed these to a separate alert channel with different SLAs than gradual stress detection.

Business Impact

The system we built produced measurable outcomes:

• Detection speed: Forest disturbances detected within 2 weeks instead of 3-5 years (field sampling cycle)
• Coverage: 100% of managed area monitored continuously vs. ~5% annual field sampling
• Early disease detection: Disease outbreaks identified an average of 4 months earlier than field detection, when treatment is still viable
• Cost reduction: Monitoring cost reduced from ~$12/hectare/year (field teams) to ~$1.50/hectare/year (satellite + ML)
• Yield improvement: Early intervention on stressed areas improved yield by an estimated 8-12% over 3 years

The ROI was clear within the first year. A single detected disease outbreak that would have gone unnoticed for 6 months under the old system — affecting approximately 2,000 hectares — would have cost more than the entire ML system’s development cost.

Lessons Learned

1. Cloud cover is the dominant challenge in tropical forestry. Plan for SAR data integration from the start. Optical-only systems have unacceptable gaps.

2. Temporal features beat spectral features for health classification. A single image can’t reliably distinguish stressed from healthy. Time series analysis is essential.

3. Class imbalance requires multi-pronged mitigation. No single technique (class weights, oversampling, threshold tuning) is sufficient alone. Combine them all.

4. Field validation data is the bottleneck. The model is only as good as the training labels. Invest heavily in field data collection protocols and quality control.

5. Operational integration matters more than model accuracy. A 75% accurate model that integrates with the client’s workflow beats a 90% accurate model that requires manual data transfer. Build the end-to-end pipeline first.

6. Tree crown detection in natural forest remains unsolved at scale. Plantation forest with regular spacing is tractable. Natural forest with irregular spacing, overlapping crowns, and multi-layered canopy is still an open research problem for satellite-based approaches.

Precision forestry is one of the most impactful applications of geospatial ML. The combination of free satellite data, mature processing tools, and clear business metrics makes it an ideal domain for data science teams looking to deliver measurable value. But the path from prototype to production requires the same engineering discipline as any production ML system — and a healthy respect for the complexity of natural systems.