GEE实践应用|热岛效应（一）地表温度计算

本文主要是介绍GEE实践应用|热岛效应（一）地表温度计算，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1.学习目标

2.理论介绍

3.从MODIS获得地表温度

4.从Landsat卫星获得地表温度

1.学习目标

①了解如何使用GEE计算地表温度

2.理论介绍

城市化涉及用建筑物、道路和停车场等建筑结构取代自然景观。这种土地覆盖的改变也改变了土地表面的特性。这些变化的范围从表面反射和吸收的辐射量到热量如何从表面消散（例如，为城市发展而去除植被会减少蒸发冷却）。这些地表特性的变化可以改变当地的天气和气候。由于城市化而引起的当地气候变化被研究最多的是城市热岛（UHI）效应。城市热岛效应是指一个城市比其周围环境或同等的未城市化地表温度更高的现象，我们对城市热岛效应的了解已有近 200 年的历史。

传统上，城市热岛指数被定义为气象站测量的城市与城外一些农村参考点之间的气温差异。这种方法的一个问题是，城市的不同部分可能有不同的气温，因此很难捕获整个城市的城市热岛。使用热带中的卫星观测使我们能够获得另一种温度测量方法：the radiometric skin temperature，通常称为陆地温度表面温度（LST）。我们可以使用LST来计算地表UHI（SUHI）强度，包括它在城市内如何在像素尺度上变化。这里需要强调的是，卫星观测到的UHI值和使用空气计算出的UHI值温度测量结果可能会有很大差异。

3.从MODIS获得地表温度

地表温度可以从 MODIS Terra 和 Aqua 卫星产品等衍生产品中提取，也可以直接根据热波段中的卫星测量结果进行估算。我们将使用中国重庆市作为感兴趣区域来探索这两种选择（图 1）。我们将从 MODIS LST开始。

我们首先加载重庆市的边界矢量数据，使用的代码如下，加载结果如图1所示：


// Load feature collection of chongqing from user assets.
var region = ee.FeatureCollection(table);
// Get dissolved feature collection using an error margin of 50 meters.
var region = region.union(50);
// Set map center and zoom level (Zoom level varies from 1to 20).
//Map.setCenter(region, 12);
// Add layer to map.
Map.addLayer(region, {}, 'chongqing boundary');

图1|重庆的边界

接下来，我们加载 MODIS MYD11A2 版本 6 产品，该产品提供来自 Aqua 卫星的八天 LST 合成数据。这对应于大约下午 1:30 的穿越赤道时间，白天和夜间凌晨1:30。相比之下，Terra 平台上的 MODIS 传感器（MOD11A2 版本 6）的通过时间约为当地时间上午10:30 和晚上10:30。


// Load MODIS image collection from the Earth Engine datacatalog.
var modisLst = ee.ImageCollection('MODIS/006/MYD11A2');
// Select the band of interest (in this case: Daytime LST).
var landSurfTemperature = modisLst.select('LST_Day_1km');

我们仅想关注夏季 SUHI，因此我们将使用一年中的日期过滤器组合每年6月1日（第152天）到8月31 日（第243天）的图像来创建LST的五年夏季合成：


// Create a summer filter.
var sumFilter = ee.Filter.dayOfYear(152, 243);
// Filter the date range of interest using a date filter.
var lstDateInt = landSurfTemperature
.filterDate('2014-01-01', '2019-01-01').filter(sumFilter);
// Take pixel-wise mean of all the images in thecollection.
var lstMean = lstDateInt.mean();

现在，我们将此图像转换为以摄氏度为单位的 LST，并遮盖所有水像素（水的高比热容会影响 LST，我们关注于陆地像素）。对于水掩模，我们使用全球地表水数据集，为了转换像素值，我们使用数据提供者提供的波段的缩放因子，然后减去 273.15 以从开尔文转换为摄氏度，比例因子可以在Earth Engine 数据摘要页面中找到（图2）。

图2 数据总结中的比例因子

最后，我们使用城市边界裁剪图像并将图层添加到地图中，如图3所示。

// Multiply each pixel by scaling factor to get the LST values.
var lstFinal = lstMean.multiply(0.02);
// Generate a water mask.
var water = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select('occurrence');
var notWater = water.mask().not();
// Clip data to region of interest, convert to degree Celsius, and mask water pixels.
var lstNewHaven = lstFinal.clip(region).subtract(273.15)
.updateMask(notWater);
// Add layer to map.
Map.addLayer(lstNewHaven, {
palette: ['blue', 'white', 'red'],
min: 25,
max: 38
},
'LST_MODIS');

图3|中国重庆白天 MODIS Aqua LST 的五年夏季合成图。红色像素显示较高的 LST 值，蓝色像素显示较低的值

4.从Landsat卫星获得地表温度

使用 MODIS LST 相对简单，因为数据已经由 NASA 团队处理。我们还可以从Landsat获取LST，它的原始分辨率（根据卫星的不同，在60m到~120 m之间）比1kmMODIS像素要精细得多。然而，我们需要根据热波段中的测量结果自行推导出LST，这通常还涉及对表面发射率的一些估算。

