diff --git a/GEE_Scripts/S1andS2_download.js b/GEE_Scripts/S1andS2_download.js new file mode 100644 index 0000000..7ee8193 --- /dev/null +++ b/GEE_Scripts/S1andS2_download.js @@ -0,0 +1,238 @@ +/** + * Sentinel-1 & Sentinel-2 哨兵一号与二号长时序数据下载 —— 适用于小范围区域年度数据获取 + * + * @author CVEO Team + * @date 2026-01-05 + * + * 1. 加载 Sentinel-1, Sentinel-2 数据与 Cloud Score+ 云掩膜, 以及 HLS 数据 + * 2. 合成年度Sentinel-1, Sentinel-2, HLS无云影像 + * 3. 使用 HLS 影像作为 Sentinel-2 影像的补充 + * 4. 合并Sentinel-1和Sentinel-2影像 + * 5. 导出COG云优化并填补缺失值的GeoTIFF影像 (小区域不分块) + * + * 注: 截止至 2026-01-05, GEE 仍未集成 2015-2017 年的 L2A 数据. + */ + +// 加载研究区域和影像 +var name_list = [ + "台北中山区大直要塞区", + "花莲县新城乡佳山基地", + "高雄左营区左营军港", + "高雄旗山区陆军第八军团指挥部", +]; +var target_index = 2; // 从 0 开始计数, 0-3 +var region_name = name_list[target_index]; +var region = ee.FeatureCollection( + demoPoints + .filter(ee.Filter.eq("Name", region_name)) + .geometry() + .buffer(2500) + .bounds() +); +var start_year = 2025; +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. +var CLEAR_THRESHOLD = 0.65; +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); +} + +// 加载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 + " S1dataset", 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 S2dataset = 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 targetS2 = sortedByPriority(S2dataset, region); +var s2_img = ee.Image(targetS2.mosaic()); +print(start_date + " - " + end_date + " S2dataset", S2dataset); +var s2_img_rgb = s2_img.select(["Red", "Green", "Blue"]); +var s2_img_frgb = s2_img.select(["NIR", "Red", "Green"]); + +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 HLS = ee.ImageCollection(HLSL30.merge(HLSS30)); +var targetHLS = sortedByPriority(HLS, region, "CLOUD_COVERAGE"); +var hls_img = ee.Image(targetHLS.mosaic()); + +var s1_vis = { + min: -25, + max: 5, +}; + +var true_rgb_vis = { + bands: ["Red", "Green", "Blue"], + min: 0.0, + max: 0.3, +}; + +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"); +Map.addLayer(s1_vh_img, s1_vis, start_year + " Sentinel-1 GRD VH"); +Map.addLayer(hls_img, false_rgb_vis, start_year + " HLS False RGB", false); +Map.addLayer(hls_img, true_rgb_vis, start_year + " HLS True RGB"); +Map.addLayer( + s2_img_frgb, + false_rgb_vis, + start_year + " Sentinel-2 False RGB", + false +); +Map.addLayer(s2_img_rgb, true_rgb_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; +Export.image.toDrive({ + image: processed_img, + description: "S1_S2_" + start_year + "_Target" + target_index, + folder: "Sentinel", + region: region, // 添加后会自动裁剪 + scale: 10, + maxPixels: 1e13, // GEE 最多支持 1e8 像素, 当超过时会自动分块 + fileFormat: "GeoTIFF", + // 导出COG云优化的GeoTIFF影像, 并明确设置缺失值为-9999.0 + formatOptions: { + cloudOptimized: true, + noData: -9999.0, + }, +});