在GEE中如何优雅的计算植被指数

在GEE中，基于Landsat8影像，实现计算任意多个植被指数的函数，并且能够方便的添加或更改植被指数，同时这个函数可以在不同影像上复用。

我才是刘双喜

475人浏览 · 2024-12-27 21:13:30

我才是刘双喜 · 2024-12-27 21:13:30 发布

摘要：在GEE中，基于Landsat8影像，实现计算任意多个植被指数的函数，并且能够方便的添加或更改植被指数，同时这个函数可以在不同影像上复用。

GEE作为一个强大的云计算平台，提供了丰富的遥感数据和强大的计算能力，使得计算植被指数变得更加高效和便捷。大多数教程中会为根据计算公式为每一个植被指数实现一个函数，然后通过map调用函数。

//NDVI 
function NDVI(image) { 
  return image.addBands( 
    image.normalizedDifference(["B5", "B4"]) 
         .rename("NDVI")); 
} 
//NDWI 
function NDWI(image) { 
  return image.addBands( 
    image.normalizedDifference(["B3", "B5"]) 
         .rename("NDWI")); 
} 
//AWEI 
function AWEI(image) { 
  return image.addBands( 
    image.expression('4*(green-SWIR1)-(0.25*NIR+2.75*SWIR2)',
            {
              green:image.select('B3'),
              NIR:image.select('B5'),
              SWIR1:image.select('B6'),
              SWIR2:image.select('B7'),
            })
          .float()
         .rename("AWEI")); 
}
//加载遥感影像 
var l8Col = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR") 
              .filterBounds(roi) 
              .filterDate(dataRage)
              .map(NDVI)
              .map(NDWI)
              .map(AWEI);

我承认这种方式比较容易实现，但是假如，我是说假如哦，你导师交给你一个任务，让你用Modis、Landsat、Sentinel分别计算100个不同的植被指数，你又该如何应对呢？

在软件设计流程中需要考虑代码的可移植性、可复用性、可维护性、可读性、可测试性、健壮性等等。因此，我们需要编写一个用来计算任意多个植被指数的函数，并且能够方便的添加或更改植被指数，同时这个函数可以在不同影像上复用。

为了实现在不同影像上的复用，在加载ImageCollection时应当对影像的波段按照各自的描述进行统一命名。

var imageCol = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                  .select(
                    ["SR_B2", "SR_B3", "SR_B4", "SR_B5", "SR_B6", "SR_B7", "QA_PIXEL"],
                    ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2", "QA_PIXEL"]
                  )
                  .filterBounds(roi)
                  .filterDate('2021-01-01', '2021-12-31')
                  .map(maskclouds)
                  .map(applyScaleFactors);

将各个指数名称和计算公式写在一个字典里，通过维护这个字典来实现添加或删除任意多个以及公式的修改。

var index_dict = {
  'EVI': '2.5 * ((NIR - Red) / (NIR + 6 * Red - 7.5 * Blue + 1))',
  'NDVI': '(NIR - Red) / (NIR + Red)',
  'DVI': 'NIR - Red',
  'RVI': 'NIR / Red',
  'NDWI': '(Green - NIR) / (Green + NIR)',
  'MNDWI': '(Green - SWIR1) / (Green + SWIR1)',
  'NDBI': '(SWIR1 - NIR) / (SWIR1 + NIR)',
  'NDSI': '(Green - SWIR1) / (Green + SWIR1)',
  'LSWI': '(NIR-SWIR1)/(NIR+SWIR1)',
  'AWEI_nsh': '4*(Green-SWIR1)-(0.25*NIR+2.75*SWIR2)',
  'AWEI_sh': 'Blue+2.5*Green-1.5*(NIR+SWIR1)-0.25*SWIR2',
  'NWI': '(Blue-(NIR+SWIR1+SWIR2)) / (Blue+(NIR+SWIR1+SWIR2))',
  'RRI': '(NIR - Red) / (NIR + Red)',
  'BSI': '(Red - Blue) / (Red + Blue)',
  'HUE': '(Green - NIR) / (Green + NIR)'
};

然后将该字典和需要计算植被指数的ImageCollection作为参数传入culIndexs函数，culIndexs函数将遍历这个字典，按照给定的公式计算植被指数，并将波段名称修改为指数名称。最后返回一个包含字典中所有指数的ImageCollection。

