Thermal Sharpening of CLMS Land Surface Temperature with Landsat

//VERSION=3

/**
 * Thermal sharpening of CLMS LST (3 km hourly V3) with Landsat thermal imagery.
 * By András Zlinszky, Sinergise (azlinszky.bsky.social) and Claude Code
 * Based on the Landsat Land Surface Temperature custom script by Mohor Gartner
 * https://custom-scripts.sentinel-hub.com/custom-scripts/landsat-8/land_surface_temperature_mapping/
 * Data fusion of three sources; see README.md for the method, caveats and setup.
 *
 * Quick setup (Copernicus Browser -> "use additional datasets (advanced)"):
 *   - Primary dataset      Landsat 8-9 L1        -> alias  LANDSAT-OT-L1
 *   - Additional dataset 1 CLMS LST 3km hourly V3 -> alias  CLMS_CAL
 *   - Additional dataset 2 CLMS LST 3km hourly V3 -> alias  CLMS_TGT
 *   Customize timespans: LANDSAT-OT-L1 + CLMS_CAL = calibration date (clear-sky
 *   Landsat overpass); CLMS_TGT = the target date you want to sharpen. Set CLMS_CAL
 *   to the hourly slot CONTAINING the overpass (~09:30 UTC over central Europe, so
 *   the 09:00 slot) - a time gap between Landsat and CLMS_CAL corrupts the residual.
 *
 * Per pixel, in Kelvin:
 *   residual      = LST_landsat_cal - LST_clms_cal
 *   LST_sharpened = LST_clms_target + residual + biasOffset
 *
 * Optional scaling of the CLMS change (changeScale = s, default 1 = the formula above):
 *   LST_sharpened = LST_landsat_cal + s * (LST_clms_target - LST_clms_cal)
 *                   + (1 - s) * deltaTMean + biasOffset
 */

//// ---- USER OPTIONS ---------------------------------------------------------

// choose which Landsat TIRS band drives brightness temperature (B10 or B11)
var band = "B10";

// scene-wide bias offset (see header). Units match the LST (Kelvin degrees).
// 0 = pure additive Landsat detail; set to mean(CLMS_cal)-mean(Landsat_cal)
// to anchor the sharpened scene-mean to CLMS.
var biasOffset = 0;

// scale s of the 3 km CLMS change between the calibration and the target date.
// 1 = additive sharpening, the full CLMS change is applied (default). Values below 1
// damp the change: use the median slope of Landsat on CLMS over several calibration
// dates (companion notebook, Section 12.5); over Budapest any value 0.25-0.5 worked.
var changeScale = 1;

// only used when changeScale != 1: mean(CLMS_TGT) - mean(CLMS_CAL) over your area,
// in Kelvin (the same number in degC). Keeps the sharpened scene mean on the CLMS level.
var deltaTMean = 0;

// physical sanity range for the Landsat-derived calibration LST (degC).
// Pixels whose derived LST falls outside this band are masked out.
var minValidC = -20;
var maxValidC = 60;

// surface roughness term, added to the emissivity of mixed soil/vegetation pixels.
// 0.005 follows Sobrino et al. (2008) and Avdan & Jovanovska (2016).
var C = 0.005;

var applyCloudMask = true;

// color-map range for the RGBA "default" output, in Kelvin. The defaults
// 240-330 K (-33 to +57 degC) reproduce the CLMS LST 3km hourly V3 palette
// exactly. Use the same values in landsat_lst_reference.js when comparing the two layers.
// degC = K - 273.15
var paletteMinK = 240;
var paletteMaxK = 330;

//// ---- LANDSAT LST CONSTANTS ------------------------------------------------

// NDVI thresholds for bare soil / full vegetation
var NDVIs = 0.2;
var NDVIv = 0.8;

// emissivity
var waterE = 0.991;
var soilE = 0.966;
var vegetationE = 0.973;

// central/mean wavelength in meters, B10 or B11
var bCent = band == "B10" ? 0.000010895 : 0.000012005;

// rho = h*c/sigma = PlanckC*velocityLight/BoltzmannC
var rho = 0.01438; // m K

//// ---- CLMS DECODE ----------------------------------------------------------

// raw CLMS LST -> Kelvin
var CLMS_FACTOR = 1 / 100;
var CLMS_OFFSET = 273.15;

//// ---- VISUALIZATION (CLMS LST 3km hourly V3 magma-style ramp) --------------

// the 10 magma colors of the CLMS ramp, spread evenly over [paletteMinK, paletteMaxK]
const RampColors = [
  [0, 0, 4],
  [27, 12, 65],
  [76, 12, 107],
  [120, 28, 109],
  [165, 45, 96],
  [206, 68, 70],
  [237, 105, 37],
  [251, 154, 7],
  [247, 208, 60],
  [252, 255, 164],
];
const ColorBar = RampColors.map(function (rgb, i) {
  var t = i / (RampColors.length - 1);
  return [paletteMinK + t * (paletteMaxK - paletteMinK), rgb];
});
const visualizer = new ColorRampVisualizer(ColorBar);

//// ---- SETUP ----------------------------------------------------------------

function setup() {
  return {
    input: [
      {
        datasource: "LANDSAT-OT-L1",
        bands: ["B03", "B04", "B05", "B10", "B11", "BQA", "dataMask"],
      },
      { datasource: "CLMS_CAL", bands: ["LST", "dataMask"] },
      { datasource: "CLMS_TGT", bands: ["LST", "dataMask"] },
    ],
    output: [
      { id: "default", bands: 4, sampleType: "UINT8" }, // RGBA color map
      { id: "index", bands: 1, sampleType: "FLOAT32" }, // sharpened LST (K)
      { id: "dataMask", bands: 1 }, // validity mask
    ],
    mosaicking: "SIMPLE", // one clear scene per layer; access samples.X[0]
  };
}

//// ---- LANDSAT EMISSIVITY (unchanged from source script) --------------------

function LSEcalc(NDVI, Pv) {
  var LSE;
  if (NDVI < 0) {
    // water
    LSE = waterE;
  } else if (NDVI < NDVIs) {
    // soil
    LSE = soilE;
  } else if (NDVI > NDVIv) {
    // vegetation
    LSE = vegetationE;
  } else {
    // mixtures of vegetation and soil
    LSE = vegetationE * Pv + soilE * (1 - Pv) + C;
  }
  return LSE;
}

