/**
 * Sentinel-1 & Sentinel-2 哨兵一号与二号长时序数据下载 —— 适用于小范围区域年度数据获取
 *
 * @author CVEO Team
 * @date 2026-01-07
 *
 * 1. 加载 Sentinel-1 GRD, Sentinel-2 L2A SR数据与 Cloud Score+ 云掩膜, 以及 HLS 数据
 * 2. 合成年度Sentinel-1, Sentinel-2, HLS无云影像
 * 3. 使用 Sentinel-2 L1C TOA 影像与 HLS 影像作为 Sentinel-2 L2A SR 影像的补充
 * 4. 合并 Sentinel-1 和 Sentinel-2 影像
 * 5. 导出 COG 云优化并填补缺失值的 GeoTIFF 影像 (小区域不分块)
 *
 * 注: 截止至 2026-01-07, GEE 仍未集成 2015-2017 年的 Sentinel-2 L2A SR数据, 且 2018 年的 L2A 数据仍不完整.
 */

// 加载研究区域和影像
var name_list = [
  "台北中山区大直要塞区",
  "花莲县新城乡佳山基地",
  "高雄左营区左营军港",
  "高雄旗山区陆军第八军团指挥部",
];
var target_index = 1; // 从 0 开始计数, 0-3
var region_name = name_list[target_index];
var region = ee.FeatureCollection(
  demoPoints
    .filter(ee.Filter.eq("Name", region_name))
    .geometry()
    .buffer(3000)
    .bounds()
);
var start_year = 2015;
var end_year = 2025;
var start_date = start_year + "-" + "01-01";
var end_date = start_year + "-" + "12-31";
var cloud_threshold = 100; // 最大云量阈值
// 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.
// 内陆建议设置为 0.8, 沿海或大江大湖旁建议设置为 0.6
var CLEAR_THRESHOLD = 0.8;
var L30Bands = ["B2", "B3", "B4", "B5", "B6", "B7"];
// 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 region_geo = region.geometry();
var bounds = region_geo.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
    .divide(10000)
    .copyProperties(image)
    .copyProperties(image, ["system:time_start", "system:index", "system:id"]);
}

/**
 * Function to mask clouds using the HLS Fmask band
 * @param {ee.Image} image HLS image
 * @return {ee.Image} cloud masked HLS image
 */
function maskHLSclouds(image) {
  var qa = image.select("Fmask");
  // Bits 1 and 2 are cloud and shadow, respectively.
  var cloudBitMask = 1 << 1;
  var shadowBitMask = 1 << 2;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa
    .bitwiseAnd(cloudBitMask)
    .eq(0)
    .and(qa.bitwiseAnd(shadowBitMask).eq(0));
  image = image.updateMask(mask);
  return image
    .copyProperties(image)
    .copyProperties(image, ["system:time_start", "system:index", "system:id"]);
}

/**
 * 按照属性值对影像集合进行排序
 * @param {ee.ImageCollection} collection 影像集合
 * @param {ee.Geometry} region 目标区域
 * @param {string} cloud_field 云量属性字段名, 默认为 "CLOUDY_PIXEL_PERCENTAGE"
 * @return {ee.ImageCollection} 排序后的影像集合
 */
function sortedByPriority(collection, region, cloud_field) {
  cloud_field = cloud_field || "CLOUDY_PIXEL_PERCENTAGE";
  // 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(cloud_field, false).sort("coverage_area");
  return ee.ImageCollection(sorted);
}

// 加载HLS L30和L30数据
var HLSL30 = ee
  .ImageCollection("NASA/HLS/HLSL30/v002")
  .filter(common_filter)
  .filter(ee.Filter.lt("CLOUD_COVERAGE", cloud_threshold))
  .map(maskHLSclouds)
  .select(L30Bands, commonBands);
var HLSS30 = ee
  .ImageCollection("NASA/HLS/HLSS30/v002")
  .filter(common_filter)
  .filter(ee.Filter.lt("CLOUD_COVERAGE", cloud_threshold))
  .map(maskHLSclouds)
  .select(s2Bands, commonBands);
var HLSdataset = ee.ImageCollection(HLSL30.merge(HLSS30));
print(start_date + " - " + end_date + " HLS dataset", HLSdataset);
var targetHLSCol = sortedByPriority(HLSdataset, region, "CLOUD_COVERAGE");
// Mosaic 处理后会丢失投影信息, 需重新设置投影
var proj = targetHLSCol.first().select(1).projection();
var hls_img = ee.Image(targetHLSCol.mosaic()).setDefaultProjection(proj);
var hls_10m_img = hls_img.resample("bicubic").reproject({
  crs: proj,
  scale: 10,
});

// 加载Sentinel-1 GRD数据
var S1dataset = ee
  .ImageCollection("COPERNICUS/S1_GRD")
  .filter(common_filter)
  .filter(ee.Filter.eq("instrumentMode", "IW"));
