Handling Edge Effects in Raster Index Generation

TL;DR: Pad input raster bands with numpy.pad(mode='reflect') before any spatial filter, run the filter on the padded array, crop back to the original extent, then reapply the validity mask — this eliminates the artificial boundary gradients that corrupt threshold classifications at field edges.

Why Edge Artifacts Appear in Agricultural Rasters

Vegetation indices such as NDVI, NDRE, and SAVI are mathematically stable at the single-pixel level, but decision-grade precision agriculture pipelines rarely stop there. Teams routinely apply spatial smoothing, rolling-window statistics, or convolution-based texture filters to reduce per-pixel sensor noise and align spectral values with agronomic management zones before threshold mapping for crop health is applied.

When those neighborhood operations reach the raster boundary, the algorithm lacks neighboring pixel data. Default behaviors vary by library: scipy.ndimage defaults to reflect; numpy.convolve zero-fills; plain slicing simply truncates. All three choices produce artificial spectral gradients along field edges, flight-line seams, and drainage features that have nothing to do with actual crop stress.

The agronomic consequences are real. A reflectance gradient created by implicit zero-padding can push NDVI below a critical classification cutoff, triggering false stress alerts, unnecessary scouting events, or misallocated variable-rate inputs. The problem compounds on orthomosaics processed from multiple drone flights — MicaSense RedEdge-MX and DJI P4 Multispectral captures each — where seam lines already carry minor radiometric offsets before any filter is applied. Explicitly managing boundary conditions is therefore a prerequisite for any drone imagery processing & vegetation index workflow that produces actionable output.

The diagram below shows what happens at the raster boundary under each padding strategy and why reflect is the correct default for crop canopy data.

Prerequisites

This page assumes you have the parent cluster’s environment in place — rasterio>=1.3, numpy>=1.24, and scipy>=1.11. The only additional requirement specific to this technique is that your input raster carries a valid nodata value or mask band (produced by any standard photogrammetry pipeline such as Pix4Dmapper, Agisoft Metashape, or DroneDeploy). A projected CRS — for example EPSG:32754 (WGS 84 / UTM zone 54S) — is required because spatial filter parameters such as window_size map to physical distances at the sensor’s ground sampling distance (GSD).

Step-by-Step Implementation

The following numbered steps produce a spatially smoothed NDVI raster that is free of boundary artifacts. The entire implementation lives in one self-contained function; paste it directly into your pipeline.

Step 1 — Open the raster and read bands plus the validity mask.

rasterio.read_masks() returns 255 for valid pixels and 0 for nodata. Reading this mask before any computation gives you the authoritative record of which pixels are real.

Step 2 — Pad input arrays using numpy.pad with mode='reflect'.

The pad_width is window_size // 2, which ensures the filter kernel always has a full neighborhood for every real pixel — including those on the first and last row or column.

Step 3 — Apply the spatial filter on the padded arrays.

scipy.ndimage.uniform_filter accepts a mode argument that mirrors the padding strategy. Pass mode='reflect' here too so the filter’s own internal boundary handling stays consistent with the explicit padding.

Step 4 — Crop back to the original extent.

Slice [pad_width:-pad_width, pad_width:-pad_width] from both smoothed arrays. This discards the artificially generated boundary rows and columns and restores the output to the exact spatial footprint of the input raster.

Step 5 — Compute the vegetation index with a safe division guard.

numpy.where(denominator == 0, 0.0, ...) prevents ZeroDivisionError and NaN propagation on bare soil or shadowed pixels where NIR + Red is near zero.

Step 6 — Reapply the original validity mask and write output.

The mask is applied after computation — not before — so the filter always operates on complete neighborhoods. Setting ndvi[~valid_mask] = np.nan ensures padded pixels never appear in the final layer.

import rasterio
import numpy as np
from scipy.ndimage import uniform_filter

def compute_ndvi_edge_safe(input_path: str, output_path: str, window_size: int = 5) -> None:
    """
    Compute spatially smoothed NDVI with explicit boundary padding.

    Band layout assumed: Band 3 = Red, Band 4 = NIR
    (standard MicaSense RedEdge-MX and DJI P4 Multispectral ordering).
    window_size must be an odd integer for symmetric padding.
    """
    if window_size % 2 == 0:
        raise ValueError(f"window_size must be odd for symmetric padding; got {window_size}.")

    pad_width = window_size // 2

    with rasterio.open(input_path) as src:
        # Step 1: read bands and validity mask
        red = src.read(3).astype(np.float32)
        nir = src.read(4).astype(np.float32)
        valid_mask = src.read_masks(1) == 255  # 255 = valid, 0 = nodata

        assert red.shape == nir.shape == valid_mask.shape, (
            f"Shape mismatch: red={red.shape}, nir={nir.shape}, mask={valid_mask.shape}"
        )

        # Step 2: pad both bands symmetrically using reflect mode
        red_pad = np.pad(red, pad_width, mode="reflect")
        nir_pad = np.pad(nir, pad_width, mode="reflect")

        # Step 3: spatial smoothing on padded arrays
        red_smooth = uniform_filter(red_pad, size=window_size, mode="reflect")
        nir_smooth = uniform_filter(nir_pad, size=window_size, mode="reflect")

        # Step 4: crop back to original extent
        red_c = red_smooth[pad_width:-pad_width, pad_width:-pad_width]
        nir_c = nir_smooth[pad_width:-pad_width, pad_width:-pad_width]

        assert red_c.shape == red.shape, (
            f"Crop failed: expected {red.shape}, got {red_c.shape}"
        )

        # Step 5: safe NDVI computation
        denom = nir_c + red_c
        ndvi = np.where(denom == 0.0, 0.0, (nir_c - red_c) / denom)

        # Step 6: reapply original validity mask
        ndvi[~valid_mask] = np.nan

        # Spot-check: interior mean should be higher than border mean for healthy canopy
        interior = ndvi[pad_width:-pad_width, pad_width:-pad_width]
        border_strip = ndvi[:pad_width, :]
        print(f"Interior NDVI mean : {np.nanmean(interior):.4f}")
        print(f"Border strip mean  : {np.nanmean(border_strip):.4f}")

        # Write output
        profile = src.profile.copy()
        profile.update(dtype=rasterio.float32, count=1, nodata=np.nan)
        with rasterio.open(output_path, "w", **profile) as dst:
            dst.write(ndvi.astype(np.float32), 1)