// Landsat single-scene LST in KELVIN, or null if the pixel is invalid.
function landsatLST(s) {
  var Bi = band == "B10" ? s.B10 : s.B11;
  var B03i = s.B03;
  var B04i = s.B04;
  var B05i = s.B05;

  // raw-band guard: screens error scenes / fill (same as source script)
  if (!(Bi > 173 && Bi < 65000 && B03i > 0 && B04i > 0 && B05i > 0)) {
    return null;
  }

  // 1 Kelvin to C (brightness temperature)
  var btC = Bi - 273.15;
  // 2 NDVI
  var NDVI = (B05i - B04i) / (B05i + B04i);
  // 3 Pv - proportional vegetation
  var Pv = Math.pow((NDVI - NDVIs) / (NDVIv - NDVIs), 2);
  // 4 LSE - land surface emissivity
  var LSE = LSEcalc(NDVI, Pv);
  // 5 LST (degC): the single-channel land surface temperature calculation based on
  // (Artis & Carnahan, 1982), converting brightness temperature (btC) to Land
  // Surface Temperature (lstC) using the emissivity (LSE) and the central wavelength
  // (bCent) and rho (Planck's constant times the speed of light divided by
  // Boltzmann's constant).
  var lstC = btC / (1 + ((bCent * btC) / rho) * Math.log(LSE));

  // physical sanity range
  if (lstC < minValidC || lstC > maxValidC) {
    return null;
  }

  // return Kelvin
  return lstC + 273.15;
}

//// ---- EVALUATE PIXEL -------------------------------------------------------

function invalid() {
  return {
    default: [0, 0, 0, 0],
    index: [NaN],
    dataMask: [0],
  };
}

function evaluatePixel(samples) {
  // note: LANDSAT-OT-L1 has hyphens, so it must be accessed via bracket notation
  var landsat = samples["LANDSAT-OT-L1"][0];
  var clmsCal = samples.CLMS_CAL[0];
  var clmsTgt = samples.CLMS_TGT[0];

  if (landsat.dataMask === 0 || clmsCal.dataMask === 0 || clmsTgt.dataMask === 0) {
    return invalid();
  }

  if (applyCloudMask) {
    var qa = decodeL8C2Qa(landsat.BQA);
    if (qa.fill || qa.cloud || qa.dilatedCloud || qa.cirrus || qa.cloudShadow) {
      return invalid();
    }
  }

  // Landsat calibration LST (Kelvin), with raw-band + physical-range guards
  var lstLandsatCal = landsatLST(landsat);
  if (lstLandsatCal === null) {
    return invalid();
  }

  // CLMS LST (Kelvin)
  var lstClmsCal = clmsCal.LST * CLMS_FACTOR + CLMS_OFFSET;
  var lstClmsTgt = clmsTgt.LST * CLMS_FACTOR + CLMS_OFFSET;

  // additive residual sharpening: here we calculate the difference between Landsat
  // and the simultaneous CLMS, and apply it to the target CLMS together with the
  // bias offset.
  var residual = lstLandsatCal - lstClmsCal;
  var lstSharpened = lstClmsTgt + residual + biasOffset;

  // optional: damp the local CLMS change by changeScale, but keep the full
  // scene-mean change deltaTMean, so the scene mean stays on the CLMS level
  if (changeScale !== 1) {
    var clmsChange = lstClmsTgt - lstClmsCal;
    lstSharpened -= (1 - changeScale) * (clmsChange - deltaTMean);
  }

  var rgb = visualizer.process(lstSharpened);
  return {
    default: rgb.concat(255),
    index: [lstSharpened],
    dataMask: [1],
  };
}

//VERSION=3

/**
 * Native Landsat LST reference layer, for validating the sharpened CLMS product.
 *
 * Deliberately identical to the Landsat half of script.js in this folder: same band,
 * same emissivity model, same constants, same guards, same CLMS magma color ramp and
 * same output ids. Run it on a Landsat overpass date and compare side by side with the
 * sharpened layer for that same date - any visible difference is the sharpening, not a
 * difference in method.
 *
 * Single dataset: Landsat 8-9 L1. No data fusion, no aliases to set.
 */

//// ---- USER OPTIONS ---------------------------------------------------------

// choose which Landsat TIRS band drives brightness temperature (B10 or B11)
var band = "B10";

// keep this at 0 to show the raw Landsat retrieval. Set it to the SAME value used in
// script.js if you want the two layers on a common absolute scale: Landsat L1 B10 is
// top-of-atmosphere brightness temperature (no atmospheric correction), so it reads
// systematically colder than the atmospherically corrected CLMS LST.
var biasOffset = 0;

// physical sanity range for the derived LST (degC)
var minValidC = -20;
var maxValidC = 60;

// surface roughness term, added to the emissivity of mixed soil/vegetation pixels.
// 0.005 follows Sobrino et al. (2008) and Avdan & Jovanovska (2016). Note that the
// Landsat-8 LST Mapping script this retrieval is ported from uses 0.009 - set it to
// that value if you need results consistent with that script (see README).
// Keep this identical to the value in script.js so the two layers stay comparable.
var C = 0.005;

// mask pixels the Landsat Collection 2 QA band flags as cloud-affected
var applyCloudMask = true;

// color-map range for the RGBA "default" output, in Kelvin. The defaults
// 240-330 K (-33 to +57 degC) reproduce the CLMS LST 3km hourly V3 palette
// exactly. Narrowing the range (e.g. 295-325 K, +22 to +52 degC, for a summer
// morning) makes intra-urban contrast easier to read - but keep these identical
// to the values in script.js, or the two layers cannot be compared by color.
// degC = K - 273.15
var paletteMinK = 240;
var paletteMaxK = 330;

//// ---- LANDSAT LST CONSTANTS ------------------------------------------------

// NDVI thresholds for bare soil / full vegetation
var NDVIs = 0.2;
var NDVIv = 0.8;

// emissivity
var waterE = 0.991;
var soilE = 0.966;
var vegetationE = 0.973;

// central/mean wavelength in meters, B10 or B11
var bCent = band == "B10" ? 0.000010895 : 0.000012005;

// rho = h*c/sigma = PlanckC*velocityLight/BoltzmannC
var rho = 0.01438; // m K

//// ---- VISUALIZATION (CLMS LST 3km hourly V3 magma-style ramp) --------------

// the 10 magma colors of the CLMS ramp, spread evenly over [paletteMinK, paletteMaxK]
const RampColors = [
  [0, 0, 4],
  [27, 12, 65],
  [76, 12, 107],
  [120, 28, 109],
  [165, 45, 96],
  [206, 68, 70],
  [237, 105, 37],
  [251, 154, 7],
  [247, 208, 60],
  [252, 255, 164],
];
const ColorBar = RampColors.map(function (rgb, i) {
  var t = i / (RampColors.length - 1);
  return [paletteMinK + t * (paletteMaxK - paletteMinK), rgb];
});
const visualizer = new ColorRampVisualizer(ColorBar);

//// ---- SETUP ----------------------------------------------------------------

function setup() {
  return {
    input: [{ bands: ["B03", "B04", "B05", "B10", "B11", "BQA", "dataMask"] }],
    output: [
      { id: "default", bands: 4, sampleType: "UINT8" }, // RGBA color map
      { id: "index", bands: 1, sampleType: "FLOAT32" }, // LST (K)
      { id: "browserStats", bands: 1, sampleType: "FLOAT32" }, // LST (K) for stats
      { id: "dataMask", bands: 1 }, // validity mask
    ],
  };
}

