GEE学习之简单阈值分割

最近一直在做土地利用分类,尝试了非监督分类、随机森林,感觉制作标签好费事,就想通过光谱指数阈值偷懒。当然结果表明,这样分割只能看大致趋势!!实际分类结果不准!!对于一些趋势分析倒是勉强可以用用,而且确实省时省力。

选取典型土地利用类型,查看光谱特征

首先我在地图上选取了几种典型的土地利用类型,计算他们的NDVI、MNDWI、SI、NDBI等光谱指数,查看哪种光谱指数能分类得更好。结果发现,EVI和SAVI和NDVI结果相差不大,MNDWI效果没有NDWI效果好,IBI结果与NDBI效果差不多,随后就选择了NDVI、MNDWI、SI、NDBI这四种指标作为分割标准。

//将各土地利用类型合成一个featurecollection;table是导入的研究区边界
var regions = new ee.FeatureCollection([sand,water,crop,wetlands,city,road,barelands]);
Map.centerObject(table, 9);
Map.addLayer(ee.Image().paint(table, 0, 2), {}, 'roi');
var date = ee.DateRange('2022-06-01','2022-09-01')

//landsat影像处理与光谱指数计算
//product:LANDSAT/LT05/C02/T1_L2
function input_landsat8L2(region,date){
//缩放 
  function applyScaleFactors(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);
  }

//NDVI 
  function NDVI(image) { 
    return image.addBands( 
     image.normalizedDifference(["SR_B5", "SR_B4"]) 
           .rename("NDVI")); 
  } 

//EVI
  var EVI = function(image) {
   var evi =  image.expression(
       '2.5 * float(nir - red)/ (nir + 6*red -7.5*blue +1 )',
       {
         'nir': image.select('SR_B5'),
          'red': image.select('SR_B4'),
          'blue':image.select('SR_B2')
       }).rename('EVI');
    return image.addBands(evi);
  };

//NDWI 
  function NDWI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B3", "SR_B5"]) 
          .rename("NDWI")); 
  } 

//MNDWI 
  function MNDWI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B3", "SR_B6"]) 
          .rename("MNDWI")); 
  } 
  
//NDBI 
  function NDBI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B6", "SR_B5"]) 
          .rename("NDBI")); 
  } 

//SI
  var SI = function(image) {
   var si =  image.expression(
       'float(swir + red-nir-blue)/ (swir +red +blue +nir )',
       {
         'swir':image.select('SR_B6'),
         'nir': image.select('SR_B5'),
          'red': image.select('SR_B4'),
          'blue':image.select('SR_B2')
       }).rename('SI');
    return image.addBands(si);
  };
//SAVI
  var SAVI = function(image) {
   var savi =  image.expression(
       '1.5*float(nir-red)/ (red +nir+0.5 )',
       {
         'nir': image.select('SR_B5'),
          'red': image.select('SR_B4'),
       }).rename('SAVI');
    return image.addBands(savi);
  };
//DEM
  function DEM(image){
    var strm=ee.Image("USGS/SRTMGL1_003");
    var dem=ee.Algorithms.Terrain(strm);
    var elevation=dem.select('elevation');
    var slope=dem.select('slope');
    return image.addBands(elevation.rename("ELEVATION")
      ).addBands(slope.rename("SLOPE"))
}
//IBI
function IBI(image){
  var ibiA = image.expression('2 * SWIR1 / (SWIR1 + NIR)', {
        'SWIR1': image.select('SR_B6'),
        'NIR'  : image.select('SR_B5')
      }).rename(['IBI_A']);
     
  var ibiB = image.expression('(NIR / (NIR + RED)) + (GREEN / (GREEN + SWIR1))', {
        'NIR'  : image.select('SR_B5'),
        'RED'  : image.select('SR_B4'),
        'GREEN': image.select('SR_B3'),
        'SWIR1': image.select('SR_B6')
      }).rename(['IBI_B']);
  var ibiAB = ibiA.addBands(ibiB);
  var ibi = ibiAB.normalizedDifference(['IBI_A', 'IBI_B']);
  return image.addBands(ibi.rename("IBI"))
  }
//remove cloud by select bitwise
  function bitwiseExtract(value, fromBit, toBit) {
    if (toBit === undefined) toBit = fromBit;
    var maskSize = 1 + toBit - fromBit; //位数
    var mask = ee.Number(1 << maskSize).subtract(1); 
    return value.rightShift(fromBit).bitwiseAnd(mask);
  }

  function cloudfree_landsat(image){
    var qa = image.select('QA_PIXEL')
    var cloudState = bitwiseExtract(qa, 3)    //5
    var cloudShadowState = bitwiseExtract(qa, 4)  //3
    var mask = cloudState.eq(0) // Clear
    .and(cloudShadowState.eq(0)) // No cloud shadow
    return image.updateMask(mask)  
  }
  
  //load landsat
  var collection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterDate(date).filterBounds(region);//LANDSAT/LT05/C02/T1_L2