var S1VVdataset = S1dataset.filter(
  ee.Filter.listContains("transmitterReceiverPolarisation", "VV")
).select("VV");
var S1VHdataset = S1dataset.filter(
  ee.Filter.listContains("transmitterReceiverPolarisation", "VH")
).select("VH");
print(start_date + " - " + end_date + " Sentinel-1 dataset", S1dataset);
var s1_vv_img = ee.Image(S1VVdataset.median());
var s1_vh_img = ee.Image(S1VHdataset.median());

// 加载Sentinel-2 L2A数据与Cloud Score+云掩膜
// 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.
// 加载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 S2TOAdataset = ee
  .ImageCollection("COPERNICUS/S2_HARMONIZED")
  .filter(common_filter)
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", cloud_threshold))
  .linkCollection(S2csPlus, [QA_BAND])
  .map(maskS2cloudsByCS)
  .select(s2Bands, commonBands);
var targetS2TOACol = sortedByPriority(S2TOAdataset, region);
var s2_toa_img = ee.Image(targetS2TOACol.mosaic());
// 加载SIAC工具集, 容易爆内存
// var SIAC = require("users/marcyinfeng/utils:SIAC");
// var s2_boa_img = SIAC.get_sur(s2_toa_img.clip(region));
print(start_date + " - " + end_date + " Sentinel-2 TOA dataset", S2TOAdataset);

var S2SRdataset = ee
  .ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filter(common_filter)
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", cloud_threshold))
  .linkCollection(S2csPlus, [QA_BAND])
  .map(maskS2cloudsByCS)
  .select(s2Bands, commonBands);
var targetS2SRCol = sortedByPriority(S2SRdataset, region);
var s2_sr_img = ee.Image(targetS2SRCol.mosaic());
print(start_date + " - " + end_date + " Sentinel-2 SR dataset", S2SRdataset);

var s2_img = ee.Image();
// 2015-2018 年间若 Sentinel-2 SR 数据为空, 则使用 Sentinel-2 TOA 数据
// 若仍然没有, 则使用重采样到 10m 的 HLS 数据
if (S2SRdataset.size().getInfo() === 0 && S2TOAdataset.size().getInfo() !== 0) {
  s2_img = s2_toa_img;
} else if (S2TOAdataset.size().getInfo() === 0 || start_year === 2015) {
  s2_img = hls_10m_img;
} else {
  s2_img = s2_sr_img;
}

// 针对试点区域的特殊设置
if (start_year === 2018 && (target_index === 0 || target_index === 3)) {
  s2_img = s2_toa_img;
}
if (start_year === 2016 && target_index === 1) {
  s2_img = hls_10m_img;
}

var s2_img_rgb = s2_img.select(["Red", "Green", "Blue"]);
var s2_img_frgb = s2_img.select(["NIR", "Red", "Green"]);

var s1_vis = {
  min: -25,
  max: 5,
};

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

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

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

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

Map.centerObject(region, 14);
Map.addLayer(s1_vv_img, s1_vis, start_year + " Sentinel-1 GRD VV", false);
Map.addLayer(s1_vh_img, s1_vis, start_year + " Sentinel-1 GRD VH", false);
Map.addLayer(hls_img, true_rgb_vis, start_year + " HLS True RGB 30m");
Map.addLayer(hls_10m_img, true_rgb_vis, start_year + " HLS True RGB 10m");
Map.addLayer(
  s2_img_frgb,
  false_rgb_vis,
  start_year + " Sentinel-2 False RGB",
  false
);
Map.addLayer(s2_img_rgb, s2_vis, start_year + " Sentinel-2 True RGB");
Map.addLayer(region.style(styling), {}, region_name);

// 导出合并后的影像
// 明确设置数据类型为Float32, 否则默认类型为Float64, 会占用更多内存
// 并对缺失值进行填充, 否则默认为nan不便于后续本地处理
var processed_img = s2_img
  .addBands(s1_vv_img)
  .addBands(s1_vh_img)
  .toFloat()
  .unmask(-9999.0);
print(
  "Start exporting " + start_year + " Sentinel-2 and S1 GRD VV/VH merged image",
  processed_img
);
target_index = target_index + 1;
var target_name = "Target" + target_index;
Export.image.toDrive({
  image: processed_img,
  description: "S1_S2_" + target_name + "_" + start_year,
  folder: "Sentinel",
  region: region, // 添加后会自动裁剪
  scale: 10,
  crs: "EPSG:4326",
  maxPixels: 1e13, // GEE 最多支持 1e8 像素, 当超过时会自动分块
  fileFormat: "GeoTIFF",
  // 导出COG云优化的GeoTIFF影像, 并明确设置缺失值为-9999.0
  formatOptions: {
    cloudOptimized: true,
    noData: -9999.0,
  },
});