Remote Sensing Labs
Browse all docs

Make an NDVI Map in QGIS with Open HLS Data

Build a local NDVI raster in QGIS, mask the pixels flagged for poor quality, compare two fixed areas, and save the map and results. The exercise takes about 60–90 minutes after the files are downloaded. It pairs with the NDVI Monitoring workflow, which explains how one date fits into a longer monitoring workflow.

Question

For the 18 June 2024 observation, is the vegetation signal in Central Park’s Great Lawn higher than the signal over the Jacqueline Kennedy Onassis Reservoir?

Inputs

  • Area: Central Park, Manhattan, New York City. The map area is the WGS 84 bounding box (west, south, east, north) -73.982, 40.764, -73.949, 40.800.
  • Product: NASA Harmonized Landsat and Sentinel-2 (HLS) Version 2.0, Sentinel-2 component S30, MGRS tile T18TWL, 18 June 2024. The exact granule is HLS.S30.T18TWL.2024170T154941.v2.0.
  • Files: download B04.tif (red), B8A.tif (narrow near-infrared), and Fmask.tif (quality flags).
  • Sample areas: the Great Lawn rectangle is -73.9670, 40.7804, -73.9630, 40.7834; the reservoir rectangle is -73.9640, 40.7837, -73.9605, 40.7870. Coordinates are west, south, east, north in WGS 84. Check the rectangles against a basemap; keep the reservoir sample inside the shoreline.

The HLS S30 files are 30 m surface-reflectance Cloud Optimized GeoTIFFs. An Earthdata Login is required for direct downloads; creating an account is free. You can also use NASA Earthdata Search: open the HLS S30 v2.0 collection, search for the full granule ID, and download only the three files listed above. No paid data or software is needed. NASA’s data access guide lists other access routes. NASA asks users to cite its data and material; see its data use policy.

HLS stores surface reflectance as signed 16-bit values with a scale factor of 0.0001 and fill value -9999. The Fmask layer is an unsigned byte and uses 255 as fill. In S30, B04 is red and B8A is narrow near-infrared. The QA bits mark cloud (bit 1), cloud/shadow adjacency (bit 2), cloud shadow (bit 3), snow/ice (bit 4), water (bit 5), and aerosol level (bits 6–7). HLS reserves bit 0. HLS dilates cloud and shadow flags by five pixels, or 150 m, and labels the added area as adjacent to cloud/shadow. The HLS Version 2.0 User Guide defines the scaling and quality bits.

Prepare the project

  1. Make a folder named central-park-ndvi with input and output subfolders. Put the three downloaded TIFF files in input.
  2. Open QGIS 4.2 or later. Drag the three TIFFs onto the map. In the Layers panel, rename them red, nir, and qa so they are easy to assign to calculator inputs A, B, and C. Save the project as central-park-ndvi/central-park-ndvi.qgz.
  3. Open each raster’s Properties → Information. Confirm that each is 3,660 columns by 3,660 rows, has 30 m pixels and uses EPSG:32618. The rasters must share the same extent and grid. Each input has 13,395,600 pixels (3,660 × 3,660); the red and NIR fill value should be -9999, and the QA fill value should be 255. These checks catch an incorrect band, date, or product version before processing.

Calculate NDVI and mask low-quality pixels

In the Processing Toolbox, open GDAL → Raster miscellaneous → Raster calculator (gdal:rastercalculator). Set input A to red, B to nir, and C to qa, using band 1 for each. Set output NoData to -9999, extent handling to Fail, and output type to Float32. Save the result as output/central-park-ndvi.tif. The earlier grid check matters here: GDAL’s calculator requires equal dimensions and does not check that projections match. Paste this expression:

numpy.where(
  logical_and(
    logical_and(C != 255, (C & 30) == 0),
    logical_and(
      (C & 192) != 192,
      logical_and(
        logical_and(A != -9999, B != -9999),
        logical_and(
          logical_and(A != 12000, B != 12000),
          (A.astype(numpy.float64) + B.astype(numpy.float64)) != 0
        )
      )
    )
  ),
  (B.astype(numpy.float64) - A.astype(numpy.float64)) /
    numpy.where(
      (B.astype(numpy.float64) + A.astype(numpy.float64)) != 0,
      B.astype(numpy.float64) + A.astype(numpy.float64),
      1
    ),
  -9999
)