//// ---- LANDSAT EMISSIVITY (unchanged from source script) --------------------

function LSEcalc(NDVI, Pv) {
  var LSE;
  if (NDVI < 0) {
    // water
    LSE = waterE;
  } else if (NDVI < NDVIs) {
    // soil
    LSE = soilE;
  } else if (NDVI > NDVIv) {
    // vegetation
    LSE = vegetationE;
  } else {
    // mixtures of vegetation and soil
    LSE = vegetationE * Pv + soilE * (1 - Pv) + C;
  }
  return LSE;
}

// Landsat single-scene LST in KELVIN, or null if the pixel is invalid.
function landsatLST(s) {
  var Bi = band == "B10" ? s.B10 : s.B11;
  var B03i = s.B03;
  var B04i = s.B04;
  var B05i = s.B05;

  // raw-band guard: screens error scenes / fill (same as source script)
  if (!(Bi > 173 && Bi < 65000 && B03i > 0 && B04i > 0 && B05i > 0)) {
    return null;
  }

  // 1 Kelvin to C (brightness temperature)
  var btC = Bi - 273.15;
  // 2 NDVI
  var NDVI = (B05i - B04i) / (B05i + B04i);
  // 3 Pv - proportional vegetation
  var Pv = Math.pow((NDVI - NDVIs) / (NDVIv - NDVIs), 2);
  // 4 LSE - land surface emissivity
  var LSE = LSEcalc(NDVI, Pv);
  // 5 LST (degC)
  var lstC = btC / (1 + ((bCent * btC) / rho) * Math.log(LSE));

  // physical sanity range
  if (lstC < minValidC || lstC > maxValidC) {
    return null;
  }

  // return Kelvin
  return lstC + 273.15;
}

//// ---- EVALUATE PIXEL -------------------------------------------------------

function invalid() {
  return {
    default: [0, 0, 0, 0],
    index: [NaN],
    browserStats: [NaN],
    dataMask: [0],
  };
}

function evaluatePixel(sample) {
  if (sample.dataMask === 0) {
    return invalid();
  }

  if (applyCloudMask) {
    var qa = decodeL8C2Qa(sample.BQA);
    if (qa.fill || qa.cloud || qa.dilatedCloud || qa.cirrus || qa.cloudShadow) {
      return invalid();
    }
  }

  var lst = landsatLST(sample);
  if (lst === null) {
    return invalid();
  }

  lst = lst + biasOffset;

  var rgb = visualizer.process(lst);
  return {
    default: rgb.concat(255),
    index: [lst],
    browserStats: [lst],
    dataMask: [1],
  };
}

//VERSION=3

/**
 * Landsat - CLMS LST residual, for measuring the biasOffset that script.js needs.
 *
 * Data fusion of two sources. Outputs the per-pixel difference
 *   residual = LST_landsat_cal - LST_clms_cal
 * in Kelvin degrees. The scene-wide mean of that field is the bias script.js needs:
 *
 *   biasOffset = -mean(residual)
 *
 * Note the sign. Landsat L1 B10 is top-of-atmosphere brightness temperature with no
 * atmospheric correction, so it reads colder than CLMS: the mean residual is negative
 * and the biasOffset that cancels it is positive.
 *
 * KNOWN COPERNICUS BROWSER LIMITATION: the statistics / histogram panel does not work
 * on data fusion layers. It fails with "Dataset with id: <alias> not found" even though
 * the layer renders correctly, and declaring a browserStats output (as this script
 * does) makes no difference. This has been reported to CDSE. Until it is fixed, read
 * the mean off the map by tuning residualCenterK (see USER OPTIONS), or compute it with
 * the Statistical API. The browserStats output is kept so that the histogram starts
 * working here as soon as the fix lands.
 *
 * The spread of the residual is the other number worth having, and it needs the
 * Statistical API. It says how tightly the residual clusters around the offset - a
 * narrow spread means the field is mostly a constant bias plus stable spatial
 * structure, which is the assumption the sharpening rests on.
 *
 * Quick setup (Copernicus Browser -> "use additional datasets (advanced)"):
 *   - Primary dataset      Landsat 8-9 L1        -> alias  LANDSAT-OT-L1
 *   - Additional dataset 1 CLMS LST 3km hourly V3 -> alias  CLMS_CAL
 *   Set both timespans to the same calibration date, with CLMS_CAL on the hourly slot
 *   CONTAINING the Landsat overpass (~09:30 UTC over central Europe, so the 09:00
 *   slot). A time gap between the two corrupts the residual.
 *
 * Keep band, C, minValidC and maxValidC identical to script.js, or the offset measured
 * here is not the offset script.js needs.
 */

//// ---- USER OPTIONS ---------------------------------------------------------

// choose which Landsat TIRS band drives brightness temperature (B10 or B11)
var band = "B10";

// physical sanity range for the Landsat-derived LST (degC).
// Keep identical to script.js.
var minValidC = -20;
var maxValidC = 60;

// surface roughness term, added to the emissivity of mixed soil/vegetation pixels.
// 0.005 follows Sobrino et al. (2008) and Avdan & Jovanovska (2016).
// Keep identical to script.js.
var C = 0.005;

// mask pixels the Landsat Collection 2 QA band flags as cloud-affected. Leave this on
// when measuring the offset: cloud tops read cold and drag the mean residual down.
var applyCloudMask = true;

// also drop water pixels (NDVI < 0). Thermal inertia decouples water temperature from
// the land signal, so water residuals are large and noisy and a river or lake inside
// the area can shift the scene mean. Turn on if your area contains open water.
var excludeWater = false;

// color-map for the RGBA "default" output, in Kelvin degrees. This is a DIVERGING
// ramp spanning residualCenterK +/- residualSpreadK: blue = Landsat colder than CLMS,
// neutral = the residual equals residualCenterK, red = Landsat warmer.
//
// These two are how you READ THE OFFSET OFF THE MAP, by bisection - currently the only
// way to obtain it in Copernicus Browser (see header). Set residualCenterK to a guess:
// if the scene skews blue the true mean is lower, if it skews red it is higher. When
// the bulk of the land renders neutral, residualCenterK IS the mean residual, and
// biasOffset = -residualCenterK. Shrink residualSpreadK as you close in - at +/-2 K the
// mean is readable to a few tenths of a Kelvin degree.
//
// Leave the centre at 0 to see the raw sign pattern: whether Landsat reads warmer or
// colder than CLMS, and where.
var residualCenterK = 0;
var residualSpreadK = 15;

//// ---- LANDSAT LST CONSTANTS ------------------------------------------------

// NDVI thresholds for bare soil / full vegetation
var NDVIs = 0.2;
var NDVIv = 0.8;

