不透水地表（Google Earth Engine(GEE)城市不透水面提取，NDBI）
先展示做出来的效果：
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);然后就可以做出一开始的那个图啦~
下一次就讲城市形态学计算！
创心域SEO版权声明：以上内容作者已申请原创保护，未经允许不得转载，侵权必究！授权事宜、对本内容有异议或投诉，敬请联系网站管理员，我们将尽快回复您，谢谢合作！