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

Twcat_tree

发布于 2023-11-04 09:04:14

3070

发布于 2023-11-04 09:04:14

文章被收录于专栏：二猫の家二猫の家

采样点的选取

如果你采用监督学习的话，那就手动打标签或者可以了解一下非监督学习

合成多季节多波段影像

首先，制作一个包含多波段的影像，每个波段作为随机森林分类器的一个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//此处值设置大一些，防止溢出
      });

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2023-11-03，如有侵权请联系 cloudcommunity@tencent.com 删除

rename