基于Landsat 1-7提取点或区域的覆盖度

主要任务:

植被覆盖度(fractional vegetation coverage)是一些生态水文过程模型的重要输入参数,本文的主要目的是提取年均覆盖度的时间序列,以约束生态水文模型。

计算原理:

基于像元二分模型(Gao et al., 2020, Redirecting),可以得到覆盖度的计算:

FVC=\frac{NDVI-NDVI_{\infty }}{NDVI_{soil}-NDVI}

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

使用的数据集:

产品名字        GEE数据产品

时间范围

(实际用的图像的范围)

计算覆盖度的波段
Landsat 1LANDSAT/LM01/C02/T21972-1978 (1974-1974)'B6', 'B5'
Landsat 2LANDSAT/LM02/C02/T21975-1982(1975-1977)'B6', 'B5'
Landsat 3LANDSAT/LM03/C02/T21978-1983(1978)'B6', 'B5'
Landsat 4LANDSAT/LT04/C02/T1_L2

1982-1993

(-)

'SR_B4', 'SR_B3'
Landsat 5LANDSAT/LT05/C02/T1_L21984-2012(1986-2010)'SR_B4', 'SR_B3'
Landsat 7LANDSAT/LE07/C02/T1_L21999-2021(2011-2012)'SR_B4', 'SR_B3'
Landsat 8LANDSAT/LC08/C02/T1_L21975-1982(2013-2023)'SR_B5', 'SR_B4'

成果:

代码

代码修改的话,只需要修改第一行的point坐标。

// 定义要查找的经纬度点
var point = ee.Geometry.Point([108.8955, 37.6935]); // Mu Us
//var point = ee.Geometry.Point([116.6821648, 43.55310465]); // Xilingol
// 定义边界
var bounds = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(point).first().geometry();


// Define the time range
var start = '1972-01-01';
var end = '2024-01-01';

// Filter Landsat data
var LC1 = ee.ImageCollection("LANDSAT/LM01/C02/T2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC1);

var LC2 = ee.ImageCollection("LANDSAT/LM02/C02/T2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC2);

var LC3 = ee.ImageCollection("LANDSAT/LM03/C02/T2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC3);

var LC4 = ee.ImageCollection("LANDSAT/LM04/C02/T2")
                          .filterBounds(point)
                          .filterDate(start, end);
//print(LC4);

var LC4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC4);


var LC5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC5);

var LC7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC7);

var LC8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                          .filterBounds(point)
                          .filterDate(start, end);
print(LC8);
// Display the result
Map.centerObject(point, 10);



// NDVI计算函数
function calcNDVI(image, bands) {
  var ndvi = image.normalizedDifference(bands);
  return image.addBands(ndvi.rename('NDVI'));
}

// FVC计算函数
function calFVC(image, minNDVI, maxNDVI) {
  var FVC = image.subtract(minNDVI).divide(maxNDVI.subtract(minNDVI))
                 //.clamp(0, 1); // 保证 FVC 值在 0 到 1 之间
  return FVC.rename('FVC');
}

function rmCloud123(image) {  
  var cloudBitMask = (1 << 3); 
  var qa = image.select("QA_PIXEL");  
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)  
  return image.updateMask(mask);  
} 

function rmCloud4578(image) {  
  var cloudBitMask = (1 << 3); 
  var cloudShadowBitMask = (1 << 4)
  var cloudsBitMask1 = (1 << 5);  
  var qa = image.select("QA_PIXEL");  
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)  
                 .and(qa.bitwiseAnd(cloudBitMask).eq(0))
                 .and(qa.bitwiseAnd(cloudsBitMask1).eq(0));  
  return image.updateMask(mask);  
} 
  
