🔧 Fix Urban heat Satellite imagery Urban Heat Mapping · step 19 Notebook By · Updated

LST Scenes Are Not Comparable

On this page

Problem statement

Two Landsat scenes of the same district, five years apart, show the surfaces 4 °C hotter in the later one. It is tempting to call that warming, or the effect of the new development, or the failure of a greening programme. It is almost certainly the weather. Surface temperature at one moment depends on that day's air temperature, sunshine, wind, the soil moisture left by the last rain and the time of overpass. Those differ between any two scenes by more than any plausible change on the ground.

Absolute LST from single scenes is not comparable across dates. Patterns within a scene are.

Quick answer

Compare relative values — each pixel against its own scene's mean or a reference class — or seasonal composites, never raw single-scene temperatures across dates.

import numpy as np, pandas as pd, rasterio, pystac_client, planetary_computer
from rasterio.warp import transform_bounds, reproject, Resampling

AOI = (-3.75, 40.38, -3.63, 40.47)                                # central Madrid
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
scenes = [i for i in cat.search(collections=["landsat-c2-l2"], bbox=list(AOI),
                                datetime="2019-06-01/2023-08-31",
                                query={"eo:cloud_cover": {"lt": 20}}).items()
          if i.properties["platform"] in ("landsat-8", "landsat-9") and i.datetime.month in (6, 7, 8)]
path_row = max({(i.properties["landsat:wrs_path"], i.properties["landsat:wrs_row"]) for i in scenes},
               key=lambda pr: sum((i.properties["landsat:wrs_path"], i.properties["landsat:wrs_row"]) == pr
                                  for i in scenes))
scenes = sorted((i for i in scenes if (i.properties["landsat:wrs_path"], i.properties["landsat:wrs_row"]) == path_row),
                key=lambda i: i.datetime)

def read(item, asset):
    with rasterio.open(planetary_computer.sign(item.assets[asset].href)) as src:
        win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", src.crs, *AOI), src.transform).round_offsets().round_lengths()
        return src.read(1, window=win), src.window_transform(win), src.crs

qa0, tf, crs = read(scenes[0], "qa_pixel")                      # one path/row: one pixel grid for all scenes
wc = np.zeros(qa0.shape, dtype="uint8")
wc_item = next(cat.search(collections=["esa-worldcover"], bbox=list(AOI), datetime="2021").items())
with rasterio.open(planetary_computer.sign(wc_item.assets["map"].href)) as src:
    reproject(rasterio.band(src, 1), wc, dst_transform=tf, dst_crs=crs, resampling=Resampling.mode)
rows = []
for it in scenes:
    try:
        dn, _, _ = read(it, "lwir11")
        qa, _, _ = read(it, "qa_pixel")
    except rasterio.errors.RasterioIOError:                      # an asset the archive cannot serve
        print(f"skipped {it.id}: asset not readable")
        continue
    cloudy = ((qa >> 1) & 1 | (qa >> 3) & 1 | (qa >> 4) & 1).astype(bool)
    t = np.where((dn == 0) | cloudy, np.nan, dn * 0.00341802 + 149.0 - 273.15)
    if np.isnan(t).mean() > 0.1:
        continue
    rows.append({"date": it.datetime.date(), "built": np.nanmedian(t[wc == 50]),
                 "trees": np.nanmedian(t[wc == 10])})
df = pd.DataFrame(rows).set_index("date")
df["gap"] = df.built - df.trees
print(f"{len(df)} clear summer scenes, WRS path/row {path_row}")
print(df.round(1).to_string())
print(f"\nrange of built-up median across scenes {df.built.max() - df.built.min():.1f} °C; "
      f"range of the built–tree gap {df.gap.max() - df.gap.min():.1f} °C")
skipped LC09_L2SP_201032_20220718_02_T1: asset not readable
23 clear summer scenes, WRS path/row ('201', '032')
            built  trees  gap
