Agricultural Data Science

---

Introduction: Why Every Agricultural Scientist Should Know NDVI

If you work in agriculture, chances are you’ve heard the term NDVI. Normalized Difference Vegetation Index — a mouthful, but one of the most powerful and widely-used tools in modern agricultural science.

Think of NDVI as a crop health report card, derived entirely from satellite images. In one number — ranging from -1 to +1 — it tells you whether your field is thriving, stressed, or bare.

But here’s the gap I see all the time among agricultural researchers in India: they understand NDVI conceptually, but they depend entirely on paid software like ERDAS Imagine or commercial platforms to compute it. That’s expensive, slow, and inflexible.

With Python, you can calculate NDVI for any field, anywhere in India, for free — in under 20 lines of code.

In this post, I’ll show you exactly how, using real Sentinel-2 satellite data. By the end, you’ll be able to:

• Understand what NDVI measures and why it matters for agricultural research
• Download and read multispectral satellite imagery in Python
• Calculate NDVI using the rasterio and numpy libraries
• Visualise NDVI maps with matplotlib
• Apply NDVI analysis to real agricultural problems (crop monitoring, drought assessment, yield estimation)

I use this workflow regularly in my own research. Let’s get into it.

---

What is NDVI and Why Does It Matter?

NDVI stands for Normalized Difference Vegetation Index. It was developed by NASA scientists in the 1970s and remains the most widely used vegetation index in the world today.

The formula is simple:

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

Where:

• NIR = Near-Infrared band reflectance
• Red = Red band reflectance

Healthy green vegetation absorbs red light for photosynthesis and strongly reflects near-infrared light. Stressed or sparse vegetation does the opposite. This contrast is what NDVI captures.

NDVI value interpretation:

NDVI Range	What It Means	Agricultural Implication
0.8 to 1.0	Dense, healthy vegetation	Peak crop canopy, excellent health
0.6 to 0.8	Moderate-high vegetation	Good crop growth
0.4 to 0.6	Moderate vegetation	Average crop health, possible stress
0.2 to 0.4	Sparse vegetation	Stress, early season, or thin canopy
0.0 to 0.2	Very sparse/bare soil	Bare ground, very early stage
Below 0.0	Non-vegetation	Water bodies, built-up areas, clouds

In Indian agricultural contexts, NDVI is used for:

1. Crop condition monitoring — Is the paddy crop healthy in July? Are wheat fields uniform in November?
2. Drought and stress detection — Which districts show early signs of moisture stress?
3. Yield forecasting — NDVI at critical crop growth stages correlates strongly with final yield
4. Kharif/Rabi crop mapping — Identifying crop types across a region based on NDVI time series
5. Assessment of PMFBY claims — Insurance companies and state governments use NDVI to validate crop loss claims

We use NDVI-based analyses to support everything from crop condition reports to state-level advisories. Once you know how to compute it in Python, you unlock all of this yourself.

---

What Satellite Data Will We Use?

We’ll use Sentinel-2 imagery from the European Space Agency (ESA). It’s free, it covers all of India, and it has a 10-metre spatial resolution for the bands we need — fine enough to see individual farm fields.

Key Sentinel-2 bands for NDVI:

Band	Name	Wavelength	Resolution	Role in NDVI
Band 4	Red	665 nm	10 m	Denominator + difference
Band 8	NIR	842 nm	10 m	Numerator + difference

You can download Sentinel-2 data freely from:

• Copernicus Open Access Hub: https://scihub.copernicus.eu/
• Google Earth Engine (for larger areas or time series)
• AWS Open Data Registry: free access, no sign-up needed for many datasets
• USGS EarthExplorer: also carries some ESA products

For this tutorial, I’ll walk you through reading a pre-downloaded GeoTIFF. If you’d like a separate post on how to download Sentinel-2 data from Copernicus Hub using Python, let me know in the comments.

---

Setting Up Your Python Environment

Before writing any code, let’s make sure the right libraries are installed.

pip install rasterio numpy matplotlib geopandas

If you’re working on a government server or a restricted environment, you can also run this in Google Colab for free — no installation needed.

Libraries we’ll use:

• rasterio — read and write geospatial raster files (GeoTIFFs)
• numpy — numerical operations, including the NDVI formula
• matplotlib — visualisation and NDVI map plotting
• geopandas — optional, for masking to a specific study area

