Applying a Terrain Illumination Correction

Compute how directly each pixel faces the sun, then normalise reflectance against it:

import numpy as np

cos_i = (np.cos(sun_zen) * np.cos(slope)
         + np.sin(sun_zen) * np.sin(slope) * np.cos(sun_az - aspect))   # local illumination

# SCS+C: bounded correction that accounts for diffuse skylight
corrected = refl * (np.cos(slope) * np.cos(sun_zen) + c) / (cos_i + c)

The same formula as a hillshade, but with the real sun of the acquisition rather than a conventional light. This page belongs to terrain analysis and DEM-derived products in Satellite Processing Workflows & Index Pipelines.

The Problem, Quantified

Reflectance tracks illumination, not forest Pixels of a single forest type plotted by near-infrared reflectance against local illumination fall along a rising line: well-lit slopes are bright and shaded slopes dark. A correction should flatten that line so reflectance no longer depends on illumination. The slope of the line and its intercept give the C parameter used by the correction. One forest type, before correction local illumination cos(i) NIR reflectance slope m, intercept b C = b / m

For a single cover type, reflectance should not depend on how the slope faces the sun, yet uncorrected pixels line up along a clear trend. That trend is both the problem and the tool: its intercept and slope give the empirical C parameter that the correction uses to account for diffuse light.

Environment & Setup

Package Version pin Used for
numpy >=1.23 The illumination geometry and correction
rioxarray >=0.15 Aligning DEM derivatives with the imagery
scikit-learn >=1.3 Robust regression for the C parameter 
pip install "numpy>=1.23" "rioxarray>=0.15" "scikit-learn>=1.3"

Complete Working Example

import numpy as np
from sklearn.linear_model import HuberRegressor


def illumination(slope_deg: np.ndarray, aspect_deg: np.ndarray,
                 sun_zenith_deg: float, sun_azimuth_deg: float) -> np.ndarray:
    s, a = np.radians(slope_deg), np.radians(aspect_deg)
    z, phi = np.radians(sun_zenith_deg), np.radians(sun_azimuth_deg)
    return (np.cos(z) * np.cos(s) + np.sin(z) * np.sin(s) * np.cos(phi - a)).astype("float32")


def fit_c(band: np.ndarray, cos_i: np.ndarray, cover_mask: np.ndarray,
          min_pixels: int = 2_000) -> float:
    """C parameter from reflectance ~ cos_i over one cover type."""
    ok = cover_mask & np.isfinite(band) & np.isfinite(cos_i) & (cos_i > 0)
    if ok.sum() < min_pixels:
        raise ValueError("too few pixels of the reference cover to fit C")
    m = HuberRegressor().fit(cos_i[ok, None], band[ok])
    slope, intercept = float(m.coef_[0]), float(m.intercept_)
    return intercept / slope if slope > 0 else np.inf


def scs_c(band: np.ndarray, slope_deg: np.ndarray, cos_i: np.ndarray,
          sun_zenith_deg: float, c: float) -> np.ndarray:
    s = np.radians(slope_deg)
    z = np.radians(sun_zenith_deg)
    out = band * (np.cos(s) * np.cos(z) + c) / (cos_i + c)
    out[cos_i <= 0] = np.nan                       # self-shadowed: no direct light at all
    return out.astype("float32")

Fitting C over a single cover type — mature forest is the usual choice in mountains — is essential. The regression assumes the only reason reflectance varies is illumination; fitting over mixed cover folds real land-cover differences into C and produces a correction that is subtly wrong everywhere. The sun angles come from the scene’s catalog properties, found with the tools in querying STAC catalogs programmatically.

Three Corrections Compared

How each correction treats shaded slopes The plain cosine correction over-brightens shaded slopes dramatically because it assumes no diffuse light. The C-correction adds an empirical diffuse term and keeps shaded slopes plausible. SCS+C additionally preserves the sun-canopy-sensor geometry of forest stands, which makes it the usual choice for vegetated mountains. Shaded-slope reflectance after correction true uncorrected cosine C SCS+C Cosine divides by small illumination values and blows up; the C term keeps the correction bounded.