// emissivity
var waterE = 0.991;
var soilE = 0.966;
var vegetationE = 0.973;

// central/mean wavelength in meters, B10 or B11
var bCent = band == "B10" ? 0.000010895 : 0.000012005;

// rho = h*c/sigma = PlanckC*velocityLight/BoltzmannC
var rho = 0.01438; // m K

//// ---- CLMS DECODE ----------------------------------------------------------

// raw CLMS LST -> Kelvin
var CLMS_FACTOR = 1 / 100;
var CLMS_OFFSET = 273.15;

//// ---- VISUALIZATION (diverging ramp, ColorBrewer RdBu reversed) -------------

// 11 stops spread evenly over the ramp. Blue = Landsat colder, the neutral middle stop
// = residualCenterK, red = Landsat warmer.
const RampColors = [
  [5, 48, 97],
  [33, 102, 172],
  [67, 147, 195],
  [146, 197, 222],
  [209, 229, 240],
  [247, 247, 247],
  [253, 219, 199],
  [244, 165, 130],
  [214, 96, 77],
  [178, 24, 43],
  [103, 0, 31],
];
var residualMinK = residualCenterK - residualSpreadK;
var residualMaxK = residualCenterK + residualSpreadK;
const ColorBar = RampColors.map(function (rgb, i) {
  var t = i / (RampColors.length - 1);
  return [residualMinK + t * (residualMaxK - residualMinK), rgb];
});
const visualizer = new ColorRampVisualizer(ColorBar);

//// ---- SETUP ----------------------------------------------------------------

function setup() {
  return {
    input: [
      {
        datasource: "LANDSAT-OT-L1",
        bands: ["B03", "B04", "B05", "B10", "B11", "BQA", "dataMask"],
      },
      { datasource: "CLMS_CAL", bands: ["LST", "dataMask"] },
    ],
    output: [
      { id: "default", bands: 4, sampleType: "UINT8" }, // RGBA color map
      { id: "index", bands: 1, sampleType: "FLOAT32" }, // residual (K)
      { id: "browserStats", bands: 1, sampleType: "FLOAT32" }, // residual, for stats
      { id: "dataMask", bands: 1 }, // validity mask
    ],
    mosaicking: "SIMPLE", // one clear scene per layer; access samples.X[0]
  };
}

//// ---- LANDSAT EMISSIVITY ---------------------------------------------------

function LSEcalc(NDVI, Pv) {
  var LSE;
  if (NDVI < 0) {
    // water
    LSE = waterE;
  } else if (NDVI < NDVIs) {
    // soil
    LSE = soilE;
  } else if (NDVI > NDVIv) {
    // vegetation
    LSE = vegetationE;
  } else {
    // mixtures of vegetation and soil
    LSE = vegetationE * Pv + soilE * (1 - Pv) + C;
  }
  return LSE;
}

// Landsat single-scene LST in KELVIN, or null if the pixel is invalid.
function landsatLST(s) {
  var Bi = band == "B10" ? s.B10 : s.B11;
  var B03i = s.B03;
  var B04i = s.B04;
  var B05i = s.B05;

  // raw-band guard: screens error scenes / fill (same as source script)
  if (!(Bi > 173 && Bi < 65000 && B03i > 0 && B04i > 0 && B05i > 0)) {
    return null;
  }

  // 1 Kelvin to C (brightness temperature)
  var btC = Bi - 273.15;
  // 2 NDVI
  var NDVI = (B05i - B04i) / (B05i + B04i);
  // 3 Pv - proportional vegetation
  var Pv = Math.pow((NDVI - NDVIs) / (NDVIv - NDVIs), 2);
  // 4 LSE - land surface emissivity
  var LSE = LSEcalc(NDVI, Pv);
  // 5 LST (degC): the single-channel land surface temperature calculation based on
  // (Artis & Carnahan, 1982), converting brightness temperature (btC) to Land
  // Surface Temperature (lstC) using the emissivity (LSE) and the central wavelength
  // (bCent) and rho (Planck's constant times the speed of light divided by
  // Boltzmann's constant).
  var lstC = btC / (1 + ((bCent * btC) / rho) * Math.log(LSE));

  // physical sanity range
  if (lstC < minValidC || lstC > maxValidC) {
    return null;
  }

  // return Kelvin
  return lstC + 273.15;
}

//// ---- EVALUATE PIXEL -------------------------------------------------------

function invalid() {
  return {
    default: [0, 0, 0, 0],
    index: [NaN],
    browserStats: [NaN],
    dataMask: [0],
  };
}

function evaluatePixel(samples) {
  // note: LANDSAT-OT-L1 has hyphens, so it must be accessed via bracket notation
  var landsat = samples["LANDSAT-OT-L1"][0];
  var clmsCal = samples.CLMS_CAL[0];

  if (landsat.dataMask === 0 || clmsCal.dataMask === 0) {
    return invalid();
  }

  if (applyCloudMask) {
    var qa = decodeL8C2Qa(landsat.BQA);
    if (qa.fill || qa.cloud || qa.dilatedCloud || qa.cirrus || qa.cloudShadow) {
      return invalid();
    }
  }

  if (excludeWater) {
    var ndvi = (landsat.B05 - landsat.B04) / (landsat.B05 + landsat.B04);
    if (ndvi < 0) {
      return invalid();
    }
  }

  // Landsat calibration LST (Kelvin), with raw-band + physical-range guards
  var lstLandsatCal = landsatLST(landsat);
  if (lstLandsatCal === null) {
    return invalid();
  }

  // CLMS LST (Kelvin)
  var lstClmsCal = clmsCal.LST * CLMS_FACTOR + CLMS_OFFSET;

  // the residual field. Its mean over a large clear-sky area, negated, is the
  // biasOffset to enter in script.js; its standard deviation says how well a single
  // scene-wide constant describes the Landsat-CLMS difference.
  var residual = lstLandsatCal - lstClmsCal;

  var rgb = visualizer.process(residual);
  return {
    default: rgb.concat(255),
    index: [residual],
    browserStats: [residual],
    dataMask: [1],
  };
}

General description

Copernicus Land Monitoring Service (CLMS) Land Surface Temperature (LST, 3 km hourly V3) collects data every hour and is published with 3-hour latency, but only at 3 km resolution. Landsat TIRS resolves LST at 100 m but has an 8–16 day revisit. Alone, neither of these can deliver daily high resolution land surface temperature for urban planning and monitoring of conditions at local scale.

This data fusion script sharpens the coarse CLMS LST using the fine spatial detail of a Landsat thermal image, following a simple additive-residual (statistical downscaling) approach.

First we calculate the difference between Landsat and CLMS temperature at a calibration date where we have both. This can be months or even years old, but should be from the same season and approximate temperature range as the target date. Then we use the difference to sharpen the CLMS temperature dataset for the time of interest. As a result, thermal conditions can be estimated at high resolution and daily frequency, suitable for complementing sensor-based heat monitoring and creating detailed maps to inform urban planning or disaster intervention.

