使用GEE计算多种时序植被指数/建筑指数/水体指数

星踪迷影

已于 2023-08-24 12:29:08 修改

阅读量2.2k

点赞数 7

文章标签： javascript 云计算

于 2023-02-20 22:25:30 首次发布

本文链接：https://blog.csdn.net/m0_71243797/article/details/129132868

版权

以Landsat5 2009年-2010年数据为源数据，使用 JavaScript 分别计算区域内逐月的 NDVI、EVI、RVI、DVI、ARVI、SAVI、NDBI、NDWI、MNDWI 九种指数。

（注：受限于笔者创作时的水平，该文章流程较为繁琐，更精简的GEE时序遥感数据批量下载流程可参考以下文章：）

GEE联合多源时序遥感数据下载：以Landsat-8和Sentinel-2为例https://blog.csdn.net/m0_71243797/article/details/132417873
实现步骤：

        1、划分卫星影像数据，生成一个列表，分别存储2009-2010年逐月的影像数据。

        2、每一张影像去云、去雪、去除水域，计算上述九种指数，并筛除环境掩膜异常的影像，每月所有处理后的影像取中位数，合成为一张影像。

        3、使用前后两个月的数据平均值，填补去云后的空白。

        4、逐月下载数据。

详细代码：

划分具体的研究区域：中国东北平原、波德平原、美国大平原：

var ThreeRegions = ee.FeatureCollection("users/FengZA/PlainRegion_2022_09_12");
var DongBei = ThreeRegions.filterMetadata('ID', 'equals', 0);
var BoDe = ThreeRegions.filterMetadata('ID', 'equals', 1);
var DaPingYuan = ThreeRegions.filterMetadata('ID', 'equals', 2);

输入起始日期和结束日期，日期格式 yyyy-MM-dd：

var date_start = ee.String('2009-01-01');
var date_end = ee.String('2011-01-01');

自定义时间片段，如3为季度，6为半年，12为1年。本次研究为逐月生成影像，值为1：

var date_interval = ee.Number(1);
var monthList = ee.List.sequence(1, 12, date_interval);

函数一：完成第一步和第二步，生成初步处理后的逐月影像列表：

