GEE+本地XGboot分类

GEE+本地XGboot分类

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

但这个代码直接在云端上进行分类，GEE会爆内存，因此我准备把数据下载到本地，使用GPU加速进行XGboot提取耕地。

董老师的代码涉及到了100多个波段特征，我删减到了45个波段，然后分块进行了数据下载：

数据下载代码：

// ========================================
// 1. 初始化与区域选择
// ========================================// 合并训练数据
var trainTable = trainTable_crop.merge(trainTable_other);// 选择第一个区域作为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 = ee.List(['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 colStr = ee.Number(col).int();var rowStr = ee.Number(row).int();var tileId = ee.String(colStr).cat('_').cat(ee.String(rowStr));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');
print('Filtered tile count:', grid.size());// 打印tile ID以验证格式
print('Tile IDs:', grid.aggregate_array('tile'));// ========================================
// 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) {// 计算NDVIvar ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');// 计算EVIvar 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');// 计算GNDVIvar gndvi = image.normalizedDifference(['B8', 'B3']).rename('GNDVI');// 计算SAVIvar savi = image.expression('((NIR - RED) / (NIR + RED + 0.5)) * 1.5', {'NIR': image.select('B8'),'RED': image.select('B4')}).rename('SAVI');// 计算MSAVI2var 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');// 计算NDWIvar ndwi = image.normalizedDifference(['B3', 'B8']).rename('NDWI');// 计算NDSIvar ndsi = image.normalizedDifference(['B3', 'B11']).rename('NDSI');// 计算NDSVIvar ndsvi = image.normalizedDifference(['B11', 'B4']).rename('NDSVI');// 计算NDTIvar ndti = image.normalizedDifference(['B11', 'B12']).rename('NDTI');// 计算RENDVIvar rendvi = image.normalizedDifference(['B8', 'B5']).rename('RENDVI');// 计算REPvar 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(monthlyMaxNDVI.get(0));for (var i = 1; i < 12; i++) {monthlyMaxNDVIImage = monthlyMaxNDVIImage.addBands(ee.Image(monthlyMaxNDVI.get(i)));}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 = annualStats.clip(region);print(finalImage,"finalImage")                // 可视化示例（可选）// Map.addLayer(finalImage.select('NDVI_seasonal'), {min: 0, max: 1, palette: ['white', 'green']}, 'NDVI Seasonal');// 10. 导出数据到Google Drivevar 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_HENAN_FENQU1",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 xgboost as xgb
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report, roc_auc_score
from tqdm import tqdm  # 用于进度指示# 读取矢量数据
CROP_FILE = r"耕地样本点.shp"
OTHERS_FILE = r"非耕地样本点.shp"
TIF_PATH = r"mosaic_result_gdal.tif"cropland = gpd.read_file(CROP_FILE)
non_cropland = gpd.read_file(OTHERS_FILE)
cropland['crop'] = 1
non_cropland['crop'] = 0
samples = pd.concat([cropland, non_cropland], ignore_index=True)with rasterio.open(TIF_PATH) as src:band_count = src.countcoords = [(point.x, point.y) for point in samples.geometry]pixel_values = list(src.sample(coords))pixel_values = np.array(pixel_values)feature_columns = [f'band_{i+1}' for i in range(band_count)]
features = pd.DataFrame(pixel_values, columns=feature_columns)
features['crop'] = samples['crop'].values# 保存特征名称以供预测阶段使用
feature_names = feature_columns.copy()# 数据预处理
features.dropna(inplace=True)
X = features.drop('crop', axis=1)
y = features['crop']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y
)# 训练模型
dtrain = xgb.DMatrix(X_train, label=y_train, feature_names=feature_names)
dtest = xgb.DMatrix(X_test, label=y_test, feature_names=feature_names)params = {'objective': 'binary:logistic','tree_method': 'hist',         # 修改为 'hist''device': 'gpu',               # 添加 'device' 参数'eval_metric': 'auc','eta': 0.1,'max_depth': 10,'subsample': 0.8,'colsample_bytree': 0.8,'seed': 42
}evallist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 100print("开始训练模型...")
bst = xgb.train(params, dtrain, num_round, evallist, early_stopping_rounds=10, verbose_eval=True)
print("模型训练完成。\n")# 评估模型
print("开始评估模型...")
y_pred_prob = bst.predict(dtest)
y_pred = (y_pred_prob > 0.5).astype(int)
accuracy = accuracy_score(y_test, y_pred)
auc = roc_auc_score(y_test, y_pred_prob)
conf_matrix = confusion_matrix(y_test, y_pred)
report = classification_report(y_test, y_pred)print(f'Accuracy: {accuracy}')
print(f'AUC: {auc}')
print('Confusion Matrix:')
print(conf_matrix)
print('Classification Report:')
print(report)
print("模型评估完成。\n")# 应用模型进行栅格分类
print("开始进行栅格分类...")
with rasterio.open(TIF_PATH) as src:profile = src.profile.copy()profile.update(dtype=rasterio.uint8,count=1,compress='lzw')# 计算窗口总数用于进度指示windows = list(src.block_windows(1))total_windows = len(windows)with rasterio.open('classified.tif', 'w', **profile) as dst:for ji, window in tqdm(windows, total=total_windows, desc="栅格分类进度"):data = src.read(window=window)# data.shape = (bands, height, width)bands, height, width = data.shapedata = data.reshape(bands, -1).transpose()  # shape: (num_pixels, bands)# 创建 DataFrame 并赋予特征名称df = pd.DataFrame(data, columns=feature_names)# 创建 DMatrixdmatrix = xgb.DMatrix(df, feature_names=feature_names)# 预测predictions = bst.predict(dmatrix)predictions = (predictions > 0.5).astype(np.uint8)# 重塑为原窗口形状out_image = predictions.reshape(height, width)# 写入输出栅格dst.write(out_image, 1, window=window)
print("栅格分类完成。")

训练完成后，就开始分类了，就出结果了：

自此，从数据下载到分类处理完毕。

样本数据多的话，也可以考虑用CNN，但分类速度比不上xgboot。

参考：

You N , Dong J , Huang J ,et al.The 10-m crop type maps in Northeast China during 2017–2019[J].Scientific Data, 2021, 8(1).DOI:10.1038/s41597-021-00827-9.