Map algebra, NDVI and reclassification

A field reflects 0.0494 of red light and 0.4214 of near-infrared: (0.4214 − 0.0494) ÷ (0.4214 + 0.0494) = 0.790, the NDVINDVI (normalised difference vegetation index): (NIR − red) ÷ (NIR + red), from −1 to 1: negative over water and rising with green leaves, which absorb red light and reflect near-infrared strongly. of dense, healthy crop. Map algebra repeats that formula for every cell of aligned rasters, 480,000 times for the lesson raster, and writes a new raster of the answers.

Here you combine bands cell by cell with a mask that keeps missing cells missing, compute NDVI and the two indices that share the name NDWINDWI (normalised difference water index): A name shared by two indices: McFeeters' (green − NIR) ÷ (green + NIR), which outlines open water, and Gao's (NIR − SWIR) ÷ (NIR + SWIR), which senses water in leaves. A map should say which bands it used., and read what their values say about water, soil, towns and crops. Then you reproduce band maths that goes wrong without an error, and turn NDVI into classes with thresholds you can defend.

Map algebra applies one formula to every cell of aligned grids

Map algebramap algebra: Computing a new raster from one or more aligned rasters with a formula applied cell by cell, such as a band ratio, a sum or a comparison with a threshold. Kinds that also read a cell's neighbours or a whole zone are moving-window and zonal operations. computes a raster from one or more rasters with a formula that works on each cell by itself: add two bands, divide one by another, compare a band with a threshold. The inputs must be aligned, with the same CRS, pixel size, origin and number of rows and columns, so that row 418, column 100 is the same patch of ground in each; the next lesson shows how to get there. Here is a corner of the lesson raster, three rows by three columns, through the NDVI formula:

Two cells are crop, 0.790 and 0.789, six are bare soil, about 0.12, and one lies in the hole without data. A cell missing in either input must stay missing in the output, whatever its nodata valuenodata value: A value, such as -9999, that marks a raster cell with no measurement. It is a code, not data, and must be left out of sums, means and counts. More in “Vector and raster data models”. Follow both kinds of cell through the steps:

The program applies those steps to all 480,000 pixels and averages the result for each kind of ground:

"""The lesson's raster: two bands, red and near-infrared, of 800 by 600 pixels of 10 m in UTM zone 43N.

Each value is a reflectance times 10,000, stored as an unsigned 16-bit integer, and 0 marks a pixel without data.
The scene is made by whole-number formulas, so every program that imports it, here or in a browser, gets the
same numbers: a river, fields of crops and bare soil, a town, and a hole of 2,400 pixels with no data.
"""

WIDTH, HEIGHT = 800, 600  # columns, rows
TRANSFORM = (10.0, 0.0, 370_000.0, 0.0, -10.0, 2_046_000.0)  # a, b, c, d, e, f
CRS = "EPSG:32643"  # WGS 84 / UTM zone 43N
NODATA = 0


def noise(row, col):
    """A fixed whole number from -40 to 40 for each row and column, the same on every computer."""
    h = ((row * 73_856_093) ^ (col * 19_349_663)) & 0xFFFFFFFF
    h = ((h ^ (h >> 13)) * 1_274_126_177) & 0xFFFFFFFF
    return h % 81 - 40


def wobble(row, col):
    """Texture that changes smoothly over 16 pixels, as real ground does, plus a little noise of its own."""
    r, c, fr, fc = row // 16, col // 16, row % 16, col % 16
    top = noise(r, c) * (16 - fc) + noise(r, c + 1) * fc
    bottom = noise(r + 1, c) * (16 - fc) + noise(r + 1, c + 1) * fc
    return (top * (16 - fr) + bottom * fr) // 256 + noise(row, col) // 8


