GEE：基于Landsat8计算陆地表面温度（Land Surface Temperature，LST）

作者：CSDN @ _养乐多_

本文将介绍在Google Earth Engine（GEE）平台上使用 Landsat 8 卫星影像数据计算陆地表面温度（Land Surface Temperature，LST）的代码。

结果如下图所示，

---

---

一、参考内容

单窗算法：Google Earth Engine Open-Source Code for Land Surface Temperature Estimation from the Landsat Series

博客：https://medium.com/@ridhomuh002/analyzing-land-surface-temperature-lst-with-landsat-8-data-in-google-earth-engine-f4dd7ca28e70

主要公式包括：

植被覆盖率（FV）计算公式：$$((NDVI-NDV{I}_{m}in)/(NDV{I}_{m}ax-NDV{I}_{m}in){)}^2$$
辐射率（EM）计算公式：$$0.004\ast FV+0.986$$
陆地表面温度（LST）计算公式：$$(TB/(1+(\lambda \ast (TB/1.438))\ast ln(em)))-273.15$$
其中，TB 代表热红外波段的亮度温度，em 代表辐射率，λ 代表波长（对应 Landsat 8 中的波段 10）。

二、示例代码链接

https://code.earthengine.google.com/9ce60d3992ea71ddd75c57e7de3acfa8?noload=true

三、完整示例代码

//#####################################################################################
//#                                                                                 #\\
//#                    Land Surface Temperature(LST) MAPPING                        #\\
//#                                                                                 #\\
//#####################################################################################

// date: 2024-03-11
// author: YangLeDuo

// 自定义研究区矢量边界的显示样式
var roi = table
var styling={color:'green',fillColor:'00000000'}
Map.addLayer(roi.style(styling),{},"roi");
Map.centerObject(roi, 8);

// 设置单波段配色方案         
var Single_bandVisParam = {  
  min: 8,  
  max: 54,  
  palette: ['C2523C','F2B60E','77ED00','1BAA7D','0B2C7A',]
};

//利用QA波段生成云掩模
function maskclouds(image) {                
         var cloudShadowBitMask = (1 << 3);
         var cloudsBitMask = (1 << 5);
         var qa = image.select('QA_PIXEL');
         var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

function L8applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}
//#############################################################################################
//定义数据集
var L8_image = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                    .filterDate('2020-01-01', '2020-12-31')
                    .filter(ee.Filter.calendarRange(1,12,"month")) //选择月份
                    .filterBounds(roi)
                    .filter(ee.Filter.lt('CLOUD_COVER', 20))
                    .map(maskclouds).map(L8applyScaleFactors)
                    .median()
                    .clip(roi);
print('L8_dataset',L8_image)

//NDVI
var L8_NDVI = L8_image.expression('(NIR-Red)/(NIR+Red)', {
      'NIR': L8_image.select('SR_B5'),
      'Red': L8_image.select('SR_B4'),
    })
    
// 计算AOI内的最小NDVI值
var ndviMin = ee.Number(L8_NDVI.reduceRegion({
  reducer   : ee.Reducer.min(),
  geometry  : roi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

// 计算AOI内的最大NDVI值
var ndviMax = ee.Number(L8_NDVI.reduceRegion({
  reducer   : ee.Reducer.max(),
  geometry  : roi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

// 植被覆盖率（FV）计算
// 公式: ((NDVI - NDVI_min) / (NDVI_max - NDVI_min))^2

// 使用指定范围内的NDVI值计算植被覆盖率（FV）。
// NDVI_min表示最小NDVI值，NDVI_max表示最大NDVI值
var fv = ((L8_NDVI.subtract(ndviMin)).divide(ndviMax.subtract(ndviMin)))
          .pow(ee.Number(2))
          .rename('FV');
        
  
// 辐射率计算
// 公式: 0.004 * FV + 0.986

// 使用植被覆盖率（FV）计算陆地表面辐射率（EM）。
// 0.004系数表示由于植被引起的辐射率变化，
// 0.986表示其他表面的基本辐射率。

var em = fv.multiply(ee.Number(0.004)).add(ee.Number(0.986)).rename('EM');

// 选择热带波段（波段10）并重命名
var thermal = L8_image.select('ST_B10').rename('thermal');

// 现在，让我们计算陆地表面温度（LST）

// 公式: (TB / (1 + (λ * (TB / 1.438)) * ln(em))) - 273.15
var L8_LST = thermal.expression(
  '(TB / (1 + (0.00115 * (TB / 1.438)) * log(em))) - 273.15', {
    'TB': thermal.select('thermal'), // 选择热带波段（TB）
    'em': em // 分配辐射率（em）
  }).rename('LST');
  
var lstMin = ee.Number(L8_LST.reduceRegion({
  reducer   : ee.Reducer.min(),
  geometry  : roi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

// 计算AOI内的最大NDVI值
var lstMax = ee.Number(L8_LST.reduceRegion({
  reducer   : ee.Reducer.max(),
  geometry  : roi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

print("最小陆地表面温度:", lstMin);
print("最大陆地表面温度:", lstMax);

Map.addLayer(L8_LST, Single_bandVisParam, 'LST');