利用GEE对季节性地物进行分类的代码实现

爱被打了一巴掌 2024-02-20 13:02 45阅读 0赞

### 采样点的选取 ###

如果你采用监督学习的话，那就手动打标签  
或者可以了解一下非监督学习  
![在这里插入图片描述][1291067ece504aad88e1e94cc9cd4335.png]

### 合成多季节多波段影像 ###

首先，制作一个包含多波段的影像，每个波段作为随机森林分类器的一个feature输入，提升feature的丰富度以保证分类精度。

**1、对landsat5用的云掩膜函数：**

// cloud mask
    var cloudMaskL457 = function(image) {
        
      var qa = image.select('pixel_qa');
      // If the cloud bit (5) is set and the cloud confidence (7) is high
      // or the cloud shadow bit is set (3), then it's a bad pixel.
      var cloud = qa.bitwiseAnd(1 << 5)
                      .and(qa.bitwiseAnd(1 << 7))
                      .or(qa.bitwiseAnd(1 << 3));
      // Remove edge pixels that don't occur in all bands
      var mask2 = image.mask().reduce(ee.Reducer.min());
      return image.updateMask(cloud.not()).updateMask(mask2);
    };

**2、构建正向指标和反指标**

// indicators
    // landsat 5
    function NDWI(img){
        
      var nir = img.select('B4');
      var green = img.select('B2');
      var ndwi = img.expression(
        '(B2-B4)/(B2+B4)',
        {
        
          'B4':nir,
          'B2':green
        });
        return ndwi;
    }
    
    function MNDWI(img){
        
      var swir1 = img.select('B5');
      var green = img.select('B2');
      var mndwi = img.expression(
        '(B2-B5)/(B2+B5)',
        {
        
          'B5':swir1,
          'B2':green
        });
        return mndwi;
    }
    
    function EWI(img){
        
      var swir1 = img.select('B5')
      var nir = img.select('B4');
      var green = img.select('B2');
      var ewi = img.expression(
        '(B2-B4-B5)/(B2+B4+B5)',
        {
        
          'B5':swir1,
          'B4':nir,
          'B2':green
        });
        return ewi;
    }
    
    function NDBI(img){
        
      var swir1 = img.select('B5');
      var nir = img.select('B4');
      var ndbi = img.expression(
        '(B5-B4)/(B5+B4)',
        {
        
          'B5':swir1,
          'B4':nir
        });
        return ndbi;
    }
    
    function NDVI(img){
        
      var nir = img.select('B4');
      var red = img.select('B3');
      var ndvi = img.expression(
        '(B4-B3)/(B4+B3)',
        {
        
          'B4':nir,
          'B3':red
        });
        return ndvi;
    }

**3、选取我们的影像数据源，这里用的是landsat5，如果更换数据源的话云掩膜函数和指标计算也要一并调整。影像源代码如下：**

// satellite data
    var l8_col = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR");

**4、开始合成季节影像。**  
这里我做的分类是1990年的，但是因为landsat的重返周期太长了，尺度稍微大点的话会碰上很多云，掩膜的话又都掩膜没了，所以我使用1989到1991三年的合成，其中冬是12月1日到3月1日、春是3月1日到6月1日、夏是6月1日到9月1日、秋是9月1日到12月1日，可以根据自己的需求调整。另外，为了避免云掩膜出现的孔洞，所以我用unmask函数把孔洞赋值为0以填补出现的孔洞。具体代码如下：

// winter
    var img_winter = ee.List([])
    for (var i=1989; i<1991; i++){
        
      var winter = ee.Image(l8_col.filterBounds(roi)
                                  .filterDate(i+"-12-1", (i+1)+"-3-1")
                                  .filter(ee.Filter.lt('CLOUD_COVER', 100))
                                  .map(cloudMaskL457)
                                  .mean());
    }
    img_winter = img_winter.add(winter);
    img_winter = ee.ImageCollection(img_winter).mean().unmask(0);
    // spring
    var img_spring = ee.List([])
    for (var i=1989; i<1991; i++){
        
      var spring = ee.Image(l8_col.filterBounds(roi)
                                  .filterDate(i+"-3-1", i+"-6-1")
                                  .filter(ee.Filter.lt('CLOUD_COVER', 100))
                                  .map(cloudMaskL457)
                                  .mean());
    }
    img_spring = img_spring.add(spring);
    img_spring = ee.ImageCollection(img_spring).mean().unmask(0);
    // summer
    var img_summer = ee.List([])
    for (var i=1989; i<1991; i++){
        
      var summer = ee.Image(l8_col.filterBounds(roi)
                                  .filterDate(i+"-6-1", i+"-9-1")
                                  .filter(ee.Filter.lt('CLOUD_COVER', 100))
                                  .map(cloudMaskL457)
                                  .mean());
    }
    img_summer = img_summer.add(summer);
    img_summer = ee.ImageCollection(img_summer).mean().unmask(0);
    // fall
    var img_fall = ee.List([])
    for (var i=1989; i<1991; i++){
        
      var fall = ee.Image(l8_col.filterBounds(roi)
                                  .filterDate(i+"-9-1", i+"-12-1")
                                  .filter(ee.Filter.lt('CLOUD_COVER', 100))
                                  .map(cloudMaskL457)
                                  .mean());
    }
    img_fall = img_fall.add(fall);
    img_fall = ee.ImageCollection(img_fall).mean().unmask(0);

