/**
 * 土地覆盖分类数据下载
 *
 * @author CVEO Team
 * @date 2025-10-27
 *
 * 1. 加载 ESA WorldCover, Dynamic World 等公开土地覆盖数据
 * - ESA WorldCover 2021 V200 版本: 11 类地物
 * - Dynamic World V1 版本: 9 类地物
 * 2. 加载 Sentinel-2 数据
 * 3. 导出 COG 云优化并填补缺失值的 GeoTIFF 影像 (大区域 GEE 自动分块下载)
 */

// 加载研究区域和影像
// var region_name = "武汉市";
// var region_name_en = "Wuhan";
// var region = allCites.filter(ee.Filter.eq("市", region_name));
var region_name = "49REL";
var region_name_en = region_name;
var region = grid;
var year = 2025;
var start_date = year + "-01-01";
var end_date = year + "-12-31";
var target_month = 5; // 指定月份
var cloud_threshold = 60;
// Use 'cs' or 'cs_cdf', depending on your use case; see docs for guidance.
var QA_BAND = "cs_cdf";
// The threshold for masking; values between 0.50 and 0.65 generally work well.
// Higher values will remove thin clouds, haze & cirrus shadows.
var CLEAR_THRESHOLD = 0.7;
// Sentinel-2的B8A波段(20m)可能更适合代表近红外波段
// var s2Bands = ["B2", "B3", "B4", "B8", "B11", "B12"];
var s2Bands = ["B2", "B3", "B4", "B8A", "B11", "B12"];
var commonBands = ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"];
var bounds = region.geometry().bounds();
var common_filter = ee.Filter.and(
  ee.Filter.bounds(bounds),
  ee.Filter.date(start_date, end_date)
);

/**
 * 基于 Cloud Score+ 云掩膜
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2cloudsByCS(image) {
  var cloudMask = image.select(QA_BAND).gte(CLEAR_THRESHOLD);
  image = image.updateMask(cloudMask);
  return image
    .copyProperties(image)
    .copyProperties(image, ["system:time_start", "system:index", "system:id"]);
}

/**
 * 按照属性值对影像集合进行排序
 * @param {ee.ImageCollection} collection 影像集合
 * @param {ee.Geometry} region 目标区域
 * @return {ee.ImageCollection} 排序后的影像集合
 */
function sortedByPriority(collection, region) {
  // 1. 计算每张影像的有效像素 (非空像素) 覆盖面积
  var withArea = collection.map(function (image) {
    var area = image
      .mask()
      .reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: region,
        scale: 10,
        maxPixels: 1e9,
      })
      .values()
      .get(0);
    return image.set("coverage_area", area);
  });

  // 2. 剔除覆盖面积为 0 的影像
  withArea = withArea.filter(ee.Filter.gt("coverage_area", 0));

  // 3. 根据云量 (降序) 和覆盖率 (升序) 对影像进行综合排序
  // 考虑影像集合是以栈的形式 (先入后出) 组织的, 在进行 mosaic 合成时, 从最后一张影像开始合并, 逆转排序
  // sort 默认升序 True
  var sorted = withArea
    .sort("CLOUDY_PIXEL_PERCENTAGE", false)
    .sort("coverage_area");
  return ee.ImageCollection(sorted);
}

// 加载Sentinel-2 L2A数据
// Cloud Score+ image collection.
// Note Cloud Score+ is produced from Sentinel-2 level 1C data and can be applied to either L1C or L2A collections.
var S2csPlus = ee.ImageCollection("GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED");
var S2dataset = ee
  .ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filter(common_filter)
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", cloud_threshold))
  .select(s2Bands, commonBands)
  .linkCollection(S2csPlus, [QA_BAND])
  .map(maskS2cloudsByCS);
var target_s2 = sortedByPriority(
  S2dataset.filter(
    ee.Filter.calendarRange(target_month, target_month, "month")
  ),
  region
);
var s2_img = ee.Image(target_s2.median());

// 加载ESA WorldCover数据
var ESA_worldcover = ee.ImageCollection("ESA/WorldCover/v200").first();
var wc_img = ee.Image(ESA_worldcover.clip(region));
var wc_class_names = {
  "Tree cover": 10,
  Shrubland: 20,
  Grassland: 30,
  Cropland: 40,
  "Built-up": 50,
  "Bare / sparse vegetation": 60,
  "Snow and ice": 70,
  "Permanent water bodies": 80,
  "Herbaceous wetland": 90,
  Mangroves: 95,
  "Moss and lichen": 100,
};
// 统计各类别覆盖面积
var wc_class_area = wc_img
  .reduceRegion({
    reducer: ee.Reducer.frequencyHistogram(),
    geometry: region,
    scale: 10,
    maxPixels: 1e13,
  })
  .get("Map");
