Remote Sensing Labs
Browse all docs

Reproduce the Central Park HLS NDVI Exercise in Python

This Python companion uses the same scene and two sample rectangles as the QGIS HLS exercise. It follows one date from reflectance bands through a masked NDVI raster to two reproducible summaries. For the index’s meaning and limits, start with spectral bands and indices and NDVI Monitoring.

The result describes spectral greenness on one date. It does not measure plant health, biomass, or a trend.

What you need

  • Python 3.11 or later.
  • The three HLS GeoTIFFs listed below. Direct downloads from NASA LP DAAC require a free Earthdata Login. The files and software are free; no Google account, Earth Engine project, or paid GIS software is needed.
  • About 150 MB of disk space for inputs and outputs. The script loads the 3,660 × 3,660 bands into arrays, so allow several hundred megabytes of available memory.

Make a working folder and a small virtual environment. These commands use macOS/Linux:

mkdir -p central-park-python-ndvi/input central-park-python-ndvi/output
cd central-park-python-ndvi
python3 -m venv .venv
source .venv/bin/activate
python -m pip install numpy rasterio

In Windows PowerShell, use:

New-Item -ItemType Directory -Force central-park-python-ndvi
Set-Location central-park-python-ndvi
New-Item -ItemType Directory -Force input, output
py -3.11 -m venv .venv
./.venv/Scripts/python.exe -m pip install numpy rasterio

Download the fixed scene

Use NASA Earthdata Search to find the HLS S30 Version 2.0 granule HLS.S30.T18TWL.2024170T154941.v2.0, observed on 18 June 2024. Select only:

Save them in input/ without renaming. The red and near-infrared rasters should be signed 16-bit reflectance with fill -9999; Fmask should be an unsigned byte with fill 255. All three should have 3,660 rows and columns, 30 m pixels, EPSG:32618, and the same grid. If one differs, stop and check the granule and band before continuing.

HLS reflectance uses a scale factor of 0.0001; the same factor applies to both bands and cancels in their normalized ratio. The S30 product uses 12000 as a saturation flag. Fmask bits 1–4 flag cloud, adjacency, shadow, and snow/ice; bit 5 flags water; bits 6–7 encode aerosol level. This workflow excludes cloud, adjacency, shadow, snow/ice, fill, saturation, and the highest aerosol code. It keeps water pixels so the reservoir remains a comparison sample. NASA’s current S30 product page and HLS Version 2.0 User Guide define the band values and quality bits.

Create the analysis script

Save the following as analyze.py in central-park-python-ndvi/. It checks that the inputs share one grid, writes a Float32 NDVI GeoTIFF, and records input hashes, processing choices, software versions, and sample results in output/run.json.

from __future__ import annotations

import hashlib
import json
import platform
from pathlib import Path

import numpy as np
import rasterio
from rasterio.features import geometry_mask
from rasterio.warp import transform_geom


ROOT = Path(__file__).resolve().parent
INPUT = ROOT / "input"
OUTPUT = ROOT / "output"
GRANULE = "HLS.S30.T18TWL.2024170T154941.v2.0"
FILES = {
    "red": INPUT / f"{GRANULE}.B04.tif",
    "nir": INPUT / f"{GRANULE}.B8A.tif",
    "fmask": INPUT / f"{GRANULE}.Fmask.tif",
}

SAMPLES = {
    "Great Lawn": {
        "type": "Polygon",
        "coordinates": [[
            [-73.9670, 40.7804], [-73.9630, 40.7804],
            [-73.9630, 40.7834], [-73.9670, 40.7834],
            [-73.9670, 40.7804],
        ]],
    },
    "Reservoir": {
        "type": "Polygon",
        "coordinates": [[
            [-73.9640, 40.7837], [-73.9605, 40.7837],
            [-73.9605, 40.7870], [-73.9640, 40.7870],
            [-73.9640, 40.7837],
        ]],
    },
}

EXPECTED = {
    "Great Lawn": {"valid_pixels": 101, "median_ndvi": 0.7441},
    "Reservoir": {"valid_pixels": 100, "median_ndvi": 0.1403},
}


