Rasterio with Real Satellite Data: Landsat 8 NDVI and Zonal Statistics for Bengaluru

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
import matplotlib.pyplot as plt

# Step 1: open the red band and read its basic facts (the metadata)
with rasterio.open("bengaluru_landsat/SR_B4.TIF") as src:
    print("size:", src.width, "x", src.height, "cells")
    print("cell size:", src.res, "metres")
    print("data type:", src.dtypes[0])
    print("CRS:", src.crs.to_epsg())
    print("NoData value:", src.nodata)
    print("bounds:", [round(v) for v in src.bounds])

# Step 2: the same for the near infrared band. Both bands must match: same size, same cells, same CRS.
with rasterio.open("bengaluru_landsat/SR_B5.TIF") as src:
    print("NIR size:", src.width, "x", src.height, "| same grid as red:", src.transform == rasterio.open("bengaluru_landsat/SR_B4.TIF").transform)

size: 1304 x 1274 cells
cell size: (30.0, 30.0) metres
data type: uint16
CRS: 32643
NoData value: 0.0
bounds: [764925, 1418145, 804045, 1456365]
NIR size: 1304 x 1274 | same grid as red: True

# Step 1: read both bands as they are stored
with rasterio.open("bengaluru_landsat/SR_B4.TIF") as src:
    red_dn = src.read(1)
    profile = src.profile                      # keep the file settings, for saving a result later
with rasterio.open("bengaluru_landsat/SR_B5.TIF") as src:
    nir_dn = src.read(1)

print("red digital numbers:", red_dn.min(), "to", red_dn.max())

# Step 2: apply the scale factor WITHOUT checking the valid range
raw_red = red_dn.astype("float32") * 0.0000275 - 0.2
print("reflectance without the valid range:", round(float(raw_red.min()), 3), "to", round(float(raw_red.max()), 3))

red digital numbers: 6187 to 55947
reflectance without the valid range: -0.03 to 1.339

# Step 3: a helper that turns digital numbers into reflectance, and blanks out invalid cells with NaN
def to_reflectance(dn):
    valid = (dn >= 7273) & (dn <= 43636)          # the producer's valid range for surface reflectance
    reflectance = dn.astype("float32") * 0.0000275 - 0.2
    return np.where(valid, reflectance, np.nan)   # NaN means "not a number": a cell we do not trust

red = to_reflectance(red_dn)
nir = to_reflectance(nir_dn)

# Step 4: check again
print("red reflectance:", round(float(np.nanmin(red)), 3), "to", round(float(np.nanmax(red)), 3))
print("nir reflectance:", round(float(np.nanmin(nir)), 3), "to", round(float(np.nanmax(nir)), 3))
print("cells outside the valid range:", int(np.isnan(red).sum()))

red reflectance: 0.006 to 0.984
nir reflectance: 0.027 to 0.872
cells outside the valid range: 11

# Step 1: read the quality band
with rasterio.open("bengaluru_landsat/QA_PIXEL.TIF") as src:
    qa = src.read(1)


# Step 2: a helper that tells us, for every cell, whether one bit is switched on
def bit_is_set(array, bit):
    return ((array >> bit) & 1) == 1       # shift right by "bit" places, then keep only the last binary digit


# Step 3: count each flag
for bit, name in [(0, "fill"), (1, "dilated cloud"), (3, "cloud"), (4, "cloud shadow"), (7, "water")]:
    print(f"{name:14s} cells: {int(bit_is_set(qa, bit).sum()):>9,}")

# Step 4: a mask for everything we do not want: fill, dilated cloud, cloud and shadow
bad = bit_is_set(qa, 0) | bit_is_set(qa, 1) | bit_is_set(qa, 3) | bit_is_set(qa, 4)
print("cells to mask:", int(bad.sum()), "of", bad.size, "=", round(float(bad.mean()) * 100, 2), "percent")

fill           cells:         0
dilated cloud  cells:        84
cloud          cells:        70
cloud shadow   cells:        83
water          cells:    14,072
cells to mask: 195 of 1661296 = 0.01 percent

# Step 1: set the flagged cells to NaN in both bands
red[bad] = np.nan
nir[bad] = np.nan

# Step 2: the NDVI formula. Cells that are NaN in either band come out as NaN.
ndvi = (nir - red) / (nir + red)

# Step 3: a summary. nanmin, nanmax and nanmean ignore the NaN cells.
print("NDVI range:", round(float(np.nanmin(ndvi)), 2), "to", round(float(np.nanmax(ndvi)), 2))
print("mean NDVI:", round(float(np.nanmean(ndvi)), 3))
print("share of valid cells with NDVI above 0.4:", round(float(np.nanmean(ndvi > 0.4)) * 100, 1), "percent")

NDVI range: -0.43 to 0.95
mean NDVI: 0.339
share of valid cells with NDVI above 0.4: 29.1 percent

