专栏名称: 锐多宝
遥感技术教程、资讯与前沿论文
目录
相关文章推荐
51好读  ›  专栏  ›  锐多宝

GEE+本地XGboot分类

锐多宝  · 公众号  ·  · 2024-12-18 23:50

正文

我想做提取耕地提取,想到了一篇董金玮老师的一篇论文,这个论文是先提取的耕地,再做作物分类,论文代码是开源的。基于这个代码下载数据。

但这个代码直接在云端上进行分类,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, 11);
  var endDate = ee.Date.fromYMD(year, 1231);
  
  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(112);
  
  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'






请到「今天看啥」查看全文