先展示做出来的效果：
1.数据导入
以10年为间隔，同时考虑Landsat卫星的运行时间和型号，设置1986、1995、2005、2015和2022为采样年份，研究1986-2022年粤港澳大湾区的建成区的变化。
对于所有年份的Landsat数据，我们都使用地表反射率数据，并根据云量进行筛选，使用median函数进行去云，因为不是本次实验的重点，具体细节可以参考我写的这篇博客：http://t.csdn.cn/QeD1k。
var first_year = 1986;var second_year =1995;var third_year =2005;var fourth_year = 2015;var fifth_year =2022;//---------------------------Part 2 load study area and data-------------------------------------------------------//load the study areavar roi = table;Map.addLayer(roi, {}, 'roi',false);Map.centerObject(roi,8);// Applies scaling factors landsat 5.function applyScaleFactors5(image) {var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0);return image.addBands(opticalBands, null, true).addBands(thermalBand, null, true);}// Applies scaling factors for landsat 8.function applyScaleFactors8(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 landsat1 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(roi).filterDate('1986-01-01', '1986-12-28').map(applyScaleFactors5).filter(ee.Filter.lt('CLOUD_COVER',12)).median();var landsat2 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(roi).filterDate('1995-01-01', '1995-12-31').map(applyScaleFactors5).sort('CLOUD_COVER').filter(ee.Filter.lt('CLOUD_COVER',12)).median();var landsat3 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(roi).filterDate('2005-01-01', '2005-12-31').map(applyScaleFactors5).sort('CLOUD_COVER').filter(ee.Filter.lt('CLOUD_COVER',5)).median();var landsat4 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterBounds(roi).filterDate('2015-01-01', '2015-12-31').map(applyScaleFactors8).sort('CLOUD_COVER').filter(ee.Filter.lt('CLOUD_COVER',19)).median();var landsat5 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterBounds(roi).filterDate('2022-04-01', '2022-12-30').map(applyScaleFactors8).filter(ee.Filter.lt('CLOUD_COVER',2)).median();2.计算NDBI
本次实验比较偷懒，不需要考虑非常准确的提取效果，使用NDBI（建筑物指数）可以将不透水面提取出来，NDBI = (SWIR – NIR) / (SWIR + NIR)。基于多次尝试，设置NDBI阈值为：-0.15。下面放了一些关键步骤。其中的SR开头的波段都是地表反射率数据的波段，比如SR_B5就是地表反射率的第五波段。注意Landsat5和Landsat8计算的公式是不一样的。normalizedDifference函数就是计算归一化的各种指数的，你还会在计算NDVI的时候用到它。
解释一下代码：首先提取2022年的不透水面，为其他年份做一个掩膜，如果不做掩膜，则可能出现每一年城市的像元变化较大的问题。built5就是2022年不透水面二值化的结果，1为不透水面，0为其他。然后剩下的年份全用2022的结果做一个mask，就得到了剩下4年的不透水面。同时，这里假设不透水面像元不会转成其他像元（对应到现实生活中就是城市不会转成乡村），所以这里对1995，2005，2015和2022年都做了一个并运算（or函数），与前一年的结果进行合并，其实这样也减少了误差。后面6句话是为了后面绘制城市化过程做准备。将1986年不透水面设置为1，1995年不透水面设置为2，以此类推，然后将全部的其他像元都设置为6。这样将影像集合取min的时候，就可以得到每一次新增的不透水面。此时注意其他像元的值一定要设成6而不能设置成别的值，否则在出图的时候就会导致颜色存在不确定性，设置成6是因为想把出图所用到的数据锁定在1，2，3，4，5，6这6个连续整数之间，颜色就可以被完全确定下来，而不会被线性拉伸了（我之前检查了好久检查不出问题）。里面那个where函数其实就有点像是其他语言里面的if函数。
//-----------Part 3 detect urban renewal----------------------------//calculate ndbi// print(landsat1)var ndbi1 = landsat1.normalizedDifference(['SR_B5','SR_B4']);var ndbi2 = landsat2.normalizedDifference(['SR_B5','SR_B4']);var ndbi3 = landsat3.normalizedDifference(['SR_B5','SR_B4']);var ndbi4 = landsat4.normalizedDifference(['SR_B6','SR_B5']);var ndbi5 = landsat5.normalizedDifference(['SR_B6','SR_B5']);//set thresholdsvar ndbi_t = -0.15;var built5 = ndbi5.gt(ndbi_t).clipToCollection(roi);built5 = built5.mask(built5);var built1 = ndbi1.gt(ndbi_t).mask(built5);var built2 = ndbi2.gt(ndbi_t).mask(built5).or(built1);var built3 = ndbi3.gt(ndbi_t).mask(built5).or(built2);var built4 = ndbi4.gt(ndbi_t).mask(built5).or(built3);built5 = built5.or(built4);built1 = built1.where(built1.gt(0),1).where(built1.eq(0),9);built2 = built2.where(built2.gt(0),2).where(built2.eq(0),9);built3 = built3.where(built3.gt(0),3).where(built3.eq(0),9);built4 = built4.where(built4.gt(0),4).where(built4.eq(0),9);built5 = built5.where(built5.gt(0),5).where(built5.eq(0),9);var collection = ee.ImageCollection([built1,built2,built3,built4,built5])NDBI会将部分的裸土错分为城市像元，但问题不是很严重，毕竟城市内部的裸土不多（图1）。同时沿海地区和入海口泥沙较多的地方都容易被提取成城市像元，但是缩小阈值又提不全不透水面。比较好的方式还可以用别的指数筛选一遍水体，这可以作为改进方向。
出图的代码如下所示：首先把主图画出来，也就是前两句话。然后添加图例。这些代码都可以套用，只需要根据自己的需要改一下min和max的值和palette里面的颜色，还有调整图例的名称即可，具体的细节就不展开说了（因为我也不是很懂functio，基本都是copy别人的，然后自己小改一下嘿嘿）。
var result=collection.min();Map.addLayer(result,{'min':1,'max':6,'palette':['EA047E','FF6D28','FCE700','00F5FF','00C897','00000000']},'result');//----------------Add the legend!------ -----------//添加图例函数 function addLegend(palette, names) {//图例的底层Panelvar legend = ui.Panel({style: {position: 'bottom-right', padding: '5px 10px'}});//图例标题var title = ui.Label({value:'城市化年份',style: {fontWeight: 'bold', color: "red", fontSize: '16px'}});legend.add(title);//添加每一列图例颜色以及说明var addLegendLabel = function(color, name) {var showColor = ui.Label({style: {backgroundColor: '#' + color, padding: '8px', margin: '0 0 4px 0'}});var desc = ui.Label({value: name, style: {margin: '0 0 4px 8px'}});return ui.Panel({widgets: [showColor, desc], layout: ui.Panel.Layout.Flow('horizontal')});};//添加所有的图例列表for (var i = 0; i < palette.length; i++) {var label = addLegendLabel(palette[i], names[i]);legend.add(label);}Map.add(legend);}//color table & legend namevar palette = ['EA047E','FF6D28','FCE700','00F5FF','00C897'];var names = ["1986","1995","2005",'2015','2022'];//添加图例 addLegend(palette, names); 然后就可以做出一开始的那个图啦~
下一次就讲城市形态学计算！