Upscaling Landsat Thermal Bands with Bilinear

Landsat 8/9 delivers its thermal band (ST_B10) on a 100 m grid while the reflectance bands sit on 30 m, so any per-pixel analysis that mixes surface temperature with a spectral index first needs the thermal band resampled onto the 30 m grid. Use rasterio.warp.reproject with Resampling.bilinear, driving the output from the reflectance band’s transform and shape:

import rasterio
from rasterio.warp import reproject, Resampling
import numpy as np

with rasterio.open("reflectance_B4.tif") as ref:
    dst_profile = ref.profile
    dst_shape = (ref.height, ref.width)

with rasterio.open("thermal_ST_B10.tif") as th:
    thermal = th.read(1, masked=True).astype("float32")
    dst = np.empty(dst_shape, dtype="float32")
    reproject(
        source=thermal, destination=dst,
        src_transform=th.transform, src_crs=th.crs,
        dst_transform=ref.transform, dst_crs=ref.crs,
        resampling=Resampling.bilinear,
    )

This page is a concrete Landsat thermal job inside Advanced Resampling & Upscaling Techniques. If you are still deciding which interpolation kernel a given band deserves, start with Choosing the Right Resampling Method for Sentinel-2 — this page assumes bilinear is already the correct call and focuses on getting the grid alignment and radiometry right.

Why This Arises in Remote Sensing Workflows

Landsat 8 and 9 carry two sensors: the OLI reflectance instrument at 30 m and the TIRS thermal instrument that acquires at 100 m and is delivered resampled. In Collection 2 Level-2 the surface temperature product ST_B10 is distributed on the 100 m grid, so when you want land surface temperature next to NDVI, an albedo estimate, or a burn severity index, the two rasters do not share a pixel geometry. Stacking them naively raises a shape mismatch, and reprojecting with the wrong kernel corrupts the temperature field.

Two things make the thermal case distinct from a routine reflectance upscale. First, the payload is a physical quantity — brightness or surface temperature in Kelvin after scaling — so an interpolation kernel that overshoots produces temperatures that never existed in the scene. Second, the fill value (0 in Collection 2 scaled integers) sits far below any real temperature, so if it leaks into the interpolation stencil it pulls valid edge pixels down by tens of degrees. Both problems are avoidable with bilinear resampling plus explicit nodata handling.

This job belongs to Advanced Resampling & Upscaling Techniques within the broader Satellite Processing Workflows & Index Pipelines section. The same grid-matching idea drives Resampling Sentinel-2 20m Bands to 10m, but the thermal payload changes the radiometric bookkeeping.

How Bilinear Maps 100 m Cells onto the 30 m Grid

Each destination 30 m cell centre falls somewhere inside a 100 m source cell. Bilinear resampling takes the four source cell centres surrounding that point and forms a distance-weighted average, so the output is smooth and every interpolated value stays bounded by its four contributors.

Bilinear weighting from four 100 m thermal cells to a 30 m target cell A destination 30 metre cell centre lies between four 100 metre source cell centres. Bilinear interpolation forms a distance-weighted average of the four surrounding temperatures, keeping the result bounded between them. 100 m source cells 298.2K 299.0K 297.5K 298.8K 30 m output cell weighted avg = 298.4K Result stays bounded by the four inputs

Nearest-neighbour would copy a single 100 m value into nine 30 m cells, producing visible blocking and a staircased temperature field. Cubic convolution samples a 4×4 stencil and can overshoot near sharp thermal boundaries (a lake edge, a fire front), yielding temperatures outside the source range. Bilinear is the middle ground that keeps the surface physically defensible.

Environment & Setup

Package Minimum version Why required
rasterio 1.3.0 rasterio.warp.reproject and Resampling enum
numpy 1.24 Masked arrays and scale/offset arithmetic 
pip install "rasterio>=1.3.0" "numpy>=1.24"

The examples assume Landsat 8/9 Collection 2 Level-2 assets, where ST_B10 is a uint16 grid with fill 0, scale 0.00341802, and additive offset 149.0 to recover Kelvin.

Complete Working Example

The function reads a 30 m reference band, resamples the 100 m thermal band onto its exact grid with bilinear interpolation, applies the Collection 2 scale and offset, and writes an aligned Kelvin GeoTIFF.

import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling

# Collection 2 Level-2 surface temperature scaling
ST_SCALE = 0.00341802
ST_OFFSET = 149.0
FILL = 0