---

Step 1: Import Libraries and Load the Satellite Bands

import rasterio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import warnings
warnings.filterwarnings('ignore')

# Define paths to your Sentinel-2 band files
# Replace these with your actual file paths
red_band_path = "T44QKF_20231015_B04_10m.tif"   # Band 4 (Red)
nir_band_path  = "T44QKF_20231015_B08_10m.tif"   # Band 8 (NIR)

# Open the Red band
with rasterio.open(red_band_path) as red_src:
    red_band = red_src.read(1).astype(float)    # Read as float for division
    profile = red_src.profile                    # Save metadata for output file
    print(f"Red band shape: {red_band.shape}")
    print(f"CRS: {red_src.crs}")
    print(f"Transform: {red_src.transform}")

# Open the NIR band
with rasterio.open(nir_band_path) as nir_src:
    nir_band = nir_src.read(1).astype(float)
    print(f"NIR band shape: {nir_band.shape}")

What this does: We read both bands as 2D numpy arrays. The .astype(float) is important — integer arrays would give wrong results during division.

Note: Sentinel-2 filenames follow the format T[tile]_[date]_[band]_[resolution].tif. The tile code for Assam and northeast India is typically T46RBN, T46RBM, or T46QBL depending on your study area.

---

Step 2: Calculate NDVI

This is the core calculation. Two lines of Python.

# Avoid division by zero — set to NaN where both bands are zero
denominator = nir_band + red_band
denominator[denominator == 0] = np.nan

# Calculate NDVI
ndvi = (nir_band - red_band) / denominator

# Quick summary statistics
print(f"NDVI Statistics:")
print(f"  Min:  {np.nanmin(ndvi):.4f}")
print(f"  Max:  {np.nanmax(ndvi):.4f}")
print(f"  Mean: {np.nanmean(ndvi):.4f}")
print(f"  Std:  {np.nanstd(ndvi):.4f}")

Sample output for an agricultural area in Assam (October, post-Kharif):

NDVI Statistics:
  Min:  -0.2341
  Max:   0.8923
  Mean:  0.4512
  Std:   0.1876

A mean NDVI of 0.45 in October for this area is consistent with late Kharif paddy — the crop is still standing but past its peak greenness.

---

Step 3: Visualise the NDVI Map

A raw array of numbers isn’t useful for agricultural interpretation. We need a colour-coded map.

# Set up a professional NDVI colour map
# RdYlGn: Red (low/stressed) → Yellow (moderate) → Green (healthy)
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.suptitle("NDVI Analysis — Agricultural Area, Assam (October 2023)",
             fontsize=14, fontweight='bold')

# Panel 1: Red Band (visible light)
ax1 = axes[0]
im1 = ax1.imshow(red_band, cmap='Reds_r', vmin=0, vmax=3000)
ax1.set_title("Red Band (Band 4)", fontsize=11)
ax1.axis('off')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04, label="Reflectance")

# Panel 2: NIR Band
ax2 = axes[1]
im2 = ax2.imshow(nir_band, cmap='YlOrRd', vmin=0, vmax=5000)
ax2.set_title("NIR Band (Band 8)", fontsize=11)
ax2.axis('off')
plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04, label="Reflectance")

# Panel 3: NDVI Map
ax3 = axes[2]
im3 = ax3.imshow(ndvi, cmap='RdYlGn', vmin=-0.3, vmax=0.9)
ax3.set_title("NDVI — Vegetation Health Index", fontsize=11)
ax3.axis('off')
cbar = plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04)
cbar.set_label("NDVI Value", rotation=270, labelpad=15)
cbar.set_ticks([-0.3, 0.0, 0.2, 0.4, 0.6, 0.8])
cbar.set_ticklabels(['Water/Cloud', 'Bare Soil', 'Sparse', 'Moderate', 'Good', 'Dense'])

plt.tight_layout()
plt.savefig("ndvi_analysis_assam.png", dpi=200, bbox_inches='tight')
plt.show()
print("Map saved as ndvi_analysis_assam.png")

---

Step 4: Classify NDVI into Agricultural Categories

Raw NDVI values are useful, but for field reports and policy briefs, classified maps are more interpretable. Let’s create a crop health classification.