function ImagesDicForMonths(geometry)
{
  //Landsat5去云、去雪、去除水域
  function maskL5sr(image){
    var DilatedCloud = 1 << 1; //云层边缘
    var Cloud = 1 << 3; //云层
    var CloudShadowBitMask = 1 << 4; //云影
    var SnowBitMask = 1 << 5; //积雪
    var Water = 1 << 7; //水域
    var qa = image.select('QA_PIXEL'); //Landsat5的QA_PIXEL波段自带地物和环境分类信息，但不完全准确，需要后续处理
    var mask = qa.bitwiseAnd(DilatedCloud).eq(0)
          .and(qa.bitwiseAnd(Cloud).eq(0))
          .and(qa.bitwiseAnd(CloudShadowBitMask).eq(0))
          .and(qa.bitwiseAnd(Water).eq(0))
          .and(qa.bitwiseAnd(SnowBitMask).eq(0));
    return image.updateMask(mask);
  }
  //生成影像底图波段,利于后续填补去云后的空白
  var backGr = function backGr(image){
    var imgMask = image.select('QA_PIXEL').mask();
    var bG = ee.Image.constant(1)
                     .updateMask(imgMask);
    return image.addBands(bG.rename('background'));
  };
  //以下为指数计算函数
  //计算归一化指数 通用方程
  var Norm_ID = function Norm_ID(B1, B2, name){
    var norm_id = function(img){
      var temp_id = img.normalizedDifference([B1, B2]);
      return img.addBands(temp_id.rename(name));
    };
    return norm_id;
  };
  //计算比例指数 通用方程
  var Two_Divide = function Two_Divide(B1, B2, name){
    var two_divide = function(img){
      var temp_td = img.expression(
        'b1/b2',
        {
          b1: img.select(B1),
          b2: img.select(B2)
        });
      return img.addBands(temp_td.rename(name));
    };
    return two_divide;
  };
  //EVI计算方程
  var EVI = function EVI(image){
    var EVI_image = image.expression(
      'float(G*(nir - red) / (nir + C1*red - C2*blue + L))',
      {
        'G': 2.5,
        'C1': 6.0,
        'C2': 7.5,
        'L': 1.0,
        'nir': image.select('SR_B4'),
        'red': image.select('SR_B3'),
        'blue': image.select('SR_B1')
      });
    return image.addBands(EVI_image.rename('EVI'));
  };
  //DVI计算方程
  var DVI = function DVI(image){
    var DVI_image = image.expression(
      'float((nir - red)/10000)',
      {
        'nir': image.select('SR_B4'),
        'red': image.select('SR_B3')
      });
    return image.addBands(DVI_image.rename('DVI'));
  };
  //ARVI计算方程
  var ARVI = function ARVI(image){
    var ARVI_image = image.expression(
      'float((nir - (2*red) + blue) / (nir + (2*red) + blue))',
      {
        'nir': image.select('SR_B4'),
        'red': image.select('SR_B3'),
        'blue': image.select('SR_B1')
      });
    return image.addBands(ARVI_image.rename('ARVI'));
  };
  //SAVI计算方程
  var SAVI = function SAVI(image){
    var SAVI_image = image.expression(
      'float((nir - red) / (nir + red + L) * (1 + L))',
      {
        'nir': image.select('SR_B4'),
        'red': image.select('SR_B3'),
        'L': 0.5
      });
    return image.addBands(SAVI_image.rename('SAVI'));
  };
  //Landsat5地表反射影像提取,完成区域筛选,去云
  var Landsat_5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
           .merge(ee.ImageCollection("LANDSAT/LT05/C02/T2_L2"))
           .filterBounds(geometry)
           .map(maskL5sr);
  //划分卫星影像数据，生成一个列表，分别存储2009-2010年逐月的影像数据
  var startYear = ee.Number.parse(date_start.slice(0, 4));
  var endYear = ee.Number.parse(date_end.slice(0, 4)).subtract(1);
  var L5_imageList = 
  monthList.map(function(startMonth){
    startMonth = ee.Number(startMonth);
    var endMonth = startMonth.add(date_interval).subtract(1);
    var images = Landsat_5
    .filter(ee.Filter.calendarRange(startYear, endYear, "year"))
    .filter(ee.Filter.calendarRange(startMonth ,endMonth, "month"));
    return images.filterBounds(geometry);
    });
  //计算上述九种指数，每月所有计算后的影像取中位数，合成为一张影像
  var range = ee.List.sequence(0, null, 1, 12);
  L5_imageList = range.map(function(number){
    var image = ee.ImageCollection(L5_imageList.get(ee.Number(number)))
      .map(backGr)
      .map(Norm_ID('SR_B4','SR_B3','NDVI'))
      .map(Norm_ID('SR_B2','SR_B5','MNDWI'))
      .map(Norm_ID('SR_B2','SR_B4','NDWI'))
      .map(Norm_ID('SR_B5','SR_B4','NDBI'))
      .map(Two_Divide('SR_B4','SR_B3','RVI'))
      .map(SAVI)
      .map(EVI)
      .map(DVI)
      .map(ARVI)
      .select(['background','NDVI','EVI','RVI','DVI','SAVI','ARVI','NDBI','NDWI','MNDWI'])
      //在查看影像时，发现有一类影像由于QA_PIXEL波段的问题，出现积雪和云层未正常掩膜去除的现象，需要处理
      //处理方法：每一张影像添加一个新属性：区域内所有像元NDVI的平均值。平均值低于0.05的为异常影像，使用过滤器去除
      .map(function(image){
        var mean = image.select('NDVI').reduceRegion({
          reducer: ee.Reducer.mean(),
          scale: 90,
          });
        return image.set("mean", mean.get('NDVI'));
      })
      .filter(ee.Filter.gt("mean", 0.05))
      .median()
      .clip(geometry);
    return image;
  });
  //添加key字段：月份，生成逐月影像字典
  monthList = monthList.map(function(month){
    var Date = ee.Date.fromYMD(1, month, 1).format("MM");
    return Date;
  });
  var VI_ImageDict = ee.Dictionary.fromLists(monthList, L5_imageList);
  return VI_ImageDict;
}

函数二：完成第三步和第四步：

