Computing EVI and NDWI from Sentinel-2 Bands

To compute EVI and NDWI from Sentinel-2 L2A bands, load the blue, green, red, and NIR bands as xarray DataArrays, scale the integer digital numbers to reflectance, then apply each index formula — EVI needs the blue band and its four coefficients, NDWI is a green/NIR normalized difference:

import xarray as xr

s = {b: xr.open_dataarray(f"{b}.tif") * 1e-4 for b in ("B02", "B03", "B04", "B08")}
evi = 2.5 * (s["B08"] - s["B04"]) / (s["B08"] + 6 * s["B04"] - 7.5 * s["B02"] + 1)
ndwi = (s["B03"] - s["B08"]) / (s["B03"] + s["B08"])

These two indices go beyond the single red/NIR ratio and show why coefficient and band handling matters in Band Math Operations with xarray.

Why This Arises in Remote Sensing Workflows

NDVI is the default vegetation index, but it has two well-known weaknesses: it saturates over dense canopy and it is sensitive to atmospheric aerosols and soil background. The Enhanced Vegetation Index was designed to fix both — it adds a blue-band aerosol correction and a soil-adjustment term, which keeps it responsive in high-biomass forests where NDVI flattens out. The Normalized Difference Water Index answers a different question entirely: where is open water. Together they cover vegetation vigour and surface water in a single band-math pass over the same Sentinel-2 scene.

This page deliberately does not re-derive the plain red/NIR ratio — that is covered in Calculating NDVI directly from xarray DataArrays. The focus here is the two things that trip people up when they move past NDVI: EVI’s four coefficients (G, C1, C2, L) with its extra blue band, and picking the correct NDWI formulation. Both depend on getting Sentinel-2’s band identifiers and reflectance scaling right first.

The work sits inside Band Math Operations with xarray and the wider Core Raster Fundamentals & STAC Mapping section. When you productionise these formulas across many scenes and indices, they graduate into a Spectral Index Calculation Pipelines workflow.

Sentinel-2 Band Mapping and the Data Flow

Both indices draw from the 10 m Sentinel-2 bands. The diagram maps each band identifier to its role and shows which inputs feed each formula.

Sentinel-2 band inputs to EVI and NDWI Sentinel-2 bands B02 blue, B03 green, B04 red, and B08 near-infrared are scaled to reflectance, then routed to two index formulas. EVI consumes blue, red, and NIR with coefficients G, C1, C2, L. NDWI consumes green and NIR as a normalized difference. Sentinel-2 L2A B02 · Blue B03 · Green B04 · Red B08 · NIR scale x 1e-4 → reflectance EVI G(NIR−Red) ÷ (NIR+C1·Red−C2·Blue+L) G=2.5 C1=6 C2=7.5 L=1 NDWI (Green−NIR) ÷ (Green+NIR) McFeeters, open water Blue is unique to EVI; NDWI reuses Green and NIR only.

Note that B08 (NIR) is shared by both formulas, blue (B02) is exclusive to EVI, and green (B03) is exclusive to NDWI. Loading all four bands once and reusing the NIR array keeps the computation lean.

Environment & Setup

Package Minimum version Why required
xarray 2023.1 DataArray broadcasting for the index arithmetic
rioxarray 0.13 Opens the band GeoTIFFs as georeferenced DataArrays with CRS/nodata
numpy 1.23 Backing array library and NaN handling 
pip install "xarray>=2023.1" "rioxarray>=0.13" "numpy>=1.23"

Sentinel-2 L2A bands are 16-bit integers. Open them with masked=True so the product’s nodata (0) is promoted to NaN, which stops fill pixels from contaminating the ratios.

Complete Working Example

The script loads the four required bands, applies the L2A offset and scale, and computes both indices with proper NaN propagation. EVI is clipped to its meaningful range because the aerosol term can overshoot on noisy pixels.

import numpy as np
import rioxarray
import xarray as xr

# EVI coefficients (Huete et al., MODIS/Sentinel-2 standard)
G, C1, C2, L = 2.5, 6.0, 7.5, 1.0
# L2A processing baseline >= 04.00 adds a -1000 DN offset before scaling
BOA_ADD_OFFSET = -1000.0
SCALE = 1e-4


def load_reflectance(path: str) -> xr.DataArray:
    """Open a Sentinel-2 L2A band as 0-1 surface reflectance with NaN nodata."""
    da = rioxarray.open_rasterio(path, masked=True).squeeze("band", drop=True)
    # Apply offset then scale; masked=True already made nodata NaN
    return (da + BOA_ADD_OFFSET) * SCALE


