Sentinel2 NDWI+NDVI gee水体提取

Question

在GEE平台上分别利用Landsat和哨兵系列遥感影像对数据进行统一预处理预处理，然后提取水体面积结合DEM反演水位并构建反演模型来估算断面流量

这是之前尝试过的代码，但有部分出错且结果不对，需要更改帮助
Sentinel2 NDWI+NDVI 水体提取

// 1. 初始化检查
if (typeof roi === 'undefined' || typeof dam === 'undefined') {
  print('🛑 【严重错误】: 请先在地图上画 roi 和 dam！');
} else {
  Map.centerObject(roi, 11);
  Map.setOptions("TERRAIN");
}
 
// 2. 参数设置
var START_DATE = '2016-01-01'; 
var END_DATE = '2025-12-29';   
var MAX_CLOUD_PERCENT = 20;    // 改MIN为MAX，语义更对
 
// 3. 核心函数（优化阈值+去云）
function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  
  return image.updateMask(mask).divide(1e4)
      .select("B.*")
      .copyProperties(image, ["system:time_start"]); 
}
 
var addNDWI = function(image) {
  return image.addBands(image.normalizedDifference(['B3', 'B8']).rename('NDWI'));
};
 
var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
};
 
// 优化组合掩膜：NDWI≥0.1 + NDVI≤0.2（更贴合三峡库区）
var CombinedMask = function(image) {
  var ndwi = image.select('NDWI');
  var ndvi = image.select('NDVI');
  var waterMask = ndwi.gte(0.1).and(ndvi.lte(0.2)).rename('Water');
  return image.addBands(waterMask);
};
 
// 4. 数据加载与按日合成
var s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
  .filterBounds(roi)
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', MAX_CLOUD_PERCENT))
  .filterDate(START_DATE, END_DATE)
  .map(function(image) { return image.clip(roi); });
 
function mosaicBy(imcol){
  var imlist = imcol.toList(imcol.size());
  var all_dates = imlist.map(function(im){ return ee.Image(im).date().format("YYYY-MM-dd"); });
  var unique_dates = all_dates.distinct();
  
  var mosaic_imlist = unique_dates.map(function(d){
    var date1 = ee.Date(d);
    var im = imcol
      .filterDate(date1, date1.advance(1, "day"))
      .mosaic(); 
    return im.set(
        "system:time_start", date1.millis(),
        "date_str", d 
    );
  });
  return ee.ImageCollection(mosaic_imlist);
}
 
var s2day = mosaicBy(s2);
var S2_Processed = s2day
    .map(maskS2clouds)
    .map(addNDWI)      
    .map(addNDVI)      
    .map(CombinedMask)
    .map(function(img) {
        var count = img.select('Water').reduceRegion({
            reducer: ee.Reducer.count(), 
            geometry: roi, 
            scale: 200, 
            maxPixels: 1e13
        }).get('Water');
        return img.set('valid_pixels', count);
    })
    .filter(ee.Filter.gt('valid_pixels', 100)); 
 
print('Sentinel2 有效影像数量:', S2_Processed.size());
 
// 5. 面积与水位计算
var Water_Calc = function(img) {
  var water = img.select('Water');
  var srtm = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation'); 
  
  var area = water.multiply(ee.Image.pixelArea()).divide(1e6)
             .reduceRegion({
               reducer: ee.Reducer.sum(),
               geometry: roi,
               scale: 30, 
               bestEffort: true,
               maxPixels: 1e13
             }).get('Water');
 
  var srtm_masked = srtm.updateMask(water).updateMask(srtm.gte(145).and(srtm.lte(175)));
  var mean_elev = srtm_masked.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: dam,
      scale: 30,
      bestEffort: true,
      maxPixels: 1e13
    }).get('Elevation');
    
  var water_level = ee.Algorithms.If(mean_elev, mean_elev, 0);
  var dateStr = img.get('date_str');
 
  return ee.Feature(null, {
    'date': dateStr,
    'area_km2': area,
    'water_level_m': water_level
  });
};
var resultTable = S2_Processed.map(Water_Calc);
 
// 6. 导出（修复语法错误）
Export.table.toDrive({
  collection: resultTable,
  description: 'ThreeGorges_Water_Data_Sentinel2',
  fileFormat: 'CSV',
  selectors: ['date', 'area_km2', 'water_level_m']
});

Answer 1

晚上好🌙🌙🌙
本答案参考通义千问

你提供的代码是用于在 Google Earth Engine (GEE) 平台上利用 Sentinel-2 数据提取水体面积，并结合 DEM（SRTM） 反演水位，进而估算断面流量的流程。然而，你在描述中提到存在“部分出错且结果不对”，需要进行修改和优化。

---

✅ 问题分析与关键错误点

1. Water_Calc 函数中的 ee.Algorithms.If 使用不当

   • ee.Algorithms.If 需要一个条件表达式、真值和假值，但当前写法没有正确传递参数。
   • 应该使用 ee.Number(mean_elev).gt(0) 来判断是否有效。

2. reduceRegion 的 geometry 参数不一致

   • 在 area 计算时使用了 roi，而在 mean_elev 中使用了 dam，这可能导致计算范围不一致或结果异常。

3. srtm_masked 未正确设置掩膜

   • updateMask(water) 是正确的，但 srtm.gte(145).and(srtm.lte(175)) 会导致所有像素都被过滤掉，除非 DEM 值在此范围内。

4. water_level 没有正确赋值

   • 如果 mean_elev 为 null 或无效值，water_level 会被设为 0，但应更安全地处理这种情况。

5. 导出数据格式可能有问题

   • Export.table.toDrive 的 selectors 会从 FeatureCollection 中提取字段，需确保字段名称正确。

