GEE+本地XGboot分类

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

但这个代码直接在云端上进行分类，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', 'NDSI', 'NDSVI', 'NDTI', 'RENDVI', 'REP'];
  
  var annualStats = s2.select(Year_Bands)
      .reduce(ee.Reducer.mean()  

        .combine(ee.Reducer.max(), null, true)
        .combine(ee.Reducer.stdDev(), null, true));

  
  // 重命名年度统计特征的波段
  var statNames = ['mean', 'max', 'stdDev'];
  var newBandNames = [];
  Year_Bands.forEach(function(band) {
    statNames.forEach(function(stat) {
      newBandNames.push(band + '_' + stat);
    });
  });
  annualStats = annualStats.rename(newBandNames);
  
  // 将月度NDVI最大值和年度统计特征合并
  annualStats = ee.Image.cat([annualStats, monthlyMaxNDVIImage]);
  
 
  
  // 9. 合并所有特征
  var finalImage = ee.Image.cat([annualStats])
                  .clip(region);
                  
  // 可视化示例（可选）
  // Map.addLayer(finalImage.select('NDVI_seasonal'), {min: 0, max: 1, palette: ['white', 'green']}, 'NDVI Seasonal');
  
  // 10. 导出数据到Google Drive
  var output_name='tile_' + tileID.getInfo()
  var name2=output_name.replace('.', '').replace('.', '')
  print(finalImage.toFloat())
  Export.image.toDrive({
    image: finalImage.toFloat(),
    description:  name2,
    scale: 10,
    folder: "download_tiles",
    region: region,  
    maxPixels: 1e13
  });
}

// ========================================
// 4. 应用函数到每个块
// ========================================

// 注意：Google Earth Engine 同时只能运行有限的Export任务（通常为3个）。
// 因此，建议分批次运行或手动触发每个块的导出任务。

// 将网格转换为特征集合列表
var gridFeatures = grid.toList(grid.size());

// 获取总块数
var totalTiles = grid.size().getInfo();

// 定义每批次导出的数量（如果需要批量控制，可以在这里调整）
var batchSize = 1;

// 处理并导出每个块
// 注意：Google Earth Engine 不支持并行启动大量导出任务，请手动管理导出任务
gridFeatures.evaluate(function(list) {
  list.forEach(function(feature) {
    processAndExport(ee.Feature(feature));
  });
});

// 打印总块数和导出说明
print('Total tiles:', totalTiles);
print('导出已启动。请在任务管理器中检查导出状态。');

然后下载完成后，用gdal做一下镶嵌（设置tile为256，LZW压缩），波段太多，导致数据非常大。最好再做一个金字塔

import os
from osgeo import gdal

# 输入和输出路径
input_dir = r"几十个波段数据"
output_file = "mosaic_result_gdal.tif"

# 获取所有tif文件
tif_files = []
for file in os.listdir(input_dir):
    if file.endswith('.tif'):
        tif_files.append(os.path.join(input_dir, file))

# 构建VRT
vrt = gdal.BuildVRT("temp.vrt", tif_files)
vrt = None

# 转换VRT为GeoTiff
gdal.Translate(
    output_file,
    "temp.vrt",
    format="GTiff",
    creationOptions=[
        "COMPRESS=LZW",
        "TILED=YES",
        "BLOCKXSIZE=256",
        "BLOCKYSIZE=256",
        "BIGTIFF=YES"
    ]
)

镶嵌完，可以放进GIS软件中查看一下。

数据分类

在此之前，需要先准备点数据，我是准备了两个点数据矢量（耕地矢量和非耕地矢量），字段属性crop为1代表耕地，0代表非耕地。如果你是做多类别，你可以多做几个矢量。

然后开始安装环境：

(1)安装CUDA，用GPU加速运行，也可以CPU，都差不多，xgboot计算量不大；

(2)安装conda，然后使用下面的命令安装环境：

conda create --prefix D:\conda_ENV\xgboot_env python=3.10
conda activate D:\conda_ENV\xgboot_env
conda install -c conda-forge numpy pandas geopandas rasterio scikit-learn tqdm

然后就可以开始分类了，代码如下：

import geopandas as gpd
import rasterio
from rasterio.sample import sample_gen
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import