NDVI From Free Satellite Data: A Working Pipeline for Indian Crops, Step by Step
The vegetation map looked wrong in a way nobody could name. Crops were green, as they should be, but every pond on the map read as faintly vegetated, and the whole district seemed a little less green than the same district a year earlier. The numbers were close to right everywhere, and slightly off everywhere, which is the worst kind of error.
The cause was a single constant. Since January 2022, the surface reflectance products from the Sentinel-2 satellites have carried an additive offset of minus 1,000 counts, and the pipeline was dividing by 10,000 without subtracting it. Water that should read about minus 0.33 on the vegetation index read about minus 0.04. A healthy crop that should read 0.78 read 0.54. A normalised difference is immune to a shared scale factor and sensitive to an added offset, so the map was biased in a way that varied with brightness.
This is a field guide to building a vegetation-index pipeline from free satellite data that does not fall into that and the other usual traps. It follows eight steps in a fixed format: what you do, the number to check, and how it goes wrong. The code runs and the numbers in it are computed, and where I rely on documentation or public discussion I say which. The target is a small Indian field, where a 10-metre pixel is a coarse thing and the monsoon is a long cloud.
What NDVI Is, and What It Is Not
The normalised difference vegetation index compares how much a surface reflects near-infrared light with how much it reflects red light. Healthy green leaves absorb red light for photosynthesis and reflect near-infrared strongly, so the difference between the two is large. Bare soil reflects both more evenly, and water reflects near-infrared very little.
For Sentinel-2, the red band is band 4 at about 665 nanometres and the near-infrared band is band 8 at about 842, both at 10-metre resolution. The index is the near-infrared reflectance minus the red, divided by their sum, and it runs from minus one to plus one.
Surface | Typical NDVI | What it tells you |
Water | Below zero | Not vegetation |
Bare soil, built-up | About 0.05 to 0.2 | Little or no green cover |
Sparse or young crop | About 0.2 to 0.4 | Early growth, or poor cover |
Established crop | About 0.4 to 0.7 | Active growth |
Dense healthy canopy | About 0.7 to 0.9 | Full canopy, with the index starting to saturate |
The ranges are rough guides and vary by crop, soil and season. Three things NDVI is not deserve saying at the outset. It is not yield: it measures greenness, which correlates with yield at some stages for some crops. It is not a diagnosis: a low value says something is wrong or early, and not what. And it is not precise at high canopy density, where the index flattens and stops distinguishing a good crop from a very good one.
The Pipeline in Eight Steps
Step | What you do | The number to check |
1 | Define the field and the question | Pure pixels per field |
2 | Find scenes in an archive | Clear observations per month |
3 | Read only the bands you need | Bytes read per field per scene |
4 | Convert counts to reflectance | Offset applied, no-data respected |
5 | Mask clouds and shadows | Valid pixel share per field |
6 | Average over the field | Pixels used, buffer applied |
7 | Fill gaps and smooth | Gap length filled |
8 | Extract season metrics | Start, peak, end |
The pipeline is eight small steps. A mistake in any one of them produces a map that looks plausible, so each step needs its own number to check.
Step 1: Define the Field and the Question
What you do. Obtain the field boundary as a polygon from a survey, a land record, a mobile app or your own digitising. Decide the question: Was the crop sown, and when? Is it growing as well as its neighbours? Did it reach a good peak? Each question needs a different set of observations.
The number to check. The count of pure pixels, meaning pixels wholly inside the field and not touching the boundary, where neighbouring land, bunds and roads contaminate the signal. With 10-metre pixels and a 10-metre inner buffer, the count falls fast as the field shrinks.
import math
def pure_pixels(area_ha, buffer_m=10, pix=10):
side = math.sqrt(area_ha * 10000) # treat the field as a square
inner = max(0.0, side - 2 * buffer_m)
return round((inner / pix) ** 2)
for a in (0.2, 0.5, 1.08, 2.0, 5.0):
print(a, pure_pixels(a), round(a * 10000 / 100))
# 0.2 ha: 6 pure of 20 0.5: 26 of 50 1.08: 70 of 108 2.0: 147 of 200 5.0: 415 of 500
Field size | Pixels in total | Pure pixels after a 10 m buffer |
0.2 hectare | 20 | 6 |
0.5 hectare | 50 | 26 |
1.08 hectares | 108 | 70 |
2 hectares | 200 | 147 |
5 hectares | 500 | 415 |
A field of the size of the national average holding, about 1.08 hectares, still has about 70 pure pixels, which is enough for a stable average. A field of a fifth of a hectare has six, which is thin, and a single cloud-edge pixel can move the average. Fields that small need an additional approach, such as combining observations from a group of neighbouring fields of the same crop or treating the result with wider error bars.
How it goes wrong. Using the polygon without checking it. A boundary drawn from an old map, or from a phone's single GPS fix, can be shifted by 10 or 20 metres, which at this resolution moves the field into its neighbour.
Choosing a Data Route: Official Archive, Cloud Mirror or Earth Engine
Three routes lead to the same satellite data, and the choice shapes the rest of the pipeline.
The official Copernicus Data Space Ecosystem is the authoritative source. It holds the full archive and, now, a reprocessed collection in which the radiometric offset is applied consistently across dates. It offers catalogue searches and partial reads, and it is the right choice when provenance matters. Public cloud mirrors expose the same data in cloud-optimised files through a standard catalogue interface, which makes reading a small window from a large archive fast and cheap, and they often add convenient metadata. Google Earth Engine offers a harmonised copy together with a computing environment, which is quick for exploration and for one-off analyses and ties you to its tools and terms.
Route | Strength | Watch out for |
Official archive | Authoritative, reprocessed collection, full metadata | Rate limits and account set-up |
Cloud mirror with a standard catalogue | Fast partial reads, easy to script | Check which processing baseline each item carries |
Earth Engine harmonised collection | Fast exploration, built-in computing | Offset already shifted, so do not apply it twice, and lock-in |
Whichever route you pick, write down which collection you used and record the scene identifier and processing baseline with every observation. Mixing sources without noting this is how the offset error returns, and it is invisible in the output.
Choose the route for provenance first and speed second. A pipeline that cannot say which processing baseline produced a number cannot be audited.
Step 2: Find Scenes: Archives, Cloud and the Monsoon
What you do. Search an archive for scenes that cover your field in your season. The main routes are the official Copernicus Data Space Ecosystem, public cloud mirrors that expose the data through a standard catalogue interface, and Google Earth Engine, which holds a harmonised copy. Sentinel-2 revisits the same ground about every five days with its satellites, so a month contains about six observations of any tile.
The number to check. The expected number of clear observations per month, and the chance of none.
def p_no_clear(p_cloud, n):
return round(p_cloud ** n * 100, 1)
for pc in (0.3, 0.6, 0.8):
print(pc, p_no_clear(pc, 6), round(6 * (1 - pc), 1))
# cloud probability 0.3: 0.1 percent of months with no clear scene, 4.2 clear on average
# 0.6: 4.7 percent, 2.4 clear 0.8: 26.2 percent, 1.2 clear
Chance a scene is cloudy | Clear scenes per month on average | Months with no clear scene |
30 percent | 4.2 | 0.1 percent |
60 percent | 2.4 | 4.7 percent |
80 percent | 1.2 | 26.2 percent |
In a monsoon month where four scenes in five are cloudy over your area, about one month in four has no clear observation at all, and the crop season is made of such months. The curve you extract will have gaps, and the gaps will fall exactly when the crop is growing fastest.
How it goes wrong. Filtering scenes by the cloud percentage in the metadata. That number describes the whole tile, which is over 100 kilometres across, and your field may be clear in a scene the metadata calls 70 percent cloudy. Filter on the field, not the tile, using the per-pixel classification in Step 5. For the cloudiest weeks, radar is the fallback, as our post on satellite radar explains, because it sees through cloud, though it measures structure and moisture and not greenness.
Step 3: Read Only the Bands You Need
What you do. Request only the red band, the near-infrared band and the scene classification, and only for the window around your field. The data is stored in cloud-friendly files that allow partial reads, so a small field needs a few kilobytes, not a gigabyte tile.
The number to check. Bytes read per field per scene. A pipeline that downloads entire tiles for a one-hectare field is doing a thousand times the work, and will not scale.
Know the resolutions. Bands 4 and 8 are at 10 metres. The scene classification is at 20 metres, so each of its pixels covers four of the others, and it must be resampled by nearest neighbour, never by interpolation, since a class code is a label and not a quantity.
How it goes wrong. Resampling the classification with an averaging method, which turns class 4 beside class 9 into class 6 or 7, a meaningless blend that sometimes happens to look like water.
Step 4: Convert to Reflectance: The Offset That Breaks Pipelines
What you do. Convert the stored integer counts to reflectance, using the scale and the offset in the product metadata. From the processing baseline 04.00, operational from 25 January 2022, the level-2A products carry a band-dependent additive offset so that very dark surfaces can have negative reflectance. The offset is initially minus 1,000 counts for all bands. The official conversion is reflectance equals counts plus the offset, divided by the quantification value of 10,000. A count of zero means no data.
The official documentation advises reading the offset and the quantification value from the metadata of each product and not hard-coding them. The processing baselines have changed over the years, and the Copernicus Data Space archive now holds a reprocessed collection in which the offset is present consistently across dates. Collections that harmonise the data, such as the one in Google Earth Engine, shift newer scenes back into the old range, in which case you must not subtract the offset a second time.
The number to check. The same pixel, computed both ways.
import numpy as np
def to_reflectance(dn, offset=-1000.0, quant=10000.0):
dn = np.asarray(dn, dtype=float)
r = (dn + offset) / quant
return np.where(dn == 0, np.nan, r) # a count of zero is no data
def ndvi(red_dn, nir_dn, offset=-1000.0):
r, n = to_reflectance(red_dn, offset), to_reflectance(nir_dn, offset)
return (n - r) / (n + r)
print(np.round(ndvi([1500], [5000]), 3), np.round(ndvi([1500], [5000], 0.0), 3)) # crop: [0.778] with offset, [0.538] without
print(np.round(ndvi([1200], [1100]), 3), np.round(ndvi([1200], [1100], 0.0), 3)) # water: [-0.333] with, [-0.043] without
print(ndvi([0], [5000])) # [nan]
The crop pixel, with true reflectances of 0.05 in the red and 0.40 in the near-infrared, gives an index of 0.778. Ignoring the offset gives 0.538, a 31 percent underestimate that looks like a plausible weak crop. The water pixel gives minus 0.333 correctly and minus 0.043 incorrectly, which turns open water into a very sparse surface. The error is largest where the surface is darkest, and a public discussion of the issue reported the median index over a scene shifting from 0.39 to 0.58 when the offset was applied properly to the same pixels.
A normalised difference ignores a shared scale factor and is sensitive to an added offset. If your pipeline reads the offset from the product metadata, this class of error cannot occur.
How it goes wrong. Three ways. Hard-coding the old convention, so scenes from before 2022 and after are treated differently and a time series shows a step at January 2022. Subtracting the offset from a harmonised collection that has already applied it. And applying the offset to the count of zero, which turns no-data into a plausible negative number.
Step 5: Mask Clouds and Shadows With the Scene Classification
What you do. Use the scene classification band, which assigns each 20-metre pixel a class, to drop pixels that are not clear land.
Class value | Meaning | Keep for vegetation work |
0 | No data | No |
1 | Saturated or defective | No |
2 | Dark area pixels | Usually no |
3 | Cloud shadows | No |
4 | Vegetation | Yes |
5 | Bare soils | Yes |
6 | Water | No, unless monitoring paddy flooding |
7 | Unclassified or low-probability cloud | With care |
8 | Medium-probability cloud | No |
9 | High-probability cloud | No |
10 | Thin cirrus | No |
11 | Snow or ice | No |
The number to check. The valid pixel share per field per date. A date with under about a third of the field's pure pixels valid should be dropped, because the few pixels left are likely to be the ones the classifier got wrong.
GOOD = {4, 5}
def masked_ndvi(red_dn, nir_dn, scl, keep=GOOD):
v = ndvi(red_dn, nir_dn)
return np.where(np.isin(scl, list(keep)), v, np.nan)
scl = np.array([4, 4, 9, 3, 5, 10])
red = np.array([1500, 1600, 1500, 1500, 2500, 1500])
nir = np.array([5000, 5200, 5000, 5000, 3500, 5000])
print(np.round(masked_ndvi(red, nir, scl), 3)) # [0.778 0.75 nan nan 0.25 nan]
Of six pixels, three are dropped: one high-probability cloud, one cloud shadow and one thin cirrus. The remaining three describe vegetation and bare soil.
How it goes wrong. Relying on the older cloud-mask band, which was not produced for a period from January 2022 to February 2024, and which stopped carrying polygons in the form older tools expect. Use the scene classification, and widen the mask by one pixel around clouds and shadows, since classifiers miss the faint edges where reflectance is most contaminated.
Step 6: Average Over the Field, Not Over the Tile
What you do. For each clear date, take the pure pixels inside the field and compute one value. Use the median where there may be stray contamination, and the mean where the field is homogeneous. Record the number of pixels used.
The number to check. Pixels used, and the share of the field's pure pixels they represent. A field value built on eight pixels from a field of 70 is a weak value, and the pipeline should say so.
Choice | Use when | Risk |
Mean | A uniform field with clean masking | One contaminated pixel moves it |
Median | Stray cloud edges or a few weedy patches | Hides real within-field variation |
Mean plus spread | You want to see variability | Needs enough pixels to be meaningful |
How it goes wrong. Averaging over the polygon without the buffer, which includes the bund and the neighbouring field, and so blends two crops with different calendars into a curve that belongs to neither.
Step 7: Fill the Gaps and Smooth the Curve
What you do. The series of field values has gaps where no clear observation exists. Fill short gaps by interpolation and smooth the result, so that noise does not look like a crop event.
The sketch below interpolates across missing values and applies a three-point median filter. More elaborate methods fit a smooth seasonal curve, such as a logistic or harmonic model, and suit longer series. A simple method is easier to audit, and is a good start.
def fill_and_smooth(series):
s = np.array(series, dtype=float)
idx = np.arange(len(s)); ok = ~np.isnan(s)
s = np.interp(idx, idx[ok], s[ok])
pad = np.pad(s, 1, mode="edge")
return np.array([np.median(pad[i:i + 3]) for i in range(len(s))])
ser = [0.20, 0.25, np.nan, 0.50, 0.55, 0.60, np.nan, np.nan, 0.40, 0.30]
print(np.round(fill_and_smooth(ser), 3))
# [0.2 0.25 0.375 0.5 0.55 0.55 0.533 0.467 0.4 0.3 ]
The number to check. The longest gap filled. A gap of one observation, about five days, is a safe fill. A gap of four or five observations, such as a cloudy stretch in July, is a guess, and the output should be flagged as interpolated so that nobody mistakes it for a measurement.
How it goes wrong. Smoothing so hard that the peak is flattened. In the example, the median filter reduced a true peak of 0.60 to 0.55. For season timing that matters little, and for peak greenness it matters, so choose the smoothing for the question.
Step 8: Turn the Curve Into Season Metrics
What you do. From the smoothed curve, extract a few numbers that answer the question: the date when the index first rises above a threshold, taken as the start of the season; the date and value of the peak; and the date it last stays above the threshold, taken as the end.
def season_metrics(days, values, thr=0.4):
v = np.array(values); d = np.array(days)
above = np.where(v >= thr)[0]
if len(above) == 0:
return None
return {"start_day": int(d[above[0]]), "peak_day": int(d[int(np.argmax(v))]),
"peak": round(float(v.max()), 3), "end_day": int(d[above[-1]])}
days = [165, 175, 185, 195, 205, 215, 225, 235, 245, 255] # day of year
vals = fill_and_smooth([0.18, 0.22, np.nan, 0.48, 0.62, 0.70, np.nan, np.nan, 0.45, 0.30])
print(season_metrics(days, vals)) # {'start_day': 195, 'peak_day': 205, 'peak': 0.62, 'end_day': 245}
For this field, the season rose above 0.4 on day 195, peaked on day 205 and fell below on day 245. Day numbers are used so that the code is not tied to a calendar. Compare these dates with the farmer's sowing and harvest records and the pattern of neighbouring fields.
The number to check. The agreement between the extracted start date and any known sowing date, over a sample of fields. A consistent lag of two or three weeks is normal, because the index rises only once the canopy covers enough soil, and the lag can be calibrated per crop.
How it goes wrong. Fixed thresholds across crops and regions. A threshold of 0.4 suits a crop that reaches 0.8 and fails for one that peaks at 0.5. Where it matters, set thresholds as a fraction of each field's own peak.
Seven Numbers to Print for Every Field
A pipeline that prints a few checks per field catches most failures before they reach a customer. Add these to the output table and review them weekly.
Number | Healthy range | If it is outside |
Pure pixels | At least 30 for a normal field | Treat the field as low confidence |
Valid pixel share on the date | Above about a third | Drop the observation |
Observations in the season | Enough for a curve, at least 8 | Flag the field as under-observed |
Longest filled gap | A few observations at most | Mark the interpolated stretch |
NDVI on a known water body in the same scene | Below zero | Suspect the offset or the mask |
Step at January 2022 in long series | None | Check processing baselines |
Agreement of start date with sowing records | Within a consistent lag | Recalibrate the threshold for the crop |
The water-body check deserves special mention. Place a small reference polygon on a stable lake or reservoir in each tile and print its NDVI each date. It should sit below zero every time. A drift toward zero is the earliest warning of the offset error, a haze problem or a masking fault, and it costs nothing to run.
Validate Against the Ground
A curve that looks like a crop is not evidence that it is one. Validation turns the pipeline from a plausible picture into a measurement.
Collect ground truth for a sample of fields, chosen to span crops, soils and management: sowing date, harvest date, crop type and variety, irrigation, and a few photographs through the season. Prefer a few well-documented fields to many poorly documented ones. Compare the extracted start, peak and end with the records and compute the error and its spread by crop. If the start date is wrong by 30 days for a quarter of fields, find out why before trusting the rest.
Watch for systematic causes of disagreement. Double-cropped fields, where a second crop follows quickly, produce two peaks that a single-season method merges. Intercropped and tree-crop fields never reach the same bare-soil baseline. Flooded paddy has a negative index before the seedlings emerge. Weeds can produce a green signal before the crop. Each is a reason to refine the method for a crop, and none is a reason to abandon it.
A vegetation curve is a hypothesis about a field. A sowing date from the farmer is evidence. Compare them on a sample before you scale.
What NDVI Cannot Tell You
Four limits matter in practice.
It saturates at high canopy density, so a good crop and an excellent one look alike. Indices that use the red-edge bands, which sit between red and near-infrared, remain sensitive to higher chlorophyll and are often better at the top of the curve.
It responds to the soil background at low cover. A pale soil under a sparse crop gives a different index from a dark soil under the same crop, which makes comparisons between farms unfair early in the season.
It confuses causes. A low value may mean late sowing, poor germination, water stress, nutrient shortage, pests or simply a different variety. Telling them apart needs other information: weather, soil moisture, field notes and, where available, other bands and indices.
And at 10 metres it mixes pixels in small fields. The unit of analysis is the field, and the field is often barely larger than a handful of pixels.
Question | Better served by |
Was the crop sown, and when | Season metrics from the NDVI curve |
Is the canopy dense and healthy at peak | Red-edge based indices |
Is the crop under water stress | Indices using shortwave infrared, with weather data |
Is the field flooded | Water indices, or radar |
What will the yield be | A calibrated model with local data, not NDVI alone |
Scaling From One Field to a Million
A script that works for one field will not serve a district without three changes.
Process by tile and by date, not by field. Read the bands once for a tile, then extract every field in it, which avoids repeated reads. Cache intermediate results, so that a new date adds one observation per field without recomputing the past. And store the output as a table of one row per field per date, with the value and the quality information, not as images.
Record provenance. For each observation, keep the scene identifier, the processing baseline, the offset applied, the number of valid pixels and the masking rules. A year later, when a processing baseline changes or a bug is found, this is what lets you tell which values to recompute.
Per-observation field | Why it is stored |
Field identifier and date | The key |
NDVI value and pixel count | The result and its support |
Valid pixel share | To filter weak observations |
Scene identifier and processing baseline | To trace and recompute after a change |
Offset and scale applied | To catch the offset error |
Interpolated flag | To separate measurement from fill |
Compute grows with the number of fields and the number of dates, and storage is modest if you keep only the table. The main cost is engineering time to make the pipeline reliable, which is why the checks in each step earn their keep.
From Curve to Advice: What a Farmer Can Use
Season metrics become useful when they turn into messages a farmer or an agronomist can act on. The step needs more care than the pipeline, because a wrong message costs more than a wrong number.
Three kinds of message work well. The first is a record: this field's crop emerged around this date, peaked at this level and began to decline on this date, which is a useful history for insurers, lenders and buyers. The second is a comparison: this field is two weeks behind its neighbours, or its peak is lower than the group, which invites a visit. The third is an alert: the index has fallen faster than expected in the last three observations, which asks someone to look.
Three kinds of message are risky. A yield estimate from NDVI alone is rarely reliable enough to act on. A diagnosis, such as saying a field has a particular disease, claims knowledge the index does not have. And a prescription, such as a quantity of fertiliser, needs soil and crop information that the satellite does not hold.
Message | Safe to send from NDVI alone | Needs more |
Emergence and season dates | Yes, with a stated lag | Sowing record to calibrate |
Behind or ahead of neighbours | Yes | Reason needs a visit |
Sudden decline alert | As a prompt to inspect | Cause needs ground information |
Yield forecast | No | Calibrated model with local data |
Disease or nutrient diagnosis | No | Field inspection and tests |
Keep the human in the loop. The satellite is a good way to decide which of a thousand fields to visit this week, and a poor substitute for the visit.
A Composite: Two Neighbouring Fields, Two Different Stories
This is a representative composite, genericised, from several agricultural monitoring projects. The numbers are illustrative.
Two neighbouring fields of about a hectare each grew the same crop on adjoining plots. Their curves, extracted with the pipeline above, looked different. One rose above the threshold on day 195 and peaked at 0.72. The other rose on day 215 and peaked at 0.58.
The first reaction was that the second farmer had a problem. The team checked before drawing any conclusion. The pure pixel counts were adequate, 68 and 71. The valid pixel share was high on most dates for both. The processing baseline was the same throughout, with the offset read from metadata.
The farmer's records settled it. The second field had been sown 18 days later because irrigation reached it later, and the lag in the curve matched the lag in sowing almost exactly. Its lower peak, however, was a real difference. Field notes showed a nitrogen top-dressing missed in the second field's early stage, and a visit confirmed pale leaves in one corner.
The curves had told a truthful story, and it took two kinds of ground information to read it: a sowing date to explain the delay and a visit to explain the peak. Neither alone would have done. The project's lesson was to keep the farmer's calendar and a short field note beside every curve.
Pitfalls That Look Like Results
Pitfall | What you see | The check |
Offset ignored | Water slightly vegetated, whole district too pale | Read the offset from metadata and test on a water pixel |
Step at January 2022 | A jump in every long series | Compare processing baselines across the series |
Harmonised data corrected twice | Values too high | Check the collection's documentation |
Class codes interpolated | Strange masks | Resample classification by nearest neighbour |
Tile-level cloud filter | Good scenes thrown away | Filter on the field |
Unbuffered boundary | Blended curves | Apply an inner buffer |
Heavy smoothing | Flattened peaks | Compare raw and smoothed |
Fixed thresholds | Wrong dates for some crops | Use a fraction of each field's peak |
A Four-Week Build Plan
Week one: collect 20 to 30 fields with boundaries and ground records, and compute pure pixel counts.
Week two: build steps 2 to 6 for one tile and one season, with the offset read from metadata and the checks printed for each date.
Week three: add the fill, the smoothing and the season metrics, and compare against the ground records.
Week four: write the provenance table, run the pitfall checks, and process a second season to confirm that the pipeline behaves the same way.
The Constant You Cannot Forget
The vegetation map that was wrong everywhere by a little was wrong by one number. The fix took a line of code, and finding it took a week, because nothing in the output announced it.
Build the pipeline so that each step has its own check, store the facts that would let you find the next such error, and compare every curve with a sowing date from the ground. Then the map will be right for the reasons you can explain.
Related service: AgriTech and Precision Agriculture
Related reading: InSAR Explained: What Satellite Radar Can and Cannot Measure on the Ground and Precision Agriculture in India: What to Measure First on a 1.08-Hectare Farm.
EurthTech is an engineering-first product design and embedded systems company with a development center in Hyderabad. To discuss a project, write to connect@eurthtech.com or see https://www.eurthtech.com/case-studies.
Comments