Google Earth Engine(GEE)农作物种植结构提取

Posted 星星和你都在我怀里

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Google Earth Engine(GEE)农作物种植结构提取相关的知识,希望对你有一定的参考价值。

目录

写在前面

前段时间因为考研的原因一直没能更新,已经完成了农作物种植结构的提取,现在给大家分享一下。
主要也是结合前面写过的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 = useimage
Map.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

以上是关于Google Earth Engine(GEE)农作物种植结构提取的主要内容,如果未能解决你的问题,请参考以下文章

Google Earth Engine(GEE)——Tasks升级,实现RUN ALL可以一键下载任务类型中的所有影像

Google Earth Engine基础使用方法

Google Earth Engine基础使用方法

Google Earth Engine 批量点击RUN任务

QGIS配置Google Earth Engine插件问题解决

QGIS配置Google Earth Engine插件问题解决