# Create NDVI classification for agricultural interpretation
def classify_ndvi_agriculture(ndvi_array):
    """
    Classify NDVI into agricultural health categories.
    Suitable for Kharif and Rabi crop monitoring in India.
    """
    classified = np.full(ndvi_array.shape, np.nan)
    
    classified[ndvi_array < 0.0]                       = 0  # Water / Non-veg
    classified[(ndvi_array >= 0.0) & (ndvi_array < 0.2)] = 1  # Bare/Very sparse
    classified[(ndvi_array >= 0.2) & (ndvi_array < 0.4)] = 2  # Sparse / stressed crop
    classified[(ndvi_array >= 0.4) & (ndvi_array < 0.6)] = 3  # Moderate crop health
    classified[(ndvi_array >= 0.6) & (ndvi_array < 0.8)] = 4  # Good crop health
    classified[ndvi_array >= 0.8]                       = 5  # Excellent / dense canopy
    
    return classified

ndvi_classified = classify_ndvi_agriculture(ndvi)

# Create a custom colour map for the classified map
class_colors = ['#2166AC',  # 0 — Blue: Water/Non-veg
                '#D7191C',  # 1 — Red: Bare/Very sparse
                '#FDAE61',  # 2 — Orange: Stressed crop
                '#FEE08B',  # 3 — Yellow: Moderate
                '#A6D96A',  # 4 — Light green: Good
                '#1A9641']  # 5 — Dark green: Excellent

cmap_custom = mcolors.ListedColormap(class_colors)
bounds = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]
norm = mcolors.BoundaryNorm(bounds, cmap_custom.N)

fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(ndvi_classified, cmap=cmap_custom, norm=norm)
ax.set_title("NDVI Crop Health Classification\nAgricultural Area, Assam (October 2023)",
             fontsize=13, fontweight='bold')
ax.axis('off')

cbar = plt.colorbar(im, ax=ax, fraction=0.03, pad=0.04)
cbar.set_ticks([0, 1, 2, 3, 4, 5])
cbar.set_ticklabels(['Water / Non-veg',
                      'Bare / Very Sparse',
                      'Stressed Crop',
                      'Moderate Health',
                      'Good Health',
                      'Excellent / Dense'])
cbar.set_label("Crop Health Category", rotation=270, labelpad=20)

plt.tight_layout()
plt.savefig("ndvi_classified_map.png", dpi=200, bbox_inches='tight')
plt.show()

---

Step 5: Calculate Area Statistics by Category

For an agricultural report, you’ll want to know how many hectares fall in each category.

# Calculate pixel counts and area for each class
# Sentinel-2 at 10m resolution: each pixel = 10m × 10m = 0.01 hectares

pixel_area_ha = 0.01  # hectares per pixel (10m × 10m)

categories = {
    0: "Water / Non-veg",
    1: "Bare / Very Sparse",
    2: "Stressed Crop",
    3: "Moderate Health",
    4: "Good Health",
    5: "Excellent / Dense"
}

print("=" * 55)
print("NDVI AREA STATISTICS — Assam Agricultural Area")
print("=" * 55)
print(f"{'Category':<25} {'Pixels':>10} {'Area (ha)':>12} {'%':>8}")
print("-" * 55)

total_valid_pixels = np.sum(~np.isnan(ndvi_classified))

for class_val, class_name in categories.items():
    pixel_count = np.sum(ndvi_classified == class_val)
    area_ha = pixel_count * pixel_area_ha
    pct = (pixel_count / total_valid_pixels) * 100 if total_valid_pixels > 0 else 0
    print(f"{class_name:<25} {pixel_count:>10,} {area_ha:>12,.1f} {pct:>7.1f}%")

print("-" * 55)
print(f"{'TOTAL':<25} {total_valid_pixels:>10,} "
      f"{total_valid_pixels * pixel_area_ha:>12,.1f} {'100.0':>8}%")
print("=" * 55)

Sample output:

=======================================================
NDVI AREA STATISTICS — Assam Agricultural Area
=======================================================
Category                  Pixels      Area (ha)        %
-------------------------------------------------------
Water / Non-veg            12,450        124.5      5.2%
Bare / Very Sparse          8,230         82.3      3.4%
Stressed Crop              18,750        187.5      7.8%
Moderate Health            64,300        643.0     26.8%
Good Health               98,760        987.6     41.2%
Excellent / Dense          37,320        373.2     15.6%
-------------------------------------------------------
TOTAL                     239,810      2,398.1    100.0%
=======================================================

