目录
- 写在前面
- 1.构建物候特征
- 2.构建光谱特征
- 3.将所有影像合并为一幅影像
- 4.构建随机森林算法进行分类
- 5.算法的存储
- 6.面积统计
写在前面
前段时间因为考研的原因一直没能更新，已经完成了农作物种植结构的提取，现在给大家分享一下。
主要也是结合前面写过的Google Earth Engine（GEE）使用土地利用数据（modis）上采样Landsat数据提取农田范围，将以上结果作为研究的基础，结合物候特征，光谱特征，地形特征，选择随机森林算法进行农作物的提取。
不想看后面的同学可以直接看代码：https://code.earthengine.google.com/3ed8a5303c610e063ae3a4517ca4e146
1.构建物候特征
在这里主要通过不同农作物NDVI值的差异来选择合适的生长时间进行影像的选择。
通过构建典型地物的NDVI时序特征曲线，不同农作物随着时间的变化，NDVI值也在跟着变化，这种变化也反映了不同农作物在不同时期的光谱特征差异。因此，选择NDVI差异较大的越冬期和成熟期作为最佳识别期进行农作物分类，可以有效进行农作物种类的区分。
var winter = ee.ImageCollection('COPERNICUS/S2_SR') .filterDate('2021-11-30', '2022-01-01') .filterBounds(roi) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) .map(maskS2clouds) .map(mndwi) .map(ndbi) .map(ndvi) .select(['B1','B2','B3','B4','B8','B12','QA60','MNDWI','NDBI','NDVI']);var winter1=winter.median() .addBands(elevation.rename("ELEVATION")) .addBands(slope.rename("SLOPE")) .clip(roi);Map.addLayer(winter1, {bands: ['B4', 'B3', 'B2'],min:0, max: 0.3}, 'winter1',false);var summer = ee.ImageCollection('COPERNICUS/S2_SR') .filterDate('2022-03-01', '2022-04-01') .filterBounds(roi) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) .map(maskS2clouds) .map(mndwi) .map(ndbi) .map(ndvi) .select(['B1','B2','B3','B4','B8','B12','QA60','MNDWI','NDBI','NDVI']);var summer1=summer.median() .addBands(elevation.rename("sum_ELEVATION")) .addBands(slope.rename("sum_SLOPE")) .clip(roi);这里我把基础影像根据农田范围进行了裁剪，只留下农田区域，能有效的排除一些非农田地物类型。
2.构建光谱特征
也就是构建NDVI、NDWI、NDBI等光谱指数，来区分植被、建筑，水体等。
//构建光谱特征function mndwi(image){ return image.addBands(image.normalizedDifference(['B3', 'B8']).rename('MNDWI'))}function ndbi(image){ return image.addBands(image.normalizedDifference(['B12', 'B8']).rename('NDBI'))}function ndvi(image){ return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'))}3.将所有影像合并为一幅影像
将上面选择的越冬期和成熟期的多幅影像合并为一幅，在这里需要更改波段名，不然会出现波段名重复的情况。
//更改波段名称,合成一幅影像var useimage = winter1.addBands([ summer1.select('B2').rename('sum_B2'), summer1.select('B3').rename('sum_B3'), summer1.select('B4').rename('sum_B4'), summer1.select('B8').rename('sum_B8'), summer1.select('B12').rename('sum_B12'), summer1.select('MNDWI').rename('sum_MNDWI'), summer1.select('NDBI').rename('sum_NDBI'), summer1.select('NDVI').rename('sum_NDVI'), summer1.select('sum_ELEVATION'), summer1.select('sum_SLOPE') ]) print(useimage,'useimage')//求交集，对影像进行掩膜var useimage1 = useimage.updateMask(esamask)// var useimage1 = useimageMap.addLayer(useimage1, {bands: ['sum_B4', 'sum_B3', 'sum_B2'],min:0, max: 0.3}, 'useimage1');var usebands = ['B2','B3','B4','B8','B12','MNDWI','NDBI','NDVI','ELEVATION','SLOPE', 'sum_B2','sum_B3','sum_B4','sum_B8','sum_B12','sum_MNDWI','sum_NDBI','sum_NDVI','sum_ELEVATION','sum_SLOPE']4.构建随机森林算法进行分类
构建算法也就是一些老套路了，无非就是建立样本，将样本分成训练样本和验证样本。
在这里想说一点，随机森林的决策树棵树选择不能随便选择，需要进行不同棵数的尝试，可以使用以下代码进行分析，选择准确率最高的颗数进行算法构建。
//选取森林棵树 var numTrees = ee.List.sequence(5, 50, 5); var accuracies = numTrees.map(function(t) { var classifier = ee.Classifier.smileRandomForest(t) .train({ features: trainingPartition, classProperty: 'landcover', inputProperties: usebands }); return testingPartition .classify(classifier) .errorMatrix('landcover', 'classification') .accuracy(); }); print(ui.Chart.array.values({ array: ee.Array(accuracies), axis: 0, xLabels: numTrees }));得到的结果大概是这样，选择准确率最高的就可以了。
5.算法的存储
因为我们上面已经训练好了算法，所以这个算法是可以运用在其他年份的，我们将它保存在个人assets里，当缺少其他年份的样本时，可以使用该算法进行分类，也就是算法复用。
var trees = ee.List(ee.Dictionary(testclassifier.explain()).get('trees'))var dummy = ee.Feature(roi.geometry())var col = ee.FeatureCollection(trees.map(function(x){return dummy.set('tree',x)}))print(col)Export.table.toAsset(col,'save_classifier','projects/dyb/assets/2022chuzhou')6.面积统计
上面已经做好了农作物的分类，接下来可以对农作物的面积做一个统计。
var areaImage = ee.Image.pixelArea().addBands(classified) var clsdArea = areaImage.reduceRegion({ 'reducer': ee.Reducer.sum().group({ 'groupField': 1, 'groupName': 'class' }), 'geometry': roi.geometry(), 'scale': 10, 'maxPixels': 1e13 }) print('Kerala landcover clsArea:', clsdArea.getInfo())最后把所有代码放在这里
https://code.earthengine.google.com/3ed8a5303c610e063ae3a4517ca4e146