// 获取跨所有传感器的 NDVI 5% 和 95% 分位数
function getAllSensorsNDVIPercentiles(startYear, endYear) {
  var start = ee.Date.fromYMD(startYear, 1, 1);
  var stop = ee.Date.fromYMD(endYear, 12, 31);

  // 所有传感器的集合
  var LM01 = ee.ImageCollection("LANDSAT/LM01/C02/T2")
                          .filterBounds(point)
                          .filterDate('1973-01-01', '1974-12-31')
                          .map(function(image) {
                            return rmCloud123(image);
                          })
                          .map(function(image) {
                            var bands = ['B6', 'B5'];
                            return calcNDVI(image, bands);
                          });                         
  
  var LM02 = ee.ImageCollection("LANDSAT/LM02/C02/T2")
                          .filterBounds(point)
                          .filterDate('1975-01-01', '1977-12-31')
                          .map(function(image) {
                            return rmCloud123(image);
                          })
                          .map(function(image) {
                            var bands = ['B6', 'B5'];
                            return calcNDVI(image, bands);
                          });   
  
  var LM03 = ee.ImageCollection("LANDSAT/LM03/C02/T2")
                          .filterBounds(point)
                          .filterDate(start, stop)
                          .map(function(image) {
                            return rmCloud123(image);
                          })
                          .map(function(image) {
                            var bands = ['B6', 'B5'];
                            return calcNDVI(image, bands);
                          });  
  
  var LT04 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2")
                          .filterBounds(point)
                          .filterDate(start, stop)
                          .map(function(image) {
                            return rmCloud4578(image);
                          })
                          .map(function(image) {
                            var bands = ['SR_B4', 'SR_B3'];
                            return calcNDVI(image, bands);
                          });

  var LT05 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
                          .filterBounds(point)
                          .filterDate(start, stop)
                          .map(function(image) {
                            return rmCloud4578(image);
                          })
                          .map(function(image) {
                            var bands = ['SR_B4', 'SR_B3'];
                            return calcNDVI(image, bands);
                          });
                  
  var LE07 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
                          .filterBounds(point)
                          .filterDate(start, stop)
                          .map(function(image) {
                            return rmCloud4578(image);
                          })
                          .map(function(image) {
                            var bands = ['SR_B4', 'SR_B3'];
                            return calcNDVI(image, bands);
                          });
  
  var LC08 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                          .filterBounds(point)
                          .filterDate(start, stop)
                          .map(function(image) {
                            return rmCloud4578(image);
                          })
                          .map(function(image) {
                            var bands = ['SR_B5', 'SR_B4'];
                            return calcNDVI(image, bands);
                          });
                
  var allSensors = LM01.merge(LM02).merge(LM03).merge(LT04).merge(LT05).merge(LE07).merge(LC08);
  var ndvi = allSensors.select('NDVI');
  var percentiles = ndvi.reduce(ee.Reducer.percentile([2, 98]));
  
  return percentiles;
}

function applyNDVI(image, sensor) {
  var bands;
  // Set bands based on sensor type
  if (sensor === 'LM01' || sensor === 'LM02' || sensor === 'LM03') {
    bands = ['B6', 'B5'];  // Landsat 1, 2, 3
  } else if (sensor === 'LT05' || sensor === 'LT04') {
    bands = ['SR_B4', 'SR_B3'];  // Landsat 4, 5, 7
  } else if (sensor === 'LC08') {
    bands = ['SR_B5', 'SR_B4'];  // Landsat 8
  } else {
    bands = ['SR_B4', 'SR_B3'];  // Default to Landsat 4, 5, 7 bands or handle other cases
  }

  return calcNDVI(image, bands);
}