Method

The underlying assumption is that the fine-scale spatial structure of urban LST anomalies (street canyons, parks, industrial zones) is approximately constant in time within a city. Given a cloud-free calibration date on which a Landsat overpass and a CLMS acquisition coincide, the script computes, per pixel and in Kelvin:

residual      = LST_landsat_cal − LST_clms_cal          (the fine-scale correction field)
LST_sharpened = LST_clms_target + residual + biasOffset

The residual captures the sub-3-km detail that CLMS cannot see. Adding it back to the CLMS LST for any target date produces a CLMS-scale temperature field that carries Landsat-scale spatial detail.

Inputs (three datasources)

Alias	Collection	Timespan	Bands
LANDSAT-OT-L1	Landsat 8-9 L1	calibration date (clear-sky overpass)	B03, B04, B05, B10, B11, BQA, dataMask
CLMS_CAL	CLMS LST 3km hourly V3	the hour matching the Landsat overpass	LST, dataMask
CLMS_TGT	CLMS LST 3km hourly V3	target date/hour to sharpen	LST, dataMask

The CLMS collection is deliberately added twice, with two different customized timespans — once for calibration (CLMS_CAL) and once for the target (CLMS_TGT). CLMS LST is a BYOC collection documented here: Collection ID 12225aec-26dd-4e2c-bbfd-994c253c1ba8, Location CDSE.

Landsat LST retrieval

The Landsat land surface temperature is derived using the mono-window / emissivity method from the Landsat-8 Land Surface Temperature Mapping script by Mohor Gartner: brightness temperature from TIRS band B10 (or B11), NDVI from B04/B05, a proportional-vegetation emissivity model, and the single-channel LST equation.

The method itself comes from Sobrino & Raissouni (2000) and Sobrino et al. (2008): pixels are classified by NDVI into water, bare soil, full vegetation and mixed, with the mixed class interpolated by proportional vegetation cover Pv (Carlson & Ripley, 1997) plus a surface roughness term C. The emissivity values themselves (water 0.991, soil 0.966, vegetation 0.973) originate from Wang et al. (2015), Table 8, computed from the ASTER Spectral Library (Baldridge et al., 2009). The NDVI thresholds are NDVIs = 0.2 (Sobrino & Raissouni, 2000) and NDVIv = 0.8 (Jiménez-Muñoz et al., 2009, who report 0.801 ± 0.012). The LST equation follows Artis & Carnahan (1982).

A note on the roughness constant C. These scripts use C = 0.005, the value given by Sobrino et al. (2008) and used by Avdan & Jovanovska (2016). The Landsat-8 LST Mapping script uses C = 0.009, a value that appears in Salih et al. (2018); we were unable to trace it to a primary measurement, so we kept the more widely cited figure. The practical difference is small — of order 0.2–0.3 K, and only on mixed soil/vegetation pixels — but it is a real difference, so C is exposed in the user options block of all three scripts. Set it to 0.009 if you need results directly comparable with the original Landsat-8 script.

CLMS LST decoding

Raw CLMS LST is decoded to Kelvin as LST_K = rawLST × 0.01 + 273.15, matching the CLMS LST 3km hourly V3 default script.

biasOffset — treating CLMS as truth

If you treat the coarse CLMS product as the stronger (but less detailed) truth, the scene-mean of the sharpened result should match CLMS rather than Landsat. A single evalscript cannot compute a per-CLMS-cell mean (evaluatePixel is strictly per-pixel, with no spatial aggregation), so the script exposes a single scene-wide constant biasOffset as the best scalar approximation.

There are two ways to obtain it, both external to script.js.

Measure the residual directly (recommended). Run bias_offset.js (see below) on your calibration date. It outputs LST_landsat_cal − LST_clms_cal per pixel, so biasOffset = −mean(residual). Because both terms come out of a single request, the same pixels are masked on both sides and the two retrievals are matched in time by construction. Open the Copernicus Browser statistics / histogram panel for the resulting layer over your area of interest. For a narrow, roughly symmetric distribution, visually estimate its centre as an approximation to the mean residual. Negate that estimate and enter it as biasOffset in script.js; see “Measuring biasOffset” below. The Browser does not report the mean or directly export histogram values; use the companion notebook or Statistical API when you need a calculated mean rather than a visual estimate.

Subtract two layer means. Run landsat_lst_reference.js and the CLMS_CAL layer over your area of interest, obtain each layer’s mean with the Statistical API or companion notebook, and set biasOffset = mean(CLMS_cal) − mean(Landsat_cal). Both are single-dataset layers. The weakness is that the two layers do not mask the same pixels — Landsat drops cloud-affected and out-of-range pixels that CLMS keeps — so you are differencing means taken over different pixel sets.

Leave biasOffset at 0 to keep pure additive Landsat detail.

Expect this offset to be sizeable (<10 K). Landsat L1 band B10 is top-of-atmosphere brightness temperature with no atmospheric correction, so it reads systematically colder than the atmospherically corrected CLMS LST — over Budapest on the example calibration date the difference is around 2.2 K. What matters for the method is not that the bias is small but that it is stable, which is exactly what a single scalar can absorb.

changeScale and deltaTMean — damping the CLMS change (optional)

The additive formula passes the full 3 km change that CLMS sees between the calibration and the target date on to every 100 m pixel. That change is partly real and partly retrieval noise and resampling. In a leave-one-date-out test over Budapest (28 Landsat dates, companion notebook Section 12.4), the sharpened map came closer to the target-date Landsat when only part of the change was used. The median pattern error was 1.78 K with the additive model and 1.53 K with the change scaled by the calibration-date slope, against 2.16 K for plain CLMS. This is one city and a small box, and the reference is Landsat itself, so the additive model stays the default. The scaling is offered as an option:

LST_sharpened = LST_landsat_cal + s × (LST_clms_target − LST_clms_cal) + (1 − s) × deltaTMean + biasOffset

• changeScale is s. At 1 (the default) the formula is exactly the additive one above and deltaTMean is ignored. At 0 the Landsat calibration pattern is used unchanged and only the scene-mean change is added.
• A sensible s is the slope of a linear regression of Landsat LST on CLMS LST over your area: how many Kelvin degrees Landsat changes per Kelvin degree of CLMS, spatially, within one date. Take the median over several calibration dates, because the slope of a single date is noisy (over Budapest the example calibration date gave 0.89, against a median of 0.54 over 28 dates). The companion notebook computes it (Section 12.5). Over Budapest a sweep of fixed values gave nearly the same error for any s between 0.25 and 0.5, rising steadily towards s = 1. In other cities the median per-date slope ranged from about 0.25 to 0.8, so measure it rather than borrow it. A fixed s for a city worked as well as a per-date one.
• deltaTMean is mean(CLMS_TGT) − mean(CLMS_CAL) over your area. Without this term an s below 1 would also damp the scene mean, and the sharpened map would drift away from the CLMS level (by 5 K for s = 0.5 and a 10 K warmer target). Obtain both means with the Statistical API or companion notebook. The Browser histograms can provide approximate centres only if both distributions are narrow and roughly symmetric. The difference is the same in °C and in Kelvin.

