基于GEE实现MOD09A1总悬浮物TSS的计算并可视化

这段代码在Google Earth Engine (GEE) 平台上实现了从MODIS卫星影像数据集中提取特定时间范围内的影像，进行去云处理、裁剪和总悬浮物 (Total Suspended Solids, TSS) 的计算，最终将处理后的影像数据可视化并导出。

总体功能：
该代码实现了对MODIS卫星影像数据进行去云处理、裁剪、计算总悬浮物（TSS）并进行可视化和导出的功能。通过对每年的影像数据进行处理，可以获得研究区内各年份的TSS变化情况，并导出处理后的影像数据供进一步分析使用。

数据集：使用的是MODIS卫星影像数据集，具体为 `MODIS/061/MOD09A1`。

大体流程和功能

1. 设置年份和研究区：
   -设置需要提取的年份（2010年）。
   -定义研究区 `ROI`。

2. 定义去云处理函数：
   - `MOD09A1(image)` 函数用于去除影像中的云和云阴影。

3. 定义像素筛选函数：
   - `filter(image)` 函数用于筛选像素数量大于50000的影像。

4. 筛选和处理MODIS影像数据：
   - 从MODIS数据集中筛选出给定时间范围内的影像数据。
   - 将筛选出的影像数据按年份进行分类，并计算每年的平均影像。

5. 计算三刺激值和TSS：
   - `getTSS(image)` 函数用于计算影像的三刺激值（X, Y, Z），色度角（hue），修正色度角（hueAngle），并最终计算TSS。

色度角的计算按照以下论文：

参考文献：李慧真，王雨辰，段高雨，等. 基于 GEE 的黄河口表层悬浮泥沙浓度时空分布及其影响因素分析[J]. 海洋学报，2023，45(8)：178–190， doi:10.12284/hyxb2023090

TSS公式是一个经验模型，根据色度角（hueAngle）计算TSS浓度。公式的形式和系数来源于遥感反演TSS浓度的研究和校准实验。

公式解释：

1. 色度角（hueAngle）：
   • 色度角是从遥感影像的反射率计算得到的一个参数，用于反映水体的光学特性。色度角的变化与水体中悬浮颗粒物的浓度有相关性。
2. 公式中的各项：
   • hueAngle.multiply(-0.07043)：对色度角进行线性变换，-0.07043是根据经验或校准数据得出的系数。
   • .exp()：对上述结果取指数。
   • .multiply(1160.51475)：对指数结果进行线性缩放，1160.51475是根据经验或校准数据得出的系数。
   • .add(1.2273)：加上一个偏置值，1.2273是根据经验或校准数据得出的常数。
   • .rename('TSS')：将计算结果重命名为'TSS'。

6. 可视化参数设置：
   - 定义了TSS的可视化参数 `visTSS`，包括色带和范围。

7. 裁剪和可视化：
   - 对计算出的TSS影像进行裁剪，并添加到地图上进行可视化显示。

8. 统计和导出：
   - 计算裁剪后影像的TSS均值。
   - 打印均值结果。
   - 将处理后的影像导出到Google Drive。

结果
- 地图显示：在GEE地图上显示经过TSS计算和去云处理后的影像。
- 统计输出：打印出裁剪后影像的TSS均值。
- 影像导出：将处理后的TSS影像导出到Google Drive，文件夹名称为"TSS"，文件名为“年份_TSS_MOD09A1”。

完整代码：

var year=2010; //设置需要提取的year
var ROI=table
get_ym_TSS(year);

function get_ym_TSS(year){
  
//去云
    function MOD09A1(image) {
    
      var cloudShadowBitMask = 1 << 3;
      var cloudsBitMask = 1 << 5;
      var qa = image.select('pixel_qa');
      var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
          .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
       return image.updateMask(mask);
    }
    
    function filter(image) {
  var PixelCount = image.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: ROI,
  scale: 500,
  maxPixels: 1e13
}).get('sur_refl_b01');
  return ee.Algorithms.If(ee.Number(PixelCount).gt(50000), 
                          image);
}
      //筛选 MOD09A1 数据
    var TSS09=ee.ImageCollection("MODIS/061/MOD09A1")
                 .filterBounds(ROI)
                 .filterDate('2020-09-01', '2020-09-30')
                 .map(function(image){
                    return image.set('year', ee.Image(image).date().get('year'))                           
                  })
    
