防风固沙（八）植被因子的计算

嘿嘿灵活的小胖子

已于 2024-06-03 16:53:46 修改

阅读量472

点赞数 22

分类专栏： GEE 防风固沙生态系统服务功能文章标签： arcgis

于 2024-05-26 22:51:07 首次发布

本文链接：https://blog.csdn.net/m0_53209421/article/details/139222546

版权

生态系统服务功能同时被 3 个专栏收录

9 篇文章 0 订阅

订阅专栏

防风固沙

6 篇文章 0 订阅

订阅专栏

GEE

3 篇文章 0 订阅

订阅专栏

大家好，我们今天学习植被因子C的计算，我是通过GEE计算的，下面让我们开始吧！！！

（1）植被因子C的计算公式

$C=e^{aSC}$

式中：当植被类型为林地时，a取-0.1535，草地：a取-0.1151，耕地：a取-0.0438，水域与建设用地：a取0，裸地：a取-0.0768.

（2）GEE计算植被因子C

var start_year = 2020
var end_year = 2020

Map.addLayer(roi, {'color':'blue'}, 'StudyArea');
Map.centerObject(roi, 6);


for(var i = start_year; i <= end_year; i++){
    var start_date = (i) + '-05-01';
    var end_date   = (i + 1) + '-10-01';
    
    var dataset = ee.ImageCollection("MODIS/061/MOD13A1")
                    .filterBounds(roi)
                    .filterDate(start_date, end_date)
                    .select('NDVI');
                    
    var bestEVI = dataset.mosaic().clip(roi).multiply(0.0001);

    var EVIVis = {
      min: -0.2,
      max: 0.5,
      palette: ['bbe029', '0a9501', '074b03'],
    };
    Map.addLayer(bestEVI, EVIVis, 'EVI'+i);

    //calculate the 5% and 95% num of the image
    var num =bestEVI.reduceRegion({
          reducer:ee.Reducer.percentile([5,95]),
          geometry:roi,
          scale:500,
          maxPixels:1e13
        });
    
    var min = ee.Number(num.get("NDVI_p5"));
    var max = ee.Number(num.get("NDVI_p95"));
    print("min",min);
    print("max",max)
    
    //quantile and combine
    var greaterPart = bestEVI.gt(max);
    var lessPart    = bestEVI.lt(min);
    var middlePart  =ee.Image(1).subtract(greaterPart).subtract(lessPart); // get the middle part
    Map.addLayer(greaterPart, EVIVis, "greaterPart", false);
    Map.addLayer(lessPart, EVIVis, "lessPart", false);
    Map.addLayer(middlePart, EVIVis, "middlePart", false);
    //calculate FVC
    var tempf1=bestEVI.subtract(min).divide(max.subtract(min));
    //ee.Image(constant): creat a constatnt image
    var FVC=ee.Image(1).multiply(greaterPart)
                  .add(ee.Image(0).multiply(lessPart))
                  .add(tempf1.multiply(middlePart));
    Map.addLayer(FVC,{},"FVC" + i);
    
    var MODIS_LULC = ee.ImageCollection("MODIS/061/MCD12Q1")
                   .filterBounds(roi)
                   .filterDate(i+'-01-01',i+'-12-31')
                   .select('LC_Type1');
    var image_LULC = MODIS_LULC.first().clip(roi);
    var LULC = image_LULC.where(image_LULC.gte(1).and(image_LULC.lte(8)), 2)
                               .where(image_LULC.eq(9).or(image_LULC.eq(10)), 3)
                               .where(image_LULC.eq(11).or(image_LULC.eq(15)).or(image_LULC.eq(17)), 4)
                               .where(image_LULC.eq(12).or(image_LULC.eq(14)), 1)
                               .where(image_LULC.eq(13), 5)
                               .where(image_LULC.eq(16), 6);
    
    
    var a = LULC.where(LULC.gte(1), -0.0438)
                .where(LULC.gte(2), -0.1535)
                .where(LULC.gte(3), -0.1151)
                .where(LULC.gte(4), 0)
                .where(LULC.gte(5), 0)
                .where(LULC.gte(6), -0.0768);
    
    var C = a.multiply(FVC).exp();
    
    var CVis = {
      min: 0.2,
      max: 1,
      palette: ['#86d68a', 'bbe029', '0a9501', '074b03', '#21cfd6', '#2d0ed6','#ca09d6'],
    };
    Map.addLayer(C,CVis,"C");
    
    // export to drive
    Export.image.toDrive(
      {
        image: C,
        description: 'C_' + i,
        folder: 'C',
        region: roi,
        scale: 500,
        crs: "EPSG:4326",
        maxPixels: 1e13,
      })
}

结果如下：