基于 GEE 可视化分析植被覆盖度 FVC 的时序变化

目录

1 植被覆盖度（FVC）

2 完整代码

3 运行结果

---

1 植被覆盖度（FVC）

植被覆盖度指某一地域植物垂直投影面积与该地域面积之比，用百分数表示。森林覆盖率,亦称森林覆被率指一个国家或地区森林面积占土地面积的百分比，是反映一个国家或地区森林面积占有情况或森林资源丰富程度及实现绿化程度的指标，又是确定森林经营和开发利用方针的重要依据之一。

长序列的时间分析，一个卫星系列的数据可能不够，就需要将多个数据整合，提取出需要的波段。目前已经发展了很多利用遥感测量植被覆盖度的方法，较为实用的方法是利用植被指数近似估算植被覆盖度，常用的植被指数为NDVI。下面是像元二分模型的模型：

其中, NDVIsoil 为完全是裸土或无植被覆盖区域的NDVI值，NDVIveg 则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。两个值的计算公式为：

利用这个模型计算植被覆盖度的关键是计算NDVIsoil和NDVIveg。这里有两种假设：

1） 当区域内可以近似取VFCmax=95%，VFCmin=5%时，公式（1）可变为：

NDVImax 和NDVImin分别为区域内最大和最小的NDVI值。由于不可避免存在噪声，NDVImax 和NDVImin一般取一定置信度范围内的最大值与最小值，置信度的取值主要根据图像实际情况来定。

2） 当区域内不能近似取VFCmax=100%，VFCmin=0%时，在有实测数据的情况下，取实测数据中的植被覆盖度的最大值和最小值作为VFCmax和 VFCmin，这两个实测数据对应图像的NDVI作为NDVImax 和NDVImin。在没有实测数据的情况下，取一定置信度范围内的NDVImax 和NDVImin。VFCmax和 VFCmin根据经验估算。

2 完整代码

// 将地图中心定位到指定几何对象，并设置缩放级别为 7
Map.centerObject(table, 7);

// 去云函数
function maskcloud(image) {
    // 定义云阴影位掩码
    var cloudShadowBitMask = 1 << 3;
    // 定义云位掩码
    var cloudsBitMask = 1 << 5;
    // 选择影像的 QA_PIXEL 波段
    var qa = image.select('QA_PIXEL');
    // 创建掩码，去除云阴影和云覆盖区域
    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
    // 更新影像的掩码，并复制系统时间和云覆盖度属性
    return image.updateMask(mask).copyProperties(image, ["system:time_start", 'CLOUD_COVER']);
}

// 计算植被覆盖度函数 
function calFVC(BestVI, region, scale) {
    // 计算指定区域内指定波段的第 5 和第 95 百分位数
    var num = BestVI.reduceRegion({
        reducer: ee.Reducer.percentile([5, 95]),
        geometry: region,
        scale: scale,
        maxPixels: 1e13
    });
    // 获取第 5 百分位数
    var min = ee.Number(num.get("NDVI_p5"));
    // 获取第 95 百分位数
    var max = ee.Number(num.get("NDVI_p95"));

    // 确定百分位区域  
    var greaterPart = BestVI.gt(max);
    var lessPart = BestVI.lt(min);
    var middlePart = ee.Image(1).subtract(greaterPart).subtract(lessPart);

    // 大于百分位 95 的是 1
    var tempf1 = BestVI.subtract(min).divide(max.subtract(min));
    var FVC = ee.Image(1).multiply(greaterPart)
      .add(ee.Image(0).multiply(lessPart))
      .add(tempf1.multiply(middlePart));
    return FVC.rename('FVC');
}

// 引入植被覆盖度函数，导图  
var year_start = 2017;
var year_end = 2022;
var LC8_BANDS = ['SR_B4', 'SR_B5']; // Landsat 8

function FVC_year(year, roi) {
    // 定义起始日期
    var Date_start = ee.Date.fromYMD(year, 6, 1);
    // 定义结束日期
    var Date_end = ee.Date.fromYMD(year, 9, 30);
    // 定义影像集合
    var l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
      .filterDate(Date_start, Date_end)
      .filterBounds(roi)
      .filter(ee.Filter.lte('CLOUD_COVER', 10))
      .select(LC8_BANDS);

    // 计算中位合成影像
    var LandsatMean = l8.median().clip(roi);
    // 释放 l8 变量所占用的内存空间，提高脚本的性能和内存管理效率
    l8 = null;
    // 计算 NDVI
    var NDVI = LandsatMean.normalizedDifference(["SR_B5", "SR_B4"]).rename("NDVI");
    // 计算 FVC
    var FVC = calFVC(NDVI, roi, 60);
    // 创建直方图图表
    var chart = ui.Chart.image.histogram({
        image: FVC,
        region: table,
        scale: 60,
    });
    // 打印图表
    // print(chart);
    // 导出影像到 Google Drive
    Export.image.toDrive({
        image: FVC,
        description: "fvc" + year,
        scale: 60,
        region: table,
        maxPixels: 1e13
    });
    return FVC.rename("FVC");
}

// 定义计算覆盖度的函数
function calculateCoverage(image) {
    // 计算覆盖度
    var coverage = image.reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: table,
        scale: 30,
        maxPixels: 1e13
    }).get('FVC');
    return coverage;
}

// 循环时间年份，执行函数
var coverageData = [];
for (var year = year_start; year < year_end; year++) {
    var ImgFVC = FVC_year(year, table);
    var coverage = calculateCoverage(ImgFVC); // 计算覆盖度均值
    coverageData.push([year, coverage]);// 添加每年平均覆盖度
    var FVC_rank = ImgFVC.where(ImgFVC.lt(0.2), 1)
      .where((ImgFVC.gte(0.2).and(ImgFVC.lt(0.4))), 2)
      .where((ImgFVC.gte(0.4).and(ImgFVC.lt(0.6))), 3)
      .where((ImgFVC.gte(0.6).and(ImgFVC.lt(0.8))), 4)
      .where(ImgFVC.gt(0.8), 5);
    // 添加 FVC 等级图层到地图
    Map.addLayer(FVC_rank, { min: 1, max: 5, palette: ["DCDCDC", "FFEBCD", "99B718", "529400", "011301"] }, "FVC_rank" + year, false);
}

var years = [];
var coverages = [];
for (var i = 0; i < coverageData.length; i++) {
    years.push(coverageData[i][0]);
    coverages.push(coverageData[i][1]);
}
// 打印每年的覆盖度数据
print(coverageData);

// 将二维数组转换为 Earth Engine 图表数据格式
var chartData01 = ee.Array(ee.List(coverageData));

// 以下两行定义的 x 重复，可删除之前的定义
// var y = ee.Array([0.4, 0.5, 0.6, 0.7, 0.8]);
// var x = ee.Array([2017, 2018, 2019, 2020, 2021]);
var x = chartData01.slice(1, 0, 1, 1); // 数组切片函数 
var y = chartData01.slice(1, 1);

// 创建 Earth Engine 图表
var coverageChart = ui.Chart.array.values({
    array: y,
    axis: 0,
    xLabels: x
})
  .setChartType('LineChart')
  .setOptions({
        title: 'Yearly Average Coverage',
        hAxis: { title: 'Year' },
        vAxis: { title: 'Coverage' },
        lineWidth: 1,
        pointSize: 3,
        series: {
            0: { color: 'blue' }
        }
    });
// 显示图表
print(coverageChart);

3 运行结果

研究区植被覆盖度可视化结果

点击RUN即可下载数据

研究区植被覆盖度时序变化折线图