for(var a = year; a < 2021; a++){
   var mean = TSS09.filterMetadata('year', 'equals', a).mean().clip(ROI)
 }
 

    //计算三刺激值
   var getTSS = function(image){
      var X = image.expression('1.1302 * b3 + 1.7517 * b4 + 2.7689 * b1',{
          b1: image.select('sur_refl_b01'),
          //b2: image.select('B2'),
          b3: image.select('sur_refl_b03'),
          b4: image.select('sur_refl_b04')
      }).rename('X');
    
      var Y = image.expression(' 0.0601 * b3 + 4.5907 * b4 + b1',{
          b1: image.select('sur_refl_b01'),
          //b2: image.select('B2'),
          b3: image.select('sur_refl_b03'),
          b4: image.select('sur_refl_b04')
    }).rename('Y');
    
      var Z = image.expression(' 5.5934 * b3 + 0.0565 * b4 ',{
          b1: image.select('sur_refl_b01'),
          //b2: image.select('B2'),
          b3: image.select('sur_refl_b03'),
          b4: image.select('sur_refl_b04')
    }).rename('Z');
    
    
    
     //计算色度角
      var chrx = X.divide((X.add(Y).add(Z))).rename('chrx');
      var chry = Y.divide((X.add(Y).add(Z))).rename('chry');
      var atan2 = chrx.subtract(1/3).atan2((chry).subtract(1/3)).rename('atan2');
      var hue_pre = atan2.multiply(ee.Number(180)).divide(Math.PI).rename('hue_pre');
      var hue = hue_pre < 0 ? (hue_pre + 360) : hue_pre.rename('hue');
      var a = hue.divide(100).rename('a');
      
     //修正色度角
      var cor = a.pow(5).multiply(33.086).add(a.pow(4).multiply(-221.6))
        .add(a.pow(3).multiply(512.65)).add(a.pow(2).multiply(-468.4))
        .add(a.multiply(137.27)).add(4.6374).rename('cor');
      var hueAngle = hue.add(cor).rename('hueAngle');
      
      
     //计算 
      var TSS = hueAngle.multiply(-0.07043).exp().multiply(1160.51475) .add(1.2273).rename('TSS');
      
      return image.addBands(TSS);
      
   }
  
      // 将计算的可视化，添加到图层
      var visTSS = {
        min: 0,
        max: 40,
        palette: [
               
              'A45308','A4600A','AF8A44','B3A053','AE9F5C', 'A8A965','ADB55F','AAB86D',
              'A5BC75', '94B660','95B645', '7DAE38', '7BA654','759E72','698C86','6D9298',
              '568F96','4B80A0','327CBB', '316DC5','2158BC',
        ]
    };
    
    
    
    //文件名 
    var year_folder = year//+''+month;
    
    var clipedTSS = getTSS(mean).select('TSS').multiply(1000).toInt().divide(1000);
    Map.addLayer(clipedTSS.clip(ROI),visTSS,year_folder+'TSS');
  
    // 减少地区，区域参数是Feature几何体
    var meanDictionary = clipedTSS.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry:ROI,
      scale: 30,
      maxPixels: 1e13
    });
    
    
     // 打印结果
    print(meanDictionary,year_folder);
    
    //输出到文件夹：FUI，文件名为：year_folder+"_CIE_MOD09"
    Export.image.toDrive({
      image:clipedTSS,
      scale:30,
      crs:'EPSG:32650',
      maxPixels : 1e13,
      region:ROI,
      folder:"TSS",
      description:year_folder+"_TSS_MOD09A1"
    });
 

}

运行结果：

工作台运行结果：