Google Earth Engine（GEE）——利用sentinel-2数据进行农作物提取分析

我们利用哨兵数据给农作物进行分类提取，主要有用到得时间节点是春夏秋冬四个季节，通过阈值法和updateMask来实现对农作物得提取。

最近发现了一个巨牛的人工智能学习网站，点击跳转到网站：前言 – 床长人工智能教程

代码：

//加载哨兵2号数据
var s2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED"),
    roi = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[114.27025202404155, 33.89965438737615],
          [114.27025202404155, 33.27732637743472],
          [115.25352839122905, 33.27732637743472],
          [115.25352839122905, 33.89965438737615]]], null, false);

Map.centerObject(roi)
Map.addLayer(roi)

//进行QA60去云操作函数
function maskS2clouds(image) {
  var qa = image.select('QA60')
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0))
  return image.updateMask(mask).divide(10000)
      .select("B.*")
      .copyProperties(image, ["system:time_start"])
}


//定义春夏秋冬的四季变换的时间节点
//春
var startDate_c = ee.Date('2021-02-03'); 
var endDate_c = ee.Date('2021-05-04'); 
//夏
var startDate_x = ee.Date('2021-05-05'); 
var endDate_x = ee.Date('2021-08-06'); 
//秋
var startDate_q = ee.Date('2021-08-07'); 
var endDate_q = ee.Date('2021-11-06'); 
//冬
var startDate_d = ee.Date('2021-11-07'); 
var endDate_d = ee.Date('2022-02-03'); 

//先将影像进行研究区过滤然后分别对不同的时间节点去筛选和去云操作
//这里代码可以优化，分别将筛选的四季的时间加入进来，然后直接进行整体的NDVI和NDWI的map
//在map完之后就可以直接进行。select().qualityMosaic()进行高质量镶嵌了
var collection = s2.filterBounds(roi)//位置过滤
var image_c=collection.filterDate(startDate_c, endDate_c)
                      .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                      .map(maskS2clouds)
                      .mosaic().clip(roi);
var image_x=collection.filterDate(startDate_x, endDate_x)
                      .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                      .map(maskS2clouds)
                      .mosaic().clip(roi);
var image_q=collection.filterDate(startDate_q, endDate_q)
                      .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                      .map(maskS2clouds)
                      .mosaic().clip(roi);
var image_d=collection.filterDate(startDate_d, endDate_d)
                      .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                      .map(maskS2clouds)
                      .mosaic().clip(roi);

var image1 = collection.filterDate('2021-01-01','2021-12-31')
                        .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                        .map(maskS2clouds);

var image2 = image1;
var image3 = image1;
var image4 = image1;
var image5 = image1;
var image6 = image1;
var image7 = image1;
var image8 = image1;

//这里可以只设定一个函数，然后分别在  return image.addBands(ndwi_c);将括号中的值进行改变即可

//NDWI
var addNDWI_c = function(image) {
  var ndwi_c = image_c.normalizedDifference(['B3', 'B8']).rename('NDWI_c');
  return image.addBands(ndwi_c);
};
var addNDWI_x = function(image) {
  var ndwi_x = image_x.normalizedDifference(['B3', 'B8']).rename('NDWI_x');
  return image.addBands(ndwi_x);
};
var addNDWI_q = function(image) {
  var ndwi_q = image_q.normalizedDifference(['B3', 'B8']).rename('NDWI_q');
  return image.addBands(ndwi_q);
};
var addNDWI_d = function(image) {
  var ndwi_d = image_d.normalizedDifference(['B3', 'B8']).rename('NDWI_d');
  return image.addBands(ndwi_d);
};

//NDVI
var addNDVI_c = function(image) {
  var ndvi_c = image_c.normalizedDifference(['B8', 'B4']).rename('NDVI_c');
  return image.addBands(ndvi_c);
};
var addNDVI_x = function(image) {
  var ndvi_x = image_x.normalizedDifference(['B8', 'B4']).rename('NDVI_x');
  return image.addBands(ndvi_x);
};
var addNDVI_q = function(image) {
  var ndvi_q = image_q.normalizedDifference(['B8', 'B4']).rename('NDVI_q');
  return image.addBands(ndvi_q);
};
var addNDVI_d = function(image) {
  var ndvi_d = image_d.normalizedDifference(['B8', 'B4']).rename('NDVI_d');
  return image.addBands(ndvi_d);
};