def sha256(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as source:
        for chunk in iter(lambda: source.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def read_band(path: Path) -> tuple[np.ndarray, dict, dict]:
    if not path.is_file():
        raise FileNotFoundError(f"Missing input: {path}")
    with rasterio.open(path) as src:
        if src.count != 1:
            raise ValueError(f"Expected one band in {path.name}; found {src.count}")
        array = src.read(1)
        profile = src.profile.copy()
        info = {
            "file": path.name,
            "bytes": path.stat().st_size,
            "sha256": sha256(path),
            "dtype": src.dtypes[0],
            "nodata": src.nodata,
            "width": src.width,
            "height": src.height,
            "crs": str(src.crs),
            "transform": list(src.transform),
        }
    return array, profile, info


red, profile, red_info = read_band(FILES["red"])
nir, nir_profile, nir_info = read_band(FILES["nir"])
fmask, qa_profile, qa_info = read_band(FILES["fmask"])

grid = (profile["width"], profile["height"], profile["crs"], profile["transform"])
for name, other in (("B8A", nir_profile), ("Fmask", qa_profile)):
    other_grid = (other["width"], other["height"], other["crs"], other["transform"])
    if other_grid != grid:
        raise ValueError(f"{name} does not match the B04 grid")
if red.shape != (3660, 3660) or red.dtype != np.int16 or nir.dtype != np.int16:
    raise ValueError("Inputs do not match the expected HLS S30 grid and reflectance type")
if fmask.dtype != np.uint8 or profile["crs"].to_epsg() != 32618:
    raise ValueError("Inputs do not match the expected HLS Fmask type or CRS")
if profile["nodata"] != -9999 or nir_profile["nodata"] != -9999 or qa_profile["nodata"] != 255:
    raise ValueError("Input fill values do not match the HLS S30 files")

red64 = red.astype(np.float64)
nir64 = nir.astype(np.float64)
denominator = red64 + nir64
valid = (
    (fmask != 255)
    & ((fmask & 30) == 0)          # cloud, adjacency, shadow, snow/ice
    & ((fmask & 192) != 192)       # aerosol code 3
    & (red != -9999)
    & (nir != -9999)
    & (red != 12000)
    & (nir != 12000)
    & (denominator != 0)
)

ndvi = np.full(red.shape, -9999, dtype=np.float32)
ndvi[valid] = ((nir64[valid] - red64[valid]) / denominator[valid]).astype(np.float32)

OUTPUT.mkdir(exist_ok=True)
ndvi_path = OUTPUT / "central-park-ndvi.tif"
output_profile = profile.copy()
output_profile.update(
    driver="GTiff",
    count=1,
    dtype="float32",
    nodata=-9999.0,
    compress="deflate",
    predictor=3,
)
with rasterio.open(ndvi_path, "w", **output_profile) as dst:
    dst.write(ndvi, 1)

results = {}
for name, polygon_wgs84 in SAMPLES.items():
    polygon = transform_geom("EPSG:4326", profile["crs"], polygon_wgs84)
    inside = geometry_mask(
        [polygon],
        out_shape=red.shape,
        transform=profile["transform"],
        all_touched=False,
        invert=True,
    )
    sample = ndvi[inside & valid]
    median = float(np.median(sample))
    row = {"valid_pixels": int(sample.size), "median_ndvi": round(median, 4)}
    if row != EXPECTED[name]:
        raise ValueError(f"{name} does not match the reference: {row}")
    results[name] = {
        **row,
        "geometry_wgs84": polygon_wgs84,
        "pixel_rule": "pixel centers inside plus Bresenham-selected polygon boundary pixels (all_touched=false)",
    }
    print(f"{name}: {row['valid_pixels']} valid pixels; median NDVI {row['median_ndvi']:.4f}")

report = {
    "product": "NASA HLS S30 Version 2.0",
    "granule": GRANULE,
    "date": "2024-06-18",
    "bands": {"red": "B04", "narrow_near_infrared": "B8A", "quality": "Fmask"},
    "reflectance_scale_factor": 0.0001,
    "ndvi": "(B8A - B04) / (B8A + B04); common scale factor cancels",
    "mask": {
        "fmask_fill": 255,
        "exclude_fmask_bits_1_to_4": True,
        "keep_water_bit_5": True,
        "exclude_aerosol_code_3": True,
        "reflectance_fill": -9999,
        "saturation_value": 12000,
        "zero_denominator": "excluded",
    },
    "grid": {
        "crs": str(profile["crs"]),
        "width": profile["width"],
        "height": profile["height"],
        "transform": list(profile["transform"]),
        "valid_pixels": int(valid.sum()),
        "ndvi_range": [float(ndvi[valid].min()), float(ndvi[valid].max())],
    },
    "inputs": {"B04": red_info, "B8A": nir_info, "Fmask": qa_info},
    "software": {
        "python": platform.python_version(),
        "numpy": np.__version__,
        "rasterio": rasterio.__version__,
        "gdal": rasterio.__gdal_version__,
    },
    "samples": results,
    "outputs": [ndvi_path.name, "run.json"],
}
(OUTPUT / "run.json").write_text(json.dumps(report, indent=2) + "\n")
print(f"Wrote {ndvi_path} and {OUTPUT / 'run.json'}")

Run it from the working folder:

python analyze.py

In Windows PowerShell, use the environment’s Python executable:

./.venv/Scripts/python.exe analyze.py

Expected output:

Great Lawn: 101 valid pixels; median NDVI 0.7441
Reservoir: 100 valid pixels; median NDVI 0.1403

geometry_mask uses Rasterio’s default selection when all_touched=False: pixels selected by their centers and by the Bresenham line algorithm along polygon boundaries. For these fixed rectangles and this grid, its counts and medians reproduce the QGIS reference. If you change the polygons or grid, recheck edge-pixel counts against the zonal-statistics method you choose. The script fails if the grids, fill values, counts, or four-decimal medians differ; check the exact granule, mask, and sample coordinates before changing an expected value. For a visual check of the lawn and reservoir rectangles, follow the QGIS exercise’s map step.

Interpret the result

The Great Lawn sample has the higher median NDVI in this scene. That is a single-date spectral comparison, not evidence of a long-term vegetation change or a measure of plant condition. Review NASA’s data use and citation guidance when sharing results, and name the product, granule, date, bands, mask, and sample areas with the values.

References