Overlaying a 13-band Sentinel-2 image on a classification raster fails for one of three reasons: the two files use different coordinate reference systems, they have different pixel sizes, or their extents do not line up. Sentinel-2 alone mixes 10, 20 and 60 metre bands, so the resampling step is unavoidable. This guide works through the whole alignment in rasterio.
Setting up#
python -m pip install rasterio numpy matplotlib
Rasterio bundles GDAL, so this is usually all you need. On Windows, conda-forge is the smoother route if pip struggles with the binary dependencies.
Step 1: inspect both files#
import rasterio
def describe(path):
with rasterio.open(path) as src:
print(path)
print(" CRS :", src.crs)
print(" size :", src.width, "x", src.height)
print(" bands :", src.count, src.dtypes[0])
print(" resolution:", src.res)
print(" bounds :", src.bounds)
print(" nodata :", src.nodata)
describe("sentinel2_13band.tif")
describe("classification.tif")
Compare four things: the CRS, the resolution, the bounds and the nodata value. Any difference in the first three means you cannot simply read both into arrays and index them together — pixel (0, 0) will not be the same place on the ground.
Step 2: the Sentinel-2 resolution problem#
| Bands | Resolution | Typical use |
|---|---|---|
| 2, 3, 4, 8 | 10 m | blue, green, red, near infrared |
| 5, 6, 7, 8A, 11, 12 | 20 m | red edge, shortwave infrared |
| 1, 9, 10 | 60 m | atmospheric correction |
A single file holding all thirteen has already been resampled to one grid by whoever produced it. If you are reading the original JP2 files, you must resample yourself before stacking:
import numpy as np
import rasterio
from rasterio.enums import Resampling
TARGET = "B02.jp2" # a 10 m band defines the grid
with rasterio.open(TARGET) as ref:
height, width = ref.height, ref.width
profile = ref.profile
def read_resampled(path, height, width):
with rasterio.open(path) as src:
return src.read(
1,
out_shape=(height, width),
resampling=Resampling.bilinear,
)
bands = [read_resampled(p, height, width) for p in band_paths]
stack = np.stack(bands, axis=0)
print(stack.shape) # (13, height, width)
Step 3: reproject the classification onto the image grid#
import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling
def align_to(source_path, reference_path):
"""Return the source raster resampled onto the reference grid."""
with rasterio.open(reference_path) as ref:
dst = np.zeros((ref.height, ref.width), dtype=np.uint8)
with rasterio.open(source_path) as src:
reproject(
source=rasterio.band(src, 1),
destination=dst,
src_transform=src.transform,
src_crs=src.crs,
dst_transform=ref.transform,
dst_crs=ref.crs,
resampling=Resampling.nearest, # categories!
src_nodata=src.nodata,
dst_nodata=0,
)
return dst
classes = align_to("classification.tif", "sentinel2_13band.tif")
print(classes.shape, np.unique(classes))
Step 4: read the image and check the shapes#
with rasterio.open("sentinel2_13band.tif") as src:
image = src.read().astype("float32") # (13, height, width)
profile = src.profile
transform = src.transform
print(image.shape, classes.shape)
assert image.shape[1:] == classes.shape, "grids still do not match"
That assertion is worth keeping. Everything downstream assumes the two arrays index the same ground positions, and a silent mismatch produces plausible-looking nonsense.
Step 5: mask by class#
import numpy as np
FOREST = 3
mask = classes == FOREST
print("forest pixels:", mask.sum())
# Mean reflectance per band, forest only
for band_index in range(image.shape[0]):
values = image[band_index][mask]
values = values[np.isfinite(values)]
print(f"band {band_index + 1:2d}: mean {values.mean():8.1f}")
A vegetation index for the same area, using bands 8 and 4 (index 7 and 3 in a zero-based array):
nir = image[7]
red = image[3]
with np.errstate(divide="ignore", invalid="ignore"):
ndvi = (nir - red) / (nir + red)
ndvi = np.where(np.isfinite(ndvi), ndvi, np.nan)
for class_id in np.unique(classes):
if class_id == 0:
continue
values = ndvi[classes == class_id]
values = values[np.isfinite(values)]
if values.size:
print(f"class {class_id}: NDVI median {np.median(values):.3f}")
The errstate context suppresses the division warnings you get where both bands are zero — over water, or outside the swath. Replacing those with NaN rather than leaving infinities is what stops them poisoning later statistics.
Step 6: plot the overlay#
import matplotlib.pyplot as plt
import numpy as np
def stretch(band, low=2, high=98):
"""Percentile stretch for display only."""
lo, hi = np.nanpercentile(band, [low, high])
return np.clip((band - lo) / (hi - lo), 0, 1)
rgb = np.dstack([stretch(image[3]), stretch(image[2]), stretch(image[1])])
fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(rgb)
ax.imshow(np.ma.masked_where(classes == 0, classes),
cmap="tab10", alpha=0.45, interpolation="nearest")
ax.set_title("Sentinel-2 true colour with classification overlay")
ax.axis("off")
plt.tight_layout()
plt.show()
Two details matter for the display. The percentile stretch is needed because raw reflectance values are not in the 0 to 1 range imshow expects. And masked_where hides the nodata class so the underlying image shows through rather than being covered by a solid block.
Step 7: save the aligned result#
profile.update(count=1, dtype="uint8", nodata=0, compress="lzw")
with rasterio.open("classification_aligned.tif", "w", **profile) as dst:
dst.write(classes, 1)
Writing the aligned classification once saves repeating the reprojection every time you run an analysis, and gives you something you can check in a GIS viewer.
Questions people ask#
Which raster should define the grid?
Usually the imagery, since resampling categorical data with nearest-neighbour loses less than resampling reflectance. If the classification is much coarser, consider whether the analysis is meaningful at the finer resolution at all.
What if the two files barely overlap?
Compute the intersection of their bounds and window both to it. Reprojecting into a grid the source barely covers gives you a mostly empty array and no warning.
Should I use rioxarray instead?
It is a pleasant layer over rasterio that keeps coordinates and metadata attached, and its reproject_match does this whole alignment in one call. Worth adopting if you are already using xarray.
Why are my reflectance values in the thousands?
Sentinel-2 products store reflectance scaled by 10000 as integers. Divide by the scale factor before computing indices, or at least be consistent about it.
Where to go next#
- Exploring np.stack() — building the band stack correctly.
- Data visualization — matplotlib beyond a single imshow.
- Python and data science — NumPy fundamentals.