def cover(row, col):
    """What lies under a pixel: water, town, soil, crop, or None in the hole without data."""
    if 420 <= row < 460 and 100 <= col < 160:
        return None
    if abs(row - (260 + col * col // 4000)) < 18:
        return "water"
    if 60 <= row < 200 and 520 <= col < 700:
        return "town"
    field_row, field_col = row // 80, col // 100
    if (field_row * 3 + field_col) % 5 == 0:
        return "soil"
    return "crop"


def pixel(row, col):
    """(red, nir) of one pixel, reflectance times 10,000."""
    kind = cover(row, col)
    if kind is None:
        return NODATA, NODATA
    w = wobble(row, col)
    if kind == "water":
        return 350 + w, 200 + w // 2
    if kind == "town":
        return 1_700 + w, 2_200 + w
    if kind == "soil":
        return 2_200 + w, 2_800 + w
    vigour = ((row // 80) * 7 + (col // 100) * 13) % 10  # 0 to 9, one value per field
    return 800 - vigour * 40 + w, 3_000 + vigour * 150 + (col % 100) * 3 + w


def bands():
    """The two bands as lists of rows: red, then near-infrared."""
    red, nir = [], []
    for row in range(HEIGHT):
        pairs = [pixel(row, col) for col in range(WIDTH)]
        red.append([p[0] for p in pairs])
        nir.append([p[1] for p in pairs])
    return red, nir

"""NDVI cell by cell: a 3 x 3 corner of the lesson raster, then every pixel, with nodata kept out."""
from scene import HEIGHT, NODATA, WIDTH, cover, pixel


def ndvi(red, nir):
    """(NIR - red) / (NIR + red) in floating point, or None where either band has no data."""
    if red == NODATA or nir == NODATA:
        return None
    return (nir - red) / (nir + red)


def show(value):
    return "nodata" if value is None else f"{value:.3f}"


rows, cols = range(418, 421), range(98, 101)
print("row  band" + "".join(f"  col {c:>3}" for c in cols))
for r in rows:
    print(f"{r}  red " + "".join(f"{pixel(r, c)[0]:>9}" for c in cols))
    print(f"{r}  nir " + "".join(f"{pixel(r, c)[1]:>9}" for c in cols))
    print(f"{r}  ndvi" + "".join(f"{show(ndvi(*pixel(r, c))):>9}" for c in cols))

totals, counts, missing = {}, {}, 0
for r in range(HEIGHT):
    for c in range(WIDTH):
        value = ndvi(*pixel(r, c))
        if value is None:
            missing += 1
            continue
        kind = cover(r, c)
        totals[kind] = totals.get(kind, 0.0) + value
        counts[kind] = counts.get(kind, 0) + 1

valid = sum(counts.values())
print(f"pixels with data: {valid:,}; nodata: {missing:,}")
print(f"mean NDVI of the pixels with data: {sum(totals.values()) / valid:.4f}")
for kind in ("water", "town", "soil", "crop"):
    print(f"  {kind:>5}: {counts[kind]:>7,} pixels, mean NDVI {totals[kind] / counts[kind]:.3f}")

row  band  col  98  col  99  col 100
418  red      2217     2215      494
418  nir      2817     2815     4214
418  ndvi    0.119    0.119    0.790
419  red      2210     2209      496
419  nir      2810     2809     4216
419  ndvi    0.120    0.120    0.789
420  red      2209     2210        0
420  nir      2809     2810        0
420  ndvi    0.120    0.120   nodata
pixels with data: 477,600; nodata: 2,400
mean NDVI of the pixels with data: 0.5046
  water:  28,000 pixels, mean NDVI -0.273
   town:  25,200 pixels, mean NDVI 0.128
   soil:  94,900 pixels, mean NDVI 0.120
   crop: 329,500 pixels, mean NDVI 0.710

477,600 pixels have data, and their mean NDVI is 0.5046: water −0.273, town 0.128, soil 0.120 and crop 0.710. On whole arrays, NumPy’s masked arrays do the same bookkeeping: a True in the mask marks a value as invalid, and computations such as a mean leave it out (NumPy, numpy.ma). rasterio returns such an array from read(masked=True).

NDVI and NDWI measure the contrast between two bands

Each index here is a normalised difference of two bands, (a − b) ÷ (a + b), which lies between −1 and 1 whenever neither band is negative. NDVI compares near-infrared with red: the chlorophyll of green leaves takes in most of the visible light for photosynthesis, while the cells of the leaf send much of the near-infrared back (NASA Earth Observatory). It was described for monitoring vegetation on the Great Plains with ERTS images (Rouse and others). NASA reads 0.1 and below as rock, sand or snow, 0.2 to 0.3 as shrub and grassland, and 0.6 to 0.8 as temperate and tropical rainforest.

NDWI names two indices. McFeeters defined one from green and near-infrared to outline open water (McFeeters); Gao defined another from near-infrared and shortwave infrared to sense the liquid water in leaves (Gao). On Sentinel-2, green is B3, red B4, near-infrared B8 and shortwave infrared B11, at about 1,610 nm (SentiWiki). The program computes all three for four surfaces, with illustrative reflectances:

"""Three normalised differences for four surfaces: NDVI and the two indices called NDWI."""

# Reflectance, 0 to 1, in four bands: illustrative values for this lesson, shaped like real surfaces
# (green 560 nm, red 665 nm, near-infrared 833 nm, shortwave infrared 1610 nm, about Sentinel-2's B3, B4, B8, B11).
SURFACES = {
    "water": {"green": 0.060, "red": 0.035, "nir": 0.020, "swir": 0.010},
    "crop": {"green": 0.080, "red": 0.050, "nir": 0.400, "swir": 0.200},
    "soil": {"green": 0.150, "red": 0.220, "nir": 0.280, "swir": 0.350},
    "town": {"green": 0.140, "red": 0.170, "nir": 0.220, "swir": 0.250},
}


def normalised_difference(a, b):
    """(a - b) / (a + b): from -1 to 1, the contrast between two bands."""
    return (a - b) / (a + b)


print(f"{'surface':>8} {'NDVI':>7} {'NDWI (green, NIR)':>18} {'NDWI (NIR, SWIR)':>17}")
for name, s in SURFACES.items():
    ndvi = normalised_difference(s["nir"], s["red"])
    ndwi_water = normalised_difference(s["green"], s["nir"])  # McFeeters (1996): open water
    ndwi_leaves = normalised_difference(s["nir"], s["swir"])  # Gao (1996): water in leaves
    print(f"{name:>8} {ndvi:>7.3f} {ndwi_water:>18.3f} {ndwi_leaves:>17.3f}")

 surface    NDVI  NDWI (green, NIR)  NDWI (NIR, SWIR)
   water  -0.273              0.500             0.333
    crop   0.778             -0.667             0.333
    soil   0.120             -0.302            -0.111
    town   0.128             -0.222            -0.064

Only water has a positive green and near-infrared NDWI, 0.500, and only water a negative NDVI. The near-infrared and shortwave index gives water and crop the same 0.333: it was made for water in leaves, not for finding ponds. Name the two bands of every NDWI you publish.

Reclassificationreclassification: Turning a raster of continuous values into classes with a table of thresholds, such as NDVI below 0 as water and 0.3 or more as vegetation. The table and the image's date belong with the result. turns such an index into classes with a table of thresholds, such as NDVI below 0 as water and 0.3 or more as vegetation. Its thresholds are decisions about one place and one date, so a reclassification is only as good as the table kept with it.

Band maths that fails without an error, reproduced and fixed Pro: opens with a Pro, Premium or Ultimate pass

Reclassification with thresholds you can defend and record Pro: opens with a Pro, Premium or Ultimate pass

Practise: NDVI with a mask

The exercise computes NDVI with a mask on the lesson’s corner and averages it, and the self-check covers the ideas of the first two sections, from cell-by-cell algebra to the two NDWIs:

def ndvi_grid(red, nir, nodata):
    """NDVI of every pixel, a list of rows, with None where either band has no data or the sum is 0."""
    # Replace this line: go through the rows and columns of both bands together.
    return []


def mean_ndvi(grid):
    """(mean rounded to 4 decimals, count) of the values that are not None; (None, 0) if there are none."""
    # Replace this line: leave out every None before you sum.
    return None, 0

from masked import mean_ndvi, ndvi_grid

RED = [[2217, 2215, 494], [2210, 2209, 496], [2209, 2210, 0]]
NIR = [[2817, 2815, 4214], [2810, 2809, 4216], [2809, 2810, 0]]


def test_lesson_corner():
    """computes NDVI cell by cell and keeps the nodata cell as None"""
    grid = ndvi_grid(RED, NIR, 0)
    assert round(grid[0][2], 3) == 0.790
    assert round(grid[1][0], 3) == 0.120
    assert grid[2][2] is None


def test_shape_is_kept():
    """returns a grid of the same shape as the bands"""
    grid = ndvi_grid(RED, NIR, 0)
    assert len(grid) == 3 and all(len(row) == 3 for row in grid)


def test_mean_leaves_nodata_out():
    """averages the eight cells with data, not nine"""
    assert mean_ndvi(ndvi_grid(RED, NIR, 0)) == (0.287, 8)


def test_zero_sum_and_other_nodata_value():
    """masks a pixel whose bands are both 0, and honours another nodata value"""
    grid = ndvi_grid([[100, 0], [-9999, 300]], [[300, 0], [500, -9999]], -9999)
    assert grid[0][0] == 0.5
    assert grid[0][1] is None and grid[1][0] is None and grid[1][1] is None
    assert mean_ndvi(grid) == (0.5, 1)


def test_empty():
    """returns (None, 0) when no cell has data"""
    assert mean_ndvi([[None, None]]) == (None, 0)


def test_equal_bands():
    """gives 0, not nodata, to a pixel whose two bands are equal"""
    assert ndvi_grid([[200, 0]], [[200, 0]], 0) == [[0.0, None]]

Interview questions

Three longer questions start from this lesson, from an index above 1 to a threshold that moves an answer; the bank gives what each tests and a first step:

Key takeaways

• Map algebra applies one formula to the same cell of aligned rasters, and a cell missing in any input stays missing in the output.
• NDVI, (NIR − red) ÷ (NIR + red), runs from −1 to 1: negative for water, about 0.1 for bare ground and 0.6 to 0.8 for rainforest.
• Two indices are called NDWI: green with near-infrared outlines open water, near-infrared with shortwave infrared senses water in leaves.
• Reclassification turns an index into classes through a table of thresholds, which must be kept with the result.