def compute_evi_ndwi(
    blue_path: str, green_path: str, red_path: str, nir_path: str
) -> xr.Dataset:
    """
    Compute EVI and NDWI from Sentinel-2 B02/B03/B04/B08 GeoTIFFs.

    Returns an xr.Dataset with two variables: 'evi' and 'ndwi'.
    """
    blue = load_reflectance(blue_path)   # B02
    green = load_reflectance(green_path) # B03
    red = load_reflectance(red_path)     # B04
    nir = load_reflectance(nir_path)     # B08

    # EVI needs blue for the aerosol term and L for soil adjustment
    evi = G * (nir - red) / (nir + C1 * red - C2 * blue + L)
    # Physically meaningful EVI is within [-1, 1]; clip aerosol overshoot
    evi = evi.clip(-1.0, 1.0)

    # NDWI (McFeeters) for open water; where() avoids divide-by-zero
    denom = green + nir
    ndwi = xr.where(denom != 0, (green - nir) / denom, np.nan)

    evi.name, ndwi.name = "evi", "ndwi"
    ds = xr.merge([evi, ndwi])
    # Preserve georeferencing from an input band
    ds = ds.rio.write_crs(red.rio.crs)
    return ds


if __name__ == "__main__":
    result = compute_evi_ndwi(
        "B02.tif", "B03.tif", "B04.tif", "B08.tif"
    )
    print("EVI  range:", float(result.evi.min()), "→", float(result.evi.max()))
    print("NDWI range:", float(result.ndwi.min()), "→", float(result.ndwi.max()))
    # Water pixels are where NDWI > 0
    water_fraction = float((result.ndwi > 0).mean())
    print(f"water fraction: {water_fraction:.3%}")
    result.evi.rio.to_raster("evi.tif")
    result.ndwi.rio.to_raster("ndwi.tif")

The (da + BOA_ADD_OFFSET) * SCALE line is the step most tutorials omit. Since Sentinel-2 processing baseline 04.00 (early 2022), L2A products carry a BOA_ADD_OFFSET of −1000 DN; skipping it biases reflectance high and pushes EVI values out of range.

Variant Patterns

1. Reading bands from a single stacked Dataset

If your bands already live in one multi-band file or an xarray.Dataset (for example from a STAC stack), select bands by name instead of opening separate files.

import xarray as xr

# ds has data variables named by band id: B02, B03, B04, B08
ds = xr.open_dataset("s2_stack.nc")
refl = ds[["B02", "B03", "B04", "B08"]].astype("float32") * 1e-4

evi = 2.5 * (refl.B08 - refl.B04) / (refl.B08 + 6 * refl.B04 - 7.5 * refl.B02 + 1)
ndwi = (refl.B03 - refl.B08) / (refl.B03 + refl.B08)

2. Applying a scene classification mask before indexing

Compute indices only over valid land/water using the L2A Scene Classification (SCL) band, so clouds and shadows do not pollute statistics.

import rioxarray
import xarray as xr

scl = rioxarray.open_rasterio("SCL.tif").squeeze("band", drop=True)
# Keep vegetation(4), bare(5), water(6), unclassified(7); drop cloud/shadow classes
valid = scl.isin([4, 5, 6, 7])

blue = rioxarray.open_rasterio("B02.tif", masked=True).squeeze("band", drop=True) * 1e-4
red = rioxarray.open_rasterio("B04.tif", masked=True).squeeze("band", drop=True) * 1e-4
nir = rioxarray.open_rasterio("B08.tif", masked=True).squeeze("band", drop=True) * 1e-4

evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)
evi_clean = evi.where(valid)  # cloud/shadow pixels become NaN

3. Vectorised multi-index function

Return a dictionary of indices so a downstream pipeline can request any subset without recomputing shared bands.

import xarray as xr

def spectral_indices(refl: xr.Dataset) -> dict[str, xr.DataArray]:
    b02, b03, b04, b08 = refl.B02, refl.B03, refl.B04, refl.B08
    return {
        "evi": (2.5 * (b08 - b04) / (b08 + 6 * b04 - 7.5 * b02 + 1)).clip(-1, 1),
        "ndwi": (b03 - b08) / (b03 + b08),
        "ndvi": (b08 - b04) / (b08 + b04),  # cheap to include; shares B08/B04
    }

This shape is exactly what a configurable Spectral Index Calculation Pipelines workflow consumes.

Common Errors

EVI values cluster near zero or explode past ±1

The bands were not scaled to reflectance, so raw digital numbers (0–10000) entered the formula and the + L soil term became negligible. Multiply by 1e-4 (and apply the −1000 offset for baseline ≥ 04.00) before computing EVI, then clip to [-1, 1].

ValueError: operands could not be broadcast together

The bands have different shapes — usually a 20 m band mixed with 10 m bands. Resample them to a common grid first; see Resampling Sentinel-2 20m Bands to 10m.

NDWI flags vegetation as water

You used the Gao (NIR − SWIR)/(NIR + SWIR) formula, which measures canopy water content, not the McFeeters (Green − NIR)/(Green + NIR) formula for open water. Use B03 and B08 for open-water delineation.