This kind of output goes directly into a district crop situation report or a field survey planning document.

---

Step 6: Save the NDVI Raster as a GeoTIFF

For use in GIS software (QGIS, ArcGIS) or further analysis, save the NDVI output as a proper georeferenced GeoTIFF.

# Save NDVI as GeoTIFF with geospatial metadata preserved
output_path = "ndvi_output_assam.tif"

# Update profile for float32 single-band output
profile.update({
    'dtype': rasterio.float32,
    'count': 1,
    'nodata': -9999
})

# Replace NaN with the nodata value before writing
ndvi_to_save = ndvi.copy()
ndvi_to_save[np.isnan(ndvi_to_save)] = -9999

with rasterio.open(output_path, 'w', **profile) as dst:
    dst.write(ndvi_to_save.astype(rasterio.float32), 1)

print(f"NDVI GeoTIFF saved successfully: {output_path}")
print(f"File size: {os.path.getsize(output_path) / (1024*1024):.1f} MB")

You can now open ndvi_output_assam.tif directly in QGIS for further GIS analysis, or share it with colleagues who don’t use Python.

---

Real-World Application: Monitoring Paddy Crop Health in Assam

Let me show you how this workflow plays out in a real agricultural scenario.

Scenario: It’s late October. The Kharif paddy crop is in its grain-filling stage across Kamrup district, Assam. There were reports of drought stress in some blocks. The district agriculture officer wants to know which blocks are showing low NDVI before planning the relief package.

Using the workflow above, here’s what you can do in a single afternoon:

1. Download Sentinel-2 images for October 2023 for Kamrup district (free from Copernicus Hub)
2. Calculate NDVI for each 10m pixel across the district
3. Clip the raster to administrative block boundaries using geopandas
4. Compute mean NDVI per block
5. Map the blocks by average NDVI — immediately showing which blocks are stressed

import geopandas as gpd
from rasterio.mask import mask
import json

# Load block-level shapefile for Kamrup district
# (Available from ICRISAT VDSA or state GIS portals)
blocks_gdf = gpd.read_file("kamrup_blocks.shp")

# Ensure same CRS as the raster
blocks_gdf = blocks_gdf.to_crs("EPSG:32646")  # UTM Zone 46N for Assam

# Compute mean NDVI per block
block_ndvi_results = []

with rasterio.open(nir_band_path) as src:
    for idx, row in blocks_gdf.iterrows():
        geom = [json.loads(row.geometry.json())]
        
        try:
            # Open raster and mask to block boundary
            with rasterio.open("ndvi_output_assam.tif") as ndvi_src:
                masked_ndvi, _ = mask(ndvi_src, geom, crop=True)
                valid_vals = masked_ndvi[masked_ndvi != -9999]
                
                if len(valid_vals) > 0:
                    mean_ndvi = float(np.nanmean(valid_vals))
                else:
                    mean_ndvi = np.nan
                    
        except Exception:
            mean_ndvi = np.nan
        
        block_ndvi_results.append({
            'block_name': row['BLOCK_NAME'],
            'mean_ndvi': mean_ndvi
        })

# Add results back to GeoDataFrame
import pandas as pd
ndvi_df = pd.DataFrame(block_ndvi_results)
blocks_gdf = blocks_gdf.merge(ndvi_df, left_on='BLOCK_NAME', right_on='block_name')

# Identify stressed blocks (mean NDVI < 0.4 in October = concern)
stressed_blocks = blocks_gdf[blocks_gdf['mean_ndvi'] < 0.4]
print(f"\nBlocks showing possible crop stress (NDVI < 0.40):")
print(stressed_blocks[['BLOCK_NAME', 'mean_ndvi']].sort_values('mean_ndvi'))

This analysis — which would take days in traditional software — runs in minutes with Python. And it’s fully reproducible, documented, and shareable.

---

NDVI Time Series: Tracking Crop Growth Through the Season

One of the most powerful applications of NDVI is monitoring how it changes across the crop season. Here’s a simple example using multi-date NDVI values:

import matplotlib.pyplot as plt
import numpy as np

# NDVI values for a paddy field in Assam — Kharif 2023
# (You would compute these from actual multi-date satellite images)
dates = ['June 15', 'July 01', 'July 15', 'Aug 01', 'Aug 15',
         'Sep 01', 'Sep 15', 'Oct 01', 'Oct 15', 'Nov 01']