**5、为了确保分类成功率，我引入GEE中现成的几个产品。**  
这些产品可以根据需要灵活调整，这里就不多介绍了。依然地，使用unmask函数填补掩膜孔洞，代码如下：

// referring image data
    var impervious = ee.Image("Tsinghua/FROM-GLC/GAIA/v10").clip(roi).unmask(-1);
    var water = ee.Image('JRC/GSW1_2/GlobalSurfaceWater').clip(roi).unmask(-1);

**6、之前合成了landsat5季节影像，别忘记上面我们写下了各种指数**  
这里就要使用这些合成影像计算指数生成各指数的波段，并进行命名，代码如下：

var ndwi_wi = NDWI(img_winter).rename('ndwi_wi');
    var mndwi_wi = MNDWI(img_winter).rename('mndwi_wi');
    var ewi_wi = EWI(img_winter).rename('ewi_wi');
    var ndbi_wi = NDBI(img_winter).rename('ndbi_wi');
    var ndvi_wi = NDVI(img_spring).rename('ndvi_wi');
    var ndwi_sp = NDWI(img_spring).rename('ndwi_sp');
    var mndwi_sp = MNDWI(img_spring).rename('mndwi_sp');
    var ewi_sp = EWI(img_spring).rename('ewi_sp');
    var ndbi_sp = NDBI(img_spring).rename('ndbi_sp');
    var ndvi_sp = NDVI(img_spring).rename('ndvi_sp');
    var ndwi_su = NDWI(img_summer).rename('ndwi_su');
    var mndwi_su = MNDWI(img_summer).rename('mndwi_su');
    var ewi_su = EWI(img_summer).rename('ewi_su');
    var ndbi_su = NDBI(img_summer).rename('ndbi_su');
    var ndvi_su = NDVI(img_summer).rename('ndvi_su');
    var ndwi_fa = NDWI(img_fall).rename('ndwi_fa');
    var mndwi_fa = MNDWI(img_fall).rename('mndwi_fa');
    var ewi_fa = EWI(img_fall).rename('ewi_fa');
    var ndbi_fa = NDBI(img_fall).rename('ndbi_fa');
    var ndvi_fa = NDVI(img_fall).rename('ndvi_fa');

**7、上面的指数是作为波段的形式。然后我们把引用的产品中的波段也提出并进行重命名，代码如下：**

var imperchange = impervious.select('change_year_index').rename('imperchange');
    var waterocc = water.select('occurrence').unmask(-1).rename('waterocc');
    var waterchange = water.select('change_norm').unmask(-150).rename('waterchange');

**8、然后，我们从每幅季节合成影像中提取波段1-7，并且加入同一幅影像，代码如下：**

// add bands into images
    var image = img_winter.addBands([ndwi_wi, mndwi_wi, ewi_wi, ndbi_wi, ndvi_wi,
                                      img_spring.select('B1').rename('B1_sp'), img_spring.select('B2').rename('B2_sp'), img_spring.select('B3').rename('B3_sp'), img_spring.select('B4').rename('B4_sp'), img_spring.select('B5').rename('B5_sp'),
                                      ndwi_sp, mndwi_sp, ewi_sp, ndbi_sp, ndvi_sp,
                                      img_summer.select('B1').rename('B1_su'), img_summer.select('B2').rename('B2_su'), img_summer.select('B3').rename('B3_su'), img_summer.select('B4').rename('B4_su'), img_summer.select('B5').rename('B5_su'),
                                      ndwi_su, mndwi_su, ewi_su, ndbi_su, ndvi_su,
                                      img_fall.select('B1').rename('B1_fa'), img_fall.select('B2').rename('B2_fa'), img_fall.select('B3').rename('B3_fa'), img_fall.select('B4').rename('B4_fa'), img_fall.select('B5').rename('B5_fa'),
                                      ndwi_fa, mndwi_fa, ewi_fa, ndbi_fa, ndvi_fa,
                                      imperchange, waterocc, waterchange]);