//全部采用高质量合成各个波段

//NDWI
image1 = image1.map(addNDWI_c).qualityMosaic('NDWI_c')
image2 = image2.map(addNDWI_x).qualityMosaic('NDWI_x')
image3 = image3.map(addNDWI_q).qualityMosaic('NDWI_q')
image4 = image4.map(addNDWI_d).qualityMosaic('NDWI_d')
//NDVI
image5 = image5.map(addNDVI_c).qualityMosaic('NDVI_c')
image6 = image6.map(addNDVI_x).qualityMosaic('NDVI_x')
image7 = image7.map(addNDVI_q).qualityMosaic('NDVI_q')
image8 = image8.map(addNDVI_d).qualityMosaic('NDVI_d')



var ndwi_c = image1.select('NDWI_c')
var ndwi_x = image2.select('NDWI_x')
var ndwi_q = image3.select('NDWI_q')
var ndwi_d = image4.select('NDWI_d')

var ndvi_c = image5.select('NDVI_c')
var ndvi_x = image6.select('NDVI_x')
var ndvi_q = image7.select('NDVI_q')
var ndvi_d = image8.select('NDVI_d')

//单独计算一个没有时间差异的NDBI
var image = collection.filterDate('2021-01-01','2021-12-31')
                        .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT',10))
                        .map(maskS2clouds)
                        .mosaic().clip(roi);

var ndbi = image.normalizedDifference(['B12', 'B8']).rename('NDBI');

//将所有波段累加在一起
var image = image.addBands([ndvi_c, ndvi_x, ndvi_q, ndvi_d, ndwi_c, ndwi_x, ndwi_q, ndwi_d, ndbi])

Map.addLayer(image, {bands: ['B4', 'B3', 'B2'],min:0.0, max: 0.3}, 'true');


//阈值分类
var ndbi = ndbi.lt(-0.1)
var ndwi_c = ndwi_c.lt(0.2).gt(-0.5)
var ndwi_x = ndwi_x.lt(0.15).gt(-0.5)
var ndwi_q = ndwi_q.lt(-0.1)
var ndwi_d = ndwi_d.lt(0).gt(-0.5)
var ndvi_c = ndvi_c.gt(0.5)
var ndvi_x = ndvi_x.gt(0.6)//0.3
var ndvi_q = ndvi_q.gt(0.6)//0.3
var ndvi_d = ndvi_d.lt(0.5).gt(0.3)
var crop = ndwi_d.updateMask(ndbi)
                 .updateMask(ndwi_c).updateMask(ndwi_x).updateMask(ndwi_q).updateMask(ndwi_d)
                 .updateMask(ndvi_c).updateMask(ndvi_x).updateMask(ndvi_q).updateMask(ndvi_d)
Map.addLayer(crop, {min: 0, max: 1,palette: ['white','yellow']}, 'crop');

总之，这个代码在看起来有很多重复的地方，后面改进的空间还有很大，大家可以按照我上面所说的注释进行相应的改进，然后尝试。

经过阈值法求得得最后庄稼生长范围，当然阈值法在使用得时候一般可以按照最大值和最小值确定，也就是使用 

ee.Reducer.max(numInputs)

Creates a reducer that outputs the maximum value of its (first) input. If numInputs is greater than one, also outputs the corresponding values of the additional inputs.

ee.Reducer.min(numInputs)

Creates a reducer that outputs the minimum value of its (first) input. If numInputs is greater than one, also outputs the corresponding values of the additional inputs.

Arguments:

numInputs (Integer, default: 1):

The number of inputs.

Returns: Reducer

或者

ee.Reducer.minMax()

Returns a Reducer that computes the minimum and maximum of its inputs.

No arguments.

Returns: Reducer