GEE: KNDVI代码

参考内容

Github
Group群组
核植被指数（kNDVI）

NDVI/EVI优缺点

目前，NDVI依然是比较流行的植被指数。但由于NDVI自身是“非线性”的，当植被覆盖持续增加时，比值NIR/RED持续增加，这种非线性的转换过程也使得NDVI对RED波段过度敏感，因此，在植被覆盖较高的或者生物含量较高地区，NDVI不可避免的伴有不同程度的饱和现象。此外，NDVI对红光（叶绿素）高度敏感，当叶绿素浓度超过一定浓度时，红光对叶绿素不再具有敏感性，很快出现饱和现象。为了调整NDVI以适应大气和土壤噪声，特别是在植被茂密地区，减缓NDVI饱和性问题，在NDVI中加入了蓝光波段，提出了增强型植被指数EVI，但是EVI中饱和现象依旧存在。

kNDVI

2021年Gustau在SCIENCE ADVANCES发表了【A unified vegetation index for quantifying the terrestrial biosphere】一文，提出了核NDVI，该研究采用机器学习的核方法理论，使NDVI线性化，避免了植被指数非线性及饱和性问题。

代码

GitHub中有详细代码，其中GEE代码进入后能够在Reader部分看到 “users/ispguv/kNDVI”，共有三部分代码。

其中一个代码如下：

var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

// -------------------------------------------------------------------------------------- 
// 1- MODIS cloud masking.
// -------------------------------------------------------------------------------------- 

/*
 * A function that returns an image containing just the specified QA bits.
 *
 * Args:
 *   image - The QA Image to get bits from.
 *   start - The first bit position, 0-based.
 *   end   - The last bit position, inclusive.
 *   name  - A name for the output image.
 */
var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out pixels that did not have observations.
var maskEmptyPixels = function(image) {
  // Find pixels that had observations.
  var withObs = image.select('num_observations_1km').gt(0);
  return image.updateMask(withObs);
};

// mask mod15 function
var maskmod15 = function(image1) {
  var maskbad = ((image1.select('FparLai_QC')).eq(0))
  .or((image1.select('FparLai_QC')).eq(2))
  .or((image1.select('FparLai_QC')).eq(32))
  .or((image1.select('FparLai_QC')).eq(34));
  return image1.mask(maskbad);

};

// A function to mask out cloudy pixels.
var maskClouds = function(image) {
  // Select the QA band.
  var QA = image.select('state_1km');
  // Get the internal_cloud_algorithm_flag bit.
  // var internalCloud = getQABits(QA, 10, 10, 'internal_cloud_algorithm_flag');
  // 10 is cloud shadow!
  // Get the cloud_state bits and find cloudy areas.
  var internalCloud = getQABits(QA, 0, 1, 'cloud_state').expression("b(0) == 1 || b(0) == 2");
  image = image.addBands(internalCloud);
  return image.updateMask(internalCloud.eq(0));
};

// A function to mask out no Land pixels.
var masknoLand = function(image) {
  // Select the QA band.
  var QA = image.select('state_1km');
  // Get the internal_cloud_algorithm_flag bit.
  // var internalCloud = getQABits(QA, 10, 10, 'internal_cloud_algorithm_flag');
  // 10 is cloud shadow!
  // Get the cloud_state bits and find cloudy areas.
  var internalCloud = getQABits(QA, 3, 5, 'Land/water flag').expression("b(0) > 1 || b(0) == 0 ");
  image = image.addBands(internalCloud);
  return image.updateMask(internalCloud.eq(0));
};

// -------------------------------------------------------------------------------------- 
// 2- Select your image collection
// -------------------------------------------------------------------------------------- 
var collection = ee.ImageCollection('MODIS/006/MOD09GA')
                   .filterDate('2018-01-01', '2019-01-11'); 

var MODIS_LAI = ee.ImageCollection("MODIS/006/MCD15A3H")
  .filterDate('2018-01-01', '2019-01-11').map(maskmod15).select('Lai');

// print(MODIS_LAI);

// Clip the geometry using a map function
function clipImage(image) {
  return image.clip(Albufera);
}
var collection = collection.map(clipImage);