function interpolation(dictionary, geometry)
{
  //自定义加法和减法函数，使12+1=1， 1-1=12，便于时间插值填补去云后的空白
  var add = function add(date)
  {
    var number = ee.Number.parse(date);
    var add_number = ee.Algorithms.If(
      number.add(date_interval).gt(12),
      ee.Number(1),
      number.add(date_interval)
    );
    var add_date = ee.Date.fromYMD(1, add_number, 1).format("MM");
    return add_date;
  };
  var sub = function sub(date)
  {
    var number = ee.Number.parse(date);
    var sub_number = ee.Algorithms.If(
      number.eq(1),
      ee.Number(12),
      number.subtract(date_interval)
    );
    var sub_date = ee.Date.fromYMD(1, sub_number, 1).format("MM");
    return sub_date;
  };
  //修改月份列表日期为标准格式
  var DateList = ee.List.sequence(1, 12, date_interval)
                   .map(function(date){
                     return ee.Date.fromYMD(1, date, 1).format("MM");
                   }).getInfo();
  var backGround = ee.Image.constant(1).clip(geometry);
  var imageAPList = DateList.map(function(mem){
    var date_1 = add(mem);
    var date_2 = sub(mem);
    var image = ee.Image(elements.get(mem));
    var image_VI = image.select(['NDVI','EVI','RVI','DVI','SAVI','ARVI','NDBI','NDWI','MNDWI']);
    var imgMask = image.select('background').mask();
    //生成目标影像的去云空白区域
    var blankArea = backGround.updateMask(imgMask.not().clip(geometry));
    //生成空白区域的填补影像 
    var image_1 = ee.Image(elements.get(date_1))
                    .select(['NDVI','EVI','RVI','DVI','SAVI','ARVI','NDBI','NDWI','MNDWI'])
                    .updateMask(blankArea);
    var image_2 = ee.Image(elements.get(date_2))
                    .select(['NDVI','EVI','RVI','DVI','SAVI','ARVI','NDBI','NDWI','MNDWI'])
                    .updateMask(blankArea);
    //影像整合，完成空白区域补充 
    var imageSrList = ee.List([image_1, image_2]).flatten();
    var imageAP = ee.ImageCollection.fromImages(imageSrList.add(image_VI)).median();
    //平滑滤波，用半径为3的矩形范围内像元平均值来代替目标像元
    imageAP = imageAP.reduceNeighborhood(ee.Reducer.median(),ee.Kernel.square(3)).rename(['NDVI']);
    //使用循环逐月下载数据
    Export.image.toDrive({
      image: imageAP,
      folder: 'VI_L5_DongBei',
      fileNamePrefix: 'BWM_L5_DongBei_'+ mem,
      description: 'BWM_L5_DongBei_'+ mem,
      region: geometry.geometry().bounds(),
      scale: 90, //重采样使空间分辨率为90m
      crs: "EPSG:4326", //使用WGS84坐标系
      maxPixels: 1e13
    });
    return imageAP;
  });
  
  return ee.Dictionary.fromLists(DateList, imageAPList);
}

最后，调用上述两个函数生成影像，以中国东北平原为例，可修改参数研究区域为波德平原或美国大平原等其他研究区域：

var elements = ImagesDicForMonths(DongBei);
print("elements:", elements);
var VI_2009_2010 = interpolation(elements, DongBei);
print("VI_2009_2010:", VI_2009_2010);

在GEE地图上在线查看数据：

//NDVI显示方案
vis = {
    "min":0,
    "max":1,
    "palette":[
        "FFFFFF","CE7E45","DF923D","F1B555",
        "FCD163","99B718","74A901","66A000",
        "529400","3E8601","207401","056201",
        "004C00","023B01","012E01","011D01",
        "011301"
    ]
};

Map.addLayer(ee.Image(VI_2009_2010.get('01')).select('NDVI'), vis, "1 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('02')).select('NDVI'), vis, "2 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('03')).select('NDVI'), vis, "3 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('04')).select('NDVI'), vis, "4 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('05')).select('NDVI'), vis, "5 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('06')).select('NDVI'), vis, "6 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('07')).select('NDVI'), vis, "7 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('08')).select('NDVI'), vis, "8 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('09')).select('NDVI'), vis, "9 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('10')).select('NDVI'), vis, "10 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('11')).select('NDVI'), vis, "11 NDVI", false);
Map.addLayer(ee.Image(VI_2009_2010.get('12')).select('NDVI'), vis, "12 NDVI", false);