print("ESA WorldCover 2021 v200");
print("各类别名称:", wc_class_names);
print("各类别覆盖面积:", wc_class_area);

// 加载Dynamic World数据
var DWdataset = ee
  .ImageCollection("GOOGLE/DYNAMICWORLD/V1")
  .filter(common_filter)
  .select("label");
// mode() 取出现次数最多的类别作为合成后的类别
var dw_img = ee.Image(DWdataset.mode().clip(region));
var dw_class_names = {
  water: 0,
  trees: 1,
  grass: 2,
  flooded_vegetation: 3,
  crops: 4,
  shrub_and_scrub: 5,
  built: 6,
  bare: 7,
  snow_and_ice: 8,
};
var built_img = dw_img.updateMask(dw_img.eq(dw_class_names.built));
var crop_img = dw_img.updateMask(dw_img.eq(dw_class_names.crops));
var flood_img = dw_img.updateMask(dw_img.eq(dw_class_names.flooded_vegetation));
// 统计各类别覆盖面积
var dw_class_area = dw_img
  .reduceRegion({
    reducer: ee.Reducer.frequencyHistogram(),
    geometry: region,
    scale: 10,
    maxPixels: 1e13,
  })
  .get("label");
print("Dynamic World v1");
print("各类别名称:", dw_class_names);
print("各类别覆盖面积:", dw_class_area);

var VIS_PALETTE = [
  "419bdf",
  "397d49",
  "88b053",
  "7a87c6",
  "e49635",
  "dfc35a",
  "c4281b",
  "a59b8f",
  "b39fe1",
];
var DW_vis = {
  min: 0,
  max: 8,
  palette: VIS_PALETTE,
};

var true_rgb_vis = {
  bands: ["Red", "Green", "Blue"],
  min: 0.0,
  max: 3000,
};

var false_rgb_vis = {
  min: 0.0,
  max: 4000,
  gamma: 1.4,
  bands: ["NIR", "Red", "Green"],
};

var styling = {
  color: "red",
  fillColor: "00000000",
};

var WC_vis = {
  bands: ["Map"],
};

Map.centerObject(region, 9);
Map.addLayer(region.style(styling), {}, region_name);
Map.addLayer(wc_img, WC_vis, "ESA 2021 Landcover");
Map.addLayer(dw_img, DW_vis, "Dynamic World LULC");
Map.addLayer(
  s2_img,
  true_rgb_vis,
  year + "-" + target_month + " Sentinel-2 True RGB",
  false
);
Map.addLayer(
  s2_img,
  false_rgb_vis,
  year + "-" + target_month + " Sentinel-2 False RGB",
  false
);
// Map.addLayer(built_img, DW_vis, 'Built DW LULC');
// Map.addLayer(crop_img, DW_vis, "Crop DW LULC");
// Map.addLayer(flood_img, DW_vis, "Flood DW LULC");

// 导出裁剪处理后的土地覆盖数据
var fileName = region_name_en + "_ESA_WorldCover_2021_v200";
// 区域过大, 需要先将 wc_img 划分为四块网格分开下载
var bounds = region.geometry().bounds();
var coords = ee.List(bounds.coordinates().get(0));

var minX = ee.Number(ee.List(coords.get(0)).get(0));
var minY = ee.Number(ee.List(coords.get(0)).get(1));
var maxX = ee.Number(ee.List(coords.get(2)).get(0));
var maxY = ee.Number(ee.List(coords.get(2)).get(1));

var midX = minX.add(maxX).divide(2);
var midY = minY.add(maxY).divide(2);

var grid1 = ee.Geometry.Rectangle([minX, minY, midX, midY]);
var grid2 = ee.Geometry.Rectangle([midX, minY, maxX, midY]);
var grid3 = ee.Geometry.Rectangle([minX, midY, midX, maxY]);
var grid4 = ee.Geometry.Rectangle([midX, midY, maxX, maxY]);

var grids = [grid1, grid2, grid3, grid4];

for (var i = 0; i < grids.length; i++) {
  var grid_region = grids[i];
  var clipped_img = wc_img.clip(grid_region);
  var grid_fileName = fileName + "_grid_" + (i + 1);
  Export.image.toDrive({
    image: clipped_img,
    description: grid_fileName,
    folder: "LULC",
    region: grid_region,
    scale: 10,
    maxPixels: 1e13, // GEE 最多支持 1e8 像素, 当超过时会自动分块
    fileFormat: "GeoTIFF",
    // 导出COG云优化的GeoTIFF影像
    formatOptions: {
      cloudOptimized: true,
    },
  });
}