# Healthy paddy field (good season)
ndvi_healthy = [0.12, 0.28, 0.45, 0.63, 0.72, 0.78, 0.74, 0.65, 0.52, 0.35]

# Stress-affected paddy field (drought mid-season)
ndvi_stressed = [0.10, 0.22, 0.38, 0.45, 0.42, 0.40, 0.37, 0.30, 0.22, 0.18]

fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(dates, ndvi_healthy, 'g-o', linewidth=2.5, markersize=8,
        label='Healthy Paddy (Kamrup Block A)')
ax.plot(dates, ndvi_stressed, 'r--s', linewidth=2.5, markersize=8,
        label='Stress-Affected Paddy (Kamrup Block B)')

# Add reference lines
ax.axhline(y=0.6, color='gray', linestyle=':', alpha=0.7, label='Moderate-Good threshold (0.6)')
ax.axhline(y=0.4, color='orange', linestyle=':', alpha=0.7, label='Stress warning threshold (0.4)')

# Shade the growing season
ax.axvspan(2, 7, alpha=0.05, color='green', label='Active growing season')

ax.set_xlabel('Date (Kharif 2023)', fontsize=12)
ax.set_ylabel('Mean NDVI', fontsize=12)
ax.set_title('NDVI Time Series — Paddy Crop Monitoring, Kamrup District, Assam\n'
             'Kharif Season 2023', fontsize=13, fontweight='bold')
ax.legend(loc='lower right', fontsize=10)
ax.set_ylim(0, 0.95)
ax.grid(True, alpha=0.3)
plt.xticks(rotation=30)

plt.tight_layout()
plt.savefig("ndvi_time_series_paddy.png", dpi=200, bbox_inches='tight')
plt.show()

The divergence between the two curves — starting from August — clearly shows where drought stress began. This kind of analysis supports early warning systems for district agriculture offices.

---

Complete Code: All Steps in One Script

Here’s the full workflow in a clean, ready-to-run script:

"""
NDVI Calculation with Python — Complete Workflow
Author: Data Science with DEB
Website: https://dibyendudeb.com
Use case: Agricultural crop health monitoring using Sentinel-2 data

Requirements:
    pip install rasterio numpy matplotlib geopandas
"""

import os
import numpy as np
import rasterio
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import warnings
warnings.filterwarnings('ignore')

# ─── CONFIG ────────────────────────────────────────────────────────────────
RED_BAND_PATH = "T44QKF_20231015_B04_10m.tif"   # Sentinel-2 Band 4 (Red)
NIR_BAND_PATH = "T44QKF_20231015_B08_10m.tif"   # Sentinel-2 Band 8 (NIR)
OUTPUT_NDVI   = "ndvi_output.tif"
PIXEL_AREA_HA = 0.01  # 10m x 10m pixel = 0.01 hectares

# ─── LOAD BANDS ────────────────────────────────────────────────────────────
print("Loading satellite bands...")
with rasterio.open(RED_BAND_PATH) as red_src:
    red = red_src.read(1).astype(float)
    profile = red_src.profile

with rasterio.open(NIR_BAND_PATH) as nir_src:
    nir = nir_src.read(1).astype(float)

# ─── CALCULATE NDVI ────────────────────────────────────────────────────────
print("Calculating NDVI...")
denom = nir + red
denom[denom == 0] = np.nan
ndvi = (nir - red) / denom

print(f"NDVI range: {np.nanmin(ndvi):.4f} to {np.nanmax(ndvi):.4f}")
print(f"Mean NDVI:  {np.nanmean(ndvi):.4f}")

# ─── VISUALISE ─────────────────────────────────────────────────────────────
print("Creating NDVI map...")
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(ndvi, cmap='RdYlGn', vmin=-0.3, vmax=0.9)
ax.set_title("NDVI — Vegetation Health Index", fontsize=14, fontweight='bold')
ax.axis('off')
cbar = plt.colorbar(im, ax=ax, fraction=0.03, pad=0.04)
cbar.set_label("NDVI Value", rotation=270, labelpad=20)
plt.tight_layout()
plt.savefig("ndvi_map.png", dpi=200, bbox_inches='tight')
plt.close()