Validity and masking

A pixel is masked (transparent, NaN index) when any of the following holds: the Landsat pixel fails the raw-band guard (B10 outside 173–65000, or B03/B04/B05 ≤ 0); the Landsat-derived LST falls outside a configurable physical range (default −20 °C … +60 °C, catching undetected cloud / emissivity blow-ups); or any Landsat or CLMS input has dataMask === 0.

With applyCloudMask = true (the default), Landsat Collection 2 QA flags decoded from BQA with decodeL8C2Qa also mask fill, cloud, dilated cloud, cirrus and cloud shadow. This is the same cloud mask used by both companion scripts and the notebook. Masking happens before the residual is calculated, so calibration-date clouds cannot propagate into target-date maps as false cool areas. Those pixels remain no-data even if the target day is clear; use a clearer calibration scene to recover coverage. Keep applyCloudMask identical across all three scripts when comparing their results.

Outputs

Three named outputs: default (RGBA color map on the CLMS magma-style ramp), index (FLOAT32 sharpened LST in Kelvin), and dataMask.

The color ramp range is exposed as paletteMinK / paletteMaxK. The defaults, 240–330 K (−33 … +57 °C), reproduce the CLMS LST 3km hourly V3 palette exactly, which is what makes the sharpened layer directly comparable with the native CLMS layer in the Browser. But that range is built to cover the whole globe across a whole year, so a single city on a single morning occupies only a narrow slice of it and most of the ramp goes unused. Narrowing the range — 295–325 K (+22 … +52 °C) suits a European summer morning — spreads the same ten colors over the temperatures actually present and makes intra-urban contrast much easier to read. Both values are in Kelvin (°C = K − 273.15). As with biasOffset, use the same values in both scripts when comparing layers; different ranges make two layers look different for purely cosmetic reasons.

There is deliberately no browserStats output. That output exists to feed the Copernicus Browser time series tool, which works by re-running the request over a sequence of time intervals — but this script pins CLMS_CAL to a fixed calibration date, so a swept time range has no coherent meaning and the tool returns an error.

The Copernicus Browser statistics / histogram panel supports data fusion layers. Use it to inspect the sharpened LST from script.js or the calibration residual from bias_offset.js over your area of interest. A histogram describes the spatial distribution for the configured dates; it does not require sweeping the calibration timespan as a time series would.

landsat_lst_reference.js does declare browserStats: it is a single-dataset script with no pinned timespans, so both the time series tool and the statistics panel behave as expected.

The three scripts

File	Datasets	Purpose
script.js	Landsat + CLMS ×2 (data fusion)	the sharpening itself
landsat_lst_reference.js	Landsat 8-9 L1 only	native Landsat LST, to check the sharpening against
bias_offset.js	Landsat + CLMS ×1 (data fusion)	the Landsat − CLMS residual, to measure biasOffset

landsat_lst_reference.js is deliberately a copy of the Landsat half of script.js — same band, same emissivity model, same constants, same guards, same color ramp. Because the retrieval is identical, any difference between the two layers is attributable to the sharpening rather than to a difference in method. It applies the same default-on Collection 2 QA cloud mask as script.js. It additionally declares a browserStats output, so the Browser time series tool works on it (see “Outputs” above).

bias_offset.js outputs the residual field itself — LST_landsat_cal − LST_clms_cal in Kelvin degrees — rather than a sharpened temperature. It uses the same Landsat retrieval and the same calibration-date pairing as script.js, so the value you measure with it is the value script.js needs. It applies the same Collection 2 QA cloud mask: cloud tops read cold and would bias the mean in the same direction as the atmospheric effect being measured, making the two indistinguishable. It can additionally drop water pixels (excludeWater), whose thermal inertia decouples them from the land signal so that a river or lake shifts the scene mean. It renders on a diverging blue–neutral–red ramp rather than the CLMS magma ramp, because the residual is a signed quantity whose whole interpretive content is the sign, and a sequential ramp gives the reader no perceptual anchor at zero.

How to use in Copernicus Browser

Sharpening (script.js)

The CLMS LST product is not in the standard Copernicus Browser dataset list — it is added as a BYOC collection, so you will need its Collection ID:

Collection ID	12225aec-26dd-4e2c-bbfd-994c253c1ba8
Type	BYOC
Location	CDSE

(from the CDSE collection documentation, which also documents the bands and the LST × 0.01 + 273.15 decoding)

1. Find a clear, cloud-free Landsat 8-9 L1 image over your target area, preferably in the same season as your target date. Use the Search panel of the Browser to find the data product and extract the exact sensing time. Note this - you will be using it later. Click visualize.

2. On the left panel, click Custom → Custom script, and tick “use additional datasets (advanced)”. This makes Landsat the primary dataset.

3. Primary dataset — Landsat 8-9 L1. This is the calibration scene that supplies the spatial detail.
   • Datasource alias: LANDSAT-OT-L1 (the default alias Copernicus Browser offers)
   • Customize timespan: the calibration date — a clear-sky Landsat overpass
4. Additional dataset 1 — CLMS LST 3km hourly V3 (BYOC, Collection ID above, Location CDSE). This is the calibration CLMS field, subtracted from Landsat to form the residual.
   • Datasource alias: CLMS_CAL
   • Customize timespan: the same day as the Landsat overpass, and the hourly slot containing it — e.g. 09:00–09:59 UTC for a 09:32 overpass. This must match; see “Matching the hour”.
5. Additional dataset 2 — CLMS LST 3km hourly V3 (same Collection ID, added a second time). This is the field you actually want to sharpen.
   • Datasource alias: CLMS_TGT
   • Customize timespan: the target date/hour, normally the same hour of day as the calibration

6. Paste script.js into the script editor and run.

7. Optionally adjust the user options at the top of the script — paletteMinK / paletteMaxK to suit the temperature range of your scene, biasOffset once you have measured it, and changeScale / deltaTMean if you want to damp the CLMS change (see above).

The aliases are case-sensitive and must match the setup() datasources exactly, or the request will not bind. Note that the CLMS collection is deliberately added twice, with two different timespans — this is what lets a single request see both the calibration and the target field.

Self-consistency check: set CLMS_TGT to the same date as CLMS_CAL. The sharpened output must then equal the pure Landsat LST (the residual and CLMS_cal cancel), which verifies the math.

Matching the hour

Two separate timing choices matter here, and they are not equally strict.

CLMS_CAL must match the Landsat overpass hour — same day, same hourly slot. This one is not optional. The residual is the difference between two supposedly simultaneous observations of the same surface; any time gap injects real thermal change into it. Check the acquisition time of your own scene rather than assuming — the Browser shows it in the scene metadata.

