我想做提取耕地提取,想到了一篇董金玮老师的一篇论文,这个论文是先提取的耕地,再做作物分类,论文代码是开源的。基于这个代码下载数据。
但这个代码直接在云端上进行分类,GEE会爆内存,因此我准备把数据下载到本地,使用GPU加速进行XGboot提取耕地。
董老师的代码涉及到了100多个波段特征,我删减到了45个波段,然后分块进行了数据下载:
数据下载代码:
// ======================================== // 1. 初始化与区域选择 // ======================================== // 选择第一个区域作为AOI var aoiFeature = fenqu.first();var aoi = aoiFeature.geometry();// 可视化AOI(可选) Map .addLayer(aoi, {color : 'blue' }, 'AOI' );// 中心定位到AOI,缩放级别10(可选) Map .centerObject(aoi, 10 );// ======================================== // 2. 划分AOI为16个块 // ======================================== // 定义划分块数(4x4网格) var numCols = 4 ;var numRows = 4 ;// 获取AOI的边界和范围 var aoiBounds = aoi.bounds();var coords = ee.List(aoiBounds.coordinates().get(0 ));var xMin = ee.Number(ee.List(coords.get(0 )).get(0 ));var yMin = ee.Number(ee.List(coords.get(0 )).get(1 ));var xMax = ee.Number(ee.List(coords.get(2 )).get(0 ));var yMax = ee.Number(ee.List(coords.get(2 )).get(1 ));// 计算AOI的宽度和高度 var aoiWidth = xMax.subtract(xMin);var aoiHeight = yMax.subtract(yMin);// 计算每个块的宽度和高度 var tileWidth = aoiWidth.divide(numCols);var tileHeight = aoiHeight.divide(numRows);// 要排除的块的ID var excludeTiles = ['0_3' , '0_2' , '3_0' ]; // 左上角、第二行第一个、右下角 // 生成4x4网格,但排除特定块 var grid = ee.FeatureCollection( ee.List.sequence(0 , numCols - 1 ).map(function (col ) { return ee.List.sequence(0 , numRows - 1 ).map(function (row ) { var tileId = ee.String(col).cat('_' ).cat(ee.String(row)); var xmin = xMin.add(tileWidth.multiply(ee.Number(col))); var ymin = yMin.add(tileHeight.multiply(ee.Number(row))); var xmax = xmin.add(tileWidth); var ymax = ymin.add(tileHeight); var rectangle = ee.Geometry.Rectangle([xmin, ymin, xmax, ymax]); return ee.Feature(rectangle, { 'tile' : tileId }); }); }).flatten() ).filter(ee.Filter.inList('tile' , excludeTiles).not());// 可视化网格(可选) Map .addLayer(grid, {color : 'red' }, 'Grid' );// ======================================== // 3. 定义数据处理和导出函数 // ======================================== function processAndExport (tileFeature ) { var tileID = ee.String(tileFeature.get('tile' )); print('Processing Tile:' , tileID); var region = tileFeature.geometry(); // 2. 定义时间范围、波段及区域 var year = 2023 ; var startDate = ee.Date.fromYMD(year, 1 , 1 ); var endDate = ee.Date.fromYMD(year, 12 , 31 ); var bands = ['B2' , 'B3' , 'B4' , 'B8' ]; // 蓝、绿、红、近红外 // 3. 云掩膜函数:基于SCL波段 function maskS2clouds (image ) { var scl = image.select('SCL' ); // SCL分类值: 3(云)、8(阴影云) var cloudMask = scl.neq(3 ).and(scl.neq(8 )); return image.updateMask(cloudMask) .clip(region) .copyProperties(image, ["system:time_start" ]); } // 4. 添加光谱指数函数 function addSpectralIndices (image ) { // 计算NDVI var ndvi = image.normalizedDifference(['B8' , 'B4' ]).rename('NDVI' ); // 计算EVI var evi = image.expression( '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))' , { 'NIR' : image.select('B8' ), 'RED' : image.select('B4' ), 'BLUE' : image.select('B2' ) } ).rename('EVI' ); // 计算GNDVI var gndvi = image.normalizedDifference(['B8' , 'B3' ]).rename('GNDVI' ); // 计算SAVI var savi = image.expression( '((NIR - RED) / (NIR + RED + 0.5)) * 1.5' , { 'NIR' : image.select('B8' ), 'RED' : image.select('B4' ) } ).rename('SAVI' ); // 计算MSAVI2 var msavi2 = image.expression( '0.5 * (2 * NIR + 1 - sqrt((2 * NIR + 1)**2 - 8 * (NIR - RED)))' , { 'NIR' : image.select('B8' ), 'RED' : image.select('B4' ) } ).rename('MSAVI2' ); // 计算NDWI var ndwi = image.normalizedDifference(['B3' , 'B8' ]).rename('NDWI' ); // 计算NDSI var ndsi = image.normalizedDifference(['B3' , 'B11' ]).rename('NDSI' ); // 计算NDSVI var ndsvi = image.normalizedDifference(['B11' , 'B4' ]).rename('NDSVI' ); // 计算NDTI var ndti = image.normalizedDifference(['B11' , 'B12' ]).rename('NDTI' ); // 计算RENDVI var rendvi = image.normalizedDifference(['B8' , 'B5' ]).rename('RENDVI' ); // 计算REP var rep = image.expression( '(705 + 35 * ((0.5 * (B6 + B4) - B2) / (B5 - B2))) / 1000' , { 'B2' : image.select('B2' ), 'B4' : image.select('B4' ), 'B5' : image.select('B5' ), 'B6' : image.select('B6' ), 'B8' : image.select('B8' ) } ).rename('REP' ); // 添加所有计算的波段 return image.addBands([ndvi, evi, gndvi, savi, msavi2, ndwi, ndsi, ndsvi, ndti, rendvi, rep]); } // 5. 加载并预处理Sentinel-2 L2A影像集合 var sentinel = ee.ImageCollection("COPERNICUS/S2_SR" ); // 确保使用正确的Sentinel-2影像集合 var s2 = sentinel .filterBounds(region) .filterDate(startDate, endDate) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE' , 20 )) // 初步云量过滤 .map(maskS2clouds) .map(addSpectralIndices) .select(['B2' , 'B3' , 'B4' , 'B8' , 'NDVI' , 'EVI' , 'GNDVI' , 'SAVI' , 'MSAVI2' , 'NDWI' , 'NDSI' , 'NDSVI' , 'NDTI' , 'RENDVI' , 'REP' ]); // 6. 计算月度NDVI最大值 var months = ee.List.sequence(1 , 12 ); var monthlyMaxNDVI = months.map(function (month ) { var monthStart = ee.Date.fromYMD(year, month, 1 ); var monthEnd = monthStart.advance(1 , 'month' ); var monthlyNDVI = s2 .filterDate(monthStart, monthEnd) .select('NDVI' ) .max(); // 使用 ee.String 和 .cat() 正确拼接字符串 var bandName = ee.String('NDVI_month_' ).cat(ee.Number(month).format('%02d' )); return monthlyNDVI.rename(bandName); }); print(monthlyMaxNDVI,"monthlyMaxNDVI" ) // 将所有月份的最大NDVI合并为一个图像 var monthlyMaxNDVIImage = ee.Image.cat.apply(null , monthlyMaxNDVI) print(monthlyMaxNDVIImage,"monthlyMaxNDVIImage" ) // 7. 提取年度统计特征 var Year_Bands = ['B2' , 'B3' , 'B4' , 'B8' , 'NDVI' , 'EVI' , 'GNDVI' , 'SAVI' , 'MSAVI2' , 'NDWI'