What is NDVI? Formula, Values, and Uses

NDVI, the normalized difference vegetation index, measures how much healthy green vegetation a pixel contains, calculated from the red and near-infrared light it reflects. Values run from −1 to +1. Dense, healthy vegetation scores high, bare ground scores near zero, and water scores below zero. NDVI is computed from Sentinel-2, Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS), drone, and aircraft imagery. This article tests one claim. A single NDVI date cannot identify a crop. NDVI supports a field decision only once it is summarized per field across a season.

The Havasu catalog holds a Sentinel-2 scene index of L2A (Level-2A, corrected for the atmosphere) imagery, seasonal mosaics at 10 m, and Fields of the World field boundaries.

Key takeaways

• NDVI compares NIR (near-infrared light) and red reflectance: NDVI = (NIR − Red) / (NIR + Red), computed per pixel on surface reflectance.
• Sentinel-2 NDVI uses bands B8 and B4 at 10 m. Landsat 8 and 9 NDVI uses bands B5 and B4 at 30 m.
• Values near +1 mean dense green vegetation, values near 0 mean bare ground, and negative values mean water.
• NDVI stops rising over dense canopy and shifts with soil brightness and haze. The enhanced vegetation index (EVI), soil-adjusted vegetation index (SAVI), and normalized difference red edge index (NDRE) correct for these cases.

The NDVI formula

The NDVI formula divides the difference between near-infrared (NIR) and red reflectance by their sum:

NDVI = (NIR − Red) / (NIR + Red)

Chlorophyll absorbs red light for photosynthesis. The cell structure of leaves reflects near-infrared light strongly. NASA's Earth Observatory explains NDVI through this contrast. Dividing by the sum keeps the index between −1 and +1 at any scene brightness.

Each sensor records red and near-infrared in different numbered bands:

Sensor	Red band	NIR band	Pixel size	NDVI formula
Sentinel-2 MultiSpectral Instrument (MSI)	B4 (665 nm)	B8 (842 nm)	10 m	(B8 − B4) / (B8 + B4)
Landsat 8 and 9 Operational Land Imager (OLI)	B4 (0.64 to 0.67 µm)	B5 (0.85 to 0.88 µm)	30 m	(B5 − B4) / (B5 + B4)
Landsat 4 to 7 Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+)	B3	B4	30 m	(B4 − B3) / (B4 + B3)

The U.S. Geological Survey (USGS) publishes the Landsat band formulas. The European Space Agency (ESA) documents the Sentinel-2 bands. Sentinel-2 also records a narrow NIR band, B8A at 20 m, which some workflows use in place of B8.

What NDVI values mean

NDVI values describe land cover along a gradient from water to dense forest. NASA's Earth Observatory gives these ranges, and notes that values of 0.8 to 0.9 indicate the highest possible density of green leaves:

NDVI value	Typical surface
Below 0	Water
0.1 and below	Barren rock, sand, or snow
0.2 to 0.3	Shrubs and grassland
0.6 to 0.8	Temperate and tropical rainforest

A healthy wheat field at peak growth can score as high as a forest, and the same field reads near zero after harvest. Sensor, season, and processing level shift the numbers. NDVI values compare best against the same place at another time, or against neighbors imaged on the same date.

NDVI vs EVI, NDRE, NDWI, and SAVI

Each related index uses different bands or adds terms for a specific weakness:

Index	Formula	Bands (Sentinel-2)	Best for
NDVI	(NIR − Red) / (NIR + Red)	B8, B4	General vegetation greenness and cover
EVI	2.5 × (NIR − Red) / (NIR + 6 × Red − 7.5 × Blue + 1)	B8, B4, B2	Dense canopy, hazy scenes
NDRE	(NIR − Red Edge) / (NIR + Red Edge)	B8A or B8, B5	Chlorophyll and nitrogen in mid- to late-season crops
NDWI	(Green − NIR) / (Green + NIR)	B3, B8	Open water mapping
SAVI	(NIR − Red) / (NIR + Red + L) × (1 + L), L = 0.5	B8, B4	Sparse vegetation over bright soil

