Map Algebra

1. What is Map Algebra?¶

In the previous chapters, you learned how to read spatial grids into memory and visualize them using cartographic techniques. Now, we transition from visualization to calculation.

Map algebra is the process of performing mathematical operations pixel by pixel across one or multiple raster grids. Because Rasterio seamlessly converts geographic bands into standard NumPy arrays, performing map algebra in Python is identical to performing matrix mathematics. If you have two perfectly aligned arrays representing different parts of the light spectrum, you can add, subtract, multiply, or divide their values to generate a brand new surface. This is how you can calculate derived environmental indices, track deforestation, and map flood boundaries.

2. Calculating Indices¶

Raw satellite data records the amount of light reflected by the Earth. By mathematically combining specific bands, we can isolate and highlight distinct environmental features.

The Normalized Difference Snow Index¶

To map snow and ice coverage, scientists use the Normalized Difference Snow Index (NDSI). Snow reflects visible light very strongly (which is why it looks bright white to our eyes) but absorbs Shortwave Infrared (SWIR) light. By comparing the Green band and the SWIR band, we can easily isolate snow from clouds and bare rock.

The mathematical formula for NDSI is:

However, calculating the raw index is only the first step. Because water also features high NDSI values (it has higher reflectance in the visible range than in the NIR and SWIR), we must exclude water bodies to prevent false positives. We do this by utilizing the NIR band. Water strongly absorbs NIR light, so we can isolate dry land and snow by only keeping pixels with NIR values greater than 0.2. Finally, to clearly separate snow and ice from everything else, we filter for NDSI values higher than 0.6 to generate a clean, binary snow map.

To calculate this in Python, we first read the Green, NIR, and SWIR2 bands into separate, perfectly aligned NumPy arrays.

Suppressing Division Warnings¶

When performing division across millions of pixels, you will inevitably encounter “divide by zero” or “invalid value” warnings. This occurs because the background areas outside the satellite image footprint contain zeros (NoData values). Dividing zero by zero produces a warning and returns NaN (NaN).

We can handle this cleanly using a NumPy context manager, which temporarily suppresses these expected warnings during the calculation.

While the raw NDSI calculation yields a continuous array ranging from -1 to 1 (where values closer to 1 represent confident snow cover), our conditional filters converted this into a binary mask. The final snow map contains only values of 0 (no snow/ice) and 1 (snow/ice), which makes it incredibly easy to quantify the total snow-covered area in subsequent analyses.

Concept Check: The Division Dilemma¶

Scenario: When calculating a spectral index like the NDSI across millions of pixels, you wrap your mathematical formula inside with np.errstate(divide='ignore', invalid='ignore'):. What is the primary purpose of this context manager?

A) To automatically delete any pixels that contain clouds or shadows from the dataset.

B) To prevent Python from filling the console with warnings or crashing when the formula inevitably attempts to divide by zero in the empty (NoData) background areas of the raster.

C) To mathematically force all negative index values to become positive numbers.

3. Change Detection¶

Map algebra is not limited to combining bands from a single image. If you have two rasters representing the exact same location at different times, you can subtract one from the other to calculate environmental change.

Let us look at an example of deforestation in Indonesia using Landsat 8 data from 2019 and 2023.

The Enhanced Vegetation Index (EVI) is heavily used to monitor biomass and canopy structure. The formula is slightly more complex than simple normalized differences:

To detect deforestation, we calculate the EVI for the older image, calculate the EVI for the newer image, and subtract the two. To verify that our math aligns with reality, it is highly useful to visualize this calculated EVI change alongside False Color Composites (FCCs) of the original dates.

Assuming we have already loaded our arrays and calculated the evi_change, we can build a three-panel comparison plot.

Deep red areas in this output clearly identify zones where active deforestation has removed the vegetation canopy over the four-year period.

Calculating Spectral Distance¶

While examining a specific index like EVI is highly useful for targeting vegetation, sometimes you want to know the total overall change across the entire spectrum. To do this, we calculate the Euclidean distance between the spectral signatures of the two dates across all available bands.

This calculation treats each pixel’s spectral signature as a point in multidimensional space and measures exactly how far that point has moved between 2019 and 2023.

Spectral distance serves as a powerful general change mask, highlighting any pixel that looks fundamentally different today than it did in the past, regardless of whether that change was driven by deforestation, urbanization, or flooding.