The cosine correction is theoretically clean and practically dangerous: it divides by the illumination term, which approaches zero on slopes facing away from the sun, so shaded slopes are brightened far beyond their true reflectance. Adding the C term to both numerator and denominator bounds the correction, and the SCS variant additionally models how forest canopies are lit, which is why SCS+C has become the default for vegetated terrain.

Self-shadowed pixels — those facing so far away from the sun that cos(i) is negative — receive no direct light at all, and no illumination correction can recover them. Masking them, as the example does, is the honest treatment; they are not observable in that acquisition.

When to Correct, and When Not To

Illumination correction matters most for single-band analysis, absolute reflectance, and classification in steep terrain, where illumination can exceed the spectral difference between classes. It matters least for normalised ratios such as NDVI, because illumination multiplies both bands roughly equally and largely cancels in the ratio — though diffuse light and adjacency break that cancellation on strongly shaded slopes.

It also matters for change detection between dates with different sun elevation: a lower winter sun lengthens shadows and changes the illumination on every slope, which a differencing product will report as change. Correcting both dates, or comparing only anniversary dates, removes that false signal, as discussed in separating real change from illumination differences.

Verification

The trend should disappear Before correction, reflectance of a single cover type rises steeply with illumination. After a successful correction the same regression is nearly flat, meaning reflectance no longer depends on how the slope faces the sun. A remaining strong slope indicates C was fitted on mixed cover or the aspect convention is wrong. Regression slope on cos(i), one cover type before: strong dependence after: nearly flat 
import numpy as np

def trend(band, cos_i, mask):
    ok = mask & np.isfinite(band) & np.isfinite(cos_i)
    return float(np.polyfit(cos_i[ok], band[ok], 1)[0])

before = trend(nir, cos_i, forest)
after = trend(nir_corrected, cos_i, forest)
print(f"slope before {before:.3f}, after {after:.3f}")
assert abs(after) < 0.25 * abs(before), "illumination dependence remains"

The regression check is the definitive test because it measures exactly the effect the correction is meant to remove. A substantial remaining slope means C was fitted over mixed cover, the aspect convention is inverted, or the solar geometry came from the wrong scene.

Common Errors

Shaded slopes become extremely bright

The plain cosine correction was used. Switch to C or SCS+C.

The correction makes things worse

Aspect is measured in the wrong convention, so illumination is computed for mirrored slopes. Verify aspect on a known slope.

C comes out negative or infinite

The regression slope is non-positive, usually because C was fitted over too little terrain or mixed cover. Fit over a larger area of one cover type.

Correction is inconsistent between scenes

Each scene needs its own solar angles and its own C. Never reuse either across acquisitions.

Frequently Asked Questions

Q: Why not use the simple cosine correction? Because it assumes all light is direct. Shaded slopes also receive diffuse skylight, so dividing by a small illumination value massively over-brightens them. The C-correction adds an empirical term for diffuse light that keeps the correction bounded.

Q: Where do the solar angles come from? From the scene metadata: the mean solar zenith and azimuth at acquisition time, usually present in the catalog item as view and sun angle properties. For large scenes, per-pixel angle grids exist in some products and are slightly more accurate.

Q: Should I correct before computing indices? For indices that are ratios of bands, much of the illumination effect already cancels, so correction matters less. For absolute reflectance, single bands and classification in steep terrain, correct before anything else.

Q: Does the DEM resolution matter for the correction? Yes. The DEM must resolve the slopes that the imagery resolves; a 90 metre DEM under 10 metre imagery smooths away small ridges, leaving their illumination uncorrected. Use a DEM at or near the imagery resolution, derived at native resolution and resampled onto the image grid.

Q: Can one C parameter serve a whole region? Only if the region shares one dominant cover and one acquisition. C depends on the band, the sun geometry and the surface, so fit it per scene and per band, and refit wherever the reference cover changes substantially across a large area.