# Inline verification — run immediately after calling the function:
# with rasterio.open("ndvi_edge_safe.tif") as check:
#     arr = check.read(1)
#     assert not np.any(np.isnan(arr[50:-50, 50:-50])), "Unexpected NaNs in interior"
#     print(f"Output shape: {arr.shape}, dtype: {arr.dtype}, nodata fraction: {np.isnan(arr).mean():.3%}")

Selecting the Right Boundary Mode

For convolutional kernels beyond a uniform filter, the boundary mode must match the kernel’s mathematical assumptions. Gaussian and Laplacian kernels behave correctly with reflect; kernels with asymmetric weights may require nearest to avoid spectral inversion at boundaries.

Scaling to Regional Orthomosaics

Full-field processing works well with the array-based approach above for orthomosaics up to roughly 500 MB. For regional mosaics — multi-field MicaSense captures processed in Pix4Dmapper or Agisoft Metashape — use windowed I/O with an overlap margin:

1. Calculate tile windows with a stride of tile_size - 2 * pad_width so adjacent tiles overlap by exactly pad_width on each side.
2. Read bands and the mask for each window including the overlap margin.
3. Apply padding, filtering, index computation, and cropping as above.
4. Write only the interior (non-overlapping) region to the corresponding output window.
5. Flush the output dataset’s write buffer every 10–20 tiles to bound RAM use.

This pattern keeps peak memory consumption proportional to the tile size rather than the full mosaic, which is essential when running multi-flight campaigns with combined orthomosaics exceeding 10 GB.

Gotchas and Edge Cases

• Band index offset: rasterio read() uses 1-based band indices, but the returned NumPy array is 0-indexed. src.read(3) returns the third raster band as array[0] — not array[2]. Confusing these is the most common cause of silent band swaps that produce negative NDVI values over healthy vegetation.
• Odd window size is mandatory: An even window_size produces an asymmetric pad_width, leaving the cropped output one row or column shorter than the input. The assert red_c.shape == red.shape guard in the code above catches this before the write step.
• Mask must be read before filtering, applied after: Reading valid_mask before any computation preserves the original nodata boundary. Applying it after the filter (not before) is what allows the filter to see complete neighborhoods at field edges. Reversing this order re-introduces the artifact you are trying to eliminate.
• Implicit nodata in drone outputs: Some photogrammetry pipelines output orthomosaics with no explicit nodata value, instead using a per-band alpha channel or a zero-filled border. If src.nodata is None, read the alpha band (src.read(src.count)) or construct a binary mask from (red > 0) | (nir > 0) before proceeding.
• Radiometric seam lines in multi-flight mosaics: reflect padding cannot correct a genuine radiometric offset between adjacent flight strips. Normalize each strip to a shared calibration panel reflectance value before merging, or use a histogram-matching step in your cloud masking for agricultural imagery preprocessing stage.

Validation

After running compute_ndvi_edge_safe, confirm that boundary artifact elimination succeeded with two quick statistical checks:

import rasterio
import numpy as np

with rasterio.open("ndvi_edge_safe.tif") as src:
    arr = src.read(1)

pad = 5  # should match window_size // 2 used during computation
interior = arr[pad:-pad, pad:-pad]
border   = np.concatenate([arr[:pad, :].ravel(), arr[-pad:, :].ravel(),
                            arr[:, :pad].ravel(), arr[:, -pad:].ravel()])

int_mean, int_std = np.nanmean(interior), np.nanstd(interior)
bdr_mean, bdr_std = np.nanmean(border),   np.nanstd(border)

print(f"Interior  mean={int_mean:.4f}  std={int_std:.4f}")
print(f"Border    mean={bdr_mean:.4f}  std={bdr_std:.4f}")

# In a correctly padded raster the means should agree within ~2%
assert abs(int_mean - bdr_mean) / (int_mean + 1e-9) < 0.02, (
    f"Border mean deviates from interior by more than 2% — check padding mode or mask application order."
)

A step-change in mean or a sudden variance drop at the perimeter indicates padding misalignment or incorrect mask application order. Clusters of NaN at the perimeter that extend deeper than pad_width rows usually mean nodata pixels are being propagated from the padded region, which confirms the mask was applied before rather than after the filter.

This guide is part of Threshold Mapping for Crop Health — see there for the full pipeline context, including threshold determination strategies and variable-rate export integration.

Related

• Threshold Mapping for Crop Health — the parent pipeline that consumes edge-corrected index rasters to generate discrete management zones
• Calculating NDVI and NDRE with Rasterio Step by Step — band extraction and safe division patterns that feed directly into the edge-handling workflow here
• Python Script for SAVI Calculation on Drone Tiles — the same pad-filter-crop pattern applied to soil-adjusted indices on tiled drone datasets
• Cloud Masking for Agricultural Imagery — upstream masking step that must run before spatial filtering to prevent cloud shadow pixels from polluting neighborhood windows