date                         
2019-07-02   44.8   42.1  2.6
2019-07-18   44.3   41.8  2.5
2019-08-03   44.2   41.9  2.3
2020-06-18   40.5   37.8  2.7
2020-07-04   44.7   41.9  2.8
2020-07-20   45.0   42.5  2.5
2020-08-05   44.3   41.8  2.5
2021-07-07   38.3   35.8  2.5
2021-07-23   40.6   39.0  1.6
2021-08-08   39.9   38.1  1.8
2022-06-24   37.9   35.7  2.2
2022-07-02   45.3   42.9  2.4
2022-07-10   48.2   45.7  2.5
2022-07-26   45.3   43.0  2.3
2022-08-03   44.6   41.4  3.2
2022-08-11   43.7   41.8  1.9
2023-06-27   42.6   40.5  2.1
2023-07-13   45.8   43.5  2.3
2023-07-29   44.2   42.4  1.9
2023-08-06   42.9   41.1  1.8
2023-08-14   42.4   40.8  1.7
2023-08-22   45.9   44.0  1.9
2023-08-30   35.6   33.9  1.6

range of built-up median across scenes 12.6 °C; range of the built–tree gap 1.6 °C

Across 23 clear summer scenes, the median temperature of the same built-up area ranges over 12.6 °C, from 35.6 to 48.2 °C. The difference between built-up surfaces and tree canopy stays between 1.6 and 3.2 °C. (One scene in the archive could not be read and is skipped; with dozens of scenes, losing one changes nothing, but the code has to allow for it.) The pattern is stable; the level is weather.

Never subtract single scenes, normalise, composite

1. Never subtract single scenes from different dates

The difference mixes land change with weather. Unless the weather was identical, the difference is mostly weather.

2. Normalise within each scene

Subtract the scene's median (or the median of a stable reference class), or convert to z-scores. Relative values compare better across dates.

3. Use composites for change

Median composites over whole summers dampen single-day weather; comparing summers still carries each summer's weather.

4. Hold the sensor and time constant

Different satellites pass at different times; Landsat 7's drift makes late scenes early-morning scenes.

5. Relate to air temperature

For absolute comparisons, pair each scene with the air temperature at overpass and model the relationship.

Median built-up and tree-canopy surface temperature in central Madrid on six clear summer Landsat scenes, showing the level moving and the gap stable. built-up median (gap to trees in the note)2 Jul 201944.8 °C · gap 2.67 Jul 202138.3 · gap 2.524 Jun 202237.9 · gap 2.210 Jul 202248.2 · gap 2.522 Aug 202345.9 · gap 1.930 Aug 202335.6 · gap 1.6025°50°Across 23 scenes: level range 12.6 °C, gap range 1.6 °C spatialworkflow.io
The level moves with the weather; the gap between surfaces holds.

A false trend and how much is the air

Example 1 — the false trend

first, last = df.iloc[0], df.iloc[-1]
print(f"{df.index[0]} → {df.index[-1]}: built-up {first.built:.1f} → {last.built:.1f} °C "
      f"({last.built - first.built:+.1f})")
print(f"trees over the same pair: {first.trees:.1f} → {last.trees:.1f} °C ({last.trees - first.trees:+.1f})")
print(f"gap: {first.gap:.1f} → {last.gap:.1f} °C ({last.gap - first.gap:+.1f})")
2019-07-02 → 2023-08-30: built-up 44.8 → 35.6 °C (-9.2)
trees over the same pair: 42.1 → 33.9 °C (-8.2)
gap: 2.6 → 1.6 °C (-1.0)

Picking the first and last scenes and differencing them produces a "cooling" of 9.2 °C in built-up land and 8.2 °C in tree canopy between 2019 and 2023 — neither of which changed much. The last scene was simply a cooler day. The gap between them changes by 1 °C.

Example 2 — how much of the scene level is the air

years = {y: pd.read_csv(f"https://noaa-global-hourly-pds.s3.amazonaws.com/{y}/08221099999.csv",
                         usecols=["DATE", "TMP"], low_memory=False) for y in {d.year for d in df.index}}
