先展示做出来的效果：
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 area
var 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 thresholds
var 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) {
//图例的底层Panel
var 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 name
var palette = ['EA047E','FF6D28','FCE700','00F5FF','00C897'];
var names = ["1986","1995","2005",'2015','2022'];
//添加图例
addLegend(palette, names); 然后就可以做出一开始的那个图啦~
下一次就讲城市形态学计算！
文章出处登录后可见！