4. Exercise: The Water Index¶

In this exercise, you will apply map algebra to map a major flood event related to the filling of the Grand Ethiopian Renaissance Dam (GERD).

You are provided with two 3 band composite images: GERD_before_composite_SWIR1_NIR_Green.tif and GERD_after_composite_SWIR1_NIR_Green.tif. The filenames indicate the band order. Band 1 is Shortwave Infrared (SWIR1), Band 2 is Near Infrared (NIR), and Band 3 is Green.

Tasks:

1. Open both datasets and read the Green and NIR bands as float arrays.

2. Calculate the Normalized Difference Water Index (NDWI) for both the “before” and “after” periods. The formula is:

3. Calculate the NDWI difference by subtracting the “before” array from the “after” array.

4. Create a binary flood mask. A pixel represents new water (a flood) if the NDWI difference is greater than 0.2.

5. To visually verify your math, generate False Color Composites (using the SWIR1, NIR, and Green bands) for both time periods. Plot these alongside your binary flood mask in a three panel side by side comparison.

Sample Solution

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# 1. Read the Green and NIR data for both periods
with rasterio.open("data/GERD_before_composite_SWIR1_NIR_Green.tif") as src_before:
    # Band 2 is NIR, Band 3 is Green
    nir_before = src_before.read(2).astype(float)
    green_before = src_before.read(3).astype(float)

with rasterio.open("data/GERD_after_composite_SWIR1_NIR_Green.tif") as src_after:
    nir_after = src_after.read(2).astype(float)
    green_after = src_after.read(3).astype(float)

# 2 & 3. Calculate NDWI and the Difference (suppressing NoData warnings)
with np.errstate(divide='ignore', invalid='ignore'):
    ndwi_before = (green_before - nir_before) / (green_before + nir_before)
    ndwi_after = (green_after - nir_after) / (green_after + nir_after)

ndwi_diff = ndwi_after - ndwi_before

# 4. Create the binary flood mask
flood_mask = ndwi_diff > 0.2

# 5. Prepare False Color Composites to verify the results visually
def normalize_fcc(array, vmin=0.0, vmax=0.3):
    """Normalize and clip array to a specific range for display."""
    clipped = np.clip(array, vmin, vmax)
    return (clipped - vmin) / (vmax - vmin)

with rasterio.open("data/GERD_before_composite_SWIR1_NIR_Green.tif") as src:
    fcc_before = np.dstack((
        normalize_fcc(src.read(1).astype(float)),
        normalize_fcc(src.read(2).astype(float)),
        normalize_fcc(src.read(3).astype(float))
    ))

with rasterio.open("data/GERD_after_composite_SWIR1_NIR_Green.tif") as src:
    fcc_after = np.dstack((
        normalize_fcc(src.read(1).astype(float)),
        normalize_fcc(src.read(2).astype(float)),
        normalize_fcc(src.read(3).astype(float))
    ))

# Plotting the three panel comparison
fig, axes = plt.subplots(1, 3, figsize=(15, 6))

axes[0].imshow(fcc_before)
axes[0].set_title("Before FCC", fontsize=15)
axes[0].axis('off')

axes[1].imshow(fcc_after)
axes[1].set_title("After FCC", fontsize=15)
axes[1].axis('off')

axes[2].imshow(flood_mask, cmap='Blues')
axes[2].set_title("Extracted Flood Extent", fontsize=15)
axes[2].axis('off')

plt.tight_layout()
plt.show()

5. Summary: Raster Math¶

Map algebra relies on the fundamental connection between spatial layers and NumPy matrices. By mastering these operations, you unlock the ability to generate completely new data from raw imagery.

• Pixel-by-Pixel Mathematics: When spatial datasets share the same extent and resolution, you can apply basic arithmetic operators directly across their arrays.

• Derived Indices: Formulas like NDSI, EVI, and NDWI combine different parts of the electromagnetic spectrum to isolate specific environmental variables like snow, biomass, and water.

• Warning Management: Always convert raw integer data to floats before division, and use np.errstate to suppress unavoidable NoData division warnings.

• Temporal Change Detection: Subtracting indices from two different time periods reveals precise geographic patterns of loss and gain.

• Spectral Distance: Calculating the Euclidean distance across multiple bands provides a robust, index-agnostic measure of overall environmental change.