Alfredo Huete introduced the soil-adjusted vegetation index in 1988. Liu and Huete corrected for atmosphere and canopy background together in a feedback-based modification of the NDVI (1995), the precursor of EVI. Huete and colleagues documented EVI as the MODIS companion to NDVI in 2002. The USGS publishes the Landsat versions of EVI and SAVI. The normalized difference water index (NDWI) here is McFeeters' water index (1996). A second index with the same name, Gao's NDWI from the same year, uses NIR and shortwave infrared to measure leaf water content. The USGS calls its Landsat version the normalized difference moisture index (NDMI). Hyperspectral imagery supports hundreds of narrow-band indices beyond these five.

Saturation over dense canopy. Once leaves cover the ground, red reflectance drops close to zero. Further leaf growth then barely changes the ratio. NDVI flattens near its upper range over forests and mature crops. The red and near-infrared plot above shows where saturation starts in the Haskell County scene. Pixels with NDVI above 0.8 cluster at a median red reflectance of 0.030 against a near-infrared of 0.421. EVI and NDRE stay sensitive there. Gitelson, Kaufman, and Merzlyak showed that a green band in place of red stays sensitive to chlorophyll (1996). Gitelson's Wide Dynamic Range Vegetation Index (2004) weights NIR by 0.1 to 0.2. It is at least three times more sensitive than NDVI at leaf area indices from 2 to 6 (square meters of leaf per square meter of ground). Mutanga and Skidmore showed that narrow-band indices overcome the saturation problem in biomass estimation (2004).

Soil background. In sparse vegetation, bright, wet, and dark soils shift NDVI by different amounts at the same plant cover. SAVI's soil factor corrects for this.

Atmosphere. Haze, smoke, and thin cloud scatter red light more than NIR and lower NDVI. Surface reflectance correction reduces the effect. EVI's blue term corrects for haze particles that remain.

How NDVI is calculated from satellite imagery

To calculate NDVI from satellite imagery:

1. Start from surface reflectance. Sentinel-2 Level-2A and Landsat Collection 2 Level-2 products are corrected for the atmosphere. Both store reflectance as scaled integers. Check whether the file or its SpatioTemporal Asset Catalog (STAC) record carries a scale and offset before computing the index. The USGS lists the Landsat factors in its Collection 2 documentation. A scale factor cancels in the ratio. An offset does not cancel. Subtracting an offset the file does not need, or skipping one it needs, shifts every NDVI value, and no error is raised. The RS_FromPath documentation spells out when WherobotsDB applies the scale and offset stored in a file. The Sentinel-2 Cloud Optimized GeoTIFF (COG) files in the Havasu scene index store reflectance as integers with no scale or offset tags. A COG is a tiled image file a reader can fetch in parts. The queries below use the stored values directly.
2. Mask clouds, shadows, and snow. Clouds score near zero or below and would read as vegetation loss. Sentinel-2 L2A ships a Scene Classification Layer (SCL) that labels cloud, cloud shadow, cirrus, snow, and water per pixel. Landsat Level-2 ships a QA_PIXEL band with the same flags as bits.
3. Align the bands. Red and NIR must share one grid. On Sentinel-2, B4 and B8 are both 10 m, and the 20 m SCL is resampled to match.
4. Apply the formula pixel by pixel. The result is single-band NDVI imagery, usually saved as a Cloud Optimized GeoTIFF and colored from red to green as an NDVI map.

Sentinel-2 revisits each location every five days at the equator with two satellites (Drusch and colleagues, 2012). Landsat extends the record back to 1972. Both are open Earth observation data.

Stacking every cloud-free observation of a pixel through a year traces a curve. NDVI rises at green-up, peaks at full canopy, and declines as leaves age and dry (senescence) or the crop is harvested. Analysts read the start, peak, and end of season from that curve. The difference from the same date in prior years is an NDVI anomaly. They fill cloud gaps with composites such as NASA's 16-day, 250 m MODIS vegetation index product (MOD13Q1). Where clouds persist for weeks, synthetic aperture radar images through them. NDVI differencing between two dates is one of the simplest forms of change detection.

NDVI examples