CLMS_TGT should normally be the same hour of day, but this is a soft constraint. The reason is that the residual field is a morning pattern: it encodes how each surface deviates from its 3 km neighbourhood at mid-morning, under that day’s sun angle. Applying it to an afternoon or night-time CLMS field assumes intra-urban thermal contrast keeps the same spatial shape through the diurnal cycle, which is a stronger assumption than the temporal stability the method already rests on. Surfaces with different thermal inertia — water, asphalt, tree canopy — diverge over the day, and after sunset the pattern can partially invert as dense built-up areas release stored heat while parks and water cool faster. Keeping both CLMS layers at the overpass hour avoids the question entirely.

That said, sharpening a different hour is the more useful application — an afternoon heat peak is what a planner actually cares about, and CLMS is hourly precisely so you can look at it. So treat the same-hour rule as the default, and a different hour as a deliberate experiment worth validating before trusting: run the sharpening for the target hour and compare against whatever independent thermal observation you can get for that hour. How far the assumption stretches across the diurnal cycle is an open question this script is well suited to test, not a settled result.

Measuring biasOffset (bias_offset.js)

Set up exactly as for the sharpening, but with only one additional dataset:

• Primary dataset: Landsat 8-9 L1 → alias LANDSAT-OT-L1
• Additional dataset 1: CLMS LST 3km hourly V3 (Collection ID above) → alias CLMS_CAL

Put both timespans on the calibration date, with CLMS_CAL on the hourly slot containing the Landsat overpass. There is no CLMS_TGT: this script measures the calibration pair, it does not sharpen anything. Keep band, C, minValidC, maxValidC and applyCloudMask identical to script.js — the offset is only valid for the retrieval that produced it. If your area contains open water, set excludeWater = true before estimating the offset; this estimates a land-only offset rather than one for all valid pixels.

The layer renders on a diverging ramp: blue where Landsat reads colder than CLMS, neutral where the residual equals residualCenterK, red where Landsat reads warmer. With the default residualCenterK = 0 a typical scene comes out uniformly blue — that is the top-of-atmosphere cold bias, and it is the thing you are about to quantify.

1. Paste bias_offset.js into the custom script editor and run the calibration pair.
2. Select the area of interest you intend to sharpen and open the statistics / histogram panel for the residual layer. Use the same area when comparing results or measuring another date.
3. Inspect the full histogram, including its tails. If it has one narrow, roughly symmetric cluster, visually estimate the residual value at its centre from the horizontal axis. This approximates the mean; the tallest bin is the mode, which is near the mean only for a roughly symmetric, single-peaked distribution. Do not use the midpoint of the color ramp as the estimate.
4. Set biasOffset to the negative of that estimated centre in script.js. For example, a histogram narrowly clustered around −2.2 K suggests biasOffset = 2.2. The values are temperature differences in Kelvin, numerically identical to differences in degrees Celsius; do not subtract 273.15. Re-run the sharpening with that value.

This is an approximate Browser-only calibration: the histogram does not report a mean or directly export bin values. Do not infer precision finer than the histogram bins support. Narrow residual distributions are expected when a scene-wide bias dominates, but check rather than assume this. A broad, skewed or multi-peaked distribution makes a visual mean estimate unreliable and suggests spatial differences that a constant cannot remove. In that case, or when numerical accuracy matters, use the companion notebook or Statistical API to calculate biasOffset = −mean(residual) over the same area and with the same masking settings.

Masked pixels are excluded from the statistics. residualCenterK and residualSpreadK only control the display colors; they do not change the residual values or the estimated offset.

Validating against native Landsat (landsat_lst_reference.js)

The sharpening is only useful if it reproduces real thermal structure on a date it was not calibrated on. To check that yourself:

1. Pick a second Landsat overpass, different from the calibration date, that is also clear-sky.
2. Run script.js with CLMS_TGT set to the hour of that second overpass, keeping CLMS_CAL and LANDSAT-OT-L1 on the original calibration date. This is a genuine prediction — no information from the second overpass enters the calculation.
3. In a second layer, run landsat_lst_reference.js on that second overpass date. This is the truth to compare against.
4. Compare the two visually and use the Browser statistics / histogram panel to compare their distributions over the same area. Use the Statistical API or companion notebook for calculated means. Statistics are available for both the native Landsat reference and the sharpened data fusion layer; different valid-pixel masks can still affect the comparison.

Keep biasOffset consistent between the two scripts. With biasOffset = 0 in both, the layers share the Landsat top-of-atmosphere scale and the difference you see is purely spatial sharpening — this is the right setting for judging whether the fine structure is reproduced. With the same non-zero offset in both, they share the CLMS absolute scale, which is the right setting for judging absolute temperatures. Mixing the two produces a scene-wide difference that says nothing about the method. The same applies to paletteMinK / paletteMaxK and to C: any comparison between the two layers is only meaningful if all of these match.

Expect the sharpened layer to look smoother than native Landsat even when the method works well: the residual field is frozen on the calibration date, so thermal structure that genuinely changed between the two dates cannot be reproduced. That limitation is the assumption being tested.

Description of representative images

Sharpened land surface temperature over Budapest, Hungary. The calibration date is a clear-sky Landsat overpass on 13 June 2025 (09:32:47 UTC, matched to the 09:00 CLMS slot); the target CLMS acquisition is 10 August 2026, 09:00 UTC — the same hour of day, during an ongoing heatwave in Budapest, and processed the same day thanks to the ~3 hour CLMS latency. The output shows CLMS-scale temperatures carrying Landsat-scale spatial detail — streets, parks, the Danube, and built-up districts become distinguishable at a resolution the native 3 km CLMS product cannot provide.

Sharpened CLMS LST over Budapest

The Danube stands out as a cool ribbon, the Buda hills on the west bank read cooler than the dense Pest districts on the east, and individual parks, rail yards and industrial roofs are resolved. None of this structure exists in the native 3 km CLMS field, where the whole city spans only a handful of pixels.

Native CLMS LST (left) and sharpened CLMS LST (right), both over a VHR true color mosaic

Native CLMS LST (3 km)	Landsat-sharpened CLMS LST

Displayed semi-transparently over a 2024 VHR true color mosaic, the thermal pattern can be read against the urban fabric that produces it: warm anomalies line up with large impervious surfaces and industrial areas, cool ones with the river, the wooded hills, and the larger parks.

The left-hand panel is the input the sharpening starts from — same view, same hour, same color ramp. At 3 km the entire metropolitan area is covered by a handful of pixels: the Danube is invisible, the Buda hills and the Pest districts fall inside the same cells, and the only discernible signal is a broad west-to-east gradient with hard block boundaries. This is the resolution gap the method closes — and it is also why the sharpened product remains useful despite its assumptions, since the alternative is not a better thermal map but essentially no intra-urban information at all.

