GEE随记（七）：结合GIS计算土壤侵蚀模型RULSE

Laurence529

已于 2024-06-02 14:21:46 修改

阅读量5.4k

点赞数 16

分类专栏： GEE随记 ArcGIS笔记文章标签：云计算

于 2024-06-02 14:16:37 首次发布

本文链接：https://blog.csdn.net/weixin_48012947/article/details/139390578

版权

GEE随记同时被 2 个专栏收录

7 篇文章

订阅专栏

ArcGIS笔记

3 篇文章

订阅专栏

1. RULSE模型介绍

RUSLE(Revised Universal Soil Loss Equation)，是修正的通用土壤流失方程。修正的通用土壤流失方程(RUSLE)是一种公认的、已得到广泛应用的经验土壤侵蚀估算模型。基于5个因子即：(1)降雨侵蚀力因子(R)、(2)土壤侵蚀力因子(K)、(3)坡长和陡度因子(LS)、(4)植被覆盖管理因子(C)、(5)水土保持实践因子(P)。

$A = R\times K\times LS\times C\times P$

式中：1）降雨侵蚀力因子(R)以动能的形式评价降雨影响的效果，并预测与该降水事件直接相关的径流量和径流量；

2）土壤可蚀性因子(K)表示在一个标准条件下，土壤抵抗雨滴影响所产生的侵蚀的能力，以及降雨影响所产生的径流量和径流量；

3）斜率长度(L)和陡度(S)因素反映了区域地形对土壤侵蚀率的影响；

4）作物经营管理因子(C)直接受植被类型、植被生长阶段和植被覆盖度的影响；

5）土壤保持措施因子(P)显示了实施措施将减少径流量和径流量的效果，从而降低了土壤侵蚀的数量和速率。

接下来将利用GEE计算R和C因子。

2. GEE计算R和C因子并导出计算LS因子所需的高程数据

此处参考推文：基于GEE的RULSE

//导入研究区
var roi = ee.FeatureCollection('projects/your_roi')
Map.centerObject(roi,12)

// 导入Landsat影像、降雨及高程数据集
//Landsat image
var Landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2");

//CHIRPS precipitation dataset -- daily (mm/day)
var CHIRPS_daily = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
           .select('precipitation') //select precipitation
           .map(function(image){return image.clip(roi.geometry())});

//ALOS DEM
var ALOS = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select('DSM')
             .map(function(image){return image.clip(roi.geometry())});
//Map.addLayer(ALOS.median().randomVisualizer())


//Landsat 8去云及波段定标
var cloudMaskL8 = function (image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 4);
  var cloudsBitMask = (1 << 3);
  var dilatedcloud = (1 << 1);
  // Get the pixel QA band.
  var qa = image.select('QA_PIXEL');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
               .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
               .and(qa.bitwiseAnd(dilatedcloud).eq(0));
  return image.addBands(opticalBands,null,true)
              .updateMask(mask).clip(roi.geometry());
}

//Calculate NDVI
function findndvi_l8(image){
  var vi = image.normalizedDifference(['SR_B5','SR_B4']).rename('ndvi')
  return vi.updateMask(vi.gt(-1).and(vi.lt(1))) 
}
//选择时间段                                           
var img2020 = Landsat8.filterDate('2019-12-31','2021-12-31')
                      .filterBounds(roi)
                      .map(cloudMaskL8).median();
Map.addLayer(img2020.select(['SR_B5','SR_B4','SR_B3']),{min:0,max:0.3})
var ndvi_img2020 = findndvi_l8(img2020).select(['ndvi']);
//Map.addLayer(ndvi_img2020.randomVisualizer())


//因子C的计算
// C ------- Vegetation Cover Factor, a=2, b=1,
// The higher the C, the worst the growth
var cf_2020 = ((ndvi_img2020.multiply(-2)).divide(ndvi_img2020.multiply(-1).add(1))).exp();



//R因子的计算
// Calculate R factor based on Flampouris' 
//phd thesis "Study of the effect of rainfall factor R on Rusle's law (in Greek)" 
//(2008) (No. GRI-2008-1712),Aristotle University of Thessaloniki.
//He used gauge stations and created some stable values 
//Retrieve CHIRPS percipitation data
var jan_chirps2019 = CHIRPS_daily.filterDate('2020-01-01','2020-02-01');
var feb_chirps2019 = CHIRPS_daily.filterDate('2020-02-01','2020-03-01');
var mar_chirps2019 = CHIRPS_daily.filterDate('2020-03-01','2020-04-01');
var apr_chirps2019 = CHIRPS_daily.filterDate('2020-04-01','2020-05-01');
var may_chirps2019 = CHIRPS_daily.filterDate('2020-05-01','2020-06-01');
var jun_chirps2019 = CHIRPS_daily.filterDate('2020-06-01','2020-07-01');
var jul_chirps2019 = CHIRPS_daily.filterDate('2020-07-01','2020-08-01');
var aug_chirps2019 = CHIRPS_daily.filterDate('2020-08-01','2020-09-01');
var sep_chirps2019 = CHIRPS_daily.filterDate('2020-09-01','2020-10-01');
var oct_chirps2019 = CHIRPS_daily.filterDate('2020-10-01','2020-11-01');
var nov_chirps2019 = CHIRPS_daily.filterDate('2020-11-01','2020-12-01');
var dec_chirps2019 = CHIRPS_daily.filterDate('2020-12-01','2021-01-01');     
       
// Calculate annual precipitation (mm) for 2020 year
var annual_precip = (jan_chirps2019.sum()).add(feb_chirps2019.sum())
                    .add(mar_chirps2019.sum()).add(apr_chirps2019.sum())
                    .add(may_chirps2019.sum()).add(jun_chirps2019.sum())
                    .add(jul_chirps2019.sum()).add(aug_chirps2019.sum())
                    .add(sep_chirps2019.sum()).add(oct_chirps2019.sum())
                    .add(nov_chirps2019.sum()).add(dec_chirps2019.sum());
var flap =  annual_precip.multiply(0.8);
Map.addLayer(flap.randomVisualizer())
var ALOS_roi = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select('DSM').median().clip(roi)
//Map.addLayer(ALOS_roi.randomVisualizer())


//cf_2020 for C factor; flap for R factor; ALOS_roi for elevation data
// export A value for each year (e.g. 2020);
Export.image.toDrive({
  image: ALOS_roi,
  description: '****',
  folder:'****',
  scale: 30,
  crs:'EPSG:4326',
  maxPixels: 1e13,
  region:roi.geometry()
});

3. ArcGIS中计算LS因子

此处参考推文：GEE系列（2）——GEE处理与arcgis 结合在土壤侵蚀上的应用

首先，通过3D Analysis中的Slope工具，利用上述下载的高程数据计算坡度roi_slope；其次，通过栅格计算器计算，公式如下：

$Con("roislope"< 5,10.8 * Sin("roislope"*3.14/180) + 0.03,Con("roislope" <10,16.8 * Sin("roislope"*3.14/180) -0.5,21.91* Sin("roislope"*3.14/180)-0.96 ))$