• Great Plains rangeland, 1973 to 1974. Rouse's team used the index to follow the spring "green wave" of rangeland along a Great Plains corridor. They reported it to NASA in Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation.
• Africa's land cover and the 1983 Sahel drought. Tucker, Townshend, and Goff classified African land cover from 19 months of NDVI from the Advanced Very High Resolution Radiometer (AVHRR) on National Oceanic and Atmospheric Administration (NOAA) satellites. They followed the drought as it happened (Science, 1985).
• Northern greening. Myneni and colleagues found increased plant growth between 45°N and 70°N from 1981 to 1991 in the AVHRR NDVI record, tied to earlier spring snowmelt (Nature, 1997).
• Kansas winter wheat yield. Becker-Reshef and colleagues built a MODIS NDVI regression model for winter wheat yields in Kansas and Ukraine (2010).
• Drought monitoring. Kogan's Vegetation Condition Index rescales each pixel's NDVI against its own historical range (1995). Brown and colleagues combined NDVI with climate data in the Vegetation Drought Response Index (2008).
• Wildfire fuel and urban green space. Fire analysts track how cured grass fuels are through the season, and planners measure tree canopy and park cover per neighborhood.

NDVI per field with zonal statistics over field boundaries

Zonal statistics return the mean, standard deviation, or count of the NDVI pixels inside each polygon, such as a field boundary. A field whose mean NDVI falls below its neighbors, or below its own prior-year curve, is a candidate for crop stress. A high standard deviation within a field points to uneven growth from drainage, soil, or irrigation problems.

Parcel data ties NDVI to a legal unit of land and its owner. Field boundary datasets follow the planted area, which a parcel can split or merge. The choice of polygon changes the result. A polygon drawn by a model can hold two management units, parts of a field that are farmed differently. Mixed 10 m pixels along a road or a neighbor's crop also shift the mean. Wherobots used this pattern for mean AlphaEarth embeddings over about 3.5 million Iowa fields across 2023 to 2025, and NDVI works the same way.

NDVI in Wherobots

The example below uses one cloud-free Sentinel-2 scene of Haskell County, Kansas, from 28 July 2024 (S2A_14SLG_20240728, 0.55% cloud). It joins that scene to the 2024 Fields of the World boundaries for a 12 km by 11 km block of center-pivot and dryland farms west of Sublette.

Source	Format	What the queries read
Sentinel-2 L2A scene index (wherobots_open_data.sentinel2.l2a_source_items)	Table of STAC items, one row per scene, with links to each band file	1 row, the 28 July scene, for the single-date query; 15 rows for the season query
Sentinel-2 bands B04 (red), B08 (NIR), and SCL (scene classification)	Cloud Optimized GeoTIFFs at 10 m, 10 m, and 20 m	The tiles inside the block, about 1.3 million 10 m pixels; 3 files per date, 45 files for the season
Fields of the World boundaries (wherobots_open_data.rasterflow_output_samples.fields_of_the_world_vector_global)	Table of field polygons, with a date for each boundary vintage	The 2024 polygons of 2 ha or more inside the block: 102 fields, 7,101 ha

WherobotsDB reads the raster files and the polygon table in one SQL query: raster and vector in one engine; most stacks pick one. RS_MapAlgebra applies the formula and the cloud mask in one pass. RS_ZonalStats returns the mean, standard deviation, and count of clear NDVI pixels inside each polygon. The scene index and the STAC reader find scenes by area, date, and cloud cover. Each band is an out-db raster: the table stores a link, and WherobotsDB reads the pixels from the COG. RasterFlow ran the Fields of the World model that produced the boundaries. The model is based on the Fields of The World benchmark by Kerner and colleagues (AAAI 2025).