# Step 1: read the corporations and check that their CRS is the same as the raster's
corp = gpd.read_file("bengaluru_landsat/corporations.gpkg")
print(corp.crs.to_epsg(), "| rows:", len(corp))
print(corp[["short", "km2", "pct"]])

# Step 2: draw NDVI with the corporation outlines on top
with rasterio.open("bengaluru_landsat/SR_B4.TIF") as src:
    left, bottom, right, top = src.bounds

fig, ax = plt.subplots(figsize=(8, 7))
image = ax.imshow(ndvi, cmap="RdYlGn", vmin=-0.2, vmax=0.8, extent=(left, right, bottom, top))
corp.boundary.plot(ax=ax, color="black", linewidth=1.2)
for _, row in corp.iterrows():
    point = row.geometry.representative_point()          # a point that is certainly inside the polygon
    ax.annotate(row["short"], (point.x, point.y), ha="center", fontsize=9, fontweight="bold")
fig.colorbar(image, ax=ax, shrink=0.75, label="NDVI")
ax.set_title("NDVI from Landsat 8, 29 March 2025")
ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)")
fig.savefig("bengaluru_ndvi.png", dpi=130, bbox_inches="tight")
print("saved")

32643 | rows: 5
     short    km2   pct
0  Central   77.9  5.63
1     East  163.7  6.67
2    North  159.8  4.75
3    South  152.3  4.73
4     West  163.6  3.20
saved

# Step 1: save NDVI as a one-band float GeoTIFF. NaN becomes the NoData value.
profile.update(dtype="float32", count=1, nodata=np.nan, compress="deflate")
with rasterio.open("ndvi_2025.tif", "w", **profile) as dst:
    dst.write(ndvi.astype("float32"), 1)

# Step 2: for each corporation, cut the NDVI raster with its shape and measure what is inside
from rasterio.mask import mask

rows = []
with rasterio.open("ndvi_2025.tif") as src:
    for _, row in corp.iterrows():
        cells, _ = mask(src, [row.geometry], crop=True, nodata=np.nan)    # cells outside the shape become NaN
        values = cells[0]
        values = values[~np.isnan(values)]                                # keep only the cells that have a value
        rows.append({
            "corporation": row["short"],
            "valid_cells": int(values.size),
            "mean_ndvi": round(float(values.mean()), 3),
            "green_share_pct": round(float((values > 0.4).mean()) * 100, 1),
        })

zonal = pd.DataFrame(rows)
print(zonal)

  corporation  valid_cells  mean_ndvi  green_share_pct
0     Central        86510      0.357             38.0
1        East       181722      0.313             23.4
2       North       177517      0.334             28.9
3       South       169201      0.321             25.7
4        West       181767      0.316             23.6

# Step 1: bring in the OpenStreetMap park and lake share, matching on the corporation's name
compare = zonal.merge(corp[["short", "pct"]], left_on="corporation", right_on="short").drop(columns="short")
compare = compare.rename(columns={"pct": "park_lake_pct"})

# Step 2: rank each corporation on both measures (1 = highest)
compare["rank_ndvi"] = compare["mean_ndvi"].rank(ascending=False).astype(int)
compare["rank_parks"] = compare["park_lake_pct"].rank(ascending=False).astype(int)
print(compare.sort_values("mean_ndvi", ascending=False).to_string(index=False))

# Step 3: how strongly do the two rankings agree? The ordinary correlation of the two rank columns
# is called the Spearman rank correlation. It runs from -1 (opposite order) to 1 (same order).
rho = compare["rank_ndvi"].corr(compare["rank_parks"])
print("Rank correlation between mean NDVI and park share:", round(float(rho), 2))

corporation  valid_cells  mean_ndvi  green_share_pct  park_lake_pct  rank_ndvi  rank_parks
    Central        86510      0.357             38.0           5.63          1           2
      North       177517      0.334             28.9           4.75          2           3
      South       169201      0.321             25.7           4.73          3           4
       West       181767      0.316             23.6           3.20          4           5
       East       181722      0.313             23.4           6.67          5           1
Rank correlation between mean NDVI and park share: 0.0

# Step 1: write the table as a CSV file that opens in Excel
compare.to_csv("ndvi_by_corporation.csv", index=False)

# Step 2: check what we saved by reading it back
print(pd.read_csv("ndvi_by_corporation.csv")[["corporation", "mean_ndvi", "rank_ndvi"]].head(3))

# Step 3: re-open the NDVI raster and confirm that its properties survived
with rasterio.open("ndvi_2025.tif") as src:
    print("ndvi_2025.tif:", src.width, "x", src.height, "|", src.dtypes[0], "| CRS", src.crs.to_epsg(), "| NoData", src.nodata)

  corporation  mean_ndvi  rank_ndvi
0     Central      0.357          1
1        East      0.313          5
2       North      0.334          2
ndvi_2025.tif: 1304 x 1274 | float32 | CRS 32643 | NoData nan