Detecting Cloud Shadows in Drone Imagery

TL;DR: Project each cloud pixel along the solar vector to its ground shadow using solar azimuth, zenith, and an estimated cloud base height, keep only dark low-NIR pixels that fall inside that projected zone, clean the result with scikit-image morphology, and dilate the union of cloud and shadow before it touches any NDVI raster.

Why Shadow Detection Is Harder Than Cloud Detection

Clouds are bright, spectrally flat, and trivial to threshold. Their shadows are the opposite problem: they are dark, and darkness is ambiguous. A cloud shadow depresses reflectance across all bands, but red and NIR drop by unequal fractions, which pushes NDVI down by 0.1–0.3 in the shaded footprint. On a crop-health map that dip is indistinguishable from nitrogen deficiency, waterlogging, or drought — so an unmasked shadow becomes a phantom stress zone that can trigger a wrong variable-rate prescription over several hectares.

A naive fix is to threshold on brightness alone, but on a real field that also catches open water, wet furrows, standing irrigation, and the deep shade inside a dense canopy. The discriminating signal is geometry: a cloud shadow must lie on the anti-solar side of the cloud that casts it, offset by a distance set by the sun’s angle and the cloud’s height. Pairing a loose spectral test (dark and low-NIR) with that hard geometric constraint is what separates true shadows from every other dark thing in the scene. This page covers the drone and small-scene case; the cloud masking for agricultural imagery guide and the Sentinel-2 time-series cloud-removal guide handle the multi-date satellite variant with the same core idea.

Solar projection geometry for cloud shadows The sun at a given zenith angle casts a cloud at height H onto the ground; the shadow is offset along the solar azimuth by a distance equal to H times the tangent of the zenith angle. Ground plane Sun Cloud H (cloud base height) cloud nadir Shadow offset = H · tan(zenith) θ zenith Direction of offset = solar azimuth + 180°

Prerequisites

Beyond the cloud masking cluster requirements, this task only adds scikit-image for the morphology and shift operations:

TEXT 
rasterio==1.3.10
numpy==1.26.4
scikit-image==0.24.0

Install with:

BASH 
pip install rasterio==1.3.10 numpy==1.26.4 scikit-image==0.24.0

Input requirements:

  • A multi-band orthomosaic with at least a red and a NIR band, projected in a metric CRS (e.g. EPSG:32615, UTM zone 15 N) so the shadow offset can be expressed in metres.
  • A boolean cloud mask raster aligned pixel-for-pixel with the orthomosaic (same transform, shape, and CRS). Produce it upstream with the brightness/blue-band test from the parent guide.
  • Solar azimuth and zenith at capture time. Drone flight logs and Sentinel-2 metadata both record these; otherwise compute them from timestamp and centroid with pvlib or astropy.

Step-by-Step

Shadow detection is geometry first, brightness second Four steps: compute the solar azimuth and elevation for the capture time, project a shadow zone from each masked cloud, find dark candidate pixels inside that zone, and verify them by checking they are also dark in near infrared, which stressed vegetation is not. Solar geometry azimuth and elevation Project from each cloud mask Candidates dark, in the projected zone Verify dark in NIR as well Brightness alone cannot separate shadow from stress; the geometry says where a shadow could be, and NIR says whether it is one.

Step 1 — Read the mask, red, and NIR bands

Open the orthomosaic and the cloud mask with rasterio and confirm they share a grid. A silent half-pixel misalignment between the mask and the imagery shifts every projected shadow, so assert the transforms match before doing any geometry.

Step 2 — Project clouds to their shadow footprint

The shadow of a cloud pixel lands on the ground offset from the cloud’s nadir by H * tan(zenith) metres, in the direction azimuth + 180° (the anti-solar direction). Convert that metre offset into a row/column pixel shift using the raster’s pixel size, then shift the entire cloud mask by that vector with scipy/skimage to obtain the projected shadow zone. Because exact cloud height is unknown, sweep a few candidate base heights and keep the shift that best overlaps the dark pixels from Step 3.

Step 3 — Flag dark, low-NIR candidate pixels

Cloud shadows are dark in the visible bands but, critically, also have depressed NIR — which is what separates them from dark healthy vegetation, where NIR stays high. Threshold on both: low visible brightness AND low NIR reflectance. This deliberately over-selects (it will also grab water and wet soil); the geometric filter in the next step removes those.

