导读：
热岛效应（Urban Heat Island Effect），是指一个地区的气温高于周围地区的现象，用两个代表性测点的气温差值（即热岛强度）表示，主要有城市热岛效应和青藏高原热岛效应两种。
可以根据卫星影像产品或者卫星热辐射波段计算热岛强度。
主要内容：
- 了解如何推导地表温度。
- 了解如何生成城市和农村参考。
- 知道如何计算地表城市热岛强度。
Section 1 ： 从Modis 获取地表温度LST （1km分辨率）现成的产品
1.MODIS/061/MYD11A2 产品
MYD11A2 V6.1 产品提供1200 x 1200公里网格内平均8天的地表温度(LST)。
MYD11A2中的每个像素值是在8天内收集的所有相应MYD11A1 LST像素的简单平均值。MYD11A2对所有每日LST值进行简单的平均，而不需要对特定的QA位进行任何过滤。每个MYD11A2的QA值都是基于任何给定像素的大多数输入每日QA值来设置的。之所以选择8天的合成周期，是因为Terra和Aqua平台的地面轨迹重复周期是8天的两倍。在本产品中，除了昼夜地表温度波段及其质量指标层外，还有MODIS波段31、32和8个观测层。
赤道过境时间约为下午1:30和凌晨1:30。
2. 代码
//========================================ROI=======================================
var regionInt = table;
//50米缓冲区
var regionInt = regionInt.union(50);
// Center the map on this area.
Map.centerObject(regionInt, 10);
// 设置地图显示选项指定了地图的底图样式。
//'HYBRID' 选项将地图显示为混合模式，即卫星图像和地图图层的组合.
//结合了卫星影像和地理信息，如道路、地名等
Map.setOptions('HYBRID');
Map.addLayer(regionInt, {}, 'ROI HAIkou');
//=======================================Modis LST====================================
var modisLst = ee.ImageCollection('MODIS/006/MYD11A2');
//选择感兴趣的波段
var landSurfTemperature = modisLst.select('LST_Day_1km');
//定义夏天时间范围6月1日（第152天）-8月31日 （第243天）
var sumFilter = ee.Filter.dayOfYear(152, 243);
//筛选14-18年5年间夏天的数据
var lstDateInt = landSurfTemperature
.filterDate('2014-01-01', '2019-01-01' )
.filter(sumFilter);
//取集合中所有图像的像素平均值。5年均值
var lstMean = lstDateInt.mean();
//将每个像素乘以缩放因子得到LST值。
var lstFinal = lstMean.multiply(0.02);
//mask out 除去所有的水像元(水的高比热容会影响地表温度)
//使用全球地表水数据集
var water = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select( 'occurrence');
var notWater = water.mask().not();
//将数据clip到感兴趣的区域，将开氏温度转换为摄氏度，并mask水像素。
var lstHaiKou = lstFinal.clip(regionInt).subtract(273.15) .updateMask(notWater);
Map.addLayer(lstHaiKou, {
palette: ['blue', 'white', 'red'],
min: 25,
max: 38
},
'LST_MODIS');
海口：
Section 2 ： 从Landsat 获取地表温度LST （60-120m分辨率） 需要通过热辐射波段计算
热辐射波段与材料的热辐射率相关。材料的热辐射率是指表征材料在温度T下热辐射能力的物理量。热辐射功率是指材料在一定温度下，单位时间内，单位表面积的热辐射强度。辐射率定义为在一定温度下该材料的热辐射功率与相同温度下黑体热辐射功率之比。实际材料的辐射率都不会大于1。如果已知材料的辐射率，通过测量热辐射功率就可以得知该物体所处的温度，这就是无接触式辐射测温的基本原理。
热辐射系数和热辐射率虽然相关，但它们描述的是不同的物理量。热辐射系数通常指的是材料的导热性能，而热辐射率则更专注于材料在特定温度下的热辐射能力。两者虽然都涉及到材料的热性能，但应用场景和计算方法有所不同。
因为卫星捕获的热辐射是地表温度和ε（surface emissivity）的函数，需要准确地规定或估计ε才能得到正确的LST。
2.1 几个概念
亮度温度BT是假设地球是一个黑体，从大气层顶部逸出的红外辐射的温度当量。它与地表温度不同。
地表温度LST需要考虑大气吸收和再发射，以及地表的发射率。地表温度(LST)是指地面的温度，太阳的热能被辐射到达地面后，一部分被反射，一部分被地面吸收，使地面增热，对地面的温度进行测量后得到的温度就是地表温度。地表温度是区域和全球尺度地表物理过程中的一个关键因子，也是研究地表和大气之间物质交换和能量交换的重要参数。许多应用如干旱、高温、林火、地质、水文、植被监测，全球环流和区域气候模型等都需要获得地表温度。
遥感反演地表温度的方法大致可归纳为5大类：单通道方法、分裂窗（双通道）方法、多通道温度-比辐射率分离方法、多角度温度反演方法和多角度与多通道相结合的方法。
2.2 代码
//========================================ROI=======================================
var regionInt = table;
//50米缓冲区
var regionInt = regionInt.union(50);
// Center the map on this area.
Map.centerObject(regionInt, 10);
// 设置地图显示选项指定了地图的底图样式。
//'HYBRID' 选项将地图显示为混合模式，即卫星图像和地图图层的组合.
//结合了卫星影像和地理信息，如道路、地名等
Map.setOptions('HYBRID');
Map.addLayer(regionInt, {}, 'ROI HAIkou');
//===============================================================================
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);
}
//定义夏天时间范围6月1日（第152天）-8月31日 （第243天）
var sumFilter = ee.Filter.dayOfYear(152, 243);
var col = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')
.filterBounds(regionInt)
.filterDate('2014-01-01', '2019-01-01')
.filter(sumFilter)
.map(cloudMask);
print('Landsat8 collection', col);
var image = col.median();
//mask out 除去所有的水像元(水的高比热容会影响地表温度)
//使用全球地表水数据集
var water = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select( 'occurrence');
var notWater = water.mask().not();
// 选择波段10(亮度温度Brightness temperature BT)
var thermal = image.select('B10')
.clip(regionInt)
.updateMask(notWater);
Map.addLayer(thermal, {
min: 280,
max: 310,
palette: ['blue', 'white', 'red']
},
'Landsat_BT');
//==================================计算NDVI=============================================
//由Landsat surface reflectance数据集计算NDVI
var ndvi = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(regionInt)
.filterDate('2014-01-01', '2019-01-01')
.filter(sumFilter)
.median()
.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI')
.clip(regionInt)
.updateMask(notWater);
Map.addLayer(ndvi, {
min: 0,
max: 1,
palette: ['blue', 'white', 'green']
},
'ndvi');
//归一化,寻找最大最小值
var minMax = ndvi.reduceRegion({
reducer: ee.Reducer.min().combine({ reducer2: ee.Reducer.max(), sharedInputs: true }),
geometry: regionInt,
scale: 30,
maxPixels: 1e9
});
print('minMax', minMax);
var min = ee.Number(minMax.get('NDVI_min'));
var max = ee.Number(minMax.get('NDVI_max'));
var fv = ndvi.subtract(min).divide(max.subtract(min)).rename('FV');
Map.addLayer(fv, {
min: 0,
max: 1,
palette: ['blue', 'white', 'green']
}, 'fv');
//==================================计算Emissivity=============================================
//使用基于这种植被覆盖度的发射率经验模型
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');
//==================================计算LST =============================================
//将该发射率与亮度温度结合起来，使用简单的单通道算法计算每个像素的LST，
//该算法是辐射传递方程的线性化近似
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: 38,
palette: ['blue', 'white', 'red'],
},
'LST_Landsat');
Section 3 ： 从GEE Landsat LST Toolbox 获取地表温度LST
使用为此目的开发的GEE 模块来计算地表温度。
//========================================ROI=======================================
var regionInt = table;
//50米缓冲区
var regionInt = regionInt.union(50);
// Center the map on this area.
Map.centerObject(regionInt, 10);
// 设置地图显示选项指定了地图的底图样式。
//'HYBRID' 选项将地图显示为混合模式，即卫星图像和地图图层的组合.
//结合了卫星影像和地理信息，如道路、地名等
Map.setOptions('HYBRID');
Map.addLayer(regionInt, {}, 'ROI HAIkou');
// Link to the module that computes the Landsat LST.
var landsatLST = require(
'projects/gee-edu/book:Part A - Applications/A1 - Human Applications/A1.5 Heat Islands/modules/Landsat_LST.js');
// Select region of interest, date range, and Landsat satellite.
var geometry = regionInt.geometry();
var satellite = 'L8';
var dateStart = '2014-01-01';
var dateEnd = '2019-01-01';
var useNdvi = true;
// Get Landsat collection with additional necessary variables.
var landsatColl = landsatLST.collection(satellite, dateStart, dateEnd, geometry, useNdvi);
//定义夏天时间范围6月1日（第152天）-8月31日 （第243天）
var sumFilter = ee.Filter.dayOfYear(152, 243);
//mask out 除去所有的水像元(水的高比热容会影响地表温度)
//使用全球地表水数据集
var water = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select( 'occurrence');
var notWater = water.mask().not();
// Create composite, clip, filter to summer, mask, and convert to degree Celsius.
var landsatComp = landsatColl
.select('LST')
.filter(sumFilter)
.median()
.clip(regionInt)
.updateMask(notWater)
.subtract(273.15);
Map.addLayer(landsatComp, {
min: 25,
max: 38,
palette: ['blue', 'white', 'red']
},
'LST_SMW');
以上3种方法算出来的结果有些差异，但从图像来看感觉差异还比较大。
Section 4 定义城市和农村参考
有了各种产品和算法对地表温度的估计，我们可以计算农村地表温度，然后减去城市地表温度，得到SUHI强度。
估算农村参考温度的方法有很多,在计算SUHI时，获得农村参考的最简单也是最常用的方法是在城市边界周围生成一个缓冲区。
//========================================ROI=======================================
var regionInt = table;
//50米缓冲区
var regionInt = regionInt.union(50);
// Center the map on this area.
Map.centerObject(regionInt, 10);
// 设置地图显示选项指定了地图的底图样式。
//'HYBRID' 选项将地图显示为混合模式，即卫星图像和地图图层的组合.
//结合了卫星影像和地理信息，如道路、地名等
Map.setOptions('HYBRID');
Map.addLayer(regionInt, {}, 'ROI HAIkou');
//====================================生成农村参考===========================================
// 方法一：缓冲固定的宽度，对于大城市来说，这种方法可能不太合适。
function bufferSubtract(feature) {
return ee.Feature(feature.geometry()
.buffer(2000)
.difference(feature.geometry()));
}
var ruralRef = regionInt.map(bufferSubtract);
Map.addLayer(ruralRef, {
color: 'green'
}, 'Buffer_ref');
print('method1:ruralRef', ruralRef);
//方法2: 使用迭代方法。该方法(请参阅下面的代码块)首先使用函数计算几何形状周围不同宽度的缓冲区，
//然后选择大小最接近原始几何形状的缓冲区域。
//定义缓冲区宽度序列
var buffWidths = ee.List.sequence(30, 3000, 30);
//定义函数生成标准化缓冲区
function bufferOptimize(feature) {
function buff(buffLength) {
var buffedPolygon = ee.Feature(feature.geometry()
.buffer(ee.Number(buffLength)))
.set({
'Buffer_width': ee.Number(buffLength)
});
var area = buffedPolygon.geometry().difference(feature .geometry()).area();
var diffFeature = ee.Feature( buffedPolygon.geometry().difference(feature .geometry()));
return diffFeature.set({
'Buffer_diff': area.subtract(feature.geometry() .area()).abs(),
'Buffer_area': area,
'Buffer_width': buffedPolygon.get('Buffer_width')
});
}
var buffed = ee.FeatureCollection(buffWidths.map(buff));
var sortedByBuffer = buffed.sort({
property: 'Buffer_diff'
});
var firstFeature = ee.Feature(sortedByBuffer.first());
return firstFeature.set({
'Urban_Area': feature.get('Area'),
'Buffer_width': firstFeature.get('Buffer_width')
});
}
var ruralRefStd = regionInt.map(bufferOptimize);
Map.addLayer(ruralRefStd, {
color: 'brown'
}, 'Buffer_ref_std');
print('method2:ruralRefStd', ruralRefStd);
//方法3: 根据土地利用数据
注：第二种方法运算量太大，建议单独运行，然后导出结果，再放到热岛效应的计算中去。
方法1的结果：
Section 5 计算地表城市热岛强度
暂时用不到不做了解。
