基于MODIS提取某个点的覆盖度

主要任务：

植被覆盖度（fractional vegetation coverage）是一些生态水文过程模型的重要输入参数。前一篇博客介绍了基于Landsat提取覆盖度的方法，本文是为了对覆盖度的结果进行多个产品之间的验证。

计算原理：

基于像元二分模型（Gao et al., 2020, Redirecting），可以得到覆盖度的计算：

​

其中NDVI是年平均NDVI，根据下表的波段计算的（前一个为nir波段，后一个为red波段），NDVIsoil和NDVI∞分别为裸土时的NDVI和植被饱和时的NDVI，有多种方法确定，这里采取NDVI的所有值的2%，98%分位数对应的NDVI计算，也有用定值计算的，如分别取0.05，0.9（可以参考Gao这篇文献。但是实操下来，还是用百分位数得到的覆盖度更符合我之前对这个地区调查的真实覆盖度。

使用的数据集：

MODIS/006/MOD09A1， 时间分辨率8day，计算NDVI的波段， 近红外nir是sur_refl_b02, 红光波段red是sur_refl_b01。

MODIS结果

​

​

两个产品结果的对比

Landsat的NDVI比MODIS的NDVI低很多，这很可能是因为MODIS分辨率是1000m更低，包括其他高NDVI的地物。

但是Landsat和MODIS的覆盖度FVC很接近，尤其是可以看出1980-2000年毛乌素沙地的平均覆盖度在0.2~0.25之间，可以较好的提供对模型参数（覆盖度）的约束。

​

​

代码

// 定义兴趣区域（用点或边界）
var point = ee.Geometry.Point([108.8955, 37.6935]); // 示例经纬度
var buffer = point.buffer(10000); // 10公里缓冲区

// 定义时间范围
var startDate = '2000-01-01';
var endDate = '2023-12-31';


// 选择 MODIS 反射率数据集（例如 MOD09A1）
var modisCollection = ee.ImageCollection('MODIS/006/MOD09A1')
  .filterBounds(buffer)
  .filterDate(startDate, endDate)
  .map(function(image){
    var ndvi = image.normalizedDifference(['sur_refl_b02', 'sur_refl_b01'])
    return image.addBands(ndvi.rename('NDVI')).clip(buffer)
  });

// 计算NDVI的分位数
var ndviCollection = modisCollection.select('NDVI');
print(ndviCollection, 'NDVI Collection');
var ndviPercentiles = ndviCollection.reduce(ee.Reducer.percentile([2, 98]));
print(ndviPercentiles, 'NDVI Percentiles');


// 获取 NDVI 分位数
var minNDVI = ndviPercentiles.select('NDVI_p2');
var maxNDVI = ndviPercentiles.select('NDVI_p98');
function calFVC(image, minNDVI, maxNDVI) {
  var FVC = image.subtract(minNDVI).divide(maxNDVI.subtract(minNDVI))
                 //.clamp(0, 1); // 保证 FVC 值在 0 到 1 之间
  return FVC.rename('FVC');
}

// 计算每年 FVC 的平均值
// 生成时间序列
function getFVCtimeSeries(startYear, endYear) {
  var years = ee.List.sequence(startYear, endYear);
  return ee.FeatureCollection(years.map(function(year) {
    var start = ee.Date.fromYMD(year, 1, 1);
    var stop = ee.Date.fromYMD(year, 12, 31);
    var modisCollection = ee.ImageCollection('MODIS/006/MOD09A1')
        .filterBounds(buffer)
        .filterDate(start, stop)
        .map(function(image){
          var ndvi = image.normalizedDifference(['sur_refl_b02', 'sur_refl_b01'])
          return image.addBands(ndvi.rename('NDVI')).clip(buffer)
        });

    var ndvi = modisCollection.select('NDVI').mean().clip(buffer); // 使用均值 NDVI
    var ndvimask = ndvi.updateMask(ndvi.gt(0));
    var FVC = calFVC(ndvimask, minNDVI, maxNDVI);
    
    var ndviValue = ndvi.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: point,
      scale: 100,
      maxPixels: 1e13
    }).set('year', year);
    
    var FVCvalue = FVC.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: point,
      scale: 100,
      maxPixels: 1e13
    }).set('year', year);
    
    return ee.Feature(null, ee.Dictionary(ndviValue).combine(FVCvalue, true));
  }));
}

var fvcAndNDVITimeSeries = ee.FeatureCollection([getFVCtimeSeries(2000, 2023),]).flatten();

// 打印时间序列数据
print('FVC 和 NDVI 时间序列:', fvcAndNDVITimeSeries);

// 可视化 FVC 和 NDVI 时间序列
var fvcChart = ui.Chart.feature.byFeature(fvcAndNDVITimeSeries, 'year', 'FVC')
  .setOptions({
    title: 'FVC 时间序列',
    hAxis: {title: '年份'},
    vAxis: {
      title: 'FVC 值',
      minValue: 0,
      maxValue: 0.8
    },
    lineWidth: 1,
    pointSize: 3
  });

var ndviChart = ui.Chart.feature.byFeature(fvcAndNDVITimeSeries, 'year', 'NDVI')
  .setOptions({
    title: 'NDVI 时间序列',
    hAxis: {title: '年份'},
    vAxis: {
      title: 'NDVI 值',
      minValue: -0.1,
      maxValue: 0.5
    },
    lineWidth: 1,
    pointSize: 3
  });

print(fvcChart);
print(ndviChart);