Calculating NDVI and NDRE with Rasterio Step by Step

Calculating NDVI and NDRE with Rasterio step by step requires loading multispectral imagery, extracting the correct spectral bands, applying normalized difference formulas using NumPy vectorization, and writing the results back to georeferenced TIFFs. The core workflow avoids Python loops entirely by leveraging rasterio for I/O and numpy for array math. You will convert raw digital numbers (DN) to floating-point arrays, mask invalid pixels, compute (NIR - Target) / (NIR + Target), and preserve the original spatial metadata during export.

1. Environment Setup & Dependency Management

Install the required stack. Rasterio 1.3.0+ and NumPy 1.20+ are recommended for stable np.divide behavior and modern GeoTIFF compression support.

pip install rasterio numpy

Compatibility Warning: Rasterio relies on GDAL under the hood. On macOS ARM64 or Windows, use conda install -c conda-forge rasterio to avoid DLL conflicts. Python 3.9–3.11 is fully supported; 3.12 requires rasterio ≥1.3.9 due to C-API changes. Consult the official Rasterio Quickstart for platform-specific build notes and GDAL wheel compatibility matrices.

2. Load Multispectral Stacks & Map Bands

Drone sensors (MicaSense, DJI P4 Multispectral, Parrot Sequoia) output band-ordered stacks or separate files. Always verify band indices before math. NDVI requires Red (typically Band 3) and NIR (Band 4). NDRE requires Red Edge (typically Band 5) and NIR.

import rasterio
import numpy as np
from pathlib import Path

def load_bands(src_path: str) -> tuple[np.ndarray, dict]:
    with rasterio.open(src_path) as src:
        # Read all bands into a 3D array: (bands, height, width)
        stack = src.read()
        meta = src.meta.copy()
        
        # Verify band count (1-based indexing in metadata vs 0-based in arrays)
        if stack.shape[0] < 5:
            raise ValueError("Input raster must contain at least 5 spectral bands.")
        return stack, meta

Band Mapping Note: Rasterio reads bands as 0-indexed arrays. If your sensor documentation lists Band 1 as Blue, Band 2 as Green, Band 3 as Red, Band 4 as NIR, and Band 5 as Red Edge, your array indices will be 2, 3, and 4 respectively. Always cross-reference your sensor’s spectral response curves before proceeding. Misaligned bands will produce inverted or ecologically meaningless vegetation indices.

3. Vectorized Band Math & Safe Division

Direct array subtraction and division trigger RuntimeWarning: invalid value encountered in divide when denominators hit zero. Use np.divide with the where parameter to safely compute indices. This approach aligns with standard practices in Band Math & Raster Algebra in Python for production-grade raster algebra.

def calculate_indices(stack: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    # Convert to float32 to prevent integer overflow during subtraction
    red = stack[2].astype(np.float32)      # Band 3
    nir = stack[3].astype(np.float32)      # Band 4
    red_edge = stack[4].astype(np.float32) # Band 5

    # Initialize output arrays with NaN to preserve nodata semantics
    ndvi = np.full_like(nir, np.nan, dtype=np.float32)
    ndre = np.full_like(nir, np.nan, dtype=np.float32)

    # Safe division: only compute where denominator > 0
    denom_ndvi = nir + red
    denom_ndre = nir + red_edge

    # See NumPy docs for out/where parameter usage: https://numpy.org/doc/stable/reference/generated/numpy.divide.html
    np.divide(nir - red, denom_ndvi, out=ndvi, where=denom_ndvi > 0)
    np.divide(nir - red_edge, denom_ndre, out=ndre, where=denom_ndre > 0)

    # Clip to valid vegetation index range [-1, 1]
    ndvi = np.clip(ndvi, -1.0, 1.0)
    ndre = np.clip(ndre, -1.0, 1.0)
    
    return ndvi, ndre

Why np.full_like + where? Initializing with NaN ensures that pixels outside the calculation mask (e.g., shadows, water, or sensor edge artifacts) remain explicitly undefined rather than defaulting to 0.0, which downstream GIS software would misinterpret as valid low-vegetation values.

4. Export Georeferenced Results

Writing indices back to disk requires updating the raster metadata. You must change the data type to float32, define a nodata value, and apply compression to keep file sizes manageable without sacrificing precision.

def export_indices(ndvi: np.ndarray, ndre: np.ndarray, meta: dict, out_dir: str):
    out_dir = Path(out_dir)
    out_dir.mkdir(exist_ok=True)

    # Update metadata for float output and compression
    meta.update({
        'count': 1,
        'dtype': 'float32',
        'compress': 'deflate',
        'nodata': np.nan
    })

    for name, data in {'ndvi.tif': ndvi, 'ndre.tif': ndre}.items():
        out_path = out_dir / name
        with rasterio.open(out_path, 'w', **meta) as dst:
            dst.write(data, 1)

Metadata Preservation: Rasterio automatically retains the original crs, transform, and width/height when you pass the copied meta dictionary. The compress: 'deflate' setting reduces file size by ~60% while maintaining lossless precision, which is critical for time-series change detection.

5. Production Considerations & Validation

For large agricultural fields or regional drone surveys, loading entire rasters into memory will cause MemoryError exceptions. Implement windowed reading using rasterio.windows to process imagery in tiles. Always validate outputs by checking histogram distributions; healthy vegetation typically shows NDVI values between 0.3 and 0.8, while NDRE ranges from 0.0 to 0.5 in mid-season crops.

When integrating these indices into broader pipelines, review the Drone Imagery Processing & Vegetation Index Workflows cluster for calibration, atmospheric correction, and time-series aggregation patterns. Proper radiometric calibration before index calculation ensures cross-flight comparability, which is critical for precision agriculture decision support.

Quick Validation Checklist:

• Verify band order against sensor spectral response curves
• Confirm np.divide masks zero-denominator pixels correctly
• Ensure output GeoTIFFs retain original CRS, transform, and resolution
• Apply NaN or -9999 nodata values consistently across exports
• Test with a small subset before scaling to full-field mosaics
• Validate against ground truth or calibrated reflectance targets