// Mask out no land areas and areas that were not observed.
collection = collection.map(maskEmptyPixels).map(masknoLand);

// Map the cloud masking function over the collection.
var collectionCloudMasked = collection.map(maskClouds);

// -------------------------------------------------------------------------------------- 
// 3- Set MODIS in reflectance values 
//    This is important to have some intuition about the sigma parameter in the RBF kernel
// -------------------------------------------------------------------------------------- 
collectionCloudMasked = collectionCloudMasked.select(
  ['sur_refl_b01', 'sur_refl_b02', 'sur_refl_b03', 'sur_refl_b04'])
  .map(function(image) {
    //return image.multiply(ee.Image.constant(1e-4))
    //            .copyProperties(image, image.propertyNames());
    return image.multiply(ee.Image.constant(1e-4))
                .set({'system:time_start': image.get('system:time_start')});
  });

var vizParams = {
  bands: ["sur_refl_b01","sur_refl_b04","sur_refl_b03"],
  min: 0, max: 0.15
};
Map.centerObject(Albufera, 9);
//Map.addLayer(collectionCloudMasked, vizParams, 'Image');

// -------------------------------------------------------------------------------------- 
// NDVI
// -------------------------------------------------------------------------------------- 

var addNDVI = function(image) {
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  var ndvi = nir.subtract(red).divide(nir.add(red)).select([0],['ndvi']);
  return image.addBands(ndvi);
};

// -------------------------------------------------------------------------------------- 
// NIRv
// -------------------------------------------------------------------------------------- 

var addNIRV = function(image) {
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  var nirv = ((nir.subtract(red).divide(nir.add(red))).subtract(0.08))
             .multiply(nir).select([0],['nirv']);
  return image.addBands(nirv);
};

// -------------------------------------------------------------------------------------- 
// kNDVI
// -------------------------------------------------------------------------------------- 

var addKNDVI = function(image) {
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  var D2 = nir.subtract(red).pow(2).select([0],['d2']); // Compute D2 and rename it to d2

 // -------------------------------------------------------------------------------------- 
  // Fix sigma parameter to a reasonable value for all pixels in the region or for your specific problem 
  // We recommend values ranging between [0.1-0.3] when working with (normalized) reflectance data.
  var sigma = 0.15;

  // kDNVI
  var kndvi = D2.divide(sigma).divide(sigma).divide(4.0).tanh().select([0], ['kndvi']);
  
  return image.addBands(kndvi);
};

// -------------------------------------------------------------------------------------- 
// Plot results
// -------------------------------------------------------------------------------------- 

var collectionNDVI = collectionCloudMasked.map(addNDVI);
collectionNDVI = collectionNDVI.map(addNIRV);
collectionNDVI = collectionNDVI.map(addKNDVI);

var mean_ndvi = collectionNDVI.select(['ndvi']).mean();
var mean_nirv = collectionNDVI.select(['nirv']).mean();
var mean_kndvi = collectionNDVI.select(['kndvi']).mean();

Map.addLayer(mean_ndvi, {'min': 0, 'max': 1, 'palette': palette}, 'NDVI');
Map.addLayer(mean_nirv, {'min': 0, 'max': 1, 'palette': palette}, 'NIRv');
Map.addLayer(mean_kndvi, {'min': 0, 'max': 1, 'palette': palette}, 'kNVDI');

// Create an image time series chart.
var chart = ui.Chart.image.series({
  imageCollection: collectionNDVI.select(['ndvi','nirv','kndvi']),
  region: Location,
  reducer: ee.Reducer.first(),
  scale: 500
}).setOptions({ title: '', hAxis: {title: 'Time'}});

// Add the chart to the map.
chart.style().set({
  position: 'bottom-right', width: '400px', height: '300px'
});
Map.add(chart);

// Create an image time series chart.
var chart2 = ui.Chart.image.series({
  imageCollection: MODIS_LAI.select('Lai'),
  region: Location,
  reducer: ee.Reducer.first(),
  scale: 500
}).setOptions({
  title: '',
  vAxis: {title: 'LAI * 10'}, hAxis: {title: 'Time'}
          
});

// Add the chart to the map.
chart2.style().set({
  position: 'bottom-left',
  width: '400px',
  height: '300px'
});
Map.add(chart2);