The expression removes cloud, adjacency, shadow, snow/ice, and high-aerosol flags while keeping water pixels. (C & 30) == 0 checks bits 1–4; (C & 192) != 192 excludes aerosol code 11 in bits 6–7. HLS bit 0 is reserved and unused; lower quality flags can be combined. Reflectance fill, QA fill, per-band saturation value 12000, and zero-denominator pixels become NoData. NASA documents 12000 as the S30 saturation flag on the HLS S30 product page. This fixed granule contains no 12000 values in B04 or B8A, checked with QGIS 4.2.3, so the extra mask leaves the reported counts and medians unchanged. The scale factor 0.0001 is identical for both reflectance bands, so it cancels in the normalized ratio; the expression uses their stored values and returns the same NDVI as scaling both bands first. QGIS documents this GDAL calculator as NumPy array algebra and lets you set the output NoData value and Float32 type in the same dialog. See the QGIS GDAL Raster calculator reference for its inputs and parameters.

Compare the Great Lawn and reservoir

Create a small polygon layer with the two fixed rectangles. In a plain-text editor, save the following as input/sample-areas.geojson:

{
  "type": "FeatureCollection",
  "name": "sample_areas",
  "features": [
    {
      "type": "Feature",
      "properties": {"site": "Great Lawn"},
      "geometry": {"type": "Polygon", "coordinates": [[
        [-73.9670, 40.7804], [-73.9630, 40.7804],
        [-73.9630, 40.7834], [-73.9670, 40.7834],
        [-73.9670, 40.7804]
      ]]}
    },
    {
      "type": "Feature",
      "properties": {"site": "Reservoir"},
      "geometry": {"type": "Polygon", "coordinates": [[
        [-73.9640, 40.7837], [-73.9605, 40.7837],
        [-73.9605, 40.7870], [-73.9640, 40.7870],
        [-73.9640, 40.7837]
      ]]}
    }
  ]
}

Add the GeoJSON with Layer → Add Layer → Add Vector Layer. In the Processing Toolbox, run Zonal statistics with sample-areas as the input polygon layer, central-park-ndvi as the raster, band 1, and prefix ndvi_. Select Count and Median, then save the output as output/ndvi-samples.gpkg. The output table should contain ndvi_count and ndvi_median fields. The QGIS Zonal statistics reference describes these settings.

Style and save the results

  1. In the NDVI layer’s Properties → Symbology, choose Singleband pseudocolor. Set the minimum to -1, the maximum to 1, and use a red–yellow–green ramp so higher values are green. Keep this range fixed for comparison.
  2. Style the sample polygons with transparent fill and contrasting outlines. Zoom to the WGS 84 map area above and check that the reservoir rectangle is over open water and the lawn rectangle is over the Great Lawn.
  3. Save the project. Open Project → New Print Layout, add a map, title, legend, scale bar, and source note (NASA HLS S30, 2024-06-18). Export the layout as output/central-park-ndvi.pdf.
  4. Right-click the output GeoPackage in the Layers panel and choose Export → Save Features As. Set Format to CSV and File name to output/ndvi-samples.csv. Under Select fields to export, keep only site, ndvi_count, and ndvi_median. Set Geometry to No geometry, then save. This gives you a two-row table without a WKT column.

Expected result and checks

The Great Lawn median should be much higher than the reservoir median. For the fixed sample rectangles and this mask, the expected results are:

SampleValid pixelsMedian NDVI
Great Lawn1010.7441
Reservoir1000.1403

These expected values were reproduced with QGIS 4.2.3’s GDAL raster calculator and Zonal statistics tools, and independently checked by reading the same COGs with Rasterio. Zonal statistics counts pixels whose centers fall inside each rectangle. Small differences in the last decimal can result from output precision. If your count differs, check the polygon coordinates, grid alignment, and mask expression before interpreting the map. The expected scene grid is 3,660 × 3,660, or 13,395,600 pixels. The masked raster has 10,639,186 valid pixels; values range from about -0.8974 to 0.9882, with none outside [-1, 1].

The map should show a strong contrast between green vegetation and open water, with mixed shoreline, trees, shadows, and atmosphere contributing local variation. One summer image describes the surface’s spectral greenness on one date; it does not establish plant health, biomass, or a trend. For those questions, see NDVI Monitoring and Cloud Masking.

Sources