From b8a617ed8d511825d89e5959a5f44044baa5fa31 Mon Sep 17 00:00:00 2001 From: gis-xh Date: Thu, 15 Jan 2026 16:47:08 +0800 Subject: [PATCH] =?UTF-8?q?feat(GEE=5FScripts):=20=E6=96=B0=E5=A2=9E?= =?UTF-8?q?=E5=9C=9F=E5=9C=B0=E8=A6=86=E7=9B=96=E5=88=86=E7=B1=BB=E6=95=B0?= =?UTF-8?q?=E6=8D=AE=E4=B8=8B=E8=BD=BD=E8=84=9A=E6=9C=AC.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GEE_Scripts/LULC_download.js | 267 +++++++++++++++++++++++++++++++++++ 1 file changed, 267 insertions(+) create mode 100644 GEE_Scripts/LULC_download.js diff --git a/GEE_Scripts/LULC_download.js b/GEE_Scripts/LULC_download.js new file mode 100644 index 0000000..1c83e78 --- /dev/null +++ b/GEE_Scripts/LULC_download.js @@ -0,0 +1,267 @@ +/** + * 土地覆盖分类数据下载 + * + * @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, + }, + }); +}