-- Mean Sentinel-2 NDVI per Fields of the World polygon, Haskell County, Kansas, 28 July 2024
WITH scene AS (
  SELECT assets.red.href AS red_href, assets.nir.href AS nir_href, assets.scl.href AS scl_href
  FROM wherobots_open_data.sentinel2.l2a_source_items
  WHERE year = 2024 AND month = 7 AND mgrs_tile = '14SLG' AND id = 'S2A_14SLG_20240728_0_L2A'
),
block AS (
  -- The 12 km by 11 km block, in the scene's UTM zone 14N
  SELECT ST_Transform(ST_SetSRID(ST_PolygonFromEnvelope(-100.98, 37.48, -100.86, 37.58), 4326),
                      'EPSG:4326', 'EPSG:32614') AS box
),
ndvi AS (
  -- One NDVI raster for the whole block. Stack red (B04), NIR (B08) and the 20 m scene
  -- classification layer resampled to 10 m; keep SCL 4 (vegetation) and 5 (not vegetated)
  SELECT RS_MapAlgebra(
           RS_AddBand(RS_AddBand(red, nir, 1, 2), RS_ReprojectMatch(scl, red, 'NearestNeighbor'), 1, 3),
           'D',
           'out = (rast[0] > 0 && rast[1] > 0 && (rast[2] == 4 || rast[2] == 5))
                  ? (rast[1] - rast[0]) / (rast[1] + rast[0]) : -9999.0;', -9999.0) AS rast
  FROM (
    -- RS_Clip(raster, band, geometry, allTouched, noDataValue, crop): reads only the COG tiles in the block
    SELECT RS_Clip(RS_FromPath(s.red_href), 1, k.box, false, 0, true) AS red,
           RS_Clip(RS_FromPath(s.nir_href), 1, k.box, false, 0, true) AS nir,
           RS_Clip(RS_FromPath(s.scl_href), 1, k.box, false, 0, true) AS scl
    FROM scene s CROSS JOIN block k
  )
),
fields AS (
  SELECT ST_AreaSpheroid(geometry) AS area_m2,
         ST_Transform(ST_SetSRID(geometry, 4326), 'EPSG:4326', 'EPSG:32614') AS geom
  FROM wherobots_open_data.rasterflow_output_samples.fields_of_the_world_vector_global
  WHERE bbox.xmin > -100.98 AND bbox.xmax < -100.86 AND bbox.ymin > 37.48 AND bbox.ymax < 37.58
    AND time < TIMESTAMP '2024-06-01 00:00:00'   -- the 2024 boundary vintage
    AND ST_AreaSpheroid(geometry) >= 20000       -- fields of 2 ha or more
)
-- Send the single NDVI raster to every field partition, then summarize each field
SELECT /*+ BROADCAST(n) */ ROUND(f.area_m2) AS area_m2,
       RS_ZonalStats(n.rast, f.geom, 1, 'mean')   AS mean_ndvi,
       RS_ZonalStats(n.rast, f.geom, 1, 'stddev') AS sd_ndvi,
       RS_ZonalStats(n.rast, f.geom, 1, 'count')  AS clear_pixels
FROM fields f CROSS JOIN ndvi n;

The map algebra runs once over about 1.3 million pixels. The broadcast hint copies that raster to every worker that holds part of the field table. Raster processing at scale: the out-of-database architecture behind WherobotsDB explains why zonal statistics read only the pixels inside each zone.

The query ran in about three minutes. On 28 July, mean NDVI per field ranged from 0.13 to 0.89, with a median of 0.44. 24 fields averaged above 0.7, most of them summer crops at full canopy. 30 averaged below 0.3, including harvested wheat and fallow ground. The median within-field standard deviation was 0.07. Four fields had a standard deviation above 0.3, a spread that points to green canopy and bare ground inside one polygon.

A full 2024 season for the same fields

The season query runs on all 15 cloud-free Sentinel-2 L2A scenes of tile 14SLG in 2024, with the same SCL mask and formula. For each field it takes the highest April or May mean (the spring peak) and the highest July or August mean (the summer peak). A field is spring peak when spring exceeds 0.6 and summer stays under 0.45. It is summer peak when summer exceeds 0.7 and spring stays under 0.45. It is low until November when both stay under 0.35. Other fields are mixed. The thresholds are set for this block. It returns the median NDVI per pattern and date and ran in under four minutes.

-- Mean NDVI per field on 15 clear Sentinel-2 dates in 2024, grouped by seasonal pattern
WITH scenes AS (
  -- REPARTITION spreads the 15 scenes over 15 tasks, so each date's reads, map algebra and
  -- zonal statistics run in parallel
  SELECT /*+ REPARTITION(15) */ id, to_date(datetime) AS scene_date,
         assets.red.href AS red_href, assets.nir.href AS nir_href, assets.scl.href AS scl_href
  FROM wherobots_open_data.sentinel2.l2a_source_items
  WHERE year = 2024 AND mgrs_tile = '14SLG'
    AND id IN ('S2A_14SLG_20240130_0_L2A','S2B_14SLG_20240224_0_L2A','S2B_14SLG_20240322_0_L2A',
               'S2B_14SLG_20240414_0_L2A','S2B_14SLG_20240514_0_L2A','S2A_14SLG_20240605_0_L2A',
               'S2A_14SLG_20240628_0_L2A','S2B_14SLG_20240713_0_L2A','S2A_14SLG_20240728_0_L2A',
               'S2A_14SLG_20240817_0_L2A','S2A_14SLG_20240903_0_L2A','S2A_14SLG_20240926_0_L2A',
               'S2A_14SLG_20241016_0_L2A','S2B_14SLG_20241110_0_L2A','S2B_14SLG_20241207_0_L2A')
),
block AS (
  SELECT ST_Transform(ST_SetSRID(ST_PolygonFromEnvelope(-100.98, 37.48, -100.86, 37.58), 4326),
                      'EPSG:4326', 'EPSG:32614') AS box
),
ndvi AS (
  -- One NDVI raster per date for the block
  SELECT s.scene_date, RS_MapAlgebra(
           RS_AddBand(RS_AddBand(r, n, 1, 2), RS_ReprojectMatch(l, r, 'NearestNeighbor'), 1, 3),
           'D',
           'out = (rast[0] > 0 && rast[1] > 0 && (rast[2] == 4 || rast[2] == 5))
                  ? (rast[1] - rast[0]) / (rast[1] + rast[0]) : -9999.0;', -9999.0) AS rast
  FROM (SELECT /*+ BROADCAST(k) */ s.scene_date,
               RS_Clip(RS_FromPath(s.red_href), 1, k.box, false, 0, true) AS r,
               RS_Clip(RS_FromPath(s.nir_href), 1, k.box, false, 0, true) AS n,
               RS_Clip(RS_FromPath(s.scl_href), 1, k.box, false, 0, true) AS l
        FROM scenes s CROSS JOIN block k) s
),
fields AS (
  -- The table has no ID column, so a 10-character geohash of each field's centroid serves as one
  SELECT ST_GeoHash(ST_Centroid(geometry), 10) AS field_id,
         ST_Transform(ST_SetSRID(geometry, 4326), 'EPSG:4326', 'EPSG:32614') AS geom
  FROM wherobots_open_data.rasterflow_output_samples.fields_of_the_world_vector_global
  WHERE bbox.xmin > -100.98 AND bbox.xmax < -100.86 AND bbox.ymin > 37.48 AND bbox.ymax < 37.58
    AND time < TIMESTAMP '2024-06-01 00:00:00'
    AND ST_AreaSpheroid(geometry) >= 20000
),
per_field AS (
  -- One row per field and date. The 102 field polygons are small, so broadcast them
  SELECT /*+ BROADCAST(f) */ f.field_id, n.scene_date,
         RS_ZonalStats(n.rast, f.geom, 1, 'mean')  AS mean_ndvi,
         RS_ZonalStats(n.rast, f.geom, 1, 'count') AS clear_pixels
  FROM ndvi n CROSS JOIN fields f
),
peaks AS (
  -- Each field's highest mean NDVI in spring (April, May) and summer (July, August)
  SELECT *,
         MAX(CASE WHEN month(scene_date) IN (4, 5) THEN mean_ndvi END) OVER (PARTITION BY field_id) AS spring_peak,
         MAX(CASE WHEN month(scene_date) IN (7, 8) THEN mean_ndvi END) OVER (PARTITION BY field_id) AS summer_peak
  FROM per_field
),
patterns AS (
  SELECT *,
         CASE WHEN spring_peak > 0.6 AND summer_peak < 0.45 THEN 'spring peak'
              WHEN summer_peak > 0.7 AND spring_peak < 0.45 THEN 'summer peak'
              WHEN spring_peak < 0.35 AND summer_peak < 0.35 THEN 'low until November'
              ELSE 'mixed' END AS pattern
  FROM peaks
)
-- One row per date and pattern, plus an all-fields row per date (pattern is NULL).
-- SELECT * FROM per_field instead returns each field's own curve.
SELECT scene_date, pattern, COUNT(*) AS fields,
       ROUND(AVG(mean_ndvi), 3) AS avg_ndvi, ROUND(percentile(mean_ndvi, 0.5), 3) AS median_ndvi,
       SUM(clear_pixels) AS clear_px
FROM patterns
GROUP BY GROUPING SETS ((scene_date, pattern), (scene_date))
ORDER BY scene_date, pattern;

All 102 fields have a value on every date. The average field mean was 0.255 on 30 January, rose to 0.384 on 14 May, and peaked at 0.479 on 13 July before falling to 0.266 by 16 October. The 28 July row reproduces the single-date result, a median field mean of 0.44.

The pattern rows split the averages into spring-peak, summer-peak, and low-until-November fields:

• Spring peak, 13 fields. Median NDVI is 0.47 in late January, 0.82 on 14 April, and 0.16 by 28 June. That matches the winter wheat calendar of southwest Kansas: sown in autumn, green through spring, harvested in early summer.
• Summer peak, 20 fields. Median NDVI stays near 0.16 through May, climbs to 0.57 on 28 June, peaks at 0.77 on 28 July, and drops to 0.22 by 16 October.
• Low until November, 14 fields. Median NDVI stays between 0.15 and 0.26 until October, then rises to 0.45 on 10 November and holds there in December. That fits fallow ground through the summer and a fall-sown crop.

The other 55 fields show double peaks, partial-field patterns, or values that fall between the thresholds. Mixed curves like these can mean a field polygon holds more than one management unit.

To cover a county or a state, the query repeats per Sentinel-2 tile. It builds one NDVI raster per tile and date and joins it to the fields it intersects. An agent connected through the Wherobots Model Context Protocol (MCP) Server can write it from a plain-language request such as "mean July NDVI for every field near Sublette". Pairing models trained on text, documents, databases, and the internet with imagery, field boundaries, and SQL is geospatial AI applied to agriculture.

NDVI research and what changes at scale

From a band ratio to a global record

The idea that leaves change the ratio of red to near-infrared light predates satellites. Carl Jordan measured light on the forest floor. He showed that the ratio of 800 nm to 675 nm light tracks leaf area, in Derivation of leaf-area index from quality of light on the forest floor (Ecology, 1969).

John Rouse, Robert Haas, John Schell, and Donald Deering at Texas A&M turned it into a normalized difference for ERTS-1, the Earth Resources Technology Satellite later renamed Landsat 1. They reported a band ratio correlated with green biomass in Monitoring vegetation systems in the Great Plains with ERTS (NASA SP-351, 1974). Compton Tucker at NASA Goddard compared combinations of red and photographic infrared bands in Red and photographic infrared linear combinations for monitoring vegetation (Remote Sensing of Environment, 1979). The normalized difference was among them. Piers Sellers then gave the index a physical basis: his radiative transfer model in Canopy reflectance, photosynthesis and transpiration (1985) linked visible and near-infrared canopy reflectance to photosynthesis and transpiration.

The NOAA AVHRR weather satellites made the index global. Pinzon and Tucker's non-stationary 1981 to 2012 AVHRR NDVI3g time series (2014) extends that calibrated record across three decades. Pettorelli and colleagues reviewed how ecologists use it in Using the satellite-derived NDVI to assess ecological responses to environmental change (2005).

NDVI time series for phenology, drought, and yield

The research that put NDVI into routine use treats each pixel or field as a time series. Zhang, Friedl, and colleagues fit S-shaped (logistic) curves to MODIS vegetation index series to monitor vegetation phenology (2003). They read green-up, maturity, senescence, and dormancy from each curve. Jönsson and Eklundh's TIMESAT (2004) became a standard tool for extracting season dates from noisy series with cloud gaps.

Wardlow, Egbert, and Kastens showed that time series of MODIS 250 m vegetation indices separate crop types in the U.S. Central Great Plains (2007). The Haskell County season tests the same principle at 10 m. On 28 July, 24 fields have a mean NDVI above 0.7, and 18 of them follow the summer-peak curve.

Adding one spring date separates the groups. Crop mapping from the timing of green-up and harvest (phenology) uses this principle.

Sensors need to agree

A long field record mixes sensors. A calibration difference between sensors can look like real change in a 40-year record. Roy and colleagues characterized the Landsat 7 to Landsat 8 NDVI difference (2016). Claverie and colleagues built the Harmonized Landsat and Sentinel-2 surface reflectance data set (2018) so the two missions can feed one time series.

Scale is a systems problem

A season of NDVI for every field in a state is billions of pixels joined to millions of polygons. Gorelick and colleagues showed planetary-scale raster analysis in Google Earth Engine (2017). Jia Yu, Jinxuan Wu, and Mohamed Sarwat introduced distributed spatial processing on Apache Spark with GeoSpark (ACM SIGSPATIAL 2015). Yu, Zhang, and Sarwat described its partitioning, indexing, and spatial join design in Spatial data management in Apache Spark: the GeoSpark perspective and beyond (GeoInformatica, 2019). GeoSpark became Apache Sedona. Its raster functions include the map algebra and zonal statistics used here, and they run the per-field summaries in parallel. Kanchan Chowdhury and Sarwat extended the stack to deep learning on raster and spatiotemporal data with GeoTorch (SIGSPATIAL 2022).

Open problems

• Greenness is not yield or health. NDVI rises with leaf area and chlorophyll. Weeds, cover crops, and a stressed but leafy canopy can all score well. Yield models add weather, soil, and crop type.
• Cloud gaps fall unevenly through the season. In the Haskell County season, the clear pixels inside the 102 fields totaled 611,103 on 7 December against 709,978 on 16 October. Cloudier regions lose more dates. Gap filling with composites, harmonized Landsat and Sentinel-2, or radar changes the curve shape.

Conclusions

1. Compute NDVI on surface reflectance and check the file's scale and offset first. A scale cancels in the ratio, and an offset shifts every value with no error raised. Mask with SCL classes 4 and 5 before the formula.
2. An NDVI value means most against a reference: the same field on earlier dates, or neighbors imaged on the same date. A wheat field at peak scores like a forest, and near zero after harvest.
3. One NDVI date does not identify a crop. Of the 24 fields above 0.7 on 28 July, 18 were summer-peak fields. Adding 14 April, when spring-peak fields read 0.82 and summer-peak fields 0.16, separates the groups. Use at least two dates for crop type or phenology.
4. NDVI flattens above about 0.8. Those pixels in the Haskell scene sit at a median red reflectance of 0.030 against 0.421 in near-infrared. Use EVI, NDRE, or the Wide Dynamic Range Vegetation Index (WDRVI) over dense canopy.
5. Compute the NDVI raster once per tile and date, repartition by scene, and broadcast the small field table. 102 fields on 15 dates ran in under four minutes.
6. Per-field NDVI curves depend on the field boundaries. 55 of the 102 fields fit none of the three peak patterns, and a mixed curve can mean a polygon holds more than one management unit. Check model-derived boundaries against parcels and the curves first.

Related terms

• Landsat: the program NDVI was first computed from, record back to 1972
• Change detection: comparing NDVI between dates to find vegetation loss
• Hyperspectral imagery: narrow-band data behind red-edge and other advanced indices
• Parcel data: land records NDVI is summarized over per property
• Earth observation data: the imagery NDVI is computed from
• Synthetic aperture radar: radar imagery unaffected by the clouds that gap NDVI time series
• Raster vs. vector data: NDVI is a raster, and field and parcel boundaries are vectors

Read more from Wherobots

• AlphaEarth Embeddings, Zonal Statistics, and PCA
• How We Delivered "Fields of The World" with RasterFlow: A Planetary-Scale GeoAI Pipeline
• Raster Data Analysis With Spatial SQL And Apache Sedona
• From the Spokane firestorm to all of Washington: real-time wildfire monitoring for under $50 a pass

Compute NDVI from Sentinel-2 and summarize it per field in SQL with a Wherobots free trial at cloud.wherobots.com.