### 分类器参数设置 ###

上面我们合成了一个含有51波段的影像，我们现在就需要选取我们分类器需要输入的波段了，代码如下：

// select bands
    var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7',
                'ndwi_wi', 'mndwi_wi', 'ewi_wi', 'ndbi_wi', 'ndvi_wi',
                'B1_sp', 'B2_sp', 'B3_sp', 'B4_sp', 'B5_sp',
                'ndwi_sp', 'mndwi_sp', 'ewi_sp', 'ndbi_sp', 'ndvi_sp',
                'B1_su', 'B2_su', 'B3_su', 'B4_su', 'B5_su',
                'ndwi_su', 'mndwi_su', 'ewi_su', 'ndbi_su', 'ndvi_su',
                'B1_fa', 'B2_fa', 'B3_fa', 'B4_fa', 'B5_fa',
                'ndwi_fa', 'mndwi_fa', 'ewi_fa', 'ndbi_fa', 'ndvi_fa',
                'imperchange', 'waterocc', 'waterchange'];

用merge函数把各分类合成一整个集合，然后它们的波段是value，然后把训练集和测试集按90%和10%分开。代码如下：

// set training data
    var SP = Water.merge(Ponds).merge(Mud).merge(Others);
    var training = image.select(bands).sampleRegions({
        
      collection: SP,
      properties:['value'],
      scale: 30
    });
    // set varified data
    var withRandom = training.randomColumn('random');
    var split = 0.9;
    var trainingP = withRandom.filter(ee.Filter.lt('random', split)); // 90% training data
    var testingP = withRandom.filter(ee.Filter.gte('random', split)); // 10% verified data

### 监督分类的实现与分类精度计算 ###

得到一个分类后的影像

// classification
    var class_img = image.select(bands).classify(classifier);

光有分类结果不行，还需要知道分类精度。把测试集也分类一下，然后用GEE自带的函数求转移矩阵并且计算overall accuracy和kappa，这些会打印在console里，代码如下：

// testing
    var test = testingP.classify(classifier);
    
    // confusion matrix
    var confusionMatrix = test.errorMatrix('value', 'classification')
    print('confusionMatrix',confusionMatrix);//面板上显示混淆矩阵
    print('overall accuracy', confusionMatrix.accuracy());//面板上显示总体精度
    print('kappa accuracy', confusionMatrix.kappa());//面板上显示kappa值

最后，当然要显示季节影像和分类后的影像，RGB波段组合还是色带什么的可以根据需要自行调整。如果选择手点的话，可以先随便点几个点，然后再根据显示的季节影像和分类影像再增加数据集，慢慢达到精度，代码如下：

// show images
    var class_color = {
        
      min: 0,
      max: 3,
      palette: ['ffffff', '0049ff', '00ffca', '97422c']
    };
    
    Map.addLayer(img_winter.clip(roi.geometry()), {
        bands: ["B4", "B3", "B2"], min:0, max:2000}, "winter");
    Map.addLayer(img_spring.clip(roi.geometry()), {
        bands: ["B4", "B3", "B2"], min:0, max:2000}, "spring");
    Map.addLayer(img_summer.clip(roi.geometry()), {
        bands: ["B4", "B3", "B2"], min:0, max:2000}, "summer");
    Map.addLayer(img_fall.clip(roi.geometry()), {
        bands: ["B4", "B3", "B2"], min:0, max:2000}, "fall");
    Map.addLayer(class_img.clip(roi.geometry()), class_color, 'classifed')
    
    Map.centerObject(roi, 7)

最后是保存在云盘里，下载

Export.image.toDrive({
        
            image:  class_img,//分类结果
            description: 'test_season_1990',//文件名
            folder: '111',
            scale: 30,//分辨率
            region: roi,//区域
            crs:"EPSG:4326",
            maxPixels:1e13//此处值设置大一些，防止溢出
          });

[1291067ece504aad88e1e94cc9cd4335.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/02/20/fe9925efb1334ee7abd86bb8bd6e415e.png