The two panels also sit at visibly different points on the color ramp, because these examples use biasOffset = 0: the sharpened field inherits the Landsat top-of-atmosphere temperature scale and so reads cooler than native CLMS. Setting biasOffset as described above brings them onto a common absolute scale.

Native Landsat LST (left) versus Landsat-sharpened CLMS LST (right)

This is the validation test described above, with both layers on the same color ramp. The faint vertical line left of the Danube is the split between the two layers, not a feature of either — and how hard it is to find is the point of the figure. Thermal features run across the boundary without a visible step: the Danube, the cool wooded blocks in the west, the warm industrial strips in the east. The sharpened half is marginally smoother, as expected, since its fine structure comes from the frozen calibration-date residual and cannot reproduce detail that changed between the two dates.

Example location and dates (Copernicus Browser):

• Location: Budapest, Hungary — zoom=13, lat=47.49523, lng=19.05913
• LANDSAT-OT-L1 timespan: 13 June 2025 (overpass at 09:32:47 UTC)
• CLMS_CAL timespan: 13 June 2025, 09:00–09:59 UTC — the hourly slot containing the overpass
• CLMS_TGT timespan: 10 August 2026, 09:00–09:59 UTC

Caveats

• Structures influencing microclimate are relatively constant in cities but can be less constant in more natural landscapes as vegetation condition or soil water content changes. Use with extra caution outside cities, keep the calibration scene as close as possible in time to the target scene.
• Clouds can affect all three datasets. Data gaps appear where either CLMS input is missing or, with applyCloudMask = true (the default), the Landsat calibration pixel is flagged as cloud-affected. Clouds missed by the QA mask can still introduce erroneous values into the calibration residual and the sharpened maps.
• The residual also absorbs the systematic bias between the two different LST retrieval algorithms (Landsat vs CLMS), not purely spatial structure.
• The additive-in-Kelvin model treats anomaly magnitude as background-independent. This means that the differences we introduce from the residuals do not scale with the actual temperature, they stay constant. It is a good idea to select a calibration scene from the same season and from the same generic conditions if possible.
• The size of the temperature change between the calibration and the target scene limits the method. In the Budapest validation (companion notebook, Section 12.2) the sharpened map was clearly better than plain CLMS when the two scene means differed by less than about 10 K, still better at 10–20 K, and worse than plain CLMS above about 20 K, e.g. a winter calibration scene applied to a summer target. Above 20 K the fine-scale pattern itself changes with the season, and no choice of biasOffset or changeScale fixes that. The time between the dates mattered much less than the temperature difference. These thresholds come from one city and should be taken as a guide.
• All three scripts apply the Landsat Collection 2 QA cloud mask by default, but it is not infallible; thin cloud or cloud edges may escape detection. Prefer a clear-sky Landsat calibration scene and inspect it visually. Flagged pixels remain no-data on every target date; the script does not fill these gaps even when the target day is clear.
• Sharpening cannot recover thermal structure that changed after the calibration date — new construction, a harvested field, a recently irrigated park. Re-calibrate periodically.
• The residual is a snapshot of one hour of the diurnal cycle. Calibration requires CLMS_CAL at the Landsat overpass hour; applying the result to a very different hour of day is an untested extrapolation (see “Matching the hour” above).

Companion notebook and citation

The method is documented and evaluated in a Jupyter notebook in the CDSE notebook samples: clms_lst_landsat_sharpening.ipynb. It reproduces the sharpening for Budapest with the Sentinel Hub Process API, then tests the method: how stable the residual pattern is over time, whether the Landsat − CLMS bias is the same in ten cities across climate zones, whether the additive model holds, and how large the error of a sharpened map is in a leave-one-date-out validation against Landsat. It runs on the CDSE JupyterLab.

If you use this method, please cite the notebook:

Zlinszky, A. (2026). Landsat-sharpened CLMS Land Surface Temperature: method, validation and error budget (Version 1.0) [Jupyter notebook]. Zenodo. https://doi.org/10.5281/zenodo.22944002

Author of the script

• Andras Zlinszky

References

Scripts and data:

• Mohor Gartner, 2019, Land Surface Temperature (LST) Mapping Script
• Copernicus Land Monitoring Service, Land Surface Temperature 3km hourly V3
• CDSE, LST global 3km hourly V3 collection documentation — BYOC Collection ID, bands and decoding
• Sentinel Hub, About Landsat 8-9 OLI-TIRS Collection 2 Level 1 Data
• Sentinel Hub, Data Fusion

Landsat LST retrieval:

• Artis, D. A., & Carnahan, W. H. (1982). Survey of emissivity variability in thermography of urban areas. Remote Sensing of Environment, 12(4), 313–329. doi:10.1016/0034-4257(82)90043-8
• Avdan, U., & Jovanovska, G. (2016). Algorithm for automated mapping of land surface temperature using LANDSAT 8 satellite data. Journal of Sensors, 2016, 1480307. doi:10.1155/2016/1480307
• Baldridge, A. M., Hook, S. J., Grove, C. I., & Rivera, G. (2009). The ASTER spectral library version 2.0. Remote Sensing of Environment, 113(4), 711–715. doi:10.1016/j.rse.2008.11.007
• Carlson, T. N., & Ripley, D. A. (1997). On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sensing of Environment, 62(3), 241–252. doi:10.1016/S0034-4257(97)00104-1
• Jiménez-Muñoz, J. C., Sobrino, J. A., Plaza, A., Guanter, L., Moreno, J., & Martínez, P. (2009). Comparison between fractional vegetation cover retrievals from vegetation indices and spectral mixture analysis. Sensors, 9(2), 768–793. doi:10.3390/s90200768
• Salih, M. M., Jasim, O. Z., Hassoon, K. I., & Abdalkadhum, A. J. (2018). Land surface temperature retrieval from LANDSAT-8 thermal infrared sensor data. International Journal of Engineering & Technology, 7(4.20), 608–612.
• Sobrino, J. A., & Raissouni, N. (2000). Toward remote sensing methods for land cover dynamic monitoring: application to Morocco. International Journal of Remote Sensing, 21(2), 353–366. doi:10.1080/014311600210876
• Sobrino, J. A., Jiménez-Muñoz, J. C., Sòria, G., Romaguera, M., Guanter, L., Moreno, J., Plaza, A., & Martínez, P. (2008). Land surface emissivity retrieval from different VNIR and TIR sensors. IEEE Transactions on Geoscience and Remote Sensing, 46(2), 316–327. doi:10.1109/TGRS.2007.904834
• Wang, F., Qin, Z., Song, C., Tu, L., Karnieli, A., & Zhao, S. (2015). An improved mono-window algorithm for land surface temperature retrieval from Landsat 8 thermal infrared sensor data. Remote Sensing, 7(4), 4268–4289. doi:10.3390/rs70404268