Land Cover Classification from Satellite Imagery: Lessons from Oil & Gas

Why Land Cover Matters Upstream

In upstream oil and gas, land cover classification from satellite imagery is not a research problem. It is an operational one. You need to know whether a given patch of land is forest, shrubland, water, or bare soil because that classification determines environmental impact assessments, community boundary demarcations, regulatory compliance, and — in some jurisdictions — whether you are legally allowed to operate there at all.

The consequence of getting it wrong is not a metric regression. It is regulatory exposure or a bad impact assessment that takes months to undo.

I built and maintained a pipeline for this problem at Accenture, working with upstream clients in Southeast Asia. What follows are the lessons that did not fit in the model cards.

The Class Imbalance Problem Is Structural

Satellite imagery of an oil field looks roughly the same everywhere: a large, relatively uniform patch of disturbed land surrounded by whatever the regional biome is. In Indonesia and Malaysia, that is usually tropical forest. In drier regions, it is shrubland or grassland.

The consequence is a dataset where the dominant class (regional vegetation) comprises 80–90% of pixels, infrastructure classes (roads, pads, facilities) comprise 1–5%, and everything else — water bodies, bare soil, urban — sits in between.

Standard data augmentation helps but does not fully compensate. What actually worked:

Stratified patch sampling. Rather than sampling patches randomly across the scene, we sampled by class label: each training batch had a guaranteed minimum representation from every class in the target taxonomy. This required knowing the class distribution before sampling, which meant running a fast, low-accuracy classifier first and using its predictions as the sampling frame.

class StratifiedPatchSampler(Sampler):
    """Guarantees each batch contains at least `min_per_class` samples per class."""
 
    def __init__(self, labels: list[int], num_classes: int,
                 min_per_class: int, batch_size: int):
        self.by_class = defaultdict(list)
        for idx, label in enumerate(labels):
            self.by_class[label].append(idx)
        self.num_classes = num_classes
        self.min_per_class = min_per_class
        self.batch_size = batch_size
 
    def __iter__(self):
        # Guarantee min_per_class slots; fill remaining slots randomly
        while True:
            batch = []
            for cls in range(self.num_classes):
                indices = self.by_class.get(cls, [])
                if indices:
                    batch.extend(random.choices(indices, k=self.min_per_class))
            # Fill the rest of the batch from all classes uniformly
            all_indices = [i for idxs in self.by_class.values() for i in idxs]
            remaining = self.batch_size - len(batch)
            batch.extend(random.choices(all_indices, k=max(0, remaining)))
            random.shuffle(batch)
            yield from batch[:self.batch_size]

Class-weighted loss. Inverse frequency weighting in the loss function was table stakes. We experimented with focal loss variants (from the RetinaNet paper) but found vanilla inverse frequency weighting performed better for our taxonomy once stratified sampling was in place. The interaction between sampling strategy and loss weighting matters more than either alone.

# Compute class weights from training set pixel counts
pixel_counts = torch.tensor([count_by_class[c] for c in range(NUM_CLASSES)], dtype=torch.float)
class_weights = (1.0 / pixel_counts) / (1.0 / pixel_counts).sum()
class_weights = class_weights * NUM_CLASSES  # scale so mean weight ≈ 1
 
criterion = nn.CrossEntropyLoss(weight=class_weights.to(device))

Refusing to treat rare classes as optional. Early in the project there was pressure to collapse the taxonomy — merge infrastructure subtypes into a single "infrastructure" class. That would have improved headline metrics. It would have also destroyed the operational utility of the output. We kept the fine-grained taxonomy and accepted lower per-class accuracy for rare classes in exchange for keeping them visible.

Resolution Trade-Offs Are Not Technical, They Are Economic

Satellite imagery comes in a range of resolutions, from 30m/pixel (Landsat, free) to 0.3m/pixel (commercial providers, expensive). The resolution you choose determines not just model accuracy but the entire cost structure of the system.

For field-level land cover classification, 10m/pixel (Sentinel-2, free) was the practical floor. Anything coarser and roads become invisible. Anything finer and the storage, preprocessing, and inference cost grows quadratically.

We used a two-stage approach: Sentinel-2 at 10m for wall-to-wall coverage and change detection, Pleiades or Maxar imagery at 0.5–0.3m for targeted inspection of flagged areas. The change detection stage was cheap; the high-resolution inspection stage was expensive but applied only to a small fraction of the total area.

The lesson here is not technical: it is that resolution selection is a product decision, not a data science decision. Someone has to decide what precision the operational consumer actually needs and what they are willing to pay for it. Data scientists who let this decision default to "use the best available data" end up with pipelines that are unusable in production.

Domain Shift Is Unavoidable; Plan for It

A model trained on imagery from one site or one season will degrade on imagery from a different site, a different season, or after the sensor is recalibrated. This is not a bug; it is the nature of remote sensing data.

The failure mode is insidious: the model continues to produce outputs without error messages. The outputs look plausible. They are wrong in ways that only a domain expert would notice — forest boundaries shift by a few pixels, a recently cleared area is misclassified as the previous vegetation class, infrastructure under construction appears as bare soil.