物体的表面发射率 (ε) 是指与相同温度下的黑体相比，物体发射热辐射的效率，范围可以从 0（对于完美反射体）到1（对于完美吸收体和发射体）。由于卫星捕获的热辐射是 LST 和 ε 的函数，因此您需要准确地规定或估计ε以获得正确的 LST。让我们考虑一种使用 Landsat8数据的简单方法。

我们将首先加载陆地卫星数据、云筛选，然后过滤到感兴趣的时间和区域。继续在同一脚本中添加以下代码：

// Function to filter out cloudy pixels.
function cloudMask(cloudyScene) {
// Add a cloud score band to the image.
var scored =ee.Algorithms.Landsat.simpleCloudScore(cloudyScene);
// Create an image mask from the cloud score band and specify threshold.
var mask = scored.select(['cloud']).lte(10);
// Apply the mask to the original image and return the masked image.
return cloudyScene.updateMask(mask);
}
// Load the collection, apply coud mask, and filter to date and region of interest.
var col = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')
.filterBounds(region)
.filterDate('2014-01-01', '2019-01-01')
.filter(sumFilter)
.map(cloudMask);
print('Landsat collection', col);

创建中值合成作为进一步减少云影响的简单方法后，我们遮盖水像素并选择亮度温度波段。

// Generate median composite.
var image = col.median();
// Select thermal band 10 (with brightness temperature).
var thermal = image.select('B10')
.clip(region)
.updateMask(notWater);
Map.addLayer(thermal, {
min: 295,
max: 310,
palette: ['blue', 'white', 'red']
},
'Landsat_BT');

亮度温度（图 4）是相当于从大气层顶部逸出的红外辐射的温度，假设地球是一个黑体。它与LST不同，LST 需要考虑大气吸收和再发射以及地表的发射率。导出像素级发射率的一种方法是作为像素植被比例的函数。为此，我们首先根据 Landsat 表面反射率数据计算归一化植被指数 (NDVI)（见图 5）。


// Calculate Normalized Difference Vegetation Index (NDVI)
// from Landsat surface reflectance.
var ndvi = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(region)
.filterDate('2014-01-01', '2019-01-01')
.filter(sumFilter)
.median()
.normalizedDifference(['SR_B5',
'SR_B4']).rename('NDVI')
.clip(region)
.updateMask(notWater);
Map.addLayer(ndvi, {
min: 0,
max: 1,
palette: ['blue', 'white', 'green']
},
'ndvi');

图4|重庆市五年夏季陆地卫星亮温合成中值，红色像素显示较高的值，蓝色像素显示较低的值

图5|重庆市上空5年夏季Landsat NDVI合成中值,白色像素显示较高的 NDVI 值，蓝色像素显示较低的值

为了将每个像素的 NDVI 映射到植被像素的实际部分（植被覆盖率），我们接下来使用基于每个像素的 NDVI 值范围的关系。

// Find the minimum and maximum of NDVI. Combine the reducers
// for efficiency (single pass over the data).
var minMax = ndvi.reduceRegion({
reducer: ee.Reducer.min().combine({
reducer2: ee.Reducer.max(),
sharedInputs: true
}),
geometry: region,
scale: 30,
maxPixels: 1e9
});
print('minMax', minMax);
var min = ee.Number(minMax.get('NDVI_min'));
var max = ee.Number(minMax.get('NDVI_max'));// Calculate fractional vegetation.
var fv = ndvi.subtract(min).divide(max.subtract(min)).rename('FV');
Map.addLayer(fv, {
min: 0,
max: 1,
palette: ['blue', 'white', 'green']
}, 'fv');

现在，我们使用基于植被覆盖率的发射率经验模型。

// Emissivity calculations.
var a = ee.Number(0.004);
var b = ee.Number(0.986);
var em =
fv.multiply(a).add(b).rename('EMM').updateMask(notWater);
Map.addLayer(em, {
min: 0.98,
max: 0.99,
palette: ['blue', 'white', 'green']
},
'EMM');

如图 6 所示，与植被上方相比，建筑结构上方的发射率较低，这是预期的。请注意，估计发射率的不同模型会导致 LST 值以及 SUHI 强度存在一些差异。

然后，我们将该发射率与亮度温度结合起来，使用简单的单通道算法计算每个像素的 LST，这是辐射传输方程的线性近似。

图6|重庆的地表发射率，基于植被比例。绿色像素显示较高的值，白色像素显示较低的值

// Calculate LST from emissivity and brightness temperature.
var lstLandsat = thermal.expression(
'(Tb/(1 + (0.001145* (Tb / 1.438))*log(Ep)))-273.15', {
'Tb': thermal.select('B10'),
'Ep': em.select('EMM')
}).updateMask(notWater);
Map.addLayer(lstLandsat, {
min: 25,
max: 35,
palette: ['blue', 'white', 'red'],
},
'LST_Landsat');