---

🔧 修正后的完整代码如下：

// 1. 初始化检查
if (typeof roi === 'undefined' || typeof dam === 'undefined') {
  print('🛑 【严重错误】: 请先在地图上画 roi 和 dam！');
} else {
  Map.centerObject(roi, 11);
  Map.setOptions("TERRAIN");
}

// 2. 参数设置
var START_DATE = '2016-01-01'; 
var END_DATE = '2025-12-29';   
var MAX_CLOUD_PERCENT = 20;    // 改MIN为MAX，语义更对

// 3. 核心函数（优化阈值+去云）
function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  
  return image.updateMask(mask).divide(1e4)
      .select("B.*")
      .copyProperties(image, ["system:time_start"]); 
}

var addNDWI = function(image) {
  return image.addBands(image.normalizedDifference(['B3', 'B8']).rename('NDWI'));
};

var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
};

// 优化组合掩膜：NDWI≥0.1 + NDVI≤0.2（更贴合三峡库区）
var CombinedMask = function(image) {
  var ndwi = image.select('NDWI');
  var ndvi = image.select('NDVI');
  var waterMask = ndwi.gte(0.1).and(ndvi.lte(0.2)).rename('Water');
  return image.addBands(waterMask);
};

// 4. 数据加载与按日合成
var s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
  .filterBounds(roi)
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', MAX_CLOUD_PERCENT))
  .filterDate(START_DATE, END_DATE)
  .map(function(image) { return image.clip(roi); });

function mosaicBy(imcol){
  var imlist = imcol.toList(imcol.size());
  var all_dates = imlist.map(function(im){ return ee.Image(im).date().format("YYYY-MM-dd"); });
  var unique_dates = all_dates.distinct();
  
  var mosaic_imlist = unique_dates.map(function(d){
    var date1 = ee.Date(d);
    var im = imcol
      .filterDate(date1, date1.advance(1, "day"))
      .mosaic(); 
    return im.set(
        "system:time_start", date1.millis(),
        "date_str", d 
    );
  });
  return ee.ImageCollection(mosaic_imlist);
}

var s2day = mosaicBy(s2);
var S2_Processed = s2day
    .map(maskS2clouds)
    .map(addNDWI)      
    .map(addNDVI)      
    .map(CombinedMask)
    .map(function(img) {
        var count = img.select('Water').reduceRegion({
            reducer: ee.Reducer.count(), 
            geometry: roi, 
            scale: 200, 
            maxPixels: 1e13
        }).get('Water');
        return img.set('valid_pixels', count);
    })
    .filter(ee.Filter.gt('valid_pixels', 100)); 

print('Sentinel2 有效影像数量:', S2_Processed.size());

// 5. 面积与水位计算
var Water_Calc = function(img) {
  var water = img.select('Water');
  var srtm = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation'); 
  
  var area = water.multiply(ee.Image.pixelArea()).divide(1e6)
             .reduceRegion({
               reducer: ee.Reducer.sum(),
               geometry: roi,
               scale: 30, 
               bestEffort: true,
               maxPixels: 1e13
             }).get('Water');

  var srtm_masked = srtm.updateMask(water).updateMask(srtm.gte(145).and(srtm.lte(175)));
  var mean_elev = srtm_masked.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: roi, // 修改为 roi，统一计算区域
      scale: 30,
      bestEffort: true,
      maxPixels: 1e13
    }).get('Elevation');

  // 安全处理 mean_elev 为 null 的情况
  var water_level = ee.Algorithms.If(
    ee.Number(mean_elev).gt(0), 
    mean_elev, 
    0
  );

  var dateStr = img.get('date_str');

  return ee.Feature(null, {
    'date': dateStr,
    'area_km2': area,
    'water_level_m': water_level
  });
};

var resultTable = S2_Processed.map(Water_Calc);

// 6. 导出（修复语法错误）
Export.table.toDrive({
  collection: resultTable,
  description: 'ThreeGorges_Water_Data_Sentinel2',
  fileFormat: 'CSV',
  selectors: ['date', 'area_km2', 'water_level_m']
});

---

✅ 重点修改说明

1. 修正 water_level 赋值逻辑

• 将原来的：

  var water_level = ee.Algorithms.If(mean_elev, mean_elev, 0);
  

  修改为：

  var water_level = ee.Algorithms.If(
    ee.Number(mean_elev).gt(0), 
    mean_elev, 
    0
  );
  

  • 确保只有当 mean_elev 不为 null 且大于 0 时才赋值，否则设为 0。

---

2. 统一 geometry 参数

• 将 mean_elev 的 geometry 从 dam 改为 roi，确保与 area 计算一致，避免因范围不同导致结果偏差。

---

3. 增加对 mean_elev 为 null 的处理

• 使用 ee.Number(mean_elev).gt(0) 判断是否为有效值，避免出现空值导致后续计算出错。

---

4. 优化 srtm_masked 的掩膜逻辑

• 如果 srtm 的高程范围（如 145~175）不合理，可以适当调整，例如：

  var srtm_masked = srtm.updateMask(water).updateMask(srtm.gte(100).and(srtm.lte(200)));
  

---

📌 后续建议

• 验证 roi 和 dam 是否正确绘制：确保两个几何对象都已定义并可用。
• 测试单个影像：在 GEE 控制台中运行 S2_Processed.first() 查看其结构和波段信息，确认 Water 波段是否生成。
• 添加日志输出：在关键步骤后添加 print() 输出，例如：

  print('Mean elevation:', mean_elev);
  print('Water level:', water_level);
  

  有助于调试。

---

如果你还有其他具体错误信息或结果不符合预期，请提供详细描述，我可以进一步帮助你定位问题。