var culIndexs = function (imageCol, index_dict) {
    for (var key in index_dict) {
      var map_function = function (image) {
        return image.addBands(
          image.expression(index_dict[key], {
              Blue: image.select("Blue"),
              Green: image.select("Green"),
              Red: image.select("Red"),
              NIR: image.select("NIR"),
              SWIR1: image.select("SWIR1"),
              SWIR2: image.select("SWIR2"),
            }).float().rename(key)
        );
      };
      imageCol = imageCol.map(map_function);
    }
    return imageCol;
};
imageCol = culIndexs(imageCol, index_dict);

完整代码如下：

var roi = ee.Geometry.Polygon(
        [[[98.76874122458284, 27.476046782884104],
          [98.76874122458284, 26.36424630156598],
          [100.16400489645784, 26.36424630156598],
          [100.16400489645784, 27.476046782884104]]], null, false);
Map.centerObject(roi, 8);                 
var maskclouds = function (img) {
    var cloudShadowBitMask = 1 << 4;
    var cloudsBitMask = 1 << 3;
    var qa = img.select("QA_PIXEL");
    var mask = qa
      .bitwiseAnd(cloudShadowBitMask)
      .eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
    return img.addBands(img.updateMask(mask), null, true);
};

var applyScaleFactors = function (image) {
    var opticalBands = image
      .select(["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"])
      .multiply(0.0000275)
      .add(-0.2).float();
    return image.addBands(opticalBands, null, true);
};

var imageCol = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                  .select(
                    ["SR_B2", "SR_B3", "SR_B4", "SR_B5", "SR_B6", "SR_B7", "QA_PIXEL"],
                    ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2", "QA_PIXEL"])
                  .filterBounds(roi)
                  .filterDate('2021-01-01', '2021-12-31')
                  .map(maskclouds)
                  .map(applyScaleFactors);
var index_dict = {
  'EVI': '2.5 * ((NIR - Red) / (NIR + 6 * Red - 7.5 * Blue + 1))',
  'NDVI': '(NIR - Red) / (NIR + Red)',
  'DVI': 'NIR - Red',
  'RVI': 'NIR / Red',
  'NDWI': '(Green - NIR) / (Green + NIR)',
  'MNDWI': '(Green - SWIR1) / (Green + SWIR1)',
  'NDBI': '(SWIR1 - NIR) / (SWIR1 + NIR)',
  'NDSI': '(Green - SWIR1) / (Green + SWIR1)',
  'LSWI': '(NIR-SWIR1)/(NIR+SWIR1)',
  'AWEI_nsh': '4*(Green-SWIR1)-(0.25*NIR+2.75*SWIR2)',
  'AWEI_sh': 'Blue+2.5*Green-1.5*(NIR+SWIR1)-0.25*SWIR2',
  'NWI': '(Blue-(NIR+SWIR1+SWIR2)) / (Blue+(NIR+SWIR1+SWIR2))',
  'RRI': '(NIR - Red) / (NIR + Red)',
  'BSI': '(Red - Blue) / (Red + Blue)',
  'HUE': '(Green - NIR) / (Green + NIR)'
};                  
var culIndexs = function (imageCol, index_dict) {
    for (var key in index_dict) {
      var map_function = function (image) {
        return image.addBands(
          image.expression(index_dict[key], {
              Blue: image.select("Blue"),
              Green: image.select("Green"),
              Red: image.select("Red"),
              NIR: image.select("NIR"),
              SWIR1: image.select("SWIR1"),
              SWIR2: image.select("SWIR2"),
            }).float().rename(key)
        );
      };
      imageCol = imageCol.map(map_function);
    }
    return imageCol;
};
imageCol = culIndexs(imageCol, index_dict);
var image_median = imageCol.median().clip(roi);
var visParam = {  
  bands: ['NDVI'],
  min: -0.2,   
  max: 0.8,  
  palette: ["FFFFFF", "CE7E45", "DF923D", "F1B555", "FCD163",   
            "99B718", "74A901", "66A000", "529400", "3E8601",   
            "207401", "056201", "004C00", "023B01", "012E01",   
            "011D01", "011301"]  
};  
Map.addLayer(image_median, visParam, "image_median");