Computing NDWI and MNDWI for Water Mapping

Two indices map surface water, and the difference between them decides how the result behaves over towns:

import xarray as xr

def nd(a, b):
    denom = a + b
    return xr.where(denom != 0, (a - b) / denom, float("nan"))

ndwi = nd(green, nir)        # B03, B08 — 10 m, confuses built-up surfaces with water
mndwi = nd(green, swir16)    # B03, B11 — 20 m, separates them cleanly

Both are index-pipeline members in the sense of Spectral Index Calculation Pipelines; the choice between them is the interesting part.

Why This Arises in Remote Sensing Workflows

Water mapping is one of the most operationally useful products remote sensing produces: flood extent, reservoir monitoring, wetland change, irrigation tracking. It is also one where a simple index gets you most of the way, which is why it is so widely used — and why its two failure modes are so widely encountered.

The first is built-up confusion. NDWI compares green against near-infrared, and water is dark in the near-infrared. Unfortunately so is asphalt, and shadowed roofs are darker still, so a city renders as a lake in an NDWI map. MNDWI substitutes the short-wave infrared band, where water remains dark but concrete and asphalt are bright, which separates the two decisively.

The second is the edge. Water boundaries are almost never aligned to pixels, so shoreline pixels are mixtures, and their index value depends on the water fraction. That makes the mapped extent sensitive to the threshold in exactly the places where the map is used — the shoreline — and it is why area figures from water maps need an uncertainty statement rather than a single number.

Why the SWIR variant separates water from concrete Water is dark in both near-infrared and short-wave infrared. Asphalt is also dark in the near-infrared, which is why NDWI confuses it with water, but it is bright in the short-wave infrared, which is what MNDWI exploits to separate the two. Reflectance by surface and band surface B03 green B08 NIR B11 SWIR open water 0.06 0.02 0.01 asphalt / concrete 0.10 0.12 — looks like water to NDWI 0.24 — clearly not water vegetation 0.08 0.35 0.18 bare soil 0.14 0.28 0.34 NDWI(asphalt) ≈ −0.09 — close enough to zero to cross the threshold in shadow. MNDWI(asphalt) ≈ −0.41 — comfortably below any sensible cut.

Environment & Setup

Package Version Why
rioxarray ≥0.15 Masked reads and grid matching
xarray ≥2023.6 Index arithmetic
rasterio ≥1.3.0 Resampling and I/O
scikit-image ≥0.22 Cleaning the binary water map 
pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0" "scikit-image>=0.22"

Complete Working Example

This function computes both indices on a common grid, thresholds MNDWI with a value derived from the scene, and cleans the result.

Water look-alikes and how each is removed Asphalt, terrain shadow, wet soil and cloud shadow all approach the water threshold for different reasons, and each has a different defence: the SWIR band, an absolute brightness floor, a minimum patch size, and the quality mask respectively. Four things that look like water, and the fix look-alike why it crosses the threshold defence asphalt and concrete dark in the near-infrared use MNDWI, not NDWI terrain shadow dark in every band exclude below a green reflectance floor wet soil after rain genuinely between soil and water minimum patch size + persistence cloud shadow dark, and not surface at all apply the quality mask first Only the first is solved by the index choice; the rest need the surrounding workflow. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling
from skimage.morphology import remove_small_holes, remove_small_objects


def nd(a: xr.DataArray, b: xr.DataArray) -> xr.DataArray:
    denom = a + b
    return xr.where(denom != 0, (a - b) / denom, np.nan).astype("float32")


def water_map(
    green_path: str,
    nir_path: str,
    swir_path: str,
    scl_path: str,
    *,
    threshold: float | None = None,
    min_object_px: int = 20,
    min_hole_px: int = 20,
) -> tuple[xr.DataArray, dict]:
    """MNDWI-based water map on the SWIR band's own 20 m grid."""
    swir = rioxarray.open_rasterio(swir_path, masked=True).squeeze(drop=True).astype("float32")
    green = rioxarray.open_rasterio(green_path, masked=True).squeeze(drop=True).astype("float32")
    nir = rioxarray.open_rasterio(nir_path, masked=True).squeeze(drop=True).astype("float32")
    scl = rioxarray.open_rasterio(scl_path).squeeze(drop=True)

    # Work on the coarser grid: downsampling 10 m green loses nothing a water map needs
    green = green.rio.reproject_match(swir, resampling=Resampling.average)
    nir = nir.rio.reproject_match(swir, resampling=Resampling.average)
    scl = scl.rio.reproject_match(swir, resampling=Resampling.nearest)

    clear = ~scl.isin([0, 1, 3, 8, 9, 10, 11])
    green, nir, swir = green.where(clear), nir.where(clear), swir.where(clear)

    mndwi = nd(green, swir)
    ndwi = nd(green, nir)

    if threshold is None:
        # Otsu-style split is reasonable here: water scenes are genuinely bimodal
        finite = mndwi.values[np.isfinite(mndwi.values)]
        try:
            from skimage.filters import threshold_otsu
            threshold = float(threshold_otsu(finite)) if finite.size else 0.0
        except Exception:
            threshold = 0.0

    water = (mndwi > threshold).fillna(False).values
    water = remove_small_objects(water, min_size=min_object_px)
    water = remove_small_holes(water, area_threshold=min_hole_px)

    out = xr.DataArray(water.astype("uint8"), coords=mndwi.coords, dims=mndwi.dims, name="water")
    out.rio.write_crs(mndwi.rio.crs, inplace=True)

    px_area_ha = abs(np.prod(mndwi.rio.resolution())) / 10_000
    stats = {
        "threshold": round(float(threshold), 4),
        "water_pixels": int(water.sum()),
        "water_area_ha": round(float(water.sum()) * px_area_ha, 2),
        "ndwi_water_pixels": int((ndwi > threshold).fillna(False).values.sum()),
        "resolution_m": abs(mndwi.rio.resolution()[0]),
    }
    return out, stats