What worked:

Temporal compositing. Rather than classifying a single image, we classified a median composite of twelve months of imagery, filtering out clouds. This reduced the effect of seasonal phenology and atmospheric variation. It also made the system useful only for annual snapshots, not real-time monitoring — a trade-off we accepted.

def temporal_composite(scenes: list[np.ndarray], cloud_masks: list[np.ndarray]) -> np.ndarray:
    """Pixel-wise median over valid (non-cloudy) observations across a year of scenes."""
    stacked = np.stack(scenes, axis=0)   # (T, H, W, C)
    masks = np.stack(cloud_masks, axis=0)  # (T, H, W), True = cloudy
    # Mask cloudy pixels as NaN before taking the median
    stacked = stacked.astype(np.float32)
    stacked[masks] = np.nan
    composite = np.nanmedian(stacked, axis=0)  # (H, W, C)
    return composite

Site-specific fine-tuning. The base model was trained on a curated multi-site dataset. For each new operational area, we collected a small set of ground-truth labels (200–500 polygons, collected via field survey or high-resolution visual interpretation) and fine-tuned the model on those labels. The fine-tuned model outperformed the base model by a significant margin on the local site and degraded gracefully on adjacent sites.

Explicit drift monitoring. We tracked the distribution of predicted class probabilities over time. A shift in the distribution — even if the hard predictions looked the same — was a signal to trigger re-evaluation. This caught one case where a sensor calibration update changed the radiometric values enough to shift the probability distribution without visibly changing the top-1 predictions.

def monitor_prediction_drift(
    baseline_probs: np.ndarray,   # (N, num_classes)
    current_probs: np.ndarray,    # (N, num_classes)
    threshold: float = 0.05,
) -> dict[str, float]:
    """Compare max-confidence distributions using KS test. Flag drift if p < threshold."""
    from scipy import stats
    results = {}
    for cls in range(baseline_probs.shape[1]):
        stat, pval = stats.ks_2samp(baseline_probs[:, cls], current_probs[:, cls])
        results[f"class_{cls}"] = pval
        if pval < threshold:
            print(f"[DRIFT] Class {cls}: KS p={pval:.4f} — retrain recommended")
    return results

Deployment: The Last 20% Takes 80% of the Time

The modeling work was straightforward once the data pipeline was stable. The hard part was operationalizing the output.

Clients do not consume GeoTIFF files. They consume shapefiles or GeoJSON polygons in their GIS systems, often with specific coordinate reference systems, attribute schemas, and classification codes mandated by their environmental reporting requirements.

We built a postprocessing layer that converted raster predictions to vector polygons (using GDAL's gdal_polygonize), applied minimum mapping unit filters (patches smaller than a threshold were merged into the dominant adjacent class), converted to the client's required CRS, and populated the attribute schema.

def apply_minimum_mapping_unit(
    class_map: np.ndarray,
    min_pixels: int = 9,           # e.g. 9 pixels × (10m)² = 900 m²
) -> np.ndarray:
    """Merge small regions into the most common adjacent class."""
    from scipy import ndimage
    result = class_map.copy()
    for cls in np.unique(class_map):
        binary = (class_map == cls).astype(int)
        labeled, num_features = ndimage.label(binary)
        for region_id in range(1, num_features + 1):
            region_mask = labeled == region_id
            if region_mask.sum() < min_pixels:
                # Find the most common neighboring class
                dilated = ndimage.binary_dilation(region_mask)
                neighbors = class_map[dilated & ~region_mask]
                if len(neighbors):
                    dominant = np.bincount(neighbors.ravel()).argmax()
                    result[region_mask] = dominant
    return result

This postprocessing layer broke more often than the model did, because the model's output format was stable but the client's requirements changed every quarter.

The other operational challenge was explaining uncertainty to non-technical stakeholders. A probability of 0.63 for "tropical forest" is meaningless to an environmental compliance officer. We developed a three-category confidence scheme (high/medium/low) based on calibrated probability thresholds and included it as an attribute in the output shapefile. Low-confidence areas were flagged for manual review. This was imperfect but it was interpretable.

What I Would Do Differently

The biggest mistake was not establishing a gold-standard validation dataset early enough. We validated on the same high-resolution imagery we used for label generation, which introduced optimistic bias. A held-out set collected independently — ideally from field surveys — would have given us a more honest picture of real-world accuracy.

The second mistake was scoping the taxonomy too late. The classification scheme evolved over the first year of the project as clients discovered what they actually needed. Taxonomy changes invalidate historical predictions, because a class that existed in year one may no longer exist in year two, or may have been split into two classes. Locking down the taxonomy early — even if it means lobbying hard for time to do it — is worth the effort.

Remote sensing AI for operational use is less about the model and more about the data pipeline, the postprocessing layer, and the validation protocol. The model is the easy part.

Companion notebook — runnable Jupyter notebook covering all four lessons (stratified sampling, class-weighted loss, domain shift demo, and postprocessing) with synthetic data: github.com/dzakybd/satellite-land-cover-classification