Step 4 — Intersect candidates with the projected zone

Take the logical AND of the spectral candidate mask and the projected shadow zone. Only dark low-NIR pixels that actually fall where a cloud could cast a shadow survive. This is the step that rejects ponds, ditches, and canopy shade sitting nowhere near a cloud’s anti-solar footprint.

Step 5 — Morphological cleanup and dilation

Remove speckle with remove_small_objects, close single-pixel gaps with a binary closing, then union the shadow mask with the original cloud mask and dilate the result by a few pixels. Penumbra — the soft partial shadow at a cloud’s edge — is not fully dark but still biases NDVI, so a small dilation buffer catches it. The complete, directly runnable script:

PYTHON 
import numpy as np
import rasterio
from skimage.morphology import (
    binary_closing, binary_dilation, remove_small_objects, disk
)

# ── Capture geometry (from flight log or Sentinel-2 metadata) ──────────────
SOLAR_AZIMUTH_DEG = 148.0   # sun bearing, degrees clockwise from north
SOLAR_ZENITH_DEG  = 38.0    # angle from vertical, degrees
CLOUD_BASE_HEIGHTS = [600, 900, 1200, 1600]  # metres AGL to sweep
RED_BAND, NIR_BAND = 3, 4   # 1-based band indices in the orthomosaic

# ── 1. Read orthomosaic + aligned cloud mask ───────────────────────────────
with rasterio.open("ortho.tif") as src:
    red = src.read(RED_BAND).astype("float32")
    nir = src.read(NIR_BAND).astype("float32")
    transform = src.transform
    ortho_shape = (src.height, src.width)
    px_x, px_y = abs(transform.a), abs(transform.e)  # pixel size in CRS units (m)
    assert src.crs.is_projected, "Ortho must be in a metric CRS (e.g. EPSG:32615)"

with rasterio.open("cloud_mask.tif") as msk:
    cloud = msk.read(1).astype(bool)
    assert msk.transform == transform, "Cloud mask grid must match the orthomosaic"
    assert (msk.height, msk.width) == ortho_shape, "Mask/ortho shape mismatch"

# Normalise reflectance to 0–1 if stored as uint DN
if red.max() > 1.5:
    red /= red.max()
    nir /= nir.max()

# ── 2. Spectral candidate test: dark AND low-NIR ───────────────────────────
brightness = (red + nir) / 2.0
DARK_T, NIR_T = 0.12, 0.18            # tune per sensor / illumination
candidates = (brightness < DARK_T) & (nir < NIR_T) & (~cloud)

# ── 3. Project the cloud mask along the anti-solar vector ───────────────────
az_rad = np.deg2rad((SOLAR_AZIMUTH_DEG + 180.0) % 360.0)  # shadow falls opposite the sun
zen_rad = np.deg2rad(SOLAR_ZENITH_DEG)

def project_cloud(cloud_mask, height_m):
    """Shift the cloud mask to where its shadow lands for a given base height."""
    offset_m = height_m * np.tan(zen_rad)
    dx_m = offset_m * np.sin(az_rad)     # easting component
    dy_m = offset_m * np.cos(az_rad)     # northing component
    col_shift = int(round(dx_m / px_x))
    row_shift = int(round(-dy_m / px_y))  # north is up → negative row direction
    return np.roll(np.roll(cloud_mask, row_shift, axis=0), col_shift, axis=1)

# ── 4. Sweep candidate heights, keep the best-overlapping projection ────────
best_zone, best_overlap = np.zeros_like(cloud), -1
for h in CLOUD_BASE_HEIGHTS:
    zone = project_cloud(cloud, h)
    overlap = int(np.logical_and(zone, candidates).sum())
    if overlap > best_overlap:
        best_zone, best_overlap = zone, overlap

# ── 5. Intersect spectral candidates with the projected shadow zone ─────────
shadow = candidates & best_zone

# ── 6. Morphological cleanup + dilate the combined cloud+shadow mask ────────
shadow = remove_small_objects(shadow, min_size=25)
shadow = binary_closing(shadow, disk(2))
combined = binary_dilation(cloud | shadow, disk(3))  # buffer for penumbra