# ─── CLASSIFY ──────────────────────────────────────────────────────────────
print("Classifying NDVI into agricultural categories...")
classified = np.full(ndvi.shape, np.nan)
classified[ndvi < 0.0]                       = 0
classified[(ndvi >= 0.0) & (ndvi < 0.2)]    = 1
classified[(ndvi >= 0.2) & (ndvi < 0.4)]    = 2
classified[(ndvi >= 0.4) & (ndvi < 0.6)]    = 3
classified[(ndvi >= 0.6) & (ndvi < 0.8)]    = 4
classified[ndvi >= 0.8]                      = 5

# ─── AREA STATS ────────────────────────────────────────────────────────────
categories = {0: "Water/Non-veg", 1: "Bare/Sparse",
              2: "Stressed Crop", 3: "Moderate Health",
              4: "Good Health",   5: "Excellent/Dense"}

print("\n" + "=" * 55)
print("NDVI AREA SUMMARY")
print("=" * 55)
total = np.sum(~np.isnan(classified))
for k, v in categories.items():
    n = np.sum(classified == k)
    print(f"{v:<22} {n:>8,} px  {n*PIXEL_AREA_HA:>10,.1f} ha  "
          f"{100*n/total:>6.1f}%")

# ─── SAVE GEOTIFF ──────────────────────────────────────────────────────────
print(f"\nSaving NDVI GeoTIFF to {OUTPUT_NDVI}...")
profile.update({'dtype': rasterio.float32, 'count': 1, 'nodata': -9999})
ndvi_save = ndvi.copy()
ndvi_save[np.isnan(ndvi_save)] = -9999
with rasterio.open(OUTPUT_NDVI, 'w', **profile) as dst:
    dst.write(ndvi_save.astype(rasterio.float32), 1)

print("✓ All outputs saved. NDVI analysis complete.")

---

Frequently Asked Questions

Q: Can I use this code with Landsat data instead of Sentinel-2?
Yes. For Landsat 8/9, the Red band is Band 4 and NIR is Band 5. The formula and code remain the same — just update the file paths.

Q: What if I don’t have downloaded satellite files? Can I compute NDVI online?
Yes — Google Earth Engine is ideal for large-area or time-series analysis. I’ll cover GEE with Python in a future post.

Q: How do I download Sentinel-2 data for free for my study area in India?
Register at Copernicus Open Access Hub and use the map interface to select your area and date. Downloads are free. I’ll write a dedicated post on this.

Q: What’s the difference between NDVI and other indices like EVI, SAVI, NDWI?
NDVI is the most widely used but has limitations over dense canopies or bare soils. SAVI (Soil Adjusted Vegetation Index) is better for sparse vegetation and is common in semi-arid agricultural research. I’ll cover these in the next post in this series.

---

What’s Next in the Agricultural Data Science Series?

This post is Part 1 of my Agricultural Data Science with Python series. Here’s what’s coming:

• Part 2: Downloading Sentinel-2 data from Copernicus Hub using Python
• Part 3: SAVI, EVI, and NDWI — Which vegetation index should you use?
• Part 4: Crop classification using Random Forest and Sentinel-2 data
• Part 5: Time series NDVI analysis for yield forecasting

---

Conclusion

NDVI is one of the most powerful tools available to modern agricultural scientists — and with Python, it’s accessible, free, and fully customisable.

In this post, you learned how to:

• Calculate NDVI from Sentinel-2 satellite data using rasterio and numpy
• Visualise NDVI as a colour-coded map with matplotlib
• Classify NDVI into agricultural health categories
• Compute area statistics for field or district-level crop reporting
• Build an NDVI time series to track crop health through the season
• Save analysis results as a georeferenced GeoTIFF

The same workflow I’ve shown you is used in real crop monitoring projects. It’s not just an academic exercise — this is how modern scinece works.

---

Want More Like This?

I write practical tutorials on data science, Python, and remote sensing for agricultural researchers — one post every week.

Next up in the Agricultural Data Science series:

• How to download free Sentinel-2 data from Copernicus Hub using Python
• SAVI vs NDVI vs EVI — which vegetation index should you use?
• Random Forest for crop classification from satellite data

If this post helped you, the best next step is to subscribe by email so you don’t miss the next one in this series.

No spam. Just one practical tutorial per week.

👉 [Subscribe here — it’s free]