if __name__ == "__main__":
    water, stats = water_map("B03.tif", "B08.tif", "B11.tif", "SCL.tif")
    print(stats)
    print("NDWI flags",
          round(stats["ndwi_water_pixels"] / max(stats["water_pixels"], 1), 2),
          "× as many pixels as MNDWI — the difference is mostly built-up surface")
    water.rio.to_raster("water_mask.tif", driver="COG", compress="DEFLATE", dtype="uint8")

Computing at 20 m rather than upsampling the SWIR band is a deliberate choice: it avoids inventing detail, and a water map at 20 m is adequate for almost every reporting purpose. If a 10 m product is required, upsample the SWIR band with bilinear and state that the effective resolution is 20 m, following Matching Landsat and Sentinel-2 Grids.

Variant Patterns

1. Choosing the threshold from the histogram

Where the threshold sits in a bimodal scene A scene containing a substantial water body produces two modes: land below zero and water above. The minimum between them is the natural cut, and it is often slightly below zero rather than exactly at zero. Turbid and shallow water sits on the left flank of the water mode, so a threshold at zero systematically under-maps it. MNDWI histogram, scene with a reservoir land water histogram minimum MNDWI = 0 −0.8 0.0 +0.8 Between the two lines sit shallow and turbid pixels — the difference is often several percent of mapped area.

2. Multi-date persistence to remove transient wet ground

import xarray as xr

stack = xr.concat([water_map(**scene)[0] for scene in scenes], dim="time")
frequency = stack.mean("time")                 # fraction of dates each pixel was water

permanent = frequency > 0.8
seasonal = (frequency > 0.2) & (frequency <= 0.8)

Water frequency is a far more useful product than any single-date extent, and it separates permanent water bodies from flooding and irrigation — the temporal reasoning developed in Temporal Aggregation and Time-Series Analysis.

3. Terrain shadow, the remaining false positive

In mountainous terrain, deep shadow produces low reflectance in every band, and the normalised difference of two small numbers is unstable. Two defences work: exclude pixels below an absolute reflectance floor in the green band, and exclude slopes above a threshold using a DEM aligned to the imagery.

too_dark = green < 0.02                      # deep shadow: index value is not meaningful
water = water & ~too_dark.values

Validating a Water Map

Three checks, in increasing order of effort.

Known water bodies must be mapped. Pick a permanent reservoir or lake, and confirm that its interior is fully classified and that the mapped area is within a few percent of its known extent. Interior holes usually mean the threshold is too high; overspill usually means shadow contamination.

Known dry areas must not be. Urban centres, airports and quarries are the classic false positives, and they are the fastest way to detect that NDWI was used where MNDWI was needed.

Edge behaviour needs a sensitivity test. Recompute the area at threshold ±0.05 and report the range: for a compact reservoir it will be small, for a shallow floodplain it can be tens of percent, and that spread is the honest uncertainty of the product — the same reporting discipline argued for in Thresholding Change Maps and Removing Noise.

Common Errors

The city is mapped as a lake

NDWI was used in a built-up scene. Switch to MNDWI; the SWIR band is what separates asphalt from water.

The water body has holes in the middle

Turbid or shallow water sits below the threshold. Lower the cut using the histogram minimum, and fill small holes morphologically.

Cloud shadow is mapped as water

The quality mask was not applied before the index. Mask first — shadow is dark in every band and will always cross a water threshold.

Frequently Asked Questions

Q: NDWI or MNDWI? MNDWI in almost any scene containing settlements, because concrete and asphalt have low near-infrared reflectance and are misclassified as water by NDWI. NDWI is adequate in purely rural or natural catchments and has the advantage of staying on the 10 m grid.

Q: Is zero the right threshold? It is the usual starting point but rarely optimal. Turbid or shallow water sits below zero and sub-pixel mixing shifts the edge, so deriving the threshold from the histogram of a scene containing both water and land is more reliable.

Q: Why does wet soil appear as water? Saturated soil after rain has a spectral response between soil and water, so it sits near the threshold. Requiring a minimum patch size and checking persistence across dates separates transient wet ground from actual water bodies.