# ── 7. Inline verification on masked fraction ───────────────────────────────
masked_frac = combined.mean()
print(f"Cloud px:  {cloud.sum():>10,}")
print(f"Shadow px: {shadow.sum():>10,}")
print(f"Combined masked fraction: {masked_frac:.1%}")
assert 0.0 <= masked_frac < 0.95, "Almost the whole scene masked — thresholds too loose"
assert shadow.sum() > 0, "No shadow pixels found — check solar geometry and band indices"

# ── 8. Persist the combined mask alongside the orthomosaic ──────────────────
with rasterio.open("ortho.tif") as src:
    profile = src.profile
profile.update(count=1, dtype="uint8", nodata=0)
with rasterio.open("cloud_shadow_mask.tif", "w", **profile) as dst:
    dst.write(combined.astype("uint8"), 1)
    dst.update_tags(1, description="1 = cloud or shadow (dilated); apply before NDVI")

print("Wrote cloud_shadow_mask.tif")

Inline verification: reopen the mask and confirm the shadow footprint sits down-sun of the cloud and covers a plausible fraction of the scene:

PYTHON 
with rasterio.open("cloud_shadow_mask.tif") as src:
    m = src.read(1).astype(bool)
    assert src.count == 1 and src.crs.is_projected
    print(f"Final masked fraction: {m.mean():.1%}")
    # A whole-field shadow > 60% almost never happens on a real drone flight
    assert m.mean() < 0.60, "Implausible shadow coverage — revisit thresholds"

Gotchas & Edge Cases

  • Azimuth sign and north-up rows. The shadow falls in the azimuth + 180° direction, and raster rows increase southward — so the northing offset maps to a negative row shift. Getting either sign wrong projects the shadow onto the sunlit side of the cloud and the height sweep silently returns near-zero overlap. Verify the direction on one obvious cloud before trusting the batch.
Four dark things, told apart by one band A table of four causes of dark patches in drone imagery — cloud shadow, crop stress, standing water and freshly worked soil — showing that all four are dark in the visible bands and only near infrared separates them. Dark patch Visible bands Near infrared Cloud shadow illumination reduced Dark Also dark — everything is dimmed Stressed crop less chlorophyll Dark Still bright — canopy still scatters Standing water specular or absorbing Dark Very dark — water absorbs NIR Freshly worked soil moisture Dark Moderate — no canopy to scatter

  • Dark healthy canopy is not a shadow. Dense, well-watered corn is dark in the visible bands but keeps NIR reflectance high (0.4–0.6). A brightness-only threshold flags it as shadow and erases exactly the vigorous crop you care about. The nir < NIR_T term is non-negotiable — it is the whole reason the spectral test works.

  • Water and irrigation mimic shadows perfectly. Ponds and wet furrows are dark with low NIR and will pass the spectral test. Only the geometric intersection rejects them; never ship the spectral candidate mask alone. If a field has extensive standing water, add a separate water mask so it is not dilated into the shadow buffer.

  • Under-dilation leaks penumbra into NDVI. The soft partial shadow at a cloud edge depresses NDVI without ever going fully dark, so a tight mask still leaves a ring of false stress. A 3–5 pixel dilation is cheap insurance; feed the buffered mask into your threshold mapping for crop health step rather than the raw one.

Frequently Asked Questions

Why do cloud shadows get mistaken for crop stress in NDVI?

A cloud shadow lowers reflectance in every band, but the red and NIR bands drop by different amounts, so the NDVI ratio shifts downward. That drop looks identical to nitrogen deficiency or drought stress in the map, so shadowed pixels must be masked out before any zone is drawn or a prescription is written.

How do I estimate cloud height to project the shadow?

For low-altitude drone flights the cloud base is usually 500 to 2000 m above ground, so a single base-height estimate is enough. Sweep a small range of candidate heights and keep the projection distance that maximises overlap between the projected zone and the detected dark pixels, which avoids needing an exact height.

Can I detect shadows without a cloud mask?

You can run the dark-pixel low-NIR spectral test alone, but it will also flag water, wet soil, and deep canopy shade. The geometric projection from known cloud locations is what removes those false positives, so a cloud mask makes the result far more reliable.

Parent Guide

This guide is part of Cloud Masking for Agricultural Imagery — see there for the full pipeline context, including the upstream cloud-detection thresholds and how the combined mask is applied before index calculation.

Parent section Cloud Masking for Agricultural Imagery Related Automating Cloud Removal in Sentinel-2 Time Series