var preinput = collection//.filterMetadata('CLOUD_COVER','less_than', 3.0).sort('CLOUD_COVER', true) //这里其实并没有用到云量控制
              .map(cloudfree_landsat).map(applyScaleFactors) 
              .map(NDVI) 
              .map(NDWI) 
              .map(NDBI)
              .map(MNDWI)
              .map(EVI).map(DEM).map(SI).map(SAVI).map(IBI)//.map(nightlight).qualityMosaic('NDVI').clip(region) //.mosaic() can be changed in different ways
var bands = [ 
  "SR_B6", "SR_B2", "SR_B3","SR_B4","SR_B5", "SR_B7", 
  "NDBI",  "NDVI" ,'MNDWI','NDWI',"SAVI",'EVI','ELEVATION','SI','SLOPE','IBI'
];
//get the base map
var input = preinput.select(bands);  
print(input)
  Map.addLayer(preinput.mean().clip(region), {min:0, max:0.3, bands: ['SR_B4', 'SR_B3', 'SR_B2']} ,"rawscenecloud");
  return input
} 
var input = input_landsat8L2(table,date)

//均值合成
var mean = input.mean().clip(table)
var ndvi_crop = mean.select('NDVI').reduceRegion({  
   reducer: ee.Reducer.mean(),  
   geometry: crop.geometry(),   
   scale: 30  
}); 
print('ndvi_mean_crop:',ndvi_crop)
var ndvi_bare = mean.select('NDVI').reduceRegion({  
   reducer: ee.Reducer.mean(),  
   geometry: barelands.geometry(),   
   scale: 30  
});
print('ndvi_mean_bareland:',ndvi_bare)
var ndvi_wetland = mean.select('NDVI').reduceRegion({  
   reducer: ee.Reducer.mean(),  
   geometry: wetlands.geometry(),   
   scale: 30  
}); 
print('ndvi_mean_wetland:',ndvi_wetland)
var ndvi_sand = mean.select('SI').reduceRegion({  
   reducer: ee.Reducer.mean(),  
   geometry: sand.geometry(),   
   scale: 30  
}); 
print('ndvi_mean_sand:',ndvi_sand)

由于影像质量每年都不确定,因此先根据每年平均值大致判断阈值(这里我为了省事直接人为设定阈值了,所以结果受主观影响很大,最好设计一套更符合每年情况的阈值确定方法)。

通过阈值进行分类

这里涉及的总思路是:(1)逐年计算6-9月样本点(田地、湿地、裸地)平均NDVI,作为阈值参考;(2)通过MNDWI分割水体(由于MNDWI对于水体分割效果很好,设置固定阈值),通过设置NDVI阈值划分湿地(阈值逐年根据NDVI均值变化);(3)通过设置NDVI阈值进一步划分田地(阈值逐年根据NDVI均值变化,考虑到夏季水面植物生长,用湿地、水域mask作为约束),通过SI阈值划分沙地(阈值不变,但用田地mask作为约束);(4)通过NDVI阈值划分裸地(阈值不变,使用水域、沙地、田地作为约束);(5)通过NDBI值划分不透水面(阈值不变,使用裸地、沙地、水域作为约束)。在各类型相互约束前提下,通过目视大致判断土地利用分布,调整阈值(湿地、田地),将各个土地类型分别赋值,不属于此类的像素mask,然后合成一张分类图,再按照灌域进行像素统计。这里还是强调,采用的方法太依赖人工经验,之后有机会的话再认真考虑阈值的选取。

//湿地——NDVI
var wetland = mean.select('NDVI').gte(0.73)
var wetland_mask = wetland.mask(wetland)
var WETLAND = wetland_mask.where(wetland_mask.select('NDVI').eq(1),1).rename('lc').set(ee.Dictionary(['lc',1])).toInt8()
var otherwl = mean.select('NDVI').lt(0.73)

//水体——MNDWI
var waterbodies = mean.select('MNDWI').gte(0.2)
var waterbodies_mask = waterbodies.mask(waterbodies)
var WATER = waterbodies_mask.where(waterbodies_mask.select('MNDWI').eq(1),3).rename('lc').set(ee.Dictionary(['lc',3])).toInt8()
var otherw = mean.select('MNDWI').lt(0.2)

//农田——NDVI
var field = mean.select('NDVI').gte(0.4)//0.35
var field_mask = field.mask(field).updateMask(otherwl.mask(otherwl)).updateMask(otherw.mask(otherw))
var FIELD = field_mask.where(field_mask.select('NDVI').eq(1),2).rename('lc').set(ee.Dictionary(['lc',2])).toInt8()
var otherf = mean.select('NDVI').lt(0.4)

//沙地——SI
var sands = mean.select('SI').gte(0.13)
var sands_mask = sands.mask(sands).updateMask(otherw.mask(otherw)).updateMask(otherf.mask(otherf))
var SANDS = sands_mask.where(sands_mask.select('SI').eq(1),4).rename('lc').set(ee.Dictionary(['lc',4])).toInt8()
var others = mean.select('SI').lt(0.13)
//var others = mean.select('SI').mask(mean.select('SI').lt(0.12)).unmask(0);

