Imagine standing in a wheat field near Okara, Punjab. From the ground, it is difficult to tell how crop growth differs across hundreds of neighboring fields. From orbit, however, Sentinel-2 satellites observe how sunlight is reflected by vegetation. With a little Python, we can turn these observations into a color-coded NDVI map and group similar pixels into spatial segments. This tutorial shows you how, one executable notebook cell at a time.
🎯 What will you build?
You will produce two maps side by side. On the left, each pixel keeps its NDVI composite value. On the right, a segmentation algorithm groups neighboring pixels and paints each group with its average NDVI. The output is useful for studying vegetation variability, planning field visits and finding areas that deserve a closer look.
Study-area note: The code uses an example agricultural box near Okara, Punjab, Pakistan. The supplied comparison figure later in this post comes from the original tutorial’s Iranian agricultural study area. It illustrates the technique and is not passed off as a Pakistani satellite image. Your Colab run will produce a fresh map for Pakistan.
What are NDVI and crop segmentation?
NDVI stands for Normalized Difference Vegetation Index. Healthy green leaves absorb much of the visible red light used in photosynthesis, while their internal structure reflects a large amount of near-infrared (NIR) light. A satellite measures both. NDVI compares these two observations to estimate vegetation greenness.
NDVI ranges mathematically from −1 to +1. Values around zero often indicate bare or sparsely vegetated surfaces; higher positive values frequently correspond to actively growing plants. NDVI does not, on its own, identify a crop species or directly measure yield. Soil, moisture, canopy structure and growth stage also matter.
Segmentation comes next. Rather than interpreting every 10 m pixel individually, an algorithm groups neighboring pixels into compact, relatively similar regions called superpixels. We calculate the average NDVI for each group, simplifying a busy mosaic of tiny colored squares into larger, interpretable patches.
| Band | What it records | Native resolution | Use |
|---|---|---|---|
| Red / B04 | Visible red reflectance | 10 m | NDVI calculation |
| NIR / B08 | Near-infrared reflectance | 10 m | NDVI calculation |
| SCL | Scene Classification Layer | 20 m | Cloud, shadow and invalid-pixel mask |
We use Sentinel-2 Level-2A surface reflectance, processed for atmospheric effects. The 20 m categorical SCL layer will be aligned to a 10 m grid using nearest-neighbor resampling, but this does not generate new 10 m cloud-classification detail.
Start a fresh Google Colab notebook
Open Google Colab, sign in, select New notebook and rename it Sentinel2_NDVI_Pakistan.ipynb. Each dark code box below belongs in its own code cell. Press the play button beside a cell or hit Shift + Enter. A normal CPU runtime is enough; you do not need a GPU.
Why so many libraries? PySTAC Client finds satellite scenes; Stackstac opens cloud-hosted bands without downloading everything first; Xarray and Dask perform the multidimensional calculations; scikit-image finds segments; SciPy summarizes them; Matplotlib draws the result.
Libraries imported successfully!
Connect Python to the satellite catalog
Downloading satellite images manually would be tedious. A STAC (SpatioTemporal Asset Catalog) API is a searchable index of satellite acquisitions. Element 84’s Earth Search is a public STAC service containing Sentinel-2 data. The following Python call connects our notebook and lets us inspect the collections that are currently published.
Catalog connected: https://earth-search.aws.element84.com/v1 Sentinel-2 collections: - sentinel-2-c1-l2a - sentinel-2-l2a ...Available catalog collections can change. Ellipses denote additional returned entries.
To stay close to the video lesson, we use sentinel-2-c1-l2a (Sentinel-2 Collection 1 Level-2A). The separately named sentinel-2-l2a is another Earth Search collection; if you switch collections, inspect its asset names and metadata before proceeding.
Choose an agricultural study area in Pakistan
Our region of interest (ROI) is described with a rectangle called a bounding box. Its four coordinates are ordered [west longitude, south latitude, east longitude, north latitude]. We will use a small sample rectangle near Okara, Punjab as our Pakistan adaptation. A beginner should use a small rectangle to prevent very large downloads.
ROI (W, S, E, N): [73.43, 30.77, 73.47, 30.8] Year: 2025 Target coordinate system: EPSG:32643
ROI. If the new location lies outside UTM zone 43N, choose the correct local projected CRS too. This is an example box, not a verified farm parcel.Search for all matching Sentinel-2 scenes in 2025
We are now asking Earth Search: Which Sentinel-2 Collection 1 Level-2A scenes overlap my Pakistani study area between 1 January and 31 December 2025? We keep the year broad and later remove clouds at the pixel level. The search returns scene metadata and asset URLs, not yet the full image data.
Items found: [a positive number] First scene: [a Sentinel-2 scene identifier] Native CRS: [an EPSG code] Available assets: ['red', 'nir', 'scl', ...]These placeholders are intentional, not invented results. The Iranian ROI used in the video returned 175 scenes for 2025. Your Pakistan result will have its own count.
Each item represents an acquisition. Its assets are the actual band files and supplementary resources. The keys red, nir and scl are the three we need. Run the last print statement before continuing and confirm those keys appear.
Stack the red, NIR and SCL bands into one data cube
Imagine a pile of 2025 satellite photographs placed on top of one another. That is essentially a spatiotemporal data cube: time along one axis, band names along a second axis, and rows/columns (y and x) across the ground. Stackstac creates such a structure and uses Dask lazy loading to avoid downloading the entire year immediately.
Data cube dimensions: {
'time': [number of scenes], 'band': 3,
'y': [pixel rows], 'x': [pixel columns]
}
Selected bands: ['red', 'nir', 'scl']
Why resolution=10? The red and near-infrared bands have 10 m native spatial resolution, meaning one working pixel represents about 10 × 10 m on the ground. We set EPSG=32643 so that 10 means meters, rather than 10 degrees of latitude. rescale=True asks Stackstac to apply any scale/offset supplied in raster metadata. Resampling.nearest protects discrete SCL classes when that 20 m layer is aligned to the 10 m output grid.
We next convert the band axis into individually named Xarray variables; this makes expressions such as ds["nir"] easy to read.
Dataset variables: ['red', 'nir', 'scl']
Dataset dimensions: {
'time': [scenes], 'y': [rows], 'x': [columns]
}
.compute() later.Remove clouds, shadows and unusable pixels
A cloud may look beautifully white in an optical image, but it is bad news for crop monitoring—it hides the ground. Sentinel-2’s Scene Classification Layer (SCL) labels pixels into predefined categories. We use it like a transparent stencil: mark contaminated pixels as missing and keep the rest for vegetation analysis.
| SCL value | What it means | Our decision |
|---|---|---|
| 0, 1 | No data; saturated/defective | Exclude |
| 3 | Cloud shadow | Exclude |
| 8, 9 | Medium- and high-probability cloud | Exclude |
| 10, 11 | Thin cirrus; snow/ice | Exclude |
| 2, 4, 5, 6, 7 | Dark surface; vegetation; soil; water; unclassified | Keep for this general land-cover example |
After: those observations become
NaN (missing), while accepted pixels keep their measured values. The computation has not been downloaded yet.Calculate NDVI for every clear observation
Now comes the vegetation index itself. For every clear pixel in every Sentinel-2 acquisition, subtract the red reflectance from NIR and divide by their sum. We explicitly avoid division by zero, because invalid values can occur near imagery boundaries.
Suppose a vegetated pixel has NIR = 0.62 and Red = 0.18. Its NDVI is (0.62 − 0.18) / (0.62 + 0.18) = 0.55.
For a less vegetated example with NIR = 0.22 and Red = 0.18, NDVI = 0.10. These are teaching examples, not measured Okara values.
At this stage we have a stack of NDVI maps, one per acquisition. A cloudy observation is missing only where the SCL mask identified a problem; clear observations from other dates can still contribute to the year-long summary.
Build monthly maximum-NDVI composites and a yearly mean
A single clear image can miss a crop’s most active stage. Instead, we aggregate across time. First we group all clear observations into 12 calendar months, selecting the maximum NDVI at each pixel within each month. Then we calculate the mean of those 12 monthly maxima, skipping missing months. This is the workflow taught in the video.
Monthly composites: 12
Annual image dimensions: ('y', 'x')
| January clear observations (example) | Monthly maximum | Contribution to annual image |
|---|---|---|
| 0.32, 0.41, 0.36 | 0.41 | January contributes 0.41 |
| February: 0.44, 0.49 | 0.49 | February contributes 0.49 |
| March: all cloudy | NaN | March is skipped at this pixel |
Imagine doing this for every month: the annual composite is the average of the available monthly peak NDVI values, rather than a simple average of all satellite scenes. This helps prevent months with many acquisitions dominating the result, but it can favor peaks and may conceal periods of poor vegetation. For genuine crop monitoring, compare multiple temporal statistics and the crop calendar.
1ME rather than 1Y? 1ME creates one end-of-month group per month. An annual resampling rule would collapse the whole year to a single value and would not reproduce the monthly-max-then-mean method.Download the computed annual NDVI pixels and inspect them
So far, the large calculations have been lazily planned. Now .compute() actually requests the relevant raster blocks from cloud storage and combines them. Dask’s progress bar shows how the job is proceeding. Loading time depends on scene count, internet access, cloud cover, image size and Colab resources.
[########################################] 100% Completed Pixels with valid annual NDVI: [positive integer] NDVI min: [value between -1 and +1] NDVI max: [value between -1 and +1]Progress-bar design and numerical results vary by environment and actual imagery.
Let us draw the first genuine map before simplifying anything. The colormap RdYlGn runs from red through yellow to green; it is familiar in vegetation analysis and easy to explain to non-specialists.
Turn thousands of pixels into SLIC superpixels
This is where the tutorial gets interesting. A normal NDVI map gives every pixel its own color. But agricultural analysis often benefits from local regions: patches of nearby pixels that have similar vegetation signals. SLIC (Simple Linear Iterative Clustering) organizes an image into spatially coherent segments using both pixel values and locations.
Think of a mosaic artist grouping small green and yellow tiles into larger shapes. n_segments=250 requests approximately 250 segments, not an exact guaranteed count; compactness=0.1 controls how strongly region shape is favored over pixel-value similarity. These are experiment settings from the source lesson, not universally optimal parameters.
Segments actually created: [approximately 250] Label-image shape: ([rows], [columns])The exact count and segment geometry depend on valid-pixel coverage and the imagery. Masked pixels have label 0.
n_segments from 250 to 100 and then 500. Fewer segments generally lead to bigger zones; more segments preserve finer image variation. Adjust compactness too, since increasing it encourages more regularly shaped superpixels. Compare the map against actual field edges in a high-resolution basemap.Assign the average NDVI to every segment
At present segment_labels contains only region IDs: 1, 2, 3, and so on. These integers tell us which pixels belong together, but not how green they are. We therefore calculate the mean NDVI in each labeled region and replace every member pixel with that region's mean.
For example, suppose segment 18 contains four NDVI values—0.44, 0.47, 0.52 and 0.49. Its mean is 0.48, so the four pixels will all be colored according to 0.48 in the segmented display. This smooths within-region variation while retaining boundaries between groups.
Output dimensions: ('y', 'x')
Example segment means: [five computed NDVI values]
The numeric array is computed from real imagery when you run the notebook; it is not pre-filled with fabricated results.
Why restore coordinates? A bare NumPy array knows its row and column positions but not where those pixels belong on Earth. By keeping the x and y coordinates from annual_ndvi, our output can be plotted with meaningful projected eastings and northings.
Make the final, publication-style two-panel NDVI map
We have reached the final deliverable. A side-by-side view is much more informative than two isolated plots: the first panel shows continuous pixel variation, while the second shows the same vegetation information summarized inside SLIC regions. We keep an identical NDVI color scale in both panels so the colors are directly comparable.
Saved: ndvi_segmentation_pakistan.pngThe notebook displays both maps. The file appears in Colab's file browser (folder icon on the left). At 360 dpi, it is suitable for high-quality on-screen or print use.
Troubleshooting: solve these common problems
❓ The STAC search says “0 items found.”
Confirm the bounding box order is west, south, east, north; inspect the area on a map; confirm the year and collection name; and try a slightly larger region. Search results may also vary when a data provider updates its archive.
❓ Colab is taking too long, freezes, or runs out of memory.
Begin with an even smaller ROI—about 1–3 km across—and consider a three-month period for a first practice run. Processing time scales with scene count and pixel dimensions. Only expand once the short experiment works.
❓ A raster request fails with 403, 404 or a connection error.
Restart the stalled cell, check internet connectivity and inspect a printed asset URL from first_item.assets["red"].href. Catalog accessibility and object storage accessibility are different things; provider-side outages and expired or changed asset addresses can occur.
❓ The NDVI plot is blank or has large empty patches.
Check the printed count of valid pixels. Too few clear images, overly strict cloud masking or a study area dominated by missing data can cause blanks. Widen the date range, inspect the SCL distribution and verify that your scene assets truly contain the desired bands.
❓ Why do I see cloud artifacts even after masking?
The SCL product is an automated classification, not perfect ground truth. Thin cloud and haze may survive. For high-stakes studies, inspect RGB imagery, consider cloud probability and compare multiple screening strategies rather than hiding warnings.
❓ SLIC throws an error about NaN values or channels.
Keep the valid mask, replace missing image inputs with zeros only for algorithm compatibility, and pass channel_axis=None because NDVI is a single-channel 2D image. The mask prevents missing pixels being treated as genuine crop observations.
❓ My map is less detailed than the figure.
Different geographies naturally contain different field shapes. Also compare the ROI extent, crop season, sensor resolution, n_segments and compactness. SLIC returns approximately the number of segments requested, and the optimal value depends on the actual landscape.
Practical uses in Pakistani agriculture
Once your notebook is working, apply the technique to an agricultural region relevant to your own work. Possible research and extension applications include:
For a better real-world agriculture project, repeat the workflow for different months or seasons, compare it with local rainfall/irrigation information and add ground observations. Wheat, rice, sugarcane and other crops have different growth calendars; a full-year maximum-NDVI average can hide those distinctions.
Original video walkthrough
This post develops a Pakistan-oriented, beginner-friendly implementation from the supplied English captions and figure for “PySTAC Tutorial-13 Sentinel-2 Crop Segmentation using NDVI in Python”, presented by remote-sensing data scientist Amirhossein Aghaee. The video uses a study area in Iran; we have adapted the notebook steps and added safe missing-data handling, explanations and a new Pakistan ROI.
🚀 Your next satellite map starts with one small rectangle
Open Colab, run the cells in order, inspect your NDVI map and try changing n_segments from 100 to 500. You do not need expensive proprietary GIS software to begin exploring satellite-based agriculture. You need a study area, a clear research question and a willingness to test what the colors really mean.
If you found this tutorial useful, share it with a student, a remote-sensing enthusiast or an agricultural researcher working in Pakistan.
Post a Comment