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

采样点的选取

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

合成多季节多波段影像

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