//裸地——NDVI
var bareland = mean.select('NDVI').gte(0.2)
var bareland_mask = bareland.mask(bareland).updateMask(otherf.mask(otherf))
.updateMask(otherw.mask(otherw)).updateMask(others.mask(others))
var BARELAND = bareland_mask.where(bareland_mask.select('NDVI').eq(1),5).rename('lc').set(ee.Dictionary(['lc',5])).toInt8()
var otherb = mean.select('NDVI').lt(0.2)

//不透水面——NDBI
var impervious = mean.select('NDBI').gt(-0.1)
var impervious_mask = impervious.mask(impervious).updateMask(others.mask(others))
.updateMask(otherb.mask(otherb)).updateMask(otherw.mask(otherw))
var IMPERVIOUS = impervious_mask.where(impervious_mask.select('NDBI').eq(1),6).rename('lc').set(ee.Dictionary(['lc',6])).toInt8();

//加载各分类情况分布
Map.addLayer(wetland_mask,{palette:["006400"]},'wetland')
Map.addLayer(field_mask,{palette:["7FFF00"]},'field')
Map.addLayer(bareland_mask,{palette:["8B4513"]},'bareland')
Map.addLayer(impervious_mask,{palette:["8A8A8A"]},'impervious')
Map.addLayer(sands_mask,{palette:["FFF68F"]},'sands')
Map.addLayer(waterbodies_mask,{palette:["1E90FF"]},'waterbodies')

这里用到了mask,mask会将符合条件的像素赋值为1,不符合条件的像素赋值为0,然后将值为0的像素掩膜掉,这里要注意不是删除原像素,只是将其掩膜。updateMask可以应用其他的mask,还有一个unmask()的命令,表示将掩膜掉的像素赋值为括号内的数字,如-9999等。

在上面的阈值分割过程中,使用了多种光谱指数对6种土地利用类型进行分割,需要注意的是不同土地利用类型之间不能有重叠,因此除了强特征的土地类型(湿地、沙地、水体),余下的(农田、裸地、不透水面)在分割的过程中需要对已分好的类型进行掩膜处理。这一步很关键,不然会影响后面面积统计的结果。类型之间的相互约束真的很重要。

在分类结束后,通过where()函数将掩膜出来的像素赋予不同的数字作为土地类型的区分,并将波段重命名为lc(landcover)。set函数是将每个影像添加一个property,名字为lc, 这是之前失败的尝试,实际上可加可不加。

合并各掩膜,进行面积统计

各土地类型分类完成后,将其集合为imagecollection,此时每张影像的波段都是lc。集合后,使用mean()函数将各个影像合成为一张土地分类图。如果之前土地掩膜没有弄好的话,这里合成就会出大问题。之前跳错说type不统一,因此在上面rename后将其转换为int(toInt8()),之后就不报错了。

//合成为一张土地分类图
var classification = ee.ImageCollection.fromImages([WETLAND,FIELD,WATER,SANDS,BARELAND,IMPERVIOUS])
//print(classification2)
Map.addLayer(classification2.mean(),{min:1,max:6,palette:["006400","7FFF00","1E90FF","FFF68F","8B4513","8A8A8A"]},'mosaic')

//提取研究区的小区域名称
var point = table.aggregate_array('Name1')

//for循环统计各小区域的土地类型面积
for(var i = 0; i<6 ; i++) {  
  var roi = table.filterMetadata('Name1','equals',point.get(i));
//count area
  var areaImage = ee.Image.pixelArea().addBands(classification.mean())
  var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'type',
      }),
      geometry: roi,
      scale: 30,
      maxPixels: 1e13
      }); 
  var classAreas = ee.List(areas.get('groups'))
  var classAreaLists = classAreas.map(function(item) {
  var areaDict = ee.Dictionary(item)
  var classNumber = ee.Number(areaDict.get('type')).format()
  var area = ee.Number(
      areaDict.get('sum')).divide(1e6)//.round()
  return ee.List([classNumber, area])
    })
  var result1 = ee.Dictionary(classAreaLists.flatten())
  print('area:'+ i ,result1)
  //}
}

//统计整个区域的各土地利用类型面积
var image_area = ee.Image.pixelArea().addBands(classification.mean())
var Area = image_area.reduceRegion({
    'reducer': ee.Reducer.sum().group({
        'groupField': 1,
        'groupName': 'class'
    }),
    'geometry': table,
    'scale': 30,
    'maxPixels': 1e13 
  })
 
  var classAreas1 = ee.List(Area.get('groups'))
  var classAreaLists1 = classAreas1.map(function(item) {
    var areaDict1 = ee.Dictionary(item)
    var classNumber1 = ee.Number(areaDict1.get('class')).format()
    var area1 = ee.Number(
      areaDict1.get('sum')).divide(1e6)//.round()
    return ee.List([classNumber1, area1])
  })
  var result2 = ee.Dictionary(classAreaLists1.flatten())
  print(result2)

后记

这个原理虽然比较简单,但是也到处找了很久代码,在设置阈值时也做了很多的前处理分析,虽然最后表明人工阈值不太可行,但是在大致的趋势上面还是足够明显的。不过果然偷懒还是不太行,老老实实地回去打标签做随机森林分类了。。

 

 

评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值