// 生成时间序列
function getFVCtimeSeries(sensor, startYear, endYear) {
  //var ndviPercentiles = getNDVIPercentiles(2012, 2023);
  var ndviPercentiles = getAllSensorsNDVIPercentiles(1972, 2019);
  
  var minNDVI = ndviPercentiles.select('NDVI_p2');
  var maxNDVI = ndviPercentiles.select('NDVI_p98');
  //var minNDVI = ee.Image.constant(0.01).rename('minNDVI');
  //var maxNDVI = ee.Image.constant(0.3).rename('maxNDVI');

  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 collection = ee.ImageCollection(ee.Algorithms.If(
      sensor === 'LM01' || sensor === 'LM02' || sensor === 'LM03' , 
      ee.ImageCollection('LANDSAT/' + sensor + '/C02/T2')
                      .filterBounds(point)
                      .filterDate(start, stop)
                      .map(function(img) {
                        var img = rmCloud123(img);
                        var img = applyNDVI(img, sensor)
                        return img
                      }),
      ee.ImageCollection('LANDSAT/' + sensor + '/C02/T1_L2')
                .filterBounds(point)
                .filterDate(start, stop)
                .map(function(img) {
                  var img = rmCloud4578(img);
                  var img = applyNDVI(img, sensor)
                  return img
                })
                ));

    var ndvi = collection.select('NDVI').mean().clip(bounds); // 使用均值 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));
  }));
}

// 获取 FVC 时间序列
var fvcAndNDVITimeSeries = ee.FeatureCollection([
  //getFVCtimeSeries('LM01', 1972, 1978),
  //getFVCtimeSeries('LM02', 1975, 1982),
  //getFVCtimeSeries('LM03', 1978, 1983),
  getFVCtimeSeries('LM01', 1974, 1974),
  getFVCtimeSeries('LM02', 1975, 1977),
  getFVCtimeSeries('LM03', 1978, 1978),
  //getFVCtimeSeries('LT04', 1982, 1993),
  getFVCtimeSeries('LT05', 1986, 2010),
  getFVCtimeSeries('LE07', 2011, 2012),
  getFVCtimeSeries('LC08', 2013, 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);

// 绘制每年 FVC 图像到地图
function drawFVCImages(sensor, startYear, endYear) {
  //var ndviPercentiles = getNDVIPercentiles(2012, 2023);
  var ndviPercentiles = getAllSensorsNDVIPercentiles(1972, 2019);
  
  var minNDVI = ndviPercentiles.select('NDVI_p2');
  var maxNDVI = ndviPercentiles.select('NDVI_p98');
  
  var years = ee.List.sequence(startYear, endYear);
  years.evaluate(function(yearsList) {
    yearsList.forEach(function(year) {
      var start = ee.Date.fromYMD(year, 1, 1);
      var stop = ee.Date.fromYMD(year, 12, 31);
      
      var collection = ee.ImageCollection('LANDSAT/' + sensor + '/C02/T1_L2')
                        .filterBounds(point)
                        .filterDate(start, stop)
                        .map(function(img) {
                          return rmCloud4578(img, sensor);
                        })
                        .map(function(image) {
                          var bands = sensor === 'LC08' ? ['SR_B5', 'SR_B4'] : ['SR_B4', 'SR_B3'];
                          return calcNDVI(image, bands);
                        });
      
      var ndvi = collection.qualityMosaic('NDVI').clip(bounds);
      var ndvimask = ndvi.select('NDVI').updateMask(ndvi.select('NDVI').gt(0));
      var FVC = calFVC(ndvimask, minNDVI, maxNDVI);
      
      // 定义颜色映射
      var visParams = {
        min: 0,
        max: 0.8,
        palette: [
          'blue',    // 0-0.1
          'cyan',    // 0.1-0.2
          'green',   // 0.2-0.3
          'yellow',  // 0.3-0.4
          'orange',  // 0.4-0.5
          'red',     // 0.5-0.6
          'darkred'  // 0.6-0.8
        ]};
        
      Map.addLayer(FVC, visParams, year.toString() + ' FVC');
    });
  });
}

// 绘制每年的 FVC 图像
//drawFVCImages('LT05', 1986, 2010);
drawFVCImages('LE07', 2011, 2012);
//drawFVCImages('LC08', 2013, 2023);

// 设置地图视图到点的位置
Map.centerObject(point, 10);
Map.addLayer(ee.Image().paint(point.buffer(1000).bounds(), 1, 2).visualize({palette: ['blue']}), {}, '正方形标记');

  • 21
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值