air = []
for d in df.index:                                           # Barajas, the hours around the 10:50 UTC overpass
    s = years[d.year]
    s = s[s.DATE.str.startswith(str(d)) & s.DATE.str[11:13].isin(["10", "11"])]
    v = pd.to_numeric(s.TMP.str.split(",").str[0], errors="coerce") / 10
    air.append(v[v.abs() < 90].mean())
df["air"] = air
print(f"correlation across scenes: built-up LST with air temperature {df.built.corr(df.air):.2f}; "
      f"gap with air temperature {df.gap.corr(df.air):.2f}")
correlation across scenes: built-up LST with air temperature 0.81; gap with air temperature -0.03

The scene-to-scene level of surface temperature tracks the air temperature at overpass (r = 0.81), while the built–tree gap does not track it at all (r = −0.03). Normalising removes the shared weather signal and keeps the surface signal.

Example 3 — normalised values

df["built_rel"] = df.built - (df.built + df.trees) / 2
print(f"standard deviation across scenes: raw built-up {df.built.std():.1f} °C, "
      f"relative to the scene {df.built_rel.std():.1f} °C")
standard deviation across scenes: raw built-up 3.0 °C, relative to the scene 0.2 °C

Expressed relative to the scene, built-up surfaces vary by 0.2 °C between dates instead of 3.0 °C — close to the level at which real change could be detected.

Rules for comparing land surface temperature across dates. before any comparison across datesnormalise within each scenebuilt-up sd 3.0 → 0.2 °Ccomposite many scenes per periodweather averages outhold sensor and time constantno drifted Landsat 7relate the level to air temperaturer = 0.81 across scenesdifference two single scenes"−9.2 °C" that was the weatherThe gap between surfaces does not track the weather (r = −0.03) spatialworkflow.io
Before any comparison across dates.

Why weather dominates

Why weather dominates

Midday surface temperature responds within hours to sunshine, wind and air temperature, and within days to soil moisture. Between summer days those vary more than a city's surfaces change in years. A single scene measures the weather of that morning as much as the land.

Why normalising works

Weather raises or lowers all surfaces together, roughly. Subtracting the scene's own median removes that common shift and leaves how each surface stands relative to the rest — the part governed by materials, vegetation and form.

Why composites are not enough on their own

A composite of one summer averages over that summer's weather, not over climate. A hotter summer still raises every pixel. Change detection between summers needs normalisation too, or many summers per period.

The first and last clear scenes of the series compared, raw and as a gap. a false trend from two scenes2 Jul 201930 Aug 2023changebuilt-up median44.8 °C35.6 °C−9.2 °Ctree canopy median42.1 °C33.9 °C−8.2 °Cbuilt − trees2.6 °C1.6 °C−1.0 °CThe last scene was a cooler day, not a cooler city spatialworkflow.io
Two scenes, a nine-degree "change", and almost none of it in the gap between surfaces.

Stable references, soil moisture and phenology

  • Normalisation assumes the reference is stable; if the reference class itself changes, the relative values shift.
  • Soil moisture after rain cools some surfaces more than others; scenes after rain are poor for comparison.
  • Phenology: vegetation in early June and late August differs; compare the same weeks.
  • Mixed paths: scenes from adjacent paths are acquired on different days and different grids.
  • Sensors: Landsat 8 and 9 are consistent; mixing with Landsat 7 after 2017 is not.

FAQ

Can I compare LST between two dates?

Not in absolute terms from single scenes; compare normalised values or composites.

Why did my city get hotter between two scenes?

Most likely the weather on the two days differed.

How do I normalise LST?

Subtract the scene's median or a stable reference class, or convert to z-scores within each scene.

Are summer composites comparable between years?

Better than single scenes, but each still carries its summer's weather.

Does the satellite matter?

Yes; different overpass times give different temperatures, and Landsat 7 drifted.

How do I detect real change in surface temperature?

Normalise within scenes, use many scenes per period, and relate to air temperature.