def upscale_thermal_to_30m(
    thermal_path: str,
    reference_path: str,
    out_path: str,
) -> None:
    """
    Resample a Landsat 100 m thermal band onto a 30 m reflectance grid
    with bilinear interpolation and write Kelvin surface temperature.
    """
    # 1. Reference grid defines the destination geometry
    with rasterio.open(reference_path) as ref:
        dst_transform = ref.transform
        dst_crs = ref.crs
        dst_h, dst_w = ref.height, ref.width

    # 2. Read the raw thermal DN as float; keep fill out of the stencil
    with rasterio.open(thermal_path) as th:
        src = th.read(1).astype("float32")
        src[src == FILL] = np.nan  # fill must not average into neighbours
        src_transform = th.transform
        src_crs = th.crs

    # 3. Reproject raw DN with bilinear onto the 30 m grid
    dst_dn = np.full((dst_h, dst_w), np.nan, dtype="float32")
    reproject(
        source=src,
        destination=dst_dn,
        src_transform=src_transform,
        src_crs=src_crs,
        dst_transform=dst_transform,
        dst_crs=dst_crs,
        src_nodata=np.nan,
        dst_nodata=np.nan,
        resampling=Resampling.bilinear,
    )

    # 4. Apply scale + offset AFTER resampling to recover Kelvin
    kelvin = dst_dn * ST_SCALE + ST_OFFSET  # NaN propagates through fill gaps

    # 5. Write an output that shares the 30 m reflectance grid exactly
    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "count": 1,
        "height": dst_h,
        "width": dst_w,
        "crs": dst_crs,
        "transform": dst_transform,
        "nodata": np.nan,
        "compress": "deflate",
        "tiled": True,
    }
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(kelvin.astype("float32"), 1)


if __name__ == "__main__":
    upscale_thermal_to_30m(
        thermal_path="LC09_ST_B10.tif",   # 100 m surface temperature DN
        reference_path="LC09_SR_B4.tif",    # 30 m red reflectance grid
        out_path="LC09_ST_30m_kelvin.tif",
    )

The key sequencing decision is that resampling happens on the raw digital numbers and scaling happens afterward. Because the scale is a linear transform, the order is mathematically equivalent, but interpolating the compact integer field before conversion avoids threading float rounding through the bilinear weights and keeps the fill mask crisp.

Reprojection here is doing two jobs at once even though the two sensors share a CRS and only differ in resolution. reproject warps from the 100 m source transform to the 30 m destination transform, which is a pure upsampling in this case, and it does so by evaluating the bilinear kernel at each destination cell centre. Because the destination grid comes verbatim from the reflectance band — same transform, same width and height, same CRS — the output is guaranteed to align pixel-for-pixel with the 30 m bands, with no half-pixel offset to chase later. That alignment guarantee is the whole point: a temperature raster that is merely “close” to the reflectance grid silently corrupts any per-pixel index that mixes the two.

Variant Patterns

1. Reprojecting to the reflectance grid with WarpedVRT

If you would rather resample lazily and read the aligned thermal band on demand — for example inside a tiled processing loop — wrap the source in a WarpedVRT targeting the 30 m grid instead of materialising a full array.

import rasterio
from rasterio.vrt import WarpedVRT
from rasterio.warp import Resampling

with rasterio.open("LC09_SR_B4.tif") as ref:
    grid = dict(crs=ref.crs, transform=ref.transform, width=ref.width, height=ref.height)

with rasterio.open("LC09_ST_B10.tif") as th:
    with WarpedVRT(th, resampling=Resampling.bilinear, **grid) as vrt:
        aligned = vrt.read(1)  # already on the 30 m grid, read any window from here

2. Stacking the upscaled thermal band with reflectance

Once the thermal band shares the reflectance grid, the arrays stack cleanly for a combined index such as a temperature–vegetation feature space. The alignment guarantee is what makes the stack safe.

import numpy as np
import rasterio

with rasterio.open("LC09_SR_B4.tif") as r, \
     rasterio.open("LC09_SR_B5.tif") as n, \
     rasterio.open("LC09_ST_30m_kelvin.tif") as t:
    red, nir, lst = r.read(1), n.read(1), t.read(1)

ndvi = (nir - red) / (nir + red)          # 30 m
stack = np.stack([ndvi, lst])              # same shape, pixel-aligned

3. Preserving the mask through the pipeline

Carry a boolean validity mask alongside the temperature so downstream steps never treat interpolated fill gaps as real. Any 30 m cell whose bilinear stencil touched a NaN comes back NaN, so a finite check reconstructs the mask exactly.

import numpy as np

valid = np.isfinite(kelvin)   # False where fill leaked or edge had no data
kelvin_filled = np.where(valid, kelvin, -9999.0).astype("float32")

Common Errors

ValueError: operands could not be broadcast together

The thermal and reflectance arrays are still on different grids — you skipped the reproject step or read the thermal band directly. Always drive the destination shape and transform from the 30 m reference band before stacking.

Interpolated coastlines or lake edges read tens of degrees too cold

The fill value bled into the bilinear stencil. Set src[src == FILL] = np.nan before reprojecting and pass src_nodata/dst_nodata, so fill cells are excluded from the weighted average instead of being averaged in as a real temperature.

Output temperatures look uniformly wrong by a constant

The scale and offset were applied to the wrong product or in the wrong order. Confirm you are using the Collection 2 Level-2 surface-temperature scale 0.00341802 and offset 149.0, and apply them after resampling the raw DN, not before.